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 their 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{\mathtt{X_1}} &\xrightleftharpoons[\mathtt{k_2}]{\mathtt{k_1}} \mathrm{\mathtt{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 OrdinaryDiffEqDefault, 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} \mathtt{X_1}\left( t \right)}{\mathrm{d}t} &= - \mathtt{k_1} \mathtt{X_1}\left( t \right) + \mathtt{k_2} \mathtt{X_2}\left( t \right) \\ \frac{\mathrm{d} \mathtt{X_2}\left( t \right)}{\mathrm{d}t} &= \mathtt{k_1} \mathtt{X_1}\left( t \right) - \mathtt{k_2} \mathtt{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} \mathtt{X_1}\left( t \right)}{\mathrm{d}t} &= - \mathtt{k_1} \mathtt{X_1}\left( t \right) + \mathtt{k_2} \left( - \mathtt{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} \mathtt{X_2}\left( t \right) &= - \mathtt{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 (vector) parameter, Γ has been added to the system:

parameters(osys)

3-element Vector{SymbolicUtils.BasicSymbolic}:
 k₁
 k₂
 Γ

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 OrdinaryDiffEqDefault, 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)

Here, while Γ 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.39450301544777
 67.61060152222583
 74.28429170304489
 78.61346954164246
 81.00276197209632
 82.28947021914873
 82.908533284148
 83.18215977995723
 83.28750769219297
 83.32210994898414
 83.331212055774
 83.33290877374532

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} &= \mathtt{X_1}\left( t \right) + \mathtt{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} \mathtt{X_2}\left( t \right) &= - \mathtt{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.

Updating conservation law values directly

Previously we noted that conservation law elimination adds the conservation law as a system parameter (named $Γ$). E.g. here we find it among the parameters of the ODE model generated by the previously used two-state system:

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

3-element Vector{SymbolicUtils.BasicSymbolic}:
 k₁
 k₂
 Γ

An advantage of this is that we can set conserved quantities's values directly in simulations. Here we simulate the model while specifying that $X₁ + X₂ = 10.0$ (and while also specifying an initial condition for $X₁$, but none for $X₂$).

u0 = [:X₁ => 6.0]
ps = [:k₁ => 1.0, :k₂ => 2.0, :Γ => [10.0]]
oprob = ODEProblem(rs, u0, 1.0, ps; remove_conserved = true)
sol = solve(oprob)

Generally, for each conservation law, one can omit specifying either the conservation law constant, or one initial condition (whichever quantity is missing will be computed from the other ones).

Just like when we create a problem, if we update the species (or conservation law parameter) values of oprob, the remaining ones will be recomputed to generate an accurate conservation law. E.g. here we create an ODEProblem, check the value of the conservation law, and then confirm that its value is updated with $X₁$.

u0 = [:X₁ => 6.0, :X₂ => 4.0]
ps = [:k₁ => 1.0, :k₂ => 2.0]
oprob = ODEProblem(rs, u0, 10.0, ps; remove_conserved = true)
oprob.ps[:Γ][1]

10.0

oprob = remake(oprob; u0 = [:X₁ => 16.0])
oprob.ps[:Γ][1]

20.0

It is also possible to update the value of $Γ$. Here, as $X₂$ is the species eliminated by the conservation law (which can be checked using conservedequations), $X₂$'s value will be modified to ensure that $Γ$'s new value is correct:

oprob = remake(oprob; p = [:Γ => [30.0]])
oprob[:X₂]

14.0

Generally, for a conservation law where we have:

• One (or more) non-eliminated species.
• One eliminated species.
• A conservation law parameter.

If the value of the conservation law parameter $Γ$'s value has never been specified, its value will be updated to accommodate changes in species values. If the conservation law parameter ($Γ$)'s value has been specified (either when the ODEProblem was created, or using remake), then the eliminated species's value will be updated to accommodate changes in the conservation law parameter or non-eliminated species's values. Furthermore, in this case, the value of the eliminated species cannot be updated.

Extracting the conservation law parameter's symbolic variable

Catalyst represents its models using the Symbolics.jl computer algebraic system, something which allows the user to form symbolic expressions of model components. If you wish to extract and use the symbolic variable corresponding to a model's conserved quantities, you can use the following syntax:

Γ = Catalyst.get_networkproperties(rs).conservedconst

\[ \begin{equation} \Gamma \end{equation} \]