Effects of the membrane AHP on the Lateral Superior Olive (LSO) (Zhou & Colburn 2010)

 Download zip file   Auto-launch 
Help downloading and running models
Accession:143114
This simulation study investigated how membrane afterhyperpolarization (AHP) influences spiking activity of neurons in the Lateral Superior Olive (LSO). The model incorporates a general integrate-and-fire spiking mechanism with a first-order adaptation channel. Simulations focus on differentiating the effects of GAHP, tauAHP, and input strength on (1) spike interval statistics, such as negative serial correlation and chopper onset, and (2) neural sensitivity to interaural level difference (ILD) of LSO neurons. The model simulated electrophysiological data collected in cat LSO (Tsuchitani and Johnson, 1985).
Reference:
1 . Zhou Y, Colburn HS (2010) A modeling study of the effects of membrane afterhyperpolarization on spike interval statistics and on ILD encoding in the lateral superior olive. J Neurophysiol 103:2355-71 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism: Auditory brainstem;
Cell Type(s): Lateral Superior Olive (LSO) cell; Abstract integrate-and-fire leaky neuron;
Channel(s): I_AHP;
Gap Junctions:
Receptor(s): GabaA; Glutamate;
Gene(s):
Transmitter(s):
Simulation Environment: NEURON; MATLAB;
Model Concept(s): Action Potential Initiation; Simplified Models; Spike Frequency Adaptation; Depolarization block; Audition;
Implementer(s): Zhou, Yi [yizhou at bu.edu];
Search NeuronDB for information about:  GabaA; Glutamate; I_AHP;
// *****************  LSO model 1.1  **********************
// ------------  Introduction   --------------------------------  
// 01/01/2011
// This model simulates spike interval statistics and ILD responses in the Lateral Superio Olive (LSO) of mammalian brainstem. 
// Model configurations and simulation methods are described in 
// Zhou and Colburn (2010) J. Neurophysiology, 103:2355-2371   

// ----- model  -----
// Leaky Integrate&Fire cell model with Afterhyperpolarization (AHP) chanel

// ----- current input  -----
// currents with Gaussian noise amplitude, described by mean and std

// ----- output  -----
// rate, ISI, gAHP, PSTH, and serial correlation statistics 






// ------------  file structure   ------------------------
// Part I: Membrane and synapse setup
// Part II: Data processing and spike-time analyses
// Part III:User interface setup 
//---------------------------------------------------------

// ---------- Code start here ----------

seed_x=xred("what is the seed for today", 1, 1, 1000) //string,default,min,max

strdef inputDir,outputDir
outputDir="output"

load_file("nrngui.hoc")
cvode.active(0)


// Part I: Membrane and synapse setup

//****  Cell Model*****
create IF_cell
access IF_cell

