Dendritica (Vetter et al 2001)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:7907
Dendritica is a collection of programs for relating dendritic geometry and signal propagation. The programs are based on those used for the simulations described in: Vetter, P., Roth, A. & Hausser, M. (2001) For reprint requests and additional information please contact Dr. M. Hausser, email address: m.hausser@ucl.ac.uk
Reference:
1 . Vetter P, Roth A, Häusser M (2001) Propagation of action potentials in dendrites depends on dendritic morphology. J Neurophysiol 85:926-37 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Dendrite;
Brain Region(s)/Organism:
Cell Type(s): Hippocampus CA1 pyramidal cell; Hippocampus CA3 pyramidal cell; Neocortex U1 L2/6 pyramidal intratelencephalic cell; Cerebellum Purkinje cell;
Channel(s): I Na,t; I L high threshold; I p,q; I K; I M; I K,Ca;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: NEURON;
Model Concept(s): Action Potential Initiation; Bursting; Active Dendrites; Influence of Dendritic Geometry; Detailed Neuronal Models; Axonal Action Potentials; Action Potentials;
Implementer(s): Hausser, M [M.Hausser at ucl.ac.uk];
Search NeuronDB for information about:  Hippocampus CA1 pyramidal cell; Hippocampus CA3 pyramidal cell; Cerebellum Purkinje cell; Neocortex U1 L2/6 pyramidal intratelencephalic cell; I Na,t; I L high threshold; I p,q; I K; I M; I K,Ca;
Files displayed below are from the implementation
/* Relating dendritic geometry and signal propagation   */
/* Philipp Vetter, Arnd Roth and Michael Hausser        */
/* geometry.hoc Version 1.0                  11.10.1999 */

func fdistance() {	local frc,n

		if (!electrotonicL) 	{
				if ($1==0)    return distance(0)
				if ($1==1)    return distance(1)
				frc = (.5+int($1*nseg))/nseg    // midpoint of seg 

				
				return distance(frc)+($1-frc)*L*sign_pk
					}

		frc = (nseg*$1)-int(mx($1*nseg-1e-6,0))
		if (!sign_pk) frc = 1-frc

		return Xo_pk + Xfrc_pk($1) + frc*Xsec_pk($1)
		}

/**************************************************************************/

func fL()	 {  if (!electrotonicL) return L return Xlen_pk }

/**************************************************************************/

func segL() 	{ if (!electrotonicL) return L/nseg  return Xsec_pk($1) }

/**************************************************************************/

func mindist() {local i,min

		if (electrotonicL) return Xo_pk  

		ifsec "soma"  { // soma is always tricky!!
			min=fdistance(1)
			for i = 0,nseg-1 { if (min>fdistance(i/nseg)) min=fdistance(i/nseg) }
			return min
				}
		 
		return mn(distance(0),distance(1))
	       }

/**************************************************************************/

func maxdist() {  local i,max
		
		  ifsec "soma"  { 		// if origin is moved, then soma maxdist is non-trivial
			max=fdistance(1)
			for i = 0,nseg-1 { if (max<fdistance(i/nseg)) max=fdistance(i/nseg) }
			return max
				}

		return mx(fdistance(0),fdistance(1))
		 }
 
/**************************************************************************/

func farea() { local I
		/* return membrane area within distance $1 from soma */
		/* this one is used for the geometry calculations   */
		/* includes spinescale corrections                  */
 
		I = 0
		forsec all  		I+= sectionarea($1)*spine_pk
		return I
	   }

/**************************************************************************/

func sectionarea() { local areavar,i,c,min,max,I
	     /* return membrane area of a given section */ 
	     /* within distance $1 from origin(originx) */

	  	if(mindist()>$1)     		            return 0           /* completely out */
	  	areavar = 0
	  	if(maxdist()<$1){ for(x) areavar +=area(x)  return areavar   } /* completely in */
										/* in between */
		if (nseg==1) 				    return area(.5)*($1-mindist())/fL()

		for i=0,nseg-1 {
				c=(i+.5)/nseg
				min = fdistance(c)-segL(c)/2
				max = min+segL(c)
				if (max<=$1) { areavar += area(c) } else {
				if (min<=$1) { areavar += area(c)*($1-min)/segL(c) }}

				}
							    return areavar
	         }

