Merge pull request #1082 from isaacsas/update_hodgkin_huxley_ex

Update hodgkin huxley example

docs/Project.toml

@@ -35,7 +35,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
-StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
+#StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
@@ -73,5 +73,5 @@ SpecialFunctions = "2.4"
 StaticArrays = "1.9"
 SteadyStateDiffEq = "2.2"
 StochasticDiffEq = "6.65"
-StructuralIdentifiability = "0.5.8"
+#StructuralIdentifiability = "0.5.8"
 Symbolics = "5.30.1, 6"

docs/src/inverse_problems/optimization_ode_param_fitting.md

@@ -32,6 +32,9 @@ data_vals = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end]
 using Plots
 plot(true_sol; idxs = :P, label = "True solution", lw = 8)
 plot!(data_ts, data_vals; label = "Measurements", seriestype=:scatter, ms = 6, color = :blue)
+plt = plot(true_sol; idxs = :P, label = "True solution", lw = 8) # hide
+plot!(plt, data_ts, data_vals; label = "Measurements", seriestype=:scatter, ms = 6, color = :blue) # hide
+Catalyst.PNG(plot(plt; fmt = :png, dpi = 200)) # hide
 ```
 
 Next, we will use DiffEqParamEstim to build a loss function to measure how well our model's solutions fit the data.
@@ -75,6 +78,8 @@ We can now simulate our model for the corresponding parameter set, checking that
 oprob_fitted = remake(oprob; p = optsol.u)
 fitted_sol = solve(oprob_fitted, Tsit5())
 plot!(fitted_sol; idxs = :P, label = "Fitted solution", linestyle = :dash, lw = 6, color = :lightblue)
+plot!(plt, fitted_sol; idxs = :P, label = "Fitted solution", linestyle = :dash, lw = 6, color = :lightblue) # hide
+Catalyst.PNG(plot(plt; fmt = :png, dpi = 200)) # hide
 ```
 
 !!! note
@@ -96,6 +101,10 @@ data_vals_P = (0.8 .+ 0.4*rand(10)) .* data_sol[:P][2:end]
 plot(true_sol; idxs=[:S, :P], label=["True S" "True P"], lw=8)
 plot!(data_ts, data_vals_S; label="Measured S", seriestype=:scatter, ms=6, color=:blue)
 plot!(data_ts, data_vals_P; label="Measured P", seriestype=:scatter, ms=6, color=:red)
+plt2 = plot(true_sol; idxs=[:S, :P], label=["True S" "True P"], lw=8)  # hide
+plot!(plt2, data_ts, data_vals_S; label="Measured S", seriestype=:scatter, ms=6, color=:blue)  # hide
+plot!(plt2, data_ts, data_vals_P; label="Measured P", seriestype=:scatter, ms=6, color=:red)  # hide
+Catalyst.PNG(plot(plt2; fmt = :png, dpi = 200)) # hide
 ```
 
 In this case we would have to use the `L2Loss(data_ts, hcat(data_vals_S, data_vals_P))` and `save_idxs=[1, 4]` arguments in `loss_function`:
@@ -112,6 +121,8 @@ optsol_S_P = solve(optprob_S_P, Optim.NelderMead())
 oprob_fitted_S_P = remake(oprob; p = optsol_S_P.u)
 fitted_sol_S_P = solve(oprob_fitted_S_P)
 plot!(fitted_sol_S_P; idxs=[:S, :P], label="Fitted solution", linestyle = :dash, lw = 6, color = [:lightblue :pink])
+plot!(plt2, fitted_sol_S_P; idxs=[:S, :P], label="Fitted solution", linestyle = :dash, lw = 6, color = [:lightblue :pink]) # hide
+Catalyst.PNG(plot(plt2; fmt = :png, dpi = 200)) # hide
 ```
 
 ## [Setting parameter constraints and boundaries](@id optimization_parameter_fitting_constraints)

docs/src/model_creation/examples/hodgkin_huxley_equation.md

@@ -1,26 +1,29 @@
 # [Hodgkin-Huxley Equation](@id hodgkin_huxley_equation)
 
