# 2-Compartment Motoneuron Model - Sharmila Venugopal (svenugopal10@gmail.com)
# Cite: S Venugopal et al., J Neurophysiology, 2011

# Applied current pulse from holding current i0 with magnitude ip
#    turns on at pon and off at poff
# Ca concentration in micromoles/liter
#
# version for computing bifurcation diagrams
#
#soma equations
dvs/dt=(iapp+ins(vs,snah,scam,scah)+outs(vs,sn,sca)-gl*(vs-(vl+vrest))+gc*(vd-vs)/parea)/cm
dsnah/dt=(snahinf(vs)-snah)/snahtau(vs)
dsn/dt=(sninf(vs)-sn)/sntau(vs)
dscam/dt=(scaminf(vs)-scam)/camtau
dscah/dt=(scahinf(vs)-scah)/cahtau
dsca/dt=-f*alpha*sica(vs,scam,scah)-f*kca*sca

ds/dt=alphamn*(1-s)*sinf(vs)-betamn*s

#dendrite equations
#dvd/dt=(ind(vd,dcam,dcah,dcas)+outd(vd,dn,dca)-INaP-gl*(vd-(vl+vrest))+gc*(vs-vd)/(1.0-parea))/cm
dvd/dt=(ind(vd,dcas)+outd(vd,dn,dca)-pernap*gNaP*mnap*hnap*(Vd-(vna+vrest))-gl*(vd-(vl+vrest))-Ir+gc*(vs-vd)/(1.0-parea))/cm
#ddcam/dt=(dcaminf(vd)-dcam)/camtau
#ddcah/dt=(dcahinf(vd)-dcah)/cahtau
ddcas/dt=(dcasinf(vd)-dcas)/dcastau
ddn/dt=(dninf(vd)-dn)/dntau(vd)
#ddca/dt=-f*alpha*(dicas(vd,dcas)+dicaf(vd,dcam,dcah))-f*kca*dca
ddca/dt=-f*alpha2*(dicas(vd,dcas))-f*kca*dca

#-----------Persistent Sodium-----------------
#dmnap/dt=(minfnap(Vd)-mnap)/taumnap(Vd)
dmnap/dt=(minfnap(Vd)-mnap)/taunap
p taunap=40

#dhnap/dt=(hinfnap(vd)-hnap)/tauhnap(Vd)
dhnap/dt=(hinfnap(vd)-hnap)/naptau
p naptau=1000
#INaP=gNaP*minfnap(Vd)*hnap*(Vd-(vna+vrest))
#INaP=gNaP*mnap*hnap*(Vd-(vna+vrest))
#INaP=gNaP*mnap*(Vd-55)

minfnap(Vd)=1/(1+exp(-(vd-thetamnap)/Kmnap))
hinfnap(vd)=1/(1+exp((vd-thetahnap)/Khnap))
#taunap(Vd)=taumxnap/cosh((vd-thetahnap)/(2*Khnap))
taumnap(Vd)=Tmnap/cosh((vd-thetamnap)/(2*Khnap))
tauhnap(Vd)=taumxnap/cosh((vd-thetahnap)/(2*Khnap))

p thetamnap=-48,Kmnap=3,thetahnap=-35,Khnap=6,taumxnap=10000,Tmnap=5000
#
#---------------RC equations-----------------
dvr/dt=(-gna*minf(vr)^3*h(nr)*(vr-ena)-gk*nr^4*(vr-ek)-gl*(vr-el)+Idr)/cr
dnr/dt=(ninf(vr)-nr)/ntau(nr)


dsr/dt=alphar*(1-sr)*sinf(vr)-betar*sr

sinf(v)=1/(1+exp(-(v-thetas)/ks))
ntau(v)=28./(exp((v+50)/20)+exp(-(v+60)/10))
ninf(v)=1.0/(1+exp((nVh-v)/nk))
minf(v)=1.0/(1+exp((mVh-v)/mk))
h(n) = .7 - 1.1*n
#--------Synaptic Currents------------
Imn=gmn*s*(vr-ve)
p Ipool=0
p Idr=6.2

Ir=gr*sr*(vd-vi)

