fix

docs/figs/ch1-07.jl

@@ -7,8 +7,7 @@ For Figures 1.7, 7.13, 7.14, 7.15
 ===#
 
 using DifferentialEquations
-using LabelledArrays
-using SimpleUnPack
+using ModelingToolkit
 using Plots
 Plots.default(linewidth=2)
 
@@ -20,37 +19,43 @@ hil(x, k) = x / (x + k)
 hil(x, k, n) = hil(x^n, k^n)
 exprel(x) = x / expm1(x)
 
-# Collins switch model
-function collins!(D, u, p, t)
-    @unpack a1, a2, β, γ, i1, i2 = p
-    @unpack s1, s2 = u
-    D.s1 = a1 * hil(1 + i2, s2, β) - s1
-    D.s2 = a2 * hil(1 + i1, s1, γ) - s2
-    return nothing
+# Define the problem
+function build_collins(;name)
+    @parameters begin
+        a1=3.0
+        a2=2.5
+        β=4.0
+        γ=4.0
+    end
+
+    @variables begin
+        t
+        s1(t)=0.075
+        s2(t)=2.5
+        i1(t)
+        i2(t)
+    end
+
+    D = Differential(t)
+
+    eqs = [
+        D(s1) ~ a1 * hil(1 + i2, s2, β) - s1
+        D(s2) ~ a2 * hil(1 + i1, s1, γ) - s2
+        i2 ~ 10 * (10 < t) * (t < 20)
+        i1 ~ 10 * (30 < t) * (t < 40)
+    ]
+    sys = ODESystem(eqs; name)
+    return structural_simplify(sys)
 end
 
-# events
-on_10!(integrator) = integrator.p.i2 = 10.
-on_20!(integrator) = integrator.p.i2 = 0.
-on_30!(integrator) = integrator.p.i1 = 10.
-on_40!(integrator) = integrator.p.i1 = 0.
+#---
+@named sys = build_collins()
 
-events = CallbackSet(
-    PresetTimeCallback(10., on_10!),
-    PresetTimeCallback(20., on_20!),
-    PresetTimeCallback(30., on_30!),
-    PresetTimeCallback(40., on_40!),
-)
-
-# Prepare
-ps = LVector(a1=3.0, a2=2.5, β=4.0, γ=4.0, i1=0.0, i2=0.0)
-u0 = LVector(s1=0.075, s2=2.5)
-tend = 50.0
-
-# solve and visualize
-prob = ODEProblem(collins!, u0, tend, ps)
-sol = solve(prob, callback=events)
+# Solve the problem
+tspan = (0.0, 50.0)
+prob = ODEProblem(sys, [], tspan, [])
+sol = solve(prob, tstops=10:10:40)
 
+# Visual
 fig = plot(sol, legend=:right, xlabel = "Time", ylabel="Concentration", title="Fig 1.7")
-
 fig |> PNG

docs/figs/ch1-09.jl

@@ -6,8 +6,8 @@ Hodgkin-Huxley model
 ===#
 
 using DifferentialEquations
-using LabelledArrays
-using SimpleUnPack
+using ModelingToolkit
+using IfElse
 using Plots
 Plots.default(linewidth=2)
 
@@ -23,49 +23,71 @@ exprel(x) = x / expm1(x)
 _istim(t) = ifelse(20 <= t <= 21, -6.6, 0.0) + ifelse(60 <= t <= 61, -6.9, 0.0)
 
 # HH Neuron model
