diff --git a/docs/Project.toml b/docs/Project.toml index 078759b32c..a009c89b8d 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/src/inverse_problems/optimization_ode_param_fitting.md b/docs/src/inverse_problems/optimization_ode_param_fitting.md index 08cb8adf05..f0226a5a14 100644 --- a/docs/src/inverse_problems/optimization_ode_param_fitting.md +++ b/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) diff --git a/docs/src/model_creation/examples/hodgkin_huxley_equation.md b/docs/src/model_creation/examples/hodgkin_huxley_equation.md index a2091035b1..fd6313647a 100644 --- a/docs/src/model_creation/examples/hodgkin_huxley_equation.md +++ b/docs/src/model_creation/examples/hodgkin_huxley_equation.md @@ -1,8 +1,8 @@ # [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 @@ -10,17 +10,20 @@ 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. \ No newline at end of file +We observe the same solutions as from our original model. \ No newline at end of file