Mathematical Models Catalyst can Generate

We now describe the types of mathematical models that Catalyst can generate from chemical reaction networks (CRNs), corresponding to reaction rate equation (RRE) ordinary differential equation (ODE) models, Chemical Langevin equation (CLE) stochastic differential equation (SDE) models, and stochastic chemical kinetics (jump process) models. For each we show the abstract representations for the models that Catalyst can support, along with concrete examples. Note that we restrict ourselves to models involving only chemical reactions, and do not consider more general models that Catalyst can support such as coupling in non-reaction ODEs, algebraic equations, or events. Please see the broader documentation for more details on how Catalyst supports such functionality.

General Chemical Reaction Notation

Suppose we have a reaction network with $K$ reactions and $M$ species, with the species labeled by $S_1$, $S_2$, $\dots$, $S_M$. We denote by

\[\mathbf{X}(t) = \begin{pmatrix} X_1(t) \\ \vdots \\ X_M(t)) \end{pmatrix}\]

the state vector for the amount of each species, i.e. $X_m(t)$ represents the amount of species $S_m$ at time $t$. This could be either a concentration or a count (i.e. "number of molecules" units), but for consistency between modeling representations we will assume the latter in the remainder of this introduction.

The $k$th chemical reaction is given by

\[\alpha_1^k S_1 + \alpha_2^k S_2 + \dots \alpha_M^k S_M \to \beta_1^k S_1 + \beta_2^k S_2 + \dots \beta_M^k S_M\]

with $\alpha^k = (\alpha_1^k,\dots,\alpha_M^k)$ its substrate stoichiometry vector, $\beta^k = (\beta_1^k,\dots,\beta_M^k)$ its product stoichiometry vector, and $\nu^k = \beta^k - \alpha^k$ its net stoichiometry vector. $\nu^k$ corresponds to the change in $\mathbf{X}(t)$ when reaction $k$ occurs, i.e. $\mathbf{X}(t) \to \mathbf{X}(t) + \nu^k$. Along with the stoichiometry vectors, we assume each reaction has a reaction rate law (ODEs/SDEs) or propensity (jump process) function, $a_k(\mathbf{X}(t),t)$.

As explained in the Catalyst introduction, for a mass action reaction where the preceding reaction has a fixed rate constant, $k$, this function would be the rate law

\[a_k(\mathbf{X}(t)) = k \prod_{m=1}^M \frac{(X_m(t))^{\alpha_m^k}}{\alpha_m^k!},\]

for RRE ODE and CLE SDE models, and the propensity function

\[a_k(\mathbf{X}(t)) = k \prod_{m=1}^M \frac{X_m(t) (X_m(t)-1) \dots (X_m(t)-\alpha_m^k+1)}{\alpha_m^k!},\]

for stochastic chemical kinetics jump process models.

Rate Law vs. Propensity Example:

For the reaction $2A + B \overset{k}{\to} 3 C$ we would have

\[\mathbf{X}(t) = (A(t), B(t), C(t))\]

with $\alpha_1 = 2$, $\alpha_2 = 1$, $\alpha_3 = 0$, $\beta_1 = 0$, $\beta_2 = 0$, $\beta_3 = 3$, $\nu_1 = -2$, $\nu_2 = -1$, and $\nu_3 = 3$. For an ODE/SDE model we would have the rate law

\[a(\mathbf{X}(t)) = \frac{k}{2} A^2 B\]

while for a jump process model we would have the propensity function

\[a(\mathbf{X}(t)) = \frac{k}{2} A (A-1) B.\]

Note, if the combinatoric factors are already included in one's rate constants, the implicit rescaling of rate constants can be disabled through use of the combinatoric_ratelaws = false argument to Base.convert or whatever Problem is being generated, i.e.

rn = @reaction_network ...
osys = convert(ODESystem, rn; combinatoric_ratelaws = false)
oprob = ODEProblem(osys, ...)

# or
oprob = ODEProblem(rn, ...; combinatoric_ratelaws = false)

In this case our ODE/SDE rate law would be

\[a(\mathbf{X}(t)) = k A^2 B\]

while the jump process propensity function is

\[a(\mathbf{X}(t)) = k A (A-1) B.\]

One can also specify during system construction that by default combinatoric scalings should be disabled, i.e.

