RMD.ode

#Nicoletti M, Loppini A, Chiodo L, Folli V, Ruocco G, et al. (2019) 
#Biophysical modeling of C. elegans neurons: Single ion currents and whole-cell dynamics of AWCon and RMD. 
#PLOS ONE 14(7): e0218738. https://doi.org/10.1371/journal.pone.0218738. PMID: 31260485


#=========================== RMD simulation =====================================================================

# potassium,calcium, leakge ad sodium reversal potentials
par ek=-80,eca=60,eleak=-80,ena=30

#===============================================================================================================


#                               VOLTAGE-GATED POTASSIUM CHANNELS 



# ==============================================================================================================

#==================== SHL-1 CHANNELS ===========================================================================
#conductance in nS
par gshal=2.48

par shalsfhit=18

# activation
par vashal=11.2,kashal=14.1,kishal=8.3,vishal=-33.1
par ptmshal1=13.8,ptmshal2=-17.5165,ptmshal3=12.9213,ptmshal4=-3.7082,ptmshal5=6.4876,ptmshal6=1.8849
par cashal=0.1

minf_shal=1/(1+exp(-(v-vashal+shalsfhit)/(kashal)))
tm_shal=(ptmshal1/(exp(-(v-ptmshal2)/ptmshal3)+exp((v-ptmshal4)/ptmshal5))+ptmshal6)*cashal

#inactivation
par pthfshal1=539.1584, pthfshal2=-28.1990, pthfshal3=4.9199, pthfshal4=27.2811
par pthsshal1=8422 pthsshal2=-37.7391 pthsshal3=6.3785 pthsshal4=118.8983,cshal=0.1

hinf_shal=1/(1+exp((v-vishal+shalsfhit)/(kishal)))
ths_shal=(pthsshal1/(1+exp((v-pthsshal2)/pthsshal3))+pthsshal4)*cshal
thf_shal=(pthfshal1/(1+exp((v-pthfshal2)/pthfshal3))+pthfshal4)*cshal


m_shal'=(minf_shal-m_shal)/tm_shal
hf_shal'=(hinf_shal-hf_shal)/thf_shal
hs_shal'=(hinf_shal-hs_shal)/ths_shal
init m_shal=0,hf_shal=1,hs_shal=1

Ishal=gshal*(m_shal^3)*(0.7*hf_shal+0.3*hs_shal)*(v-ek)

#=================== SHK1 CHANNELS ==============================================================================

#conductance
par gshak=1.1

#activation
par vashak=20.4,kashak=7.7
par ptmshak1=26.571450568169027 ,ptmshak2=-33.741611800716130,ptmshak3=15.757936311607475 
par ptmshak4=15.364937728953288,ptmshak5=1.990037272604829 ,shiftV05=0

minf_shak=1/(1+exp(-(v-vashak+shiftV05)/kashak))
tm_shak=ptmshak1/(exp(-(v-(ptmshak2+shiftV05))/ptmshak4)+exp((v-(ptmshak2+shiftV05))/ptmshak3))+ptmshak5

#inactivation
par kishak=5.8,vishak=-6.95
par pthshak=1400

hinf_shak=1/(1+exp((v-vishak+shiftV05)/kishak))
th_shak=pthshak

m_shak'=(minf_shak-m_shak)/tm_shak
h_shak'=(hinf_shak-h_shak)/th_shak

init m_shak=0,h_shak=1

I_shak=gshak*m_shak*h_shak*(v-ek)


#================  EGL-36 CHANNELS ==============================================================================

#conductance
par gegl36=1.3

#activation
par va_egl36=63,ka_egl36=28.5,t1_egl36=355,t2_egl36=63,t3_egl36=13
par a1=0.31,a2=0.36,a3=0.39

min1_egl36=1/(1+exp(-(v-va_egl36)/(ka_egl36)))
min2_egl36=1/(1+exp(-(v-va_egl36)/(ka_egl36)))
min3_egl36=1/(1+exp(-(v-va_egl36)/(ka_egl36)))

tm1_egl36=t1_egl36
tm2_egl36=t2_egl36
tm3_egl36=t3_egl36

m1_egl36'=(min1_egl36-m1_egl36)/tm1_egl36
m2_egl36'=(min2_egl36-m2_egl36)/tm2_egl36
m3_egl36'=(min3_egl36-m3_egl36)/tm3_egl36

init m1_egl36=0,m2_egl36=0,m3_egl36=0

I_egl36=gegl36*((a1*m1_egl36)+(a2*m2_egl36)+(a3*m3_egl36))*(v-ek)

#======================= IRK CHANNELS =================================================================================

par gkir=0.2

par va_kir=-52, ka_kir=13 
par p1tmkir=17.0752,p2tmkir=-17.8258, p3tmkir=20.3154, p4tmkir=-43.4414, p5tmkir=11.1691, p6tmkir=3.8329

