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
#
#
# README
#
# This code computes time series like in Fig. 5. SD is 
# triggered either by current stimulation (dynamics of variable
# "stim") or by pump interruption (dynamics of pump coefficient
# "z"). Two regulation schemes for potassium are implemented:
# glial buffering and diffusive coupling to the vascular system
# (with a switch parameter "s", see below). For diffusive 
# coupling with elevated "k_bath" values the time series of Fig. 7
# can be simulated.
#

##################################
# RATE EQUATIONS FOR 5D SD MODEL #
##################################
v'   = 1000. * v_DOT
n'   = 1000. * n_DOT
ki'  = 1000. * ki_DOT
cli' = 1000. * cli_DOT
dk'  = 1000. * dk_DOT

#########################################################
# STIMULATE MEMBRANE WITH A 0.5sec DEPOLARIZING CURRENT #
#########################################################
stim' = 0
global 1 1.000000000001*t -  20.0        {stim=Amp}
global 1 1.000000000001*t - (20.0 + 0.5) {stim=0}
par Amp=0.
#par Amp=100.

#####################################
# SWITCH PUMP OFF FOR delta SECONDS #
#####################################
z' = 0
global 1 1.000000000001*t -  20.0          {z=0.0}
global 1 1.000000000001*t - (20.0 + delta) {z=1.0}
#par delta=0.
par delta=9.5


######################
# INITIAL CONDITIONS #
######################
init v=-67.193253
init n=0.069410823
init ki=129.25764
init cli=9.900239
init dk=0

init stim=0
init z=1

#########################
# REGULATION PARAMETERS #
#########################
# s=0: glial buffering
# s=1: diffusive coupling
# k_bath is the bath concentration
par k_bath=4
par s=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 ki0=129.25764
par nai0=25.231485
par cli0=9.900239
par ke0=4.
par nae0=125.30555
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)

############
# 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 NA_PUMP=25
par K_PUMP=5.5
par MAX_PUMP=6.8

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)) * z

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

########################
# POTASSIUM REGULATION #
########################
par lambda=3.0e-05
par b0=500
par K_cr=15.
par k1=5.e-8

J_DIFF   = lambda * (k_bath - ke)
dk_DOT_1 = J_DIFF
k2       = k1 / (1 + exp((K_cr - ke) / 1.09))
dk_DOT_0 = -k1 * dk - k2 * ke * (B0 + dk)

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

v_DOT    = -1. / C * (I_CL_L + I_NA + I_K - stim)
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
dk_DOT   = (1 - s) * dk_DOT_0 + s * dk_DOT_1

###############
# 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

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

done

Loading data, please wait...