/**************************************************************************/

func fseg() {   local max, min /* calculate sum of diams of sections at dist $1 from soma */

	diameterlist = new Vector()

	forsec all {if ((maxdist()>=$1)&&(mindist()<$1)) {
			diameterlist.append(diam(($1-mindist())/fL()))
							 }
		   }

	if (diameterlist.size==0) { print "patch for distance ", $1 
				    return no_of_sections }

	no_of_sections     	= diameterlist.size()
	sum_of_diameters	= diameterlist.sum()
	mean_diameter		= diameterlist.mean()
	equiv_crossdiameter	= (diameterlist.sumsq)^(0.5)
	diameterlist.pow(3/2)
	equiv_ralldiameter      = diameterlist.sum()^(2/3)
	return no_of_sections

       	    }

/**************************************************************************/

func fbranch() { 	local i, dmax
			/* return branchnumber for section
			   1 = 2 daughter branches
			   0 = 1 daugheter branch
			  -1 = termination */

			branch_int 	= $2
			branch 		= 0
			termination 	= 0

			origin.sec distance(0,originx)

			forsec all {	if (maxdist()<($1+branch_int) && maxdist()>=$1)  { 
						i = 0
						subtree_traverse("{i+=1 }")
						if (i >= 2) branch += i-1
						if (i == 0) termination +=1  	   }		    
			              }			
			return branch
		}

/**************************************************************************/

proc get_parent() { /* make Parent a SectionRef to the parent section */

		root = 0
		if (parent_section()) {
					push_section(parent_section())
					this_section() Parent = new SectionRef()
					pop_section()
					// Parent.sec print secname()

					} else { 
					this_section() Parent = new SectionRef()
					root = 1 
					}
		}

/**************************************************************************/

func pbranchpoint() { /* if section is a branchpoint, return 1 */
			j = 0
		        Parent.sec subtree_traverse("j+=1")
			branchnumber = j - 1
			if (j>=2) return 1 
			if (root==1) return 1
			return 0			
		   }

/**************************************************************************/

func nextparent() { Parent.sec get_parent()
		    return pbranchpoint() }

/**************************************************************************/

func branchpoint() {	/* if section is a branchpoint, return 1 */
			j = 0
		        subtree_traverse("j+=1")
			branchnumber = j - 1
			if (j>=2) return 1
			return 0			
		   }

/**************************************************************************/

proc get_root()		{ /* put root into Parent  */

			  Soma.sec get_parent()
			  Parent.sec	{ while (!root) { get_parent() } } 				
			}

/**************************************************************************/

proc ubranch() {local a 
		this_section() Parent = new SectionRef()
		a = 0
		while (a==0) {  Parent.sec connection = parent_connection()
				a = nextparent()
				Parent.sec pdistance = fdistance(connection)
			     }
		udistance = abs(mindist() - pdistance)
		}

/**************************************************************************/

proc get_rall() {	local i, k /* calculate the relative diameters of the parent to branch */

			dchild.resize(0)	
			dparent.resize(0)	
			dY.resize(0)		
			ch.resize(0)		

			branchpointd.resize(0)		
			relaread.resize(0)		
			ralld.resize(0)		


			/* create 3 Vectors [ dchild dparent dY ] */
			get_children()
			forsec children {  Y               = parent_connection()
					   X               = section_orientation()
					   X2              = X + diam_int/L - (X)*2*diam_int/L
					   X2              = mn(X2,1)
					   X2          	   = mx(X2,0)
					   diameter_child  = (diam(X)+diam(X2))/2
			       Parent.sec  diameter_parent = diam(Y)
					   dY.append(Y)
					   dchild.append(diameter_child)
					   dparent.append(diameter_parent)
					  }
			
			while (dY.size() >= 1) {

				ref = dY.indvwhere("==",dY.x[0])

				if (dY.x[0] != 1)    ch.append(dparent.x[0])

				if (ref.size()==1) { ch.append(dchild.x[0])
						     if (dY.x[0] != 1) rall_calc(dparent.x[0],ch) }

				if (ref.size()>1)  { for i = 0,ref.size()-1 {ch.append(dchild.x[ref.x[i]])}
						     rall_calc(dparent.x[0],ch)}
									
				ref  	= dY.indvwhere("!=",dY.x[0]) /* clean up */
				dY   	= dY.ind(ref)
				dchild  = dchild.ind(ref)
				dparent = dparent.ind(ref)
					     }
		}