-This tutorial shows how to programmatically construct a
+This tutorial shows how to construct a
 [Catalyst](http://docs.sciml.ai/Catalyst/stable/) [`ReactionSystem`](@ref) that
-is coupled to a constraint ODE, corresponding to the [Hodgkin–Huxley
+includes a coupled ODE, corresponding to the [Hodgkin–Huxley
 model](https://en.wikipedia.org/wiki/Hodgkin%E2%80%93Huxley_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.
+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:
 ```@example hh1
-using ModelingToolkit, Catalyst, NonlinearSolve
-using OrdinaryDiffEq, Symbolics
-using Plots
-t = default_t()
-D = default_time_deriv()
+using ModelingToolkit, Catalyst, NonlinearSolve, Plots, OrdinaryDiffEq
 ```
 
-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
+## 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\}$$
@@ -49,49 +52,42 @@ end
 nothing # hide
 ```
 
-We now declare the symbolic variable, `V(t)`, that will represent the
-transmembrane potential
-
+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. 
 ```@example hh1
-@variables V(t)
-nothing # hide
+Iapp(t,I₀) = I₀ * sin(2*pi*t/30)^2
 ```
 
-and a `ReactionSystem` that models the opening and closing of receptors
+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.
 
 ```@example hh1
-hhrn = @reaction_network hhmodel begin
-    (αₙ($V), βₙ($V)), n′ <--> n
-    (αₘ($V), βₘ($V)), m′ <--> m
-    (αₕ($V), βₕ($V)), h′ <--> h
+hhmodel = @reaction_network hhmodel begin
+    @parameters begin
+        C = 1.0 
+        ḡNa = 120.0 
+        ḡK = 36.0 
+        ḡ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 * (ḡK*n^4*(V-EK) + ḡNa*m^3*h*(V-ENa) + ḡL*(V-EL)) + Iapp(t,I₀)
+    end
 end
-nothing # hide
-```
-
-Next we create a `ModelingToolkit.ODESystem` to store the equation for `dV/dt`
-
-```@example hh1
-@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ₜ = default_time_deriv()
-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)
-nothing # hide
-```
-
-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
-
-```@example hh1
-@named hhmodel = extend(voltageode, hhrn)
-hhmodel = complete(hhmodel)
-nothing # hide
 ```
+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
@@ -111,14 +107,16 @@ From the artificial initial condition we specified, the solution approaches the
 physiological steady-state via firing one action potential:
 
 ```@example hh1
-plot(hhsssol, idxs = V)
+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.
+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:
 
 ```@example hh1
-u_ss = hhsssol.u[end]
+u_ss = unknowns(hhmodel) .=> hhsssol(tspan[2], idxs = unknowns(hhmodel))
 nothing # hide
 ```
 
@@ -127,10 +125,71 @@ amplitude of the stimulus is non-zero and see if we get action potentials
 
 ```@example hh1
 tspan = (0.0, 50.0)
-@unpack I₀ = hhmodel
+@unpack V,I₀ = hhmodel
 oprob = ODEProblem(hhmodel, u_ss, tspan, [I₀ => 10.0])
 sol = solve(oprob)
-plot(sol, vars = V, legend = :outerright)
+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:
+```@example hh1
+IKmodel = @reaction_network IKmodel begin
+    @parameters ḡK = 36.0 EK = -82.0 
+    @variables V(t) Iₖ(t)
+    (αₙ(V), βₙ(V)), n′ <--> n
+    @equations Iₖ ~ ḡK*n^4*(V-EK)
+end
+
+INamodel = @reaction_network INamodel begin
+    @parameters ḡNa = 120.0 ENa = 45.0 
+    @variables V(t) Iₙₐ(t)
+    (αₘ(V), βₘ(V)), m′ <--> m
+    (αₕ(V), βₕ(V)), h′ <--> h
+    @equations Iₙₐ ~ ḡNa*m^3*h*(V-ENa) 
+end
+
+ILmodel = @reaction_network ILmodel begin
+    @parameters ḡL = .3 EL = -59.0 
+    @variables V(t) Iₗ(t)
+    @equations Iₗ ~ ḡL*(V-EL)
+end
+nothing # hide
+```
+
+We next define the voltage dynamics with unspecified values for the currents
+```@example hh1
+hhmodel2 = @reaction_network 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
+nothing # hide
+```
+Finally, we extend the `hhmodel` with the systems defining the ion channel currents
+```@example hh1
+@named hhmodel2 = extend(IKmodel, hhmodel2)
+@named hhmodel2 = extend(INamodel, hhmodel2)
+@named hhmodel2 = extend(ILmodel, hhmodel2)
+hhmodel2 = complete(hhmodel2)
+```
+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`:
+
+```@example hh1
+@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 three action potentials due to the steady applied current.
+We observe the same solutions as from our original model.