/* Implementation of BCM-like synaptic plasticity model in a
* conductance-based model. The model contains Na, KDR, KA and h currents,
* and employs GHK-based implementations of AMPA and NMDA receptors. A
* calcium-based synaptic plasticity rule (from Shouval et al., PNAS, 2002)
* is the incorporated into the model to update the weight w. All weight
* updates are performed in the MOD file "Wghkampa.mod". This file helps
* recreate Fig. 2A and Fig. 2B of the following paper:
* Narayanan R, Johnston D. The h current is a candidate mechanism for
* regulating the sliding modification threshold in a BCM-like synaptic
* learning rule. J Neurophysiol. 2010 Aug;104(2):1020-33. doi:
* 10.1152/jn.01129.2009. PubMed PMID: 20554832. PubMed Central PMCID:
* PMC2934916.
* Implemented by Rishikesh Narayanan. Contact rishi.n@gmail.com.
*/
create dend
access dend
// All parameters are set as if this is an apical dendritic
// section sitting around 200um from the soma.
objref cvode
cvode = new CVode()
cvode.active(1)
objref ampa, nmda, s
objref ncl, weight
double isilist[1]
//**************************************************************//
// The Parameters
//**************************************************************//
tstop = 30000
steps_per_ms = 40
dt = 0.025
// --------------------------------------------------------------
// Resistances and Capacitances
// --------------------------------------------------------------
ra = 150
rm = 28000
rmdend = rm/2
c_m = 0.75
cmdend = c_m*2
v_init = -65
celsius = 34
// --------------------------------------------------------------
// Equilibrium potentials
// --------------------------------------------------------------
Eh=-30
Ek = -90
Ena = 55
//--------------------------------------------------------------
// Active conductance densities and other related parameters
//--------------------------------------------------------------
gna=0.03
nash=10
gkdr=0.005
gka=0.044 // Base value of gka=0.022*(1+xdist/100)
gh=0.0042 // Base value of gh=0.0006*(1+3*xdist/100)
//--------------------------------------------------------------
// NetStim parameters
//--------------------------------------------------------------
tsyn=20
intrvl=10
nmbr=900
noisevar=0.0
//--------------------------------------------------------------
// Synapse parameters
//--------------------------------------------------------------
P=1e-6 // AMPA Permeability in cm/s.
NAR=1.5 // Poirazi et al., 2003 use NAR of 0.6 for Trunk and 2.5 for
// non trunk - this is somewhere in between!
w=0.25
//--------------------------------------------------------------
// The code.
//--------------------------------------------------------------
dend {
L=50 // Take a 50um dendrite
diam=1 // It is 1um thick
}
// Insert passive elements.
insert pas
e_pas = v_init
Ra=ra
cm=cmdend
g_pas=1/rmdend
// Insert active elements.
dend {
insert nas gbar_nas=gna ar_nas=0.8
insert kdr gkdrbar_kdr=gkdr
insert kad gkabar_kad=0
insert hd
insert cad
ghdbar_hd = gh
sh_nas=nash-4 // xdist=200; nash-8*(xdist-100)/200
vhalfl_hd=-73-2 //xdist=200; -73-4*(xdist-100)/200
gkabar_kad = gka
ek = Ek
ena = Ena
ehd_hd=Eh
finitialize(v_init)
fcurrent()
for (x) {
e_pas(x)=v(x)
e_pas(x)=e_pas(x)+(ina(x)+ik(x)+i_hd(x))/g_pas(x)
}
}
s = new NetStim(0.5)
s.interval=intrvl // ms (mean) time between spikes
s.number=nmbr // (average) number of spikes
s.start=tsyn // ms (mean) start time of first spike
s.noise=0 // range 0 to 1. Fractional randomness.
ampa=new Wghkampa(0.5)
ncl=new List()
dend ncl.append(new NetCon(s, ampa, 0, 0, 0))
ampa.taur = 2
ampa.taud = 10
ampa.Pmax = P
ampa.winit = w
nmda=new ghknmda(0.5)
dend ncl.append(new NetCon(s, nmda, 0, 0, 0))
nmda.taur = 5
nmda.taud = 50
nmda.Pmax = P*NAR
/********************************************************************/
proc update_Gh() {
dend{
ghdbar_hd = gh
}
update_init()
}
/********************************************************************/
proc update_init(){
finitialize(v_init)
fcurrent()
for (x) {
e_pas(x)=v(x)
e_pas(x)=e_pas(x)+(ina(x)+ik(x)+i_hd(x))/g_pas(x)
}
}
/********************************************************************/
// This function generates the % change in synaptic weight as a
// function of induction frequency for various values of gh
// (Similar to Fig. 2A of the paper). The output is saved in a file named
// as "Fig2A_Output.txt". For each value of gh (which is printed first in
// the file), there are two columns of output: the induction frequency
// used and the % change in synaptic weight for that induction frequency.An
// example output, as generated by this function, is also attached with the
// same filename "Fig2A_Output.txt".
proc Fig2A() {
gh=0.0042
weight= new Vector()
double isilist[5000]
isicount=0
for (i=0.5;i<25.5; i+=0.5) {
isilist[isicount]=1000/i
isicount=isicount+1
}
ampa.winit=0.25
w=0.25
wopen("Fig2A_Output.txt")
i_init=0
for (gh=0.00; gh<=0.064; gh*=2) { // For each value of Ih density
fprint ("%f\n", gh)
print gh
update_Gh()
for (i=i_init; i<isicount; i+=1) { // For each ISI in the list
s.interval=isilist[i]
weight.record(&a.w)
update_init()
tstop=s.interval*900 // 900 pulses of that ISI
while (t < tstop){
fadvance()
}
sum=0
right=weight.size()-1000
if (right<0) right=0
count=0
for (j=weight.size()-1; j>right; j-=1){
sum=sum+weight.x[j]
count=count+1
}
sum=sum/count
fprint("%f\t%f\n",1000/isilist[i],100*(sum/w-1))
print 1000/isilist[i], "\t", sum, "\t", 100*(sum/w-1)
}
if (gh==0.0) {
gh=0.00025
}
}
wopen()
}
/********************************************************************/
// Computes the sliding modification threshold for each value of gh and
// saves in a file called "Fig2B_Output.txt". The output file has two
// columns: the first is the gh, and the second is the sliding modification
// threshold. An example output, as generated by this function, is also
// attached with the same filename "Fig2B_Output.txt".
proc Fig2B() {
weight= new Vector()
double isilist[1000]
isicount=0
for (i=0.2;i<25; i+=0.05) {
isilist[isicount]=1000/i
isicount=isicount+1
}
ampa.winit=0.25
w=0.25
wopen("Fig2B_Output.txt")
i_init=0
for (gh=0.0; gh<=0.012; gh+=0.00025) { // For each value of Ih density
print gh
update_Gh()
for (i=i_init; i<isicount; i+=1) { // For each ISI in the list
s.interval=isilist[i]
weight.record(&a.w)
update_init()
tstop=s.interval*900 // 900 pulses of that ISI
while (t < tstop){
fadvance()
}
sum=0
right=weight.size()-1000
if (right<0) right=0
count=0
for (j=weight.size()-1; j>right; j-=1){
sum=sum+weight.x[j]
count=count+1
}
sum=sum/count
print 1000/isilist[i], "\t", sum, "\t", sum/w
if (sum/w > 1){
fprint ("%f\t%f\n", gh, 1000/isilist[i])
break // Break if LTP threshold has been reached.
}
// Increasing gh shifts the sliding modification threshold to the right, so
// the following step to have i_init to start from 10 steps behind the
// current i is a time-saving exercise, so that i_init does not have to
// start from 0 always.
i_init=i-10
if (i_init<0) i_init=0
}
}
wopen()
}
/********************************************************************/
// Run Fig2A () for obtaining data for Fig. 2A of the paper.
// Run Fig2B () for obtaining data for Fig. 2B of the paper.
xpanel("Fig 2")
xbutton("Run Fig2A", "Fig2A()")
xbutton("Run Fig2B", "Fig2B()")
xlabel("Select one of the two simulations above")
xpanel()
//Fig2A()
//Fig2B()