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        */
/* electrophysiology.hoc Version 1.0         11.10.1999 */

proc rest() { 	local i,n /* set neuron to steady-state resting potential and t = 0 */

		n=numarg()
		inactive vclamp.loc(0)  // take out electrode until rest is achieved
		inactive st.loc(0)
		inactive synapse.loc(0)


		finitialize(-70.494513)  // actual resting - already perfect setting for act0
		t   = - Rest_dur
		dt  = Rest_dur/10
       	        while (t<0) fadvance()
		finitialize()


		if (!n) { // usually put point processes back in
			Iclamp.sec st.loc(Iclampx) 
			Vclamp.sec vclamp.loc(Vclampx) 
			Synapse.sec synapse.loc(Synapsex) 
			}

		if (movie) { doEvents()  splot.flush() }
		}

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

proc simcore()  {/* simulation-core (equivalent to run) */

		rest("leave point processes off")

		if (simMode==0)	Iclamp.sec st.loc(Iclampx) 
		if (simMode==1) Vclamp.sec vclamp.loc(Vclampx) 
		if (simMode==5) Synapse.sec synapse.loc(Synapsex) 


		dt = sim_dtI
		if (sim_durI>0) while (t<sim_durI) { fadvance() }
		dt = sim_dt
		while (t<(sim_durI+sim_dur))  { fadvance()
						if (movie) { doEvents() splot.flush() } 
					      }

		forall for(x)   { // 20.8.98
		  vhalf_pk(x)     = vpeak_pk(x)/2 + vonset_pk(x)/2
		  onset_ref_pk(x) = dvdtpeak_pk(x)*sqrt(abs(1.44/(vpeak_pk(x)+73)-1/93)) 
		  if (forward)	        vonset_pk(x) = vrest_pk(x)
		  if (diffvonset)	vonset_pk(x) = vrest_pk(x)
		  // 8.4.99 after discussion with MH AR PV
				}
		  }

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

proc sim()	{ local n, st_intensitysave, vclamp_dursave,mn,mx
	/* this is the core simulation procedure */
	/* composed of Rest and then actual simulation */
	/* specify, whether voltage clamp or current clamp */
	/*****************************************************/
	/* sim runs the simulation */
	/* it is composed of rest() and simcore() */
	/* 					*/
	/* simulations come in four flavours */
	/* simMode 0	IClamp at originx */
	/* simMode 1 	VClamp at originx */
	/* simMode 2	IClamp at originx g_na = 0 */
	/* simMode 3	VClamp at originx g_na = 0 */


	n = numarg()
	if (movie) { initPlot()  }

	// save locations of point processes

	Vclampx    = vclamp.get_loc()
	sectionname(sec)
	sec Vclamp = new SectionRef()  
	pop_section()

	Iclampx    = st.get_loc()
	sectionname(sec)
	sec Iclamp = new SectionRef()  
	pop_section()

	Synapsex    = synapse.get_loc()
	sectionname(sec)
	sec Synapse = new SectionRef()  
	pop_section() 


	rec_set()       // set recording options


	if (simMode==0) {
			/* Current clamp */
			simcore()
			if (!movie)	simcore()
			}

	if (simMode==1) {
			/* Voltage clamp */
			/* play p18 somatic AP waveform */
			/* scale waveform by $1 (optional) */
			/* specify waveform vector with $o2 (optional) */


			if (n==2) set_waveform($o2.c)
			if (n==3) set_waveform($s2,1)
			if (waveform.size()==0) {  set_waveform() }
			wave = waveform.c
			mn = wave.min
			mx = wave.max
			if (n>0) wave.scale(mn,mn+$1*(mx-mn))
			wave.play(&vclamp.amp1)		

			simcore()
			if (!movie) simcore()
			}

	if (simMode==2) {	
			/* Vclamp passive geometry but use IClamp waveform */
			simMode=0
			simcore()     /* get the vector to play */
			simMode=1
			active_set(0,g_kv,0)
			ref2 = vsoma.c
			ref2.play(&vclamp.amp1)		
			simMode=2
			simcore()
			}

	if (simMode==3) {	
			/* do passive AP  with VClamp */
			simMode=1
			active_set(0)

			if (n==2) set_waveform($o2.c)
			if (n==3) set_waveform($s2,1)
			if (waveform.size()==0) {  set_waveform() }
			wave = waveform.c
			mn = wave.min
			mx = wave.max
			if (n>0) wave.scale(mn,mn+$1*(mx-mn))
			wave.play(&vclamp.amp1)		

			simcore()
			if (!movie) simcore()
			active_set()
			simMode=3
			}

	if (simMode==5) {	
			/* do Synapse */
			simcore()
			if (!movie) simcore()
			}

			Iclamp.sec st.loc(Iclampx)
			Vclamp.sec vclamp.loc(Vclampx)
			Synapse.sec synapse.loc(Synapsex)

	}

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








