BPF.hoc

/* Implementation of Short-term plasticity at CA3-CA1 synapses in a conductance 
 * based model in the paper:
 * 
 * Mukunda, C. L. and Narayanan, R. (2017), Degeneracy in the regulation of 
 * short‐term 
 * plasticity and synaptic filtering by presynaptic mechanisms. J Physiol, 595: 
 * 2611-2637. doi:10.1113/JP273482
 * 
 * The model contains Na, KA, KDr, CaN, CaL and h currents in the terminal, 
 * while other compartments have Na and KDr currents only. AMPA receptor in the CA1 spine is
 * implemented using GHK equations. A detailed reaction-diffusion model for 
 * calcium in the presynaptic terminal, with calcium ON and OFF mechanisms and radial 
 * diffusion has been adopted from Ashhad & Narayanan (2013). The model is tuned to obtain 
 * short-term plasticity profile as observed by Dittman et al. (J Neuroscience, 2000).
 * 
 * This code generates data for Figure 1. of the paper.
 * 
 * Implemented by Chinmayee L Mukunda. Contact: chinmayeelm@gmail.com
 *  
 */

load_file("stdlib.hoc")
//load_file("nrngui.hoc")

file_num = 1

nms=50
nnor=50
celsius     =   34



create soma, AH, AIS, MS[nms], NOR[nnor], terminal
access soma

objref serialsec, NORs, MSs
serialsec = new SectionList()
NORs      = new SectionList()
MSs       = new SectionList()

/***********************************/
/*Parameters*/
/**********************************/


// --------------------------------------------------------------
// Resistances and Capacitances
// --------------------------------------------------------------
Rm          =   40000
Cm          =   1
RA          =   120
RmMS        =   Rm*2
CmMS        =   Cm/2

//--------------------------------------------------------------
// Active conductance densities and other related parameters
//--------------------------------------------------------------
gNaBar      =   0.053866
gKdrBar     =   0.016384
gCalBar     =   0.000938
gCanBar     =   0.0008
gHdBar      =   0.000188
gKaBar      =   0.015907
gNaBarNOR   =   18*gNaBar
gKdrBarNOR  =   56.667*gKdrBar
v_init      =   -65         //Resting membrane potential

//--------------------------------------------------------------
// Parameters related to Calcium stores
//--------------------------------------------------------------
Gamma       =   196916014533.504303
Alpha       =   .9
Beta        =   1.3
Cathresh    =   2e-4
VmaxSERCA   =   4e-4
Fca         =   5
kds         =   4
TotBufs     =   0.522645
TotBufm     =   0
Jmax        =   3.5e-3

//--------------------------------------------------------------
//Neurotransmitter kinetics related parameters
//--------------------------------------------------------------

ttinf       =   0.0001
TauTT       =   125
TauTC       =   0.0001
ft          =   11e9

P           =   1e-7   //Permeability of AMPA receptor in cm/s

//-------------------------------------------------------------
//Physical size parameters of MS and NOR compartments
//-------------------------------------------------------------
diameterNOR =   1
lengthNOR   =   4
lengthMS    =   200
diameterMS  =   2



//--------------------------------------------------------------
// The code.
//--------------------------------------------------------------

soma {
    nseg = 1
    diam = 100
    L    = 100
    Ra   = 120
    cm   = Cm
    
    //Initialize Passive elements
    insert pas
    g_pas=1/Rm
    e_pas=v_init
    
    //Initialize Active elements
    insert nahh
    gnabar_nahh=gNaBar
    
    insert borgkdr
    gkdrbar_borgkdr=gKdrBar
    
    
    insert minmax
    serialsec.append()
    
    ena         =   55
    ek          =   -90
}

AH {
    diam = 3
    L = 15
    Ra = 120
    cm=Cm
    
    insert pas
    g_pas=1/Rm
    e_pas=v_init
    nseg=int((L/(0.1*lambda_f(100))+0.9)/2)*2+1
    
    insert nahh
    gnabar_nahh= 5*gNaBar
    
    insert borgkdr
    gkdrbar_borgkdr= gKdrBar
    
    insert minmax
    serialsec.append()
    
    ena         =   55
    ek          =   -90
    
}

AIS {
    diam = 2//1.5
    L = 25
    Ra = 120
    cm=Cm
    
    insert pas
    g_pas=1/Rm
    e_pas=v_init
    nseg=int((L/(0.1*lambda_f(100))+0.9)/2)*2+1
    
    insert nahh
    gnabar_nahh = 16*gNaBar
    
    
    
    insert borgkdr
    gkdrbar_borgkdr = gKdrBar
    
    
    insert minmax
    serialsec.append()
    
    ena         =   55
    ek          =   -90
    
}

