Library of Basic Chemical Reaction Network Models

using Pkg
Pkg.activate(".")
Pkg.add("Catalyst")
Pkg.add("OrdinaryDiffEqDefault")
Pkg.add("StochasticDiffEq")
Pkg.add("JumpProcesses")
Pkg.add("Plots")

Below we will present various simple and established chemical reaction network (CRN) models. Each model is given some brief background, implemented using the @reaction_network DSL, and basic simulations are performed.

Birth-death process

The birth-death process is one of the simplest possible CRN models. It consists of a single component ($X$) which is both produced and degraded at linear rates:

using Catalyst
bd_process = @reaction_network begin
    (p,d), ∅ <--> X
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{p} \mathrm{X} \end{align*} \]

Next, we define simulation conditions. Note that the initial condition is integer-valued (more natural than decimal numbers for jump simulations).

u0 = [:X => 1]
tspan = (0.0, 10.0)
ps = [:p => 1.0, :d => 0.2]

We can now simulate our model using all three interpretations. First, we perform a reaction rate equation-based ODE simulation:

using OrdinaryDiffEqDefault
oprob = ODEProblem(bd_process, u0, tspan, ps)
osol = solve(oprob)

Next, a chemical Langevin equation-based SDE simulation:

using StochasticDiffEq
sprob = SDEProblem(bd_process, u0, tspan, ps)
ssol = solve(sprob, STrapezoid())

Next, a stochastic chemical kinetics-based jump simulation:

using JumpProcesses
jinput = JumpInputs(bd_process, u0, tspan, ps)
jprob = JumpProblem(jinput)
jsol = solve(jprob)

Finally, we plot the results:

using Plots
oplt = plot(osol; title = "Reaction rate equation (ODE)")
splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
plot(oplt, splt, jplt; lw = 3, size=(800,700), layout = (3,1))

Two-state model

The two-state model describes a component (here called $X$) which can exist in two different forms (here called $X₁$ and $X₂$). It switches between these forms at linear rates. First, we simulate the model using both ODEs and SDEs:

using Catalyst, OrdinaryDiffEqDefault, StochasticDiffEq, Plots
two_state_model = @reaction_network begin
    (k₁,k₂), X₁ <--> X₂
end

u0 = [:X₁ => 50.0, :X₂ => 50.0]
tspan = (0.0, 1.0)
ps = [:k₁ => 2.0, :k₂ => 3.0]

oprob = ODEProblem(two_state_model, u0, tspan, ps)
osol = solve(oprob)
oplt = plot(osol; title = "Reaction rate equation (ODE)")

sprob = SDEProblem(two_state_model, u0, tspan, ps)
ssol = solve(sprob, STrapezoid())
splt = plot(ssol; title = "Chemical Langevin equation (SDE)")

plot(oplt, splt; lw = 3, size = (800,550), layout = (2,1))

What is interesting about this model is that it has a conserved quantity, where $X₁ + X₂$ remains constant throughout the simulation (both in deterministic and stochastic cases). We can show this by instead plotting this conserved quantity.

@unpack X₁, X₂ = two_state_model
oplt = plot(osol; idxs = X₁ + X₂, title = "Reaction rate equation (ODE)")
splt = plot(ssol; idxs = X₁ + X₂, title = "Chemical Langevin equation (SDE)")
plot(oplt, splt; lw = 3, ylimit = (99,101), size = (800,450), layout = (2,1))

Michaelis-Menten enzyme kinetics

Michaelis-Menten enzyme kinetics is a simple description of an enzyme ($E$) transforming a substrate ($S$) into a product ($P$). Under certain assumptions, it can be simplified to a single function (a Michaelis-Menten function) and used as a reaction rate. Here we instead present the full system model:

using Catalyst
mm_system = @reaction_network begin
    kB, S + E --> SE
    kD, SE --> S + E
    kP, SE --> P + E
end

\[ \begin{align*} \mathrm{S} + \mathrm{E} &\xrightleftharpoons[\mathtt{kD}]{\mathtt{kB}} \mathrm{\mathtt{SE}} \\ \mathrm{\mathtt{SE}} &\xrightarrow{\mathtt{kP}} \mathrm{P} + \mathrm{E} \end{align*} \]

Next, we perform ODE, SDE, and jump simulations of the model:

u0 = [:S => 301, :E => 100, :SE => 0, :P => 0]
tspan = (0., 100.)
ps = [:kB => 0.00166, :kD => 0.0001, :kP => 0.1]

using OrdinaryDiffEqDefault
oprob = ODEProblem(mm_system, u0, tspan, ps)
osol  = solve(oprob)

