Hodgkin-Huxley single compartment with ModelingToolkit.jl

hh_mtk.jl

using DifferentialEquations, ModelingToolkit, Plots

# Channel gating functions
αn(V) = (0.01*(V + 55))/(1 - exp(-(V + 55)/10))
βn(V) = 0.125 * exp(-(V + 65)/80)
αm(V) = (0.1*(V + 40))/(1 - exp(-(V + 40)/10))
βm(V) = 4*exp(-(V + 65)/18)
αh(V) = 0.07*exp(-(V+65)/20)
βh(V) = 1/(1 + exp(-(V + 35)/10))

# Steady-state channel gating
m∞(V) = αm(V)/(αm(V) + βm(V))
h∞(V) = αh(V)/(αh(V) + βh(V))
n∞(V) = αn(V)/(αn(V) + βn(V))

@parameters t gl gna gk El Ena Ek Cm Iapp
@variables V(t) m(t) h(t) n(t)
@derivatives D'~t

# HH equations
eqs = [D(V) ~ (Iapp - gl*(V - El) - gna*m^3*h*(V - Ena) - gk*n^4*(V - Ek))/Cm,
       D(m) ~ αm(V) * (1 - m) - βm(V)*m,
       D(h) ~ αh(V) * (1 - h) - βh(V)*h,
       D(n) ~ αn(V) * (1 - n) - βn(V)*n]

sys = ODESystem(eqs, t, [V, m, h, n], [gl, gna, gk, El, Ena, Ek, Cm, Iapp])

# Callback for applying a step current
function current_step!(integrator)
    if integrator.p[8] == zero(Float64)
      integrator.p[8] += 1.0e-3
    else
      integrator.p[8] = zero(Float64)
    end
end

iapp_cb = PresetTimeCallback([100.0, 200.0], current_step!)

Vinit = -65.0
radius = 0.0025 # given in cm 
area = 4*pi*(radius)^2 # area in cm²

# Set applied current to zero initially
u0 = [V => Vinit
      m => m∞(Vinit)
      h => h∞(Vinit)
      n => n∞(Vinit)]

p = [gl  => 0.3*area
     gna => 120.0*area
     gk  => 36.0*area
     El  => -54.4
     Ena => 50.0
     Ek  => -77.0
     Cm  => 1.0e-3*area
     Iapp => zero(Float64)]

tspan = (0.0, 250.0)
prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, Tsit5(), callback=iapp_cb)
plot(sol, vars=(0,1))