Fig 6.14

Fig 6.14#

Model of E. coli chemotaxis signalling pathway

using ModelingToolkit
using Catalyst
using OrdinaryDiffEq
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]]) |> complete

┌ 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);
sol = solve(prob)

┌ 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/8S2W1/src/systems/diffeqs/abstractodesystem.jl:1441

retcode: Success
Interpolation: 3rd order Hermite
t: 105-element Vector{Float64}:
  0.0
  0.0
  0.05968737770312525
  0.08662460033554424
  0.12694391787581313
  0.16950791095146842
  0.2229946265995895
  0.28241446330146397
  0.3413393099804022
  0.39629908274591924
  ⋮
 34.00753385281169
 34.97233280943664
 36.03648174376251
 37.33758933001082
 38.83290672234507
 40.68813264306904
 42.956499116505604
 45.95257092480216
 50.0
u: 105-element Vector{Vector{Float64}}:
 [0.036, 0.0595, 1.5593, 0.7356]
 [0.036, 0.0595, 1.5593, 0.7356]
 [0.03595027525385366, 0.05953953220886678, 1.5593219305751664, 0.7355997935085642]
 [0.03601289608672456, 0.05947514303222826, 1.5593068605221378, 0.7355997452996131]
 [0.036022073388818514, 0.05946463801131664, 1.5593088555837615, 0.7355996283853471]
 [0.03602437773678273, 0.05946174455183917, 1.5593129524853817, 0.7355995010070753]
 [0.036023848727154105, 0.05946199844888866, 1.5593188688752477, 0.7355993390621877]
 [0.036020335741447994, 0.05946541683014946, 1.559326004755733, 0.7355991575254756]
 [0.0360109087395074, 0.059474816689258594, 1.5593345280134596, 0.7355989742805211]
 [0.03599997865292204, 0.05948573915286399, 1.559342857705766, 0.7355988021675578]
 ⋮
 [0.033118745750131566, -0.008363190194611653, 3.33627050096727, 0.737046439075847]
 [0.03381860741903603, -0.009063051863577583, 3.4259384369778125, 0.7371243010757769]
 [0.03433985395185947, -0.009584298396307667, 3.492655671737599, 0.7371858177820777]
 [0.03474484195035795, -0.009989286399439662, 3.5444349293865143, 0.7372383500371893]
 [0.03501673769086339, -0.010261182135087733, 3.57914417946151, 0.7372797060979512]
 [0.035195198330194535, -0.010439642774645616, 3.601863363932512, 0.7373152235911977]
 [0.03529536939888693, -0.010539813843335709, 3.6145359767544063, 0.7373466319031647]
 [0.0353461348131183, -0.010590579257581754, 3.620844276609218, 0.7373792751630395]
 [0.035367616730279924, -0.010612061174729173, 3.62335021004333, 0.7374171909033692]

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.