using Catalyst
rn = @reaction_network begin
    @combinatoric_ratelaws false
    k, 3A + 2B --> A + 3D
end
osys = convert(ODESystem, rn)

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= - 2 \left( A\left( t \right) \right)^{3} \left( B\left( t \right) \right)^{2} k \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} &= - 2 \left( A\left( t \right) \right)^{3} \left( B\left( t \right) \right)^{2} k \\ \frac{\mathrm{d} D\left( t \right)}{\mathrm{d}t} &= 3 \left( A\left( t \right) \right)^{3} \left( B\left( t \right) \right)^{2} k \end{align} \]

Reaction Rate Equation (RRE) ODE Models

The RRE ODE models Catalyst creates for a general system correspond to the coupled system of ODEs given by

\[\frac{d X_m}{dt} =\sum_{k=1}^K \nu_m^k a_k(\mathbf{X}(t),t), \quad m = 1,\dots,M.\]

These models can be generated by creating ODEProblems from Catalyst ReactionSystems, and solved using the solvers in OrdinaryDiffEq.jl. Similarly, creating NonlinearProblems or SteadyStateProblems will generate the coupled algebraic system of steady-state equations associated with a RRE ODE model, i.e.

\[0 =\sum_{k=1}^K \nu_m^k a_k(\bar{\mathbf{X}}), \quad m = 1,\dots,M\]

for a steady-state $\bar{\mathbf{X}}$. Note, here we have assumed the rate laws are autonomous so that the equations are well-defined.

RRE ODE Example

Let's see the generated ODEs for the following network

using Catalyst, ModelingToolkit, Latexify
rn = @reaction_network begin
    k₁, 2A + B --> 3C
    k₂, A --> 0
    k₃, 0 --> A
end
osys = convert(ODESystem, rn)

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= \mathtt{k{_3}} - \mathtt{k{_2}} A\left( t \right) - \left( A\left( t \right) \right)^{2} \mathtt{k{_1}} B\left( t \right) \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} &= - \frac{1}{2} \left( A\left( t \right) \right)^{2} \mathtt{k{_1}} B\left( t \right) \\ \frac{\mathrm{d} C\left( t \right)}{\mathrm{d}t} &= \frac{3}{2} \left( A\left( t \right) \right)^{2} \mathtt{k{_1}} B\left( t \right) \end{align} \]

Likewise, the following drops the combinatoric scaling factors, giving unscaled ODEs

osys = convert(ODESystem, rn; combinatoric_ratelaws = false)

\[ \begin{align} \frac{\mathrm{d} A\left( t \right)}{\mathrm{d}t} &= \mathtt{k{_3}} - \mathtt{k{_2}} A\left( t \right) - 2 \left( A\left( t \right) \right)^{2} \mathtt{k{_1}} B\left( t \right) \\ \frac{\mathrm{d} B\left( t \right)}{\mathrm{d}t} &= - \left( A\left( t \right) \right)^{2} \mathtt{k{_1}} B\left( t \right) \\ \frac{\mathrm{d} C\left( t \right)}{\mathrm{d}t} &= 3 \left( A\left( t \right) \right)^{2} \mathtt{k{_1}} B\left( t \right) \end{align} \]

Chemical Langevin Equation (CLE) SDE Models

The CLE SDE models Catalyst creates for a general system correspond to the coupled system of SDEs given by

\[d X_m = \sum_{k=1}^K \nu_m^k a_k(\mathbf{X}(t),t) dt + \sum_{k=1}^K \nu_m^k \sqrt{a_k(\mathbf{X}(t),t)} dW_k(t), \quad m = 1,\dots,M,\]

where each $W_k(t)$ represents an independent, standard Brownian Motion. Realizations of these processes can be generated by creating SDEProblems from Catalyst ReactionSystems, and sampling the processes using the solvers in StochasticDiffEq.jl.

CLE SDE Example

Consider the same network as above,

rn = @reaction_network begin
    k₁, 2A + B --> 3C
    k₂, A --> 0
    k₃, 0 --> A
end

We obtain the CLE SDEs