minf_kir=1/(1+exp((v-va_kir+30)/ka_kir))

tm_kir=p1tmkir/(exp(-(v-p2tmkir)/p3tmkir)+exp((v-p4tmkir)/p5tmkir))+p6tmkir

m_kir'=(minf_kir-m_kir)/tm_kir

I_kir=gkir*m_kir*(v-ek)
aux I_kir=I_kir


#=========================================================================================================================


#                                      CALCIUM CHANNELS 


#=========================================================================================================================


#=================== UNC-2 CHANNELS ======================================================================================= 
#conductance
par gunc2=0.9

#activation
par va_unc2=-12.17,ka_unc2=3.97,stm2=25
par p1tmunc2=1.4969,p2tmunc2=-8.1761,p3tmunc2=9.0753,p4tmunc2=15.3456,p5tmunc2=0.1029
par shiftmunc2=30,constmunc2=3

minf_unc2=1/(1+exp(-(v-va_unc2+stm2)/(ka_unc2)))
tm_unc2=(p1tmunc2/(exp(-(v-p2tmunc2+shiftmunc2)/(p3tmunc2))+exp((v-p2tmunc2+shiftmunc2)/(p4tmunc2)))+p5tmunc2)*constmunc2

#inactivation
par vi_unc2=-52.47,ki_unc2=5.6,sth2=25
par p1thunc2=83.8037, p2thunc2=52.8997,p3thunc2=3.4557,p4thunc2=72.0995,p5thunc2=23.9009,p6thunc2=3.5903
par consthunc2=1.7, shifthunc2=30


hinf_unc2=1/(1+exp((v-vi_unc2+sth2)/(ki_unc2)))
th_unc2=(p1thunc2/(1+exp((v-p2thunc2+shifthunc2)/(p3thunc2)))+p4thunc2/(1+exp(-(v-p5thunc2+shifthunc2)/(p6thunc2))))*consthunc2

# equations
m_unc2'=(minf_unc2-m_unc2)/tm_unc2
h_unc2'=(hinf_unc2-h_unc2)/th_unc2

init m_unc2=0,h_unc2=1

I_unc2=gunc2*m_unc2*h_unc2*(v-eca)

#alpha and beta for slo1 ans slo2 channels 
alpha=minf_unc2/tm_unc2
beta=(1/tm_unc2)-alpha


#================== EGL-19 CHANNELS ===========================================================================================
#conductance
par gegl19=0.99

#activation
par va_egl19=5.6, ka_egl19=7.50,stm19=10,sth19=10
par pdg1=2.3359,pdg2=2.9324,pdg3=5.2357,pdg4=6.0,pdg5=1.8739,pdg6=1.3930,pdg7=30.0,stau19=10

minf_egl19=1/(1+exp(-(v-va_egl19+stm19)/ka_egl19))
tm_egl19=(pdg1+(pdg2*exp(-(v-pdg3+stau19)^2/(pdg4)^2))+(pdg5*exp(-(v-pdg6+stau19)^2/(pdg7)^2)))

#inactivation
par p1hegl19=1.4314,p2hegl19=24.8573,p3hegl19=11.9541,p4hegl19=0.1427,p5hegl19=5.9589,p6hegl19=-10.5428,p7hegl19=8.0552,p8hegl19=0.6038
par pds1=0.4,pds2=0.55,pds3=81.1179,pds4=-22.9723,pds5=5,pds6=43.0937,pds7=0.9,pds8=40.4885,pds9=28.7251,pds10=3.7125,pds11=0,shiftdps=10

hinf_egl19=((p1hegl19/(1+exp(-(v-p2hegl19+sth19)/p3hegl19))+p4hegl19)*(p5hegl19/(1+exp((v-p6hegl19+sth19)/p7hegl19))+p8hegl19))
ths_egl19=pds1*(((pds2*pds3)/(1+exp((v-pds4+shiftdps)/pds5)))+pds6+((pds7*pds8)/(1+exp((v-pds9+shiftdps)/pds10)))+pds11)

# equations
m_egl19'=(minf_egl19-m_egl19)/tm_egl19
hs_egl19'=(hinf_egl19-hs_egl19)/ths_egl19

init m_egl19=0,hs_egl19=1

I_egl19=gegl19*m_egl19*hs_egl19*(v-eca)

# alpha and beta for slo1 and slo2  channels
alpha1=minf_egl19/tm_egl19
beta1=(1/tm_egl19)-alpha1

#===============  CCA-1 channels ==================================================================================================
#conductance
par gcca1=3.1

par  va_cca1=-42.65,ka_cca1=1.7,sscca1=15
par stmcca1=30,sthcca1=15,sshcca1=15,constmcca1=0.5
par fcca=1.4,f2cca1=1.15,f3ca=1.7,f4ca=1.1
par p1tmcca1=40,p2tmcca1=-62.5393,p3tmcca1=-12.4758,p4tmcca1=0.6947

