diff --git a/README.md b/README.md index 1eb5ddc5d2..93d97be52c 100644 --- a/README.md +++ b/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/#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). +- [Conservation laws can be detected and utilized](https://docs.sciml.ai/Catalyst/stable/model_creation/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. diff --git a/docs/pages.jl b/docs/pages.jl index d12215f4d1..7f6181a11d 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -11,6 +11,7 @@ pages = Any[ "model_creation/programmatic_CRN_construction.md", "model_creation/compositional_modeling.md", "model_creation/constraint_equations.md", + "model_creation/conservation_laws.md", "model_creation/parametric_stoichiometry.md", "model_creation/model_file_loading_and_export.md", "model_creation/model_visualisation.md", diff --git a/docs/src/index.md b/docs/src/index.md index 5781e5b5d0..06d93dca9f 100644 --- a/docs/src/index.md +++ b/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_conservation_laws) 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 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. diff --git a/docs/src/inverse_problems/structural_identifiability.md b/docs/src/inverse_problems/structural_identifiability.md index e723673a7e..649f68ab40 100644 --- a/docs/src/inverse_problems/structural_identifiability.md +++ b/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_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. +contain [conservation laws](@ref 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 diff --git a/docs/src/model_creation/conservation_laws.md b/docs/src/model_creation/conservation_laws.md new file mode 100644 index 0000000000..383b85812c --- /dev/null +++ b/docs/src/model_creation/conservation_laws.md @@ -0,0 +1,91 @@ +# [Working with conservation laws](@id 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](@ref homotopy_continuation_conservation_laws)), but is also available for users to utilise directly (e.g. to potentially [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 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 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 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 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 conservation_laws +observed(osys) +``` +Using the `unknowns` function we can confirm that the ODE only has a single unknown variable: +```@example conservation_laws +unknowns(osys) +``` +Next, using `parameters` we note that an additional parameter, `Γ[1]` has been added to the system: +```@example 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 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 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 parameters, `Γ`, 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 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`. + +!!! 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. + +## [Conservation law accessor functions](@id conservation_laws_accessors) + +For any given `ReactionSystem` model, we can use `conservationlaw_constants` to compute all of a system's conserved quantities: +```@example 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 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 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. \ No newline at end of file diff --git a/docs/src/model_creation/network_analysis.md b/docs/src/model_creation/network_analysis.md index 138cb1a799..3f77a99f55 100644 --- a/docs/src/model_creation/network_analysis.md +++ b/docs/src/model_creation/network_analysis.md @@ -306,7 +306,9 @@ we find that there are five independent species. Let's check this is correct: using LinearAlgebra rank(netstoichmat(rn)) == 5 ``` -So we know that the rank of our reaction network is five. +So we know that the rank of our reaction network is five. An extended section discussing how to +utilise conservation law elimination during chemical reaction network modelling can be found +[here](@ref conservation_laws). The deficiency, ``\delta``, of a reaction network is defined as ```math @@ -500,97 +502,6 @@ 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) diff --git a/docs/src/model_simulation/ode_simulation_performance.md b/docs/src/model_simulation/ode_simulation_performance.md index 0b6d70e314..11934526ff 100644 --- a/docs/src/model_simulation/ode_simulation_performance.md +++ b/docs/src/model_simulation/ode_simulation_performance.md @@ -172,7 +172,7 @@ Generally, the use of preconditioners is only recommended for advanced users who ## [Elimination of system conservation laws](@id ode_simulation_performance_conservation_laws) -Previously, we have described how Catalyst, when it generates ODEs, is able to [detect and eliminate conserved quantities](@ref network_analysis_conservation_laws). In certain cases, doing this can improve performance. E.g. in the following example we will eliminate the single conserved quantity in a [two-state model](@ref basic_CRN_library_two_states). This results in a differential algebraic equation with a single differential equation and a single algebraic equation (as opposed to two differential equations). However, as the algebraic equation is fully determined by the ODE solution, Catalyst moves it to be an observable and our new system therefore only contains one ODE that must be solved numerically. Conservation laws can be eliminated by providing the `remove_conserved = true` option to `ODEProblem`: +Previously, we have described how Catalyst, when it generates ODEs, is able to [detect and eliminate conserved quantities](@ref conservation_laws). In certain cases, doing this can improve performance. E.g. in the following example we will eliminate the single conserved quantity in a [two-state model](@ref basic_CRN_library_two_states). This results in a differential algebraic equation with a single differential equation and a single algebraic equation (as opposed to two differential equations). However, as the algebraic equation is fully determined by the ODE solution, Catalyst moves it to be an observable and our new system therefore only contains one ODE that must be solved numerically. Conservation laws can be eliminated by providing the `remove_conserved = true` option to `ODEProblem`: ```@example ode_simulation_performance_conservation_laws using Catalyst, OrdinaryDiffEq diff --git a/docs/src/spatial_modelling/lattice_reaction_systems.md b/docs/src/spatial_modelling/lattice_reaction_systems.md index 675dad48cf..32df70adcb 100644 --- a/docs/src/spatial_modelling/lattice_reaction_systems.md +++ b/docs/src/spatial_modelling/lattice_reaction_systems.md @@ -239,6 +239,6 @@ edge_parameters(lrs) ``` ## [Spatial modelling limitations](@id spatial_lattice_modelling_intro_limitations) -Many features which are supported for non-spatial `ReactionSystem`s are currently unsupported for [`LatticeReactionSystem`](@ref)s. This includes [observables](@ref dsl_advanced_options_observables), [algebraic and differential equations](@ref constraint_equations), [hierarchical models](@ref compositional_modeling), and [events](@ref constraint_equations_events). It is possible that these features will be supported in the future. Furthermore, [removal of conserved quantities](@ref network_analysis_conservation_laws) is not supported when creating spatial `ODEProblem`s. +Many features which are supported for non-spatial `ReactionSystem`s are currently unsupported for [`LatticeReactionSystem`](@ref)s. This includes [observables](@ref dsl_advanced_options_observables), [algebraic and differential equations](@ref constraint_equations), [hierarchical models](@ref compositional_modeling), and [events](@ref constraint_equations_events). It is possible that these features will be supported in the future. Furthermore, [removal of conserved quantities](@ref conservation_laws) is not supported when creating spatial `ODEProblem`s. If you are using Catalyst's features for spatial modelling, please give us feedback on how we can improve these features. Additionally, just letting us know that you use these features is useful, as it helps inform us whether continued development of spatial modelling features is worthwhile. diff --git a/docs/src/steady_state_functionality/homotopy_continuation.md b/docs/src/steady_state_functionality/homotopy_continuation.md index d1589d7903..85b7274cea 100644 --- a/docs/src/steady_state_functionality/homotopy_continuation.md +++ b/docs/src/steady_state_functionality/homotopy_continuation.md @@ -47,7 +47,7 @@ two_state_model = @reaction_network begin (k1,k2), X1 <--> X2 end ``` -Catalyst allows the conservation laws of such systems to be [computed using the `conservationlaws` function](@ref network_analysis_conservation_laws). By using these to reduce the dimensionality of the system, as well as specifying the initial amount of each species, homotopy continuation can again be used to find steady states. Here we do this by designating such an initial condition (which is used to compute the system's conserved quantities, in this case $X1 + X2$): +Catalyst allows the conservation laws of such systems to be [computed using the `conservationlaws` function](@ref conservation_laws). By using these to reduce the dimensionality of the system, as well as specifying the initial amount of each species, homotopy continuation can again be used to find steady states. Here we do this by designating such an initial condition (which is used to compute the system's conserved quantities, in this case $X1 + X2$): ```@example hc_claws import HomotopyContinuation # hide ps = [:k1 => 2.0, :k2 => 1.0] diff --git a/docs/src/steady_state_functionality/nonlinear_solve.md b/docs/src/steady_state_functionality/nonlinear_solve.md index bfc4d29530..fb4ae52f0a 100644 --- a/docs/src/steady_state_functionality/nonlinear_solve.md +++ b/docs/src/steady_state_functionality/nonlinear_solve.md @@ -65,7 +65,7 @@ two_state_model = @reaction_network begin (k1,k2), X1 <--> X2 end ``` -It has an infinite number of steady states. To make steady state finding possible, information of the system's conserved quantities (here $C = X1 + X2$) must be provided. Since these can be computed from system initial conditions (`u0`, i.e. those provided when performing ODE simulations), designating an `u0` is often the best way. There are two ways to do this. First, one can perform [forward ODE simulation-based steady state finding](@ref steady_state_solving_simulation), using the initial condition as the initial `u` guess. Alternatively, any conserved quantities can be eliminated when the `NonlinearProblem` is created. This feature is supported by [Catalyst's conservation law finding and elimination feature](@ref network_analysis_conservation_laws). +It has an infinite number of steady states. To make steady state finding possible, information of the system's conserved quantities (here $C = X1 + X2$) must be provided. Since these can be computed from system initial conditions (`u0`, i.e. those provided when performing ODE simulations), designating an `u0` is often the best way. There are two ways to do this. First, one can perform [forward ODE simulation-based steady state finding](@ref steady_state_solving_simulation), using the initial condition as the initial `u` guess. Alternatively, any conserved quantities can be eliminated when the `NonlinearProblem` is created. This feature is supported by [Catalyst's conservation law finding and elimination feature](@ref conservation_laws). To eliminate conservation laws we simply provide the `remove_conserved = true` argument to `NonlinearProblem`: ```@example steady_state_solving_claws