/**************************************************************************/

proc rall_calc() { local rall_children /* input parent, children */

			parent = $1
			rall_children = 0
			for i = 0,$o2.size()-1 { rall_children += $o2.x[i]^(3/2) }
			
			branchpointd.append($o2.sum()/parent)
			relaread.append( $o2.sumsq()/(parent^2))
			ralld.append(rall_children/(parent^(3/2)))
		}

/**************************************************************************/

proc ename() { 	
		str1 = $s1
		if (electrotonicL) sprint(str1,"a%s",str1)
		sprint(str1,"%s = %s",str1,$s2)
		execute(str1)

		}
		
/**************************************************************************/

func gstep()		{ 
	/* returns the increment that is used for either el or physical lengths */

	if (!electrotonicL) {

if (strob.substr(cellname,"Purkinje")>-1) {stepdist=purkinjestep}else{stepdist=1}
	 if (synthetic) stepdist = 5
	
	 return stepdist
			    }
	stepdist = 0.004
	if (synthetic) stepdist = 0.02

	return stepdist
			}

/**************************************************************************/

proc make_dAr()         { local i
		dAr.resize(Ar.size())
		dAr.indgen(gstep(),gstep())
		dAr.x[0] = Ar.x[0]
	    	for i = 1,dAr.size()-1  { dAr.x[i] = (Ar.x[i]-Ar.x[i-1])/gstep() }
			}

/**************************************************************************/

proc get_gdist()	{ /* sets size of gdist, dAr & Ar */

		L_switch() /* make sure to get adist.max */
		if (numarg()>0) L_switch($1)
		{gdist.indgen(0,adist.max,gstep())   gdist.x[0] = gstep()/40 }
		if (Ar.size>1) make_dAr()
			}

/**************************************************************************/