minf_cca1=1/(1+exp(-(v-va_cca1+sscca1)/(ka_cca1*fcca)))
tm_cca1=((p1tmcca1/(1+exp(-(v-p2tmcca1+stmcca1)/(p3tmcca1*f3ca))))+p4tmcca1)*constmcca1

par vi_cca1=-58,ki_cca1=7,consthcca1=0.08
par p1thcca1=280,p2thcca1=-60.7312,p3thcca1=8.5224,p4thcca1=19.7456

hinf_cca1=1/(1+exp((v-vi_cca1+sshcca1)/(ki_cca1*f2cca1)))
th_cca1=((p1thcca1/(1+exp((v-p2thcca1+sthcca1)/(p3thcca1*f4ca))))+p4thcca1)*consthcca1

# equations
m_cca1'=(minf_cca1-m_cca1)/tm_cca1
h_cca1'=(hinf_cca1-h_cca1)/th_cca1

init m_cca1=0,h_cca1=1

I_cca1=gcca1*m_cca1^2*h_cca1*(v-eca)



#======================================================================================================================================


#                                      CALCIUM-REGULATED POTASSIUM CHANNELS 


#======================================================================================================================================


#================== Ca nanodomain calculation (Only for SLO1 and SLO2)=================================================================

# r: nanodomain radius
# d: diffusion coefficient
# F: Faraday constant C/mol
# kb:
# b:
# gsc: conductance of calcium channels

par r=13e-9,d=250e-12,F=96485,kb=500e6,b=30e-6,backgr=0.05,gsc=40e-12

cao_nano=(((abs(gsc*(v-eca)*1e-3)/(8*pi*r*d*F))*exp(-r/sqrt(d/(kb*b))))*1e6*1e-3 )+backgr
cac_nano=backgr

#================= SLO1 CHANNELS  ====================================================================================================
# no clustering 1:1  stoichiometry with EGL-19 and UNC-2 channels 

# activation 
# (parameters from genetic algorithm)
par wom=3.152961,wyx=0.012643,kyx=34.338784,nyx=0.000100
par wop=0.156217,wxy=-0.027527,kxy=55.726816,nxy=1.299198

kcm=wom*exp(-wyx*v)*(1/(1+((cac_nano)/kyx)^nyx))
kom=wom*exp(-wyx*v)*(1/(1+((cao_nano)/kyx)^nyx))
kop=wop*exp(-wxy*v)*(1/(1+(kxy/(cao_nano))^nxy))
kcp=0

# SLO1-UNC2 COMPLEX

#conductance 
par gbk=0.3

minf_bk=(m_unc2*kop*(alpha+beta+kcm))/((kop+kom)*(kcm+alpha)+(beta*kcm))
tm_bkunc2=((alpha+beta+kcm)/((kop+kom)*(kcm+alpha)+(beta*kcm)))

mbk'=(minf_bk-mbk)/tm_bkunc2
init mbk=0

I_bk=gbk*mbk*h_unc2*(v-ek)

# SLO1-EGL19 COMPLEX

#conductance
par gslo1=0.3

#activation
# same activation parameters of SLO1-UNC2 complex

kcm2=wom*exp(-wyx*v)*(1/(1+((cac_nano)/kyx)^nyx))
kom2=wom*exp(-wyx*v)*(1/(1+((cao_nano)/kyx)^nyx))
kop2=wop*exp(-wxy*v)*(1/(1+(kxy/(cao_nano))^nxy))
kcp2=0

minf_slo1=(m_egl19*kop*(alpha1+beta1+kcm2))/((kop2+kom2)*(kcm2+alpha1)+(beta1*kcm2))
tm_slo1=((alpha1+beta1+kcm2)/((kop2+kom2)*(kcm2+alpha1)+(beta1*kcm2)))

mslo1'=(minf_slo1-mslo1)/tm_slo1

init mslo1=0

I_slo1=gslo1*mslo1*hs_egl19*(v-ek)


#================ SLO-2 CHANNELS =============================================================================================
# no clustering, 1:1 stoichiometry with EGL-19 and UNC-2


# activation
# parameters from genetic algorithm

par wom1=0.896395,wyx1=0.019405,kyx1=3294.553404,nyx1=0.000010
par wop1=0.026719,wxy1=-0.024123,kxy1=93.449423,nxy1=1.835067

kcm1=wom1*exp(-wyx1*v)*(1/(1+(cac_nano/kyx1)^nyx1))
kom1=wom1*exp(-wyx1*v)*(1/(1+((cao_nano)/kyx1)^nyx1))
kop1=wop1*exp(-wxy1*v)*(1/(1+(kxy1/(cao_nano))^nxy1))
kcp1=0