proc initsimvclamp() {
			ActiveModel = $s1
			get(p18.loc,ActiveModel,"no read")
			threshold_find()
			st_set(st_intensity)
			simMode = 0
			sim()
       			sprint(cellvar,"../data/%s/simvclamp_p18",ActiveModel)
			fileob.wopen(cellvar)
			vsoma.c.vwrite(fileob)
			}





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

proc dvdr_calc() {	local i, x, V1
			forall  {get_parent()
				Parent.sec dr  = .5*L/nseg
				Parent.sec V1  =  vpeak_pk((nseg-.5)/nseg) 
					   dr += .5*L/nseg
					    x  = .5/nseg
			dvdr_pk(x) = (vpeak_pk(x)-V1)/(.5*(vpeak_pk(x)+V1) - vonset_pk(x))/dr
					
				if (nseg > 1) { 
						dr = L/nseg
						V1 = vpeak_pk(x) 
						for i = 2,nseg {  x = (i-.5)/nseg 
			dvdr_pk(x) = (vpeak_pk(x)-V1)/(.5*(vpeak_pk(x)+V1) - vonset_pk(x))/dr
	 				        V1 = vpeak_pk(x) }
						}  } 				
		}

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

func sim_err() {
		/* run simulation of active model, calculate error of settings */
		/* -> used for fitting routine	*/

		active_set($&2[0],$&2[1],$&2[2])    
		finitialize()
		active_run()
		return err
		}

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

proc sim_fit() {    
		/* fit active parameters */
		/* $1 gbar_kv/100 */
		/* $2 gbar_na/100 */
		/* $3/10'000 gbar_na (node) */

		param[0] = $1
		param[1] = $2
		param[2] = $3
		act_init()
		attr_praxis(1e-5,5000,0)
		fit_praxis(3,"sim_err",&param[0])	
	          }

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

proc rinput_calc()     { local vx,sx,syx
			/* calculate input resistance */
			/* put result in vector rin  */

	rest("don't put back")

	Soma.sec 	impR = new Impedance()
	Soma.sec	impR.loc(0.5)
	Soma.sec	impR.compute(0)
	Soma.sec	rin.x[0] = impR.input(0.5)
		   print "Rin = ", rin.x[0] 
	rest()

		      }

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