proc geometry_calc()	{ local a, i, stemdiam, n, Distance, Diameter, bp, j
		/* calculate all values necessary for geometry_save */
		/* writes to vector "areas" */
		/* for branchpoints finds distance to soma [SDist] */
		/* upstream branch distance [DDist] */
		/* parent diameter [Diam_abs] */
		/* relative branch areas [Diam_relarea] */
 		/* rall ratio [Diam_rall] */

		n = numarg()

if (!equiv) {	

	/*********************/
	/* Rall values 	     */
	/*********************/

	somadist.resize(0)
	updst.resize(0)
	branchpointdiam.resize(0)
	ralldiam.resize(0)
	relareadiam.resize(0)
	ralldiam.resize(0)
	somadist_noend.resize(0)
 	updst_noend.resize(0)
	branchpointdiam_noend.resize(0)
	ralldiam_noend.resize(0)
	relareadiam_noend.resize(0)
	ralldiam_noend.resize(0)
	terminationdist.resize(0)

	forsec branchpoints {		ubranch()
					updst.append(udistance)
		if (!isterminal())	updst_noend.append(udistance)

					somadist.append(fdistance(0))
		if (!isterminal())	somadist_noend.append(fdistance(0))

					get_rall()
					branchpointdiam.append(branchpointd)
					relareadiam.append(relaread)
					ralldiam.append(ralld)
		if (!isterminal()){	branchpointdiam_noend.append(branchpointd)
					relareadiam_noend.append(relaread)
					ralldiam_noend.append(ralld)	}}

	forsec terminations   terminationdist.append(fdistance(1))

	branchpoints_num       = ralldiam.size()

	ename("branchdensity","mean(updst)")  
	ename("branchdensity_noend","mean(updst_noend)")

	length  = 0
	forsec all length += fL() 
	ename("branchdensityII","div(ralldiam.size(),length)")

	length  = 0
	forsec all_noend length += fL()	
	ename("branchdensityII_noend","div(ralldiam_noend.size(),length)")

	rallratio_peak		= 0
	rallratio_mean		= 0
	diamratio_peak		= 0
	diamratio_mean 		= 0
	rallratio_noend_peak	= 0
	rallratio_noend_mean	= 0
	diamratio_noend_peak	= 0
	diamratio_noend_mean	= 0

	if (ralldiam.size>0) {
				  gauss(ralldiam,0.01,0.01) /* calculate ralldiam "peak" */
	rallratio_peak		= gaussx.x[gaussy.indwhere("==",gaussy.max)]
	rallratio_mean		= ralldiam.mean()
				  gauss(branchpointdiam,0.01,0.01) 
	diamratio_peak		= gaussx.x[gaussy.indwhere("==",gaussy.max)]
	diamratio_mean 		= branchpointd.mean
			}
	

	if (ralldiam_noend.size>0) {
				  gauss(ralldiam_noend,0.01,0.01) 
	rallratio_noend_peak	= gaussx.x[gaussy.indwhere("==",gaussy.max)]
	rallratio_noend_mean	= ralldiam_noend.mean()
				  gauss(branchpointdiam_noend,0.01,0.01)
	diamratio_noend_peak	= gaussx.x[gaussy.indwhere("==",gaussy.max)]
	diamratio_noend_mean	= branchpointdiam_noend.mean 
			}




	/********************************************************************/
	/* calculate mean stem dendrite diameter (ref.mean) */
	/********************************************************************/

	ref = new Vector()
	forsec "soma" { 
			get_children()
			forsec children { if (!issection("som*")) { if (!issection("hill*")) {
						stemdiam = mn(1,5/L) 
						ref.append(diam(stemdiam)) }}}}

	ename("mean_stem_dendrite_diam","mean(ref)")


	/********************************************************************/
	/* taper = tap_ddiam.sum()/tap_ddist.sum()                          */
	/********************************************************************/

	tap_dist  = new Vector()
	tap_ddist = new Vector()
	tap_ddiam = new Vector()
	forsec all {
 		 get_parent()
		 Parent.sec bp = branchpoint()

		if (!bp) {      i = .5/nseg
			 	Distance = segL(i)
				Diameter = diam(i)

		      		Parent.sec	Distance += segL(i)
				Parent.sec 	Diameter -= diam(1-i)

				tap_ddist.append(Distance) 
				tap_dist.append(fdistance(i))
				tap_ddiam.append(Diameter)
				}

		if (nseg > 1) {

		for i = 2,nseg { j = (i-.5)/nseg
				 tap_ddist.append(segL(j))
				 tap_dist.append(fdistance(j))
				 tap_ddiam.append(diam(j) - diam(j-(1/nseg)))}
		}}

	ename("taper","div(tap_ddiam.sum,tap_ddist.sum)")

	/********************************************************************/

	}

	if (n==0) {

	/***************************************/
	/* calculate area as a fxn of distance */
	/***************************************/

	{gdist.indgen(0,adist.max,gstep())   gdist.x[0] = 0.0001 }
	Ar = gdist.c
   	Ar.apply("farea")		
	make_dAr()
	if (!electrotonicL) temp = 25 else temp = 0.12

	area_max	        = Ar.max
	val = mn(temp/stepdist,dAr.size()-1) // patch to make sure no false copy
	mdarea_max = dAr.c(val,dAr.size()-1).max
	mdarea_maxdist = dAr.indwhere("==",mdarea_max)*gstep()

	ename("darea_max","mdarea_max")
	ename("darea_maxdist","mdarea_maxdist") 

	writevec_el("areas",Ar)

		}	

}


/* functions to avoid NaNs */

func mean() { if ($o1.size()>0) return $o1.mean() else return 9999 }
func div()  { if ($2!=0)        return $1/$2      else return 9999 }



/**************************************************************************/
/**************************************************************************/
/**************************************************************************/
/**************************************************************************/
/**************************************************************************/
/**************************************************************************/
/**************************************************************************/
/**************************************************************************/
/**************************************************************************/
/**************************************************************************/
/**************************************************************************/
/**************************************************************************/
/**************************************************************************/
/**************************************************************************/

/**************************************************************************/
/* equivalent cylinder calculations                                       */
/**************************************************************************/

