// ***************** 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 -----------------------------------------
//------------------------------------------------------------
