/*
Stimulus protocol as used by Sjöström et al. (2001)
Demonstrating the synaptic plasticity rule by Ebner et al. (2019)
https://doi.org/10.1016/j.celrep.2019.11.068
Made to work with the Hay et al. (2011) L5b pyramidal cell model
*/
load_file("nrngui.hoc")
// ## Object references
objref freqs, loc_prox, loc_dist
objref trec, vrec_prox, vrec_dist
objref pre_post, post_pre
objref syn, stim, nc, ic
// ## Parameters
REPS = 5 // Number of pre- and postsynaptic spikes
DT = 0.025 // ms, integration step.
AMP = 2.7 // nA, amplitude of current injection to trigger postsynaptic spikes
DUR = 5.0 // ms, duration of the current injection
WARM_UP = 2300 // ms, silence phase before stimulation
COOL_DOWN = 100 // ms, silence phase after stimulation
WEIGHT = 0.0035 // µS, conductance of (single) synaptic potentials
LOCATION = 1 // 1 ... proximal (90 µm) / 2 ... distal (669 µm)
freqs = new Vector() // stimulation frequencies (Hz)
// freqs.append(0.1, 10, 20, 40, 50) // use this to include 0.1 Hz condition. This will take some time.
freqs.append(10, 20, 40, 50)
// ## Create cell
load_file("import3d.hoc")
objref L5PC
strdef morphology_file
morphology_file = "morphologies/cell1.asc"
load_file("models/L5PCbiophys4.hoc")
load_file("models/L5PCtemplate.hoc")
L5PC = new L5PCtemplate(morphology_file)
// ## Create synapse
L5PC.apic[1] loc_prox = new SectionRef()
L5PC.apic[50] loc_dist = new SectionRef()
if(LOCATION == 1) {loc_prox.sec syn = new Syn4P(0.5)
} else {loc_dist.sec syn = new Syn4P(0.5)}
syn.tau_a = 0.2 // time constant of EPSP rise
syn.tau_b = 2 // time constant of EPSP decay
syn.e = 0 // reversal potential
syn.w_pre_init = 0.5 // pre factor initial value
syn.w_post_init = 2.0 // post factor initial value
syn.s_ampa = 0.5 // contribution of AMPAR currents
syn.s_nmda = 0.5 // contribution of NMDAR currents
syn.tau_G_a = 2 // time constant of presynaptic event G (rise)
syn.tau_G_b = 50 // time constant of presynaptic event G (decay)
syn.m_G = 10 // slope of the saturation function for G
syn.A_LTD_pre = 3e-3 // amplitude of pre-LTD
syn.A_LTP_pre = 33e-4 // amplitude of pre-LTP
syn.A_LTD_post = 36e-5 // amplitude of post-LTD
syn.A_LTP_post = 2e-1 // amplitude of post-LTP
syn.tau_u_T = 10 // time constant for filtering u to calculate T
syn.theta_u_T = -60 // voltage threshold applied to u to calculate T
syn.m_T = 1.7 // slope of the saturation function for T
syn.theta_u_N = -30 // voltage threshold applied to u to calculate N
syn.tau_Z_a = 1 // time constant of presynaptic event Z (rise)
syn.tau_Z_b = 15 // time constant of presynaptic event Z (decay)
syn.m_Z = 6 // slope of the saturation function for Z
syn.tau_N_alpha = 7.5 // time constant for calculating N-alpha
syn.tau_N_beta = 30 // time constant for calculating N-beta
syn.m_N_alpha = 2 // slope of the saturation function for N_alpha
syn.m_N_beta = 10 // slope of the saturation function for N_beta
syn.theta_N_X = 0.2 // threshold for N to calculate X
syn.theta_u_C = -68 // voltage threshold applied to u to calculate C
syn.theta_C_minus = 15 // threshold applied to C for post-LTD (P activation)
syn.theta_C_plus = 35 // threshold applied to C for post-LTP (K-alpha activation)
syn.tau_K_alpha = 15 // time constant for filtering K_alpha to calculate K_alpha_bar
syn.tau_K_gamma = 20 // time constant for filtering K_beta to calculate K_gamma
syn.m_K_alpha = 1.5 // slope of the saturation function for K_alpha
syn.m_K_beta = 1.7 // slope of the saturation function for K_beta
syn.s_K_beta = 100 // scaling factor for calculation of K_beta
// ## Custom init procedure
proc init(){ local i
finitialize(-70)
dt = 10
i = 0
for i = 1,224 fadvance() // Run warm-up phase using big steps
print t
dt = DT
}
// ## Current trace vector
obfunc current_trace() { local freq, delta_t, total_t, start_t, end_t, counter localobj trace
freq = $1
delta_t = $2
total_t = $3
trace = new Vector(total_t/DT, 0)
for counter = 0, REPS - 1 {
start_t = (0 + delta_t + counter * (1000.0 / freq) + WARM_UP)
end_t = (start_t + DUR) // this is the actual duration of the current injection
trace.fill(AMP, start_t/DT, end_t/DT)
}
return trace
}
// ## Stimulation protocol; goes through all frequencies at a particular timing
obfunc stimulation() { local delta_t, counter localobj delta_w, ic_trace
delta_t = $1
delta_w = new Vector()
// Create objects
stim = new NetStim()
stim.number = REPS
stim.start = WARM_UP
stim.noise = 0
nc = new NetCon(stim, syn, 0, 0, WEIGHT)
L5PC.soma ic = new IClamp(0.5)
ic.del = 0
ic.dur = 1e9 // not the actual duration
for counter = 0, freqs.size() - 1 {
print "\n*Loop ", counter + 1, "of ", freqs.size, "// f = ", freqs.x[counter], "Hz"
stim.interval = 1000.0 / freqs.x[counter] // calculate interval
total_time = WARM_UP + REPS * (1000.0 / freqs.x[counter]) + COOL_DOWN
tstop = total_time
ic_trace = new Vector()
ic_trace = current_trace(freqs.x[counter], delta_t, total_time)
ic_trace.play(&ic.amp, DT) // play the trace vector into the IClamp protocol
// Recording
trec = new Vector()
trec.record(&t)
vrec_prox = new Vector()
vrec_prox.record(&loc_prox.sec.v(0.5)) // voltage at proximal location
vrec_dist = new Vector()
vrec_dist.record(&loc_dist.sec.v(0.5)) // voltage at distal location
init()
run()
w_init = syn.w_pre_init * syn.w_post_init
delta_w_pre = syn.w_pre - syn.w_pre_init // Difference of present and initial values
delta_w_post = syn.w_post - syn.w_post_init
if(freqs.x[counter] == 0.1) {
w_final_pre = syn.w_pre_init + delta_w_pre * 10 // 10 repetitions at 0.1 Hz
w_final_post = syn.w_post_init + delta_w_post * 10
} else {
w_final_pre = syn.w_pre_init + delta_w_pre * 15 // 15 repetitions at 10, 20, 40 and 50 Hz
w_final_post = syn.w_post_init + delta_w_post * 15
}
if (w_final_pre > 1.0) { w_final_pre = 1.0 }
if (w_final_pre < 0.0) { w_final_pre = 0.0 }
if (w_final_post > 5.0) { w_final_post = 5.0 }
if (w_final_post < 0.0) { w_final_post = 0.0 }
w_final = w_final_pre * w_final_post
outcome = (w_final - w_init) / w_init * 100
print "dw(pre) = ", delta_w_pre, "--> ", (w_final_pre - syn.w_pre_init) / syn.w_pre_init * 100, "%"
print "dw(post) = ", delta_w_post, "--> ", (w_final_post - syn.w_post_init) / syn.w_post_init * 100, "%"
print "w(start) = ", w_init
print "w(end) = ", w_final
print "normalized weight change at f = ", freqs.x[counter], " --> ", outcome, "%"
delta_w.append(outcome)
}
return delta_w
}
// ## Main
pre_post = new Vector()
post_pre = new Vector()
print "\n*** Starting Pre-Post ***"
pre_post = stimulation(10)
print "\n*** Starting Post-Pre ***"
post_pre = stimulation(-10)
// ## Visualization
objref voltplot
voltplot = new Graph(0) // proximal & distal voltage over time (final run only)
voltplot.size(0, 10625, -85, 10)
voltplot.view(-0.3, -85, 10625, 105, 65, 105, 300, 200)
vrec_prox.line(voltplot, 1, 2)
vrec_dist.line(voltplot, 3, 2)
objref compplot, compvec1, compvec2
compvec1 = new Vector()
compvec1.append(14, 29, 53, 56) // Experimental data (Sjöström et al., 2001) pre-post
compvec2 = new Vector()
compvec2.append(-41, -35, 56, 75) // Experimental data (Sjöström et al., 2001) post-pre
compplot = new Graph(0)
compplot.size(0, 3, -50, 80)
compplot.view(-0.3, -50, 3.3, 130, 65, 105, 300, 200)
compvec1.line(compplot, 9, 2)
pre_post.line(compplot, 2, 2)
compvec2.line(compplot, 9, 2)
post_pre.line(compplot, 2, 2)
