istim as time-dep param

docs/figs/ch1-09.jl

@@ -1,8 +1,6 @@
 #===
-
 # Fig 1.09
 Hodgkin-Huxley model
-
 ===#
 using DifferentialEquations
 using ModelingToolkit
@@ -16,6 +14,7 @@ exprel(x) = x / expm1(x)
 
 # HH Neuron model
 function build_hh(;name)
+    @independent_variables t
     @parameters begin
         E_N = 55.0       ## Reversal potential of Na (mV)
         E_K = -72.0      ## Reversal potential of K (mV)
@@ -24,10 +23,8 @@ function build_hh(;name)
         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))
+        iStim(t) = 0.0   ## stimulation current
     end
-
-    @independent_variables t
-
     @variables begin
         mα(t)
         mβ(t)
@@ -38,7 +35,6 @@ function build_hh(;name)
         iNa(t)
         iK(t)
         iLeak(t)
-        iStim(t)
         v(t) = -59.8977
         m(t) = 0.0536
         h(t) = 0.5925
@@ -57,27 +53,27 @@ function build_hh(;name)
         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α,
     ]
+    discrete_events = [[20] => [iStim ~ -6.6], [21] => [iStim ~ 0], [60] => [iStim ~ -6.9], [61] => [iStim ~ 0]]
 
-    return ODESystem(eqs, t; name)
+    return ODESystem(eqs, t; name, discrete_events)
 end
 
 #---
 tend = 100.0
 @mtkbuild sys = build_hh()
-prob = ODEProblem(sys, [], tend, [])
+prob = ODEProblem(sys, [], tend)
 
 #---
-sol = solve(prob, tstops=[20., 60.])
+sol = solve(prob)
 
 #---
 @unpack v, m, h, n, iStim = sys
-p1 = plot(sol, idxs=[v], ylabel="Membrane potential (mV)", xlabel="", legend=false, title="Fig 1.9")
+p1 = plot(sol, idxs = v, ylabel="Voltage (mV)", xlabel="",  labels="Membrane potential", title="Fig 1.9", legend=:topleft)
 p2 = plot(sol, idxs = [m, h, n], xlabel="")
-p3 = plot(sol, idxs = [iStim], xlabel="Time (ms)", ylabel="Current", labels="Stimulation current")
+p3 = plot(sol, idxs = iStim, xlabel="Time (ms)", ylabel="Current (uA/cm^2)", labels="Stimulation current", legend=:left)
 plot(p1, p2, p3, layout=(3, 1), size=(600, 900), leftmargin=5*Plots.mm)