using StochasticDiffEq
sprob = SDEProblem(mm_system, u0, tspan, ps)
ssol = solve(sprob, STrapezoid())

using JumpProcesses
jinput = JumpInputs(mm_system, u0, tspan, ps)
jprob = JumpProblem(jinput)
jsol = solve(jprob)

using Plots
oplt = plot(osol; title = "Reaction rate equation (ODE)")
splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1))
fullplt = plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1), fmt = :png,

SIR infection model

The SIR model is the simplest model of the spread of an infectious disease. While the real system is very different from the chemical and cellular processes typically modelled with CRNs, it (and several other epidemiological systems) can be modelled using the same CRN formalism. The SIR model consists of three species: susceptible ($S$), infected ($I$), and removed ($R$) individuals, and two reaction events: infection and recovery.

using Catalyst
sir_model = @reaction_network begin
    α, S + I --> 2I
    β, I --> R
end

\[ \begin{align*} \mathrm{S} + \mathrm{I} &\xrightarrow{\alpha} 2 \mathrm{I} \\ \mathrm{I} &\xrightarrow{\beta} \mathrm{R} \end{align*} \]

First, we perform a deterministic ODE simulation:

using OrdinaryDiffEqDefault, Plots
u0 = [:S => 99, :I => 1, :R => 0]
tspan = (0.0, 500.0)
ps = [:α => 0.001, :β => 0.01]

# Solve ODEs.
oprob = ODEProblem(sir_model, u0, tspan, ps)
osol = solve(oprob)
plot(osol; title = "Reaction rate equation (ODE)", size=(800,350))

Next, we perform 3 different Jump simulations. Note that for the stochastic model, the occurrence of an outbreak is not certain. Rather, there is a possibility that it fizzles out without a noteworthy peak.

using JumpProcesses
jinput = JumpInputs(sir_model, u0, tspan, ps)
jprob = JumpProblem(jinput)

jsol1 = solve(jprob)
jsol2 = solve(jprob)
jsol3 = solve(jprob)

jplt1 = plot(jsol1; title = "Outbreak")
jplt2 = plot(jsol2; title = "Outbreak")
jplt3 = plot(jsol3; title = "No outbreak")
plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1))

Chemical cross-coupling

In chemistry, cross-coupling is when a catalyst combines two substrates to form a product. In this example, the catalyst ($C$) first binds one substrate ($S₁$) to form an intermediary complex ($S₁C$). Next, the complex binds the second substrate ($S₂$) to form another complex ($CP$). Finally, the catalyst releases the now-formed product ($P$). This system is an extended version of the Michaelis-Menten system presented earlier.

using Catalyst
cc_system = @reaction_network begin
    k₁, S₁ + C --> S₁C
    k₂, S₁C + S₂ --> CP
    k₃, CP --> C + P
end

\[ \begin{align*} \mathrm{\mathtt{S_1}} + \mathrm{C} &\xrightarrow{\mathtt{k_1}} \mathrm{\mathtt{S_{1}C}} \\ \mathrm{\mathtt{S_{1}C}} + \mathrm{\mathtt{S_2}} &\xrightarrow{\mathtt{k_2}} \mathrm{\mathtt{CP}} \\ \mathrm{\mathtt{CP}} &\xrightarrow{\mathtt{k_3}} \mathrm{C} + \mathrm{P} \end{align*} \]

Below, we perform a simple deterministic ODE simulation of the system. Next, we plot both:

• The concentration of the substrates and the product.
• The concentration of the catalyst and the intermediaries.

In two separate plots.

using OrdinaryDiffEqDefault, Plots
u0 = [:S₁ => 1.0, :C => 0.05, :S₂ => 1.2, :S₁C => 0.0, :CP => 0.0, :P => 0.0]
tspan = (0., 15.)
ps = [:k₁ => 5.0, :k₂ => 5.0, :k₃ => 100.0]

# solve ODEs
oprob = ODEProblem(cc_system, u0, tspan, ps)
osol  = solve(oprob)

plt1 = plot(osol; idxs = [:S₁, :S₂, :P], title = "Substrate and product dynamics")
plt2 = plot(osol; idxs = [:C, :S₁C, :CP], title = "Catalyst and intermediaries dynamics")
plot(plt1, plt2; lw = 3, size = (800,600), layout = (2,1))

The Wilhelm model

The Wilhelm model was introduced in Wilhelm (2009) as the smallest CRN model (with constant rates) that exhibits bistability.

