# -*- 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)