# SLO2- EGL19 COMPLEX 

#conductance

par gbk2=0.3

minf_bk2=(m_egl19*kop1*(alpha1+beta1+kcm1))/((kop1+kom1)*(kcm1+alpha1)+(beta1*kcm1))
tm_bkegl19=((alpha1+beta1+kcm1)/((kop1+kom1)*(kcm1+alpha1)+(beta1*kcm1)))

mbk2'=(minf_bk2-mbk2)/tm_bkegl19

init mbk2=0

I_bk2=gbk2*mbk2*hs_egl19*(v-ek)


# SLO2-UNC2 COMPLEX

#conductance
par gslo2=0.3

kcm3=wom1*exp(-wyx1*v)*(1/(1+(cac_nano/kyx1)^nyx1))
kom3=wom1*exp(-wyx1*v)*(1/(1+((cao_nano)/kyx1)^nyx1))
kop3=wop1*exp(-wxy1*v)*(1/(1+(kxy1/(cao_nano))^nxy1))
kcp3=0

minf_slo2=(m_unc2*kop3*(alpha+beta+kcm3))/((kop3+kom3)*(kcm3+alpha)+(beta*kcm3))
tm_slo2=((alpha+beta+kcm3)/((kop3+kom3)*(kcm3+alpha)+(beta*kcm3)))

mslo2'=(minf_slo2-mslo2)/tm_slo2



I_slo2=gslo2*mslo2*h_unc2*(v-ek)

													
#================================= Intracellular calcium calculation =========================================================
# fd: Faraday constant in C/mmol

par fd=96.485,fca=0.001,t_ca=50,rcell=1
I_ca=I_unc2+I_egl19+I_cca1
aux I_ca=I_ca

# RMD total volume in um^3 (from Neuromorpho)
par vol=5.65


alpha_ca=1/(2*vol*fd)

J_ca1=fca*alpha_ca*I_ca
aux J_ca1=J_ca1

backgr2=0.05e-3

rs=if(I_ca<0)then(-fca*(alpha_ca*I_ca)-((ca_intra1-backgr2)/t_ca))else((backgr2-ca_intra1)/t_ca)
ca_intra1'=rs

init ca_intra1=5e-05


#============================= KCNL CHANNELS  ==================================================================================
# Rinzel 

#conductance
par gca=0.06

#activation
par k_sk2=0.00033

minf_sk=ca_intra1/(k_sk2+ca_intra1)
t_sk=6.3

#equations
m_sk'=(minf_sk-m_sk)/t_sk

init m_sk=0.13563


I_sk=gca*m_sk*(v-ek)



#==================================================================================================================================


#                                      LEAKAGE CHANNELS 


#===================================================================================================================================



#================= LEAKAGE CURRENT =================================================================================================
par gleak=0.4
I_leak=gleak*(v-eleak)

#================== Nca current ====================================================================================================
par gnca=0.05
I_nca=gnca*(v-ena)



#===================================================================================================================================


#                                                      SIMULATIONS


#====================================================================================================================================


# Total current


I_tot=Ishal+I_bk+I_bk2+I_shak+I_cca1+I_egl19+I_unc2+I_sk+I_leak+I_egl36+I_nca+I_slo1+I_slo2+I_kir

aux Itot=I_tot


#===================== CURRENT CLAMP ================================================================================================

# IH: holding current
# c: membrane capacitance pF

par ih=0
par c=1.2


v'=(Istim-I_tot)/c
init v=-70
#current clamp protocol
Istim=if(t>ton)then(if(t<toff)then(iclamp)else(if(t>t2)then(if(t<t3)then(iclamp2)else(ih))else(ih)))else(ih)
par iclamp=10,ton=310,toff=360,t2=410,iclamp2=-15,t3=430
aux prot=if(t>ton)then(if(t<toff)then(iclamp)else(if(t>t2)then(if(t<t3)then(iclamp2)else(ih))else(ih)))else(ih)



### XPP SETTINGS
@ meth=stiff,trans=200,total=400,dt=0.01,maxstor=1000000,bound=100000000000000,atol=1e-8,tol=1e-8
@ xp=t,yp=v,xlo=200,xhi=1000,yhi=30,ylo=-100
done




#===================== VOLTAGE CLAMP =================================================================================================

# voltage clamp protocol

#par vk=-70

#v=if(t>ton)then(if(t<toff)then(vclamp)else(vk))else(vk)
#par vclamp=1100,ton=1100,toff=2300

### XPP SETTINGS
#@ meth=stiff,trans=1000,total=2400,dt=0.01,maxstor=1000000,bound=100000000000000,atol=1e-8,tol=1e-8
#@ xp=t,yp=Itot,xlo=1000,xhi=2400,yhi=800,ylo=-100
#done