Updated docs

Author: TorkelE

docs/src/model_creation/conservation_laws.md

@@ -71,6 +71,9 @@ sol[:X₂]
 !!! 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.
 
+!!! warn
+    Users of the [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl) package might be familiar with the `structural_simplify` function. While it can be applied to Catalyst models as well, generally, this should be avoided (as `remove_conserved` performs a similar role, but is better adapted to these models). Furthermore, applying `structural_simplify` will interfere with conservation law removal, preventing users from accessing eliminated quantities.
+
 ## [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:
@@ -88,12 +91,52 @@ 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.
 
 ## [Updating conservation law values directly](@id conservation_laws_prob_updating)
-Previously we noted that creation law elimination adds the conservation law as a parameter of the system, e.g. as in the case of this dimerisation system:
-```@example conservation_laws_prob_updating
-using Catalyst # hide
-rn = @reaction_network begin
-    (kD,kU), 2X <--> X2
-end
-parameters(convert(ODESystem, rn; remove_conserved = true))
+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:
+```@example conservation_laws
+parameters(convert(ODESystem, rs; remove_conserved = true))
+```
+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₂$).
+```@example conservation_laws
+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)
+nothing # hide
 ```
-Ad advantage of this is that we can set the conserved quantity#s value directly in simulations. E.g. 
+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). 
+!!! note
+    As previously mentioned, the conservation law parameter $Γ$ is a [*vector-valued* parameter](@ref dsl_advanced_options_vector_variables) with one value for each conservation law. That is why we above provide its value as a vector (`:Γ => [5.0]`). If there had been multiple conservation laws, we would provide the `:Γ` vector with one value for each of them (e.g. `:Γ => [10.0, 15.0]`).
+!!! warn
+    If you specify the value of a conservation law parameter, you *must not* specify the value of all species of that conservation law (this can result in an error). Instead, the value of exactly one species must be left unspecified.
+
+Just like when we create a problem, if we [update the species (or conservation law parameter) values of `oprob`](@ref simulation_structure_interfacing_problems), 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₁$.
+```@example conservation_laws
+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]
+```
+```@example conservation_laws
+oprob = remake(oprob; u0 = [:X₁ => 16.0])
+oprob.ps[:Γ][1]
+```
+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:
+```@example conservation_laws
+oprob = remake(oprob; p = [:Γ => [30.0]])
+oprob[:X₂]
+```
+
+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*. 
+
+!!! warn
+    When updating the values of problems with conservation laws it is additionally important to use `remake` (as opposed to direct indexing, e.g. setting `oprob[:X₁] = 16.0`).
+
+### [Extracting the conservation law parameter's symbolic variable](@id conservation_laws_prob_updating_symvar)
+Catalyst represents its models using the [Symbolics.jl](https://github.com/JuliaSymbolics/Symbolics.jl) computer algebraic system, something which allows the user to [form symbolic expressions of model components](@ref simulation_structure_interfacing_symbolic_representation). If you wish to extract and use the symbolic variable corresponding to a model's conserved quantities, you can use the following syntax:
+```@example conservation_laws
+Γ = Catalyst.get_networkproperties(rs).conservedconst
+```