Fig 1.09 - MMSB

Hodgkin-Huxley model

using OrdinaryDiffEq
using ComponentArrays: ComponentArray
using DiffEqCallbacks
using SimpleUnPack
using CairoMakie

HH Neuron model

function hh_neuron(u, p, t)
    exprel(x) = x / expm1(x)
    @unpack E_N, E_K, E_LEAK, G_N_BAR, G_K_BAR, G_LEAK, C_M, iStim = p
    @unpack v, m, h, n = u
    ma = exprel(-0.10 * (v + 35))
    mb = 4.0 * exp(-(v + 60) / 18.0)
    ha = 0.07 * exp(-(v + 60) / 20)
    hb = 1 / (exp(-(v + 30) / 10) + 1)
    na = 0.1 * exprel(-0.1 * (v + 50))
    nb = 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)
    return (;
        dv = -(iNa + iK + iLeak + iStim) / C_M,
        dm = -(ma + mb) * m + ma,
        dh = -(ha + hb) * h + ha,
        dn = -(na + nb) * n + na,
        ma, mb, ha, hb, na, nb, iNa, iK, iLeak
    )
end

hh_neuron (generic function with 1 method)

Inplace version of the HH neuron model

function hh_neuron!(du, u, p, t)
    @unpack dv, dm, dh, dn = hh_neuron(u, p, t)
    du.v = dv
    du.m = dm
    du.h = dh
    du.n = dn
    nothing
end

hh_neuron! (generic function with 1 method)

Problem definition

ps = ComponentArray(
    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))
    iStim=0.0       ## stimulation current
)

u0 = ComponentArray(
    v=-59.8977,
    m=0.0536,
    h=0.5925,
    n=0.3192,
)

tend = 100.0
prob = ODEProblem(hh_neuron!, u0, tend, ps)

ODEProblem with uType ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(v = 1, m = 2, h = 3, n = 4)}}} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 100.0)
u0: ComponentVector{Float64}(v = -59.8977, m = 0.0536, h = 0.5925, n = 0.3192)

Callbacks for external current and solve the problem

cbs = let
    affect_stim_on1!(integrator) = integrator.p.iStim = -6.6
    affect_stim_off1!(integrator) = integrator.p.iStim = 0.0
    affect_stim_on2!(integrator) = integrator.p.iStim = -6.9
    affect_stim_off2!(integrator) = integrator.p.iStim = 0.0
    cb_stim_on1 = PresetTimeCallback(20.0, affect_stim_on1!)
    cb_stim_off1 = PresetTimeCallback(21.0, affect_stim_off1!)
    cb_stim_on2 = PresetTimeCallback(60.0, affect_stim_on2!)
    cb_stim_off2 = PresetTimeCallback(61.0, affect_stim_off2!)
    cbs = CallbackSet(cb_stim_on1, cb_stim_off1, cb_stim_on2, cb_stim_off2)
end

@time sol = solve(prob, Tsit5(), callback=cbs)

  1.392225 seconds (5.85 M allocations: 384.636 MiB, 7.97% gc time, 99.98% compilation time)

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 189-element Vector{Float64}:
   0.0
   0.28649552814014345
   0.6282062598888134
   1.1292780059864103
   1.788539543713246
   2.7280332499045175
   3.8225722202572276
   4.8381191258330505
   5.700042620198738
   6.449974938460642
   ⋮
  94.07807125059809
  94.83081709374544
  95.58381904806036
  96.33693945892486
  97.09012174895203
  97.84335317992398
  98.59662808732999
  99.34992900744201
 100.0
u: 189-element Vector{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(v = 1, m = 2, h = 3, n = 4)}}}}:
 ComponentVector{Float64}(v = -59.8977, m = 0.0536, h = 0.5925, n = 0.3192)
 ComponentVector{Float64}(v = -59.896772589630615, m = 0.053584781785332665, h = 0.592500685682472, n = 0.3192027439091019)
 ComponentVector{Float64}(v = -59.896078417587496, m = 0.05358354985022893, h = 0.5925003853539793, n = 0.31920656870895786)
 ComponentVector{Float64}(v = -59.895380703866564, m = 0.05358732070491404, h = 0.592498563192231, n = 0.3192127047916002)
 ComponentVector{Float64}(v = -59.894827103417214, m = 0.05359154097918554, h = 0.5924946675523771, n = 0.3192210724000593)
 ComponentVector{Float64}(v = -59.894615228409705, m = 0.05359428457937094, h = 0.5924881300109393, n = 0.31923235771422587)
 ComponentVector{Float64}(v = -59.89505890343335, m = 0.05359755851450737, h = 0.5924816552294706, n = 0.31924307874932606)
 ComponentVector{Float64}(v = -59.896269814551104, m = 0.053625090215204584, h = 0.5924780715825779, n = 0.3192501479046132)
 ComponentVector{Float64}(v = -59.897715332885426, m = 0.053671317462470676, h = 0.5924772857413374, n = 0.3192539371264494)
 ComponentVector{Float64}(v = -59.89809941543937, m = 0.05366424957886711, h = 0.5924790001540866, n = 0.3192552276259311)
 ⋮
 ComponentVector{Float64}(v = -59.98623587619844, m = 0.05302161756509676, h = 0.5922846705578787, n = 0.31921639393388995)
 ComponentVector{Float64}(v = -59.96859305380443, m = 0.05312410654815457, h = 0.5925413551051383, n = 0.3190635018269061)
 ComponentVector{Float64}(v = -59.949158949319276, m = 0.05324337078057165, h = 0.5927208261133716, n = 0.31896734932290344)
 ComponentVector{Float64}(v = -59.930163912602, m = 0.053364270459756426, h = 0.5928277036944383, n = 0.3189218752241481)
 ComponentVector{Float64}(v = -59.91320983696415, m = 0.05347549038166708, h = 0.5928719324601499, n = 0.31891799228867124)
 ComponentVector{Float64}(v = -59.899299395534406, m = 0.05356949564317245, h = 0.5928664831171276, n = 0.3189452118913437)
 ComponentVector{Float64}(v = -59.88890577995788, m = 0.05364220381357519, h = 0.5928254089335973, n = 0.3189929498913998)
 ComponentVector{Float64}(v = -59.882066431720574, m = 0.053692471965090964, h = 0.5927623163142357, n = 0.31905149706887587)
 ComponentVector{Float64}(v = -59.878518366198804, m = 0.053701898037931545, h = 0.5926997287437964, n = 0.31910420561357977)

Visualization

fig, ax, sp = lines(0 .. tend, t -> sol(t).v, axis=(title="Fig 1.9\nHodgkin-Huxley Neuron", xlabel="Time (ms)", ylabel="Membrane potential (mV)"))

This notebook was generated using Literate.jl.