proc sim_calc() {	local i, j, k
			/* calculate values from simulation run     			*/
			/* vectors vpk, amp, vmax, plat, olat, half			*/
			/* site of initiation						*/
			/* err - deviation from experimental data in Pyramidal cells	*/

		err_amp = err_vmax = err_plat = err_olat = initiation = 0

		dvdr_calc()

		i=0
		vpk.resize(0)
		amp.resize(0)
		vmax.resize(0)
		plat.resize(0)
		olat.resize(0)
		half.resize(0)
		dvdr.resize(0)
				
	   	forsec distlist { for j=0,nseg-1 {
						k = (j+.5)/nseg
						vpk.append(vpeak_pk(k))
						amp.append(vpeak_pk(k)-vonset_pk(k))
						vmax.append(dvdtpeak_pk(k))
						plat.append(tpeak_pk(k) - Soma.sec.tpeak_pk)
						olat.append(onset_pk(k)-Soma.sec.onset_pk)
						dvdr.append(dvdr_pk(k))
						half.append(halfwidth_pk(k))
						 } }


		/* calculate measure of error for active parameter settings */

		g = 0 
		forall ifsec "iseg[1]" g = 1
		
		if ((numarg()<1)&&(g)){
		ref = dist.c
		err_amp  += amp.meansqerr(ref.div(-amp_lmd).apply("exp")\
			   .mul(amp_max).add(amp_bas))/amp.max
		ref = dist.c
		err_vmax += vmax.meansqerr(ref.div(-vmax_lmd).apply("exp")\
			   .mul(vmax_amp))/vmax.max
		ref = dist.c
		err_plat += plat.meansqerr(ref.div(plat_s))/1.2
		ref = dist.c
		err_olat += olat.meansqerr(ref.div(olat_s))/0.5

		isg = iseg.tpeak_pk 
		sma = Soma.sec.tpeak_pk 
		nde = node[0].tpeak_pk 
		
		if (sma <  isg && isg <= nde) print i, "initiation soma -> iseg -> node"
		if (sma >  isg && isg >= nde) print i, "initiation node -> iseg -> soma"
		if (sma >  isg && isg <  nde) {
			if (sma < nde)	    { print i, "initiation iseg -> soma -> node"}\
			else	            { print i, "initiation iseg -> node -> soma"} }
		if (sma ==isg && isg <= nde)  print i, "initiation soma == iseg -> node"
		if (sma ==isg && isg > nde)   print i, "initiation node -> iseg == soma"

		if (iseg.tpeak_pk > Soma.sec.tpeak_pk) initiation += 1000000
		err = err_amp*2 + err_vmax + err_plat + err_olat + initiation }
              


		sprint(str1,"g_kv %g  g_na %g  g_na(node) %g  err %g", Soma.sec.gbar_kv,\
								       Soma.sec.gbar_na,\
								       node[0].gbar_na,err)
		if (numarg() < 1) print str1
		
		}







/*****************************************************************/
/* threshold_calc() -> threshold_find() -> threshold() -> thresh */
/*****************************************************************/


proc threshold_calc() {
	/* calculate various thresholds, which quantify the amount */
	/* of backpropagation in the cell. */
	/* - g_na threshold when all sections are depolarized >0 mV */ 
	/* - g_na threshold when all sections are depolarized >0 mV (VClamp p18 somatic AP) */ 
	/* - g_na threshold when terminations are depolarized >0 mV */ 
	/* - g_na threshold when terminations are depolarized >0 mV (VClamp p18 somatic AP) */ 
	/* - AP amplitude 200 um from soma relatic to somatic AP amplitude */ 
	/* - AP amplitude 200 um from soma relatic to somatic AP amplitude (passive cell) */ 
	/* - AP amplitude 200 um dependence on g_na */ 

		rinput_calc()
		input_resistance = rin.x[0]			
		threshold_find(1,0) /* (all>0mV) */
		remove_axon()
		threshold_find(1,1) /* vlcamp p18 somatic AP waveform (all>0mV) */
		threshold_find(3,1) /* vlcamp p18 somatic AP waveform (terminations>0mV)*/
		APdecay_sensitivity()
		APdecay(1) // active 
		APdecay(3) // passive 

			}

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

proc forwardthreshold_calc() {

	/* for forward propagation, only VClamp commands make sense */
	/* just look for depolarizations > injection  at the soma */
	/* it won't work for all sections, since they are far & there is decay */
	/* finds g_na thresh with AP200 waveform (soma) */

if (!equiv) threshold_find(0,1) else {if (diffvonset) threshold_find(5,1) else threshold_find(1,1) }

	APdecay(1)  			// active 
	APdecay(3)  			// passive

	rinput_calc()
	input_resistance = rin.x[0]			

			}

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

