|
| 1 | +# [Mathematical Models Catalyst can Generate](@id math_models_in_catalyst) |
| 2 | +We now describe the types of mathematical models that Catalyst can generate from |
| 3 | +chemical reaction networks (CRNs), corresponding to reaction rate equation (RRE) |
| 4 | +ordinary differential equation (ODE) models, Chemical Langevin equation (CLE) |
| 5 | +stochastic differential equation (SDE) models, and stochastic chemical kinetics |
| 6 | +(jump process) models. For each we show the abstract representations for the |
| 7 | +models that Catalyst can support, along with concrete examples. Note that we |
| 8 | +restrict ourselves to models involving only chemical reactions, and do not |
| 9 | +consider more general models that Catalyst can support such as coupling in |
| 10 | +non-reaction ODEs, algebraic equations, or events. Please see the broader |
| 11 | +documentation for more details on how Catalyst supports such functionality. |
| 12 | + |
| 13 | +!!! note |
| 14 | + This documentation assumes you have already read the [Introduction to Catalyst](@ref introduction_to_catalyst) tutorial. |
| 15 | + |
| 16 | +## General Chemical Reaction Notation |
| 17 | +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 |
| 18 | +```math |
| 19 | +\mathbf{X}(t) = \begin{pmatrix} X_1(t) \\ \vdots \\ X_M(t)) \end{pmatrix} |
| 20 | +``` |
| 21 | +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. |
| 22 | + |
| 23 | +The $k$th chemical reaction is given by |
| 24 | +```math |
| 25 | +\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 |
| 26 | +``` |
| 27 | +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)$. |
| 28 | + |
| 29 | +As explained in [the Catalyst introduction](@ref introduction_to_catalyst), for a mass action reaction where the preceding reaction has a fixed rate constant, $k$, this function would be the rate law |
| 30 | +```math |
| 31 | +a_k(\mathbf{X}(t)) = k \prod_{m=1}^M \frac{(X_m(t))^{\sigma_m^k}}{\sigma_m^k!}, |
| 32 | +``` |
| 33 | +for RRE ODE and CLE SDE models, and the propensity function |
| 34 | +```math |
| 35 | +a_k(\mathbf{X}(t)) = k \prod_{m=1}^M \frac{X_m(t) (X_m(t)-1) \dots (X_m(t)-\sigma_m^k+1)}{\sigma_m^k!}, |
| 36 | +``` |
| 37 | +for stochastic chemical kinetics jump process models. |
| 38 | + |
| 39 | +### Rate Law vs. Propensity Example: |
| 40 | +For the reaction $2A + B \overset{k}{\to} 3 C$ we would have |
| 41 | +```math |
| 42 | +\mathbf{X}(t) = (A(t), B(t), C(t)) |
| 43 | +``` |
| 44 | +with $\sigma_1 = 2$, $\sigma_2 = 1$, $\sigma_3 = 0$, $\beta_1 = 0$, $\beta_2 = |
| 45 | +0$, $\beta_3 = 3$, $\nu_1 = -2$, $\nu_2 = -1$, and $\nu_3 = 3$. For an ODE/SDE |
| 46 | +model we would have the rate law |
| 47 | +```math |
| 48 | +a(\mathbf{X}(t)) = \frac{k}{2} A^2 B |
| 49 | +``` |
| 50 | +while for a jump process model we would have the propensity function |
| 51 | +```math |
| 52 | +a(\mathbf{X}(t)) = \frac{k}{2} A (A-1) B. |
| 53 | +``` |
| 54 | + |
| 55 | +Note, if the combinatoric factors are already included in one's rate constants, |
| 56 | +the implicit rescaling of rate constants can be disabled through use of the |
| 57 | +`combinatoric_ratelaws = false` argument to [`Base.convert`](@ref) or whatever |
| 58 | +Problem is being generated, i.e. |
| 59 | +```julia |
| 60 | +rn = @reaction_network ... |
| 61 | +osys = convert(ODESystem, rn; combinatoric_ratelaws = false) |
| 62 | +oprob = ODEProblem(osys, ...) |
| 63 | + |
| 64 | +# or |
| 65 | +oprob = ODEProblem(rn, ...; combinatoric_ratelaws = false) |
| 66 | +``` |
| 67 | +In this case our ODE/SDE rate law would be |
| 68 | +```math |
| 69 | +a(\mathbf{X}(t)) = k A^2 B |
| 70 | +``` |
| 71 | +while the jump process propensity function is |
| 72 | +```math |
| 73 | +a(\mathbf{X}(t)) = k A (A-1) B. |
| 74 | +``` |
| 75 | + |
| 76 | +## Reaction Rate Equation (RRE) ODE Models |
| 77 | +The RRE ODE models Catalyst creates for a general system correspond to the coupled system of ODEs given by |
| 78 | +```math |
| 79 | +\frac{d X_m}{dt} =\sum_{k=1}^K \nu_m^k a_k(\mathbf{X}(t),t), \quad m = 1,\dots,M. |
| 80 | +``` |
| 81 | +These models can be generated by creating `ODEProblem`s from Catalyst `ReactionSystem`s, and solved using the solvers in [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl). Similarly, creating `NonlinearProblem`s or `SteadyStateProblem`s will generate the coupled algebraic system of steady-state equations associated with a RRE ODE model, i.e. |
| 82 | +```math |
| 83 | +0 =\sum_{k=1}^K \nu_m^k a_k(\bar{\mathbf{X}}), \quad m = 1,\dots,M |
| 84 | +``` |
| 85 | +for a steady-state $\bar{\mathbf{X}}$. Note, here we have assumed the rate laws are [autonomous](https://en.wikipedia.org/wiki/Autonomous_system_(mathematics)) so that the equations are well-defined. |
| 86 | + |
| 87 | +### RRE ODE Example |
| 88 | +Let's see the generated ODEs for the following network |
| 89 | +```@example math_examples |
| 90 | +using Catalyst, ModelingToolkit, Latexify |
| 91 | +rn = @reaction_network begin |
| 92 | + k₁, 2A + B --> 3C |
| 93 | + k₂, A --> 0 |
| 94 | + k₃, 0 --> A |
| 95 | +end |
| 96 | +osys = convert(ODESystem, rn) |
| 97 | +``` |
| 98 | +Likewise, the following drops the combinatoric scaling factors, giving unscaled ODEs |
| 99 | +```@example math_examples |
| 100 | +osys = convert(ODESystem, rn; combinatoric_ratelaws = false) |
| 101 | +``` |
| 102 | + |
| 103 | +## Chemical Langevin Equation (CLE) SDE Models |
| 104 | +The CLE SDE models Catalyst creates for a general system correspond to the coupled system of SDEs given by |
| 105 | +```math |
| 106 | +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, |
| 107 | +``` |
| 108 | +where each $W_k(t)$ represents an independent, standard Brownian Motion. Realizations of these processes can be generated by creating `SDEProblem`s from Catalyst `ReactionSystem`s, and sampling the processes using the solvers in [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl). |
| 109 | + |
| 110 | +### CLE SDE Example |
| 111 | +Consider the same network as above, |
| 112 | +```julia |
| 113 | +rn = @reaction_network begin |
| 114 | + k₁, 2A + B --> 3C |
| 115 | + k₂, A --> 0 |
| 116 | + k₃, 0 --> A |
| 117 | +end |
| 118 | +``` |
| 119 | +We obtain the CLE SDEs |
| 120 | +```math |
| 121 | +\begin{align} |
| 122 | +dA(t) &= \left(- k_1 A^{2} B - k_2 A + k_3 \right) dt |
| 123 | + - 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) |
| 124 | +\\ |
| 125 | +dB(t) &= - \frac{k_1}{2} A^{2} B \, dt - \sqrt{\frac{k_1}{2} A^{2} B} \, dW_1(t) \\ |
| 126 | +dC(t) &= \frac{3}{2} k_1 A^{2} B \, dt + 3 \sqrt{\frac{k_1}{2} A^{2} B} \, dW_1(t). |
| 127 | +\end{align} |
| 128 | +``` |
| 129 | + |
| 130 | +## Stochastic Chemical Kinetics Jump Process Models |
| 131 | +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 |
| 132 | +```math |
| 133 | +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. |
| 134 | +``` |
| 135 | +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 `JumpProblem`s from Catalyst `ReactionSystem`s, and sampling the processes using the solvers in [JumpProcesses.jl](https://github.com/SciML/JumpProcesses.jl). |
| 136 | + |
| 137 | +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 |
| 138 | +```math |
| 139 | +\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]. |
| 140 | +``` |
| 141 | +While Catalyst does not currently support generating and solving for $P(\mathbf{x},t)$, for sufficiently small models the [FiniteStateProjection.jl](https://github.com/SciML/FiniteStateProjection.jl) package can be used to generate and solve such models directly from Catalyst [`ReactionSystem`](@ref)s. |
| 142 | + |
| 143 | +### Stochastic Chemical Kinetics Jump Process Example |
| 144 | +Consider the same network as above, |
| 145 | +```julia |
| 146 | +rn = @reaction_network begin |
| 147 | + k₁, 2A + B --> 3C |
| 148 | + k₂, A --> 0 |
| 149 | + k₃, 0 --> A |
| 150 | +end |
| 151 | +``` |
| 152 | +The time change process representation would be |
| 153 | +```math |
| 154 | +\begin{align*} |
| 155 | +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) \\ |
| 156 | +B(t) &= B(0) - Y_1\left( \frac{k_1}{2} \int_0^t A(s^-)(A(s^-)-1) B(s^-) \, ds \right) \\ |
| 157 | +C(t) &= C(0) + 3 Y_1\left( \frac{k_1}{2} \int_0^t A(s^-)(A(s^-)-1) B(s^-) \, ds \right), |
| 158 | +\end{align*} |
| 159 | +``` |
| 160 | +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 |
| 161 | +```math |
| 162 | +\begin{align*} |
| 163 | +\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] \\ |
| 164 | +&\phantom{=} + \left[k_2 (a+1) P(a+1,b,c,t) - k_2 a P(a,b,c,t)\right] \\ |
| 165 | +&\phantom{=} + \left[k_3 P(a-1,b,c,t) - k_3 P(a,b,c,t)\right]. |
| 166 | +\end{align*} |
| 167 | +``` |
| 168 | +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$. |
0 commit comments