Fig 6.14

Fig 6.14#

Model of E. coli chemotaxis signalling pathway

using ModelingToolkit
using Catalyst
using DifferentialEquations
using Plots
Plots.default(linewidth=2)
rn = @reaction_network begin
    @parameters L(t)
    mm(Am, k1 * BP, KM1), Am => A
    mm(AmL, k2 * BP, KM2), AmL => AL
    km1 * R , A => Am
    km2 * R , AL => AmL
    (k3 * L, km3), Am <--> AmL
    (k4 * L, km4), A <--> AL
    (k5 * Am, km5), B <--> BP
end
\[\begin{split} \begin{align*} \mathrm{\mathtt{Am}} &\xrightarrow{\frac{\mathtt{Am} \mathtt{BP} \mathtt{k1}}{\mathtt{Am} + \mathtt{KM1}}} \mathrm{A} \\ \mathrm{\mathtt{AmL}} &\xrightarrow{\frac{\mathtt{AmL} \mathtt{BP} \mathtt{k2}}{\mathtt{AmL} + \mathtt{KM2}}} \mathrm{\mathtt{AL}} \\ \mathrm{A} &\xrightarrow{R \mathtt{km1}} \mathrm{\mathtt{Am}} \\ \mathrm{\mathtt{AL}} &\xrightarrow{R \mathtt{km2}} \mathrm{\mathtt{AmL}} \\ \mathrm{\mathtt{Am}} &\xrightleftharpoons[\mathtt{km3}]{\mathtt{k3} L\left( t \right)} \mathrm{\mathtt{AmL}} \\ \mathrm{A} &\xrightleftharpoons[\mathtt{km4}]{\mathtt{k4} L\left( t \right)} \mathrm{\mathtt{AL}} \\ \mathrm{B} &\xrightleftharpoons[\mathtt{km5}]{\mathtt{Am} \mathtt{k5}} \mathrm{\mathtt{BP}} \end{align*} \end{split}\]
setdefaults!(rn, [
    :Am => 0.0360,
    :AmL => 1.5593,
    :A => 0.0595,
    :AL => 0.3504,
    :B => 0.7356,
    :BP => 0.2644,
    :k1 => 200,
    :k2 => 1,
    :k3 => 1,
    :km1 => 1,
    :km2 => 1,
    :km3 => 1,
    :k5 => 0.05,
    :km5 => 0.005,
    :k4 => 1,
    :km4 => 1,
    :KM1 => 1,
    :KM2 => 1,
    :R => 1,
    :L => 20
])

