Hodgkin-Huxley Equation

This tutorial shows how to construct a Catalyst ReactionSystem that includes a coupled ODE, corresponding to the Hodgkin–Huxley model for an excitable cell. The Hodgkin–Huxley model is a mathematical model that describes how action potentials in neurons are initiated and propagated. It is a continuous-time dynamical system given by a coupled system of nonlinear differential equations that model the electrical characteristics of excitable cells such as neurons and muscle cells.

We'll present two different approaches for constructing the model. The first will show how it can be built entirely within a single DSL model, while the second illustrates another work flow, showing how separate models containing the chemistry and the dynamics of the transmembrane potential can be combined into a complete model.

We begin by importing some necessary packages:

using ModelingToolkit, Catalyst, NonlinearSolve, Plots, OrdinaryDiffEqRosenbrock

Building the model via the Catalyst DSL

Let's build a simple Hodgkin-Huxley model for a single neuron, with the voltage, $V(t)$, included as a coupled ODE. We first specify the transition rates for three gating variables, $m(t)$, $n(t)$ and $h(t)$.

\[s' \xleftrightarrow[\beta_s(V(t))]{\alpha_s(V(t))} s, \quad s \in \{m,n,h\}\]

Here each gating variable is used in determining the fraction of active (i.e. open) or inactive ($m' = 1 - m$, $n' = 1 -n$, $h' = 1 - h$) sodium ($m$ and $h$) and potassium ($n$) channels.

The transition rate functions, which depend on the voltage, $V(t)$, are given by

function αₘ(V)
    theta = (V + 45) / 10
    ifelse(theta == 0.0, 1.0, theta/(1 - exp(-theta)))
end
βₘ(V) = 4*exp(-(V + 70)/18)

αₕ(V) = .07 * exp(-(V + 70)/20)
βₕ(V) = 1/(1 + exp(-(V + 40)/10))

function αₙ(V)
    theta = (V + 60) / 10
    ifelse(theta == 0.0, .1, .1*theta / (1 - exp(-theta)))
end
βₙ(V) = .125 * exp(-(V + 70)/80)

We also declare a function to represent an applied current in our model, which we will use to perturb the system and create action potentials.

Iapp(t,I₀) = I₀ * sin(2*pi*t/30)^2

Iapp (generic function with 1 method)

We now declare a ReactionSystem that encompasses the Hodgkin-Huxley model. Note, we will also give the (default) values for our parameters as part of constructing the model to avoid having to specify them later on via parameter maps.

hhmodel = @reaction_network hhmodel begin
    @parameters begin
        C = 1.0
        ḡNa = 120.0
        ḡK = 36.0
        ḡL = .3
        ENa = 45.0
        EK = -82.0
        EL = -59.0
        I₀ = 0.0
    end

    @variables V(t)

    (αₙ(V), βₙ(V)), n′ <--> n
    (αₘ(V), βₘ(V)), m′ <--> m
    (αₕ(V), βₕ(V)), h′ <--> h

    @equations begin
        D(V) ~ -1/C * (ḡK*n^4*(V-EK) + ḡNa*m^3*h*(V-ENa) + ḡL*(V-EL)) + Iapp(t,I₀)
    end
end

\[ \begin{align*} \mathrm{\mathtt{n\prime}} &\xrightleftharpoons[0.125 e^{\frac{1}{80} \left( -70 - V\left( t \right) \right)}]{\mathrm{ifelse}\left( 0.1, \frac{0.010000000000000002 \left( 60 + V\left( t \right) \right)}{1 - e^{ - \frac{1}{10} \left( 60 + V\left( t \right) \right)}}; \frac{1}{10} \left( 60 + V\left( t \right) \right) = 0.0 \right)} \mathrm{n} \\ \mathrm{\mathtt{m\prime}} &\xrightleftharpoons[4 e^{\frac{1}{18} \left( -70 - V\left( t \right) \right)}]{\mathrm{ifelse}\left( 1.0, \frac{\frac{1}{10} \left( 45 + V\left( t \right) \right)}{1 - e^{ - \frac{1}{10} \left( 45 + V\left( t \right) \right)}}; \frac{1}{10} \left( 45 + V\left( t \right) \right) = 0.0 \right)} \mathrm{m} \\ \mathrm{\mathtt{h\prime}} &\xrightleftharpoons[\frac{1}{1 + e^{\frac{1}{10} \left( -40 - V\left( t \right) \right)}}]{0.07 e^{\frac{1}{20} \left( -70 - V\left( t \right) \right)}} \mathrm{h} \\ \frac{\mathrm{d} V\left( t \right)}{\mathrm{d}t} &= \frac{ - \left( - \mathtt{EL} + V\left( t \right) \right) \mathtt{{\bar{g}}L} - \left( n\left( t \right) \right)^{4} \left( - \mathtt{EK} + V\left( t \right) \right) \mathtt{{\bar{g}}K} - \left( m\left( t \right) \right)^{3} \left( - \mathtt{ENa} + V\left( t \right) \right) h\left( t \right) \mathtt{{\bar{g}}Na}}{C} + \sin^{2}\left( 0.20944 t \right) \mathtt{I{_0}} \end{align*} \]

For now we turn off the applied current by setting its amplitude, I₀, to zero.

hhmodel is now a ReactionSystem that is coupled to an internal constraint ODE for $dV/dt$. Let's now solve to steady-state, as we can then use these resting values as an initial condition before applying a current to create an action potential.

tspan = (0.0, 50.0)
u₀ = [:V => -70, :m => 0.0, :h => 0.0, :n => 0.0,
	  :m′ => 1.0, :n′ => 1.0, :h′ => 1.0]
oprob = ODEProblem(hhmodel, u₀, tspan)
hhsssol = solve(oprob, Rosenbrock23())

From the artificial initial condition we specified, the solution approaches the physiological steady-state via firing one action potential:

plot(hhsssol, idxs = hhmodel.V)

We now save this steady-state to use as the initial condition for simulating how a resting neuron responds to an applied current. We save the steady-state values as a mapping from the symbolic variables to their steady-states that we can later use as an initial condition:

u_ss = unknowns(hhmodel) .=> hhsssol(tspan[2], idxs = unknowns(hhmodel))

Finally, starting from this resting state let's solve the system when the amplitude of the stimulus is non-zero and see if we get action potentials

tspan = (0.0, 50.0)
@unpack V,I₀ = hhmodel
oprob = ODEProblem(hhmodel, u_ss, tspan, [I₀ => 10.0])
sol = solve(oprob)
plot(sol, idxs = V, legend = :outerright)

We observe three action potentials due to the steady applied current.

Building the model via composition of separate systems for the ion channel and transmembrane voltage dynamics

As an illustration of how one can construct models from individual components, we now separately construct and compose the model components.

We start by defining systems to model each ionic current. Note we now use @network_component instead of @reaction_network as we want the models to be composable and not marked as finalized.

IKmodel = @network_component IKmodel begin
    @parameters ḡK = 36.0 EK = -82.0
    @variables V(t) Iₖ(t)
    (αₙ(V), βₙ(V)), n′ <--> n
    @equations Iₖ ~ ḡK*n^4*(V-EK)
end

INamodel = @network_component INamodel begin
    @parameters ḡNa = 120.0 ENa = 45.0
    @variables V(t) Iₙₐ(t)
    (αₘ(V), βₘ(V)), m′ <--> m
    (αₕ(V), βₕ(V)), h′ <--> h
    @equations Iₙₐ ~ ḡNa*m^3*h*(V-ENa)
end

ILmodel = @network_component ILmodel begin
    @parameters ḡL = .3 EL = -59.0
    @variables V(t) Iₗ(t)
    @equations Iₗ ~ ḡL*(V-EL)
end

We next define the voltage dynamics with unspecified values for the currents

hhmodel2 = @network_component hhmodel2 begin
    @parameters C = 1.0 I₀ = 0.0
    @variables V(t) Iₖ(t) Iₙₐ(t) Iₗ(t)
    @equations D(V) ~ -1/C * (Iₖ + Iₙₐ + Iₗ) + Iapp(t,I₀)
end

Finally, we extend the hhmodel with the systems defining the ion channel currents

@named hhmodel2 = extend(IKmodel, hhmodel2)
@named hhmodel2 = extend(INamodel, hhmodel2)
@named hhmodel2 = extend(ILmodel, hhmodel2)
hhmodel2 = complete(hhmodel2)

\[ \begin{align*} \mathrm{\mathtt{n\prime}} &\xrightleftharpoons[0.125 e^{\frac{1}{80} \left( -70 - V\left( t \right) \right)}]{\mathrm{ifelse}\left( 0.1, \frac{0.010000000000000002 \left( 60 + V\left( t \right) \right)}{1 - e^{ - \frac{1}{10} \left( 60 + V\left( t \right) \right)}}; \frac{1}{10} \left( 60 + V\left( t \right) \right) = 0.0 \right)} \mathrm{n} \\ \mathrm{\mathtt{m\prime}} &\xrightleftharpoons[4 e^{\frac{1}{18} \left( -70 - V\left( t \right) \right)}]{\mathrm{ifelse}\left( 1.0, \frac{\frac{1}{10} \left( 45 + V\left( t \right) \right)}{1 - e^{ - \frac{1}{10} \left( 45 + V\left( t \right) \right)}}; \frac{1}{10} \left( 45 + V\left( t \right) \right) = 0.0 \right)} \mathrm{m} \\ \mathrm{\mathtt{h\prime}} &\xrightleftharpoons[\frac{1}{1 + e^{\frac{1}{10} \left( -40 - V\left( t \right) \right)}}]{0.07 e^{\frac{1}{20} \left( -70 - V\left( t \right) \right)}} \mathrm{h} \\ \frac{\mathrm{d} V\left( t \right)}{\mathrm{d}t} &= \frac{ - \mathtt{I_k}\left( t \right) - \mathtt{I_{n a}}\left( t \right) - \mathtt{I_l}\left( t \right)}{C} + \sin^{2}\left( 0.20944 t \right) \mathtt{I_0} \\ \mathtt{I_k}\left( t \right) &= \left( n\left( t \right) \right)^{4} \left( - \mathtt{EK} + V\left( t \right) \right) \mathtt{{\bar{g}}K} \\ \mathtt{I_{n a}}\left( t \right) &= \left( m\left( t \right) \right)^{3} \left( - \mathtt{ENa} + V\left( t \right) \right) h\left( t \right) \mathtt{{\bar{g}}Na} \\ \mathtt{I_l}\left( t \right) &= \left( - \mathtt{EL} + V\left( t \right) \right) \mathtt{{\bar{g}}L} \end{align*} \]

Let's again solve the system starting from the previously calculated resting state, using the same applied current as above (to verify we get the same figure). Note, we now run structural_simplify from ModelingToolkit to eliminate the algebraic equations for the ionic currents when constructing the ODEProblem:

@unpack I₀,V = hhmodel2
oprob = ODEProblem(hhmodel2, u_ss, tspan, [I₀ => 10.0]; structural_simplify = true)
sol = solve(oprob)
plot(sol, idxs = V, legend = :outerright)

We observe the same solutions as from our original model.