proc threshold_find() { local n, value, where, dx, hi, lo
			/* find stimulation threshold to initiate nodal AP while g_na = 0 */
			/* time limit is 15 msec 					  */
			/* With this stimulus find g_na treshold that depolarizes all	  */
			/* sections beyond 0 mV during simulation			  */

		n = numarg()
		if (n>0) where  = $1 else where = 1  /* 0 soma; 1 all; 3 only terminations */
		if (n>1) simMode = $2
		if (n>2) dx 	= $3 else dx = 0.005 /* accuracy */
		if (n>3) hi 	= $4 else hi = 3000      /* max */
		if (n>4) lo 	= $5 else lo = 0	/* min */

		/* stimulus threshold (if necessary) */

		if (simMode!=1) {
		active_set(0)                    /* g_na = 0 */
		st_intensity = threshold(0,2,0.001) /* find stimulus threshold node */
		print "stimulus threshold ", st_amp, "nA\t(nodal AP within 10 ms; g_na = 0)" 
		st_set(st_intensity)		/* set amplitude */
				}


		/* g_na threshold for the different conditions */

		value = threshold(1,where,dx,hi,lo)   /* calculate the threshold */

		if (where==1&&simMode==0) { 	nathreshold        = value
						print "nathreshold = ", nathreshold}
		if (where==1&&simMode==1) { 	nathresholdvclamp  = value
						print "nathresholdvclamp = ",nathresholdvclamp}
		if (where==3&&simMode==0) { 	nathreshold2       = value
						print "nathreshold2 = ", nathreshold2}
		if (where==3&&simMode==1) { 	nathresholdvclamp2 = value
						print "nathresholdvclamp2 = ",nathresholdvclamp2}


		if (where==0&&simMode==0) { 	nathreshold2       = value
						print "nathreshold = ", nathreshold}
		if (where==0&&simMode==1) { 	nathresholdvclamp = value
						print "nathresholdvclamp = ",nathresholdvclamp}

		active_set() // set all the conductances back to standard
		}

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

func threshold() 	{ local n,what,lo,hi,dx,where
		/* find the na channel density threshold at which the
		   distal dendrite is depolarized > 0 mV 
		   search interval [$1 .. $2] */
		/* does stimulus *and* g_na thresholds */

		n = numarg()
		if (n>0) type	= $1 else type  = 1     /* threshold: 0=stim, 1=g_na */
		if (n>1) where	= $2 else where = 1     /* where: 0 soma, 1 all, 2 node */ 
		if (n>2) dx 	= $3 else dx 	= 0.005 /* accuracy */
		if (n>3) hi 	= $4 else if (type==1) hi = 300  else hi = 2  /* max */
		if (n>4) lo 	= $5 else lo 	= 0	   /* min */


		/*** help */

		if (type==0)  sprint(str1,"stimulus threshold")
		if (type==1)  sprint(str1,"g_na threshold")
		if (where==0) sprint(str1,"%s at soma",str1) 
		if (where==1) sprint(str1,"%s at all sections",str1) 
		if (where==2) sprint(str1,"%s at node",str1) 
		if (where==3) sprint(str1,"%s at terminations",str1) 
	       		      sprint(str1,"%s [%0.0f-%0.0f]",str1,lo,hi)
		if (type==0)  sprint(str1,"%s nA; accuracy %g",str1,dx)
		if (type==1)  sprint(str1,"%s pS/um2; accuracy %g",str1,dx)
		print str1

		if (type==1)  sim_set(threshold_dur) else sim_set(10)

	while (abs(lo-hi)>dx) {if(thresh((lo+hi)/2,where)) hi=(lo+hi)/2 else lo=(lo+hi)/2}

		sim_set()
		return hi
		      }

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


func thresh() { local where
		/* set na channel density or stimulus intensity ($1), then run simulation,   
		   if all sections are depolarized > 0 mV return 1 else 0 */ 	
		ok=1

		if (type==1) { /* nathresh */
				g_na	= $1
				active_set(g_na)}

		if (type==0) { /* stthresh */
				st_amp	= $1
				st_set(st_amp) }

		where 	= $2

		sim()

 if (where==0) if(Soma.sec.vpeak_pk(0.5)<waveform.max-0.01) return 0 
 if (where==1) forsec all for(x) if (vpeak_pk(x) < 0)              return 0   
 if (where==2) if (node.vpeak_pk < 10)                    return 0
 if (where==3) forsec terminations if(vpeak_pk(1) < 0)  	   return 0
 if (where==4) forsec all for(x) if(vpeak_pk(x)<waveform.max-0.01)  return 0   
 if (where==5) forsec terminations if (vpeak_pk(1)<waveform.max-0.01)  return 0   
		return 1 	

	      }







