## Steps to generate the bifurcation curves shown in Figure 4:
## 1) Load this .ode file into XPPAUT and 
## 2) Set gtonic to 0.3 for panel A, 0.12 for panel B, 0.22 for panel C, and 0.18 for panel D
## 3) Keystrokes to run simulation and compute bifurcation curve: (i)(g)(i)(l)(f)(a)(r)(s)
## 4) Keystrokes to save the data: (f)(w), and then click on Ok 

# BRS model fast subsystem with h as bifurcation parameter
par h=-2

# Set gtonic as described above for accordingly for panels A-D
par gtonic=0.3 
#par gtonic=0.1247
#par gtonic=0.2194
#par gtonic=0.1806

# parameters
par i=0,c=21,etonic=0
par gl=2.8,el=-65
par gna=28,ena=50
par gk=11.2,ek=-85
par gnap=2.8

# gating functions
xinf(v,vt,sig)=1/(1+exp((v-vt)/sig))
taux(v,vt,sig,tau)=tau/cosh((v-vt)/(2*sig))

# persistent sodium
pinf(v)=xinf(v,-40,-6)
inap=gnap*pinf(v)*h*(v-ena)

# transient sodium 
minf(v)=xinf(v,-34,-5)
ina=gna*minf(v)^3*(1-n)*(v-ena)

# potassium
ninf(v)=xinf(v,-29,-4)
taun(v)=taux(v,-29,-4,10)
ik=gk*n^4*(v-ek)

# leak
il=gl*(v-el)

# tonic
itonic=gtonic*(v-etonic)

# differential equations
v' = (i-il-ina-ik-inap-itonic)/c
n'=(ninf(v)-n)/taun(v)

init v=-60

# XPP settings
@ total=40000,dt=.1,meth=cvode,maxstor=10000000
@ tol=1e-8,atol=1e-8
@ xlo=0,xhi=40000,ylo=-80,yhi=20

# AUTO settings
@ parmin=-100,parmax=100,autoxmin=-2,autoxmax=2,autoymin=-80,autoymax=20
@ dsmin=1e-4,dsmax=.1,nmax=500

done