func equivalent_calc() { local i,temp, F

	/* geometric transformation of cell to incorporate spines */

	if (equiv) { print "equivalent_calc() cannot be applied to equivalent cables"
			   return 1 }

	//forsec "soma"  		
	forsec dendrites	spinetransform(spine_pk)

	print "geometry changed; reload cell to continue working"

	/* mean diameter [mdiam] and # sections as fxn(r) */
	/**************************************************/

   	origin.sec distance(0,originx)
	sections.resize(gdist.size())
	mdiam.resize(gdist.size())
	cdiam.resize(gdist.size())
	rdiam.resize(gdist.size())
	for i=0,gdist.size()-1 {     
			fseg(gdist.x[i])
			sections.x[i]   = no_of_sections
			mdiam.x[i]  	= mean_diameter
			cdiam.x[i] 	= equiv_crossdiameter
			rdiam.x[i]  	= equiv_ralldiameter  }	




	if (!electrotonicL) { temp = 25 } else {temp=0.12}
	val = mn(temp/stepdist,rdiam.size()-1) // patch to make sure no false copy
	mdeq_max 	  = rdiam.c(val,rdiam.size()-1).max
	ename("deq_max","mdeq_max")
	ename("deq_maxdist","rdiam.indwhere(\"==\",mdeq_max)*gstep()")
	ename("sections_mean","mean(sections)")
	ename("sections_max","sections.max")
	ename("sections_maxdist","sections.indwhere(\"==\",sections.max)*gstep()")
	ename("diam_mean","mean(mdiam)")
	ename("taper_mean","lintaper()")  // uses mdiam

	writevec_el("cdiam",cdiam)
	writevec_el("rdiam",rdiam)
	writevec_el("mdiam",mdiam)
	writevec_el("sections",sections)

	return 1
			}

/**************************************************************************/
double tmpx[1000]
double tmpy[1000]
double tmpz[1000]
double tmpdiam[1000]


proc spinetransform() { local tmpn,i 
			tmpn = n3d()
    			F = $1
    			L *= F^(2/3)
   
			for i=0,tmpn-1 {
					tmpx[i] = x3d(i)
					tmpy[i] = y3d(i)
					tmpz[i] = z3d(i)
					tmpdiam[i] = diam3d(i)*F^(1/3)
					}

			pt3dclear()
 			for i=0,tmpn-1 { pt3dadd(tmpx[i], tmpy[i], tmpz[i], tmpdiam[i])  }
			define_shape()
  			}

/**************************************************************************/

proc make_equivalent_cable() { local i,j

	/* need gstep & rdiam in the electrotonic domain */
	/* create soma, dend with physical lengths */
	
	L_switch(1)
	dX = gstep()

	sprint(celladdress,"%sequiv",celladdress)	

	writeveca("equiv")
	fprint("create soma\n")	
	fprint("create dend[%d]\n\n",rdiam.size()-1)	

	fprint("soma.diam = %2.4f\n",  mx(rdiam.x[0],0.0001))
	fprint("soma.L    = %2.8f\n\n",mx(100*dX*sqrt(0.25*Rm*rdiam.x[0]/Rax),0.00001))

	fprint("dend[0].diam = %2.4f\n",  rdiam.x[1])
	fprint("dend[0].L    = %2.4f\n",  100*dX*sqrt(0.25*Rm*rdiam.x[1]/Rax))
	fprint("soma connect dend[0](0), 1\n\n")


	for i = 2,rdiam.size()-1 { j = i-1

		fprint("dend[%d].diam = %2.4f\n",  j, rdiam.x[i])
		fprint("dend[%d].L    = %2.4f\n",  j, 100*dX*sqrt(0.25*Rm*rdiam.x[i]/Rax))
		fprint("dend[%d] connect dend[%d](0), 1\n\n",j-1,j)

				}
	wopen()
	cell_name(cell)  /* restore names */

}



/**************************************************************************/
/* these procedures need to be checked by eye for outliers                */
/**************************************************************************/