-function hh!(D, u, p, t)
-    @unpack G_N_BAR, E_N, G_K_BAR, E_K, G_LEAK, E_LEAK, C_M = p
-    @unpack v, m, h, n = u
-    mα = exprel(-0.10 * (v + 35))
-    mβ  = 4.0 * exp(-(v + 60) / 18.0)
-    hα  = 0.07 * exp(- ( v + 60) / 20)
-    hβ  = 1 / (exp(-(v+30)/10) + 1)
-    nα  = 0.1 * exprel(-0.1 * (v+50))
-    nβ  = 0.125 * exp( -(v+60) / 80)
-    iNa = G_N_BAR * (v - E_N) * (m^3) * h
-    iK  = G_K_BAR * (v - E_K) * (n^4)
-    iLeak = G_LEAK * (v - E_LEAK)
-    iStim = _istim(t)
-    D.v = -(iNa + iK + iLeak + iStim) / C_M
-    D.m = -(mα + mβ) * m + mα
-    D.h = -(hα + hβ) * h + hα
-    D.n = -(nα + nβ) * n + nα
-    return nothing
+function build_hh(;name)
+    @parameters begin
+        E_N = 55.0       ## Reversal potential of Na (mV)
+        E_K = -72.0      ## Reversal potential of K (mV)
+        E_LEAK = -49.0   ## Reversal potential of leaky channels (mV)
+        G_N_BAR = 120.0  ## Max. Na channel conductance (mS/cm^2)
+        G_K_BAR = 36.0   ## Max. K channel conductance (mS/cm^2)
+        G_LEAK = 0.30    ## Max. leak channel conductance (mS/cm^2)
+        C_M = 1.0        ## membrane capacitance (uF/cm^2))
+    end
+
+    @variables begin
+        t
+        mα(t)
+        mβ(t)
+        hα(t)
+        hβ(t)
+        nα(t)
+        nβ(t)
+        iNa(t)
+        iK(t)
+        iLeak(t)
+        iStim(t)
+        v(t) = -59.8977
+        m(t) = 0.0536
+        h(t) = 0.5925
+        n(t) = 0.3192
+    end
+
+    D = Differential(t)
+
+    eqs = [
+        mα ~ exprel(-0.10 * (v + 35)),
+        mβ ~ 4.0 * exp(-(v + 60) / 18.0),
+        hα ~ 0.07 * exp(- ( v + 60) / 20),
+        hβ ~ 1 / (exp(-(v+30)/10) + 1),
+        nα ~ 0.1 * exprel(-0.1 * (v+50)),
+        nβ ~ 0.125 * exp( -(v+60) / 80),
+        iNa ~ G_N_BAR * (v - E_N) * (m^3) * h,
+        iK  ~ G_K_BAR * (v - E_K) * (n^4),
+        iLeak ~ G_LEAK * (v - E_LEAK),
+        iStim ~ -6.6 * (20 < t) * (t < 21) + (-6.9) * (60 < t) * (t < 61),
+        D(v) ~ -(iNa + iK + iLeak + iStim) / C_M,
+        D(m) ~ -(mα + mβ) * m + mα,
+        D(h) ~ -(hα + hβ) * h + hα,
+        D(n) ~ -(nα + nβ) * n + nα,
+    ]
+
+    sys = ODESystem(eqs; name)
+    return structural_simplify(sys)
 end
 
 #---
-ps = (
-    E_N = 55.0,       ## Reversal potential of Na (mV)
-    E_K = -72.0,      ## Reversal potential of K (mV)
-    E_LEAK = -49.0,   ## Reversal potential of leaky channels (mV)
-    G_N_BAR = 120.0,  ## Max. Na channel conductance (mS/cm^2)
-    G_K_BAR = 36.0,   ## Max. K channel conductance (mS/cm^2)
-    G_LEAK = 0.30,    ## Max. leak channel conductance (mS/cm^2)
-    C_M = 1.0         ## membrane capacitance (uF/cm^2))
-)
-u0 = LVector(v=-59.8977, m=0.0536, h=0.5925, n=0.3192)
 tend = 100.0
-
-#---
-prob = ODEProblem(hh!, u0, tend, ps)
+@named sys = build_hh()
+prob = ODEProblem(sys, [], tend, [])
 
 #---
 sol = solve(prob, tstops=[20., 60.])
 
 #---
-p1 = plot(sol, idxs=[:v], ylabel="Membrane potential (mV)", xlabel="", legend=false, title="Fig 1.9")
-p2 = plot(sol, idxs = [:m, :h, :n], xlabel="")
-p3 = plot(_istim, sol.t, xlabel="Time (ms)", ylabel="Current", labels="Stimulation current")
+@unpack v, m, h, n, iStim = sys
+p1 = plot(sol, idxs=[v], ylabel="Membrane potential (mV)", xlabel="", legend=false, title="Fig 1.9")
+p2 = plot(sol, idxs = [m, h, n], xlabel="")
+p3 = plot(sol, idxs = [iStim], xlabel="Time (ms)", ylabel="Current", labels="Stimulation current")
 fig = plot(p1, p2, p3, layout=(3, 1), size=(600, 900), leftmargin=5*Plots.mm)
 
 fig |> PNG