using Catalyst
wilhelm_model = @reaction_network begin
    k1, Y --> 2X
    k2, 2X --> X + Y
    k3, X + Y --> Y
    k4, X --> 0
end

\[ \begin{align*} \mathrm{Y} &\xrightarrow{\mathtt{k1}} 2 \mathrm{X} \\ 2 \mathrm{X} &\xrightarrow{\mathtt{k2}} \mathrm{X} + \mathrm{Y} \\ \mathrm{X} + \mathrm{Y} &\xrightarrow{\mathtt{k3}} \mathrm{Y} \\ \mathrm{X} &\xrightarrow{\mathtt{k4}} \varnothing \end{align*} \]

We can simulate the model for two different initial conditions, demonstrating the existence of two different stable steady states.

using OrdinaryDiffEqDefault, Plots
u0_1 = [:X => 1.5, :Y => 0.5]
u0_2 = [:X => 2.5, :Y => 0.5]
tspan = (0., 10.)
ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5]

oprob1 = ODEProblem(wilhelm_model, u0_1, tspan, ps)
oprob2 = ODEProblem(wilhelm_model, u0_2, tspan, ps)
osol1 = solve(oprob1)
osol2 = solve(oprob2)
plot(osol1; lw = 4, idxs = :X, label = "X(0) = 1.5")
plot!(osol2; lw = 4, idxs = :X, label = "X(0) = 2.5", yguide = "X", size = (800,350))

Simple self-activation loop

The simplest self-activation loop consists of a single species (here called $X$) which activates its own production. If its production rate is modelled as a Hill function with $n>1$, the system may exhibit bistability.

using Catalyst
sa_loop = @reaction_network begin
    v₀ + hill(X,v,K,n), ∅ --> X
    d, X --> ∅
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{\mathtt{v_0} + \frac{X^{n} v}{K^{n} + X^{n}}} \mathrm{X} \end{align*} \]

A simple example of such a loop is a transcription factor which activates its own gene. Here, $v₀$ represents a basic transcription rate (leakage) in the absence of the transcription factor.

We simulate the self-activation loop from a single initial condition using both deterministic (ODE) and stochastic (jump) simulations. We note that while the deterministic simulation reaches a single steady state, the stochastic one switches between two different states.

using JumpProcesses, OrdinaryDiffEqDefault, Plots
u0 = [:X => 4]
tspan = (0.0, 1000.0)
ps = [:v₀ => 0.1, :v => 2.0, :K => 10.0, :n => 2, :d => 0.1]

oprob = ODEProblem(sa_loop, u0, tspan, ps)
osol = solve(oprob)

jinput = JumpInputs(sa_loop, u0, tspan, ps)
jprob = JumpProblem(jinput)
jsol = solve(jprob)

fplt = plot(osol; lw = 3, label = "Reaction rate equation (ODE)")
plot!(fplt, jsol; lw = 3, label = "Stochastic chemical kinetics (Jump)", yguide = "X", size = (800,350))

The Brusselator

The Brusselator is a well-known (theoretical) CRN model able to produce oscillations (its name is a portmanteau of "Brussels" and "oscillator").

using Catalyst
brusselator = @reaction_network begin
    A, ∅ --> X
    1, 2X + Y --> 3X
    B, X --> Y
    1, X --> ∅
end

\[ \begin{align*} \varnothing &\xrightarrow{A} \mathrm{X} \\ 2 \mathrm{X} + \mathrm{Y} &\xrightarrow{1} 3 \mathrm{X} \\ \mathrm{X} &\xrightarrow{B} \mathrm{Y} \\ \mathrm{X} &\xrightarrow{1} \varnothing \end{align*} \]

It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this result is based on models generated when combinatorial adjustment of rates is not performed. Since Catalyst automatically perform these adjustments, and one reaction contains a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in the second case the system can generate oscillations.

using OrdinaryDiffEqDefault, Plots
u0 = [:X => 1.0, :Y => 1.0]
tspan = (0., 50.)
ps1 = [:A => 1.0, :B => 1.0]
ps2 = [:A => 1.0, :B => 1.8]

oprob1 = ODEProblem(brusselator, u0, tspan, ps1)
oprob2 = ODEProblem(brusselator, u0, tspan, ps2)
osol1  = solve(oprob1)
osol2  = solve(oprob2)
oplt1 = plot(osol1; title = "No Oscillation")
oplt2 = plot(osol2; title = "Oscillation")

plot(oplt1, oplt2; lw = 3, size = (800,600), layout = (2,1))

The Repressilator