/*****************************************************************/
/* AP decay - other measure than thresholds                      */
/*****************************************************************/

proc APdecay() {local g,min,max,areavar,y,Max,Min,n,APdst,apo200,toggle,active,AP

	/* calculates distance of AP half-decay & AP rel. amplitude at 200 um */
	/* APhalf_pass,  AP200_pass */
	/* if numarg() == 1 then active case */

	if (electrotonicL) { L_switch(0) toggle=1 } else toggle = 0

	n = numarg()
	if (n>0) simMode=$1 else simMode=1
	sim()         // run the simulation


	/* APhalf */
	APdst = 9999
	origin.sec  AP = vpeak_pk(originx) - vonset_pk(originx)

	forsec all for(x){
if(vpeak_pk(x)-vonset_pk(x)<AP/2 && fdistance(x)<APdst) {APdst=fdistance(x)}}


	if (simMode==3) APhalf_pass=APdst   else APhalf=APdst 
	if (simMode==3)	fprint("passive ") else fprint("active  ")
	fprint("AP halfdecay distance %5.0f\t",APdst)

	/* AP200 */
	AP200v = new Vector()
	forsec all {	g = (maxdist()-200)/fL()  // 0 < g < 1 // -> section is within
			if(g>=0&&g<=1) { AP200v.append(vpeak_pk(g) - vonset_pk(g))  }}

	apo200 = mean(AP200v)/AP

	if (simMode==3) AP200_pass=apo200  else  AP200=apo200
	fprint("normalised AP200 %5.3f\n",apo200)

	if (toggle) L_switch(1)
	if (simMode==3) simMode =1
		}


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


proc APdecay_sensitivity() { local i, mult, n, mode
		
		n = numarg()
		if (n>0) simMode = $1 else simMode = 1
		if (simMode==1)  remove_axon()

		mult = 1.4

		for i = 0,2 sens[i] = new Vector()
		sens[0].indgen(0,nathresholdvclamp*mult,nathresholdvclamp/11)
		for i = 0,sens[0].size()-1 { 
						active_set(sens[0].x[i]) 
						APdecay()   
					  	sens[1].append(APhalf)
					  	sens[2].append(AP200)
					}
		sigmoidal_calc()
		active_set()
			}






func sigmoidal() { /* sigmoidal fit for AP200 titration */

		  ref 	= sens[0].c
		  shalf  = $&2[0]
		  steep = $&2[1]
		  range = $&2[2]
		  basis = $&2[3] 
		  ref.sub(shalf).mul(-1).div(steep)
		  if (ref.max > 709) return 1000
		   ref.apply("exp").add(1).div(range).apply("flip")
		   ref.add(basis)
		   return ref.meansqerr(sens[2])
		}



proc sigmoidal_calc() { /* sigmoidal fits for AP200 */

			n = numarg()

			param[0] = 15
			param[1] = 1
			param[2] = 0.5
			param[3] = sens[2].x[0]

			attr_praxis(1e-6,4,0)
			fit_praxis(4,"sigmoidal",&param[0])

			if (n>=1) fig(sens[0],sens[2])
			if (n>=1) fig(sens[0],ref,0,1,2)

			print "half ",shalf,"\tsteep ",steep , "\n"

			AP200_half  = shalf
			AP200_steep = steep
			AP200_basis = basis
			AP200_range = range
			}


proc scrappy() {
		ActiveModel = act0
		cell_name($s1)
		number = new Vector()

		readveca("nathreshold")
		for i=0,2 sens[i].resize(15)
		sigmoidal_calc()
		ref4 = ref.c
		fig(sens[0],sens[2],1,0,1)
		fig(sens[0],ref4,1,1,2)
		figlab(cellname,0.3,0.9)
		}

Loading data, please wait...