@unpack L = rn
osys = convert(ODESystem, rn; remove_conserved = true, discrete_events = [[10] => [L ~ 40], [30] => [L ~ 80]]) |> structural_simplify
┌ Warning: You are creating a system or problem while eliminating conserved quantities. Please note,
│         due to limitations / design choices in ModelingToolkit if you use the created system to
│         create a problem (e.g. an `ODEProblem`), or are directly creating a problem, you *should not*
│         modify that problem's initial conditions for species (e.g. using `remake`). Changing initial
│         conditions must be done by creating a new Problem from your reaction system or the
│         ModelingToolkit system you converted it into with the new initial condition map.
│         Modification of parameter values is still possible, *except* for the modification of any
│         conservation law constants (Γ), which is not possible. You might
│         get this warning when creating a problem directly.
│ 
│         You can remove this warning by setting `remove_conserved_warn = false`.
└ @ Catalyst ~/.julia/packages/Catalyst/48wH3/src/reactionsystem_conversions.jl:456
\[\begin{split} \begin{align} \frac{\mathrm{d} \mathtt{Am}\left( t \right)}{\mathrm{d}t} &= - mm\left( \mathtt{Am}\left( t \right), \mathtt{k1} \left( - B\left( t \right) + \Gamma_{2} \right), \mathtt{KM1} \right) + R \mathtt{km1} + \mathtt{km3} \mathtt{AmL}\left( t \right) - \mathtt{k3} \mathtt{Am}\left( t \right) L\left( t \right) \\ \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= mm\left( \mathtt{Am}\left( t \right), \mathtt{k1} \left( - B\left( t \right) + \Gamma_{2} \right), \mathtt{KM1} \right) - R \mathtt{km1} + \mathtt{km4} \left( - \mathtt{AmL}\left( t \right) - \mathtt{Am}\left( t \right) - A\left( t \right) + \Gamma_{1} \right) - \mathtt{k4} A\left( t \right) L\left( t \right) \\ \frac{\mathrm{d} \mathtt{AmL}\left( t \right)}{\mathrm{d}t} &= - mm\left( \mathtt{AmL}\left( t \right), \mathtt{k2} \left( - B\left( t \right) + \Gamma_{2} \right), \mathtt{KM2} \right) + R \mathtt{km2} - \mathtt{km3} \mathtt{AmL}\left( t \right) + \mathtt{k3} \mathtt{Am}\left( t \right) L\left( t \right) \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} &= \mathtt{km5} \left( - B\left( t \right) + \Gamma_{2} \right) - \mathtt{k5} \mathtt{Am}\left( t \right) B\left( t \right) \end{align} \end{split}\]
observed(osys)
\[\begin{split} \begin{align} \mathtt{AL}\left( t \right) &= - \mathtt{AmL}\left( t \right) - \mathtt{Am}\left( t \right) - A\left( t \right) + \Gamma_{1} \\ \mathtt{BP}\left( t \right) &= - B\left( t \right) + \Gamma_{2} \end{align} \end{split}\]
equations(osys)
\[\begin{split} \begin{align} \frac{\mathrm{d} \mathtt{Am}\left( t \right)}{\mathrm{d}t} &= - mm\left( \mathtt{Am}\left( t \right), \mathtt{k1} \left( - B\left( t \right) + \Gamma_{2} \right), \mathtt{KM1} \right) + R \mathtt{km1} + \mathtt{km3} \mathtt{AmL}\left( t \right) - \mathtt{k3} \mathtt{Am}\left( t \right) L\left( t \right) \\ \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= mm\left( \mathtt{Am}\left( t \right), \mathtt{k1} \left( - B\left( t \right) + \Gamma_{2} \right), \mathtt{KM1} \right) - R \mathtt{km1} + \mathtt{km4} \left( - \mathtt{AmL}\left( t \right) - \mathtt{Am}\left( t \right) - A\left( t \right) + \Gamma_{1} \right) - \mathtt{k4} A\left( t \right) L\left( t \right) \\ \frac{\mathrm{d} \mathtt{AmL}\left( t \right)}{\mathrm{d}t} &= - mm\left( \mathtt{AmL}\left( t \right), \mathtt{k2} \left( - B\left( t \right) + \Gamma_{2} \right), \mathtt{KM2} \right) + R \mathtt{km2} - \mathtt{km3} \mathtt{AmL}\left( t \right) + \mathtt{k3} \mathtt{Am}\left( t \right) L\left( t \right) \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} &= \mathtt{km5} \left( - B\left( t \right) + \Gamma_{2} \right) - \mathtt{k5} \mathtt{Am}\left( t \right) B\left( t \right) \end{align} \end{split}\]
tend = 50.0
prob = ODEProblem(osys, [], tend)
┌ Warning: Initialization system is overdetermined. 2 equations for 0 unknowns. Initialization will default to using least squares. `SCCNonlinearProblem` can only be used for initialization of fully determined systems and hence will not be used here. To suppress this warning pass warn_initialize_determined = false. To make this warning into an error, pass fully_determined = true
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/klLLV/src/systems/diffeqs/abstractodesystem.jl:1354
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 50.0)
u0: 4-element Vector{Float64}:
 0.036
 0.0595
 1.5593
 0.7356
sol = solve(prob)

plot(sol, idxs=[osys.Am], title="Fig 6.14", xlabel="Time", ylabel="Active CheA ([Am])", ylims=(0.01, 0.04), legend=false)
../_images/d732d6c7ebb966735d89b9e16b2d9691d0c1b769e8d3f7d2dd80a82e476d471d.png

This notebook was generated using Literate.jl.