Fig 7.7

Fig 7.7#

model of lac operon in E. coli

using Catalyst
using ModelingToolkit
using DifferentialEquations
using Plots
Plots.default(linewidth=2)

rn = @reaction_network begin
    @parameters Le(t)
    a1 / (1 + RToverK1 * (K2 / (K2 + L))^4), 0 --> M
    δM, M --> 0
    (c1 * M, δY), 0 <--> Y
    mm(Le, kL * Y, KML), 0 --> L
    mm(L, 2 * kg * (Y / 4), KMg), L => 0
    δL, L --> 0
end

\[\begin{split} \begin{align*} \varnothing &\xrightleftharpoons[\mathtt{{\delta}M}]{\frac{\mathtt{a1}}{1 + \left( \frac{\mathtt{K2}}{\mathtt{K2} + L} \right)^{4} \mathtt{RToverK1}}} \mathrm{M} \\ \varnothing &\xrightleftharpoons[\mathtt{{\delta}Y}]{M \mathtt{c1}} \mathrm{Y} \\ \varnothing &\xrightleftharpoons[\frac{\frac{1}{2} L Y \mathtt{kg}}{\mathtt{KMg} + L}]{\frac{Y \mathtt{kL} \mathrm{\mathtt{Le}}\left( t \right)}{\mathtt{KML} + \mathrm{\mathtt{Le}}\left( t \right)}} \mathrm{L} \\ \mathrm{L} &\xrightarrow{\mathtt{{\delta}L}} \varnothing \end{align*} \end{split}\]

setdefaults!(rn, [
    :δM => 0.48,
    :δY => 0.03,
    :δL => 0.02,
    :a1 => 0.29,
    :K2 => 2.92 * 1e6,
    :RToverK1 => 213.2,
    :c1 => 18.8,
    :kL => 6 * 1e4,
    :KML => 680,
    :kg => 3.6 * 1e3,
    :KMg => 7 * 1e5,
    :Le => 0.0,
    :M => 0.01,
    :Y => 0.1,
    :L => 0.0
])

@unpack Le = rn
discrete_events = [[500.0] => [Le ~ 50], [1000.0] => [Le ~ 100], [1500.0] => [Le ~ 150], [2000.0] => [Le ~ 0]]
osys = convert(ODESystem, rn; discrete_events, remove_conserved=true) |> 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} M\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{a1}}{1 + \left( \frac{\mathtt{K2}}{\mathtt{K2} + L\left( t \right)} \right)^{4} \mathtt{RToverK1}} - M\left( t \right) \mathtt{{\delta}M} \\ \frac{\mathrm{d} Y\left( t \right)}{\mathrm{d}t} &= \mathtt{c1} M\left( t \right) - Y\left( t \right) \mathtt{{\delta}Y} \\ \frac{\mathrm{d} L\left( t \right)}{\mathrm{d}t} &= - mm\left( L\left( t \right), \frac{1}{2} \mathtt{kg} Y\left( t \right), \mathtt{KMg} \right) + mm\left( \mathtt{Le}\left( t \right), \mathtt{kL} Y\left( t \right), \mathtt{KML} \right) - L\left( t \right) \mathtt{{\delta}L} \end{align} \end{split}\]

Fig 7.07 (A)#

tend = 2500.0
prob = ODEProblem(osys, [], tend)
sol = solve(prob)

@unpack M, Y, L, Le = osys
plot(sol, idxs=[Y, Le], xlabel="Time (min)", title="Fig 7.7 (A)", label=["β-galactosidase monomer" "External lactose (μM)"])

../_images/998920a9c1e443e3e6fb5a47020a00de494d9e446871600ff759aea80d573dd3.png

Fig 7.07 (B)#

Compare the original model and the modified model

rn_mod = @reaction_network begin
    @parameters Le(t)
    a1 / (1 + RToverK1 * (K2 / (K2 + L))^4), 0 --> M
    δM, M --> 0
    (c1 * M, δY), 0 <--> Y
    mm(Le, kL * 4 * Enz, KML), 0 --> L
    mm(L, 2 * kg * Enz, KMg), L ⇒ 0
    δL, L --> 0
