Chapter 1 - MMSB

Fig 1.7¶

Collins toggle switch model for Figures 1.7, 7.13, 7.14, 7.15.

using Startup

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using Plots
Plots.gr(linewidth=1.5)

Plots.GRBackend()

function collins_sys(; name=:collins)
    @parameters a1=3 a2=2.5 β=4 γ=4
    @variables i1(t) i2(t) s1(t)=0.075 s2(t)=2.5
    hil(x, k) = x / (x + k)
    hil(x, k, n) = hil(x^n, k^n)
    eqs = [
        D(s1) ~ a1 * hil(1 + i2, s2, β) - s1,
        D(s2) ~ a2 * hil(1 + i1, s1, γ) - s2,
        i1 ~ 10 * (t > 30) * (t < 40),
        i2 ~ 10 * (t > 10) * (t < 20)
    ]
    return ODESystem(eqs, t; name)
end

collins_sys (generic function with 1 method)

tend = 50.0
@time "Build system" @mtkbuild sys = collins_sys()
@time "Build problem" prob = ODEProblem(sys, [], tend)
@time "Solve problem" sol = solve(prob, Tsit5(), tstops=[10.0, 20.0, 30.0, 40.0])

Build system: 5.296717 seconds (2.18 M allocations: 155.218 MiB, 0.83% gc time, 99.71% compilation time: 2% of which was recompilation)
Build problem: 4.320197 seconds (7.87 M allocations: 570.874 MiB, 3.33% gc time, 99.24% compilation time: 48% of which was recompilation)
Solve problem: 1.942443 seconds (3.44 M allocations: 238.634 MiB, 3.72% gc time, 99.84% compilation time)

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 51-element Vector{Float64}:
  0.0
  0.3854804910906137
  1.27863135868535
  2.4645387315945353
  4.060892983567712
  6.1990651272741975
  9.163107137014425
 10.0
 10.237564173675388
 10.34590562274151
  ⋮
 36.91499716488063
 37.930797103851646
 39.07484274668416
 40.0
 41.178943649884275
 42.81986325561439
 44.70418578725867
 47.24334762487948
 50.0
u: 51-element Vector{Vector{Float64}}:
 [2.5, 0.075]
 [2.4999747263084027, 0.07496310584167311]
 [2.499943104158718, 0.07491895384143081]
 [2.499927975740838, 0.07489947360105878]
 [2.4999227222886917, 0.07489347508722309]
 [2.499921565722204, 0.07489238103870029]
 [2.4999214349680665, 0.0748923054927769]
 [2.4999213893330796, 0.07489223296465887]
 [2.372119211508157, 0.639134614676165]
 [2.3191402463267385, 0.8809050681389696]
 ⋮
 [2.4974663084645137, 0.08211174811522608]
 [2.499081896309855, 0.0776109365193797]
 [2.499707057292952, 0.07579211279533471]
 [2.499713078290432, 0.07525607577666821]
 [2.499856577005165, 0.0750131183637864]
 [2.499908252504363, 0.07491847827207444]
 [2.499919138177413, 0.07489687124757002]
 [2.499920936863924, 0.0748929461571923]
 [2.499921249118503, 0.07489232760866728]

plot(sol, title="Fig. 1.7", xlabel="Time", ylabel="Concentration")

plot(sol, idxs=[sys.i1, sys.i2], labels=["i1" "i2"], xlabel="Time", ylabel="Signal strength")

Fig 1.09¶

Hodgkin-Huxley model

using Startup
using OrdinaryDiffEq
using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using Plots
Plots.gr(linewidth=1.5)

Plots.GRBackend()

