yscience/pricklypear
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Added neuron model

Browse Source
This commit is contained in:
Yarmo Mackenbach 2020-04-01 14:41:11 +02:00
parent 6e1f465865
commit 216ccec72e
15 changed files with 3518 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
neuron/.gitignore vendored Normal file
View File

@ -0,0 +1,5 @@
*.o
*.c
*.lnk
*.tmp
*.dll

68
neuron/Alpha.mod Normal file
View File

@ -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
}

1244
neuron/MSO_Zhouetal_2005.hoc Normal file
View File

File diff suppressed because it is too large Load Diff

223
neuron/gaussstim.mod Normal file
View File

@ -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)
}
}
}

80
neuron/ih_VCN2003.mod Normal file
View File

@ -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

95
neuron/kHT_VCN2003.mod Normal file
View File

@ -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

93
neuron/kLT_VCN2003.mod Normal file
View File

@ -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

34
neuron/mosinit.hoc Normal file
View File

@ -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)
}

378
neuron/msomodel_base.hoc Normal file
View File

@ -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()

249
neuron/msomodel_default.hoc Normal file
View File

@ -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%%%

236
neuron/msomodel_passive_axon.hoc Normal file
View File

@ -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%%%

107
neuron/na_VCN2003.mod Normal file
View File

@ -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

30
neuron/readme.txt Normal file
View File

@ -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).

75
neuron/relaystim.mod Normal file
View File

@ -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)
}
}
}

601
neuron/test_inputs.hoc Normal file
View File

@ -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)
//----------------------------------------------------------------