Working with conservation laws

Catalyst can detect, and eliminate for differential-equation based models, conserved quantities, i.e. linear combinations of species which are conserved via the chemistry. This functionality is both automatically utilised by Catalyst (e.g. to remove singular Jacobians during steady state computations), but is also available for users to utilise directly (e.g. to potentially improve simulation performance).

To illustrate conserved quantities, let us consider the following two-state model:

using Catalyst
rs = @reaction_network begin
 (k₁,k₂), X₁ <--> X₂
end

\[ \begin{align*} \mathrm{X_1} &\xrightleftharpoons[k_2]{k_1} \mathrm{X_2} \end{align*} \]

If we simulate it, we note that while the concentrations of $X₁$ and $X₂$ change throughout the simulation, the total concentration of $X$ ($= X₁ + X₂$) is constant:

using OrdinaryDiffEq, Plots
u0 = [:X₁ => 80.0, :X₂ => 20.0]
ps = [:k₁ => 10.0, :k₂ => 2.0]
oprob = ODEProblem(rs, u0, (0.0, 1.0), ps)
sol = solve(oprob)
plot(sol; idxs = [rs.X₁, rs.X₂, rs.X₁ + rs.X₂], label = ["X₁" "X₂" "X₁ + X₂ (a conserved quantity)"])

This makes sense, as while $X$ is converted between two different forms ($X₁$ and $X₂$), it is neither produced nor degraded. That is, it is a conserved quantity. Next, if we consider the ODE that our model generates:

using Latexify
latexify(rs; form = :ode)

\[\begin{align} \frac{\mathrm{d} X_1\left( t \right)}{\mathrm{d}t} &= - k_1 X_1\left( t \right) + k_2 X_2\left( t \right) \\ \frac{\mathrm{d} X_2\left( t \right)}{\mathrm{d}t} &= k_1 X_1\left( t \right) - k_2 X_2\left( t \right) \end{align} \]

we note that it essentially generates the same equation twice (i.e. $\frac{dX₁(t)}{dt} = -\frac{dX₂(t)}{dt}$). By designating our conserved quantity $X₁ + X₂ = Γ$, we can rewrite our differential equation model as a differential-algebraic equation (with a single differential equation and a single algebraic equation):

\[\frac{dX₁(t)}{dt} = - k₁X₁(t) + k₂(-X₁(t) + Γ) \\ X₂(t) = -X₁(t) + Γ\]

Using Catalyst, it is possible to detect any such conserved quantities and eliminate them from the system. Here, when we convert our ReactionSystem to an ODESystem, we provide the remove_conserved = true argument to instruct Catalyst to perform this elimination:

osys = convert(ODESystem, rs; remove_conserved = true)

\[ \begin{align} \frac{\mathrm{d} X_1\left( t \right)}{\mathrm{d}t} &= - k_1 X_1\left( t \right) + k_2 \left( - X_1\left( t \right) + \Gamma_{1} \right) \end{align} \]

We note that the output system only contains a single (differential) equation and can hence be solved with an ODE solver. The second (algebraic) equation is stored as an observable, and can be retrieved using the observed function:

observed(osys)

\[ \begin{align} X_2\left( t \right) &= - X_1\left( t \right) + \Gamma_{1} \end{align} \]

Using the unknowns function we can confirm that the ODE only has a single unknown variable:

unknowns(osys)

1-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 X₁(t)

Next, using parameters we note that an additional parameter, Γ[1] has been added to the system:

parameters(osys)

3-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 k₁
 k₂
 Γ[1]

Here, Catalyst encodes all conserved quantities in a single, vector-valued, parameter called Γ.

Practically, the remove_conserved = true argument can be provided when a ReactionSystem is converted to an ODEProblem:

using OrdinaryDiffEq, Plots
u0 = [:X₁ => 80.0, :X₂ => 20.0]
ps = [:k₁ => 10.0, :k₂ => 2.0]
oprob = ODEProblem(rs, u0, (0.0, 1.0), ps; remove_conserved = true)

┌ 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 ~/work/Catalyst.jl/Catalyst.jl/src/reactionsystem_conversions.jl:456

Here, while Γ[1] becomes a parameter of the new system, it has a default value equal to the corresponding conservation law. Hence, its value is computed from the initial condition [:X₁ => 80.0, :X₂ => 20.0], and does not need to be provided in the parameter vector. Next, we can simulate and plot our model using normal syntax:

sol = solve(oprob)
plot(sol)

While X₂ is an observable (and not unknown) of the ODE, we can access it just like if remove_conserved = true had not been used:

sol[:X₂]

14-element Vector{Float64}:
 20.0
 43.56620284189201
 56.39450336401671
 67.61060115947159
 74.28429177308591
 78.61346950718746
 81.00276209445914
 82.28947026842505
 82.90853333133784
 83.18215980724716
 83.28750769854895
 83.32210994972924
 83.3312120554613
 83.33290877375747

Conservation law accessor functions

For any given ReactionSystem model, we can use conservationlaw_constants to compute all of a system's conserved quantities:

conservationlaw_constants(rs)

\[ \begin{align} \Gamma_{1} &= X_1\left( t \right) + X_2\left( t \right) \end{align} \]

Next, the conservedequations function can be used to compute the algebraic equations produced when a system's conserved quantities are eliminated:

conservedequations(rs)

\[ \begin{align} X_2\left( t \right) &= - X_1\left( t \right) + \Gamma_{1} \end{align} \]

Finally, the conservationlaws function yields a $m \times n$ matrix, where $n$ is a system's number of species, $m$ its number of conservation laws, and element $(i,j)$ corresponds to the copy number of the $j$th species that occurs in the $i$th conserved quantity:

conservationlaws(rs)

1×2 Matrix{Int64}:
 1  1

I.e. in this case we have a single conserved quantity, which contains a single copy each of the system's two species.