end

\[\begin{split} \begin{align*} \varnothing &\xrightleftharpoons[\mathtt{{\delta}M}]{\frac{\mathtt{a1}}{1 + \left( \frac{\mathtt{K2}}{\mathtt{K2} + L} \right)^{4} \mathtt{RToverK1}}} \mathrm{M} \\ \varnothing &\xrightleftharpoons[\mathtt{{\delta}Y}]{M \mathtt{c1}} \mathrm{Y} \\ \varnothing &\xrightleftharpoons[\frac{2 \mathtt{Enz} L \mathtt{kg}}{\mathtt{KMg} + L}]{\frac{4 \mathtt{Enz} \mathtt{kL} \mathrm{\mathtt{Le}}\left( t \right)}{\mathtt{KML} + \mathrm{\mathtt{Le}}\left( t \right)}} \mathrm{L} \\ \mathrm{L} &\xrightarrow{\mathtt{{\delta}L}} \varnothing \end{align*} \end{split}\]

setdefaults!(rn_mod, [
    :δM => 0.48,
    :δY => 0.03,
    :δL => 0.02,
    :a1 => 0.29,
    :K2 => 2.92 * 1e6,
    :RToverK1 => 213.2,
    :c1 => 18.8,
    :kL => 6 * 1e4,
    :KML => 680,
    :kg => 3.6 * 1e3,
    :KMg => 7 * 1e5,
    :Le => 0.0,
    :M => 0.01,
    :Y => 0.1,
    :L => 0.0,
    :Enz => 40.0
])

osys_mod = convert(ODESystem, rn_mod; remove_conserved=true) |> 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} M\left( t \right)}{\mathrm{d}t} &= \frac{\mathtt{a1}}{1 + \left( \frac{\mathtt{K2}}{\mathtt{K2} + L\left( t \right)} \right)^{4} \mathtt{RToverK1}} - M\left( t \right) \mathtt{{\delta}M} \\ \frac{\mathrm{d} Y\left( t \right)}{\mathrm{d}t} &= \mathtt{c1} M\left( t \right) - Y\left( t \right) \mathtt{{\delta}Y} \\ \frac{\mathrm{d} L\left( t \right)}{\mathrm{d}t} &= - mm\left( L\left( t \right), 2 \mathtt{Enz} \mathtt{kg}, \mathtt{KMg} \right) + mm\left( \mathtt{Le}\left( t \right), 4 \mathtt{Enz} \mathtt{kL}, \mathtt{KML} \right) - L\left( t \right) \mathtt{{\delta}L} \end{align} \end{split}\]

prob = SteadyStateProblem(rn, [])
prob_mod = SteadyStateProblem(rn_mod, [])

@unpack Le, Y = osys
lerange = range(0, 100, 101)

eprob = EnsembleProblem(prob;
    prob_func=(prob, i, repeat) -> remake(prob, p=[Le=>lerange[i]]),
    output_func=(sol, i) -> (sol[Y] / 4, false)
)

eprob_mod = EnsembleProblem(prob_mod;
    prob_func=(prob, i, repeat) -> remake(prob, p=[Le=>lerange[i]]),
    output_func=(sol, i) -> (sol[Y] / 4, false)
)

sim = solve(eprob, DynamicSS(Rodas5()); trajectories=length(lerange), abstol=1e-9, reltol=1e-9)
sim_mod = solve(eprob_mod, DynamicSS(Rodas5()); trajectories=length(lerange), abstol=1e-9, reltol=1e-9)

fig = plot(lerange, sim.u, label="Original")
fig = plot!(fig, lerange, sim_mod.u, label="Modified")
fig = plot!(fig,
    xlabel="External lactose concentration (μM)",
    ylabel="β-galactosidase",
    title="Fig 7.7 (B)"
)

../_images/25651d53eaa9bb213bffb9e9eb071dfe7f9f97efd3f89f34eea5d907530dcdcb.png

This notebook was generated using Literate.jl.

Fig 7.7

Contents

Fig 7.7#

Fig 7.07 (A)#

Fig 7.07 (B)#