# Rotenone_ChemAut.ode
# 
# Huange et al. Multiple actions of rotenone, an inhibitor of mitochondrial respiratory chain, 
# on ionic currents and minature end-plate potential in mouse hippocampal (mHippoE-14) neurons.
# Cell Physiol Biochem 2018;47:330-343.
# Chemical autaptic regulation
# 
# tau=0.5
# init gaut=0
# dgaut/dt=0
para gaut=0.00005
para tau=5, kk=0
para kkk=1, theta=0, Vsyn=2
# n=5000
wiener b
# dv/dt=(I*100/(PI*2*rsoma*lsoma)-gna*h*(v-vna)*m^3-gkdr*nkdr^3*lkdr*(v+91) \
-gka*nka*lka*(v+91)-gkc*okc*(v+91)-gcat*0.001*nt32^2*lt32*(v-120)-gl*(v-vl)-gaut*(v-delay(v, tau)))/cm
dv/dt=(I*100/(PI*2*rsoma*lsoma)-gna*h*(v-vna)*m^3-gkdr*nkdr^3*lkdr*(v+91) \
-gka*nka*lka*(v+91)-gkc*okc*(v+91)-gcat*0.001*nt32^2*lt32*(v-120)-gl*(v-vl)-gaut*(v-Vsyn)*(1/(1+exp(-kkk*(delay(v,tau)-theta)))))/cm
ica(v) = gcat*0.001*nt32*nt32*lt32*(v-120)
# channel densities (Yilmaz and Ozer, Physica A, 2015 421:455-462)
para rhona=60, rhok=18, rhoca=10, nnna=100, nnk=80, nnca=50
# Na1.1 WT
dh/dt=(hinf(v)-h)/tauh(v)+kk*sqrt(((1-2*hinf(v))*h+hinf(v))/(rhona*nnna))*b

# Cat3.1, 3.2, 3.3
# dnt31/dt=(nt31inf(v)-nt31)/taun31(v)
# dlt31/dt=(lt31inf(v)-lt31)/taul31(v)
 dnt32/dt=(nt32inf(v)-nt32)/taun32(v)+kk*sqrt(((1-2*nt32inf(v))*nt32+nt32inf(v))/(rhoca*nnca))*b
 dlt32/dt=(lt32inf(v)-lt32)/taul32(v)+kk*sqrt(((1-2*lt32inf(v))*lt32+lt32inf(v))/(rhoca*nnca))*b

# dnt33/dt=(nt33inf(v)-nt33)/taun33(v)
# dlt33/dt=(lt33inf(v)-lt33)/taul33(v)
# K dr
dnkdr/dt=andr(v)*(1-nkdr)-bndr(v)*nkdr+kk*sqrt(((1-2*(andr(v)/(andr(v)+bndr(v))))*nkdr+(andr(v)/(andr(v)+bndr(v))))/(rhok*nnk))*b
dlkdr/dt=aldr(v)*(1-lkdr)-bldr*lkdr+kk*sqrt(((1-2*(aldr(v)/(aldr(v)+bldr)))*lkdr+(aldr(v)/(aldr(v)+bldr)))/(rhok*nnk))*b
# K A
dnka/dt= ana(v)*(1-nka)-bna(v)*nka+kk*sqrt(((1-2*(ana(v)/(ana(v)+bna(v))))*nka+(ana(v)/(ana(v)+bna(v))))/(rhok*nnk))*b

dlka/dt= ala(v)*(1-lka)-bla*lka+kk*sqrt(((1-2*(ala(v)/(ala(v)+bla)))*lka+(ala(v)/(ala(v)+bla)))/(rhok*nnk))*b

# K C
  dokc/dt= ao(v)*(1-okc)-bo(v)*okc+kk*sqrt(((1-2*(ao(v)/(ao(v)+bo(v))))*okc+(ao(v)/(ao(v)+bo(v))))/(rhok*nnk))*b
# [Ca2+] dynamics
 dca/dt = 0.1*(-0.14*ica(v)-ca/5) 
