diff --git a/neuron/.gitignore b/neuron/.gitignore new file mode 100644 index 0000000..58d9170 --- /dev/null +++ b/neuron/.gitignore @@ -0,0 +1,5 @@ +*.o +*.c +*.lnk +*.tmp +*.dll \ No newline at end of file diff --git a/neuron/Alpha.mod b/neuron/Alpha.mod new file mode 100644 index 0000000..093b90d --- /dev/null +++ b/neuron/Alpha.mod @@ -0,0 +1,68 @@ +COMMENT +** This mod file is based on the alpha synapse in NEURON, except for a ratio factor con in the NET_RECEIVE block.** +ENDCOMMENT + +NEURON { + POINT_PROCESS Alpha + RANGE tau1, tau2, e, i,con + NONSPECIFIC_CURRENT i + + RANGE g + GLOBAL total +} + +UNITS { + (nA) = (nanoamp) + (mV) = (millivolt) + (umho) = (micromho) +} + +PARAMETER { + tau1=.1 (ms) <1e-9,1e9> + tau2 = 10 (ms) <1e-9,1e9> + e=0 (mV) + con=10 +} + +ASSIGNED { + v (mV) + i (nA) + g (umho) + factor + total (umho) +} + +STATE { + A (umho) + B (umho) +} + +INITIAL { + LOCAL tp + total = 0 + if (tau1/tau2 > .9999) { + tau1 = .9999*tau2 + } + A = 0 + B = 0 + tp = (tau1*tau2)/(tau2 - tau1) * log(tau2/tau1) + factor = -exp(-tp/tau1) + exp(-tp/tau2) + factor = 1/factor +} + +BREAKPOINT { + SOLVE state METHOD cnexp + g = B - A + i = g*(v - e) +} + +DERIVATIVE state { + A' = -A/tau1 + B' = -B/tau2 +} + +NET_RECEIVE(weight (umho)) { + state_discontinuity(A, A + weight*factor*con) + state_discontinuity(B, B + weight*factor*con) + total = total+weight*con +} diff --git a/neuron/MSO_Zhouetal_2005.hoc b/neuron/MSO_Zhouetal_2005.hoc new file mode 100644 index 0000000..73c211d --- /dev/null +++ b/neuron/MSO_Zhouetal_2005.hoc @@ -0,0 +1,1244 @@ +// ***************** Coincidence Detector 2.0 ********************** + +// ------------ Introduction -------------------------------- +// 01/01/2005 +// This model simulates interaural time difference (ITD) dependent +// responses in the Medial Superio Olive (MSO) of the mammalian brainstem. +// Model configurations and simulation methods are described in +// Zhou, Carney, and Colburn (2005) J. Neuroscience, 25(12):3046-3058 + +// Ionic channels: ILTK, IHTK, Ina, Ih, Il from TypeII VCN cells +// in Rothman and Manis (2003c) J Neurophysiol 89:3097-113 + +// ------------ Model file structure ------------------------ +// Part I: Membrane and synapse setup +// Part II: Data processing and rate-ITD measurements +// Part III: User interface setup +//-------------------------------------------------------------- + + + + + +// ******* Codes start here **************************************** + +seed_x=xred("Please choose a random seed number <1,1000>", 1, 1, 1000) //string,default,min,max + + + +//------------------ Part I --------------------------------------- +// *** cell membrane and synapse descriptions *** +// This section defines the geometry and membrane properties of a bipolar MSO model. +// Excitatory and inhibitory synapses are implemented by the point-process tool in NEURON. +// Mod files associated with this section include: +// kHT_VCN2003.mod kLT_VCN2003.mod na_VCN2003.mod ih_VCN2003.mod Alpha.mod + + +//**** cell template ***** +//----------------------------- + + begintemplate Cell + + public soma, axon, dend + create soma,axon,dend[1] + + proc init() { + ndend=$1 + create soma,axon,dend[ndend] + + soma { + nseg = 1 + diam = 20 // um + L = 40 // um + Ra=200 //ohm.cm + cm=1 //uF/cm2 + + insert kHT_VCN2003 + insert kLT_VCN2003 + insert na_VCN2003 + insert ih_VCN2003 + + ek=-70 //mV + ena=55 + eh=-43 + + gkbar_kHT_VCN2003 = 0 //S/cm2 + gkbar_kLT_VCN2003=0 + gnabar_na_VCN2003= 0.1 + ghbar_ih_VCN2003=0 + gl_na_VCN2003=0.002 + + + + } + axon { + nseg = 51 + diam = 2 + L = 400 //um^2 + Ra=200 //ohm.cm + cm=1 //uF/cm2 + + insert kHT_VCN2003 + insert kLT_VCN2003 + insert na_VCN2003 + insert ih_VCN2003 + ek=-70 + ena=55 + eh=-43 + + // total G[nS]=gkbar*area*10 + + + gkbar_kHT_VCN2003 = 0.02 //S/cm2 + gkbar_kLT_VCN2003=0.03 + gnabar_na_VCN2003= 0.3 + ghbar_ih_VCN2003=0.0015 + gl_na_VCN2003=0.002 + } + + for i = 0, ndend-1 dend[i] { + L = 200 //um + nseg = 20 + diam= 3 // um + Ra=200 //ohm.cm + cm=1 //uF/cm2 + + insert pas + e_pas=-65 + g_pas=0.002 + } + + connect dend[0](0), soma(0) + connect dend[1](0), soma(1) + connect axon(0), dend[1](0.225) // 4.5/20 nseg(45um) + + } + + endtemplate Cell + //------- END cell template -------------------------------- + //---------------------------------------------------------- + + + load_file("nrngui.hoc") + + objectvar maincell + nmaindend =2 // two dendrites + maincell = new Cell(nmaindend) + access maincell.soma + + + + + + //**** Excitatory and inhibitory synapses **************** + //--------------------------------------------------------- + + //---- parameter description ---------------------- + //syn0[]: E alpha synapses on contralateral dendrite, dend0 + //syn1[]: E alpha synapses on ipsilateral dendrite, dend1 + //gen0[]: E stimuli to syn0 + //gen1[]: E stimuli to syn1 + //netcon0[]: input- E synapse connection on dend0 + //netcon1[]: input- E synapse connection on dend1 + + //syn_inhi0[]: I alpha synapse on the soma + //gen_inhi0[]: I stimuli to syn_inhi0 + //netcon_inhi0[]: input- I synapse connection on the soma + + + // ************ Excitatory synapses on two dendrites ******** + objectvar syn0[10], gen0[10],syn1[10], gen1[10] + + //----- add spike, virtual spikes + for i=0, 9{ + gen0[i] = new GaussStim(0.5) + gen0[i].interval=2 // [ms] input period + gen0[i].start=0 // [ms] arrival time of input spikes + gen0[i].number=10000 + gen0[i].factor=20 // phase-locking factor F + + gen1[i] = new GaussStim(0.5) + gen1[i].interval=2 + gen1[i].start=0 + gen1[i].number=10000 + gen1[i].factor=20 + + gen0[i].seed(seed_x+i) + gen1[i].seed(seed_x+10+i) + } + + //----- add EE expsyn at dend0 + maincell.dend[0] { //contra side + for i=0,9 { + n=20 //n=nseg doesn't change after nseg + syn0[i] = new Alpha(i*1/n+1/n/2) // scatter along the firsthalf dend + + syn0[i].tau1 =0.1 //ms + syn0[i].tau2 =0.1 //ms + syn0[i].con=11 //nS- mho + syn0[i].e =0 // mV + } + } + + maincell.dend[1] { //ipsi side + for i=0,9 { + n=20 //n=nseg + syn1[i] = new Alpha(i*1/n+1/n/2) + + syn1[i].tau1 =0.1 //ms may add some jitters + syn1[i].tau2 =0.1 //ms + syn1[i].con=11 //nS-mho + syn1[i].e =0 // mV + } + } + + + //--- Connect gen[] with syn[] + objref netcon0[10],netcon1[10] + + Rave_E_contra=240 // spikes/sec + Rave_E_ipsi=240 + + // x_thres controls the prob. of passing a spike + + x_thres_contra=1-Rave_E_contra*gen0[0].interval/1000 // prob(x>1-p)=p + x_thres_ipsi=1-Rave_E_ipsi*gen1[0].interval/1000 + + if (x_thres_contra<=0) {x_thres_contra=1e-4} + if (x_thres_ipsi<=0) {x_thres_ipsi=1e-4} + + + e_delay=0 // no synaptic delay + + for i=0,9 { + netcon0[i]=new NetCon(gen0[i], syn0[i],x_thres_contra,e_delay,0.001) + //thres delay gain[uS] + netcon1[i]=new NetCon(gen1[i], syn1[i],x_thres_ipsi,e_delay,0.001) + //thres delay gain[uS] + } + + //--------------------------------------------------------- + + + // ************ Inhibitory synapses on the soma ******** + + objectvar syn_inhi0[10], gen_inhi0[10], syn_inhi1[10], gen_inhi1[10] + //---------- gen_inhi ----------------- + for i=0,9{ + gen_inhi0[i]=new GaussStim(0.5) // contralateral inputs + gen_inhi0[i].interval=2 + gen_inhi0[i].start=0 + gen_inhi0[i].number=10000 + gen_inhi0[i].factor=10 + gen_inhi0[i].seed(seed_x+30+i) + + gen_inhi1[i]=new GaussStim(0.5) // ipsilateral inputs + gen_inhi1[i].interval=2 + gen_inhi1[i].start=0 + gen_inhi1[i].number=10000 + gen_inhi1[i].factor=10 + gen_inhi1[i].seed(seed_x+50+i) + } + + + + //----------- syn_inhi at soma----------------- + + maincell.soma { + for i=0,9{ + syn_inhi0[i] = new Alpha(0.5) + syn_inhi0[i].tau1 =0.1 //ms + syn_inhi0[i].tau2 =2 //ms + syn_inhi0[i].con=0 //umho + syn_inhi0[i].e =-70 // mV + } + + for i=0,9{ + syn_inhi1[i] = new Alpha(0.5) + syn_inhi1[i].tau1 =0.1 //ms + syn_inhi1[i].tau2 =2 //ms + syn_inhi1[i].con=0 //umho + syn_inhi1[i].e =-70 // mV + } + } + + //---------- Connect gen[] with syn[] --------------- + I_delay0=0 + I_delay1=0 // no synaptic delay + + objref netcon_inhi0[10],netcon_inhi1[10] + + Rave_I_contra=240 // spikes/sec + Rave_I_ipsi=240 + + x_thres_I_contra=1-Rave_I_contra*gen_inhi0[0].interval/1000 // prob(x>1-p)=p + x_thres_I_ipsi=1-Rave_I_ipsi*gen_inhi1[0].interval/1000 + + if (x_thres_I_contra<=0) {x_thres_I_contra=1e-4} + if (x_thres_I_ipsi<=0) {x_thres_I_ipsi=1e-4} + + + for i=0,9{ + netcon_inhi0[i]=new NetCon(gen_inhi0[i], syn_inhi0[i],x_thres_I_contra,I_delay0,0.001) + netcon_inhi1[i]=new NetCon(gen_inhi1[i], syn_inhi1[i],x_thres_I_ipsi,I_delay1,0.001) + } + +//-------- END synapse setup ------------------------------------- + + + + //**** procedure "reset the resting leaky battery" **** + //----------------------------------------------------------- + + v_init=-65 + celsius=38 + +proc reset_membrane() { + finitialize(v_init) + fcurrent() // make sure that leak current balance the non-leak current + maincell.soma { + print "gl_na_VCN2003=[S/cm^2]=",gl_na_VCN2003 + print "soma el=",el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003 + print "g_rest=S/cm2 ",g_rest1=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003 + print "tau_rest=",cm/(g_rest1*1000) + } + maincell.axon { print "axon el=",el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003 + print "g_rest=S/cm2 ",g_rest3=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003 + print "tau_rest=",cm/(g_rest3*1000) + } + maincell.dend[0] { + print "dendrite el=", e_pas=v_init + print " dend g_rest=S/cm2 ",g_rest2=g_pas + print "tau_rest=",cm/(g_rest2*1000) + } + maincell.dend[1] { + e_pas=v_init + } + + fcurrent() + +} + +reset_membrane() + +//-------- END reset_membrane ------------------------------ + +//------------------ END PART I ----------------------------------- + + + + + +//------------------ PART II --------------------------------------- +//*** data analysis and rate-ITD response *** +// In this section, synaptic inputs are integrated and membrane potentials +// are calculated at discrete time steps. The rate-ITD response at one ITD +// and/or at different ITDs (-1ms~ 1ms) are measured. + +// Model responses to ITDs are tested with and without inhibition. + +// The model discharges at one stimulus ITD condition are stored in two vectors: +// v_ISI for the event time of action potentials and v_voltage for the membrane potential +// at each integration time step. Statistics of discharge patterns are processed +// at the end of the program running, such as the post-stimulus time histogram, +// the inter-spike interval histograms, and firing rates. + +// The overall rate-ITD responses and firing patterns can also be saved into binary files for further analysis. + + + + + + //---- parameter and data vector descriptions ------------ + + // v_ISI: [ms] AP event time + // v_voltage: [mV] membrane potentials at each dt + // pre_E,pre_I: [ms] E/I input event time + + tstop=1000 // [ms] stimulus duration + dt=0.025 // [ms] simulation time step + bin=0.1 // [ms] binwidth of APs + + cvode.active(0) // constant integration time step + + + access maincell.axon + + objref apc,v_ISI,v_voltage + objref pre_E,pre_I_contra,pre_I_ipsi + + v_ISI=new Vector() + v_voltage=new Vector(tstop/dt+2,0) + + pre_E=new Vector() + pre_I_contra=new Vector() + pre_I_ipsi=new Vector() + + + //==== counting APs + AP_threshold=-10 //[mV] AP detection threshold + apc = new APCount(0.5) // measure action potenials at the midpoint of the axon + apc.thresh=AP_threshold + + //==== saving event times and membrane potentials + apc.record(v_ISI) //record event time + v_voltage.record(&maincell.axon.v(0.5),dt) + + + + //==== monitoring input events + netcon0[0].record(pre_E) //record input event time + netcon_inhi0[0].record(pre_I_contra) //record input event time + netcon_inhi1[0].record(pre_I_ipsi) //record input event time + + + //-------- END --------------------------------- + //----------------------------------------------- + + + + + //**** function "psth": post-stimulus time histogram **** + //----------------------------------------------- + //---- data vector description ------------ + //v_save: bin v_voltage for each run + //v_all: summation of all v_save over trials + //v_psth: PSTH vector + + objref v_save,v_all,v_psth + objref g_psth,x_psth + objref d, xtime + + bin_psth=bin + proc psth() { + length_psth=tstop/bin //length of period points + + v_psth=new Vector(v_all.size,0) + v_psth.add(v_all) + + x_psth= new Vector(length_psth,0) + x_psth.indgen(bin) + g_psth= new Graph() + g_psth.size(0,tstop,0,10) + g_psth.align(0.5,0.2) + g_psth.label(0.9,0.2,"time") + g_psth.label(0.1,0.9,"PSTH") + v_psth.plot(g_psth,x_psth,2,2) + + } + + //--- END psth () ------------------------------ + //----------------------------------------------- + + + + + //**** function "pth": period time histogram **** + //----------------------------------------------- + //---- data vector description ------------ + // v_save: bin v_voltage for each run + // v_pth: wrap v_all with length=gen.interval/bin and normalize + // vs_real: vector strength + // phase_real: mean phase + + // get_pth(): get v_save from each run + + + objref v_pth, g_pth, x_pth + + proc get_pth() { + 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 + } + + proc pth() { + x=0 + y=0 + vs_real=0 + angel_real=0 + length=gen0[0].interval/bin //length of period points + + N_fold=v_all.size/length + v_pth= new Vector(length,0) + x_pth= new Vector(length,0) // phase vector + x_pth.indgen(1/length) + + //==== wrapping over one period + for i=0,N_fold-1 { + for j=0,length-1 { + v_pth.x[j]=v_pth.x[j]+v_all.x[length*i+j] + } + } + v_pth.div(v_pth.sum) + + //==== calculate the VS and mean phase of responses + for i=0,length-1 { + x=x+v_pth.x[i]*cos(2*PI*i/length) + y=y+v_pth.x[i]*sin(2*PI*i/length) + } + vs_real=sqrt(x^2+y^2) //compute VS + print "VS=",vs_real + + angel_real=atan(y/x) + if(x<0 ) { angel_real=angel_real+PI + } else if (y<0 && x>0) { + angel_real=angel_real+2*PI + } + print "Phase=",angel_real/(2*PI) + + + g_pth.size(0,1,0,0.3) + g_pth.align(0.5,0.2) + //g_pth.label(0.9,0.1,"Phase") + g_pth.label(0.3,0.9,"Period Hist") + g_pth.begin() + v_pth.mark(g_pth,x_pth,"+",9,1) + v_pth.line(g_pth,x_pth,1,2) + + +} + + proc VS() { + vs=exp(-PI^2/(2*gen0[0].factor^2)) + } + + + //-------------- END pth () ------------------------- + //--------------------------------------------------- + + + + + + //**** function "ISI": interspike interval **** + //----------------------------------------------- + //---- data vector description ------------ + + // v_ISI: event times + // 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 statistics for each run + + // get_ISI(): get data for each run + + objref hist_ISI,xtime,g_ISI,hist + objref mean_T, sq_T, num_T + objref temp,t1, v1,v2 // operating vectors + + mean_ISI=0 + OVF=0 // overflow + EI=0 + + mean_ISI2=0 + std_ISI=0 + sq_sum=0 + CV_ISI=0 + + proc get_ISI() { + + 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 initial interval + print "mean ISI interval ", v2.mean + print "ISI std ",v2.stdev + + temp = new Vector(v2.size,0) + temp.add(v2) + temp.sub(v2.mean) // for calculate the square sum + + 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 "**********************************" + + hist=v2.histogram(0, topbin+OVFbin, bin_ISI) + + if(hist.size-hist_ISI.size==1) { hist.resize(hist_ISI.size)} + hist_ISI.add(hist) + + } + + proc ISI() { + hist_ISI.div(hist_ISI.sum) //normalize + OVF=hist_ISI.sum(topbin/bin_ISI, (topbin+OVFbin)/bin_ISI-1) + EI=hist_ISI.sum(0,1.5*T/bin_ISI-1) // the first whole interval + + hist_ISI.resize(topbin/bin_ISI+1) + + xtime=new Vector(hist_ISI.size,0) + xtime.indgen(bin_ISI) + + g_ISI.size(0,topbin,0,0.3) + g_ISI.align(0.5,0.5) + //g_ISI.label(0.9,0.2,"T") + g_ISI.label(0.2,0.9,"ISI") + g_ISI.beginline() + hist_ISI.mark(g_ISI,xtime,"o",6,1) + hist_ISI.line(g_ISI,xtime,1,2) + + + //==== 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_ISI2=t2/num_T.sum() //overall ISI mean + + t1=new Vector(N_run,0) + t1.add(mean_T) + t1.sub(mean_ISI2) + 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_ISI2>0) { + CV_ISI=std_ISI/mean_ISI2 //coefficient of variation + } else{ CV_ISI=1 + } + + } + + + //----------- END ISI() ------------------------- + //----------------------------------------------- + + + + + + //**** procedures: rate-ITD responses with and without inhibition **** + //---------------------------------------------------------------------- + //---- data Vecter description + // r_mean: discharge rates at different ITDs + // std_mean: standard deviations of rates at different ITDs + // vs_ITD, phase-ITD: response VS and mean phases at different ITDs + + + //---- procedure description ------------ + // reset() : initializing + // rerun() : main simulation function, repeat every tstop + // integrate(): integration function + + // Rate_ITD_EI(): rate-ITD response without bilateral excitation contralateral inhibition + // getcon(): upgrade new synapse values based on syn0[0], syn1[0], syn_inhi0[0], syn_inhi1[0]print "lamda_DC=[um]" + // getstim():upgrade new stimuli values based on gen0[0], gen1[0], gen_inhi0[0], gen_inhi1[0] + // savecon(): save synapse parameters + // changelevel(): update new input parameters + + objref v_AP, mean_T, sq_T, num_T // rates and spike intervals at one ITD over N_run + objref r_mean,std_mean,vs_ITD, phase_ITD // rates at all ITDs + objref T_ITD,std_ITD,CV_ITD // spike interval statistics at all ITDs + objref ITD,g_Rate,g_VS,data + + N_run=10 // ten trials + ITD_max=1 //[ms] maximum ITD tested + + + + //**** rate-ITD responses with excitation and inhibition ******************* + // ITD is defined as the arrival time difference between excitatory inputs on two dendrites. + //------------------------------------------------------------ + + proc Rate_ITD_EI() { + N_run=10 + ITD_step=0.05 // [ms] change to finer + half_ITD=ITD_max/ITD_step + N_ITD=2*half_ITD+1+2 // [-ITD_max :0.05:ITD_max C I ] + + ITD=new Vector(N_ITD,0) + ITD.indgen(ITD_step) + ITD.sub(ITD_max) + + r_mean=new Vector(N_ITD,0) + std_mean=new Vector(N_ITD,0) + vs_ITD=new Vector(N_ITD,0) + phase_ITD=new Vector(N_ITD,0) + + T_ITD=new Vector(N_ITD,0) + std_ITD=new Vector(N_ITD,0) + CV_ITD=new Vector(N_ITD,0) + + save_con() // save all the initial syn.con + + gen0[0].start=0 + gen1[0].start=0 + gen_inhi0[0].start=contra_inhi_start + gen_inhi1[0].start=ipsi_inhi_start + + + //==== test binaural EI responses + + for k=0,half_ITD { + gen0[0].start=(half_ITD-k)*ITD_step // contra delayed + gen_inhi0[0].start=contra_inhi_start+(half_ITD-k)*ITD_step // contra inhi delay + + gen1[0].start=0 + gen_inhi1[0].start=ipsi_inhi_start // ipsi inhi delay + rerun() + + r_mean.x[k]=v_AP.mean() + if(v_AP.size >1) {std_mean.x[k]=v_AP.stdev() + } else {std_mean.x[k]=0} + vs_ITD.x[k]=vs_real + phase_ITD.x[k]=angel_real/(2*PI) + + T_ITD.x[k]=mean_ISI2 + std_ITD.x[k]=std_ISI + CV_ITD.x[k]=CV_ISI + + } + + for k=half_ITD+1,2*half_ITD { + gen0[0].start=0 + gen_inhi0[0].start=contra_inhi_start + + gen1[0].start=(k-half_ITD)*ITD_step // ipsi delayed + gen_inhi1[0].start=ipsi_inhi_start +(k-half_ITD)*ITD_step // ipsi inhi delay + rerun() + + r_mean.x[k]=v_AP.mean() + if(v_AP.size >1) { std_mean.x[k]=v_AP.stdev() + } else {std_mean.x[k]=0} + vs_ITD.x[k]=vs_real + phase_ITD.x[k]=angel_real/(2*PI) + + T_ITD.x[k]=mean_ISI2 + std_ITD.x[k]=std_ISI + CV_ITD.x[k]=CV_ISI + } + + + //==== test contra monaural + k=2*half_ITD+1 + gen0[0].start=0 + gen1[0].start=0 + gen_inhi0[0].start=contra_inhi_start // contra inhi lead 0.1ms + gen_inhi1[0].start=ipsi_inhi_start + + syn0[0].con=e0 + syn1[0].con=0 + syn_inhi0[0].con=inhi0 + syn_inhi1[0].con=0 + + + rerun() + + r_mean.x[k]=v_AP.mean() + if(v_AP.size >1) {std_mean.x[k]=v_AP.stdev() + } else {std_mean.x[k]=0} + vs_ITD.x[k]=vs_real + phase_ITD.x[k]=angel_real/(2*PI) + + T_ITD.x[k]=mean_ISI2 + std_ITD.x[k]=std_ISI + CV_ITD.x[k]=CV_ISI + + + //==== test ipsi monaural + k=2*half_ITD+2 + gen0[0].start=0 + gen1[0].start=0 + gen_inhi0[0].start=contra_inhi_start + gen_inhi1[0].start=ipsi_inhi_start + + syn0[0].con=0 + syn1[0].con=e1 + syn_inhi0[0].con=0 + syn_inhi1[0].con=inhi1 + + rerun() + r_mean.x[k]=v_AP.mean() + if(v_AP.size >1) {std_mean.x[k]=v_AP.stdev() + } else {std_mean.x[k]=0} + vs_ITD.x[k]=vs_real + phase_ITD.x[k]=angel_real/(2*PI) + + + T_ITD.x[k]=mean_ISI2 + std_ITD.x[k]=std_ISI + CV_ITD.x[k]=CV_ISI + + print "############################################" + print"ITD= -1:0.05:1 C I " + print"R=" + r_mean.printf("% 3.2f") + print"\nR_std=" + std_mean.printf("% 2.2f") + print "\nVS=" + vs_ITD.printf("% 1.4f") + print "\nphase=" + phase_ITD.printf("% 1.4f") + + //print "T=" + //T_ITD.printf("% 2.4f") + //print "std_T=" + //T_ITD.printf("% 2.4f") + //print "CV_T=" + //CV_ITD.printf("% 2.4f") + + + + //==== draw the rate-ITD curve + + data=new VBox() + data.intercept(1) + xpanel("") + xlabel("Rate") + g_Rate=new Graph() + g_Rate.beginline() + g_Rate.size(-ITD_max, ITD_max,0, 50) + r_mean.plot(g_Rate,ITD,2,2) + r_mean.ploterr(g_Rate,ITD,std_mean,10,1,2) + xpanel() + + data.intercept(0) + data.map("DATA",300,50,-1,50) + + } + +//----------- END Rate_ITD_EI() ------------------------- +//------------------------------------------------------- + + + + + //----- reset function for all simulations + //------------------------------------------ + + proc reset() { + + + //==== Reset membrane, synapse, and inputs + reset_membrane() + getcon() + getstim() + changelevel() + + + //==== Reset all vec and var + N_run=10 + + v_voltage=new Vector(tstop/dt+2,0) + v_voltage.record(&maincell.axon.v(0.5),dt) // measure the midpoint of the axon + + v_AP=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) + + T=gen0[0].interval + bin_ISI=bin + topbin=int(10*T/bin_ISI)*bin_ISI // 10 intervals and 1 for OVF + OVFbin=int(1*T/bin_ISI)*bin_ISI + + mean_ISI2=0 + std_ISI=0 + sq_sum=0 + CV_ISI=0 + + mean_ISI=0 + OVF=0 + EI=0 + + hist_ISI=new Vector((topbin+OVFbin)/bin_ISI+1,0) + v_all=new Vector(tstop/dt+2,0) + v_all.rebin(v_all,bin/dt) + + + } + + //---------- end reset() ---------------- + //--------------------------------------- + + + + + + + //----- Main simulation function -------------------- + // Measure model response at one ITD condition + //--------------------------------------------------- + + + proc rerun() { + + reset() + + //********************************* + while (N_run>0) { + v_ISI=new Vector() + v_voltage=new Vector(tstop/dt+2,0) + apc.record(v_ISI) //record event time + v_voltage.record(&maincell.axon.v(0.5),dt) // measure APs at the midpoint of the axon + + finitialize(v_init) + fcurrent() + + integrate() + get_pth() + get_ISI() + v_AP.x[N_run-1]=apc.n + N_run=N_run-1 + } + + print "R=",v_AP.mean() + if(v_AP.size >1) {print "R_std=",v_AP.stdev() + print "R_sem=",v_AP.stderr() + } + print " " + N_run=10 + + g_ISI.erase_all() + g_pth.erase_all() + ISI() + pth() + + //print "ISI mean=[ms] ", mean_ISI/N_run + //print "T= ",mean_ISI2 + //print "T_std= ",std_ISI + //print "T_CV= ",CV_ISI + + //print "EI= ",EI + //print "OVF= ",OVF + print "----------- END -----------------------" +} + + + proc integrate() { + + while (t < tstop) { + fadvance() // advance solution + } + print "N_run=",N_run // show run time + print "spike number=",apc.n + print "-------------------------" + //print "Number_input_E=",pre_E.size + //print "Number_input_I_contra=",pre_I_contra.size + //print "Number_input_I_ipsi=",pre_I_ipsi.size + //print "*************" + + if(v_ISI.size<=2) { v_ISI=new Vector(3,0) } + L_ISI=v_ISI.size + //print "L_ISI=",L_ISI + L_PTH=v_voltage.size + //print "L_PTH=",L_PTH +} + + + + +//******* update synapse and input information **** +//--------------------------------------------------- + + proc save_con() { + e0=syn0[0].con + e1=syn1[0].con + inhi0=syn_inhi0[0].con + inhi1=syn_inhi1[0].con + + contra_inhi_start=gen_inhi0[0].start + ipsi_inhi_start=gen_inhi1[0].start + } + + + proc changelevel() { + + x_thres_contra=1-Rave_E_contra*gen0[0].interval/1000 // prob(x>1-p)=p + x_thres_ipsi=1-Rave_E_ipsi*gen1[0].interval/1000 + + x_thres_I_contra=1-Rave_I_contra*gen_inhi0[0].interval/1000 // prob(x>1-p)=p + x_thres_I_ipsi=1-Rave_I_ipsi*gen_inhi1[0].interval/1000 + + if (x_thres_contra<=0) {x_thres_contra=1e-4} + if (x_thres_ipsi<=0) {x_thres_ipsi=1e-4} + if (x_thres_I_contra<=0) {x_thres_I_contra=1e-4} + if (x_thres_I_ipsi<=0) {x_thres_I_ipsi=1e-4} + + + for i=0,9 { + netcon0[i].threshold=x_thres_contra + netcon1[i].threshold=x_thres_ipsi + + netcon_inhi0[i].threshold=x_thres_I_contra + netcon_inhi1[i].threshold=x_thres_I_ipsi + } + + } + + proc getcon() { + for i=0,9 { + syn0[i].con=syn0[0].con //umho + syn1[i].con=syn1[0].con //umho + + syn_inhi0[i].con=syn_inhi0[0].con //umho + syn_inhi1[i].con=syn_inhi1[0].con //umho + + } + print "*******************" + print "syn0 con =",syn0[1].con + print "syn1 con =",syn1[1].con + print "*******************" + print "syn_inhi0 con =",syn_inhi0[1].con + print "syn_inhi1 con =",syn_inhi1[1].con + print "--------------------------" + } + + proc getstim() { + //------- EE ------------- + for i=0, 9{ + gen0[i].interval=gen0[0].interval // + gen0[i].start=gen0[0].start // identify ITD + gen0[i].factor=gen0[0].factor + + gen1[i].interval=gen1[0].interval + gen1[i].start=gen1[0].start // identify ITD + gen1[i].factor=gen1[0].factor + } + + //------- II------------- + for i=0, 9{ + gen_inhi0[i].interval=gen_inhi0[0].interval // + gen_inhi0[i].start=gen_inhi0[0].start + gen_inhi0[i].factor=gen_inhi0[0].factor + + gen_inhi1[i].interval=gen_inhi1[0].interval + gen_inhi1[i].start=gen_inhi1[0].start + gen_inhi1[i].factor=gen_inhi1[0].factor + } + + print "E0 start =",gen0[1].start + print "E0 interval =",gen0[1].interval + print "E1 start =",gen1[1].start + print "E1 interval =",gen1[1].interval + print "*******************" + + print "I0 inhi start =",gen_inhi0[1].start + print "I0 inhi interval =",gen_inhi0[1].interval + print "I1 Inhi start =",gen_inhi1[1].start + print "I1 inhi interval =",gen_inhi1[1].interval + + print "*******************" + + } + //----------- END rerun() ----------------------------------- + //----------------------------------------------------------- + + + +//------------------ END PART II ---------------------------- + + + + + + + +//------------------ PART III --------------------------------------- +// the user interface for adjusting the synapse and input parameters +// and for the program control. + + + + //**** GUI "synapse and stimuli" **** + //--------------------------------------------- + + objref syn_con, gen_start + + syn_con = new HBox() + syn_con.intercept(1) + + xpanel("") + xlabel("E0 syn contra") + xvalue("G [nS]","syn0[0].con") + xlabel("E1 syn ipsi") + xvalue("G [nS]","syn1[0].con") + xpanel() + + xpanel("") + xlabel("I0 syn contra") + xvalue("G [nS]","syn_inhi0[0].con") + xlabel("I1 syn ipsi") + xvalue("G [nS]","syn_inhi1[0].con") + xpanel() + + xpanel("") + xlabel("Na at soma") + xvalue("na [S/cm2]","maincell.soma.gnabar_na_VCN2003") + + xpanel() + syn_con.intercept(0) + syn_con.map("SYN",0,150,-1,50) + + + //-------------------------------------------------------------------------- + + gen_start = new HBox() + gen_start.intercept(1) + + xpanel("") + xlabel("E0 input contra") + xvalue("delay [ms]","gen0[0].start") + xvalue("period","gen0[0].interval") + xvalue("synchrony factor","gen0[0].factor") + xvalue("rate [spikes/sec]","Rave_E_contra") + + xlabel("E1 input ipsi") + xvalue("delay [ms]","gen1[0].start") + xvalue("period [ms]","gen1[0].interval") + xvalue("synchrony factor","gen1[0].factor") + xvalue("rate [spikes/sec]","Rave_E_ipsi") + xpanel() + + xpanel("") + xlabel("I0 input contra") + xvalue("delay [ms]","gen_inhi0[0].start") + xvalue("period","gen_inhi0[0].interval") + xvalue("synchrony factor","gen_inhi0[0].factor") + xvalue("rate [spikes/sec]","Rave_I_contra") + + xlabel("I1 input ipsi") + xvalue("delay [ms]","gen_inhi1[0].start") + xvalue("period [ms]","gen_inhi1[0].interval") + xvalue("synchrony factor","gen_inhi1[0].factor") + xvalue("rate [spikes/sec]", "Rave_I_ipsi") + xpanel() + + gen_start.intercept(0) + gen_start.map("GEN",0,350,-1,50) + + +//-------------------------------------------------------------------------- + + + objref boxv + boxv=new VBox() + boxv.intercept(1) + xpanel("") + xbutton("calculate VS of E0 input","VS()") + xvalue("vs") + xbutton("Reset","reset()") + xbutton("Single Run", "rerun()") + xbutton("Rate-ITD response with EE/II inputs","Rate_ITD_EI()") + xbutton("Quit", "quit()") + //xbutton("PSTH","psth()") + g_ISI=new Graph() + g_pth=new Graph() + xpanel() + boxv.intercept(0) + boxv.map("CONTROL",600,50,-1,50) + //---------------------------------------------------------------- + + + + + + + //**** Electrotonic properties of the model ***** + //---------------------------------------------------------------- + + proc pas_membrane() { + + + maincell.axon { + print "======================================" + area3=PI*diam*L + G_rest3=area3*g_rest3*10 //nS + Rm_rest3=1000/G_rest3 //Mohm + Ri_rest3=4*Ra*L/(100*PI*diam*diam) //Mohm + + R_inf3=2*sqrt(Ra/(gl_na_VCN2003*diam^3))/PI //Mohm + + lambda=100*sqrt(diam/(4*Ra*gl_na_VCN2003)) + R_axon=R_inf3/tanh(L/lambda) + + print "axon_area=[um^2]",area3 + print "lambda_DC=[um]",lambda + print "total axon_conductance=[nS]",G_rest3 + + print "seal-end input resistance=[Mohn]",R_axon + print "axon el=[mV]",el_na_VCN2003 + print "g_rest=[S/cm2]",g_rest3 + print "axon tau_rest=[ms]",cm/(g_rest3*1000) + print "======================" + } + + maincell.dend[0] { + area2=PI*diam*L + G_rest2=area2*g_rest2*10 //nS + Rm_rest2=1000/G_rest2 //Mohm + Ri_rest2=4*Ra*L/(100*PI*diam*diam) //Mohm + + R_inf2=2*sqrt(Ra/(g_pas*diam^3))/PI //Mohm + + lambda=100*sqrt(diam/(4*Ra*g_rest2)) + R_dend=R_inf2/tanh(L/lambda) + + print "each dend_area=[um^2]",area2 + print "lambda_DC=[um]",lambda + print "each dend_conductance=[nS]",G_rest2 + + print "seal-end input resistance=[Mohn]",R_dend + print "dend el=[mV]",e_pas + print "g_rest=[S/cm2]",g_rest2 + print "dend tau_rest=[ms]",cm/(g_rest2*1000) + print "======================" + } + + + maincell.soma { + area1=PI*diam*L + G_rest1=area1*g_rest1*10 //nS + Rm_rest1=1000/G_rest1 //Mohm + Ri_rest1=4*Ra*L/(100*PI*diam*diam) //Mohm + R_inf1=2*sqrt(Ra/(g_rest1*diam^3))/PI //Mohm + + + + print "soma_area=[um^2]",area1 + print "lamda_DC=[um]",lambda=100*sqrt(diam/(4*Ra*gl_na_VCN2003)) + + print "total soma_conductance=[nS]",G_rest1 + print "total membrane resistance=[Mohm]", Rm_rest1 + print "soma el=[mV]",el_na_VCN2003 + print "g_rest=[S/cm2] ", g_rest1 + print "soma tau_rest=[ms]",cm/(g_rest1*1000) + print "======================" + + R_soma=R_inf1*(R_dend/2+R_inf1*tanh(L/lambda))/(R_inf1+R_dend*tanh(L/lambda)/2) + R_dend_soma=R_inf1*(R_dend+R_inf1*tanh(L/lambda))/(R_inf1+R_dend*tanh(L/lambda)) + } + + + R_input=R_soma + + + + +print "input resistance=[Mohm]", R_input +print "rho=",R_dend_soma/R_inf2 + +} // end of pas_membrane() + + pas_membrane() + + //------- END Part III --------------------------------------------------------- + diff --git a/neuron/gaussstim.mod b/neuron/gaussstim.mod new file mode 100644 index 0000000..be397b6 --- /dev/null +++ b/neuron/gaussstim.mod @@ -0,0 +1,223 @@ +COMMENT +This mod file is based on netstim.mod in NEURON +modified by Yi Zhou to generate Gaussian interval distributions +ENDCOMMENT + +NEURON { + POINT_PROCESS GaussStim + RANGE x + RANGE interval, number, start, factor, refrac + RANGE N_backward, N_forward, N_normal, N_total + RANGE rnd + RANGE rand +} + +PARAMETER { + interval = 10 (ms) <1e-9,1e9> : ideal time between spikes (ms) + number = 10000 <0,1e9> : maximum number of spikes + start = 10 (ms) : delay before first spiketime + factor = 4 : fraction of mean to compute std (= mean / factor) + refrac = 0.5 <0,3> : absolute refractory period (upper limit of firing = 1 kHz) + : refrac >= dt, otherwise x can't be reset to 0 +} + +ASSIGNED { + x : if x > 0, cell is in refractory period + event (ms) : spiketime of next event in absolute time + on + end (ms) + m (ms) : mean of Gaussian + diff (ms) + N_forward : swap spike whose value exceed forwardly one interval + N_backward : swap spike whose value exceed backwardly one interval + N_normal + N_total + rnd + rand +} + +PROCEDURE seed(x) { + set_seed(x) +} + +INITIAL { + on = 0 + x = 0 + diff = 0 + m = interval / 2 : each T has normal distribution N(interval/2, interval/2/factor) + N_forward = 0 + N_backward = 0 + N_normal = 0 + N_total = 0 + + if (start >= 0 && number > 0) { + :: generate the first spiketime + event = start + invl(m) + + :: make sure first spiketime is not before 0 + if (event < 0) { + event = 0 + } + + :: skip directly to next event + net_send(event, 3) + } +} + +PROCEDURE init_sequence(t(ms)) { +:: Turn on spike generation if n_spikes_max > 0 +:: - "on" is set to 1 +:: - "event" is set to "t" (input of this call) +:: - "end" is set to "interval * n_spikes_max" + + if (number > 0) { + on = 1 + event = t + end = t + 1e-6 + interval * (number-1) + } +} + +PROCEDURE event_time() { +:: Generate the next spiketime +:: - call invl which calls normrand() + + :: define local variables + LOCAL diff2, T + + :: BAD diff gets overwritten in invl :facepalm: + diff2 = diff + + :: if n_spikes_max is larger than 0 + if (number > 0) { + :: generate a new interval until next spiketime + T = invl(m) + + :: BAD useless + rnd = T + + :: if new spiketime T is equal to an existing spiketime + if (T == 0 && diff2 == 0) { + T = T + dt + } + + :: compute absolute time of next event (relative to 0 ms) + event = event + diff2 + T + + :: augment total counter + N_total = N_total + 1 + } + + :: if next event is after the end, turn off spike generation + if (event > end) { + on = 0 :: BAD not the role of this function + } +} + +FUNCTION invl(mean (ms)) (ms) { +:: Returns a new interval invl +:: Warning: also sets the global variable diff + + :: define local variable + LOCAL std + + :: compute standard deviation + if (mean <= 0.) { + mean = .01 (ms) + } + std = mean / factor + + :: generate new interval invl + invl = normrand(mean, std) + + :: make sure invl is within 0 and max interval + if (invl >= interval) { + :: decrease invl by n*interval + invl = fmod(invl, interval) + :: augment forward counter + N_forward = N_forward + 1 + } else if (invl < 0) { + :: increase invl by n*interval + invl = fmod(invl, interval) + interval + :: augment backward counter + N_backward = N_backward + 1 + } else { + :: augment normal counter + N_normal = N_normal + 1 + } + + :: overwrite global variable diff + diff = interval - invl +} + +NET_RECEIVE(w) { + :: EXTERNAL CALL + if (flag == 0) { + if (w > 0 && on == 0) { + :: turn on spike generation + init_sequence(t) + net_send(0, 1) + } else if (w < 0 && on == 1) { + :: turn off spike generation + on = 0 + } + } + + :: MAIN CALL + if (flag == 1 && on == 1) { + if (x == 0) { :: time to spike (refractory period has passed) + :: set x to random uniform number (useless, just needs to > 0) + rand = scop_random() : [0, 1] uniform + x = rand + + :: output event to synapse + net_event(t) + + :: generate next spiketime + event_time() + + :: if next spiketime is between refrac and refrac+dt + if (event-t >= refrac && event-t <= refrac+dt) { + :: move next spiketime up by one dt (just after refrac) + event = event + dt + } + + :: if spike generation is turned on + if (on == 1) { + :: skip directly to next event + net_send(event - t, 1) + } + + :: refractory period call + net_send(refrac, 2) + } else if (x != 0) { :: time to spike (within refractory period) + :: output event to synapse + net_event(t) :: BAD even still emitted even though refractory period, should be deleted + + :: generate next spiketime + event_time() : although this spike will be ignored, still call next spike (independent cycles) + + :: if spike generation is turned on + if (on == 1) { + :: skip directly to next event + net_send(event - t, 1) + } + } + } + + :: REFRACTORY PERIOD CALL (END OF REFRACTORY PERIOD) + if (flag == 2) { + x = 0 + } + + :: INITIAL CALL + if (flag == 3) { + :: if spike generation is not yet turned on + if (on == 0) { + :: turn on spike generation + init_sequence(t) + + :: call net_receive directly + net_send(0, 1) + } + } +} diff --git a/neuron/ih_VCN2003.mod b/neuron/ih_VCN2003.mod new file mode 100644 index 0000000..3c6d52f --- /dev/null +++ b/neuron/ih_VCN2003.mod @@ -0,0 +1,80 @@ +TITLE Ih channels in VCN auditory neurons +: ih=ghmax*r*(v-eh) +: based on Rothman and Manis (2003c) +: Modifications by Yi Zhou for an MSO model + +INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)} + +NEURON { + SUFFIX ih_VCN2003 + USEION h READ eh WRITE ih VALENCE 1 + RANGE ghbar + RANGE r_inf,tau_r,r_exp + RANGE ih,gh + +} + + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) +} + +PARAMETER { + ghbar = 0.002 (mho/cm2) + eh=-43 (mV) + celsius =22 (degC) + dt (ms) + v (mV) + +} + +STATE { + r +} + +ASSIGNED { + gh (mho/cm2) + ih (mA/cm2) + r_inf + tau_r + r_exp + tadj +} + + +BREAKPOINT { + SOLVE states + gh=ghbar *r + ih = gh*(v-eh) +} + + +PROCEDURE states() { : this discretized form is more stable + evaluate_fct(v) + r = r + r_exp * (r_inf - r) + VERBATIM + return 0; + ENDVERBATIM +} + +UNITSOFF +INITIAL { +: +: Q10 was assumed to be 3 for both currents +: + tadj = 3.0 ^ ((celsius-22)/ 10 ) + evaluate_fct(v) + r= r_inf + } + +PROCEDURE evaluate_fct(v(mV)) { + + tau_r = (100000/(237*exp((v+60)/12)+17*exp(-(v+60)/14))+25)/ tadj + r_inf = 1/(1+exp((v+76)/7)) + + r_exp = 1 - exp(-dt/tau_r) + +} + +UNITSON diff --git a/neuron/kHT_VCN2003.mod b/neuron/kHT_VCN2003.mod new file mode 100644 index 0000000..6219169 --- /dev/null +++ b/neuron/kHT_VCN2003.mod @@ -0,0 +1,95 @@ +TITLE high threshold potassium channels in VCN auditory neurons + +: k_HT=ght*(rr*n^2+(1-rr)*p)*(v-Ek) +: based on Rothman and Manis (2003c) +: +: Modifications by Yi Zhou for an MSO model + +INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)} + +NEURON { + SUFFIX kHT_VCN2003 + USEION k READ ek WRITE ik + RANGE gkbar + RANGE n_inf,p_inf + RANGE tau_n,tau_p + RANGE n_exp,p_exp + RANGE ik,gk + +} + + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) +} + +PARAMETER { + gkbar = 0.03 (mho/cm2) + ek=-70 (mV) + celsius =22 (degC) + dt (ms) + v (mV) + +} + +STATE { + n p rr +} + +ASSIGNED { + gk(mho/cm2) + ik (mA/cm2) + n_inf + p_inf + tau_n + tau_p + n_exp + p_exp + tadj +} + + +BREAKPOINT { + SOLVE states + gk=gkbar *(rr*n^2+(1-rr)*p) + ik = gk*(v-ek) +} + + + +PROCEDURE states() { : this discretized form is more stable + evaluate_fct(v) + n = n + n_exp * (n_inf - n) + p = p + p_exp * (p_inf - p) + VERBATIM + return 0; + ENDVERBATIM +} + +UNITSOFF +INITIAL { +: +: Q10 was assumed to be 3 for both currents +: + tadj = 3.0 ^ ((celsius-22)/ 10 ) + evaluate_fct(v) + n= n_inf + p= p_inf + rr=0.85 + } + +PROCEDURE evaluate_fct(v(mV)) { + + tau_n = (100/(11*exp((v+60)/24)+21*exp(-(v+60)/23))+0.7)/ tadj + n_inf = 1/(1+exp(-(v+15)/5))^2 + + tau_p = (100/(4*exp((v+60)/32)+5*exp(-(v+60)/22))+5)/ tadj + p_inf = 1/(1+exp(-(v+23)/6)) + + n_exp = 1 - exp(-dt/tau_n) + p_exp = 1 - exp(-dt/tau_p) + +} + +UNITSON diff --git a/neuron/kLT_VCN2003.mod b/neuron/kLT_VCN2003.mod new file mode 100644 index 0000000..4018e91 --- /dev/null +++ b/neuron/kLT_VCN2003.mod @@ -0,0 +1,93 @@ +TITLE low threshold potassium channels in VCN auditory neurons +: k_LT=glt*w^4*z*(v-Ek) +: based on Rothman and Manis (2003c) +: Modifications by Yi Zhou for an MSO model + +INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)} + +NEURON { + SUFFIX kLT_VCN2003 + USEION k READ ek WRITE ik + RANGE gkbar + RANGE w_inf,z_inf + RANGE tau_w,tau_z + RANGE w_exp,z_exp + RANGE ik,gk + +} + + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) +} + +PARAMETER { + gkbar = 0.04 (mho/cm2) + ek=-70 (mV) + celsius =22 (degC) + dt (ms) + v (mV) + +} + +STATE { + w z +} + +ASSIGNED { + gk(mho/cm2) + ik(mA/cm2) + w_inf + z_inf + tau_w + tau_z + w_exp + z_exp + tadj +} + + +BREAKPOINT { + SOLVE states + gk=gkbar *w^4*z + ik = gk*(v-ek) +} + + + +PROCEDURE states() { : this discretized form is more stable + evaluate_fct(v) + w = w + w_exp * (w_inf - w) + z = z + z_exp * (z_inf - z) + VERBATIM + return 0; + ENDVERBATIM +} + +UNITSOFF +INITIAL { +: +: Q10 was assumed to be 3 for both currents +: + tadj = 3.0 ^ ((celsius-22)/ 10 ) + evaluate_fct(v) + w= w_inf + z= z_inf + } + +PROCEDURE evaluate_fct(v(mV)) {LOCAL h + + tau_w = (100/(6*exp((v+60)/6)+16*exp(-(v+60)/45))+1.5)/ tadj + w_inf = 1/(1+exp(-(v+48)/6))^0.25 + + tau_z = (1000/(exp((v+60)/20)+exp(-(v+60)/8))+50)/ tadj + h=0.5 + z_inf = (1-h)/(1+exp((v+71)/10))+h + + w_exp = 1 - exp(-dt/tau_w) + z_exp = 1 - exp(-dt/tau_z) + +} + +UNITSON diff --git a/neuron/mosinit.hoc b/neuron/mosinit.hoc new file mode 100644 index 0000000..9e7a354 --- /dev/null +++ b/neuron/mosinit.hoc @@ -0,0 +1,34 @@ +load_file("nrngui.hoc") + +strdef tstr + +xpanel("Figures for Zhou et al 2005") +xlabel("Press a button to create the figure") +xradiobutton("Main simulation", "restart(\"MSO_Zhouetal_2005\")") +xradiobutton("Additional simulation", "restart(\"test_inputs\")") +xlabel("After you click on one of the above buttons you will be asked") +xlabel("to enter a number which will be used to seed the random number") +xlabel("generator. After you have entered the number please press the") +xlabel("\"Single Run\" button which will appear in a Control window") +xlabel("Note: the main simulation takes 3 1/4 minutes on a 2Ghz machine") +xlabel("and the additional simulation takes less than 15 seconds.") +xpanel(5,100) + +pwmcnt = PWManager[0].count // the initial gui should not be dismissed + +proc restart() {local i + objref graphItem, save_window_ + for i=0, n_graph_lists-1 { + graphList[i].remove_all() + } + flush_list.remove_all() + fast_flush_list.remove_all() + doNotify() + for (i= PWManager[0].count-1; i >= pwmcnt; i -= 1) { + PWManager[0].close(i) + doNotify() + } + + sprint(tstr, "%s.hoc", $s1) + load_file(1, tstr) +} diff --git a/neuron/msomodel_base.hoc b/neuron/msomodel_base.hoc new file mode 100644 index 0000000..71829e5 --- /dev/null +++ b/neuron/msomodel_base.hoc @@ -0,0 +1,378 @@ +// Init measurements +objref vm_s, vm_d0p, vm_d0m, vm_d0d, vm_d1p, vm_d1m, vm_d1d, vm_ap, vm_ad +vm_s = new Vector(tstop/dt+2,0) +vm_d0p = new Vector(tstop/dt+2,0) +vm_d0m = new Vector(tstop/dt+2,0) +vm_d0d = new Vector(tstop/dt+2,0) +vm_d1p = new Vector(tstop/dt+2,0) +vm_d1m = new Vector(tstop/dt+2,0) +vm_d1d = new Vector(tstop/dt+2,0) +vm_ap = new Vector(tstop/dt+2,0) +vm_ad = new Vector(tstop/dt+2,0) +vm_s.record(&maincell.soma.v(0.5), dt) +vm_d0p.record(&maincell.dend[0].v(0.1), dt) +vm_d0m.record(&maincell.dend[0].v(0.5), dt) +vm_d0d.record(&maincell.dend[0].v(0.9), dt) +vm_d1p.record(&maincell.dend[1].v(0.1), dt) +vm_d1m.record(&maincell.dend[1].v(0.5), dt) +vm_d1d.record(&maincell.dend[1].v(0.9), dt) +vm_ap.record(&maincell.axon.v(0.1), dt) +vm_ad.record(&maincell.axon.v(0.5), dt) + + +// Variables left over from obsolete functions +// TODO find obsolete variables +bin = 0.1 // [ms] binwidth of APs +cvode.active(0) // constant integration time step + +access maincell.axon + +objref apc, v_ISI, v_voltage +objref pre_E, pre_I_contra, pre_I_ipsi + +v_ISI = new Vector() +v_voltage = new Vector(tstop/dt+2,0) + +pre_E = new Vector() +pre_I_contra = new Vector() +pre_I_ipsi = new Vector() + +//==== counting APs +AP_threshold = -10 //[mV] AP detection threshold +apc = new APCount(0.5) // measure action potenials at the midpoint of the axon +apc.thresh = AP_threshold + +//==== saving event times and membrane potentials +apc.record(v_ISI) //record event time +v_voltage.record(&maincell.axon.v(0.5), dt) + +//==== monitoring input events +netcon0[0].record(pre_E) //record input event time +netcon_inhi0[0].record(pre_I_contra) //record input event time +netcon_inhi1[0].record(pre_I_ipsi) //record input event time + +objref v_save,v_all,v_psth +objref g_psth,x_psth +objref d, xtime + +bin_psth=bin + +objref v_pth, g_pth, x_pth + +objref hist_ISI,xtime,g_ISI,hist +objref mean_T, sq_T, num_T +objref temp,t1, v1,v2 // operating vectors + +mean_ISI=0 +OVF=0 // overflow +EI=0 + +mean_ISI2=0 +std_ISI=0 +sq_sum=0 +CV_ISI=0 + +objref v_AP, mean_T, sq_T, num_T // rates and spike intervals at one ITD over N_run +objref r_mean,std_mean,vs_ITD, phase_ITD // rates at all ITDs +objref T_ITD,std_ITD,CV_ITD // spike interval statistics at all ITDs +objref ITD,g_Rate,g_VS,data + +N_run = n_runs_total +ITD_max=1 //[ms] maximum ITD tested + + +// Functions +proc reset_membrane() { + finitialize(v_init) + fcurrent() + + maincell.soma { + gl_na_VCN2003 + el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003 + g_rest1=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003 + cm/(g_rest1*1000) + } + maincell.axon { + el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003 + g_rest3=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003 + cm/(g_rest3*1000) + } + maincell.dend[0] { + e_pas=v_init + g_rest2=g_pas + cm/(g_rest2*1000) + } + maincell.dend[1] { + e_pas=v_init + } + + fcurrent() +} +// proc reset_membrane() { +// finitialize(v_init) +// fcurrent() // make sure that leak current balance the non-leak current +// maincell.soma { +// print "gl_na_VCN2003=[S/cm^2]=", gl_na_VCN2003 +// print "soma el=", el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003 +// print "g_rest=S/cm2 ", g_rest1=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003 +// print "tau_rest=", cm/(g_rest1*1000) +// } +// maincell.axon { +// print "axon el=", el_na_VCN2003=(ina+ik+ih+gl_na_VCN2003*v_init)/gl_na_VCN2003 +// print "g_rest=S/cm2 ", g_rest3=gl_na_VCN2003+gk_kLT_VCN2003+gk_kHT_VCN2003+gna_na_VCN2003+gh_ih_VCN2003 +// print "tau_rest=", cm/(g_rest3*1000) +// } +// maincell.dend[0] { +// print "dendrite el=", e_pas=v_init +// print " dend g_rest=S/cm2 ", g_rest2=g_pas +// print "tau_rest=", cm/(g_rest2*1000) +// } +// maincell.dend[1] { +// e_pas=v_init +// } + +// fcurrent() +// } + +reset_membrane() + +proc reset() { + reset_membrane() + + //==== Reset all vec and var + N_run = n_runs_total + + v_voltage=new Vector(tstop/dt+2,0) + v_voltage.record(&maincell.axon.v(0.5),dt) // measure the midpoint of the axon + + vm_s = new Vector(tstop/dt+2,0) + vm_d0p = new Vector(tstop/dt+2,0) + vm_d0m = new Vector(tstop/dt+2,0) + vm_d0d = new Vector(tstop/dt+2,0) + vm_d1p = new Vector(tstop/dt+2,0) + vm_d1m = new Vector(tstop/dt+2,0) + vm_d1d = new Vector(tstop/dt+2,0) + vm_ap = new Vector(tstop/dt+2,0) + vm_ad = new Vector(tstop/dt+2,0) + vm_s.record(&maincell.soma.v(0.5), dt) + vm_d0p.record(&maincell.dend[0].v(0.1), dt) + vm_d0m.record(&maincell.dend[0].v(0.5), dt) + vm_d0d.record(&maincell.dend[0].v(0.9), dt) + vm_d1p.record(&maincell.dend[1].v(0.1), dt) + vm_d1m.record(&maincell.dend[1].v(0.5), dt) + vm_d1d.record(&maincell.dend[1].v(0.9), dt) + vm_ap.record(&maincell.axon.v(0.1), dt) + vm_ad.record(&maincell.axon.v(0.5), dt) + + v_AP=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) +} + +objref fl +objref flPG +objref spkt_contra +objref spkt_ipsi +strdef flname +strdef cmd +fl = new File() +flPG = new File() + +ispkt_contra = 0 +ispkt_ipsi = 0 +spkt_contra = new Vector() +spkt_ipsi = new Vector() + +proc rerun() { + reset() + + //********************************* + while (N_run>0) { + v_ISI=new Vector() + v_voltage=new Vector(tstop/dt+2,0) + apc.record(v_ISI) //record event time + v_voltage.record(&maincell.axon.v(0.5),dt) // measure APs at the midpoint of the axon + + vm_s = new Vector(tstop/dt+2,0) + vm_d0p = new Vector(tstop/dt+2,0) + vm_d0m = new Vector(tstop/dt+2,0) + vm_d0d = new Vector(tstop/dt+2,0) + vm_d1p = new Vector(tstop/dt+2,0) + vm_d1m = new Vector(tstop/dt+2,0) + vm_d1d = new Vector(tstop/dt+2,0) + vm_ap = new Vector(tstop/dt+2,0) + vm_ad = new Vector(tstop/dt+2,0) + vm_s.record(&maincell.soma.v(0.5), dt) + vm_d0p.record(&maincell.dend[0].v(0.1), dt) + vm_d0m.record(&maincell.dend[0].v(0.5), dt) + vm_d0d.record(&maincell.dend[0].v(0.9), dt) + vm_d1p.record(&maincell.dend[1].v(0.1), dt) + vm_d1m.record(&maincell.dend[1].v(0.5), dt) + vm_d1d.record(&maincell.dend[1].v(0.9), dt) + vm_ap.record(&maincell.axon.v(0.1), dt) + vm_ad.record(&maincell.axon.v(0.5), dt) + + flPG.ropen("eventtimes_contra.txt") + spkt_contra.scanf(flPG) + flPG.close() + flPG.ropen("eventtimes_ipsi.txt") + spkt_ipsi.scanf(flPG) + flPG.close() + + for i = 0,9 { + gen0[i].timeToNextEvt = spkt_contra.x[0] - dt + gen1[i].timeToNextEvt = spkt_ipsi.x[0] - dt + } + + finitialize(v_init) + fcurrent() + + integrate() + // get_pth() + // get_ISI() + v_AP.x[N_run-1] = apc.n + N_run = N_run-1 + + // Export to temporary files + sprint(flname, "../data/_direct/run_%05i__vm_c.txt.tmp", iter) + fl.wopen(flname) + vm_s.printf(fl) + fl.close() + + sprint(flname, "../data/_direct/run_%05i__vm_d0p.txt.tmp", iter) + fl.wopen(flname) + vm_d0p.printf(fl) + fl.close() + + sprint(flname, "../data/_direct/run_%05i__vm_d0m.txt.tmp", iter) + fl.wopen(flname) + vm_d0m.printf(fl) + fl.close() + + sprint(flname, "../data/_direct/run_%05i__vm_d0d.txt.tmp", iter) + fl.wopen(flname) + vm_d0d.printf(fl) + fl.close() + + sprint(flname, "../data/_direct/run_%05i__vm_d1p.txt.tmp", iter) + fl.wopen(flname) + vm_d1p.printf(fl) + fl.close() + + sprint(flname, "../data/_direct/run_%05i__vm_d1m.txt.tmp", iter) + fl.wopen(flname) + vm_d1m.printf(fl) + fl.close() + + sprint(flname, "../data/_direct/run_%05i__vm_d1d.txt.tmp", iter) + fl.wopen(flname) + vm_d1d.printf(fl) + fl.close() + + sprint(flname, "../data/_direct/run_%05i__vm_ap.txt.tmp", iter) + fl.wopen(flname) + vm_ap.printf(fl) + fl.close() + + sprint(flname, "../data/_direct/run_%05i__vm_ad.txt.tmp", iter) + fl.wopen(flname) + vm_ad.printf(fl) + fl.close() + + // Rename files from temporary name to destined name + sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_c.txt.tmp run_%05i__vm_c.txt", iter, iter) + print "Cmd: ",cmd + system(cmd) + sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d0p.txt.tmp run_%05i__vm_d0p.txt", iter, iter) + print "Cmd: ",cmd + system(cmd) + sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d0m.txt.tmp run_%05i__vm_d0m.txt", iter, iter) + system(cmd) + sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d0d.txt.tmp run_%05i__vm_d0d.txt", iter, iter) + system(cmd) + sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d1p.txt.tmp run_%05i__vm_d1p.txt", iter, iter) + system(cmd) + sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d1m.txt.tmp run_%05i__vm_d1m.txt", iter, iter) + system(cmd) + sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_d1d.txt.tmp run_%05i__vm_d1d.txt", iter, iter) + system(cmd) + sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_ap.txt.tmp run_%05i__vm_ap.txt", iter, iter) + system(cmd) + sprint(cmd, "rename ..\\data\\_direct\\run_%05i__vm_ad.txt.tmp run_%05i__vm_ad.txt", iter, iter) + system(cmd) + } + + N_run = n_runs_total + +} + +proc integrate() { + // while the simulation has not yet ended + while (t < tstop) { + // CONTRA + // if t is earlier than the last spiketime + if (t < spkt_contra.x[spkt_contra.size()-1]) { + // if t is right after a spiketime has been issued, issue the next one + if (t >= spkt_contra.x[ispkt_contra-1]+dt && t < spkt_contra.x[ispkt_contra-1]+dt*2) { + if (ispkt_contra == spkt_contra.size()-1) { + for i = 0,9 { + gen0[i].timeToNextEvt = 0 + } + } else { + for i = 0,9 { + gen0[i].timeToNextEvt = spkt_contra.x[ispkt_contra+1] - spkt_contra.x[ispkt_contra] + } + } + + ispkt_contra = ispkt_contra + 1 + } else if (t >= dt && t < dt*2) { // if t is right the beginning of the simulation + for i = 0,9 { + gen0[i].timeToNextEvt = spkt_contra.x[ispkt_contra+1] - spkt_contra.x[ispkt_contra] - dt + } + ispkt_contra = ispkt_contra + 1 + } + } + + // IPSI + // if t is earlier than the last spiketime + if (t < spkt_ipsi.x[spkt_ipsi.size()-1]) { + // if t is right after a spiketime has been issued, issue the next one + if (t >= spkt_ipsi.x[ispkt_ipsi-1]+dt && t < spkt_ipsi.x[ispkt_ipsi-1]+dt*2) { + if (ispkt_ipsi == spkt_ipsi.size()-1) { + for i = 0,9 { + gen1[i].timeToNextEvt = 0 + } + } else { + for i = 0,9 { + gen1[i].timeToNextEvt = spkt_ipsi.x[ispkt_ipsi+1] - spkt_ipsi.x[ispkt_ipsi] + } + } + + ispkt_ipsi = ispkt_ipsi + 1 + } else if (t >= dt && t < dt*2) { // if t is right the beginning of the simulation + for i = 0,9 { + gen1[i].timeToNextEvt = spkt_ipsi.x[ispkt_ipsi+1] - spkt_ipsi.x[ispkt_ipsi] - dt + } + ispkt_ipsi = ispkt_ipsi + 1 + } + } + + fadvance() // advance solution + } + print "N_run=",N_run // show run time + print "spike number=",apc.n + // print "-------------------------" + //print "Number_input_E=",pre_E.size + //print "Number_input_I_contra=",pre_I_contra.size + //print "Number_input_I_ipsi=",pre_I_ipsi.size + //print "*************" + + if(v_ISI.size<=2) { v_ISI=new Vector(3,0) } + L_ISI=v_ISI.size + //print "L_ISI=",L_ISI + L_PTH=v_voltage.size + //print "L_PTH=",L_PTH +} + +rerun() \ No newline at end of file diff --git a/neuron/msomodel_default.hoc b/neuron/msomodel_default.hoc new file mode 100644 index 0000000..c0813fc --- /dev/null +++ b/neuron/msomodel_default.hoc @@ -0,0 +1,249 @@ +//---------- +// VARIABLES +//---------- + +seed_x = 1 + +iter = 1 + +n_runs_total = 1 +tstop = 2000 +dt = 0.025 + +v_init = -65 +celsius = 38 + +axon_soma_distance = 45 +dend_n = 2 +dend_n_syn = 10 +dend_n_seg = 20 +dend_length = 200 +dend_syn_spread = 0.5 +dend_syn_offset = 0.45 + +dend_0_exc_G = 11 +dend_1_exc_G = 11 + +dend_0_inh_G = 0 +dend_1_inh_G = 0 + +dend_0_exc_gain = 0.0008 +dend_1_exc_gain = 0.0008 + +dend_0_inh_gain = 0.001 +dend_1_inh_gain = 0.001 + +// Spike rate / threshold +spk_thres_exc_contra = 0 +spk_thres_exc_ipsi = 0 +spk_thres_inh_contra = 1 +spk_thres_inh_ipsi = 1 + +// Overrule variables using templating +//%%%TEMPLATE%%% + + +//---------------- +// CELL DEFINITION +//---------------- + +begintemplate Cell + public soma, axon, dend + create soma, axon, dend[1] + + proc init() { + ndend = $1 + dend_length = $2 + axon_soma_distance = $3 + create soma, axon, dend[ndend] + + soma { + nseg = 1 + diam = 20 // um + L = 40 // um + Ra = 200 // ohm.cm + cm = 1 // uF/cm2 + + insert kHT_VCN2003 + insert kLT_VCN2003 + insert na_VCN2003 + insert ih_VCN2003 + + ek = -70 // mV + ena = 55 // mV + eh = -43 // mV + + gkbar_kHT_VCN2003 = 0 // S/cm2 + gkbar_kLT_VCN2003 = 0 // S/cm2 + gnabar_na_VCN2003 = 0.1 // S/cm2 + ghbar_ih_VCN2003 = 0 // S/cm2 + gl_na_VCN2003 = 0.002 // S/cm2 + } + + axon { + nseg = 51 + diam = 2 // um + L = 400 // um^2 + Ra = 200 // ohm.cm + cm = 1 // uF/cm2 + + insert kHT_VCN2003 + insert kLT_VCN2003 + insert na_VCN2003 + insert ih_VCN2003 + + ek = -70 // mV + ena = 55 // mV + eh = -43 // mV + + // total G[nS]=gkbar*area*10 + + gkbar_kHT_VCN2003 = 0.02 // S/cm2 + gkbar_kLT_VCN2003 = 0.03 // S/cm2 + gnabar_na_VCN2003 = 0.3 // S/cm2 + ghbar_ih_VCN2003 = 0.0015 // S/cm2 + gl_na_VCN2003 = 0.002 // S/cm2 + } + + for i = 0, ndend-1 dend[i] { + nseg = 20 + diam = 3 // um + L = 200 // um + Ra = 200 // ohm.cm + cm = 1 // uF/cm2 + + insert pas + e_pas = -65 + g_pas = 0.002 + } + + connect dend[0](0), soma(0) + connect dend[1](0), soma(1) + connect axon(0), dend[1](axon_soma_distance/dend_length) + } +endtemplate Cell + + +//-------------------- +// CELL INITIALISATION +//-------------------- + +objectvar maincell +maincell = new Cell(dend_n, dend_length, axon_soma_distance) +access maincell.soma + + +//-------------------- +// EXCITATORY SYNAPSES +//-------------------- + +// Variables +objectvar syn0[10], gen0[10], syn1[10], gen1[10] + +// Event generators +for i = 0,9 { + // Contra + gen0[i] = new RelayStim(0.5) + gen0[i].seed(seed_x+i) + + // Ipsi + gen1[i] = new RelayStim(0.5) + gen1[i].seed(seed_x+10+i) +} + +// Contra synapses +maincell.dend[0] { + for i = 0,9 { + syn0[i] = new Alpha(i/((dend_n_syn-1)/dend_syn_spread)+dend_syn_offset) + + syn0[i].tau1 = 0.1 // ms + syn0[i].tau2 = 0.1 // ms + syn0[i].con = dend_0_exc_G // nS-mho + syn0[i].e = 0 // mV + } +} + +// Ipsi synapses +maincell.dend[1] { + for i = 0,9 { + syn1[i] = new Alpha(i/((dend_n_syn-1)/dend_syn_spread)+dend_syn_offset) + + syn1[i].tau1 = 0.1 // ms + syn1[i].tau2 = 0.1 // ms + syn1[i].con = dend_1_exc_G // nS-mho + syn1[i].e = 0 // mV + } +} + +// Connect generators and synapses +objref netcon0[10], netcon1[10] +e_delay = 0 // Synaptic delay (ms) + +for i = 0,9 { + netcon0[i] = new NetCon(gen0[i], syn0[i], spk_thres_exc_contra, e_delay, dend_0_exc_gain) + netcon1[i] = new NetCon(gen1[i], syn1[i], spk_thres_exc_ipsi, e_delay, dend_1_exc_gain) +} + + +//-------------------- +// INHIBITORY SYNAPSES +//-------------------- + +// Variables +objectvar syn_inhi0[10], gen_inhi0[10], syn_inhi1[10], gen_inhi1[10] + +// Event generators +for i = 0,9 { + gen_inhi0[i]=new GaussStim(0.5) // contralateral inputs + gen_inhi0[i].interval=2 + gen_inhi0[i].start=0 + gen_inhi0[i].number=10000 + gen_inhi0[i].factor=10 + gen_inhi0[i].seed(seed_x+30+i) + + gen_inhi1[i]=new GaussStim(0.5) // ipsilateral inputs + gen_inhi1[i].interval=2 + gen_inhi1[i].start=0 + gen_inhi1[i].number=10000 + gen_inhi1[i].factor=10 + gen_inhi1[i].seed(seed_x+50+i) +} + +// Somatic synapses +maincell.soma { + for i = 0,9 { + syn_inhi0[i] = new Alpha(0.5) + syn_inhi0[i].tau1 = 0.1 //ms + syn_inhi0[i].tau2 = 2 //ms + syn_inhi0[i].con = dend_0_inh_G //umho + syn_inhi0[i].e = -70 // mV + } + + for i = 0,9 { + syn_inhi1[i] = new Alpha(0.5) + syn_inhi1[i].tau1 = 0.1 //ms + syn_inhi1[i].tau2 = 2 //ms + syn_inhi1[i].con = dend_1_inh_G //umho + syn_inhi1[i].e = -70 // mV + } +} + +// Connect generators and synapses +objref netcon_inhi0[10],netcon_inhi1[10] +I_delay0 = 0 // Synaptic delay +I_delay1 = 0 // Synaptic delay + +for i = 0,9{ + netcon_inhi0[i] = new NetCon(gen_inhi0[i], syn_inhi0[i], spk_thres_inh_contra, I_delay0, dend_0_inh_gain) + netcon_inhi1[i] = new NetCon(gen_inhi1[i], syn_inhi1[i], spk_thres_inh_ipsi, I_delay1, dend_1_inh_gain) +} + + +//--------------------- +// LOAD SIMULATION CODE +//--------------------- + +load_file("msomodel_base.hoc") + +// Ending template (useful to quit automated runs) +//%%%TEMPLATE_ENDING%%% diff --git a/neuron/msomodel_passive_axon.hoc b/neuron/msomodel_passive_axon.hoc new file mode 100644 index 0000000..589351d --- /dev/null +++ b/neuron/msomodel_passive_axon.hoc @@ -0,0 +1,236 @@ +//---------- +// VARIABLES +//---------- + +seed_x = 1 + +iter = 1 + +n_runs_total = 1 +tstop = 2000 +dt = 0.025 + +v_init = -65 +celsius = 38 + +axon_soma_distance = 45 +dend_n = 2 +dend_n_syn = 10 +dend_n_seg = 20 +dend_length = 200 +dend_syn_spread = 0.5 +dend_syn_offset = 0.45 + +dend_0_exc_G = 11 +dend_1_exc_G = 11 + +dend_0_inh_G = 0 +dend_1_inh_G = 0 + +dend_0_exc_gain = 0.0008 +dend_1_exc_gain = 0.0008 + +dend_0_inh_gain = 0.001 +dend_1_inh_gain = 0.001 + +// Spike rate / threshold +spk_thres_exc_contra = 0 +spk_thres_exc_ipsi = 0 +spk_thres_inh_contra = 1 +spk_thres_inh_ipsi = 1 + +// Overrule variables using templating +//%%%TEMPLATE%%% + + +//---------------- +// CELL DEFINITION +//---------------- + +begintemplate Cell + public soma, axon, dend + create soma, axon, dend[1] + + proc init() { + ndend = $1 + dend_length = $2 + axon_soma_distance = $3 + create soma, axon, dend[ndend] + + soma { + nseg = 1 + diam = 20 // um + L = 40 // um + Ra = 200 // ohm.cm + cm = 1 // uF/cm2 + + insert kHT_VCN2003 + insert kLT_VCN2003 + insert na_VCN2003 + insert ih_VCN2003 + + ek = -70 // mV + ena = 55 // mV + eh = -43 // mV + + gkbar_kHT_VCN2003 = 0 // S/cm2 + gkbar_kLT_VCN2003 = 0 // S/cm2 + gnabar_na_VCN2003 = 0.1 // S/cm2 + ghbar_ih_VCN2003 = 0 // S/cm2 + gl_na_VCN2003 = 0.002 // S/cm2 + } + + axon { + nseg = 51 + diam = 2 // um + L = 400 // um^2 + Ra = 200 // ohm.cm + cm = 1 // uF/cm2 + + insert pas + e_pas = -65 + g_pas = 0.002 + } + + for i = 0, ndend-1 dend[i] { + nseg = 20 + diam = 3 // um + L = 200 // um + Ra = 200 // ohm.cm + cm = 1 // uF/cm2 + + insert pas + e_pas = -65 + g_pas = 0.002 + } + + connect dend[0](0), soma(0) + connect dend[1](0), soma(1) + connect axon(0), dend[1](axon_soma_distance/dend_length) + } +endtemplate Cell + + +//-------------------- +// CELL INITIALISATION +//-------------------- + +objectvar maincell +maincell = new Cell(dend_n, dend_length, axon_soma_distance) +access maincell.soma + + +//-------------------- +// EXCITATORY SYNAPSES +//-------------------- + +// Variables +objectvar syn0[10], gen0[10], syn1[10], gen1[10] + +// Event generators +for i = 0,9 { + // Contra + gen0[i] = new RelayStim(0.5) + gen0[i].seed(seed_x+i) + + // Ipsi + gen1[i] = new RelayStim(0.5) + gen1[i].seed(seed_x+10+i) +} + +// Contra synapses +maincell.dend[0] { + for i = 0,9 { + syn0[i] = new Alpha(i/((dend_n_syn-1)/dend_syn_spread)+dend_syn_offset) + + syn0[i].tau1 = 0.1 // ms + syn0[i].tau2 = 0.1 // ms + syn0[i].con = dend_0_exc_G // nS-mho + syn0[i].e = 0 // mV + } +} + +// Ipsi synapses +maincell.dend[1] { + for i = 0,9 { + syn1[i] = new Alpha(i/((dend_n_syn-1)/dend_syn_spread)+dend_syn_offset) + + syn1[i].tau1 = 0.1 // ms + syn1[i].tau2 = 0.1 // ms + syn1[i].con = dend_1_exc_G // nS-mho + syn1[i].e = 0 // mV + } +} + +// Connect generators and synapses +objref netcon0[10], netcon1[10] +e_delay = 0 // Synaptic delay (ms) + +for i = 0,9 { + netcon0[i] = new NetCon(gen0[i], syn0[i], spk_thres_exc_contra, e_delay, dend_0_exc_gain) + netcon1[i] = new NetCon(gen1[i], syn1[i], spk_thres_exc_ipsi, e_delay, dend_1_exc_gain) +} + + +//-------------------- +// INHIBITORY SYNAPSES +//-------------------- + +// Variables +objectvar syn_inhi0[10], gen_inhi0[10], syn_inhi1[10], gen_inhi1[10] + +// Event generators +for i = 0,9 { + gen_inhi0[i]=new GaussStim(0.5) // contralateral inputs + gen_inhi0[i].interval=2 + gen_inhi0[i].start=0 + gen_inhi0[i].number=10000 + gen_inhi0[i].factor=10 + gen_inhi0[i].seed(seed_x+30+i) + + gen_inhi1[i]=new GaussStim(0.5) // ipsilateral inputs + gen_inhi1[i].interval=2 + gen_inhi1[i].start=0 + gen_inhi1[i].number=10000 + gen_inhi1[i].factor=10 + gen_inhi1[i].seed(seed_x+50+i) +} + +// Somatic synapses +maincell.soma { + for i = 0,9 { + syn_inhi0[i] = new Alpha(0.5) + syn_inhi0[i].tau1 = 0.1 //ms + syn_inhi0[i].tau2 = 2 //ms + syn_inhi0[i].con = dend_0_inh_G //umho + syn_inhi0[i].e = -70 // mV + } + + for i = 0,9 { + syn_inhi1[i] = new Alpha(0.5) + syn_inhi1[i].tau1 = 0.1 //ms + syn_inhi1[i].tau2 = 2 //ms + syn_inhi1[i].con = dend_1_inh_G //umho + syn_inhi1[i].e = -70 // mV + } +} + +// Connect generators and synapses +objref netcon_inhi0[10],netcon_inhi1[10] +I_delay0 = 0 // Synaptic delay +I_delay1 = 0 // Synaptic delay + +for i = 0,9{ + netcon_inhi0[i] = new NetCon(gen_inhi0[i], syn_inhi0[i], spk_thres_inh_contra, I_delay0, dend_0_inh_gain) + netcon_inhi1[i] = new NetCon(gen_inhi1[i], syn_inhi1[i], spk_thres_inh_ipsi, I_delay1, dend_1_inh_gain) +} + + +//--------------------- +// LOAD SIMULATION CODE +//--------------------- + +load_file("msomodel_base.hoc") + +// Ending template (useful to quit automated runs) +//%%%TEMPLATE_ENDING%%% diff --git a/neuron/na_VCN2003.mod b/neuron/na_VCN2003.mod new file mode 100644 index 0000000..d0d80e4 --- /dev/null +++ b/neuron/na_VCN2003.mod @@ -0,0 +1,107 @@ +TITLE Na channels in VCN auditory neurons of guinea pig + +: na=gna*m^3*h +: based on Rothman and Manis 2003 +: Modifications by Yi Zhou for an MSO model + + +INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)} + +NEURON { + SUFFIX na_VCN2003 + USEION na READ ena WRITE ina + NONSPECIFIC_CURRENT il + RANGE gnabar + RANGE m_inf,h_inf + RANGE tau_m,tau_h + RANGE m_exp,h_exp + RANGE ina,gna + RANGE gl,el + +} + + +UNITS { + (mA) = (milliamp) + (mV) = (millivolt) +} + +PARAMETER { + gnabar = 0.2 (mho/cm2) + ena=55 (mV) + gl= 4.0e-4 (mho/cm2) + el=-57 (mV) + celsius=22 (degC) + dt (ms) + v (mV) + +} + +STATE { + m h +} + +ASSIGNED { + gna (mho/cm2) + ina (mA/cm2) + il (mA/cm2) + m_inf + h_inf + tau_m + tau_h + m_exp + h_exp + tadj3 + +} + + +BREAKPOINT { + SOLVE states + gna=gnabar * m*m*h + ina = gna * (v - ena) + il=gl*(v-el) +} + + + +PROCEDURE states() { : this discretized form is more stable + evaluate_fct(v) + m = m + m_exp * (m_inf - m) + h = h + h_exp * (h_inf - h) + VERBATIM + return 0; + ENDVERBATIM +} + +UNITSOFF +INITIAL { +: +: Q10 was assumed to be 3 for both currents +: + tadj3 = 3.0 ^ ((celsius-22)/ 10 ) + + evaluate_fct(v) + m= m_inf + h= h_inf + } + +PROCEDURE evaluate_fct(v(mV)) { + + m_inf = 1 / (1+exp(-(v + 38) / 7)) + h_inf = 1 / (1+exp((v + 65) / 6)) + + tau_m = (10 / (5*exp((v+60) / 18) + 36*exp(-(v+60) / 25))) + 0.04 + tau_h = (100 / (7*exp((v+60) / 11) + 10*exp(-(v+60) / 25))) + 0.6 + + tau_m=tau_m/tadj3 + tau_h=tau_h/tadj3 + + m_exp = 1 - exp(-dt/tau_m) + h_exp = 1 - exp(-dt/tau_h) + +} + +UNITSON + + diff --git a/neuron/readme.txt b/neuron/readme.txt new file mode 100644 index 0000000..7dd6dcb --- /dev/null +++ b/neuron/readme.txt @@ -0,0 +1,30 @@ +This model package provides NEURON codes for a multi-compartment model described +in Zhou Y, Carney LH, Colburn HS (2005). The model simulates recent physiological +results from the gerbil medial superior olive (MSO) that reveal that blocking +glycinergic inhibition can shift the tuning for the interaural time difference (ITD) +of the cell (Brand et al., 2002). + +The model has a bipolar structure and incorporates two anatomic observations +in the MSO: (1) the axon arises from the dendrite that receives ipsilateral +inputs (Smith, 1995); and (2) inhibitory synapses are located primarily +on the soma in adult animals. Specific Hodgkin-Huxley type of ionic channels +are inserted into the axon to generate action potentials: fast sodium (Ina), +low threshold potassium (ILTK), high threshold potassium (IHTK), +hyperpolarization-active inward current (Ih), and the leakage channel (Il). + +The kinetics of all channels are based on the model of Type II cells in the +anteroventral cochlear nuclei (Rothman and Manis, 2003c). The soma has only +fast sodium channels. Ionic channels and the driving function for synaptic +inputs are described in .mod files. The hoc files describe simulation procedures +for modeling neural responses to ITDs (MSO_Zhouetal_2005.hoc) and a test kit +for exploring input characteristics (test_inputs.hoc). Use NEURON command nrnivmodl +to compile *.mod files, and then execute a *.hoc file. Please refer to NEURON +documentations for further instructions on compiling mod files and +launching hoc files on different operation systems. + +The file MSO_Zhouetal_2005.hoc produces the rate-ITD responses +in Figures. 5, 6, 8, 11 in Zhou Y, Carney LH, Colburn HS (2005). + +Go to http://www.bu.edu/dbin/binaural/MSOdemo for more instructions on using the model. + +Any questions for implementing the model should be directed to Yi Zhou (yizhou@bu.edu). diff --git a/neuron/relaystim.mod b/neuron/relaystim.mod new file mode 100644 index 0000000..f0d9ecf --- /dev/null +++ b/neuron/relaystim.mod @@ -0,0 +1,75 @@ +COMMENT +This mod file is based on netstim.mod in NEURON +modified by Yarmo Mackenbach to relay spike times provided by the hoc script +ENDCOMMENT + +NEURON { + POINT_PROCESS RelayStim + RANGE x, timeToNextEvt, refrac +} + +PARAMETER { + timeToNextEvt = 10 (ms) <1e-9,1e9> : time until next spiketime (ms) + refrac = 0.1 (ms) <0,3> : absolute refractory period (ms) +} + +ASSIGNED { + x : unknown significance + on : is spike generation active? +} + +PROCEDURE seed(x) { + set_seed(x) +} + +INITIAL { + : init variables + x = 0 + on = 0 + + : wait until first event + net_send(timeToNextEvt, 3) +} + +NET_RECEIVE(w) { + : MAIN CALL + if (flag == 1 && on == 1) { + : because somehow this is a crucial step + x = 1 + + : output event to synapse + net_event(t) + + : disable spike generation if timeToNextEvt is null + if (timeToNextEvt == 0) { + on = 0 + } + + : if spike generation is turned on + if (on == 1) { + : wait until next event + net_send(timeToNextEvt, 1) + } + + : refractory period call (somehow required) + net_send(refrac, 2) + } + + : REFRACTORY PERIOD CALL (END OF REFRACTORY PERIOD) + if (flag == 2) { + : unknown signifance (but a required step) + x = 0 + } + + : INITIAL CALL + if (flag == 3) { + : if spike generation is not yet turned on + if (on == 0) { + : turn on spike generation + on = 1 + + : call net_receive directly + net_send(0, 1) + } + } +} diff --git a/neuron/test_inputs.hoc b/neuron/test_inputs.hoc new file mode 100644 index 0000000..c37363d --- /dev/null +++ b/neuron/test_inputs.hoc @@ -0,0 +1,601 @@ +// ***************** Test Inputs ********************** + +// ------------ Introduction -------------------------------- +// 01/01/2005 +// This program generats periodic input spikes that simulate afferent inputs to cells in the Medial Superio Olive (MSO) +// of the mammalian brainstem as described in +// Zhou, Carney, and Colburn (2005) J. Neuroscience, 25(12):3046-3058 + + + + + +// ******* Codes start here **************************************** + +seed_x=xred("Please choose a random seed number <1,1000>", 1, 1, 1000) //string,default,min,max + +strdef inputDir,outputDir +outputDir="D:\msomodel\" +inputDir="D:\msomodel\" + + + load_file("nrngui.hoc") + + create cell + + access cell + + + + + + //**** Excitatory and inhibitory synapses **************** + //--------------------------------------------------------- + + //---- parameter description ---------------------- + //syn0[]: E alpha synapses on contralateral dendrite, dend0 + //syn1[]: E alpha synapses on ipsilateral dendrite, dend1 + //gen0[]: E stimuli to syn0 + //gen1[]: E stimuli to syn1 + //netcon0[]: input- E synapse connection on dend0 + //netcon1[]: input- E synapse connection on dend1 + + //syn_inhi0[]: I alpha synapse on the soma + //gen_inhi0[]: I stimuli to syn_inhi0 + //netcon_inhi0[]: input- I synapse connection on the soma + + + // ************ Excitatory synapses on two dendrites ******** + objectvar syn0[10], gen0[10],syn1[10], gen1[10] + + //----- add spike, virtual spikes + for i=0, 9{ + gen0[i] = new GaussStim(0.5) + gen0[i].interval=2 // [ms] input period + gen0[i].start=0 // [ms] arrival time of input spikes + gen0[i].number=10000 + gen0[i].factor=20 // phase-locking factor F + + gen1[i] = new GaussStim(0.5) + gen1[i].interval=2 + gen1[i].start=0 + gen1[i].number=10000 + gen1[i].factor=20 + + gen0[i].seed(seed_x+i) + gen1[i].seed(seed_x+10+i) + } + + //----- add EE expsyn at dend0 + cell { //contra side + for i=0,9 { + n=20 //n=nseg doesn't change after nseg + syn0[i] = new Alpha(i*1/n+1/n/2) // scatter along the firsthalf dend + + syn0[i].tau1 =0.1 //ms + syn0[i].tau2 =0.1 //ms + syn0[i].con=10 //nS- mho + syn0[i].e =0 // mV + } + } + + cell { //ipsi side + for i=0,9 { + n=20 //n=nseg + syn1[i] = new Alpha(i*1/n+1/n/2) + + syn1[i].tau1 =0.1 //ms may add some jitters + syn1[i].tau2 =0.1 //ms + syn1[i].con=0 //nS-mho + syn1[i].e =0 // mV + } + } + + + //--- Connect gen[] with syn[] + objref netcon0[10],netcon1[10] + + Rave_E_contra=240 // spikes/sec + Rave_E_ipsi=240 + + x_thres_contra=1-Rave_E_contra*gen0[0].interval/1000 // prob(x>1-p)=p + x_thres_ipsi=1-Rave_E_ipsi*gen1[0].interval/1000 + + if (x_thres_contra<=0) {x_thres_contra=1e-4} + if (x_thres_ipsi<=0) {x_thres_ipsi=1e-4} + + + e_delay=0 // no synaptic delay + + for i=0,9 { + netcon0[i]=new NetCon(gen0[i], syn0[i],x_thres_contra,e_delay,0.001) + //thres delay gain[uS] + netcon1[i]=new NetCon(gen1[i], syn1[i],x_thres_ipsi,e_delay,0.001) + //thres delay gain[uS] + } + + + + + v_init=-65 + + + + //---- parameter and data vector descriptions ------------ + + // v_ISI: [ms] AP event time + // v_voltage: occurrence of events at each dt + + + tstop=1000 // [ms] stimulus duration + dt=0.025 // [ms] simulation time step + bin=0.1 // [ms] binwidth of APs + + cvode.active(0) // constant integration time step + + + + + objref v_ISI,v_voltage + objref pre_E + + v_ISI=new Vector() + v_voltage=new Vector(tstop/dt+2,0) + + + + //==== for input ISIs + netcon0[0].record(v_ISI)//record input event time + v_voltage.record(&gen0[0].x,dt) //record the event + AP_threshold=x_thres_contra // reset if record inputs + + + + //**** function "pth": period time histogram **** + //----------------------------------------------- + //---- data vector description ------------ + // v_save: bin v_voltage for each run + // v_pth: wrap v_all with length=gen.interval/bin and normalize + // vs_real: vector strength + // phase_real: mean phase + + // get_pth(): get v_save from each run + + + objref v_save, v_all, v_pth, g_pth, x_pth, v_stand0 + + proc get_pth() { + 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 + } + + proc pth() { + x=0 + y=0 + vs_real=0 + angel_real=0 + length=gen0[0].interval/bin //length of period points + + N_fold=v_all.size/length + v_pth= new Vector(length,0) + x_pth= new Vector(length,0) // phase vector + x_pth.indgen(1/length) + + //==== wrapping over one period + for i=0,N_fold-1 { + for j=0,length-1 { + v_pth.x[j]=v_pth.x[j]+v_all.x[length*i+j] + } + } + v_pth.div(v_pth.sum) + + v_stand0=new Vector(x_pth.size,0) // dend[0] input PSH + v_stand0.indgen(1/length) + v_stand0.apply("Gauss0") + v_stand0.rotate(gen0[0].start/bin) + + //==== calculate the VS and mean phase of responses + for i=0,length-1 { + x=x+v_pth.x[i]*cos(2*PI*i/length) + y=y+v_pth.x[i]*sin(2*PI*i/length) + } + vs_real=sqrt(x^2+y^2) //compute VS + print "VS=",vs_real + + angel_real=atan(y/x) + if(x<0 ) { angel_real=angel_real+PI + } else if (y<0 && x>0) { + angel_real=angel_real+2*PI + } + print "Phase=",angel_real/(2*PI) + + + g_pth.size(0,1,0,0.3) + g_pth.align(0.5,0.2) + //g_pth.label(0.9,0.1,"Phase") + g_pth.color(1) + g_pth.label(0.3,0.9,"Period Histogram") + g_pth.begin() + v_pth.mark(g_pth,x_pth,"+",9,1) + v_pth.line(g_pth,x_pth,1,2) + g_pth.color(2) + v_stand0.line(g_pth,x_pth,2,2) + + +} + + proc VS() { + vs=exp(-PI^2/(2*gen0[0].factor^2)) + } + + func Gauss0() { //for PST + sigma2=0.5/gen0[0].factor + return exp(-($1-0.5)^2/(2*sigma2^2))/((sqrt(2*PI)*sigma2)*length) +} + + + //-------------- END pth () ------------------------- + //--------------------------------------------------- + + + + + + //**** function "ISI": interspike interval **** + //----------------------------------------------- + //---- data vector description ------------ + + // v_ISI: event times + // 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 statistics for each run + + // get_ISI(): get data for each run + + objref hist_ISI,xtime,g_ISI,hist, v_stand + objref mean_T, sq_T, num_T + objref temp,t1, v1,v2 // operating vectors + + mean_ISI=0 + OVF=0 // overflow + EI=0 + + mean_ISI2=0 + std_ISI=0 + sq_sum=0 + CV_ISI=0 + + proc get_ISI() { + + 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 initial interval + print "mean ISI interval ", v2.mean + print "ISI std ",v2.stdev + + temp = new Vector(v2.size,0) + temp.add(v2) + temp.sub(v2.mean) // for calculate the square sum + + 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 "**********************************" + + hist=v2.histogram(0, topbin+OVFbin, bin_ISI) + + if(hist.size-hist_ISI.size==1) { hist.resize(hist_ISI.size)} + hist_ISI.add(hist) + + } + + proc ISI() { + hist_ISI.div(hist_ISI.sum) //normalize + OVF=hist_ISI.sum(topbin/bin_ISI, (topbin+OVFbin)/bin_ISI-1) + EI=hist_ISI.sum(0,1.5*T/bin_ISI-1) // the first whole interval + + hist_ISI.resize(topbin/bin_ISI+1) + + xtime=new Vector(hist_ISI.size,0) + xtime.indgen(bin_ISI) + + v_stand=new Vector(hist_ISI.size,0) + v_stand.indgen(bin_ISI) + v_stand.apply("Gauss") + + g_ISI.size(0,topbin,0,0.2) + g_ISI.align(0.5,0.5) + //g_ISI.label(0.9,0.2,"T") + g_ISI.label(0.2,0.9,"ISI") + g_ISI.beginline() + hist_ISI.mark(g_ISI,xtime,"o",6,1) + hist_ISI.line(g_ISI,xtime,1,2) + v_stand.line(g_ISI,xtime,2,2) + + + //==== 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_ISI2=t2/num_T.sum() //overall ISI mean + + t1=new Vector(N_run,0) + t1.add(mean_T) + t1.sub(mean_ISI2) + 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_ISI2>0) { + CV_ISI=std_ISI/mean_ISI2 //coefficient of variation + } else{ CV_ISI=1 + } + + } + + func Gauss() { //for ISI + sigma=sqrt(2)/2*(gen0[0].interval/gen0[0].factor) //std=sqrt(2)*T/2/factor + return exp(-($1-gen0[0].interval)^2/(2*sigma^2))*bin/(sqrt(2*PI)*sigma) +} + + + //----------- END ISI() ------------------------- + //----------------------------------------------- + + + N_run=10 // ten trials + + + //----- reset function for all simulations + //------------------------------------------ + + proc reset() { + + + //==== Reset synapse, and inputs + + getcon() + getstim() + changelevel() + + + //==== Reset all vec and var + N_run=10 + + v_voltage=new Vector(tstop/dt+2,0) + v_voltage.record(&gen0[0].x,dt) //record the event + + mean_T=new Vector(N_run,0) + sq_T=new Vector(N_run,0) + num_T=new Vector(N_run,0) + + T=gen0[0].interval + bin_ISI=bin + topbin=int(10*T/bin_ISI)*bin_ISI // 10 intervals and 1 for OVF + OVFbin=int(1*T/bin_ISI)*bin_ISI + + mean_ISI2=0 + std_ISI=0 + sq_sum=0 + CV_ISI=0 + + mean_ISI=0 + OVF=0 + EI=0 + + hist_ISI=new Vector((topbin+OVFbin)/bin_ISI+1,0) + v_all=new Vector(tstop/dt+2,0) + v_all.rebin(v_all,bin/dt) + + + } + + //---------- end reset() ---------------- + //--------------------------------------- + + + + + + + //----- Main simulation function -------------------- + // Measure model response at one ITD condition + //--------------------------------------------------- + + + proc rerun() { + + reset() + + //********************************* + while (N_run>0) { + v_ISI=new Vector() + v_voltage=new Vector(tstop/dt+2,0) + + //==== for input ISI and PS + netcon0[0].record(v_ISI)//record input event time + v_voltage.record(&gen0[0].x,dt) //record the event + + finitialize(v_init) + fcurrent() + + integrate() + get_pth() + get_ISI() + + N_run=N_run-1 + } + + + N_run=10 + + g_ISI.erase_all() + g_pth.erase_all() + ISI() + pth() + + print "----------- END -----------------------" +} + + + proc integrate() { + + while (t < tstop) { + fadvance() // advance solution + } + print "N_run=",N_run // show run time + print "-------------------------" + print "Number_input_E=",v_ISI.size + + + if(v_ISI.size<=2) { v_ISI=new Vector(3,0) } + L_ISI=v_ISI.size + //print "L_ISI=",L_ISI + L_PTH=v_voltage.size + //print "L_PTH=",L_PTH +} + + + + +//******* update synapse and input information **** + + + proc changelevel() { + + x_thres_contra=1-Rave_E_contra*gen0[0].interval/1000 // prob(x>1-p)=p + x_thres_ipsi=1-Rave_E_ipsi*gen1[0].interval/1000 + + + + if (x_thres_contra<=0) {x_thres_contra=1e-4} + if (x_thres_ipsi<=0) {x_thres_ipsi=1e-4} + + + for i=0,9 { + netcon0[i].threshold=x_thres_contra + netcon1[i].threshold=x_thres_ipsi + } + + } + + proc getcon() { + for i=0,9 { + syn0[i].con=syn0[0].con //umho + syn1[i].con=syn1[0].con //umho + + } + print "*******************" + print "syn0 con =",syn0[1].con + print "syn1 con =",syn1[1].con + + print "--------------------------" + } + + proc getstim() { + //------- EE ------------- + for i=0, 9{ + gen0[i].interval=gen0[0].interval // + gen0[i].start=gen0[0].start // identify ITD + gen0[i].factor=gen0[0].factor + + gen1[i].interval=gen1[0].interval + gen1[i].start=gen1[0].start // identify ITD + gen1[i].factor=gen1[0].factor + } + + + + print "E0 start =",gen0[1].start + print "E0 interval =",gen0[1].interval + print "E1 start =",gen1[1].start + print "E1 interval =",gen1[1].interval + + print "*******************" + + } + //----------- END rerun() ----------------------------------- + //----------------------------------------------------------- + + + + //**** GUI "synapse and stimuli" **** + //--------------------------------------------- + + objref syn_con, gen_start + + syn_con = new HBox() + syn_con.intercept(1) + + xpanel("") + xlabel("E0 syn contra") + xvalue("G [nS]","syn0[0].con") + xpanel() + + + syn_con.intercept(0) + syn_con.map("SYN",0,150,-1,50) + + + //-------------------------------------------------------------------------- + + gen_start = new HBox() + gen_start.intercept(1) + + xpanel("") + xlabel("E0 input contra") + xvalue("delay [ms]","gen0[0].start") + xvalue("period","gen0[0].interval") + xvalue("synchrony factor","gen0[0].factor") + xvalue("rate [spikes/sec]","Rave_E_contra") + xpanel() + + + + gen_start.intercept(0) + gen_start.map("GEN",0,350,-1,50) + + +//-------------------------------------------------------------------------- + + + objref boxv + boxv=new VBox() + boxv.intercept(1) + xpanel("") + xbutton("calculate VS of E0 input","VS()") + xvalue("vs") + xbutton("Reset","reset()") + xbutton("Single Run", "rerun()") + xbutton("Quit", "quit()") + + g_ISI=new Graph() + g_pth=new Graph() + xpanel() + boxv.intercept(0) + boxv.map("CONTROL",600,50,-1,50) + //---------------------------------------------------------------- + + + +