\[\begin{align} dA(t) &= \left(- k_1 A^{2} B - k_2 A + k_3 \right) dt - 2 \sqrt{\tfrac{k_1}{2} A^{2} B} \, dW_1(t) - \sqrt{k_2 A} \, dW_2(t) + \sqrt{k_3} \, dW_3(t) \\ dB(t) &= - \frac{k_1}{2} A^{2} B \, dt - \sqrt{\frac{k_1}{2} A^{2} B} \, dW_1(t) \\ dC(t) &= \frac{3}{2} k_1 A^{2} B \, dt + 3 \sqrt{\frac{k_1}{2} A^{2} B} \, dW_1(t). \end{align}\]

Stochastic Chemical Kinetics Jump Process Models

The stochastic chemical kinetics jump process models Catalyst creates for a general system correspond to the coupled system of jump processes, in the time change representation, given by

\[X_m(t) = X_m(0) + \sum_{k=1}^K \nu_m^k Y_k\left( \int_{0}^t a_k(\mathbf{X}(s^-),s) \, ds \right), \quad m = 1,\dots,M.\]

Here each $Y_k(t)$ denotes an independent unit rate Poisson counting process with $Y_k(0) = 0$, which counts the number of times the $k$th reaction has occurred up to time $t$. Realizations of these processes can be generated by creating JumpProblems from Catalyst ReactionSystems, and sampling the processes using the solvers in JumpProcesses.jl.

Let $P(\mathbf{x},t) = \operatorname{Prob}[\mathbf{X}(t) = \mathbf{x}]$ represent the probability the state of the system, $\mathbf{X}(t)$, has the concrete value $\mathbf{x}$ at time $t$. The forward equation, i.e. Chemical Master Equation (CME), associated with $\mathbf{X}(t)$ is then the coupled system of ODEs over all possible values for $\mathbf{x}$ given by

\[\frac{dP}{dt}(\mathbf{x},t) = \sum_{k=1}^k \left[ a_k(\mathbf{x} - \nu^k,t) P(\mathbf{x} - \nu^k,t) - a_k(\mathbf{x},t) P(\mathbf{x},t) \right].\]

While Catalyst does not currently support generating and solving for $P(\mathbf{x},t)$, for sufficiently small models the FiniteStateProjection.jl package can be used to generate and solve such models directly from Catalyst ReactionSystems.

Stochastic Chemical Kinetics Jump Process Example

Consider the same network as above,

rn = @reaction_network begin
    k₁, 2A + B --> 3C
    k₂, A --> 0
    k₃, 0 --> A
end

The time change process representation would be

\[\begin{align*} A(t) &= A(0) - 2 Y_1\left( \frac{k_1}{2} \int_0^t A(s^-)(A(s^-)-1) B(s^-) \, ds \right) - Y_2 \left( k_2 \int_0^t A(s^-) \, ds \right) + Y_3 \left( k_3 t \right) \\ B(t) &= B(0) - Y_1\left( \frac{k_1}{2} \int_0^t A(s^-)(A(s^-)-1) B(s^-) \, ds \right) \\ C(t) &= C(0) + 3 Y_1\left( \frac{k_1}{2} \int_0^t A(s^-)(A(s^-)-1) B(s^-) \, ds \right), \end{align*}\]

while the CME would be the coupled (infinite) system of ODEs over all realizable values of the non-negative integers $a$, $b$, and $c$ given by

\[\begin{align*} \frac{dP}{dt}(a,b,c,t) &= \left[\tfrac{k_1}{2} (a+2) (a+1) (b+1) P(a+2,b+1,c-3,t) - \tfrac{k_1}{2} a (a-1) b P(a,b,c,t)\right] \\ &\phantom{=} + \left[k_2 (a+1) P(a+1,b,c,t) - k_2 a P(a,b,c,t)\right] \\ &\phantom{=} + \left[k_3 P(a-1,b,c,t) - k_3 P(a,b,c,t)\right]. \end{align*}\]

If we initially have $A(0) = a_0$, $B(0) = b_0$, and $C(0) = c_0$ then we would have one ODE for each of possible state $(a,b,c)$ where $a \in \{0,1,\dots\}$ (i.e. $a$ can be any non-negative integer), $b \in \{0,1,\dots,b_0\}$, and $c \in \{c_0, c_0 + 1,\dots, c_0 + 3 b_0\}$. Other initial conditions would lead to different possible ranges for $a$, $b$, and $c$.