IF_cell {
	nseg = 1
       	diam = 20  //um
       	L = 50    // 1000pi um^2
       	Ra=150   //ohm.cm
       	cm=1     //uF/cm2
        
       	insert pas
	
	e_pas=-65  //mV
	g_pas=0.001 //0.0004S/cm2  //  tau=1ms->0.001uS    tau=0.2ms->0.005uS     0.25ms->0.004uS

	area_cell=PI*diam*L
		
	print "area[um^2]= ",area_cell
	print "g_pas[S/cm2]=  ",g_pas
        G_pas=area(0.5)*g_pas/100  //uS
	print "total axon_conductance[uS]= ",G_pas
	print "tau_rest[ms]= ",1/(g_pas*1000)
	print "======================"
  
    	
}

  // add absolute refractory period by voltage clamping 

  V_thres=-50   // -50: 15mV above the rest

  objref refrac_abs, netCell_abs
	refrac_abs=new Refrac_abs(0.5)
	refrac_abs.tr=2  // absolute RP, if tr>1, may cause next input with smaller tau
	refrac_abs.e=-65   // [mV] reversal potential
	refrac_abs.gr=10   //[uS] conductance during absolute refractory period  : 10 times the G_pas

	
	delays=0
	weight=0     // not used in *.mod 
	netCell_abs=new NetCon(&IF_cell.v(0.5),refrac_abs,V_thres,delays,weight)

  // put relative refractory period inside the Cell   : not used in this program -- refrac_ahp serves both
  
  objref refrac_rel, netCell_rel
	refrac_rel=new Refrac_rel(0.5)
	refrac_rel.tau=5  // the decay time constant
	refrac_rel.e=-65   // [mV] reversal potential
	refrac_rel.gr=0   //[uS] conductance after the abs refractory period

	
	refrac_rel_delay=refrac_abs.tr   // if delay=refrac_abs.tr, then send out signal after the abs refractory 
	//refrac_rel_delay=0
	weight=0		// not used in *.mod 
	netCell_rel=new NetCon(&IF_cell.v(0.5),refrac_rel,V_thres,refrac_rel_delay,weight)


  
  // add AHP channel
  
  objref refrac_AHP, netCell_AHP
	refrac_AHP=new Refrac_rel(0.5)
	refrac_AHP.tau=5  // the decay time constant
	refrac_AHP.e=-65   // [mV] reversal potential
	refrac_AHP.gr=0.08   //[uS] conductance after the abs refractory period

	
	refrac_AHP_delay=refrac_abs.tr    // if delay=refrac_abs.tr, then send out signal after the abs refractory 
	//refrac_AHP_delay=0
	weight=0        // not used in *.mod 
	netCell_AHP=new NetCon(&IF_cell.v(0.5),refrac_AHP,V_thres,refrac_AHP_delay,weight)
  


 //------- END Cell Model -------------------------------------
 //----------------------------------------------------------





  //**** procedure  "noisy current stimulus"  ****** 
  //---------------------------------------------------------
  objref current_inj
  current_inj=new current_gauss(0.5)
  current_inj.del=0
  current_inj.dur=200
  current_inj.mean=1.4
  current_inj.std0=0.4
	
  
  //------- END current stim -------------------------------------
  //----------------------------------------------------------
 

 



  //**** procedure  "stimulus and target connection"  ****** 
  //---------------------------------------------------------
  
  //----     parameter description     ----------------------
  // 
  //syn_contra[]: I alpha synapses, each group has 50 synapses 
  //syn_ipsi[]: E alpha synapses,each group has 50 synapses  
  //
  //gen_contra[]: I stimuli to syn_contra[]
  //gen_ipsi[]: E stimuli to syn_ipsi[]
  
  //netcon_contra[]: connection on contra
  //netcon_ipsi[]: connection on ipsi
  //
  
  // Synaptic inputs not used  in this model
  
  syn_number=50   

  objectvar syn_contra[syn_number], gen_contra[syn_number]
  objectvar syn_ipsi[syn_number], gen_ipsi[syn_number]

 
  access IF_cell 
 
   //----- virtual AVCN and MNTB spiking activity  
   for i=0, syn_number-1{
            gen_contra[i] = new VecStim(0.5)
			gen_ipsi[i] = new VecStim(0.5)
  	} 

	// contralateral inhibition from MNTB  
	for i=0,syn_number-1 {
		syn_contra[i] = new Alpha(0.5)  // scatter along the whold dend
		syn_contra[i].tau1 =0.1  //ms    may add some jitters
		syn_contra[i].tau2 =1   //ms
		syn_contra[i].con=0    //nS- mho
		syn_contra[i].e =-70   // mV    
		}                            //syn.i --- nA
   
	// ipsilateral excitation from AVCN
 	for i=0,syn_number-1 {
		syn_ipsi[i] = new Alpha(0.5)
		syn_ipsi[i].tau1 =0.1  //ms    may add some jitters
		syn_ipsi[i].tau2 =1   //ms
		syn_ipsi[i].con=0    //nS-mho
		syn_ipsi[i].e =0   // mV
    	}                            //syn.i --- nA

  
    objref netcon_contra[syn_number],netcon_ipsi[syn_number]

    x_threshold=0.05  // x_threshold was not used to control the input firing rate in this model 
    i_delay_contra=0  //  response delay
    e_delay_ipsi=0  

   for i=0,syn_number-1  { 
		netcon_contra[i]=new NetCon(gen_contra[i], syn_contra[i]) 
		netcon_contra[i].delay=i_delay_contra
      	netcon_contra[i].weight=0.001  // delay gain[uS]

      	netcon_ipsi[i]=new NetCon(gen_ipsi[i], syn_ipsi[i]) 
        	netcon_ipsi[i].delay=e_delay_ipsi
        	netcon_ipsi[i].weight=0.001  // delay gain[uS]
       }
    
  //-------- END connection ---------------------------------
  //--------- END PART I ------------------------------------

   
  
  
  // Part II: Data processing and spike-time analyses
    
  // *** procedure   "record AP numbers"  ***
  //-----------------------------------------------
  //----     parameter description     ------------ 
  // bin:[ms] bin width of AP
  // total:   total number of entries recorded in v_all
  // v_ISI: [ms] record the event time
  // v_voltage:  record occurrence of event at each dt
  // pre_E: [ms] record the precell event time of netcon_contra[0]
  
  //sptime_out: spike time for all trials, exporting to matlab
  
  tstop=200+e_delay_ipsi   // all inputs last 200ms, maximum stim time=300ms 
  start_ss=40+e_delay_ipsi   // steady state starts at 40 ms  after the E onset; reset in the reset()
  

  dt=0.025    // numerical integration step
  bin=0.05    // 50us time bin size, bin AP
  
  total=tstop/bin
   
  N_run=50
  


  access IF_cell

  objref apc,v_ISI,v_voltage,pre_E, sptime_out

  objref g_AHP
  
  v_ISI=new Vector()
  v_voltage=new Vector(tstop/dt+2,0)
  pre_E=new Vector()
  
  AP_threshold=V_thres  //-50 mV 
  apc = new APCount(0.5)
  apc.thresh=AP_threshold


  
  
  //==== for test input ISI and PS
  //netcon_contra[0].record(v_ISI)
  
  //==== for test output ISI and PS
  apc.record(v_ISI) //record event time
  v_voltage.record(&IF_cell.v(0.5),dt) //record the event: event time=index*dt
  
  netcon_contra[0].record(pre_E)//record input event time 	


  // === monitoring the AHP conductance
   g_AHP=new Vector(tstop/dt+2,0)   
   g_AHP.record(&refrac_AHP.g,dt)

  //-------- END record AP ------------------------
  //-----------------------------------------------
  
  
  
  

  //****  function  "PSTH"  **** 
  //-----------------------------------------------
  //----     parameter description     ------------ 
  //v_save: bin v_voltage for each run
  //v_all: summation of all v_save 
  //v_psth: output PSTH vector
  //x_psth: AN PSTH vector
  //t_psth: time axis
  //v_trigger: spike-triggered avergage of sub-membrane potential 

  
  objref v_save,v_all,v_psth
  objref v_trigger, t_trigger, v_average, v_trigger_mean,v_trigger_std,g_trigger, trigger_time
  objref g_psth,t_psth, v_stand
  objref zero_pad
  objref d, xtime
  objref v_temp


  bin_psth=10*bin  //0.5 ms 
  trigger_length=200*bin  //ms

  nrow_trigger=N_run
  ncol_trigger=trigger_length/bin

  v_average=new Matrix(nrow_trigger,ncol_trigger)

  proc get_psth() { 
        v_save=new Vector(L_PTH)
  		v_save.spikebin(v_voltage,AP_threshold) //(v_source, threshold)
		v_save.rebin(v_save,bin/dt)  // rebin the time axis
		v_all.add(v_save)   // add all trials, the number of spikes
	
		//*******   measure the pre-spike average of membrane potential during the steady-state
		v_trigger=new Vector(L_PTH,0)
		t_trigger=new Vector(L_PTH,0)
		v_trigger.add(v_voltage)
		t_trigger.spikebin(v_voltage,AP_threshold) // all-none train

		v_trigger.rebin(v_trigger,bin/dt) 
		v_trigger.div(bin/dt)
		t_trigger.rebin(t_trigger,bin/dt)  // rebin the time axis

		// only keep the steady-state phase
		start_index=start_ss/bin
		v_trigger.remove(0, start_index)
		t_trigger.remove(0, start_index)
		
		v_temp=new Vector()
		v_temp.correl(v_trigger,t_trigger)
		v_temp=v_temp.c(0,ncol_trigger-1)
		v_temp.div(t_trigger.sum)   // normalized by number of spikes
		v_average.setrow(N_run-1, v_temp)
		
  	}  



  proc psth() {
		length_psth=tstop/bin_psth  //length of period points
		
		v_psth=new Vector(v_all.size,0)   
		v_psth.add(v_all)
		
		
		// rebin the time axis    // can't compress if ratio=1
		if ( bin_psth > bin ) { v_psth.rebin(v_psth,bin_psth/bin) } 
			
		v_psth.div(N_run*bin_psth*0.001)   //  #spikes/bin-> #spikes/ms-> #spikes/sec
		
		
		t_psth=new Vector(length_psth)
		t_psth.indgen(bin_psth)
		
		
		g_psth= new Graph()
		g_psth.size(0,tstop,0,1000)
  		g_psth.align(0.5,0.2)
  		g_psth.label(0.9,0.2,"time")

  		g_psth.color(2)
  		g_psth.label(0.6,0.8,"LSO(#/sec)")
  		v_psth.plot(g_psth,t_psth,2,2)
       
		
		// average overll all trials for spike-triggered potentials
		v_temp=new Vector()
		for i=0,ncol_trigger-1 {
			v_average.getcol(i,v_temp)
			v_trigger_mean.x[i]=v_temp.mean()
			v_trigger_std.x[i]=0 //v_temp.stdev()
			}
		
		trigger_time=new Vector(ncol_trigger,0)
		trigger_time.indgen(bin)
   }
  
  //---  END function  "PSTH" ---------------------
  //-----------------------------------------------
  

 

  //****  function  "ISI"  **** 
  //-----------------------------------------------
  //----     parameter description     ------------ 
  // get_ISI(): get histogram for each run
  //v_ISI: record each event time
  //v1,v2,v1_ss,v2_ss: operating vector 
  //bin_ISI : bin width of ISI	
  //hist_ISI: histogram of ISI for all run
  //hist: histogram of ISI for each run
  //mean_T, sq_T, num_T: ISI numbers at each run
  
  objref v1,v2,v1_ss,v2_ss
  objref g_AHP_index
  objref hist_ISI,hist_ISI_ss,hist,hist_ss
  objref F,recover	
  
  objref serialISI,vector_noZero,time_noZero,meanfit_ISI, meanfit_AHP

  objref xtime,xtime_recover,g_ISI,g_serialISI
  objref mean_T, sq_T, num_T
  objref temp,t1, temp_ss, temp_count,temp_AHP
   
  objref serialAHP, AHP_count
  objref mean_AHP,sq_AHP

  objref rate_infor
  
  mean_ISI=0  // mean T for each run; only use for records
  OVF=0       // overflow out of topbins
  EI=0        // entrainment index   
  mean_ISI_total=0
  std_ISI=0
  sq_sum=0
  CV_ISI=0
  
  mean_AHP_total=0
  std_AHP=0


  serial=1  // the order of serial dependence !!! 


  proc get_ISI() { 
  		
		// ****    For the whole section of responses ******
  		v2=new Vector(L_ISI,0)
  		v1=new Vector(L_ISI,0)
  		
  		
  		v1.add(v_ISI)
  		v2.add(v_ISI)
  		v1.rotate(1) // right shift
  		v2.sub(v1)
  		v2.rotate(-1)
  		v2.resize(v2.size-1)  // cut off the last interval
  		
  		
  		//*** spike times adding all trials 
  		if (v_ISI.max!=0) sptime_out.append(v_ISI)   
  		 
  		 
  		//*** for calculate the square sum 
  		temp = new Vector(v2.size,0)
  		temp.add(v2)
  		temp.sub(v2.mean)  
  		
  		mean_T.x[N_run-1]=v2.mean
  		num_T.x[N_run-1]=v2.size
  		sq_T.x[N_run-1]=temp.sumsq()
  		
  		mean_ISI=v2.mean + mean_ISI  // mean ISI for each run
  
  		print "% rate=",(v2.size+1)/(tstop/interval)
  		print " spikes=",v2.size+1
  		print "**********************************"
  		
  		hist=v2.histogram(0, topbin+OVFbin, bin_ISI)  //bin width=?
		
  		if(hist.size-hist_ISI.size==1) { hist.resize(hist_ISI.size)} 
  		hist_ISI.add(hist)      // units N_run*spikes/bin_ISI


	// ****    For the steady-state section of responses ******
  		v2_ss=new Vector(L_ISI,0)
  		v1_ss=new Vector(L_ISI,0)
  		temp = new Vector(L_ISI,0)
		temp.indvwhere(v_ISI, "[]", start_ss, tstop)   // get the indice for spike between start_ss -> tstop
		R40=v_ISI.size-temp.size     // BONUS!! the spikes number for the 40 ms onset phase
		
		if(temp.size<=1) { temp=new Vector(2,0) }  // not saving spikes less than one

  		v2_ss = v_ISI.ind(temp)
  		v1_ss = v_ISI.ind(temp)
  		v1_ss.rotate(1) // right shift
  		v2_ss.sub(v1_ss)
  		v2_ss.rotate(-1)
  		v2_ss.resize(v2_ss.size-1)  // cut off the last interval

		//*** save the indice of g_AHP right after a spike(i.e., before refrac and the next spike) during the SS

		g_AHP_index=new Vector(L_ISI,0)
		g_AHP_index=v_ISI.ind(temp)  // the spike times during the SS
		g_AHP_index.div(dt)          // sp index
		g_AHP_index.resize(g_AHP_index.size-1)


		// *** save the mean and square sum of g_AHP during the SS
		temp = new Vector(g_AHP.size,0)
  		temp.add(g_AHP)
		temp.rotate(-start_ss/dt)  // left shift
		size_ss=g_AHP.size-start_ss/dt
		temp.resize(size_ss)
  		
		mean_AHP.x[N_run-1]=temp.mean()
 		temp.sub(temp.mean)  
  		sq_AHP.x[N_run-1]=temp.sumsq()
	
	// **** measure the serial dependence of ISI for the steady-state
        
	temp_ss=new Vector(v2_ss.size,0)
	temp_ss.add(v2_ss)
	temp_ss.resize(temp_ss.size-serial)  // cut off the last interval, which has no serial ISI following
	
	for i=0, topbin/bin_ISI {
		temp = new Vector(200,0)   // for saving the indice 
		temp.indvwhere(temp_ss,"[)",i*bin_ISI, (i+1)*bin_ISI)   // indice of intervals belongs to each bin
		temp_AHP=new Vector(200,0)   // for saving the indice of g_AHP for each bin 

		if (temp.size>=1) { 
			// the  g_AHP that starts at the current interval 
			AHP_count.x(i)=AHP_count.x(i)+temp.size
			temp_AHP=g_AHP_index.ind(temp)
			serialAHP.x(i)=serialAHP.x(i)+g_AHP.ind(temp_AHP).sum   // sum over previous run

			// the next interval
			temp.add(serial)   
			temp_count.x(i)=temp_count.x(i)+temp.size
			serialISI.x(i)=serialISI.x(i)+v2_ss.ind(temp).sum   // sum over previous run
				   } 	
		}
	
	// **** measure the histograms 
		hist_ss=v2_ss.histogram(0, topbin+OVFbin, bin_ISI)  //bin width=0.2
  		if(hist_ss.size-hist_ISI_ss.size==1) { hist_ss.resize(hist_ISI_ss.size)} 
  		hist_ISI_ss.add(hist_ss)   // units N_run*spikes on each bin
  	}	
  	
  proc ISI() { 
  		hist_ISI.div(hist_ISI.sum)  //normalize units: prob()
		OVF=hist_ISI.sum(topbin/bin_ISI+1, (topbin+OVFbin)/bin_ISI-1) // may overflow by T  		
  		hist_ISI.resize(topbin/bin_ISI+1)

	//******  cut off the first element of sptime 
		sptime_out.rotate(-1)
  		sptime_out.resize(sptime_out.size-1) 
		
	//******   steady-state   *********  
		hist_ISI_ss.div(hist_ISI_ss.sum)  //normalize
		hist_ISI_ss.resize(topbin/bin_ISI+1)		
		
		//** Hazard function  : f(x)/[1-F(x)]  : using F(x)=1-F(x)
		F=new Vector(hist_ISI_ss.size-1,0)
		for i=0, hist_ISI_ss.size-2  {
			F.x[i]=hist_ISI_ss.sum(i+1,hist_ISI_ss.size-1)	
		}

		recover=new Vector(hist_ISI_ss.size,0)
		temp = new Vector(hist_ISI_ss.size,0)
		temp.indvwhere(F, ">=", 0.05)    // save 1-F(x)  with prob > 0.05  instead of hist_ISI_ss>0.1

		recover=hist_ISI_ss.ind(temp)
		F.resize(recover.size)
		recover.div(F)

		//*** time axis  
		xtime=new Vector(hist_ISI.size,0)
  		xtime.indgen(bin_ISI)
		xtime_recover=new Vector(recover.size,0)
  		xtime_recover.indgen(bin_ISI)

		//** all vairals convert to 1/sec  units  
        hist_ISI.div(bin_ISI*0.001)
		hist_ISI_ss.div(bin_ISI*0.001)
		recover.div(bin_ISI*0.001)

	   //******   serial ISI   *********  
		serialISI.div(temp_count)
		
       //******   serial AHP   *********  
		serialAHP.div(AHP_count)		

  		g_ISI.size(0,topbin,0,1.2*recover.max)
  		g_ISI.align(0.5,0.5)
  		//g_ISI.label(0.2,0.9,"ISI")
  		g_ISI.beginline()  

  		//ISI of the whole responses (not shown)
  		//hist_ISI.mark(g_ISI,xtime,"o",6,1)    // size 6; color 1
  		//hist_ISI.line(g_ISI,xtime,1,2)
		//g_ISI.color(1)
		//g_ISI.label(0.7,0.2,"whole")  		
		
		hist_ISI_ss.mark(g_ISI,xtime,"+",6,2)
  		hist_ISI_ss.line(g_ISI,xtime,2,2)
		g_ISI.color(2)
		g_ISI.label(0.7,0.5,"ISIss")    		
		recover.line(g_ISI,xtime_recover,3,3)   // color 3; size 3
		g_ISI.color(3)
		g_ISI.label(0.7,0.9,"Recovery Func") 


    //******  plot the serial dependence of ISIs   ***********

		g_serialISI.size(0,topbin,0,1.2*serialISI.max)
  		g_serialISI.align(0.5,0.5)
  		g_serialISI.beginline()  
		serialISI.mark(g_serialISI,xtime,"o",6,3)
		g_serialISI.color(3)
		g_serialISI.label(0.7,0.5,"serial ISIs")    
		
		serialAHP.mul(100) // zoom in X100
		serialAHP.mark(g_serialISI,xtime,"+",6,5)
		g_serialISI.color(5)
		g_serialISI.label(0.7,0.9,"serial AHPs (X100)")    
		serialAHP.div(100) // zoon back to the original
	
  		//****   calculate the mean and the std of intervals  
   		t1=new Vector(N_run,0)
   		t1.add(mean_T)
   		t1.mul(num_T)
   		t2=t1.sum()
   		mean_ISI_total=t2/num_T.sum()  //overall ISI mean 
  		
  		t1=new Vector(N_run,0)
  		t1.add(mean_T)
  		t1.sub(mean_ISI_total)
  		t1.mul(t1)
  		t1.mul(num_T)
  		sq_sum=sq_T.sum()+t1.sum()
  		
  		if(num_T.sum()>1) {
  		    std_ISI=sqrt(sq_sum/(num_T.sum-1))   //Overall ISI std
  		    } else { std_ISI=sqrt(sq_sum) } 
  		
  		if(mean_ISI_total>0) {
  				CV_ISI=std_ISI/mean_ISI_total
  				} else{ CV_ISI=1}            //coefficient of variation
  
  	//****   calculate the mean and the std of g_AHP during the SS 
 
   		t1=new Vector(N_run,0)
   		t1.add(mean_AHP)
   		mean_AHP_total=t1.mean()  //overall mean(g_AHP) 
  		t1.sub(mean_AHP_total)
  		t1.mul(t1)
  		t1.mul(size_ss)
		sq_sum_AHP=sq_AHP.sum()+t1.sum()
  		std_AHP=sqrt(sq_sum_AHP/(size_ss*N_run))
  
   }

  proc meanFit() { 
		// ** linear fit for the serial_ISI data 
		temp = new Vector()
		temp.indvwhere(serialISI, ">", 0)   // fit elements >0

		vector_noZero=new Vector()
		vector_noZero=serialISI.ind(temp)

		time_noZero=new Vector()
		time_noZero=xtime.ind(temp)

		meanfit_ISI=new Vector(time_noZero.size,0)
		p1_ISI=-0.3  // don't choose 0 as the inital value for the fit function
		p0_ISI=8
		error = vector_noZero.fit(meanfit_ISI, "line", time_noZero, &p1_ISI, &p0_ISI)
		print "ISI mean square error=",error
		strdef fit_func   //p1x+p0
		sprint(fit_func, "T(n+1)=%2.2fT(n)+%2.2f",p1_ISI,p0_ISI)   

		meanfit_ISI.line(g_serialISI,time_noZero,1,2)  //color 1;  line2
		g_serialISI.color(3)
		g_serialISI.label(0.7,0.2,fit_func)       

		
		// ** linear fit for the serial_AHP data 
		temp = new Vector()
		temp.indvwhere(serialAHP, ">", 0)

		vector_noZero=new Vector()
		vector_noZero=serialAHP.ind(temp)

		time_noZero=new Vector()
		time_noZero=xtime.ind(temp)

		meanfit_AHP=new Vector(time_noZero.size,0)
		p1_AHP=0.1  // don't choose 0 as the inital value for the fit function
		p0_AHP=8
		error = vector_noZero.fit(meanfit_AHP, "line", time_noZero, &p1_AHP, &p0_AHP)
		print "AHP mean square error=",error
		strdef fit_func   //p1x+p0
		sprint(fit_func, "T(n+1)=%2.2fT(n)+%2.2f",p1_AHP,p0_AHP)   

		meanfit_AHP.line(g_serialISI,time_noZero,1,2)  //color 1;  line2
		g_serialISI.color(5)
		g_serialISI.label(0.7,0.9,fit_func)       

  }
  

  //****  function  Write to file  **** 
  //-----------------------------------------------
  //----     parameter description     ------------ 
  
  objref f_hist, f_histss,f_recover,f_serial,f_serialAHP, f_psth, f_rate,f_sptime
  objref sys_dir  // the output dir	  

    strdef stimulus
	stimulus="temp"

  proc WriteFile() { local a,b,c,d,e

	chdir(outputDir) 

	a=refrac_abs.tr
    b=refrac_AHP.gr
	c=refrac_AHP.tau
	d=current_inj.mean
	e=current_inj.std

	strdef stimulus   //abs_AHPg_AHPtau
	sprint(stimulus, "abs(%2.1f)_AHPg(%2.2f)_AHPtau(%2.1f)_mean(%2.1f)_std(%2.1f)",a,b,c,d,e)  

	strdef filename   //variable_stimulus
	sprint(filename, "hist_%s",stimulus)  
	    f_hist=new File(filename)
        f_hist.wopen()
        hist_ISI.vwrite(f_hist,4)   
        f_hist.close()

	//-----------------------------

	strdef filename
	sprint(filename, "histss_%s",stimulus)   
	    f_histss=new File(filename)
        f_histss.wopen()
        hist_ISI_ss.vwrite(f_histss,4)   
        f_histss.close()

	//-----------------------------

	strdef filename
	sprint(filename, "recover_%s",stimulus) 
	    f_recover=new File(filename)
        f_recover.wopen()
        recover.vwrite(f_recover,4)   // sptime output 
        f_recover.close()
        
        
	//-----------------------------

	strdef filename
	sprint(filename, "serial_%s",stimulus) 
	    f_serial=new File(filename)
        f_serial.wopen()
        serialISI.vwrite(f_serial,4)   
        f_serial.close()

	//-----------------------------

	strdef filename
	sprint(filename, "serialAHP_%s",stimulus) 
	    f_serialAHP=new File(filename)
        f_serialAHP.wopen()
        serialAHP.vwrite(f_serialAHP,4)   
        f_serialAHP.close()

	//-----------------------------

	strdef filename
	sprint(filename, "psth_%s",stimulus) 
	    f_psth=new File(filename)
        f_psth.wopen()
        v_psth.vwrite(f_psth,4)   
        f_psth.close()

	//-----------------------------

	strdef filename
	sprint(filename, "rate_%s",stimulus) 
	    f_rate=new File(filename)
        f_rate.wopen()
        rate_infor.vwrite(f_rate,4)   //rate_infor 
        f_rate.close()

	//-----------------------------

	strdef filename
	sprint(filename, "sptime_%s",stimulus) 
	    f_sptime=new File(filename)
        f_sptime.wopen()
        sptime_out.vwrite(f_sptime,4)   //spike times for all runs
        f_sptime.close()
        
	chdir("..")  
	}	


 //**** procedure "run N-run times, save firing time, compute mean and std of firing rates" ****
  //---------------------------------------------------
  //----     parameter description     ------------
  // N_run: number of trials 
  // v_AP: firing rates of each run
  //
  // rerun() : main function, repeat every tstop
  // integrate(): integration function
  // getcon(): upgrade the new synapse information according to syn_contra[0] and syn_ipsi[0]
  // getstim():upgrade the new stimuli information according to gen_contra[0] and gen_ipsi[0]
  
  objref v_AP,v_AP_40,mean_T, sq_T, num_T   // vec for one IPD with N_run 
  objref mean_AHP,sq_AHP
  objref rate_infor
 
  v_AP=new Vector(N_run,0)    // the whole 200 ms firing rates	
  v_AP_40=new Vector(N_run,0)  // the 40 ms transient firing rates	

