3D olfactory bulb: operators (Migliore et al, 2015)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:168591
"... Using a 3D model of mitral and granule cell interactions supported by experimental findings, combined with a matrix-based representation of glomerular operations, we identify the mechanisms for forming one or more glomerular units in response to a given odor, how and to what extent the glomerular units interfere or interact with each other during learning, their computational role within the olfactory bulb microcircuit, and how their actions can be formalized into a theoretical framework in which the olfactory bulb can be considered to contain "odor operators" unique to each individual. ..."
Reference:
1 . Migliore M, Cavarretta F, Marasco A, Tulumello E, Hines ML, Shepherd GM (2015) Synaptic clusters function as odor operators in the olfactory bulb. Proc Natl Acad Sci U S A 112:8499-504 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Realistic Network;
Brain Region(s)/Organism:
Cell Type(s): Olfactory bulb main mitral GLU cell; Olfactory bulb main interneuron granule MC GABA cell;
Channel(s): I Na,t; I A; I K;
Gap Junctions:
Receptor(s): AMPA; NMDA; Gaba;
Gene(s):
Transmitter(s): Gaba; Glutamate;
Simulation Environment: NEURON; Python;
Model Concept(s): Activity Patterns; Dendritic Action Potentials; Active Dendrites; Synaptic Plasticity; Action Potentials; Synaptic Integration; Unsupervised Learning; Sensory processing; Olfaction;
Implementer(s): Migliore, Michele [Michele.Migliore at Yale.edu]; Cavarretta, Francesco [francescocavarretta at hotmail.it];
Search NeuronDB for information about:  Olfactory bulb main mitral GLU cell; Olfactory bulb main interneuron granule MC GABA cell; AMPA; NMDA; Gaba; I Na,t; I A; I K; Gaba; Glutamate;
/
figure1eBulb3D
readme.html
ampanmda.mod *
distrt.mod *
fi.mod *
fi_stdp.mod *
kamt.mod *
kdrmt.mod *
naxn.mod *
ThreshDetect.mod *
.hg_archival.txt
all2all.py *
balance.py *
bindict.py
binsave.py
binspikes.py
BulbSurf.py
catfiles.sh
colors.py *
common.py
complexity.py *
custom_params.py *
customsim.py
destroy_model.py *
determine_connections.py
distribute.py *
falsegloms.txt
fixnseg.hoc *
g37e1i002.py
gidfunc.py *
Glom.py *
granule.hoc *
granules.py
grow.py
input-odors.txt *
loadbalutil.py *
lpt.py *
m2g_connections.py
mayasyn.py
mgrs.py
misc.py
mitral.hoc *
mkdict.py
mkmitral.py
modeldata.py *
multisplit_distrib.py *
net_mitral_centric.py
odordisp.py *
odors.py *
odorstim.py
params.py
parrun.py
realgloms.txt *
realSoma.py *
runsim.py
spike2file.hoc *
split.py *
util.py *
vrecord.py
weightsave.py *
                            
# -*- coding: cp1252 -*-
from copy import copy
import realSoma
from math import sqrt, sin, cos, pi
from misc import *
from params import *

import params

class Section:
        def __init__(self):
                self.index = -1
                self.points = []
                self.sons = []
                self.parent = None
                
class Neuron:
        def __init__(self): 
                self.dend = []
                self.apic = None
                self.tuft = []
                self.soma = None
class Extreme:
        SOMA = -1; DENDRITE = 0; APICAL = 1; TUFT = 2
        def __init__(self):
                self.sec = None
                self.phi = 0.
                self.theta = 0.
                self.dist = 0.
                self.limit = 0.
                self.extr_type = Extreme.SOMA
                self.can_bifurc = False
                self.bif_dist = 0.
                self.bif_limit = 0.
                self.depth = 0
                self.basePhi = 0.
                self.baseTheta = 0.
                



# find intersection with Ellipsoid
def EllipsoidLineIntersec(u, p, axis):
  return Ellipsoid(params.bulbCenter,axis).intersection(p,u)

# Ellipsoidd with intersection
def EllipsoidIntersec(p, axis):
  e = Ellipsoid(params.bulbCenter, axis)
  p = e.project(p)
  lamb, phi = e.toElliptical(p)[1:]
  return p, pi*0.5-lamb, phi
  



new_version = False

