Fig 3.03 - MMSB

Michaelis-Menten kinetics

using OrdinaryDiffEq
using ComponentArrays: ComponentArray
using SimpleUnPack
using CairoMakie

Enzyme kinetics full model

function model303(u, p, t)
    @unpack k1, km1, k2 = p
    @unpack S, ES, P = u
    E = p.ET - ES
    v1 = k1 * S * E - km1 * ES
    v2 = k2 * ES
    return (; dS=-v1, dES=v1-v2, dP=v2, E)
end

function model303!(du, u, p, t)
    @unpack dS, dES, dP = model303(u, p, t)
    du.S = dS
    du.ES = dES
    du.P = dP
    nothing
end

model303! (generic function with 1 method)

ps303 = ComponentArray(k1=30.0, km1=1.0, k2=10.0, ET=1.0)
u0 = ComponentArray(S=5.0, ES=0.0, P=0.0)
tend = 1.0
prob303 = ODEProblem(model303!, u0, tend, ps303)

ODEProblem with uType ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(S = 1, ES = 2, P = 3)}}} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (0.0, 1.0)
u0: ComponentVector{Float64}(S = 5.0, ES = 0.0, P = 0.0)

@time sol = solve(prob303, Tsit5())

  0.819338 seconds (4.75 M allocations: 316.245 MiB, 10.16% gc time, 100.00% compilation time)

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 40-element Vector{Float64}:
 0.0
 6.665333466693315e-6
 7.331866813362645e-5
 0.000486426499836535
 0.0012463104267083492
 0.0022559163259547044
 0.003638919897606555
 0.005455833654299618
 0.007824959148420711
 0.010840766457437317
 ⋮
 0.6498572831449281
 0.6992361771444529
 0.7455747547375974
 0.7903699863468574
 0.8349926420141908
 0.8806910765547713
 0.9286579911037935
 0.9801092501915283
 1.0
u: 40-element Vector{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(S = 1, ES = 2, P = 3)}}}}:
 ComponentVector{Float64}(S = 5.0, ES = 0.0, P = 0.0)
 ComponentVector{Float64}(S = 4.999000802751873, ES = 0.000999163942258781, P = 3.330586857606668e-8)
 ComponentVector{Float64}(S = 4.989074749914229, ES = 0.010921237104096317, P = 4.012981675171883e-6)
 ComponentVector{Float64}(S = 4.930127945660642, ES = 0.06969993841061337, P = 0.00017211592874508786)
 ComponentVector{Float64}(S = 4.832226260206192, ES = 0.1666948718607426, P = 0.001078867933065927)
 ComponentVector{Float64}(S = 4.720095602191208, ES = 0.2765713014448811, P = 0.003333096363911436)
 ComponentVector{Float64}(S = 4.593086950718923, ES = 0.3988758764938654, P = 0.008037172787212141)
 ComponentVector{Float64}(S = 4.461543546979021, ES = 0.5219980795240161, P = 0.01645837349696343)
 ComponentVector{Float64}(S = 4.332646087873498, ES = 0.6370812621468813, P = 0.030272649979621397)
 ComponentVector{Float64}(S = 4.214625749264197, ES = 0.7343117806650907, P = 0.05106247007071272)
 ⋮
 ComponentVector{Float64}(S = 0.046264116063786755, ES = 0.26918672283275225, P = 4.684549161103462)
 ComponentVector{Float64}(S = 0.020977105513513457, ES = 0.18357572785827508, P = 4.795447166628212)
 ComponentVector{Float64}(S = 0.010414950263948865, ES = 0.12369155788070174, P = 4.86589349185535)
 ComponentVector{Float64}(S = 0.005656467922722563, ES = 0.08276248765840352, P = 4.911581044418875)
 ComponentVector{Float64}(S = 0.003279673934083416, ES = 0.054842165901592244, P = 4.941878160164325)
 ComponentVector{Float64}(S = 0.001969425821800517, ES = 0.0357536784060861, P = 4.962276895772114)
 ComponentVector{Float64}(S = 0.0011916706791066738, ES = 0.0227352496031876, P = 4.976073079717707)
 ComponentVector{Float64}(S = 0.0007101807195001128, ES = 0.013958539864174157, P = 4.985331279416327)
 ComponentVector{Float64}(S = 0.0005834152171657474, ES = 0.011555343089102612, P = 4.987861241693733)

