11# [ Hodgkin-Huxley Equation] (@id hodgkin_huxley_equation)
22
3- This tutorial shows how to programmatically construct a
3+ This tutorial shows how to construct a
44[ Catalyst] ( http://docs.sciml.ai/Catalyst/stable/ ) [ ` ReactionSystem ` ] ( @ref ) that
5- is coupled to a constraint ODE, corresponding to the [ Hodgkin–Huxley
5+ includes a coupled ODE, corresponding to the [ Hodgkin–Huxley
66model] ( https://en.wikipedia.org/wiki/Hodgkin%E2%80%93Huxley_model ) for an
77excitable cell. The Hodgkin–Huxley model is a mathematical model that describes
88how action potentials in neurons are initiated and propagated. It is a
99continuous-time dynamical system given by a coupled system of nonlinear
1010differential equations that model the electrical characteristics of excitable
1111cells such as neurons and muscle cells.
1212
13- We begin by importing some necessary packages.
13+ We'll present two different approaches for constructing the model. The first
14+ will show how it can be built entirely within a single DSL model, while the
15+ second illustrates another work flow, showing how separate models containing the
16+ chemistry and the dynamics of the transmembrane potential can be combined into a
17+ complete model.
18+
19+ We begin by importing some necessary packages:
1420``` @example hh1
15- using ModelingToolkit, Catalyst, NonlinearSolve
16- using OrdinaryDiffEq, Symbolics
17- using Plots
18- t = default_t()
19- D = default_time_deriv()
21+ using ModelingToolkit, Catalyst, NonlinearSolve, Plots, OrdinaryDiffEq
2022```
2123
22- We'll build a simple Hodgkin-Huxley model for a single neuron, with the voltage,
23- V(t), included as a constraint ODE. We first specify the transition rates for
24+ ## Building the model via the Catalyst DSL
25+ Let's build a simple Hodgkin-Huxley model for a single neuron, with the voltage,
26+ $V(t)$, included as a coupled ODE. We first specify the transition rates for
2427three gating variables, $m(t)$, $n(t)$ and $h(t)$.
2528
2629$$ s' \xleftrightarrow[\beta_s(V(t))]{\alpha_s(V(t))} s, \quad s \in \{m,n,h\} $$
4952nothing # hide
5053```
5154
52- We now declare the symbolic variable, ` V(t) ` , that will represent the
53- transmembrane potential
54-
55+ We also declare a function to represent an applied current in our model, which we
56+ will use to perturb the system and create action potentials.
5557``` @example hh1
56- @variables V(t)
57- nothing # hide
58+ Iapp(t,I₀) = I₀ * sin(2*pi*t/30)^2
5859```
5960
60- and a ` ReactionSystem ` that models the opening and closing of receptors
61+ We now declare a ` ReactionSystem ` that encompasses the Hodgkin-Huxley model.
62+ Note, we will also give the (default) values for our parameters as part of
63+ constructing the model to avoid having to specify them later on via parameter
64+ maps.
6165
6266``` @example hh1
63- hhrn = @reaction_network hhmodel begin
64- (αₙ($V), βₙ($V)), n′ <--> n
65- (αₘ($V), βₘ($V)), m′ <--> m
66- (αₕ($V), βₕ($V)), h′ <--> h
67+ hhmodel = @reaction_network hhmodel begin
68+ @parameters begin
69+ C = 1.0
70+ ḡNa = 120.0
71+ ḡK = 36.0
72+ ḡL = .3
73+ ENa = 45.0
74+ EK = -82.0
75+ EL = -59.0
76+ I₀ = 0.0
77+ end
78+
79+ @variables V(t)
80+
81+ (αₙ(V), βₙ(V)), n′ <--> n
82+ (αₘ(V), βₘ(V)), m′ <--> m
83+ (αₕ(V), βₕ(V)), h′ <--> h
84+
85+ @equations begin
86+ D(V) ~ -1/C * (ḡK*n^4*(V-EK) + ḡNa*m^3*h*(V-ENa) + ḡL*(V-EL)) + Iapp(t,I₀)
87+ end
6788end
68- nothing # hide
69- ```
70-
71- Next we create a ` ModelingToolkit.ODESystem ` to store the equation for ` dV/dt `
72-
73- ``` @example hh1
74- @parameters C=1.0 ḡNa=120.0 ḡK=36.0 ḡL=.3 ENa=45.0 EK=-82.0 EL=-59.0 I₀=0.0
75- I = I₀ * sin(2*pi*t/30)^2
76-
77- # get the gating variables to use in the equation for dV/dt
78- @unpack m,n,h = hhrn
79-
80- Dₜ = default_time_deriv()
81- eqs = [Dₜ(V) ~ -1/C * (ḡK*n^4*(V-EK) + ḡNa*m^3*h*(V-ENa) + ḡL*(V-EL)) + I/C]
82- @named voltageode = ODESystem(eqs, t)
83- nothing # hide
84- ```
85-
86- 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.
87-
88- Finally, we add this ODE into the reaction model as
89-
90- ``` @example hh1
91- @named hhmodel = extend(voltageode, hhrn)
92- hhmodel = complete(hhmodel)
93- nothing # hide
9489```
90+ For now we turn off the applied current by setting its amplitude, ` I₀ ` , to zero.
9591
9692` hhmodel ` is now a ` ReactionSystem ` that is coupled to an internal constraint
9793ODE for $dV/dt$. Let's now solve to steady-state, as we can then use these
@@ -111,14 +107,16 @@ From the artificial initial condition we specified, the solution approaches the
111107physiological steady-state via firing one action potential:
112108
113109``` @example hh1
114- plot(hhsssol, idxs = V)
110+ plot(hhsssol, idxs = hhmodel. V)
115111```
116112
117113We now save this steady-state to use as the initial condition for simulating how
118- a resting neuron responds to an applied current.
114+ a resting neuron responds to an applied current. We save the steady-state values
115+ as a mapping from the symbolic variables to their steady-states that we can
116+ later use as an initial condition:
119117
120118``` @example hh1
121- u_ss = hhsssol.u[end]
119+ u_ss = unknowns(hhmodel) .=> hhsssol(tspan[2], idxs = unknowns(hhmodel))
122120nothing # hide
123121```
124122
@@ -127,10 +125,71 @@ amplitude of the stimulus is non-zero and see if we get action potentials
127125
128126``` @example hh1
129127tspan = (0.0, 50.0)
130- @unpack I₀ = hhmodel
128+ @unpack V, I₀ = hhmodel
131129oprob = ODEProblem(hhmodel, u_ss, tspan, [I₀ => 10.0])
132130sol = solve(oprob)
133- plot(sol, vars = V, legend = :outerright)
131+ plot(sol, idxs = V, legend = :outerright)
132+ ```
133+
134+ We observe three action potentials due to the steady applied current.
135+
136+ ## Building the model via composition of separate systems for the ion channel and transmembrane voltage dynamics
137+
138+ As an illustration of how one can construct models from individual components,
139+ we now separately construct and compose the model components.
140+
141+ We start by defining systems to model each ionic current:
142+ ``` @example hh1
143+ IKmodel = @reaction_network IKmodel begin
144+ @parameters ḡK = 36.0 EK = -82.0
145+ @variables V(t) Iₖ(t)
146+ (αₙ(V), βₙ(V)), n′ <--> n
147+ @equations Iₖ ~ ḡK*n^4*(V-EK)
148+ end
149+
150+ INamodel = @reaction_network INamodel begin
151+ @parameters ḡNa = 120.0 ENa = 45.0
152+ @variables V(t) Iₙₐ(t)
153+ (αₘ(V), βₘ(V)), m′ <--> m
154+ (αₕ(V), βₕ(V)), h′ <--> h
155+ @equations Iₙₐ ~ ḡNa*m^3*h*(V-ENa)
156+ end
157+
158+ ILmodel = @reaction_network ILmodel begin
159+ @parameters ḡL = .3 EL = -59.0
160+ @variables V(t) Iₗ(t)
161+ @equations Iₗ ~ ḡL*(V-EL)
162+ end
163+ nothing # hide
164+ ```
165+
166+ We next define the voltage dynamics with unspecified values for the currents
167+ ``` @example hh1
168+ hhmodel2 = @reaction_network hhmodel2 begin
169+ @parameters C = 1.0 I₀ = 0.0
170+ @variables V(t) Iₖ(t) Iₙₐ(t) Iₗ(t)
171+ @equations D(V) ~ -1/C * (Iₖ + Iₙₐ + Iₗ) + Iapp(t,I₀)
172+ end
173+ nothing # hide
174+ ```
175+ Finally, we extend the ` hhmodel ` with the systems defining the ion channel currents
176+ ``` @example hh1
177+ @named hhmodel2 = extend(IKmodel, hhmodel2)
178+ @named hhmodel2 = extend(INamodel, hhmodel2)
179+ @named hhmodel2 = extend(ILmodel, hhmodel2)
180+ hhmodel2 = complete(hhmodel2)
181+ ```
182+ Let's again solve the system starting from the previously calculated resting
183+ state, using the same applied current as above (to verify we get the same
184+ figure). Note, we now run ` structural_simplify ` from ModelingToolkit to
185+ eliminate the algebraic equations for the ionic currents when constructing the
186+ ` ODEProblem ` :
187+
188+ ``` @example hh1
189+ @unpack I₀,V = hhmodel2
190+ oprob = ODEProblem(hhmodel2, u_ss, tspan, [I₀ => 10.0]; structural_simplify = true)
191+ sol = solve(oprob)
192+ plot(sol, idxs = V, legend = :outerright)
134193```
135194
136- We observe three action potentials due to the steady applied current .
195+ We observe the same solutions as from our original model .
0 commit comments