Fig 6.16 - MMSB

Model of apoptosis signalling pathway

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

function model616!(D, u, p, t)
    @unpack k1, k2, k3, k4, k5, k6, k7, k8, k9, k10, k11, k12, k13, k14, k15, k16, k17, k18, k19, I = p
    @unpack C8, C8s, C3, C3s, BAR, IAP, C8sBAR, C3sIAP = u
    v1 = k1 - k2 * C8
    v2 = k3 * (C3s + I) * C8
    v3 = k4 * C8s
    v4 = k5 * C8s * BAR - k6 * C8sBAR
    v5 = k7 - k8 * C3
    v6 = k9 * C8s * C3
    v7 = k10 * C3s
    v8 = k11 * C3s * IAP - k12 * C3sIAP
    v9 = k13 - k14 * BAR
    v10 = k15 - (k16 + k17 * C3s) * IAP
    v11 = k18 * C8sBAR
    v12 = k19 * C3sIAP
    D.C8 = v1 - v2
    D.C8s = v2 - v3 - v4
    D.C3 = v5 - v6
    D.C3s = v6 - v7 - v8
    D.BAR = v9 - v4
    D.IAP = v10 - v8
    D.C8sBAR = v4 - v11
    D.C3sIAP = v8 - v12
    nothing
end

ps616 = ComponentArray(
    k1 = 507.0,
    k2 = 3.9e-3,
    k3 = 1e-5,
    k4 = 5.8e-3,
    k5 = 5e-4,
    k6 = 0.21,
    k7 = 81.9,
    k8 = 3.9e-3,
    k9 = 5.8e-6,
    k10 = 5.8e-3,
    k11 = 5e-4,
    k12 = 0.21,
    k13 = 40.0,
    k14 = 1e-3,
    k15 = 464.0,
    k16 = 1.16e-2,
    k17 = 3e-4,
    k18 = 1.16e-2,
    k19 = 1.73e-2,
    I = 0.0
)

u0616 = ComponentArray(
    C8 = 1.3E5,
    C8s = 0.0,
    C3 = 0.21E5,
    C3s = 0.0,
    BAR = 0.4E5,
    IAP = 0.4E5,
    C8sBAR = 0.0,
    C3sIAP = 0.0
)

ComponentVector{Float64}(C8 = 130000.0, C8s = 0.0, C3 = 21000.0, C3s = 0.0, BAR = 40000.0, IAP = 40000.0, C8sBAR = 0.0, C3sIAP = 0.0)

Event: increase I at t=100, decrease I at t=1200

affect_i1!(integrator) = integrator.p.I = 200.0
affect_i2!(integrator) = integrator.p.I = 0.0
event_i1 = PresetTimeCallback([100.0], affect_i1!)
event_i2 = PresetTimeCallback([1200.0], affect_i2!)
cbs = CallbackSet(event_i1, event_i2)
tend = 1800.0
prob = ODEProblem(model616!, u0616, tend, ps616)

ODEProblem with uType ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(C8 = 1, C8s = 2, C3 = 3, C3s = 4, BAR = 5, IAP = 6, C8sBAR = 7, C3sIAP = 8)}}} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 1800.0)
u0: ComponentVector{Float64}(C8 = 130000.0, C8s = 0.0, C3 = 21000.0, C3s = 0.0, BAR = 40000.0, IAP = 40000.0, C8sBAR = 0.0, C3sIAP = 0.0)

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

  3.416280 seconds (14.36 M allocations: 951.649 MiB, 6.48% gc time, 99.95% compilation time)

retcode: Success
Interpolation: 3rd order Hermite
t: 111-element Vector{Float64}:
    0.0
    9.999999999999999e-5
    0.0010999999999999998
    0.011099999999999997
    0.11109999999999996
    1.1110999999999995
   11.111099999999993
  100.0
  100.0
  100.06195510015208
    ⋮
 1241.561632628816
 1249.264191805789
 1290.5099275795083
 1303.911757651201
 1417.7834741983045
 1443.1408419372815
 1584.5171441145644
 1622.9073233724994
 1800.0