interval=4    // Global;  don't change, 
  
 proc reset() {
  
    N_run=50
    tstop=200+e_delay_ipsi   // all inputs last 200ms, maximum stim time=300ms 
    start_ss=40+e_delay_ipsi     // define the start of the steady-state for ISI()
    current_inj.dur=tstop

    interval=4 	// arbitrary length; Global;  don't change, 
	bin_ISI=bin*4   // bin*2=0.1ms  bin*4=0.2ms  ****!!!!!!!!!!!!!****!!!!!!!!!!!!!
  	topbin=int(5*interval/bin_ISI)*bin_ISI  // 20 ms for binning and 20 ms for OVF
  	OVFbin=int(5*interval/bin_ISI)*bin_ISI
    bin_psth=10*bin  //0.5 ms 
        
    
//***** reclaim all the vectors   *************	
    v_AP=new Vector(N_run,0)
	v_AP_40=new Vector(N_run,0)

  	mean_T=new Vector(N_run,0)
  	sq_T=new Vector(N_run,0)
  	num_T=new Vector(N_run,0)
	
	mean_AHP=new Vector(N_run,0)
  	sq_AHP=new Vector(N_run,0)
  	
	nrow_trigger=N_run
    ncol_trigger=trigger_length/bin
    v_average=new Matrix(nrow_trigger,ncol_trigger)
	v_trigger_mean=new Vector(ncol_trigger,0)
	v_trigger_std=new Vector(ncol_trigger,0)

    R=0
	std=0
    R40=0
	R40_std=0
  	mean_ISI_total=0
   	std_ISI=0
   	sq_sum=0
   	CV_ISI=0
  	mean_ISI=0  // only for check on the each run
  	OVF=0
  	EI=0
    mean_AHP_total=0
	std_AHP=0
  	
  	hist_ISI=new Vector((topbin+OVFbin)/bin_ISI+1,0)
	hist_ISI_ss=new Vector((topbin+OVFbin)/bin_ISI+1,0)
	serialISI=new Vector (topbin/bin_ISI+1,0)
	temp_count=new Vector(topbin/bin_ISI+1, 1e-4)

    serialAHP=new Vector (topbin/bin_ISI+1,0)
	AHP_count=new Vector(topbin/bin_ISI+1, 1e-4)

	v_all=new Vector(tstop/dt+1,0)
  	v_all.rebin(v_all,bin/dt)

	rate_infor=new Vector(13,0)  
  
  //==== for record the spike times for all trials
  sptime_out=new Vector(1,3000)   // the first element is spared later
  }
  
  
  
  
  proc rerun() {
        
  	reset()
    while (N_run>0) {
		// reclaim vector otherwise wrong message like  "Segmentation violation" 	
		v_ISI=new Vector()
  		v_voltage=new Vector(tstop/dt+2,0)
		g_AHP=new Vector(tstop/dt+2,0) 
  		
  		//==== for output ISI and PS
        apc.record(v_ISI) //record event time
  		v_voltage.record(&IF_cell.v(0.5),dt) 
  		
   		g_AHP.record(&refrac_AHP.g,dt)
  		
  		
  		v_init=-65      
  		finitialize(v_init) 
  		fcurrent()
  		integrate()
  		get_ISI()
		get_psth()

        v_AP.x[N_run-1]=apc.n
		v_AP_40.x[N_run-1]=R40
		N_run=N_run-1 
		
	} 

   N_run=50   // for psth() function
   
   R=v_AP.mean()
   if (v_AP.size>1) { R_std=v_AP.stdev() }
   R40=v_AP_40.mean()
   if (v_AP_40.size>1) { R40_std=v_AP_40.stdev() }
   g_ISI.erase_all()
   g_serialISI.erase_all()
   ISI()
   psth()

   
   // **** just checking  ***********
 
   print "R=",R
   print "R_std=",R_std
   
   print "R40=",R40
   print "R40_std=",R40_std
   
   print "T=", mean_ISI_total
   print "T_std=", std_ISI
   print "T_CV=", CV_ISI
   
   print "EI=",  EI
   print "OVF=", OVF

   print "AHP=", mean_AHP_total
   print "AHP_std=", std_AHP
   
//*****  save rate infor into a vector
   rate_infor.x[0]=R
   rate_infor.x[1]=R_std
   rate_infor.x[2]=R40
   rate_infor.x[3]=R40_std
   rate_infor.x[4]=mean_ISI_total
   rate_infor.x[5]=std_ISI
   rate_infor.x[6]=CV_ISI
   rate_infor.x[7]=OVF
   rate_infor.x[8]=mean_AHP_total
   rate_infor.x[9]=std_AHP
   rate_infor.x[10]=bin_ISI
   rate_infor.x[11]=bin_psth
   rate_infor.x[12]=interval
   print "\n----------- END -----------------------"
}
  
  
  proc integrate() {
    
    while (t < tstop) {
    fadvance()  // advance solution 
    
  }
  print "======= info ======="
  print "N_run=",N_run   // show run time 
  print "spike number=",apc.n 
  print "N_input_E=",pre_E.size

  if(v_ISI.size<=1) { v_ISI=new Vector(2,0) }  // not saving spikes less than one
  L_ISI=v_ISI.size
  print "L_ISI=",L_ISI
  L_PTH=v_voltage.size		
  //print "L_PTH=",L_PTH	
}
 
 //-----END rerun() ------
  //------------END PART II---------

  