proc slope_darea() {	local n	/* batch procedure */
			/* writes to "d2area" */
			/* calculates d2area_max,d2area_maxdist,*/
			/*d2Area_maxAr_ratio,d2area_maxAr_percent */ 
			/* $s1 cell */
			/* $s2 active model */
			/* n>2 print figure at distance $3 */
			/* n=4 don't write results */
			/*   e.g. slope_darea(int2,act0) */

			n = numarg()
			if (n>0) cell_name($s1)
			if (n>1) ActiveModel=$s2
			activecell = cellname
			dAr   = new Vector()
			Ar    = new Vector()
			gdist = new Vector()
			readvec("areas",Ar)
			electrotonicL=0
			make_dAr()
			{gdist.indgen(0,(dAr.size+1)*gstep(),gstep())   gdist.x[0] = 0.0001 }

			if (n<3) rolling(gdist,dAr,63)
			if (n==3) rolling(gdist,dAr,63,"print figure",$3)
			if (n==4) rolling(gdist,dAr,63,"print figure")

			fprint("d2area_max           = %2.4f\n", d2area_max)
			fprint("d2area_maxdist       = %2.4f\n", d2area_maxdist)
			fprint("d2area_maxAr_ratio   = %2.4f\n", d2area_maxAr_ratio)
			fprint("d2area_maxAr_percent = %2.4f\n", d2area_maxAr_percent)

			if (n<4) writevec("darea")
			fprint("d2area_max           = %2.4f\n", d2area_max)
			fprint("d2area_maxdist       = %2.4f\n", d2area_maxdist)
			fprint("d2area_maxAr_ratio   = %2.4f\n", d2area_maxAr_ratio)
			fprint("d2area_maxAr_percent = %2.4f\n", d2area_maxAr_percent)
			if (n==3) fprint("// modified manually\n")
			if (n<4) wopen()
		}

/**************************************************************************/

proc slope_deq() {	local n
			/* calculates ddeq_max,ddeq_maxdist,*/
			/* ddeq_maxAr_ratio,ddeq_maxAr_percent */ 
			/*   $s1 cell */
			/*   $s2 active model */
			/*   $3 manually set maximum in x (optional) */
			/*   $4 do not set slope - just check to see */
			/*   e.g. slope_deq(int2,act0) */
			n = numarg()
			if (n>0) cell_name($s1)
			if (n>1) ActiveModel=$s2
			activecell = cellname
			rdiam = new Vector()
			gdist = new Vector()
			readveca("rdiam",rdiam)
			electrotonicL=1
			{gdist.indgen(0,(rdiam.size+1)*gstep(),gstep())   gdist.x[0] = 0.0001 }
			if (n<3)  rolling2(gdist,rdiam,63)
			if (n==3)  rolling2(gdist,rdiam,63,"print figure",$3)
			if (n==4)  rolling2(gdist,rdiam,63,"print figure")

			fprint("ddeq_max           = %2.4f\n", d2area_max*100)
			fprint("ddeq_maxdist       = %2.4f\n", d2area_maxdist)
			fprint("ddeq_maxAr_ratio   = %2.4f\n", d2area_maxAr_ratio)
			fprint("ddeq_maxAr_percent = %2.4f\n", d2area_maxAr_percent)

			if (n<4)  writeveca("ddeq")
			fprint("ddeq_max           = %2.4f\n", d2area_max*100)
			fprint("ddeq_maxdist       = %2.4f\n", d2area_maxdist)
			fprint("ddeq_maxAr_ratio   = %2.4f\n", d2area_maxAr_ratio)
			fprint("ddeq_maxAr_percent = %2.4f\n", d2area_maxAr_percent)
			if (n==3) fprint("// modified manually\n")
			if (n<4) wopen()
		}



/**************************************************************************/
/* procedures requiring manual supervision */
/**************************************************************************/


