Fig 6.18 - MMSB

Model of calcium-induced calcium release in hepatocytes

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

_r618(u, p, t) = p.Rtot - u.RI - u.RIC - u.RICC
function model618!(D, u, p, t)
    hil(x, k) = x / (k + x)
    hil(x, k, n) = hil(x^n, k^n)
    @unpack k1, km1, k2, km2, k3, km3, vr, γ0, γ1, p1, p2, Cer, I, Rtot = p
    @unpack RI, RIC, RICC, C = u
    R = _r618(u, p, t)
    v1 = k1 * I * R - km1 * RI
    v2 = k2 * C * RI - km2 * RIC
    v3 = k3 * C * RIC - km3 * RICC
    v4 = vr * (γ0 + γ1 * RIC) * (Cer - C) - p1 * hil(C, p2, 4)
    D.RI = v1 - v2
    D.RIC = v2 - v3
    D.RICC = v3
    D.C = v4
    nothing
end

model618! (generic function with 1 method)

ps618 = ComponentArray(
    k1 = 12.0,
    km1 = 8.0,
    k2 = 15.0,
    km2 = 1.65,
    k3 = 1.8,
    km3 = 0.21,
    vr = 0.185,
    γ0 = 0.1,
    γ1 = 20.5,
    p1 = 8.5,
    p2 = 0.065,
    Cer = 8.37,
    I = 0.0,
    Rtot = 1.0
)

u0618 = ComponentArray(
    C = 0.0,
    RIC = 0.0,
    RICC = 0.0,
    RI = 0.0
)

ComponentVector{Float64}(C = 0.0, RIC = 0.0, RICC = 0.0, RI = 0.0)

Fig 6.18 (A)¶

ps618a = ComponentArray(ps618; I = 2.0)
tend = 25.0
prob618a = ODEProblem(model618!, u0618, tend, ps618a)
@time sol618a = solve(prob618a, TRBDF2())

  3.184190 seconds (13.84 M allocations: 917.779 MiB, 5.86% gc time, 99.92% compilation time)

retcode: Success
Interpolation: 3rd order Hermite
t: 130-element Vector{Float64}:
  0.0
  5.1031034455724376e-6
  5.613413790129681e-5
  0.0005664444824585405
  0.005669547928030976
  0.014623065166013298
  0.028865586780927013
  0.049765795640144904
  0.0801215720712506
  0.12151486555199609
  ⋮
 23.17079591207013
 23.218630914553177
 23.249313546288032
 23.34806167877875
 23.513276530671067
 23.776512335769134
 24.06164559459949
 24.541246152872603
 25.0
u: 130-element Vector{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(C = 1, RIC = 2, RICC = 3, RI = 4)}}}}:
 ComponentVector{Float64}(C = 0.0, RIC = 0.0, RICC = 0.0, RI = 0.0)
 ComponentVector{Float64}(C = 7.901900157306073e-7, RIC = 3.0682382937080804e-15, RICC = 4.612307390177343e-27, RI = 0.00012246448336776218)
 ComponentVector{Float64}(C = 8.692086072223609e-6, RIC = 3.882829629189479e-12, RICC = 7.061529274251455e-22, RI = 0.0013460101620063807)
 ComponentVector{Float64}(C = 8.77106612115255e-5, RIC = 3.944313709516432e-9, RICC = 9.145920895579543e-17, RI = 0.01347232476228388)
 ComponentVector{Float64}(C = 0.000878088016399012, RIC = 3.6784529858890115e-6, RICC = 1.0581076224824966e-11, RI = 0.12454825957643698)
 ComponentVector{Float64}(C = 0.002271104790886495, RIC = 5.1383237212217004e-5, RICC = 7.963158288960104e-10, RI = 0.280776495053278)
 ComponentVector{Float64}(C = 0.004554634626292422, RIC = 0.00033072495318715704, RICC = 1.892309788951009e-8, RI = 0.4534522094771505)
 ComponentVector{Float64}(C = 0.008312481169457603, RIC = 0.0013948381339247576, RICC = 2.412269889333395e-7, RI = 0.5989666946255119)
 ComponentVector{Float64}(C = 0.015491768820918671, RIC = 0.0047783417803678235, RICC = 2.367387481539599e-6, RI = 0.6915158390860954)
 ComponentVector{Float64}(C = 0.028041553772052436, RIC = 0.01379027202436416, RICC = 1.833409769336898e-5, RI = 0.7254375885723534)
 ⋮
 ComponentVector{Float64}(C = 0.08513135651292704, RIC = 0.17949380328498832, RICC = 0.7485456217398918, RI = 0.055553292160303795)
 ComponentVector{Float64}(C = 0.07677096815442358, RIC = 0.17517913585556258, RICC = 0.7422801031054088, RI = 0.06360420267350435)
 ComponentVector{Float64}(C = 0.07679094122343164, RIC = 0.17273980163664537, RICC = 0.7382487307697272, RI = 0.06842263172708603)
 ComponentVector{Float64}(C = 0.0747303256413803, RIC = 0.16648294751983936, RICC = 0.7253549968800836, RI = 0.08260288731612718)
 ComponentVector{Float64}(C = 0.07284200630300504, RIC = 0.16019310635658965, RICC = 0.7041355535844593, RI = 0.10297884691754)
 ComponentVector{Float64}(C = 0.0719406010433005, RIC = 0.15720620352724687, RICC = 0.6715442803553211, RI = 0.12939704617667475)
 ComponentVector{Float64}(C = 0.07259025635850241, RIC = 0.1596715057398764, RICC = 0.6382035438980028, RI = 0.15237178809131313)
 ComponentVector{Float64}(C = 0.07558132310476838, RIC = 0.17004537652992946, RICC = 0.5870531678971557, RI = 0.18275924926253315)
 ComponentVector{Float64}(C = 0.07985526458698419, RIC = 0.18333538903525337, RICC = 0.5439326625512514, RI = 0.20500187739306777)