// Part III:  User interface setup 
  cvode.active(0)
  access IF_cell
  
  objref syn_con
  
  syn_con = new VBox()
  syn_con.intercept(1)       //all following creations go into the "vbox" box

  xpanel("")
  xlabel("Cell")
  xvalue("IF_cell.g[S/cm2]","IF_cell.g_pas")
  xpanel()
  
  xpanel("")
  xlabel("Refrac_abs")
  xvalue("refrac_abs.gr[uS]","refrac_abs.gr")
  xvalue("refrac_abs.e[mV]","refrac_abs.e")
  xvalue("refrac_abs.tr[ms]","refrac_abs.tr")
  
  xlabel("AHP")
  xvalue("AHP.gr[uS]","refrac_AHP.gr")
  xvalue("AHP.e[mV]","refrac_AHP.e")
  xvalue("AHP.tau[ms]","refrac_AHP.tau")
  
  xpanel()

  xpanel("")
  xlabel("Current Stim")
  xvalue("I_mean[nA]","current_inj.mean")
 
  xvalue("I_std_4000Hz[nA]","current_inj.std0")
  //xvalue("I_std[nA]","current_inj.std")

  xvalue("I_start[ms]","current_inj.del")
  xvalue("I_duration[ms]","current_inj.dur")
  
  xpanel()

  syn_con.intercept(0)       //ends intercept mode
  syn_con.map("SYN",0,100,-1,50)              //draw the box and its content


  //----------------
  xpanel("ControlPanel")
  xbutton("Single Run", "rerun()")
  xbutton("Reset", "reset()")
  xbutton("Quit", "quit()")
  xbutton("PSTH","psth()")
  xbutton("WriteFile","WriteFile()")
  xpanel(50,100)



  //----------------
  objref boxv
  boxv=new VBox()
  boxv.intercept(1)
  xpanel("")
  xlabel("ISI")
  g_ISI=new Graph()
  xpanel()
  
  xpanel("")
  xlabel("serial ISI")
  g_serialISI=new Graph()
  xpanel()

  

  boxv.intercept(0)
  boxv.map("CONTROL",600,50,-1,50)
  //----------------------------------------------------------------




 
  //****  procedure  "reset the resting leaky battery"   ****
  //----------------------------------------------------------- 
  
  v_init=-65
  celsius=38    
  reset()
  finitialize(v_init) 
  fcurrent()  // make sure that leak current balance the non-leak current
  print "total axon_conductance[uS]= ",G_pas
  print "tau_rest[ms]= ",1/(g_pas*1000)

  print "interval=[ms]", interval
  print "bin_ISI=[ms]",bin_ISI
  print "bin_psth=[ms]",bin_psth
  
  //-------- END reset -----------------------------------------
  //------------------------------------------------------------


  

Loading data, please wait...