Hodgkin-Huxley with dynamic ion concentrations (Hubel and Dahlem, 2014)

 Download zip file 
Help downloading and running models
Accession:167714
The classical Hodgkin--Huxley (HH) model neglects the time-dependence of ion concentrations in spiking dynamics. The dynamics is therefore limited to a time scale of milliseconds, which is determined by the membrane capacitance multiplied by the resistance of the ion channels, and by the gating time constants. This model includes slow dynamics in an extended HH framework that simulates time-dependent ion concentrations, pumps, and buffers. Fluxes across the neuronal membrane change intra- and extracellular ion concentrations, whereby the latter can also change through contact to reservoirs in the surroundings. The dynamics on three distinct slow times scales is determined by the cell volume-to-surface-area ratio and the membrane permeability (seconds), the buffer time constants (tens of seconds), and the slower backward buffering (minutes to hours). The modulatory dynamics and the newly emerging excitable dynamics corresponds to pathological conditions observed in epileptiform burst activity, and spreading depression in migraine aura and stroke, respectively.
Reference:
1 . Hübel N, Dahlem MA (2014) Dynamics from seconds to hours in Hodgkin-Huxley model with time-dependent ion concentrations and buffer reservoirs. PLoS Comput Biol 10:e1003941 [PubMed]
Model Information (Click on a link to find other models with that property)
Model Type: Neuron or other electrically excitable cell;
Brain Region(s)/Organism:
Cell Type(s):
Channel(s): I Na,t; I K; I K,leak; Na/K pump; I Cl, leak; I Na, leak;
Gap Junctions:
Receptor(s):
Gene(s):
Transmitter(s):
Simulation Environment: XPP;
Model Concept(s): Homeostasis; Spreading depression;
Implementer(s):
Search NeuronDB for information about:  I Na,t; I K; I K,leak; Na/K pump; I Cl, leak; I Na, leak;
#  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

Loading data, please wait...