378 lines
11 KiB
Plaintext
378 lines
11 KiB
Plaintext
// 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() |