proc dAdr_calc() { local n
/* calculates dAdr (has to be checked manually) */
/* written into "dAdr" */
/* $s1 specifies cell */
/* n==1 calculate and save dAdr */
/* n==2 calculate and display dAdr (don't save) */
/* n==3 calculate and display dAdr; save max ($2) min ($3) */
/* n==4 calculate and display and save dAdr */

	 	  n = numarg()		
		  cell_name($s1)
		  electrotonicL=0
		  Ar    =new Vector()
		  dAr   =new Vector()
		  gdist =new Vector()
		  readvec("areas",Ar)	
		  gdist.resize(Ar.size)
		  { gdist.indgen(0,gstep()) gdist.x[0]=0.0001 }
		  make_dAr()

		  if (n==1) estcore(dAr,1) else  estcore(dAr)
		  if (n==1) dAdr_write(max,min)
		  if (n==3) dAdr_write($2,$3)
		  if (n==4) dAdr_write(max,min)

		}

/**************************************************************************/

proc dAdr_write() {     local n /* writes to "dAdr" */
			n=numarg()
			if (n>3) ActiveModel=$s4
			if (n>2) cell_name($s3)
			if ($1 < $2) print "please type values in again\n"
			fprint("dAdr_relmax	= %2.4f\n",$1)			
			fprint("dAdr_relmin	= %2.4f\n",$2)			
			fprint("dAdr_ratio	= %2.4f\n",$1/$2)			
			writevec("dAdr")
			fprint("dAdr_relmax	= %2.4f\n",$1)			
			fprint("dAdr_relmin	= %2.4f\n",$2)			
			fprint("dAdr_ratio	= %2.4f\n",$1/$2)			
			wopen()
		}

/**************************************************************************/

proc deq_calc() { local n
		/* calculates deq (has to be checked manually) */
		/* written into "deq" */
		/* $s1 specifies cell */
		/* n==1 calculate and save deq */
		/* n==2 calculate and display deq (don't save) */
		/* n==3 calculate and display deq; save max ($2) min ($3) */
		/* n==4 calculate and display and save deq */

		n = numarg()
		  cell_name($s1)

		  electrotonicL=1
		  rdiam=new Vector()
		  gdist=new Vector()
		  readveca("rdiam",rdiam)
		  gdist.resize(rdiam.size)
		  { gdist.indgen(0,gstep()) gdist.x[0]=0.0001 }

		  if (n==1) estcore(rdiam,1) else estcore(rdiam)
		  if (n==1) deq_write(max,min)
		  if (n==3) deq_write($2,$3)
		  if (n==4) deq_write(max,min)
		}


proc deq_write() {      local n /* writes to "deq" */
			n=numarg()
			if (n>3) ActiveModel=$s4
			if (n>2) cell_name($s3)
			if ($1 < $2) print "please type values in again\n"
			fprint("deq_relmax	= %2.4f\n",$1)			
			fprint("deq_relmin	= %2.4f\n",$2)			
			fprint("deq_ratio	= %2.4f\n",$1/$2)			
			writeveca("deq")
			fprint("deq_relmax	= %2.4f\n",$1)			
			fprint("deq_relmin	= %2.4f\n",$2)			
			fprint("deq_ratio	= %2.4f\n",$1/$2)			
			wopen()
		}




/**************************************************************************/

proc estcore() { local n, start
		n = numarg()

		if (electrotonicL)  start = 0.12/gstep() else start = 25/stepdist
		max =$o1.c(start, $o1.size()-1).max
		min =$o1.c(start,mx($o1.indwhere("==",max),start)).min

		maxx = gdist.x[$o1.indwhere("==",max)]
		minx = gdist.x[$o1.indwhere("==",min)]

		ref = new Vector(2,min)
		ref2 = new Vector(2,minx)
		ref.x[1] = max
		ref2.x[1] = maxx

    if (n<2) {  fig(gdist,$o1)
		figlab(cellname,0.2,0.9)
		fig(ref2,ref,1,1,2,8,"o")
		figure[figcurrent].unmap
	        figure[figcurrent].view(0,0,gdist.max,$o1.max,0,100,1000,500) }
		print "min = ",min, "\tmax = ", max
		}

/**************************************************************************/

func tap()     { /* ugly function for section taper */

		tobj 	= gdist.c
		tobj.mul(0.001*$&2[0]).add($&2[1])

		return   mdiam.meansqerr(tobj)
		}

/**************************************************************************/

func lintaper(){ /* ugly linear regression for section taper */

		if (!electrotonicL) param[0] = -1 else param[0] = -1000
		param[1] = 10

                attr_praxis(1e-5,4,0)
                fit_praxis(2,"tap",&param[0])
                return param[0] // slope in 1/mm 
                }









