Fig 1.7 - MMSB

For Figures 1.7, 7.13, 7.14, 7.15.

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

Collins toggle switch model

function collins!(du, u, p, t)
    hil(x, k) = x / (x + k)
    hil(x, k, n) = hil(x^n, k^n)
    @unpack a1, a2, β, γ, i1, i2 = p
    @unpack s1, s2 = u
    du[1] = a1 * hil(1 + i2, s2, β) - s1
    du[2] = a2 * hil(1 + i1, s1, γ) - s2
    nothing
end

collins! (generic function with 1 method)

Setup the problem

tend = 50.0
ps = ComponentArray(a1=3.0, a2=2.5, β=4.0, γ=4.0, i1=0.0, i2=0.0)
u0 = ComponentArray(s1=0.075, s2=2.5)
prob = ODEProblem(collins!, u0, (0.0, tend), ps)

ODEProblem with uType ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(s1 = 1, s2 = 2)}}} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 50.0)
u0: ComponentVector{Float64}(s1 = 0.075, s2 = 2.5)

Callbacks and solve the problem

cbs = let
    affect_i2_on!(integrator) = integrator.p.i2 = 10.0
    affect_i2_off!(integrator) = integrator.p.i2 = 0.0
    affect_i1_on!(integrator) = integrator.p.i1 = 10.0
    affect_i1_off!(integrator) = integrator.p.i1 = 0.0
    cb_i2_on = PresetTimeCallback(10.0, affect_i2_on!)
    cb_i2_off = PresetTimeCallback(20.0, affect_i2_off!)
    cb_i1_on = PresetTimeCallback(30.0, affect_i1_on!)
    cb_i1_off = PresetTimeCallback(40.0, affect_i1_off!)
    cbs = CallbackSet(cb_i2_on, cb_i2_off, cb_i1_on, cb_i1_off)
end

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

  1.387313 seconds (5.85 M allocations: 385.570 MiB, 8.58% gc time, 99.98% compilation time)

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 45-element Vector{Float64}:
  0.0
  0.38548049109061366
  1.2786313586853497
  2.464538731594535
  4.06089298356771
  6.199065127274195
  9.163107137014421
 10.0
 10.0
 10.90492160678711
  ⋮
 37.811045850454185
 38.97117380017039
 40.0
 40.0
 41.45655048835799
 43.01214443871626
 45.107570048628844
 47.65220584313297
 50.0
u: 45-element Vector{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(s1 = 1, s2 = 2)}}}}:
 ComponentVector{Float64}(s1 = 0.075, s2 = 2.5)
 ComponentVector{Float64}(s1 = 0.07496310584167308, s2 = 2.4999747263084022)
 ComponentVector{Float64}(s1 = 0.0749189538414308, s2 = 2.499943104158718)
 ComponentVector{Float64}(s1 = 0.0748994736010588, s2 = 2.49992797574084)
 ComponentVector{Float64}(s1 = 0.07489347508722315, s2 = 2.4999227222886917)
 ComponentVector{Float64}(s1 = 0.07489238103870041, s2 = 2.4999215657222047)
 ComponentVector{Float64}(s1 = 0.07489230549277733, s2 = 2.4999214349680683)
 ComponentVector{Float64}(s1 = 0.07489223296465886, s2 = 2.49992138933308)
 ComponentVector{Float64}(s1 = 0.07489223296465886, s2 = 2.49992138933308)
 ComponentVector{Float64}(s1 = 1.813887696370784, s2 = 1.607495538140783)
 ⋮
 ComponentVector{Float64}(s1 = 0.07794218751213779, s2 = 2.498990428766632)
 ComponentVector{Float64}(s1 = 0.07588614749855141, s2 = 2.499682949330218)
 ComponentVector{Float64}(s1 = 0.07525537152132933, s2 = 2.4998865900140177)
 ComponentVector{Float64}(s1 = 0.07525537152132933, s2 = 2.4998865900140177)
 ComponentVector{Float64}(s1 = 0.07497917389522031, s2 = 2.4999126610145854)
 ComponentVector{Float64}(s1 = 0.07491123383177443, s2 = 2.499919371700559)
 ComponentVector{Float64}(s1 = 0.07489514451004413, s2 = 2.4999210438808794)
 ComponentVector{Float64}(s1 = 0.07489273938700687, s2 = 2.499921301318458)
 ComponentVector{Float64}(s1 = 0.074892273010352, s2 = 2.4999213465656758)

Visualization

ts = 0:0.1:tend
data = Array(sol(ts; idxs=[1, 2]))
fig, ax, sp = series(ts, data; labels=["s1", "s2"], axis=(title="Fig. 1.7", xlabel="Time", ylabel="Concentration"))
axislegend(ax)
fig

This notebook was generated using Literate.jl.