# Hodgkin-Huxley Equation

This tutorial shows how to programmatically construct a Catalyst `ReactionSystem`

that is coupled to a constraint 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 begin by importing some necessary packages.

```
using ModelingToolkit, Catalyst, NonlinearSolve
using DifferentialEquations, Symbolics
using Plots
```

We'll build a simple Hodgkin-Huxley model for a single neuron, with the voltage, V(t), included as a constraint 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 now declare the symbolic variable, `V(t)`

, that will represent the transmembrane potential

`@variables t V(t)`

and a `ReactionSystem`

that models the opening and closing of receptors

```
hhrn = @reaction_network hhmodel begin
(αₙ($V), βₙ($V)), n′ <--> n
(αₘ($V), βₘ($V)), m′ <--> m
(αₕ($V), βₕ($V)), h′ <--> h
end
```

Next we create a `ModelingToolkit.ODESystem`

to store the equation for `dV/dt`

```
@parameters C=1.0 ḡNa=120.0 ḡK=36.0 ḡL=.3 ENa=45.0 EK=-82.0 EL=-59.0 I₀=0.0
I = I₀ * sin(2*pi*t/30)^2
# get the gating variables to use in the equation for dV/dt
@unpack m,n,h = hhrn
Dₜ = Differential(t)
eqs = [Dₜ(V) ~ -1/C * (ḡK*n^4*(V-EK) + ḡNa*m^3*h*(V-ENa) + ḡL*(V-EL)) + I/C]
@named voltageode = ODESystem(eqs, t)
```

Notice, we included an applied current, `I`

, that we will use to perturb the system and create action potentials. For now we turn this off by setting its amplitude, `I₀`

, to zero.

Finally, we add this ODE into the reaction model as

`@named hhmodel = extend(voltageode, hhrn)`

`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())
nothing # hide
```

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

`plot(hhsssol, idxs = V)`

We now save this steady-state to use as the initial condition for simulating how a resting neuron responds to an applied current.

```
u_ss = hhsssol.u[end]
nothing # hide
```

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 I₀ = hhmodel
oprob = ODEProblem(hhmodel, u_ss, tspan, [I₀ => 10.0])
sol = solve(oprob)
plot(sol, vars = V, legend = :outerright)
```

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