if new_version:

        
  ''' init the soma position '''
  def gen_soma_pos(r, glompos, mgid):
    nmg=params.Nmitral_per_glom
    glomproj, phi_base, theta_base = EllipsoidIntersec(glompos, params.somaAxis[0])
    phi_base = 0.5*pi - phi_base
    phi_phase = 2*pi/nmg*(mgid%nmg)
    
    # generate a new position
    def gen_pos():
      phi = mgid%nmg * 2*pi / nmg + r.uniform(-0.4*pi/nmg, 0.4*pi/nmg) + phi_phase
      theta = r.normal(0.5*pi, 0.05*(pi*0.5)**2)
      rho = r.uniform(0.2, 1) * params.GLOM_DIST
      phi, theta = convert_direction(phi, theta, phi_base, theta_base)
      return Spherical.xyz(rho, phi, theta, glomproj)

    # check if the position is good
    def in_soma_zone(pos):
      soma_inf = Ellipsoid(bulbCenter, somaAxis[0])
      soma_sup = Ellipsoid(bulbCenter, somaAxis[1])
      return soma_inf.normalRadius(pos) > 1. and soma_sup.normalRadius(pos) < 1.

    while True:
      pos = gen_pos()
      if in_soma_zone(pos):
        break

    return pos

  ''' init the soma '''
  def mk_soma(mgid, nrn):
    sec = Section()
    glompos = params.glomRealCoords[int(mgid / params.Nmitral_per_glom)]
    r = params.ranstream(mgid, params.stream_soma)
    somapos = gen_soma_pos(r, glompos, mgid)

    import realSoma
    index = int(r.discunif(0, realSoma.N_SOMA-1))
    sec.points = realSoma.realSoma(index, somapos)
    nrn.soma = sec

    return somapos


  def can_change_type(extr, glompos):
    if extr.extr_type == Extreme.APICAL:
      return distance(extr.sec.points[-1], glompos) < params.GLOM_RADIUS
    return False

  ''' check the barrier '''
  def feasible(p, extr, nrn, glompos):
    if extr.extr_type == Extreme.APICAL:
      return distance(p, glompos) < distance(extr.sec.points[-1], glompos)
    elif extr.extr_type == Extreme.TUFT:
      d = distance(p, glompos)/params.GLOM_RADIUS
      return d <= 1 and d >= 0.75
    return True
    
  # noise
  def noise(r, b,c=0,d=0):
    return rLaplace(r, 0., b)

  # generate a walk, contains the growing rules.
  def genWalk(r, extr, glomPos):


          
          r = r[extr.extr_type]
          rho = WALK_RHO
          phi = extr.phi
          theta = extr.theta
          diam = 1.                        
          pts = extr.sec.points
          
          if extr.extr_type == Extreme.APICAL:
                  diam = APIC_DIAM
                  dphi = noise(r, params.NS_PHI_B)
                  dtheta = noise(r, params.NS_THETA_B)
                  phi += dphi
                  theta += dtheta
                  absphi, abstheta = convert_direction(phi, theta, extr.basePhi, extr.baseTheta)
                  p = Spherical.xyz(rho, absphi, abstheta, extr.sec.points[-1])
                  p += [ diam ]
                  return p, phi, theta
          elif extr.extr_type == Extreme.TUFT:
                  #def correctTuft(norm, p):
                  #        _rho, phi, theta = Spherical.to(getP(norm, versor(glomPos, p), p), extr.sec.points[-1])
                  #        return phi, theta
                  #
                  #
                  #_p = Spherical.xyz(rho, phi, theta, pts[-1])
                  #dglom = distance(_p, glomPos)
                  #if dglom > 0.9 * GLOM_RADIUS:
                  #        phi, theta = correctTuft(dglom - .9 * GLOM_RADIUS, _p)
                  #elif dglom < 0.6 * GLOM_RADIUS:
                  #        phi, theta = correctTuft(dglom - .6 * GLOM_RADIUS, _p)
                  ##
                  diam = TUFT_DIAM
                  dphi = noise(r, params.NS_PHI_B)
                  dtheta = noise(r, params.NS_THETA_B)
                  phi += dphi
                  theta += dtheta
                  absphi, abstheta = convert_direction(phi, theta, extr.basePhi, extr.baseTheta)
                  p = Spherical.xyz(rho, absphi, abstheta, extr.sec.points[-1])
                  p += [ diam ]
          elif extr.extr_type == Extreme.DENDRITE:
                  ## diam
                  diam = gen_dend_diam(extr.dist)
                  pts = extr.sec.points

                  # first curvature
                  #if len(pts) < 5 and extr.depth == 0:
                    #      theta += (pi / 3. - pi / 15) / 4
                  #        _phi, _theta = ConvertDirection(phi, theta, extr.basePhi, extr.baseTheta, True)
                   #       p = Spherical.xyz(rho, _phi, _theta, pts[-1])

                  
                  # simulate Ellipsoid pression
                  def EllipsoidPression(p, axis, up, k):
                          e = Ellipsoid(bulbCenter, axis)
                          h = e.normalRadius(p)

                          # check if the pression if from up or down to the surface
                          if up:
                                  h_check = h > 1.
                                  h_diff = h
                          else:
                                  h_check = h < 1.
                                  h_diff = 1. - h
                                  
                          q = None
                          if h_check:
                                  q, _lamb, _phi = EllipsoidIntersec(p, axis)
                          else:
                                  _h, _lamb, _phi = e.toElliptical(p)
                                  F = h_diff * k
                                  q = [ -F * sin(_lamb) * cos(_phi) + p[0], -F * cos(_lamb) * cos(_phi) + p[1], -F * sin(_phi) + p[2] ]
                                  
                          vNew = versor(q, pts[-1])
                          _p = getP(rho, vNew, pts[-1])
                          _rho, _phi, _theta = Spherical.to(_p, pts[-1])
                          return _p, _phi, _theta

                  # check
                  _phi, _theta = convert_direction(phi, theta, extr.basePhi, extr.baseTheta)
                  p = Spherical.xyz(rho, _phi, _theta, pts[-1])
                  if not inSomaZone(pts[-1]):
                          e = Ellipsoid(bulbCenter, bulbAxis)
                          if inSomaZone(p):
                                  p, _phi, _theta = EllipsoidPression(p, somaAxis[1], True, 0.)
                                  phi, theta = convert_direction(_phi, _theta, extr.basePhi, extr.baseTheta, True)
                          elif e.normalRadius(pts[-1]) < e.normalRadius(p):
                                  #p, _phi, _theta = EllipsoidPression(p, bulbAxis, True, .6)
                                  p, _phi, _theta = EllipsoidPression(p, bulbAxis, True, GROW_RESISTANCE)
                                  phi, theta = convert_direction(_phi, _theta, extr.basePhi, extr.baseTheta, True)

                  _phi += noise(r, NS_PHI_B, NS_PHI_MIN, NS_PHI_MAX)
                  _theta += noise(r, NS_THETA_B, NS_THETA_MIN, NS_THETA_MAX)
                  p = Spherical.xyz(rho, _phi, _theta, pts[-1])                
                  p.append(diam)
                  return p, phi, theta



                  
          # add the noise        
          phi += noise(r, NS_PHI_B, NS_PHI_MIN, NS_PHI_MAX)
          theta += noise(r, NS_THETA_B, NS_THETA_MIN, NS_THETA_MAX)
          p = Spherical.xyz(rho, phi, theta, extr.sec.points[-1])
          p.append(diam)                               
          return p, phi, theta

  def canDelete(extr, flag_fail):
    if extr.extr_type != Extreme.APICAL:
      return (((extr.dist >= extr.limit or abs(extr.dist - extr.limit) < WALK_RHO) and not extr.can_bifurc) or flag_fail)
    
    return False

  def canBifurcate(extr):
    if extr.extr_type == Extreme.DENDRITE:
      return extr.can_bifurc and (extr.bif_dist >= extr.bif_limit or abs(extr.bif_dist - extr.bif_limit) < WALK_RHO)
    
    return False


                                                                                       
  def updateExtr(extr, p, phi, theta):
          d = distance(p, extr.sec.points[-1])
          extr.dist += d
          extr.bif_dist += d
          extr.sec.points += [ p ]                
          extr.phi = phi
          extr.theta = theta
          
  def mkDendrite(r, depth, parent):
          sec = Section()
          sec.parent = parent
          if parent: parent.sons.append(sec)
          ##
          extr = Extreme()
          extr.extr_type = Extreme.DENDRITE
          extr.sec = sec
          
          # bifurcation decision
          extr.can_bifurc = depth < len(BIFURC_PROB) and r.uniform(0, 1) <= BIFURC_PROB[depth] and extr.dist < (DEND_LEN_MAX - DEND_LEN_MIN)
          extr.depth = depth
         #extr.can_bifurc = False
          if extr.can_bifurc:
                  if depth <= 0:
                          min_len = BIFURC_LEN_MIN_1
                          max_len = BIFURC_LEN_MAX_1
                          mu = BIFURC_LEN_MU_1
                  else:
                          min_len = BIFURC_LEN_MIN_2
                          max_len = BIFURC_LEN_MAX_2
                          mu = BIFURC_LEN_MU_2
                  bl = r.negexp(mu)
                  while bl < min_len or bl > max_len: bl = r.repick()
                  extr.bif_limit = bl
          else:
                  minLen = DEND_LEN_MIN
                  if extr.dist > 0: minLen = extr.dist + DEND_LEN_MIN
                  
                  dendLen = r.normal(DEND_LEN_MU, DEND_LEN_VAR)
                  while dendLen < minLen or dendLen > DEND_LEN_MAX: dendLen = r.repick()
                  extr.limit = dendLen
                  
          return sec, extr




  def mkApical(mid, somaPos, glomPos, nrn, extrLs):
          sec = Section()
          sec.points = [ copy(somaPos) ]
          sec.points[-1].append(APIC_DIAM)
          sec.parent = nrn.soma
          nrn.soma.sons.append(sec)
          nrn.apic = sec
          extr = Extreme()
          extr.sec = sec
          #rho, phi, theta = Spherical.to(glomPos, somaPos)
          extr.phi = 0; extr.theta = 0
          extr.extr_type = Extreme.APICAL
          extr.limit = APIC_LEN_MAX
          extrLs.append(extr)



  def mkTuft(r, extrLs, nrn, pos, glomPos,basePhi,baseTheta):
          # orientation respect glomerulus
          _rho, gl_phi, gl_theta = Spherical.to(glomPos, pos)
          
          # make tuft        
          TUFTS = int(r.discunif(N_MIN_TUFT, N_MAX_TUFT))
          for i in range(TUFTS):
                  sec = Section()
                  sec.index = i
                  sec.parent = nrn.apic
                  nrn.apic.sons.append(sec)
                  nrn.tuft.append(sec)
                  sec.points = [ copy(nrn.apic.points[-1]) ]
                  sec.points[-1][-1] = TUFT_DIAM
                  extr = Extreme()
                  extr.sec = sec
                  extr.phi, extr.theta = convert_direction(i * 2 * pi / TUFTS * (1 + 0.75 * r.uniform(-0.5 / TUFTS, 0.5 / TUFTS)), pi / 4, gl_phi, gl_theta)
                  extr.limit = r.uniform(TUFT_MIN_LEN, TUFT_MAX_LEN)
                  extr.extr_type = Extreme.TUFT
                  extrLs.append(extr)
                  extrLs[-1].basePhi=basePhi;extrLs[-1].baseTheta=baseTheta
          


  # change extremity, especially at the end of apical                                                             
  def change_type(r, extrLs, i, nrn, glomPos):
                                   
          mkTuft(r[Extreme.TUFT], extrLs, nrn, extrLs[i].sec.points[-1], glomPos,extrLs[i].basePhi,extrLs[i].baseTheta)
          


  # make a bifurcation
  def Bifurcate(r, extrLs, i, nrn):
          r = r[Extreme.DENDRITE]
          def get_phi():
                  phi = r.negexp(BIFURC_PHI_MU)
                  while phi < BIFURC_PHI_MIN or phi > BIFURC_PHI_MAX: phi = r.repick()
                  return phi

          
          sec1, extr1 = mkDendrite(r, extrLs[i].depth + 1, extrLs[i].sec)
          sec1.points = [ copy(extrLs[i].sec.points[-1]) ]
          sec1.index = len(nrn.dend)
          nrn.dend.append(sec1)
          
          extr1.phi = extrLs[i].phi + get_phi(); extr1.theta = extrLs[i].theta
          extr1.dist = extrLs[i].dist
          extr1.basePhi = extrLs[i].basePhi; extr1.baseTheta = extrLs[i].baseTheta
          extrLs.append(extr1)

          sec2, extr2 = mkDendrite(r, extrLs[i].depth + 1, extrLs[i].sec)
          sec2.points = [ copy(extrLs[i].sec.points[-1]) ]
          sec2.index = len(nrn.dend)
          nrn.dend.append(sec2)
          extr2.basePhi = extrLs[i].basePhi; extr2.baseTheta = extrLs[i].baseTheta
          extr2.phi = extrLs[i].phi - get_phi(); extr2.theta = extrLs[i].theta
          extr2.dist = extrLs[i].dist
          extrLs.append(extr2)
          

  def inSomaZone(p):
          somaInf = Ellipsoid(bulbCenter, somaAxis[0]); somaSup = Ellipsoid(bulbCenter, somaAxis[1])
          return somaInf.normalRadius(p) >= 1. and somaSup.normalRadius(p) < 1.


  # initialization of algorithm
  def initGrow(mid, noDend):
          
          r = [ ranstream(mid, stream_dend), ranstream(mid, stream_apic), ranstream(mid, stream_tuft) ]
          extrLs = [ ]
          nrn = Neuron()
          somaPos, glomPos = mkSoma(mid, nrn)
          
          mkApical(r[Extreme.APICAL], somaPos, glomPos, nrn, extrLs)

          # make dendrites.
          if not noDend:
                  # get basic inclination
                  somaLvl = Ellipsoid(bulbCenter, somaAxis[0])
                  h, phi_base, theta_base = somaLvl.toElliptical(somaPos)
                  theta_base = pi / 2 - theta_base; phi_base = pi / 2 - phi_base
                  _r = r[Extreme.DENDRITE]
                  
                  ###
                  DENDRITES = int(_r.negexp(N_MEAN_DEND - N_MIN_DEND)) + N_MIN_DEND
                  while DENDRITES > N_MAX_DEND: DENDRITES = int(_r.negexp(N_MEAN_DEND - N_MIN_DEND)) + N_MIN_DEND
                  nmg=params.Nmitral_per_glom
                  dphi=2*pi/DENDRITES
                  for i in range(DENDRITES):
                          phiPhase = i*dphi + 2*pi/nmg*(mid%nmg)
                          
                          sec, extr = mkDendrite(_r, 0, nrn.soma)
                          p = copy(somaPos); p.append(gen_dend_diam(0.))
                          sec.points = [ p ]
                          sec.index = i              
                          nrn.dend.append(sec)
                          
                          ## initial direction
                          phi = phiPhase + _r.uniform(-0.4,0.4)*dphi
                          
                          theta = pi / 3. # - _r.uniform(pi / 15, pi / 12)
                          newPhi, newTheta = convert_direction(phi, theta, phi_base, theta_base)
                          extr.phi = phi; extr.theta = theta
                          extr.basePhi = phi_base; extr.baseTheta = theta_base
                          extrLs.append(extr)
          return r, nrn, extrLs, glomPos


  
  def Grow(mid, noDend):
          r = nrn = extrLs = glomPos = None
          r, nrn, extrLs, glomPos = initGrow(mid, noDend) # initialization        
          
          it = 0
          while it < GROW_MAX_ITERATIONS and len(extrLs) > 0:
                  i = 0
                  while i < len(extrLs):
                          j = 0
                          while j < GROW_MAX_ATTEMPTS:
                                  p, phi, theta = genWalk(r, extrLs[i], glomPos)
                                  
                                  if feasible(p, extrLs[i], nrn, glomPos):
                                          updateExtr(extrLs[i], p, phi, theta)
                                          break
                                  j += 1
                                           
                          if canDelete(extrLs[i], j >= GROW_MAX_ATTEMPTS):
                                  if len(extrLs[i].sec.points) <= 1 and extrLs[i].extr_type == Extreme.DENDRITE:
                                          if extrLs[i].sec.index < (len(nrn.dend) - 1):
                                                  for q in range(extrLs[i].sec.index + 1, len(nrn.dend)):
                                                          nrn.dend[q].index -= 1
                                          del nrn.dend[extrLs[i].sec.index] 
                                  
                                  del extrLs[i]        
                          elif can_change_type(extrLs[i], glomPos):
                                  change_type(r, extrLs, i, nrn, glomPos)
                                  del extrLs[i]
                          elif canBifurcate(extrLs[i]):
                                 Bifurcate(r, extrLs, i, nrn)
                                 del extrLs[i]
                          else:
                                  i += 1
                  it += 1
                  
          return nrn

  def genMitral(mid): return Grow(mid, False)
  def genSomaApicalTuft(mid): return Grow(mid, True)