func get_link()		{ local i,ok
		/* find sections between given section G and origin O [exclusive] */
		/* sign_pk = 1 if G is downstream of O or neutral wrt O */
		/* returns link, sign_pk and origin_c */

		link=new SectionList()
		sign_pk=1    
		origin_c   = parent_connection()
 
		ifsec originsec  { if (originx>0.5) sign_pk=-1 return 1	}
		get_parent()
		ok=0
		Parent.sec 	ifsec originsec  ok=1
		if (ok) return 1

		// search append G -> O

		while (!root) {	Parent.sec link.append()
				Parent.sec origin_c = parent_connection() 
				Parent.sec get_parent()
				Parent.sec ifsec originsec  ok = 1
				if (ok) return 1
			        }

			//  O is not root
			//  O is not upstream of G
			//  O  may be downstram of G
			//  put link > link2
			//  search O -> G 
			// Hypothesis: G upstream of O 

			sign_pk=-1 	
			link2 = new SectionList()
			forsec link { link2.append() link.remove() } 
			origin.sec   origin_c = parent_connection()
		    	origin.sec   get_parent()
			str1 = secname()

		    	while (!root) { Parent.sec link.append()
					Parent.sec get_parent()
					Parent.sec ifsec str1 	ok = 1
					if (ok) return 1  
				      }
			
			// G is neutral with respect to O
			// combine  link + link2 */
			// all sections which appear in both sectionlists must disappear
			// combine sections which appear in only one of the two sectionlists

			sign_pk=1
			forsec link2 {	ifsec link { link.remove() link2.remove()}}
			forsec link2 {  link.append() }
			forsec link2 link2.remove()
			link2.append()
			forsec link2 { ifsec link { link.remove() }}
			
			return 1
			}


proc set_electrotonic() { local Lx, i, lambdax, f, oseg, Xcum, toggle, c, d, xo

		/* set the electrotonic distances for all segments */
		/* set Xsec for each segment */
		/* Xfrc_pk is fraction of section Length before given segment */
		/* Xo_pk is base length */

		if (!electrotonicL) toggle=1 else toggle=0
		L_switch(1)

		forsec all { Lx      = L/nseg 	
		      	 Xlen_pk = 0
			 for i = 1,nseg {  f 		= (i-.5)/nseg
					  lambdax 	= sqrt(diam(f)/(4e-4*g_pas*Ra))
					  Xsec_pk(f)  	= Lx/lambdax 
					  Xlen_pk      += Xsec_pk(f)
					}

			Xcum    = 0
			get_link()  // get sign_pk
		        for i = 1,nseg {  f 		= (i-.5)/nseg
					  Xfrc_pk(f)    = Xcum
					  Xcum         += Xsec_pk(f)
					  if (!sign_pk) Xfrc_pk(f) = Xlen_pk-Xcum //!!!
				       }
			 }

		if (originx!=1) {  /* Xfrc is slightly different here */
			
		access origin.sec   

		oseg = (nseg*originx)-int(mx(originx*nseg-.0001,0))+1  // which seg is origin
		for(x) Xfrc_pk(x) = Xsec_pk((oseg-.5)/nseg)
		Xfrc_pk((oseg-.5)/nseg) = 0   // sufficient until nseg >3

		if (nseg>3) {	

				for i=1,nseg { 	

		if (i<oseg-1) for j=i+1,oseg-1 { Xfrc_pk((i-.5)/nseg) += Xsec_pk((j-.5)/nseg) }
		if (i>oseg+1) for j=oseg+1,i-1 { Xfrc_pk((i-.5)/nseg) += Xsec_pk((j-.5)/nseg) }
						}			
			     }			
				}




		forsec all  Xo_pk = 0
		forsec all  {

  		get_link()    
		origin.sec.Xo_pk = 0
		origin.sec   xo=fdistance(origin_c) // how much of origin section to be added
		forsec link  xo += Xlen_pk 
		Xo_pk = xo
		origin.sec.Xo_pk = 0

			    } 

		if (toggle)  L_switch(0) 
		}

Loading data, please wait...