Merge pull request #987 from isaacsas/add_math_docs

math model description page

HISTORY.md

@@ -39,9 +39,8 @@ briefly summarised in the following bullet points:
   new merged and/or composed `ReactionSystem`s from multiple component systems.
 
 #### General changes
-- The `default_t()` and `default_time_deriv()` functions are now the preferred
-  approaches for creating the default time independent variable and its
-  differential. i.e.
+- `default_t()` and `default_time_deriv()` functions should be used for creating
+  the default time independent variable and its differential. i.e.
   ```julia
   # do
   t = default_t()

docs/pages.jl

@@ -2,7 +2,8 @@ pages = Any[
     "Home" => "index.md",
     "Introduction to Catalyst" => Any[
         "introduction_to_catalyst/catalyst_for_new_julia_users.md",
-        "introduction_to_catalyst/introduction_to_catalyst.md"
+        "introduction_to_catalyst/introduction_to_catalyst.md",
+        "introduction_to_catalyst/math_models_intro.md"
     ],
     "Model Creation and Properties" => Any[
         "model_creation/dsl_basics.md",

docs/src/introduction_to_catalyst/introduction_to_catalyst.md

@@ -3,7 +3,11 @@ In this tutorial we provide an introduction to using Catalyst to specify
 chemical reaction networks, and then to solve ODE, jump, and SDE models
 generated from them [1]. At the end we show what mathematical rate laws and
 transition rate functions (i.e. intensities or propensities) are generated by
-Catalyst for ODE, SDE and jump process models.
+Catalyst for ODE, SDE and jump process models. The [Mathematical Models Catalyst
+can Generate](@ref math_models_in_catalyst) documentation illustrates the
+abstract mathematical models Catalyst reaction models can be converted to,
+but please note it assumes one has already read this tutorial as a
+prerequisite.
 
 We begin by installing Catalyst and any needed packages into a new environment.
 This step can be skipped if you have already installed them in your current,
@@ -350,6 +354,8 @@ and the ODE model
 \end{align*}
 ```
 
+A more detailed summary of the precise mathematical equations Catalyst can generate is available in the [Mathematical Models Catalyst can Generate](@ref math_models_in_catalyst) documentation.
+
 ---
 ## Notes
 1. For each of the preceding models we converted the `ReactionSystem` to, i.e.,

docs/src/introduction_to_catalyst/math_models_intro.md

@@ -0,0 +1,168 @@
+# [Mathematical Models Catalyst can Generate](@id math_models_in_catalyst)
+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.
+
+!!! note
+    This documentation assumes you have already read the [Introduction to Catalyst](@ref introduction_to_catalyst) tutorial.
+
+## 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
+```math
+\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
+```math
+\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](@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
+```math
+a_k(\mathbf{X}(t)) = k \prod_{m=1}^M \frac{(X_m(t))^{\sigma_m^k}}{\sigma_m^k!},
+```
+for RRE ODE and CLE SDE models, and the propensity function
+```math
+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!},
+```
+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
+```math
+\mathbf{X}(t) = (A(t), B(t), C(t))
+```
+with $\sigma_1 = 2$, $\sigma_2 = 1$, $\sigma_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
+```math
+a(\mathbf{X}(t)) = \frac{k}{2} A^2 B
+```
+while for a jump process model we would have the propensity function
+```math
+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`](@ref) or whatever
+Problem is being generated, i.e.
+```julia
+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
+```math
+a(\mathbf{X}(t)) = k A^2 B
+```
+while the jump process propensity function is
+```math
+a(\mathbf{X}(t)) = k A (A-1) B.
+```
+
+## 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
+```math
+\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 `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.
+```math
+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](https://en.wikipedia.org/wiki/Autonomous_system_(mathematics)) so that the equations are well-defined.
+
+### RRE ODE Example
+Let's see the generated ODEs for the following network
+```@example math_examples
+using Catalyst, ModelingToolkit, Latexify
+rn = @reaction_network begin
+    k₁, 2A + B --> 3C
+    k₂, A --> 0
+    k₃, 0 --> A
+end
+osys = convert(ODESystem, rn)
+```
+Likewise, the following drops the combinatoric scaling factors, giving unscaled ODEs
+```@example math_examples
+osys = convert(ODESystem, rn; combinatoric_ratelaws = false)
+```
+
+## 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
+```math
+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 `SDEProblem`s from Catalyst `ReactionSystem`s, and sampling the processes using the solvers in [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl).
+
+### CLE SDE Example
+Consider the same network as above,
+```julia
+rn = @reaction_network begin
+    k₁, 2A + B --> 3C
+    k₂, A --> 0
+    k₃, 0 --> A
+end
+```
+We obtain the CLE SDEs
+```math
+\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
+```math
+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 `JumpProblem`s from Catalyst `ReactionSystem`s, and sampling the processes using the solvers in [JumpProcesses.jl](https://github.com/SciML/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
+```math
+\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](https://github.com/SciML/FiniteStateProjection.jl) package can be used to generate and solve such models directly from Catalyst [`ReactionSystem`](@ref)s.
+
+### Stochastic Chemical Kinetics Jump Process Example
+Consider the same network as above,
+```julia
+rn = @reaction_network begin
+    k₁, 2A + B --> 3C
+    k₂, A --> 0
+    k₃, 0 --> A
+end
+```
+The time change process representation would be
+```math
+\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
+```math
+\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$.

docs/src/v14_migration_guide.md

@@ -136,18 +136,19 @@ using Catalyst
 @variables t
 nothing # hide
 ```
-MTKv9 has introduced a standard global time variable, and as such a new, preferred, interface has been developed:
+MTKv9 has introduced a standard global time variable, and as such a new interface has been developed:
 ```@example v14_migration_3
 t = default_t()
 nothing # hide
 ```
+Note that internally MTK9 actually represents `t` as a parameter, so the old approach should never be used. As the type of `t` is now considered internal to MTK, it should always be declared using the `default_t()` function.
 
 Similarly, the time differential (primarily relevant when creating combined reaction-ODE models) used to be declared through
 ```@example v14_migration_3
 D = Differential(t)
 nothing # hide
 ```
-where the preferred method is now
+where now one must use
 ```@example v14_migration_3
 D = default_time_deriv()
 nothing # hide