else:

  def inSomaZone(p):
          somaInf = Ellipsoid(bulbCenter, somaAxis[0]); somaSup = Ellipsoid(bulbCenter, somaAxis[1])
          return somaInf.normalRadius(p) >= 1. and somaSup.normalRadius(p) < 1.
      
  ''' init the soma position '''
  def gen_soma_pos(r, glompos, mgid):
    nmg=params.Nmitral_per_glom
    glomproj1, phi_base, theta_base = EllipsoidIntersec(glompos, params.somaAxis[0])
    glomproj2 = EllipsoidIntersec(glompos, params.somaAxis[1])[0]
    glomproj = centroid([ glomproj1, glomproj2 ])
    
    phi_base = 0.5*pi - phi_base
    
    # generate a new position
    def gen_pos():
      phi = 2*pi/nmg*(mgid%nmg) +2*pi/nmg*r.uniform(-0.4,0.4)
      
      theta = r.normal(0.5*pi,0.05*(0.5*pi)**2)
      rho = r.uniform(0.5, 1) * params.GLOM_DIST
      phi, theta = convert_direction(phi, theta, phi_base, theta_base)
      return Spherical.xyz(rho, phi, theta, glomproj)

    # check if the position is good
    def in_soma_zone(pos):
      soma_inf = Ellipsoid(bulbCenter, somaAxis[0])
      soma_sup = Ellipsoid(bulbCenter, somaAxis[1])
      return soma_inf.normalRadius(pos) > 1. and soma_sup.normalRadius(pos) < 1.

    while True:
      pos = gen_pos()
      if in_soma_zone(pos):
        break

    return pos

  ''' init the soma '''
  def mk_soma(mgid, nrn):
    sec = Section()
    glompos = params.glomRealCoords[int(mgid / params.Nmitral_per_glom)]
    r = params.ranstream(mgid, params.stream_soma)
    somapos = gen_soma_pos(r, glompos, mgid)

    import realSoma
    index = int(r.discunif(0, realSoma.N_SOMA-1))
    sec.points = realSoma.realSoma(index, somapos)
    nrn.soma = sec

    return somapos

  ''' check the barrier '''
  def feasible(p, extr, nrn, glompos):
    if extr.extr_type == Extreme.APICAL:
      return distance(p, glompos) < distance(extr.sec.points[-1], glompos)
    elif extr.extr_type == Extreme.TUFT:
      d = distance(p, glompos)/params.GLOM_RADIUS
      return d <= 1 and d >= 0.75
    return True

  # noise
  def noise(r, b, minval, maxval):
          ns = rLaplace(r, 0., b)
          while ns < minval or ns > maxval: ns = rLaplace(r, 0., b)
          return ns

  # generate a walk, contains the growing rules.
  def genWalk(r, extr, glomPos):


          
          r = r[extr.extr_type]
          rho = WALK_RHO; phi = extr.phi; theta = extr.theta; diam = 1.                        
          pts = extr.sec.points
          if extr.extr_type == Extreme.APICAL:
                  _rho, phi, theta = Spherical.to(glomPos, extr.sec.points[-1])
                  diam = APIC_DIAM                
          elif extr.extr_type == Extreme.TUFT:
                  def correctTuft(norm, p):
                          _rho, phi, theta = Spherical.to(getP(norm, versor(glomPos, p), p), extr.sec.points[-1])
                          return phi, theta
                  #
                  
                  _p = Spherical.xyz(rho, phi, theta, pts[-1])
                  dglom = distance(_p, glomPos)
                  if dglom > 0.9 * GLOM_RADIUS:
                          phi, theta = correctTuft(dglom - .9 * GLOM_RADIUS, _p)
                  elif dglom < 0.6 * GLOM_RADIUS:
                          phi, theta = correctTuft(dglom - .6 * GLOM_RADIUS, _p)
                  ##
                  diam = TUFT_DIAM
          elif extr.extr_type == Extreme.DENDRITE:
                  ## diam
                  diam = gen_dend_diam(extr.dist)
                  pts = extr.sec.points

                  # first curvature
                  #if len(pts) < 5 and extr.depth == 0:
                    #      theta += (pi / 3. - pi / 15) / 4
                  #        _phi, _theta = ConvertDirection(phi, theta, extr.basePhi, extr.baseTheta, True)
                   #       p = Spherical.xyz(rho, _phi, _theta, pts[-1])

                  
                  # simulate Ellipsoid pression
                  def EllipsoidPression(p, axis, up, k):
                          e = Ellipsoid(bulbCenter, axis)
                          h = e.normalRadius(p)

                          # check if the pression if from up or down to the surface
                          if up:
                                  h_check = h > 1.
                                  h_diff = h
                          else:
                                  h_check = h < 1.
                                  h_diff = 1. - h
                                  
                          q = None
                          if h_check:
                                  q, _lamb, _phi = EllipsoidIntersec(p, axis)
                          else:
                                  _h, _lamb, _phi = e.toElliptical(p)
                                  F = h_diff * k
                                  q = [ -F * sin(_lamb) * cos(_phi) + p[0], -F * cos(_lamb) * cos(_phi) + p[1], -F * sin(_phi) + p[2] ]
                                  
                          vNew = versor(q, pts[-1])
                          _p = getP(rho, vNew, pts[-1])
                          _rho, _phi, _theta = Spherical.to(_p, pts[-1])
                          return _p, _phi, _theta

                  # check
                  _phi, _theta = convert_direction(phi, theta, extr.basePhi, extr.baseTheta)
                  p = Spherical.xyz(rho, _phi, _theta, pts[-1])
                  if not inSomaZone(pts[-1]):
                          e = Ellipsoid(bulbCenter, bulbAxis)
                          if inSomaZone(p):
                                  p, _phi, _theta = EllipsoidPression(p, somaAxis[1], True, 0.)
                                  phi, theta = convert_direction(_phi, _theta, extr.basePhi, extr.baseTheta, True)
                          elif e.normalRadius(pts[-1]) < e.normalRadius(p):
                                  #p, _phi, _theta = EllipsoidPression(p, bulbAxis, True, .6)
                                  p, _phi, _theta = EllipsoidPression(p, bulbAxis, True, GROW_RESISTANCE)
                                  phi, theta = convert_direction(_phi, _theta, extr.basePhi, extr.baseTheta, True)

                  _phi += noise(r, NS_PHI_B, NS_PHI_MIN, NS_PHI_MAX)
                  _theta += noise(r, NS_THETA_B, NS_THETA_MIN, NS_THETA_MAX)
                  p = Spherical.xyz(rho, _phi, _theta, pts[-1])                
                  p.append(diam)
                  return p, phi, theta



                  
          # add the noise        
          phi += noise(r, NS_PHI_B, NS_PHI_MIN, NS_PHI_MAX)
          theta += noise(r, NS_THETA_B, NS_THETA_MIN, NS_THETA_MAX)
          p = Spherical.xyz(rho, phi, theta, extr.sec.points[-1])
          p.append(diam)                               
          return p, phi, theta

  def canDelete(extr, flag_fail):
    if extr.extr_type != Extreme.APICAL:
      return (((extr.dist >= extr.limit or abs(extr.dist - extr.limit) < WALK_RHO) and not extr.can_bifurc) or flag_fail)
    
    return False

  def canBifurcate(extr):
    if extr.extr_type == Extreme.DENDRITE:
      return extr.can_bifurc and (extr.bif_dist >= extr.bif_limit or abs(extr.bif_dist - extr.bif_limit) < WALK_RHO)
    
    return False

  def canChangeTypeOfExtr(extr, glomPos): return extr.extr_type == Extreme.APICAL and (extr.dist >= extr.limit or distance(extr.sec.points[-1], glomPos) < GLOM_RADIUS)
                                                                                       
  def updateExtr(extr, p, phi, theta):
          d = distance(p, extr.sec.points[-1])
          extr.dist += d
          extr.bif_dist += d
          extr.sec.points += [ p ]                
          extr.phi = phi
          extr.theta = theta
          
  def mkDendrite(r, depth, parent):
          sec = Section()
          sec.parent = parent
          if parent: parent.sons.append(sec)
          ##
          extr = Extreme()
          extr.extr_type = Extreme.DENDRITE
          extr.sec = sec
          
          # bifurcation decision
          extr.can_bifurc = depth < len(BIFURC_PROB) and r.uniform(0, 1) <= BIFURC_PROB[depth] and extr.dist < (DEND_LEN_MAX - DEND_LEN_MIN)
          extr.depth = depth
         #extr.can_bifurc = False
          if extr.can_bifurc:
                  if depth <= 0:
                          min_len = BIFURC_LEN_MIN_1
                          max_len = BIFURC_LEN_MAX_1
                          mu = BIFURC_LEN_MU_1
                  else:
                          min_len = BIFURC_LEN_MIN_2
                          max_len = BIFURC_LEN_MAX_2
                          mu = BIFURC_LEN_MU_2
                  bl = r.negexp(mu)
                  while bl < min_len or bl > max_len: bl = r.repick()
                  extr.bif_limit = bl
          else:
                  minLen = DEND_LEN_MIN
                  if extr.dist > 0: minLen = extr.dist + DEND_LEN_MIN
                  
                  dendLen = r.normal(DEND_LEN_MU, DEND_LEN_VAR)
                  while dendLen < minLen or dendLen > DEND_LEN_MAX: dendLen = r.repick()
                  extr.limit = dendLen
                  
          return sec, extr



  def mkApical(mid, somaPos, glomPos, nrn, extrLs):
          sec = Section()
          sec.points = [ copy(somaPos) ]
          sec.points[-1].append(APIC_DIAM)
          sec.parent = nrn.soma
          nrn.soma.sons.append(sec)
          nrn.apic = sec
          extr = Extreme()
          extr.sec = sec
          rho, phi, theta = Spherical.to(glomPos, somaPos)
          extr.phi = phi; extr.theta = theta
          extr.extr_type = Extreme.APICAL
          extr.limit = APIC_LEN_MAX
          extrLs.append(extr)



  def mkTuft(r, extrLs, nrn, pos, glomPos):
          # orientation respect glomerulus
          _rho, gl_phi, gl_theta = Spherical.to(glomPos, pos)
          
          # make tuft        
          TUFTS = int(r.discunif(N_MIN_TUFT, N_MAX_TUFT))
          for i in range(TUFTS):
                  sec = Section()
                  sec.index = i
                  sec.parent = nrn.apic
                  nrn.apic.sons.append(sec)
                  nrn.tuft.append(sec)
                  sec.points = [ copy(nrn.apic.points[-1]) ]
                  sec.points[-1][-1] = TUFT_DIAM
                  extr = Extreme()
                  extr.sec = sec
                  extr.phi, extr.theta = convert_direction(i * 2 * pi / TUFTS * (1 + 0.75 * r.uniform(-0.5 / TUFTS, 0.5 / TUFTS)), pi / 4, gl_phi, gl_theta)
                  extr.limit = r.uniform(TUFT_MIN_LEN, TUFT_MAX_LEN)
                  extr.extr_type = Extreme.TUFT
                  extrLs.append(extr)
          
          

  # change extremity, especially at the end of apical                                                             
  def changeTypeExtr(r, extrLs, i, nrn, glomPos):
          
          if distance(extrLs[i].sec.points[-1], glomPos) != GLOM_RADIUS:
                  pos = getP(distance(extrLs[i].sec.points[-1], glomPos) - GLOM_RADIUS,
                             versor(glomPos, extrLs[i].sec.points[-1]),
                             extrLs[i].sec.points[-1])
                  stretchSection(extrLs[i].sec.points, pos)
                  
                  
          mkTuft(r[Extreme.TUFT], extrLs, nrn, pos, glomPos)
          
                  
  # make a bifurcation
  def Bifurcate(r, extrLs, i, nrn):
          r = r[Extreme.DENDRITE]
          def get_phi():
                  phi = r.negexp(BIFURC_PHI_MU)
                  while phi < BIFURC_PHI_MIN or phi > BIFURC_PHI_MAX: phi = r.repick()
                  return phi
          
          sec1, extr1 = mkDendrite(r, extrLs[i].depth + 1, extrLs[i].sec)
          sec1.points = [ copy(extrLs[i].sec.points[-1]) ]
          sec1.index = len(nrn.dend)
          nrn.dend.append(sec1)
          
          extr1.phi = extrLs[i].phi + get_phi(); extr1.theta = extrLs[i].theta
          extr1.dist = extrLs[i].dist
          extr1.basePhi = extrLs[i].basePhi; extr1.baseTheta = extrLs[i].baseTheta
          extrLs.append(extr1)

          sec2, extr2 = mkDendrite(r, extrLs[i].depth + 1, extrLs[i].sec)
          sec2.points = [ copy(extrLs[i].sec.points[-1]) ]
          sec2.index = len(nrn.dend)
          nrn.dend.append(sec2)
          extr2.basePhi = extrLs[i].basePhi; extr2.baseTheta = extrLs[i].baseTheta
          extr2.phi = extrLs[i].phi - get_phi(); extr2.theta = extrLs[i].theta
          extr2.dist = extrLs[i].dist
          extrLs.append(extr2)
          



  # initialization of algorithm
  def initGrow(mid, noDend):
          
          r = [ ranstream(mid, stream_dend), ranstream(mid, stream_apic), ranstream(mid, stream_tuft) ]
          extrLs = [ ]
          nrn = Neuron()

          glomPos = params.glomRealCoords[int(mid/params.Nmitral_per_glom)]
          somaPos = mk_soma(mid, nrn)
          
          mkApical(r[Extreme.APICAL], somaPos, glomPos, nrn, extrLs)

          # make dendrites.
          if not noDend:
                  # get basic inclination
                  somaLvl = Ellipsoid(bulbCenter, somaAxis[0])
                  h, phi_base, theta_base = somaLvl.toElliptical(somaPos)
                  theta_base = pi / 2 - theta_base; phi_base = pi / 2 - phi_base
                  _r = r[Extreme.DENDRITE]
                  
                  ###
                  DENDRITES = int(_r.negexp(N_MEAN_DEND - N_MIN_DEND)) + N_MIN_DEND
                  while DENDRITES > N_MAX_DEND: DENDRITES = int(_r.negexp(N_MEAN_DEND - N_MIN_DEND)) + N_MIN_DEND
                  
                  nmg=params.Nmitral_per_glom
                  dphi=2*pi/DENDRITES
                  phi_phase=_r.uniform(0,2*pi)
                  for i in range(DENDRITES):
                          phiPhase = i*dphi + phi_phase #2*pi/nmg*(mid%nmg)
                          
                          sec, extr = mkDendrite(_r, 0, nrn.soma)
                          p = copy(somaPos); p.append(gen_dend_diam(0.))
                          sec.points = [ p ]
                          sec.index = i              
                          nrn.dend.append(sec)
                          
                          ## initial direction
                          phi = phiPhase
                          
                          theta = pi / 3.
                          newPhi, newTheta = convert_direction(phi, theta, phi_base, theta_base)
                          extr.phi = phi; extr.theta = theta
                          extr.basePhi = phi_base; extr.baseTheta = theta_base
                          extrLs.append(extr)
          return r, nrn, extrLs, glomPos

  def Grow(mid, noDend):
          r = nrn = extrLs = glomPos = None
          r, nrn, extrLs, glomPos = initGrow(mid, noDend) # initialization        
          
          it = 0
          while it < GROW_MAX_ITERATIONS and len(extrLs) > 0:
                  i = 0
                  while i < len(extrLs):
                          j = 0
                          while j < GROW_MAX_ATTEMPTS:
                                  p, phi, theta = genWalk(r, extrLs[i], glomPos)
                                  
                                  if feasible(p, extrLs[i], nrn, glomPos):
                                          updateExtr(extrLs[i], p, phi, theta)
                                          break
                                  j += 1
                                           
                          if canDelete(extrLs[i], j >= GROW_MAX_ATTEMPTS):
                                  if len(extrLs[i].sec.points) <= 1 and extrLs[i].extr_type == Extreme.DENDRITE:
                                          if extrLs[i].sec.index < (len(nrn.dend) - 1):
                                                  for q in range(extrLs[i].sec.index + 1, len(nrn.dend)):
                                                          nrn.dend[q].index -= 1
                                          del nrn.dend[extrLs[i].sec.index] 
                                  
                                  del extrLs[i]        
                          elif canChangeTypeOfExtr(extrLs[i], glomPos):
                                  changeTypeExtr(r, extrLs, i, nrn, glomPos)
                                  del extrLs[i]
                          elif canBifurcate(extrLs[i]):
                                 Bifurcate(r, extrLs, i, nrn)
                                 del extrLs[i]
                          else:
                                  i += 1
                  it += 1
                  
          return nrn

  def genMitral(mid): return Grow(mid, False)
  def genSomaApicalTuft(mid): return Grow(mid, True)


Loading data, please wait...