Fig 6.3 - MMSB

Two component pathway

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

function model603!(D, u, p, t)
    @unpack k1, km1, k2, k3, L, Rtot, Ptot = p
    @unpack RL, Ps = u
    P = Ptot - Ps
    R = Rtot - RL
    v1 = k1 * L * R - km1 * RL
    v2 = k2 * RL * P
    v3 = k3 * Ps
    D.RL = v1
    D.Ps = v2 - v3
    nothing
end

model603! (generic function with 1 method)

ps603 = ComponentArray(
    k1=5.0,
    km1=1.0,
    k2=6.0,
    k3=2.0,
    L=0.0,
    Rtot=3.0,
    Ptot=8.0
)

u0603 = ComponentArray(
    RL=0.0,
    Ps=0.0
)

ComponentVector{Float64}(RL = 0.0, Ps = 0.0)

Events

affect_L1!(integrator) = integrator.p.L = 3.0
affect_L2!(integrator) = integrator.p.L = 0.0
event_L1 = PresetTimeCallback([1.0], affect_L1!)
event_L2 = PresetTimeCallback([3.0], affect_L2!)
cbs = CallbackSet(event_L1, event_L2)

SciMLBase.CallbackSet{Tuple{}, Tuple{SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##287".affect_L1!)}, typeof(Main.var"##287".affect_L1!), DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##287".affect_L1!)}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}, SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##287".affect_L2!)}, typeof(Main.var"##287".affect_L2!), DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##287".affect_L2!)}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}}}((), (SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##287".affect_L1!)}, typeof(Main.var"##287".affect_L1!), DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##287".affect_L1!)}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}(DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##287".affect_L1!)}([1.0], true, SciMLBase.INITIALIZE_DEFAULT, Main.var"##287".affect_L1!), Main.var"##287".affect_L1!, DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##287".affect_L1!)}([1.0], true, SciMLBase.INITIALIZE_DEFAULT, Main.var"##287".affect_L1!), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing, ()), SciMLBase.DiscreteCallback{DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##287".affect_L2!)}, typeof(Main.var"##287".affect_L2!), DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##287".affect_L2!)}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}(DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##287".affect_L2!)}([3.0], true, SciMLBase.INITIALIZE_DEFAULT, Main.var"##287".affect_L2!), Main.var"##287".affect_L2!, DiffEqCallbacks.PresetTimeFunction{Vector{Float64}, typeof(SciMLBase.INITIALIZE_DEFAULT), typeof(Main.var"##287".affect_L2!)}([3.0], true, SciMLBase.INITIALIZE_DEFAULT, Main.var"##287".affect_L2!), SciMLBase.FINALIZE_DEFAULT, Bool[1, 1], nothing, ())))

Fig. 6.3 A¶

tend = 10.0
prob603a = ODEProblem(model603!, u0603, tend, ps603)
@time sol603a = solve(prob603a, KenCarp47(), callback=cbs)

  7.965807 seconds (23.90 M allocations: 1.154 GiB, 1.67% gc time, 99.98% compilation time)

retcode: Success
Interpolation: 3rd order Hermite
t: 21-element Vector{Float64}:
  0.0
  9.999999999999999e-5
  0.0010999999999999998
  0.011099999999999997
  0.11109999999999996
  1.0
  1.0
  1.125
  1.2177545069860003
  1.4620499505051638
  ⋮
  3.0
  3.0
  3.573018109346427
  4.187581809384962
  5.410454977791814
  6.633328146198666
  8.231667896488167
  9.952410929564572
 10.0
u: 21-element Vector{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(RL = 1, Ps = 2)}}}}:
 ComponentVector{Float64}(RL = 0.0, Ps = 0.0)
 ComponentVector{Float64}(RL = 0.0, Ps = 0.0)
 ComponentVector{Float64}(RL = 0.0, Ps = 0.0)
 ComponentVector{Float64}(RL = 0.0, Ps = 0.0)
 ComponentVector{Float64}(RL = 0.0, Ps = 0.0)
 ComponentVector{Float64}(RL = 0.0, Ps = 0.0)
 ComponentVector{Float64}(RL = 0.0, Ps = 0.0)
 ComponentVector{Float64}(RL = 2.4294624699557787, Ps = 5.2120810276262155)
 ComponentVector{Float64}(RL = 2.7255453099084783, Ps = 6.736841418278899)
 ComponentVector{Float64}(RL = 2.810274122624125, Ps = 7.1431341864406415)
 ⋮
 ComponentVector{Float64}(RL = 2.812499963583009, Ps = 7.152312640756655)
 ComponentVector{Float64}(RL = 2.812499963583009, Ps = 7.152312640756655)
 ComponentVector{Float64}(RL = 1.585763851826017, Ps = 6.698527880240048)
 ComponentVector{Float64}(RL = 0.8577190343012714, Ps = 5.95548691376097)
 ComponentVector{Float64}(RL = 0.25262183656552534, Ps = 3.941849387941322)
 ComponentVector{Float64}(RL = 0.07440407646242693, Ps = 1.994610675493069)
 ComponentVector{Float64}(RL = 0.01507622131397586, Ps = 0.5819613613775646)
 ComponentVector{Float64}(RL = 0.002705368615998952, Ps = 0.1215200833111412)
 ComponentVector{Float64}(RL = 0.0025796380492617435, Ps = 0.1161565025353953)

fig = Figure()
ax = Axis(fig[1, 1], xlabel="Time", ylabel="Concentration", title="Fig. 6.3 (A)")
lines!(ax, 0 .. 10, t -> sol603a(t).RL, label="RL")
lines!(ax, 0 .. 10, t -> sol603a(t).Ps, label="P*")
lines!(ax, 0 .. 10, t -> 3 * (1 <= t <= 3), label="Ligand", linestyle=:dash, color=:black, alpha=0.7)
axislegend(ax)
fig

Fig 6.3 B¶

lrange = 0:0.01:1

prob_func = (prob, i, repeat) -> remake(prob, p=ComponentArray(ps603; L=lrange[i]))
prob603b = SteadyStateProblem(model603!, u0603, ps603)
trajectories = length(lrange)
alg = DynamicSS(KenCarp47())
eprob = EnsembleProblem(prob603b; prob_func)
@time sim = solve(eprob, alg; trajectories);

pstar = map(s -> s.u.Ps, sim)
rl = map(s -> s.u.RL, sim)

fig = Figure()
ax = Axis(fig[1, 1],
    xlabel="Ligand",
    ylabel="Steady-state concentration",
    title="Fig. 6.3 (B)"
)
lines!(ax, lrange, pstar, label="P*")
lines!(ax, lrange, rl, label="RL")
limits!(ax, 0, 1, 0, 8)
axislegend(ax, position=:rc)
fig

  3.940189 seconds (11.17 M allocations: 577.274 MiB, 1.39% gc time, 183.19% compilation time)

This notebook was generated using Literate.jl.