#  Code for article by N. Huebel and M. A. Dahlem: 
# "Dynamics from seconds to hours in Hodgkin–Huxley model with 
#  time-dependent ion concentrations and buffer reservoirs" 
#  PLoS Comp Biol (2014)
#  email: niklas.huebel@gmail.com
#  Oct. 7, 2014
#
#
# 
#
# To compute the fixed point continuation of Fig. 2 run xppaut 
# with this file. To reach the exact fixed point use (I)nitial 
# and (G)o first. Then use (F)ile + (A)UTO to open the AUTO 
# interface. (R)un + (S)teady state will start the forward 
# continuation. Then change the (N)umerics parameter DS from 0.2 
# to -0.2, (G)rab (+ 'Enter') the starting point of the forward 
# continuation curve, and (R)un again.
# (Remark: The Hopf line of Fig. 4 can only be obtained by changing
# this code so that also "cli" is a parameter.) 
# 

########################################
# RATE EQUATIONS FOR 4D BISTABLE MODEL #
########################################
v'     = 1000. * v_DOT
n'     = 1000. * n_DOT
ki'    = 1000. * ki_DOT
cli'   = 1000. * cli_DOT

####################################
# PHYSIOLOGICAL RESTING CONDITIONS #
####################################
init v=-67.193253
init n=0.069410823
init ki=129.25764
init cli=9.900239

##############################
# MAIN BIFURCATION PARAMETER #
##############################
par dk=0.

##########
# GATING #
##########
A_N = 0.01 * (v + 34.0) / (1.0 - exp(-0.1 * (v + 34.0))) 
B_N = 0.125 * exp(-(v + 44.0) / 80.0)
A_M = 0.1 * (v + 30.0) / (1.0 - exp(-0.1 * (v + 30.0))) 
B_M = 4.0 * exp(-(v + 55.0) / 18.0) 
m   = A_M / (A_M + B_M)
h   = 1 - 1. / (1 + exp(-6.5*(n-0.35)))

######################
# ION CONCENTRATIONS #
######################
par vol_i=2.16
par vol_e=0.72
par nai0=25.231485
par ki0=129.25764
par cli0=9.900239
par nae0=125.30555
par ke0=4.
par cle0=123.2716

nai =  nai0 + ki0 - ki - cli0 + cli
nae = (nai0 * vol_i + nae0 * vol_e - nai * vol_i) / vol_e
cle = (cli0 * vol_i + cle0 * vol_e - cli * vol_i) / vol_e
ke  = (ki0  * vol_i + ke0  * vol_e - ki  * vol_i) / vol_e + dk

#####################
# NERNST POTENTIALS #
#####################
EK  = 26.64 * log(ke /ki)
ENA = 26.64 * log(nae/nai)
ECL =-26.64 * log(cle/cli)

#####################
# TYPES OF CURRENTS #
#####################
par G_NA_L=0.0175
par G_NA_G=100.
par G_K_L=0.05
par G_K_G=40.
par G_CL_L=0.02

par MAX_PUMP=6.8
par NA_PUMP=25
par K_PUMP=5.5

I_NA_L = G_NA_L            * (v - ENA)
I_NA_G = G_NA_G * m**3 * h * (v - ENA)
I_K_L  = G_K_L             * (v - EK)
I_K_G  = G_K_G * n**4      * (v - EK)
I_CL_L = G_CL_L            * (v - ECL)
IPUMP  = MAX_PUMP / (1.0 + exp((NA_PUMP - nai)/3.)) / (1. + exp(K_PUMP - ke))

#################
# FULL CURRENTS #
#################
I_NA = I_NA_L + I_NA_G + 3. * IPUMP
I_K  = I_K_L  + I_K_G  - 2. * IPUMP

#############################
# RATE FUNCTIONS FOR SOLVER #
#############################
par C=1
par CONV=9.55589e-05
par PHI=3

v_DOT   = -1. / C * (I_NA + I_K + I_CL_L)
n_DOT   =  PHI * (A_N * (1 - n) - B_N * n)
ki_DOT  = -1. / vol_i * CONV * I_K
cli_DOT =  1. / vol_i * CONV * I_CL_L

###############
# AUXILIARIES #
###############
aux _nai    = nai
aux _nae    = nae
aux _cle    = cle
aux _ke     = ke

aux _EK	    = EK
aux _ENA    = ENA
aux _ECL    = ECL

aux _I_NA_L = I_NA_L
aux _I_NA_G = I_NA_G
aux _I_K_L  = I_K_L
aux _I_K_G  = I_K_G
aux _I_CL_L = I_CL_L
aux _IPUMP  = IPUMP

aux _I_NA   = I_NA
aux _I_K    = I_K

########################
# INTEGRATION NUMERICS #
########################
@ meth=stiff
@ dt=5e-2
@ maxstor=10000000, bounds=10000000
@ total=500
@ bell=0

################################
# AUTO CONTINUATION PARAMETERS #
################################
@ NTST=50, NMAX=900000, NPR=100000
@ DS=0.2, DSMIN=0.1, DSMAX=0.5
@ PARMIN=-100, PARMAX=100
@ AUTOXMIN=-100, AUTOXMAX=100, AUTOYMIN=-150, AUTOYMAX=50

################
# PLOT OPTIONS #
################
@ xhi=500
@ nplot=3, yp1=v, yp2=_EK, yp3=_ENA, ylo=-150, yhi=160

done