# 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
# 
# 7 AUG 2017:
# introduce parameters pharmP, pharmB to block inhibition in pBC, BC, resp.
# change name to basenetTA.ode

# 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

foute(v)=1./(1+exp((v-V12outpute)/Koutpute))
fouti(v)=1./(1+exp((v-V12outputi)/Koutputi))

# 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

init v1=-60,v2=-60,v3=-60,v4=-60
init h1=0.35,h2=0.35,h3=0.35,h4=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)-pharmP*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)-pharmP*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)-pharmB*gsynI*(b23*fouti(v2)+b43*fouti(v4)+sdum3)*(v3-EsynI)-gsynE*(c13*Drive1)*(v3-EsynE)+Iapp3)/C
v4'=(-Inapi(v4,h4)-Iad(v4,m4)-ILI(v4)-pharmB*gsynI*(b24*fouti(v2)+b34*fouti(v3)+sdum4)*(v4-EsynI)-gsynE*(c14*Drive1+c24*Drive2)*(v4-EsynE)+Iapp4)/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)
m2'=(-m2+Kad2*fouti(v2))/Tad2
m3'=(-m3+Kad3*fouti(v3))/Tad3
m4'=(-m4+Kad4*fouti(v4))/Tad4

par pharmP=1.0,pharmB=1.0

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)

only t,fv1,v1

@ dt=0.1,total=20000,transient=10000,meth=qualrk,bound=5000,maxstor=10000001
@ range=1,rangereset=no,rangeover=b31,rangestep=14,rangelow=0.1,rangehigh=0.52,rangeoldic=yes
@ output=b31range_ton.test

done