function hh_sys(; name=:hh)
    exprel(x) = x / expm1(x)
    @discretes iStim(t)=0.0
    @parameters E_N=55 E_K=-72 E_LEAK=-49 G_N_BAR=120 G_K_BAR=36 G_LEAK=0.3 C_M=1
    @variables v(t)=-59.8977 m(t)=0.0536 h(t)=0.5925 n(t)=0.3192 iNa(t) iK(t) iLeak(t) ma(t) mb(t) ha(t) hb(t) na(t) nb(t)

    # Electrical stimulation events
    stim_on_1 = ModelingToolkit.SymbolicDiscreteCallback([20.0] => [iStim ~ -6.6], discrete_parameters = iStim, iv = t)
    stim_off_1 = ModelingToolkit.SymbolicDiscreteCallback([21.0] => [iStim ~ 0.0], discrete_parameters = iStim, iv = t)
    stim_on_2 = ModelingToolkit.SymbolicDiscreteCallback([60.0] => [iStim ~ -6.9], discrete_parameters = iStim, iv = t)
    stim_off_2 = ModelingToolkit.SymbolicDiscreteCallback([61.0] => [iStim ~ 0.0], discrete_parameters = iStim, iv = t)

    eqs = [
        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),
        D(v) ~ -(iNa + iK + iLeak + iStim) / C_M,
        D(m) ~ -(ma + mb) * m + ma,
        D(h) ~ -(ha + hb) * h + ha,
        D(n) ~ -(na + nb) * n + na
    ]
    sys = ODESystem(eqs, t; name, discrete_events=[stim_on_1, stim_off_1, stim_on_2, stim_off_2])
    return sys
end

hh_sys (generic function with 1 method)

tend = 100.0
@time "Build system" @mtkcompile sys = hh_sys()
@time "Build problem" prob = ODEProblem(sys, [], tend)
@time "Solve problem" sol = solve(prob, TRBDF2())

Build system: 2.102366 seconds (2.92 M allocations: 207.594 MiB, 2.70% gc time, 96.83% compilation time: 5% of which was recompilation)
Build problem: 3.699117 seconds (4.23 M allocations: 300.862 MiB, 1.55% gc time, 98.45% compilation time)
Solve problem: 10.037315 seconds (14.47 M allocations: 1007.853 MiB, 10.62% gc time, 99.90% compilation time)

retcode: Success
Interpolation: 3rd order Hermite
t: 86-element Vector{Float64}:
   0.0
   0.0439333801989886
   0.8603019183785635
   1.4238266341729644
   4.248524319244319
   7.256696580279309
  14.008668650102909
  20.0
  20.0
  20.523672537957506
   ⋮
  80.98994690756788
  82.75533576957824
  85.28523794002348
  87.0567951172999
  89.58671291457803
  92.11663071185616
  95.36073853933443
  98.60484636681271
 100.0
u: 86-element Vector{Vector{Float64}}:
 [0.3192, 0.5925, 0.0536, -59.8977]
 [0.31920037837610854, 0.5925001823589415, 0.05359578382735057, -59.89751905826826]
 [0.3192092568010022, 0.5924998590652717, 0.05358184355412253, -59.89566008856614]
 [0.3192163677409184, 0.5924970879076324, 0.05358943258416223, -59.8950814157433]
 [0.3192459341881293, 0.5924809146952182, 0.053592196231923135, -59.895028309864934]
 [0.31925609549509726, 0.5924817912524404, 0.0535783077758142, -59.89726393088571]
 [0.319243765883664, 0.5925213611499963, 0.053569581372286346, -59.8984121547399]
 [0.31924371445418437, 0.5925326003276697, 0.05357583828365717, -59.89745801667102]
 [0.31924371445418437, 0.5925326003276697, 0.05357583828365717, -59.89745801667102]
 [0.3215591731712101, 0.5890432009452046, 0.06551860167610046, -56.81727065722306]
 ⋮
 [0.311392477904874, 0.6031106305241054, 0.052022691780803765, -60.036248508020066]
 [0.3142669515137393, 0.6000691383062855, 0.0558651230148874, -59.487387329833965]
 [0.3190176931317738, 0.5935667912021313, 0.057101906306891825, -59.36939888583164]
 [0.3208963029340365, 0.5905559702437168, 0.05581313670700942, -59.58538988535852]
 [0.32106895880988423, 0.5898014922646461, 0.05365651750240193, -59.90978213551699]
 [0.3199015194882542, 0.5912671247973528, 0.05274804147431585, -60.032332929045424]
 [0.3188688381870651, 0.5928587296366981, 0.05307221874116379, -59.96919401789336]
 [0.3188802294577743, 0.5930203127977015, 0.05364577194975229, -59.88204595309853]
 [0.31904457349408305, 0.592822868458166, 0.05375258860521209, -59.86802351279653]

plot(sol, idxs=[sys.v], xlabel="Time (ms)", ylabel="Membrane potential (mV)", title="Fig 1.9")

This notebook was generated using Literate.jl.