Fig 1.09

Hodgkin-Huxley model

using DifferentialEquations
using ModelingToolkit
using Plots
Plots.default(linewidth=2)

Convenience functions

hil(x, k) = x / (x + k)
hil(x, k, n) = hil(x^n, k^n)
exprel(x) = x / expm1(x)

exprel (generic function with 1 method)

HH Neuron model

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

    @independent_variables t

    @variables begin
        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α,
    ]

    return ODESystem(eqs, t; name)
end

build_hh (generic function with 1 method)

tend = 100.0
@mtkbuild sys = build_hh()
prob = ODEProblem(sys, [], tend, [])

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100.0)
u0: 4-element Vector{Float64}:
 -59.8977
   0.0536
   0.5925
   0.3192

sol = solve(prob, tstops=[20., 60.])

retcode: Success
Interpolation: 3rd order Hermite
t: 150-element Vector{Float64}:
   0.0
   0.28649552814014345
   0.6282062598888134
   1.1292780059864103
   1.788539543713246
   2.7280332499045175
   3.8225722202572276
   4.8381191258330505
   5.700042620198738
   6.449974938460642
   ⋮
  83.74271111900204
  84.50452188467617
  85.2635572190637
  86.77618556649925
  89.58562532955888
  91.88341150430689
  94.43179906313678
  97.18664452190608
 100.0
u: 150-element Vector{Vector{Float64}}:
 [-59.8977, 0.0536, 0.5925, 0.3192]
 [-59.896772589630615, 0.053584781785332665, 0.592500685682472, 0.3192027439091019]
 [-59.896078417587496, 0.05358354985022893, 0.5925003853539793, 0.31920656870895786]
 [-59.895380703866564, 0.05358732070491404, 0.592498563192231, 0.3192127047916002]
 [-59.894827103417214, 0.05359154097918554, 0.5924946675523771, 0.3192210724000593]
 [-59.894615228409705, 0.05359428457937118, 0.5924881300109393, 0.31923235771422587]
 [-59.895058903433394, 0.05359755851450915, 0.5924816552294707, 0.31924307874932606]
 [-59.89626981455131, 0.05362509021521666, 0.592478071582578, 0.31925014790461326]
 [-59.897715332885845, 0.053671317462498445, 0.5924772857413373, 0.3192539371264496]
 [-59.898099415439695, 0.053664249578892986, 0.5924790001540865, 0.3192552276259314]
 ⋮
 [-59.39721366318545, 0.05681424427265232, 0.5958733330387953, 0.31764558621382655]
 [-59.42371104071792, 0.056748808979855125, 0.5940950680259607, 0.3188466034006442]
 [-59.489583305581625, 0.05639106127893273, 0.5926212702706732, 0.31978970525716416]
 [-59.673910113241966, 0.055200202081011916, 0.5907897617474163, 0.3208380897911074]
 [-59.954080791351544, 0.05332485741598256, 0.5905841814746136, 0.32061931407505107]
 [-60.014398725623025, 0.05284322076282067, 0.5918220833295753, 0.3196484716518568]
 [-59.96450150421426, 0.05311030476612361, 0.5928561412562576, 0.3189612969650616]
 [-59.89560866267993, 0.05355728120715776, 0.593058906066063, 0.31889448927416275]
 [-59.87107444726174, 0.05374103259628208, 0.5927454895103067, 0.31914299983250066]

@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")
plot(p1, p2, p3, layout=(3, 1), size=(600, 900), leftmargin=5*Plots.mm)

---

This notebook was generated using Literate.jl.