# ----------RC parameters---------------
p cr=1
p alphamn=2,betamn=0.05
p alphar=1
p betar=0.09
# other values of betar = 0.025, 0.0125
p thetas=-15, ks=.02
#p ggb=0.005,ggl=0.005
p gr=0.01	# Changed on Oct 6th 2009
p gmn=0.15
p ve=50,vi=-80
p rgk=100
p gna=80.0,rnamth=-35.0,rnamslp=7.8,rnahth=-55.0,rnahslp=7.0
par mVh=-40,mk=7,nVh=-45,nk=15
p rgna=70,el=-54,Gk=40
p Ena=50,Ek=-77,eca=80
# soma functions
ins(vs,snah,scam,scah)=-sina(vs,snah)-sica(vs,scam,scah)
outs(vs,sn,sca)=-sikdr(vs,sn)-sikca(vs,sca)
sina(vs,snah)=sgna*snaminf(vs)**3*snah*(vs-(vna+vrest))
sica(vs,scam,scah)=sgca*scam**2*scah*(vs-(vca+vrest))
sikdr(vs,sn)=sgk*sn**4*(vs-(vk+vrest))
sikca(vs,sca)=perkca*sgkca*(sca**p/(kd**p+sca**p))*(vs-(vk+vrest))
snaminf(vs)=1.0/(1.0+exp(-(vs-snamth)/snamslp))
snahinf(vs)=1.0/(1.0+exp((vs-snahth)/snahslp))
snahtau(vs)=snahc/(exp((vs-snahv)/snaha)+exp(-(vs-snahv)/snahb))
sninf(vs)=1.0/(1.0+exp(-(vs-nth)/nslp))
sntau(vs)=nc/(exp((vs-nv)/na)+exp(-(vs-nv)/nb))
scaminf(vs)=1.0/(1.0+exp(-(vs-camth)/camslp))
scahinf(vs)=1.0/(1.0+exp((vs-cahth)/cahslp))
# dendrite functions
ind(vd,dcas)=-dicas(vd,dcas)
outd(vd,dn,dca)=-dikdr(vd,dn)-dikca(vd,dca)
dicaf(vd,dcam,dcah)=dgcaf*dcam**2*dcah*(vd-(vca+vrest))
dicas(vd,dcas)=percap*dgcas*dcas*(vd-(vca+vrest))
dikdr(vd,dn)=dgk*dn**4*(vd-(vk+vrest))
dikca(vd,dca)=perkca*dgkca*(dca**p/(kd**p+dca**p))*(vd-(vk+vrest))
dcaminf(vd)=1.0/(1.0+exp(-(vd-camth)/camslp))
dcahinf(vd)=1.0/(1.0+exp((vd-cahth)/cahslp))
dcasinf(vd)=1.0/(1.0+exp(-(vd-dcasth)/dcasslp))
dninf(vd)=1.0/(1.0+exp(-(vd-nth)/nslp))
dntau(vd)=nc/(exp((vd-nv)/na)+exp(-(vd-nv)/nb))
# applied current - here a parameter iapp
#iapp(t)=i0+heav(poff-t)*heav(t-pon)*ip
#
# bifurcation parameter
p iapp=0
# common parameters
p vna=115.0,vk=-20.0,vl=0.0,vca=140.0,vrest=-60.0,cm=1.0,gl=0.51,
p p=1,kd=0.2,f=0.01,kca=2.0,alpha=0.009
# soma parameters
p sgna=80.0,snamth=-35.0,snamslp=7.8,snahth=-55.0,snahslp=7.0
p snahv=-50.0,snaha=15.0,snahb=16.0,snahc=30.0
p sgk=100.0,nth=-28.0,nslp=12.0,nv=-40.0,na=40.0,nb=50.0,nc=7.0
p sgca=14.0,camth=-30.0,camslp=5.0,camtau=4.0
p cahth=-45.0,cahslp=5.0,cahtau=40.0
#p sgkca=5.0
p sgkca=6
# dendrite parameters
p dgcaf=0.3
p dcasth=-39.0,dcasslp=7.0,dcastau=40.0
p dgk=0.0
# coupling parameters
p gc=0.1,parea=0.1
#p gc=0.08,parea=0.1
# current parameters
p i0=0,ip=10,pon=50,poff=1000
# conductance percentage factors
p perkca=1.0, perslca=1.0
# conductances sensitive to SCI
p dgkca=1.0
p gnap=0.1
p dgcas=0.25
p gpratio=1.772
p alpha2=0.009
p percap=1
p pernap=1

--------
Iramp=offsetr+scaler*(t-tonr)*(heav(t-tonr)*heav(toffr-t))+2*scaler*(tswitchr-t)*(heav(t-tswitchr)*heav(toffr-t)) 
par offsetr=0
par scaler=0.005
par tonr=0
par toffr=10000
par tswitchr=4000

# Changed on Nov 16 in order to compare the grading of sec. freq by RC inhibition
#Iramp=offsetr+scaler*(t-tonr)*(heav(t-tonr)*heav(toffr-t))
#par offsetr=0
#par scaler=0.0007
#par tonr=0
#par toffr=30000
#par tswitchr=4000
--------
# soma initial conditions
vs(0)=-57.34
snah(0)=0.5829
sn(0)=0.1239
scam(0)=0.004199
scah(0)=0.9219
sca(0)=0.0001406
# dendrite initial conditions
vd(0)=-56.64
#dcam(0)=0.00483
#dcah(0)=0.9112
dcas(0)=0.08493
dn(0)=0.1291
dca(0)=0.01724
#
#aux xiapp=iapp(t)
aux I=Iramp
aux Ipulse=Iapp
aux Ivclamp=(ins(vs,snah,scam,scah)+outs(vs,sn,sca)-gl*(vs-(vl+vrest))+gc*(vd-vs)/parea)/cm
aux Irc=-Ir
aux Ipic=ind(vd,dcas)-pernap*gNaP*mnap*hnap*(Vd-(vna+vrest))
#
#xpp formatting
@ TOTAL=10000,DT=0.05,
@ MAXSTOR=1000000,BOUNDS=600000,method=qualrk
@ xlo=0,xhi=10000,ylo=-60,yhi=40,nplot=2,yp2=vr,nplot=3,yp3=Ipulse
done