u: 111-element Vector{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(C8 = 1, C8s = 2, C3 = 3, C3s = 4, BAR = 5, IAP = 6, C8sBAR = 7, C3sIAP = 8)}}}}:
 ComponentVector{Float64}(C8 = 130000.0, C8s = 0.0, C3 = 21000.0, C3s = 0.0, BAR = 40000.0, IAP = 40000.0, C8sBAR = 0.0, C3sIAP = 0.0)
 ComponentVector{Float64}(C8 = 130000.0, C8s = 0.0, C3 = 21000.0, C3s = 0.0, BAR = 40000.0, IAP = 40000.0, C8sBAR = 0.0, C3sIAP = 0.0)
 ComponentVector{Float64}(C8 = 130000.0, C8s = 0.0, C3 = 21000.0, C3s = 0.0, BAR = 40000.0, IAP = 40000.0, C8sBAR = 0.0, C3sIAP = 0.0)
 ComponentVector{Float64}(C8 = 130000.0, C8s = 0.0, C3 = 21000.0, C3s = 0.0, BAR = 40000.0, IAP = 40000.0, C8sBAR = 0.0, C3sIAP = 0.0)
 ComponentVector{Float64}(C8 = 130000.0, C8s = 0.0, C3 = 21000.0, C3s = 0.0, BAR = 40000.0, IAP = 40000.0, C8sBAR = 0.0, C3sIAP = 0.0)
 ComponentVector{Float64}(C8 = 130000.0, C8s = 0.0, C3 = 21000.0, C3s = 0.0, BAR = 40000.0, IAP = 40000.0, C8sBAR = 0.0, C3sIAP = 0.0)
 ComponentVector{Float64}(C8 = 130000.0, C8s = 0.0, C3 = 21000.0, C3s = 0.0, BAR = 40000.0, IAP = 40000.0, C8sBAR = 0.0, C3sIAP = 0.0)
 ComponentVector{Float64}(C8 = 130000.0, C8s = 0.0, C3 = 21000.0, C3s = 0.0, BAR = 40000.0, IAP = 40000.0, C8sBAR = 0.0, C3sIAP = 0.0)
 ComponentVector{Float64}(C8 = 130000.0, C8s = 0.0, C3 = 21000.0, C3s = 0.0, BAR = 40000.0, IAP = 40000.0, C8sBAR = 0.0, C3sIAP = 0.0)
 ComponentVector{Float64}(C8 = 129983.89369370748, C8s = 9.617147330746889, C3 = 20999.960259547162, C3s = 0.02559584026504929, BAR = 39993.51096038827, IAP = 39999.977331502625, C8sBAR = 6.487210257864291, C3sIAP = 0.014138631205567447)
 ⋮
 ComponentVector{Float64}(C8 = 9228.542295342766, C8s = 74232.16035526957, C3 = 188.50587147856592, C3s = 5082.137859543312, BAR = 20.57799680348132, IAP = 268.25737372634785, C8sBAR = 3446.6455711194967, C3sIAP = 2998.947362170737)
 ComponentVector{Float64}(C8 = 9241.648230184679, C8s = 74222.79042366291, C3 = 188.5329345488147, C3s = 5080.656482071339, BAR = 20.580515516156183, IAP = 268.33538910751037, C8sBAR = 3446.630871281678, C3sIAP = 2998.945251490195)
 ComponentVector{Float64}(C8 = 9273.65907539395, C8s = 74209.90220843168, C3 = 188.57242166482806, C3s = 5073.839115259634, BAR = 20.583791454991985, IAP = 268.69479764542496, C8sBAR = 3446.5788347181797, C3sIAP = 2998.932800823995)
 ComponentVector{Float64}(C8 = 9277.464105234418, C8s = 74211.76568952105, C3 = 188.56787003690178, C3s = 5071.974677447004, BAR = 20.58322642425905, IAP = 268.79323958152617, C8sBAR = 3446.568202710099, C3sIAP = 2998.9288741950763)
 ComponentVector{Float64}(C8 = 9299.229499851712, C8s = 74226.58811520053, C3 = 188.5304024443297, C3s = 5060.697151439122, BAR = 20.578841994259438, IAP = 269.3901290433434, C8sBAR = 3446.519821054066, C3sIAP = 2998.9048716384414)
 ComponentVector{Float64}(C8 = 9302.254867060747, C8s = 74229.2807245234, C3 = 188.52355404925996, C3s = 5059.094256189493, BAR = 20.57807328542653, IAP = 269.4751841891138, C8sBAR = 3446.5158964109933, C3sIAP = 2998.9014690609456)
 ComponentVector{Float64}(C8 = 9313.255995298343, C8s = 74240.28971289778, C3 = 188.49559493196685, C3s = 5053.336167103835, BAR = 20.574963455306033, IAP = 269.78116328761047, C8sBAR = 3446.5056371168534, C3sIAP = 2998.889183314006)
 ComponentVector{Float64}(C8 = 9314.91459386362, C8s = 74242.31409224124, C3 = 188.49045218152608, C3s = 5052.4722048375825, BAR = 20.574397941842452, IAP = 269.8271348032898, C8sBAR = 3446.504845825682, C3sIAP = 2998.887350647231)
 ComponentVector{Float64}(C8 = 9319.273231255373, C8s = 74248.40502498693, C3 = 188.47497557163314, C3s = 5050.193725057417, BAR = 20.57270096263192, IAP = 269.94833378338603, C8sBAR = 3446.5031785513866, C3sIAP = 2998.882502760388)

fig = Figure()
ax = Axis(fig[1, 1], xlabel="Time", ylabel="Concentration", title="Fig 6.16\nApoptosis signalling pathway")
lines!(ax, 0..tend, t -> sol(t).C8s, label="C8s")
lines!(ax, 0..tend, t -> sol(t).C3s, label="C3s")
lines!(ax, 0..tend, t -> 100 * 200 * (100<=t<1200), label="I (×100)")
axislegend(ax, position=:rc)
fig

This notebook was generated using Literate.jl.