# parameters
# par nt32=0, okc=0
par vna=50,vk=-91,gcat=5.0, i=0.02,  vl=-65
par gna=2,gkc=0.0004,gkdr=0.08,gka=0.0001,gl=1.6667e-5,cm=1e-3
num rm=60000, ra=200e-4, F=96520,F1=96.4853, R=8.3134, Temp=303.14
par cao=2,cai=5e-5
par rsoma=15, lsoma=25.5
par qt=1.18, q10na=1.62,q10m=3
# 
# rate functions for ion channels
# SCN1A data from lossin et al., 2002
minf(v) =  1/(1+exp(-(v+26.4)/7.1))
m=minf(v)
hinf(v) =  1/(1+exp(-(v+67.5)/(-6.2)))
tauh(v) =  0.197+10.701/(1+exp((v+78.87)*0.0538))
#
# T-type Ca2+ channels
# CaT31
 nt31inf(v) = 1/(1+exp(-(v+49.3)/4.6))
 lt31inf(v) = 1/(1+exp((v+74.2)/5.5))
 taun31(v) = if(v>(-56))then((0.8+0.025*exp(-v/14.5))/qt)else((1.71*exp((v+120)/38))/qt)
 taul31(v) = (if(v>(-60))then((12.3+0.12*exp(-v/10.8))/qt)else(137/qt)) 
# CaT32
 nt32inf(v) = 1/(1+exp(-(v+48.4)/5.2))
 lt32inf(v) = 1/(1+exp((v+75.6)/6.2))
 taun32(v) = (if(v>=(-56))then((1.34+0.035*exp(-v/11.8))/qt)else((2.44*exp((v+120)/40))/qt))
 taul32(v) = (if(v>=(-60))then((18.3+0.005*exp(-v/6.2))/qt)else(500/qt))
# CaT33
 nt33inf(v) = 1/(1+exp(-(v+41.5)/6.2))
 lt33inf(v) = 1/(1+exp((v+69.8)/6.1))
 taun33(v) = (if(v>=(-60))then((7.2+0.02*exp(-v/14.7))/qt)else((0.875*exp((v+120)/41))/qt))
 taul33(v) = (if(v>=(-60))then((79.5+2.0*exp(-v/9.3))/qt)else(260/qt)) 
# K+  C
 ao(v) = ca*0.28/(ca+0.48e-3*exp(-2*0.84*v*F1/(R*Temp)))
 bo(v) = 0.48/(1+ca/(0.13e-6*exp(-2*v*F1/(R*Temp))))
# K+ DR
andr(v) = 0.03*exp(1e-3*2*(v+32)*F/(R*Temp))
bndr(v) = 0.03*exp(-1e-3*3*(v+32)*F/(R*Temp))
aldr(v) = 0.001*exp(-1e-3*2*(v+61)*F/(R*Temp))
bldr = 0.001
# K+ A
ana(v) = 0.02*exp(1e-3*1.8*(v+33.6)*F/(R*Temp))
bna(v) = 0.02*exp(-1e-3*1.2*(v+33.6)*F/(R*Temp))
ala(v) = 0.08*exp(-1e-3*4*(v+83)*F/(R*Temp))
bla = 0.08
#
init v=-65, ca=5e-5
init nkdr=0.00179769,lkdr=0.576006
init nka=0.026395,lka=0.0596601
 init okc=0.00337468
init h=0.9793
 init nt32=0.0394562
 init lt32=0.153206
# ina=gna*h*(v-vna)*m^3
# ikc=gkc*okc*(v+91)
# aux ina=ina
# aux ikc=ikc
# init nt31=0.0318903,lt31=0.158061,nt32=0.0394562,lt32=0.153206,nt33=0.0220894,lt33=0.312838
@ dt=0.001, bound=100000000, maxstor=10000000, total=100
# @ xp=gaut, xlo=0, xhi=0.002
@ ylo=-100, yhi=60, yp=v, delay=50, xlo=0, xhi=100
# @ poimap=period, poivar=v, poisgn=0, poistop=0
# @ rangeover=gaut. rangestep=200, ranglelow=0, rangehigh=0.002, rangeoldic=yes
# @ rangereset=no, lt=0
# @ output=Wuoutput.dat
done