for i=0,49 MS[i] {
    diam = diameterMS
    L = lengthMS
    Ra = 120
    cm= CmMS
    nseg=int((L/(0.1*lambda_f(100))+0.9)/2)*2+1
    
    insert pas
    g_pas= 1/RmMS
    e_pas=v_init
    
    insert minmax
    MSs.append()
    

}

for i=0,49 NOR[i] {
    diam = diameterNOR
    L = lengthNOR
    Ra = 120
    cm=Cm
    nseg=int((L/(0.1*lambda_f(100))+0.9)/2)*2+1
    
    insert pas
    g_pas=1/Rm
    e_pas=v_init
    
    insert nahh
    gnabar_nahh= gNaBarNOR
    
    insert borgkdr
    gkdrbar_borgkdr= gKdrBarNOR
    
    
    insert minmax
    NORs.append()
    
    ena         =   55
    ek          =   -90
}

for i=0,49 {
    
    NOR[i] serialsec.append()
}

terminal {
    
    nseg=1
    diam = 2
    L = 0.5
    Ra = 120
    cm=Cm
    
    insert pas
    g_pas=1/Rm
    e_pas=v_init
    
    insert nahh
    gnabar_nahh= 10*gNaBar
    
    insert borgkdr
    gkdrbar_borgkdr= gKdrBar
    
    
    insert cal
    gcalbar_cal=gCalBar
    
    insert can
    gcanbar_can=gCanBar
    
    insert hd
    ghdbar_hd = gHdBar
    
    insert kad
    gkabar_kad = gKaBar
    
    insert minmax
    
    
    insert cal4
    gamma_cal4=Gamma
    cath_cal4=Cathresh
    vmax_cal4=VmaxSERCA
    TBufs_cal4=TotBufs
    TBufm_cal4=TotBufm
    alpha_cal4=Alpha
    beta_cal4=Beta
    KDs_cal4=kds
    factor_cal4 = Fca
    jmax_cal4   = Jmax
    
    insert stp
    
    fT_stp = ft
    tauTC_stp = TauTC
    TTINF_stp = ttinf
    tauTT_stp = TauTT
    
    
    ena         =   55
    ek          =   -90
    
}

//Connecting the compartments to build neuron
connect soma(1), AH(0)
connect AH(1), AIS(0)
connect AIS(1), MS[0](0)



for(i=0;i<49;i+=1){
    connect MS[i](1), NOR[i](0)
    connect NOR[i](1), MS[i+1](0)
}

connect MS[49](1), NOR[49](0)
connect NOR[49](1), terminal(0)

access soma

stamp=2
objectvar stim
soma stim = new Pulses(0.5)

//Parameters of the stimulation pulses

stim.del = 500
stim.dur = 2
stim.amp = stamp
stim.npulses=10
stim.period=20

//tstop = 300

//Useful only during initial tuning of certain parameters for efficient propagation of spikes
/*  
access NOR[0]
newPlotV()

proc update_MS(){
    forsec MSs{
        L = lengthMS
        diam = diameterMS
        cm = CmMS
        g_pas= 1/RmMS
        nseg=int((L/(0.1*lambda_f(100))+0.9)/2)*2+1
        
    }
}

proc update_NOR(){
    forsec NORs{
        L = lengthNOR
        diam = diameterNOR
        gnabar_nahh = gNaBarNOR
        gkdrbar_borgkdr = gKdrBarNOR
        nseg=int((L/(0.1*lambda_f(100))+0.9)/2)*2+1
        
    }
}

proc update_RA(){
    forsec serialsec{
        Ra = RA
        nseg=int((L/(0.1*lambda_f(100))+0.9)/2)*2+1
        
    }
    
    forsec MSs{
        Ra = RA
        nseg=int((L/(0.1*lambda_f(100))+0.9)/2)*2+1
    }
    
    terminal{
        Ra = RA
        nseg=int((L/(0.1*lambda_f(100))+0.9)/2)*2+1
    }
    
}

proc update_RM(){
    forsec serialsec{
        g_pas = 1/Rm
        nseg=int((L/(0.1*lambda_f(100))+0.9)/2)*2+1
        
    }
    
    forsec MSs{
        g_pas = 1/(Rm*2)
        nseg=int((L/(0.1*lambda_f(100))+0.9)/2)*2+1
    }
    
    terminal{
        g_pas = 1/Rm
        nseg=int((L/(0.1*lambda_f(100))+0.9)/2)*2+1
    }
    
}
*/

/****************************************************************/
proc update_init(){
	finitialize(v_init)
	fcurrent()
	forsec serialsec {
		e_pas=v_init
		e_pas=e_pas+(ina+ik)/g_pas
	}
    
}
/***************************************************************/

proc update_init_MS(){
	finitialize(v_init)
	fcurrent()
	forsec MSs {
		e_pas=v_init
	}
    
}
/***************************************************************/