fig = Figure()
ax = Axis(fig[1, 1], xlabel="Time", ylabel="Abundance", title="Fig 6.18 (A)\nCalcium-induced calcium release in hepatocytes")
lines!(ax, 0..tend, t -> sol618a(t).RIC, label="RIC")
lines!(ax, 0..tend, t -> sol618a(t).RICC, label="RICC")
lines!(ax, 0..tend, t -> sol618a(t).C, label="C")
axislegend(ax, position=:rt)
fig

Fig 6.18 (B)¶

affect_i1!(integrator) = integrator.p.I = 0.7
affect_i2!(integrator) = integrator.p.I = 1.2
affect_i3!(integrator) = integrator.p.I = 4.0
event_1 = PresetTimeCallback(20.0, affect_i1!)
event_2 = PresetTimeCallback(60.0, affect_i2!)
event_3 = PresetTimeCallback(90.0, affect_i3!)
cbs = CallbackSet(event_1, event_2, event_3)

tend = 120.0
prob618b = ODEProblem(model618!, u0618, tend, ps618)

ODEProblem with uType ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(C = 1, RIC = 2, RICC = 3, RI = 4)}}} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 120.0)
u0: ComponentVector{Float64}(C = 0.0, RIC = 0.0, RICC = 0.0, RI = 0.0)

@time sol618b = solve(prob618b, TRBDF2(), callback=cbs)

  0.872292 seconds (3.03 M allocations: 196.618 MiB, 5.56% gc time, 99.30% compilation time)

retcode: Success
Interpolation: 3rd order Hermite
t: 415-element Vector{Float64}:
   0.0
   9.999999999999999e-5
   0.0010999999999999998
   0.011099999999999997
   0.11109999999999996
   0.1315904332337221
   0.23404259940233282
   0.27075993885074146
   0.5062606199915098
   0.6401486705792568
   ⋮
 118.08994110676943
 118.13277318787593
 118.17560526898242
 118.29028919348603
 118.46325525235876
 118.73720855433824
 119.03450585526818
 119.53670312341583
 120.0
u: 415-element Vector{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(C = 1, RIC = 2, RICC = 3, RI = 4)}}}}:
 ComponentVector{Float64}(C = 0.0, RIC = 0.0, RICC = 0.0, RI = 0.0)
 ComponentVector{Float64}(C = 1.5484485676844983e-5, RIC = 0.0, RICC = 0.0, RI = 0.0)
 ComponentVector{Float64}(C = 0.00017032776685632227, RIC = 0.0, RICC = 0.0, RI = 0.0)
 ComponentVector{Float64}(C = 0.0017185886041296257, RIC = 0.0, RICC = 0.0, RI = 0.0)
 ComponentVector{Float64}(C = 0.01605093263445068, RIC = 0.0, RICC = 0.0, RI = 0.0)
 ComponentVector{Float64}(C = 0.018348813813633527, RIC = 0.0, RICC = 0.0, RI = 0.0)
 ComponentVector{Float64}(C = 0.023578091052589883, RIC = 0.0, RICC = 0.0, RI = 0.0)
 ComponentVector{Float64}(C = 0.023819780953078767, RIC = 0.0, RICC = 0.0, RI = 0.0)
 ComponentVector{Float64}(C = 0.024001728560150547, RIC = 0.0, RICC = 0.0, RI = 0.0)
 ComponentVector{Float64}(C = 0.023969416206032417, RIC = 0.0, RICC = 0.0, RI = 0.0)
 ⋮
 ComponentVector{Float64}(C = 0.08525407450631939, RIC = 0.18492818901449287, RICC = 0.7448311978150932, RI = 0.06076581001681161)
 ComponentVector{Float64}(C = 0.07904023424425649, RIC = 0.18092066313855648, RICC = 0.7393071535472192, RI = 0.0689413283287082)
 ComponentVector{Float64}(C = 0.07832686555459364, RIC = 0.17747709091353853, RICC = 0.7337698930027249, RI = 0.0766059951889182)
 ComponentVector{Float64}(C = 0.07598081053424906, RIC = 0.1706694253163579, RICC = 0.7190439674261723, RI = 0.0949840944093115)
 ComponentVector{Float64}(C = 0.0743088323904794, RIC = 0.16535562667951, RICC = 0.697246151000266, RI = 0.11813721300124799)
 ComponentVector{Float64}(C = 0.07399363833153798, RIC = 0.16453167754148051, RICC = 0.6641048468595826, RI = 0.1471631679424781)
 ComponentVector{Float64}(C = 0.07545202459971793, RIC = 0.16958976043490814, RICC = 0.6303717188554444, RI = 0.1716825974419474)
 ComponentVector{Float64}(C = 0.08024095361322868, RIC = 0.1844866619338284, RICC = 0.5790670291869783, RI = 0.20282393770580884)
 ComponentVector{Float64}(C = 0.08687426225151251, RIC = 0.20167551095274625, RICC = 0.538189678536161, RI = 0.22308173929642008)

fig = Figure()
ax = Axis(fig[1, 1], xlabel="Time", ylabel="Ca concentration", title="Fig 6.18 (B)\nCalcium-induced calcium release in hepatocytes")
lines!(ax, 0..tend, t -> sol618b(t).C, label="C")
axislegend(ax, position=:rt)
fig

This notebook was generated using Literate.jl.