# RESPIRATORY RHYTHM MODEL
# J. Rubin, N. Shevtsova, B. Ermentrout, J. Smith, I. Rybak
# J. Neurophysiol. 101:2146-2165, 2009

# parameters/initial conditions set to give 3-phase rhythm

# edited Dec 14, 2015:  include NaP in early-I, post-I and aug-E
# 26 JAN 2017:
# changed name from rubinetal_resp_cpg_NaP_added.ode to basenet.ode
# included applied current parameters for each neuron

# 20 DEC 2017:  include excitatory population of post-I neurons (pico)
# 3 JAN 2018: v2:  parameters for synaptic output from pico differ from pbc
# 3 JAN 2018: v3:  include excitation from pre-I to pico

# general parameters
p C=20
p gLE=5,gLI=2.8
p ELE=-65,ELI=-60
p gsynE=10
p EsynE=0
p gsynI=60
p EsynI=-75
p EK=-85
p ENa=50

# parameters of NaP and Kdr currents 
p V12n=-29,kn=-4,gKdr=3,V12mp=-40,kmp=-6
p V12hp=-40,khp=6,tauhp=6000,gNaPe=5,gNaPi=2.5

# parameters of adaptive neurons
p gAD=10,Tad2=2000,Tad3=1500,Tad4=2000
p Kad2=0.9,Kad3=1.3,Kad4=0.9

# parameters of output function
p V12outpute=-45,Koutpute=-5
p V12outputi=-30,Koutputi=-4
p V12outputp=-25,Koutputp=-12

foute(v)=1./(1+exp((v-V12outpute)/Koutpute))
fouti(v)=1./(1+exp((v-V12outputi)/Koutputi))
# for pico
foutp(v)=1./(1+exp((v-V12outputp)/Koutputp))

# between neurons
p a12=0.4,b23=0.23,b24=0.3,b31=0.1,b32=0.25,b34=0.48
p b41=0.045,b42=0.4,b43=0.09
p Drive1=1,Drive2=1,Drive3=1
p c11=0,c12=0.3,c13=0.61,c14=0.38
p c21=0.07,c22=0.3,c24=0.4,c31=0.025
p nf=1
p sdum1=0,sdum2=0,sdum3=0,sdum4=0,edum=0
p Iapp1=0,Iapp2=0,Iapp3=0,Iapp4=0
p b25=0,b45=0,sdum5=0
p c15=0,c25=0,c35=0
p Iapp5=0,a53=0,a15=0

init v1=-60,v2=-60,v3=-60,v4=-60,v5=-60,
init h1=0.35,h2=0.35,h3=0.35,h4=0.35,h5=0.35,m2=0,m3=0,m4=0

# functions
ninf(v)=1./(1+exp((v-V12n)/kn))
mpinf(v)=1./(1+exp((v-V12mp)/kmp))
hpinf(v)=1./(1+exp((v-V12hp)/khp))
tauinf(v)=tauhp/cosh((v-V12hp)/(2*khp))

# currents
Inape(v,h)=gNaPe*mpinf(v)*h*(v-ENa)
Inapi(v,h)=gNaPi*mpinf(v)*h*(v-ENa)
Ikdr(v)=gKdr*ninf(v)^4*(v-EK)
Iad(v,m)=gAD*m*(v-EK)
ILE(v)=gLE*(v-ELE)
ILI(v)=gLI*(v-ELI)

# ODEs
v1'=(-Inape(v1,h1)-Ikdr(v1)-ILE(v1)-gsynI*(b31*fouti(v3)+b41*fouti(v4)+sdum1)*(v1-EsynI)-gsynE*(c11*Drive1+c21*Drive2+c31*Drive3)*(v1-EsynE)+Iapp1)/C
v2'=(-Inapi(v2,h2)-Iad(v2,m2)-ILI(v2)-gsynI*(b32*fouti(v3)+b42*fouti(v4)+sdum2)*(v2-EsynI)-gsynE*(a12*foute(v1)+c22*Drive2+c12*Drive1+edum)*(v2-EsynE)+Iapp2)/C
v3'=(-Inapi(v3,h3)-Iad(v3,m3)-ILI(v3)-gsynI*(b23*fouti(v2)+b43*fouti(v4)+sdum3)*(v3-EsynI)-gsynE*(a53*foutp(v5)+c13*Drive1)*(v3-EsynE)+Iapp3)/C
v4'=(-Inapi(v4,h4)-Iad(v4,m4)-ILI(v4)-gsynI*(b24*fouti(v2)+b34*fouti(v3)+sdum4)*(v4-EsynI)-gsynE*(c14*Drive1+c24*Drive2)*(v4-EsynE)+Iapp4)/C
v5'=(-Inape(v5,h5)-Ikdr(v5)-ILE(v5)-gsynI*(b25*fouti(v2)+b45*fouti(v4)+sdum5)*(v5-EsynI)-gsynE*(a15*foute(v1)+c15*Drive1+c25*Drive2+c35*Drive3)*(v5-EsynE)+Iapp5)/C
h1'=(hpinf(v1)-h1)/tauinf(v1)
h2'=(hpinf(v2)-h2)/tauinf(v2)
h3'=(hpinf(v3)-h3)/tauinf(v3)
h4'=(hpinf(v4)-h4)/tauinf(v4)
h5'=(hpinf(v5)-h5)/tauinf(v5)
m2'=(-m2+Kad2*fouti(v2))/Tad2
m3'=(-m3+Kad3*fouti(v3))/Tad3
m4'=(-m4+Kad4*fouti(v4))/Tad4

aux inh1=b31*fouti(v3)+b41*fouti(v4)
aux inh2=b32*fouti(v3)+b42*fouti(v4)
aux inh3=b23*fouti(v2)+b43*fouti(v4)
aux inh4=b24*fouti(v2)+b34*fouti(v3)
aux exc2=a12*foute(v1)

# new auxiliaries
aux fv1=foute(v1)
aux fv2=fouti(v2)
aux fv3=fouti(v3)
aux fv4=fouti(v4)
aux fv5=foute(v5)

@ dt=0.1,total=10000,meth=qualrk,xp=t,yp=v1,xlow=0,xhi=1000,ylo=-80,yhi=20.,bound=5000,maxstor=10000001

done