docs/figs/ch4-01.jl

@@ -60,7 +60,7 @@ fig |> PNG
 
 # ## Fig. 4.2 B (Phase plot)
 
-fig = plot( sols[1], idxs=(:a, :b),
+fig = plot( sols[1], idxs=(1, 2),
     xlabel="[A]", ylabel="[B]",
     aspect_ratio=:equal, legend=nothing,
     title="Fig. 4.2 B (Phase plot)",
@@ -86,7 +86,7 @@ fig |> PNG
 fig = plot(title="Fig. 4.3 B (Phase plot)")
 
 for sol in sols
-    plot!(fig, sol, idxs=(:a, :b), legend=nothing)
+    plot!(fig, sol, idxs=(1, 2), legend=nothing)
 end
 
 plot!(fig, xlabel="[A]", ylabel="[B]", xlims=(0., 2.), ylims=(0., 2.), aspect_ratio=:equal)
@@ -120,7 +120,7 @@ fig = plot(title="Fig. 4.4 A (Phase plot with vector field)")
 quiver!(fig, xx, yy, quiver=∂F44, line=(:lightgrey), arrow=(:closed), aspect_ratio=:equal)
 
 for sol in sols
-    plot!(fig, sol, idxs=(:a, :b), linealpha=0.7, legend=nothing)
+    plot!(fig, sol, idxs=(1, 2), linealpha=0.7, legend=nothing)
 end
 
 plot!(fig, size=(600, 600), xlims=(rxy[1], rxy[end]), ylims=(rxy[1], rxy[end]), xlabel="[A]", ylabel="[B]")
@@ -133,7 +133,7 @@ fig = plot(title="Fig. 4.5 A (Phase plot with nullclines)")
 
 ## Phase plots
 for sol in sols
-    plot!(fig, sol, idxs=(:a, :b), linealpha=0.7, lab=nothing)
+    plot!(fig, sol, idxs=(1, 2), linealpha=0.7, lab=nothing)
 end
 
 ## nullclines

docs/figs/ch4-15.jl

@@ -83,7 +83,7 @@ plot!(fig, Float64[], Float64[], line=(:black), label="A nullcline")
 contour!(fig, 0:0.01:4, 0:0.01:4, ∂B, levels=[0], cbar=false, line=(:black, :dash))
 plot!(fig, Float64[], Float64[], line=(:black, :dash), label="B nullcline")
 for sol in sols
-    plot!(fig, sol, idxs=(:a, :b), label=nothing)
+    plot!(fig, sol, idxs=(1, 2), label=nothing)
 end
 plot!(fig, xlim=(0, 4), ylim=(0, 4), legend=:topright, size=(600, 600), xlabel="[A]", ylabel="[B]")
 
@@ -140,7 +140,7 @@ plot!(fig, Float64[], Float64[], line=(:black), label="A nullcline")
 contour!(fig, 0:0.01:4, 0:0.01:4, ∂B, levels=[0], cbar=false, line=(:black, :dash))
 plot!(fig, Float64[], Float64[], line=(:black, :dash), label="B nullcline")
 for sol in sols
-    plot!(fig, sol, idxs=(:a, :b), label=nothing)
+    plot!(fig, sol, idxs=(1, 2), label=nothing)
 end
 plot!(fig, xlim=(0, 4), ylim=(0, 4), legend=:topright, size=(600, 600), xlabel="[A]", ylabel="[B]")
 
@@ -151,7 +151,7 @@ fig |> PNG
 sol = solve(ODEProblem(model415!, LVector(a=2.0, b=1.5), 10.0, ps2))
 
 fig = plot(title="Fig 4.17")
-plot!(fig, sol, idxs=(:a, :b), label=nothing, arrow=:closed)
+plot!(fig, sol, idxs=(1, 2), label=nothing, arrow=:closed)
 quiver!(fig, xx, yy, quiver=∂F416, line=(:lightgrey), arrow=(:closed), aspect_ratio=:equal)
 contour!(fig, 1:0.01:3, 1:0.01:3, ∂A, levels=[0], cbar=false, line=(:black))
 plot!(fig, identity, 0, 0, line=(:black), label="A nullcline")