Specialised conservation law section on the network analysis doc page

README.md

@@ -55,7 +55,7 @@ be found in its corresponding research paper, [Catalyst: Fast and flexible model
 - Leveraging ModelingToolkit, generated models can be converted to symbolic reaction rate equation ODE models, symbolic Chemical Langevin Equation models, and symbolic stochastic chemical kinetics (jump process) models. These can be simulated using any [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) [ODE/SDE/jump solver](https://docs.sciml.ai/Catalyst/stable/model_simulation/simulation_introduction/), and can be used within `EnsembleProblem`s for carrying out [parallelized parameter sweeps and statistical sampling](https://docs.sciml.ai/Catalyst/stable/model_simulation/ensemble_simulations/). Plot recipes are available for [visualization of all solutions](https://docs.sciml.ai/Catalyst/stable/model_simulation/simulation_plotting/).
 - Non-integer (e.g. `Float64`) stoichiometric coefficients [are supported](https://docs.sciml.ai/Catalyst/stable/model_creation/dsl_basics/#dsl_description_stoichiometries_decimal) for generating ODE models, and symbolic expressions for stoichiometric coefficients [are supported](https://docs.sciml.ai/Catalyst/stable/model_creation/parametric_stoichiometry/) for all system types.
 - A [network analysis suite](https://docs.sciml.ai/Catalyst/stable/model_creation/network_analysis/) permits the computation of linkage classes, deficiencies, reversibility, and other network properties.
-- [Conservation laws can be detected and utilized](https://docs.sciml.ai/Catalyst/stable/model_creation/network_analysis/#network_analysis_deficiency) to reduce system sizes, and to generate non-singular Jacobians (e.g. during conversion to ODEs, SDEs, and steady state equations).
+- [Conservation laws can be detected and utilized](https://docs.sciml.ai/Catalyst/stable/model_creation/network_analysis/#working_with_conservation_laws) to reduce system sizes, and to generate non-singular Jacobians (e.g. during conversion to ODEs, SDEs, and steady state equations).
 - Catalyst reaction network models can be [coupled with differential and algebraic equations](https://docs.sciml.ai/Catalyst/stable/model_creation/constraint_equations/) (which are then incorporated during conversion to ODEs, SDEs, and steady state equations).
 - Models can be [coupled with events](https://docs.sciml.ai/Catalyst/stable/model_creation/constraint_equations/#constraint_equations_events) that affect the system and its state during simulations.
 - By leveraging ModelingToolkit, users have a variety of options for generating optimized system representations to use in solvers. These include construction of [dense or sparse Jacobians](https://docs.sciml.ai/Catalyst/stable/model_simulation/ode_simulation_performance/#ode_simulation_performance_sparse_jacobian), [multithreading or parallelization of generated derivative functions](https://docs.sciml.ai/Catalyst/stable/model_simulation/ode_simulation_performance/#ode_simulation_performance_parallelisation), [automatic classification of reactions into optimized jump types for Gillespie type simulations](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#jump_types), [automatic construction of dependency graphs for jump systems](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#Jump-Aggregators-Requiring-Dependency-Graphs), and more.

docs/src/index.md

@@ -24,7 +24,7 @@ etc).
 - Leveraging ModelingToolkit, generated models can be converted to symbolic reaction rate equation ODE models, symbolic Chemical Langevin Equation models, and symbolic stochastic chemical kinetics (jump process) models. These can be simulated using any [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) [ODE/SDE/jump solver](@ref simulation_intro), and can be used within `EnsembleProblem`s for carrying out [parallelized parameter sweeps and statistical sampling](@ref ensemble_simulations). Plot recipes are available for [visualization of all solutions](@ref simulation_plotting).
 - Non-integer (e.g. `Float64`) stoichiometric coefficients [are supported](@ref dsl_description_stoichiometries_decimal) for generating ODE models, and symbolic expressions for stoichiometric coefficients [are supported](@ref parametric_stoichiometry) for all system types.
 - A [network analysis suite](@ref network_analysis) permits the computation of linkage classes, deficiencies, reversibility, and other network properties.
-- [Conservation laws can be detected and utilized](@ref network_analysis_deficiency) to reduce system sizes, and to generate non-singular Jacobians (e.g. during conversion to ODEs, SDEs, and steady state equations).
+- [Conservation laws can be detected and utilized](@ref network_analysis_conservation_laws) to reduce system sizes, and to generate non-singular Jacobians (e.g. during conversion to ODEs, SDEs, and steady state equations).
 - Catalyst reaction network models can be [coupled with differential and algebraic equations](@ref constraint_equations_coupling_constraints) (which are then incorporated during conversion to ODEs, SDEs, and steady state equations).
 - Models can be [coupled with events](@ref constraint_equations_events) that affect the system and its state during simulations.
 - By leveraging ModelingToolkit, users have a variety of options for generating optimized system representations to use in solvers. These include construction of [dense or sparse Jacobians](@ref ode_simulation_performance_sparse_jacobian), [multithreading or parallelization of generated derivative functions](@ref ode_simulation_performance_parallelisation), [automatic classification of reactions into optimized jump types for Gillespie type simulations](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#jump_types), [automatic construction of dependency graphs for jump systems](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#Jump-Aggregators-Requiring-Dependency-Graphs), and more.

docs/src/inverse_problems/structural_identifiability.md

@@ -124,7 +124,7 @@ rs = @reaction_network begin
     (k1,k2), X1 <--> X2
 end
 ```
-contain [conservation laws](@ref network_analysis_deficiency) (in this case $Γ = X1 + X2$, where $Γ = X1(0) + X2(0)$ is a constant). Because the presence of such conservation laws makes structural identifiability analysis prohibitively computationally expensive (for all but the simplest of cases), these are automatically eliminated by Catalyst. This is handled internally, and should not be noticeable to the user. The exception is the `make_si_ode` function. For each conservation law, its output will have one ODE removed, and instead have a conservation parameter (of the form `Γ[i]`) added to its equations. This feature can be disabled through the `remove_conserved = false` option.
+contain [conservation laws](@ref network_analysis_conservation_laws) (in this case $Γ = X1 + X2$, where $Γ = X1(0) + X2(0)$ is a constant). Because the presence of such conservation laws makes structural identifiability analysis prohibitively computationally expensive (for all but the simplest of cases), these are automatically eliminated by Catalyst. This is handled internally, and should not be noticeable to the user. The exception is the `make_si_ode` function. For each conservation law, its output will have one ODE removed, and instead have a conservation parameter (of the form `Γ[i]`) added to its equations. This feature can be disabled through the `remove_conserved = false` option.
 
 ## [Systems with exponent parameters](@id structural_identifiability_exp_params)
 Structural identifiability cannot currently be applied to systems with parameters (or species) in exponents. E.g. this

docs/src/model_creation/examples/basic_CRN_library.md

@@ -54,12 +54,12 @@ The two-state model describes a component (here called $X$) which can exist in t
 ```@example crn_library_two_states
 using Catalyst, OrdinaryDiffEq, StochasticDiffEq, Plots
 two_state_model = @reaction_network begin
-    (k1,k2), X₁ <--> X₂
+    (k₁,k₂), X₁ <--> X₂
 end
 
 u0 = [:X₁ => 50.0, :X₂ => 50.0]
 tspan = (0.0, 1.0)
-ps = [:k1 => 2.0, :k2 => 3.0]
+ps = [:k₁ => 2.0, :k₂ => 3.0]
 
 oprob = ODEProblem(two_state_model, u0, tspan, ps)
 osol = solve(oprob)

docs/src/model_creation/network_analysis.md

@@ -8,7 +8,7 @@ section
 Note, currently API functions for network analysis and conservation law analysis
 do not work with constant species (currently only generated by SBMLToolkit).
 
-## Network representation of the Repressilator `ReactionSystem`
+## [Network representation of the Repressilator `ReactionSystem`](@id network_analysis_repressilator_representation)
 We first load Catalyst and construct our model of the repressilator
 ```@example s1
 using Catalyst
@@ -49,7 +49,7 @@ the reaction rate equation ODE model for the repressilator is
 \end{aligned}
 ```
 
-## Matrix-vector reaction rate equation representation
+## [Matrix-vector reaction rate equation representation](@id network_analysis_matrix_vector_representation)
 In general, reaction rate equation (RRE) ODE models for chemical reaction networks can
 be represented as a first-order system of ODEs in a compact matrix-vector notation. Suppose
 we have a reaction network with ``K`` reactions and ``M`` species, labelled by the state vector
@@ -208,7 +208,7 @@ substrate complexes into product complexes where the rate is just a number or
 parameter, and red arrows indicate the conversion of substrate complexes into
 product complexes where the rate is an expression involving chemical species.
 
-## Aspects of reaction network structure
+## [Aspects of reaction network structure](@id network_analysis_structural_aspects)
 The reaction complex representation can be exploited via [Chemical Reaction
 Network Theory](https://en.wikipedia.org/wiki/Chemical_reaction_network_theory)
 to provide insight into possible steady state and time-dependent properties of
@@ -234,7 +234,7 @@ complexgraph(rn)
 
 ![network_1](../assets/complex_rn.svg)
 
-#### Linkage classes and sub-networks of the reaction network
+#### [Linkage classes and sub-networks of the reaction network](@id network_analysis_structural_aspects_linkage)
 The preceding reaction complex graph shows that `rn` is composed of two
 disconnected sub-graphs, one containing the complexes ``A+B``, ``C``, ``D+E``, and
 ``F``, the other containing the complexes ``2A``, ``B + G``, and ``H``. These sets,
@@ -273,7 +273,7 @@ and,
 
 ![subnetwork_2](../assets/complex_subnets2.svg)
 
-#### [Deficiency of the network](@id network_analysis_deficiency)
+#### [Deficiency of the network](@id network_analysis_structural_aspects_deficiency)
 A famous theorem in Chemical Reaction Network Theory, the Deficiency Zero
 Theorem [^1], allows us to use knowledge of the net stoichiometry matrix and the
 linkage classes of a *mass action* RRE ODE system to draw conclusions about the
@@ -325,7 +325,7 @@ Quoting Feinberg [^1]
 > stoichiometry vectors] are as independent as the partition of complexes into
 > linkage classes will allow.
 
-#### Reversibility of the network
+#### [Reversibility of the network](@id network_analysis_structural_aspects_reversibility)
 A reaction network is *reversible* if the "arrows" of the reactions are
 symmetric so that every reaction is accompanied by its reverse reaction.
 Catalyst's API provides the [`isreversible`](@ref) function to determine whether
@@ -374,7 +374,7 @@ isweaklyreversible(rn, subnets)
 Every reversible network is also weakly reversible, but not every weakly
 reversible network is reversible.
 
-#### Deficiency Zero Theorem
+#### [Deficiency Zero Theorem](@id network_analysis_structural_aspects_deficiency_zero_theorem)
 Knowing the deficiency and weak reversibility of a mass action chemical reaction
 network ODE model allows us to make inferences about the corresponding
 steady state behavior. Before illustrating how this works for one example, we
@@ -464,7 +464,7 @@ which we see is mass action and has deficiency zero, but is not weakly
 reversible. As such, we can conclude that for any choice of rate constants the
 RRE ODEs cannot have a positive equilibrium solution.
 
-## Caching of Network Properties in `ReactionSystems`
+## [Caching of Network Properties in `ReactionSystems`](@id network_analysis_caching_properties)
 When calling many of the network API functions, Catalyst calculates and caches
 in `rn` a variety of information. For example the first call to
 ```julia
@@ -500,6 +500,97 @@ Catalyst.reset_networkproperties!(rn)
 Network property functions will then recalculate their associated properties and
 cache the new values the next time they are called.
 
+
+## [Working with conservation laws](@id network_analysis_conservation_laws)
+Catalyst contain specialised functionality to work with *conserved quantities*, i.e. species quantities which are conserved throughout e.g. a simulation. This functionality is both automatically utilised by Catalyst (e.g. to [remove singular Jacobians during steady state computations](@ref homotopy_continuation_conservation_laws)), but also something that is available for the user to utilise directly (e.g. to [improve simulation performance](@ref ode_simulation_performance_conservation_laws)).
+
+To illustrate conserved quantities, let us consider the following [two-state](@ref basic_CRN_library_two_states) model:
+```@example network_analysis_conservation_laws
+using Catalyst
+rs = @reaction_network begin
+ (k₁,k₂), X₁ <--> X₂
+end
+```
+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:
+```@example network_analysis_conservation_laws
+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:
+```@example network_analysis_conservation_laws
+using Latexify
+latexify(rs; form = :ode)
+```
+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](https://en.wikipedia.org/wiki/Differential-algebraic_system_of_equations) (with a single differential equation and a single algebraic equation):
+```math
+\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:
+```@example network_analysis_conservation_laws
+osys = convert(ODESystem, rs; remove_conserved = true)
+```
+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*](@ref dsl_advanced_options_observables), and can be retrieved using the `observed` function:
+```@example network_analysis_conservation_laws
+observed(osys)
+```
+Using the `unknowns` function we can confirm that the ODE only has a single unknown variable:
+```@example network_analysis_conservation_laws
+unknowns(osys)
+```
+Next, using `parameters` we note that an additional parameter, `Γ[1]` has been added to the system:
+```@example network_analysis_conservation_laws
+parameters(osys)
+```
+Here, Catalyst encodes all conserved quantities in a single, [vector-valued](@ref dsl_advanced_options_vector_variables), parameter called `Γ`.
+
+Practically, the `remove_conserved = true` argument can be provided when a `ReactionSystem` is converted to an `ODEProblem`:
+```@example network_analysis_conservation_laws
+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)
+nothing # hide
+```
+Here, while `Γ[1]` becomes a parameter of the new system, it has a [default value](@ref dsl_advanced_options_default_vals) 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:
+```@example network_analysis_conservation_laws
+sol = solve(oprob)
+plot(sol)
+```
+!!! note
+    Any species eliminated using `remove_conserved = true` will not be plotted by default. However, it can be added to the plot using [the `idxs` plotting option](@ref simulation_plotting_options). E.g. here we would use `plot(sol; idxs = [:X₁, :X₂])` to plot both species.
+
+!!! danger
+    Currently, there is a bug in MTK where the values associated with conservation laws are not updated properly in response to [`remake`](@ref simulation_structure_interfacing_problems_remake) (or [other problem-updating functions, such as `getu`](@ref simulation_structure_interfacing_functions)). Hence, problems created using `remove_conserved = true` *should not* be modified, i.e. by changing initial conditions or the values of the `Γ`'s directly. Instead, to change these values new problems must be generated with the new initial condition map.
+
+While `X₂` is an observable (and not unknown) of the ODE, we can [access it](@ref simulation_structure_interfacing_problems) just like if `remove_conserved = true` had not been used:
+```@example network_analysis_conservation_laws
+sol[:X₂]
+```
+!!! note
+    Generally, `remove_conserved = true` should not change any model workflows. I.e. anything that works without this option should also work when an `ODEProblem` is created using `remove_conserved = true`.
+
+For any given `ReactionSystem` model, we can use `conservationlaw_constants` to compute all of a system's conserved quantities:
+```@example network_analysis_conservation_laws
+conservationlaw_constants(rs)
+```
+Next, the `conservedequations` function can be used to compute the algebraic equations produced when a system's conserved quantities are eliminated:
+```@example network_analysis_conservation_laws
+conservedequations(rs)
+```
+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:
+```@example network_analysis_conservation_laws
+conservationlaws(rs)
+```
+I.e. in this case we have a single conserved quantity, which contains a single copy each of the system's two species.
+
+!!! note
+    The `remove_conserved = true` option is available when creating `SDEProblem`s, `NonlinearProblem`s, and `SteadyStateProblem`s (and their corresponding systems). However, it cannot be used when creating `JumpProblem`s.
+
 ---
 ## References
 [^1]: [Feinberg, M. *Foundations of Chemical Reaction Network Theory*, Applied Mathematical Sciences 202, Springer (2019).](https://link.springer.com/book/10.1007/978-3-030-03858-8?noAccess=true)