proc update_init_terminal(){
	finitialize(v_init)
	fcurrent()
	terminal{
        e_pas=v_init
    }
}
/***************************************************************/

/*
xpanel("Parameters")
xvalue("gnanor","gNaBarNOR", 1, "update_NOR()",1,1)
xvalue("gkdrnor","gKdrBarNOR", 1, "update_NOR()",1,1)
xvalue("lennor","lengthNOR", 1, "update_NOR()",1,1)
xvalue("dianor","diameterNOR", 1, "update_NOR()",1,1)
xvalue("lenms","lengthMS", 1, "update_MS()",1,1)
xvalue("diams","diameterMS", 1, "update_MS()",1,1)
xvalue("CM","CmMS", 1, "update_MS()",1,1)
xvalue("RM","RmMS", 1, "update_MS()",1,1)
xvalue("Raxial","RA", 1, "update_RA()",1,1)
xvalue("Rmem","Rm",1, "update_RM()",1,1)
xpanel()

update_init()
update_init_terminal()
update_init_MS()

init()
run()

nrncontrolmenu()

*/

create POST

POST{
    nseg=1
    diam = 2
    L = 0.5
    Ra = 120
    cm=Cm
    
    insert pas
    g_pas=1/Rm
    e_pas=v_init
}

access POST


//Inserting voltage clamp to obtain current fluctuations
objref vclamp
POST vclamp= new SEClamp(0.5)

vclamp.dur1=10
vclamp.dur3=10
vclamp.amp1=-70
vclamp.amp2=-70
vclamp.amp3=-70

objectvar c
c = new ghkampaC()      //creating AMPA receptor point process
POST c.loc(0.5)
c.Pmax=P



objref EPSCvec,CAvec, TTvec, TCvec, APSvec, APTvec,INAvec, IKvec
objref EPSCfile, APSfile,APTfile, CAfile, TCfile,TTfile, VALfile

EPSCvec=new Vector()
CAvec=new Vector()
TTvec=new Vector()
TCvec=new Vector()
APSvec=new Vector()
APTvec=new Vector()
INAvec=new Vector()
IKvec=new Vector()

EPSCfile = new File()
APSfile = new File()
APTfile = new File()
CAfile = new File()
TTfile = new File()
TCfile = new File()
VALfile = new File()



strdef EPSC, APS, APT, CA, TT, TC, VAL



sprint(EPSC, "data/EPSC_%d.txt", file_num)
EPSCfile.wopen(EPSC)

sprint(VAL, "data/VAL_%d.txt", file_num)
VALfile.wopen(VAL)

sprint(APS, "data/APS_%d.txt", file_num)
APSfile.wopen(APS)

sprint(APT, "data/APT_%d.txt", file_num)
APTfile.wopen(APT)

sprint(CA, "data/CA_%d.txt", file_num)
CAfile.wopen(CA)

sprint(TT,"data/TT_%d.txt", file_num)
TTfile.wopen(TT)

sprint(TC, "data/TC_%d.txt", file_num )
TCfile.wopen(TC)

for(k=1; k<=50; k+=1){
    
    period = 1000/k
    access soma
    stim.period = period
    

    
    APSvec.record(&soma.v(0.5))
    CAvec.record(&terminal.cai(0.5))
    APTvec.record(&terminal.v(0.5))
    TTvec.record(&terminal.TT_stp(0.5))
    TCvec.record(&terminal.TC_stp(0.5))
    INAvec.record(&POST.ina(0.5))
    IKvec.record(&POST.ik(0.5))

    
    tstop= 500 + 10*period + 500

    
    access POST
    vclamp.dur2= tstop
    
    update_init()
    update_init_terminal()
    update_init_MS()
    
    finitialize()
    while(t<tstop){
        fadvance()
        c.C = terminal.TC_stp(0.5)
    }
    
    EPSCvec=INAvec.c.add(IKvec)
    
    ap  = APTvec.max(20000,21000) - v_init
    cal_spike = CAvec.max(20000,21000)
    epsc = EPSCvec.min(20000,21000)
    
    
    EPSCfile.aopen(EPSC)
    APSfile.aopen(APS)
    APTfile.aopen(APT)
    CAfile.aopen(CA)
    TCfile.aopen(TC)
    TTfile.aopen(TT)
    
    
    APSvec.printf(APSfile)
    APTvec.printf(APTfile)
    CAvec.printf(CAfile)
    EPSCvec.printf(EPSCfile)
    TCvec.printf(TCfile)
    TTvec.printf(TTfile)
    
    
    EPSCfile.close()
    APSfile.close()
    APTfile.close()
    CAfile.close()
    TCfile.close()
    TTfile.close()

    VALfile.aopen(VAL)
    VALfile.printf("%f\n%e\n%e\n",ap,cal_spike,epsc)
    VALfile.close()


}