The Repressilator was introduced in Elowitz & Leibler (2000) as a simple system that can generate oscillations (most notably, they demonstrated this both in a model and in a synthetic in vivo implementation in Escherichia col). It consists of three genes, repressing each other in a cycle. Here, we will implement it using three species ($X$, $Y$, and $Z$) whose production rates are (repressing) Hill functions.

using Catalyst
repressilator = @reaction_network begin
    hillr(Z,v,K,n), ∅ --> X
    hillr(X,v,K,n), ∅ --> Y
    hillr(Y,v,K,n), ∅ --> Z
    d, (X, Y, Z) --> ∅
end

\[ \begin{align*} \varnothing &\xrightarrow{\frac{K^{n} v}{K^{n} + Z^{n}}} \mathrm{X} \\ \varnothing &\xrightarrow{\frac{K^{n} v}{K^{n} + X^{n}}} \mathrm{Y} \\ \varnothing &\xrightarrow{\frac{K^{n} v}{K^{n} + Y^{n}}} \mathrm{Z} \\ \mathrm{X} &\xrightarrow{d} \varnothing \\ \mathrm{Y} &\xrightarrow{d} \varnothing \\ \mathrm{Z} &\xrightarrow{d} \varnothing \end{align*} \]

Whether the Repressilator oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model can sustain oscillations in both cases. This is an example of the phenomena of noise-induced oscillation.

using OrdinaryDiffEqDefault, StochasticDiffEq, Plots
u0 = [:X => 50.0, :Y => 15.0, :Z => 15.0]
tspan = (0., 200.)
ps1 = [:v => 10.0, :K => 20.0, :n => 3, :d => 0.1]
ps2 = [:v => 10.0, :K => 50.0, :n => 3, :d => 0.1]

oprob1 = ODEProblem(repressilator, u0, tspan, ps1)
oprob2 = ODEProblem(repressilator, u0, tspan, ps2)
osol1  = solve(oprob1)
osol2  = solve(oprob2)
oplt1 = plot(osol1; title = "Oscillation (ODE, K = 20)")
oplt2 = plot(osol2; title = "No oscillation (ODE, K = 50)")

sprob1 = SDEProblem(repressilator, u0, tspan, ps1)
sprob2 = SDEProblem(repressilator, u0, tspan, ps2)
ssol1  = solve(sprob1, STrapezoid())
ssol2  = solve(sprob2, STrapezoid())
splt1 = plot(ssol1; title = "Oscillation (SDE, K = 20)")
splt2 = plot(ssol2; title = "Oscillation (SDE, K = 50)")

plot(oplt1, oplt2, splt1, splt2; lw = 2, layout = (2,2), size = (800,600))

The Willamowski–Rössler model

The Willamowski–Rössler model was introduced in Willamowski & Rössler (1979) as an example of a simple CRN model which exhibits chaotic behaviours. This means that small changes in initial conditions can produce relatively large changes in the system's trajectory.

using Catalyst
wr_model = @reaction_network begin
    k1, 2X --> 3X
    k2, X --> 2X
    k3, Z + 2X --> 2Z
    k4, Y + X --> 2Y
    k5, Y --> ∅
    k6, 2Z --> ∅
    k7, Z --> ∅
end

\[ \begin{align*} 2 \mathrm{X} &\xrightleftharpoons[\mathtt{k2}]{\mathtt{k1}} 3 \mathrm{X} \\ \mathrm{Z} + 2 \mathrm{X} &\xrightarrow{\mathtt{k3}} 2 \mathrm{Z} \\ \mathrm{Y} + \mathrm{X} &\xrightarrow{\mathtt{k4}} 2 \mathrm{Y} \\ \mathrm{Y} &\xrightarrow{\mathtt{k5}} \varnothing \\ 2 \mathrm{Z} &\xrightarrow{\mathtt{k6}} \varnothing \\ \mathrm{Z} &\xrightarrow{\mathtt{k7}} \varnothing \end{align*} \]

Here we simulate the model for a single initial condition, showing both time-state space and phase space how it reaches a strange attractor.

using OrdinaryDiffEqDefault, Plots
u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5]
tspan = (0.0, 50.0)
p = [:k1 => 2.1, :k2 => 0.7, :k3 => 2.9, :k4 => 1.1, :k5 => 1.0, :k6 => 0.5, :k7 => 2.7]
oprob = ODEProblem(wr_model, u0, tspan, p)
sol = solve(oprob)

plt1 = plot(sol; title = "Time-state space")
plt2 = plot(sol; idxs = (:X, :Y, :Z), title = "Phase space")
plot(plt1, plt2; layout = (1,2), size = (800,400))