Chapter 1 - MMSB

Fig 1.7¶

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

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 System(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, FBDF(), tstops=[10.0, 20.0, 30.0, 40.0])

Build system: 62.225787 seconds (64.65 M allocations: 3.654 GiB, 2.81% gc time, 99.96% compilation time: 76% of which was recompilation)
Build problem: 25.358029 seconds (35.27 M allocations: 1.962 GiB, 2.75% gc time, 99.89% compilation time: 75% of which was recompilation)
Solve problem: 4.558831 seconds (6.76 M allocations: 379.062 MiB, 3.26% gc time, 99.89% compilation time)

retcode: Success
Interpolation: specialized backward-difference stiffness-aware interpolation
t: 15-element Vector{Float64}:
  0.0
  0.008511615915238276
  0.05979154576106624
  0.18909483334920232
  0.40827865077437575
  0.9585750295617196
  1.9414394399940162
  3.2076912766067505
  4.984858331839563
  7.923266680509096
 10.0
 20.0
 30.0
 40.0
 50.0
u: 15-element Vector{Vector{Float64}}:
 [2.5, 0.075]
 [2.4999993324564906, 0.07499901316443638]
 [2.4999955080415166, 0.07499337957725291]
 [2.4999869748705517, 0.074980914995141]
 [2.499975123564857, 0.07496383355214947]
 [2.4999518720065805, 0.07493083282369019]
 [2.499930098465337, 0.07490126037714417]
 [2.499920181734469, 0.07488901753396895]
 [2.4999178241877056, 0.07488707862901119]
 [2.499918909501426, 0.074889202649177]
 [2.4999205286833726, 0.07489129549531684]
 [2.4999213039219836, 0.0748923150738949]
 [2.499921381405481, 0.07489224248014598]
 [2.499921361592247, 0.07489218736862312]
 [2.4999213548157253, 0.07489218162120131]

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 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, FBDF())

Build system: 4.077274 seconds (6.34 M allocations: 350.496 MiB, 1.40% gc time, 99.17% compilation time: 20% of which was recompilation)
Build problem: 8.862363 seconds (10.07 M allocations: 548.516 MiB, 1.47% gc time, 99.54% compilation time: <1% of which was recompilation)
Solve problem: 12.087476 seconds (12.14 M allocations: 660.458 MiB, 1.63% gc time, 99.90% compilation time)

retcode: Success
Interpolation: specialized backward-difference stiffness-aware interpolation
t: 182-element Vector{Float64}:
   0.0
   0.0019301418957088378
   0.027623969322989078
   0.12760563302440844
   0.3746355725449328
   1.0201014567145164
   4.075844659345648
  13.100230938043945
  20.0
  20.0
   ⋮
  86.84190471896495
  87.9409066362059
  89.4465845030154
  91.31385483648587
  93.18112516995635
  95.04839550342682
  96.91566583689729
  99.23268816641837
 100.0
u: 182-element Vector{Vector{Float64}}:
 [-59.8977, 0.3192, 0.0536, 0.5925]
 [-59.89769165826955, 0.31920001623366245, 0.05359979475069142, 0.5925000086904005]
 [-59.89758711629514, 0.3192002388453548, 0.05359739305800681, 0.5925001129927862]
 [-59.89725054056525, 0.3192011824980739, 0.05359144076508554, 0.5925003770933467]
 [-59.896649208209425, 0.3192038129647644, 0.05358615729298744, 0.5925004163496281]
 [-59.8957403533129, 0.3192114350319489, 0.053586629477550676, 0.5924982647556389]
 [-59.89564301215512, 0.31923509915867776, 0.05358735975870319, 0.5924898618575208]
 [-59.89721595216636, 0.31924636918090815, 0.0535777453043398, 0.5925062165007702]
 [-59.89763806618148, 0.31924655525927653, 0.053574927303614905, 0.592519749649679]
 [-59.89763806618148, 0.31924655525927653, 0.053574927303614905, 0.592519749649679]
 ⋮
 [-59.96007274531345, 0.319141773266003, 0.05315989709892987, 0.5909604809560154]
 [-59.941896421356205, 0.31901473935402114, 0.05326761329957201, 0.5913717835729757]
 [-59.91487367830194, 0.31895960816856417, 0.053439854796262475, 0.5917315798769647]
 [-59.89164533369352, 0.31900597610734055, 0.053596211225019355, 0.5919474137827075]
 [-59.881096535395216, 0.3191105003933809, 0.05367433275797232, 0.5920113184568345]
 [-59.88112879690023, 0.3192112798389982, 0.05368248709829526, 0.5920230409358463]
 [-59.88707146851936, 0.3192753531079135, 0.05364852461137692, 0.5920446966491143]
 [-59.89706543028938, 0.319297441833161, 0.05358420376823441, 0.5921115642080232]
 [-59.89959333840461, 0.319290600052596, 0.053566845845859606, 0.5921480521923568]

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

This notebook was generated using Literate.jl.