fig303 = Figure()
ax303 = Axis(
    fig303[1, 1],
    xlabel="Time",
    ylabel="Concentration",
    title="Fig. 3.03 (Full model)"
)
lines!(ax303, 0 .. tend, t -> sol(t).S, label="S")
lines!(ax303, 0 .. tend, t -> sol(t).ES, label="ES")
lines!(ax303, 0 .. tend, t -> sol(t).P, label="P")
lines!(ax303, 0 .. tend, t -> model303(sol(t), ps303, t).E, label="E")
axislegend(ax303, position=:rc)
fig303

QSSA of ES complex¶

function model303mm(u, p, t)
    @unpack k1, km1, k2, ET, S0 = p
    @unpack P = u
    S = S0 - P
    ES = (k1 * S * ET) / (km1 + k2 + k1 * S)
    E = ET - ES
    return (; dP=k2 * ES, S, E, ES)
end
function model303mm!(du, u, p, t)
    @unpack dP = model303mm(u, p, t)
    du.P = dP
    nothing
end

ps303mm = ComponentArray(ps303; S0=5.0)
u0303mm = ComponentArray(P=0.0)
prob303mm = ODEProblem(model303mm!, u0303mm, tend, ps303mm)

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

@time sol303mm = solve(prob303mm, Tsit5())

  0.790327 seconds (4.74 M allocations: 315.730 MiB, 9.18% gc time, 100.00% compilation time)

retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 14-element Vector{Float64}:
 0.0
 9.999999999999999e-5
 0.0010999999999999998
 0.011099999999999997
 0.0782815192255917
 0.22560725880659876
 0.40777510082471
 0.5814448248891889
 0.6444132780246575
 0.709328257260174
 0.7626632358601979
 0.8371444047425067
 0.9162616372026308
 1.0
u: 14-element Vector{ComponentArrays.ComponentVector{Float64, Vector{Float64}, Tuple{ComponentArrays.Axis{(P = 1,)}}}}:
 ComponentVector{Float64}(P = 0.0)
 ComponentVector{Float64}(P = 0.0009316710873866917)
 ComponentVector{Float64}(P = 0.010247728722854758)
 ComponentVector{Float64}(P = 0.10334216159542618)
 ComponentVector{Float64}(P = 0.7253458014622206)
 ComponentVector{Float64}(P = 2.061211027540912)
 ComponentVector{Float64}(P = 3.608693155495224)
 ComponentVector{Float64}(P = 4.735254033847741)
 ComponentVector{Float64}(P = 4.921045864324828)
 ComponentVector{Float64}(P = 4.982819834816747)
 ComponentVector{Float64}(P = 4.9957916536455365)
 ComponentVector{Float64}(P = 4.999292966351351)
 ComponentVector{Float64}(P = 4.9998890349369844)
 ComponentVector{Float64}(P = 4.999982400534711)

fig303mm = Figure()
ax303mm = Axis(
    fig303mm[1, 1],
    xlabel="Time",
    ylabel="Concentration",
    title="Fig. 3.03 (QSSA)"
)
lines!(ax303mm, 0 .. tend, t -> sol(t).S, label="S (full)", linestyle=:dash)
lines!(ax303mm, 0 .. tend, t -> sol(t).P, label="P (full)", linestyle=:dash)
lines!(ax303mm, 0 .. tend, t -> model303mm(sol303mm(t), ps303mm, t).S, label="S (MM)")
lines!(ax303mm, 0 .. tend, t -> sol303mm(t).P, label="P (MM)")
axislegend(ax303mm, position=:rc)
fig303mm

This notebook was generated using Literate.jl.