From 361311e222addbb142d317f968e10b7e431cc267 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Sat, 29 Mar 2025 23:16:57 +0000
Subject: [PATCH 01/11] init

---
 HISTORY.md                                    |  3 +
 docs/src/model_creation/conservation_laws.md  |  3 -
 .../bifurcation_kit_extension.jl              |  3 +-
 .../homotopy_continuation_extension.jl        |  3 +-
 .../structural_identifiability_extension.jl   |  2 +-
 src/reactionsystem_conversions.jl             | 54 ++++-----------
 src/steady_state_stability.jl                 |  2 +-
 test/network_analysis/conservation_laws.jl    | 69 ++++---------------
 .../simulation_and_solving/solve_nonlinear.jl |  2 +-
 9 files changed, 37 insertions(+), 104 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index 902c6a597b..a3a3a1cbd1 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -13,6 +13,9 @@
   solver that auto-switches between explicit and implicit methods, install `OrdinaryDiffEqDefault`. To 
   use `Tsit5` install `OrdinaryDiffEqTsit5`, etc. The possible sub-libraries, each containing different solvers,
   can be viewed [here](https://github.com/SciML/OrdinaryDiffEq.jl/tree/master/lib).
+- It should now be safe to use `remake` on problems which have had conservation laws removed.
+  The warning regarding this when eliminating conservation laws have been removed (in conjecture
+  with the `remove_conserved_warn` argument for disabling this warning).
 - New formula for inferring variables from equations (declared using the `@equations` options) in the DSL. The order of inference of species/variables/parameters is now:
     (1) Every symbol explicitly declared using `@species`, `@variables`, and `@parameters` are assigned to the correct category.
     (2) Every symbol used as a reaction reactant is inferred as a species.
diff --git a/docs/src/model_creation/conservation_laws.md b/docs/src/model_creation/conservation_laws.md
index 0fe504d372..6086952d25 100644
--- a/docs/src/model_creation/conservation_laws.md
+++ b/docs/src/model_creation/conservation_laws.md
@@ -61,9 +61,6 @@ 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₂]
diff --git a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
index 5bd5355b25..6528a93041 100644
--- a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
+++ b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
@@ -37,7 +37,6 @@ function bkext_make_nsys(rs, u0)
     cons_default = [cons_eq.rhs for cons_eq in cons_eqs]
     cons_default = Catalyst.get_networkproperties(rs).conservedconst => cons_default
     defaults = Dict([u0; cons_default])
-    nsys = convert(NonlinearSystem, rs; defaults,
-        remove_conserved = true, remove_conserved_warn = false)
+    nsys = convert(NonlinearSystem, rs; defaults, remove_conserved = true)
     return complete(nsys)
 end
diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
index 13bb916dbb..4e70b619f7 100644
--- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
+++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
@@ -51,8 +51,7 @@ function steady_state_polynomial(rs::ReactionSystem, ps, u0)
     # Creates the appropriate nonlinear system, and converts parameters to a form that can
     # be substituted in later.
     rs = Catalyst.expand_registered_functions(rs)
-    ns = complete(convert(NonlinearSystem, rs;
-        remove_conserved = true, remove_conserved_warn = false))
+    ns = complete(convert(NonlinearSystem, rs; remove_conserved = true))
     pre_varmap = [symmap_to_varmap(rs, u0)..., symmap_to_varmap(rs, ps)...]
     Catalyst.conservationlaw_errorcheck(rs, pre_varmap)
     p_dict = make_p_val_dict(pre_varmap, rs, ns)
diff --git a/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl b/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl
index 3b210ea69d..b9189b68fc 100644
--- a/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl
+++ b/ext/CatalystStructuralIdentifiabilityExtension/structural_identifiability_extension.jl
@@ -169,7 +169,7 @@ function make_osys(rs::ReactionSystem; remove_conserved = true)
         error("Identifiability should only be computed for complete systems. A ReactionSystem can be marked as complete using the `complete` function.")
     end
     rs = complete(Catalyst.expand_registered_functions(flatten(rs)))
-    osys = complete(convert(ODESystem, rs; remove_conserved, remove_conserved_warn = false))
+    osys = complete(convert(ODESystem, rs; remove_conserved))
     vars = [unknowns(rs); parameters(rs)]
 
     # Computes equations for system conservation laws.
diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl
index 7211492094..b3e7ade1e5 100644
--- a/src/reactionsystem_conversions.jl
+++ b/src/reactionsystem_conversions.jl
@@ -449,21 +449,6 @@ diff_2_zero(expr) = (Symbolics.is_derivative(expr) ? 0 : expr)
 
 COMPLETENESS_ERROR = "A ReactionSystem must be complete before it can be converted to other system types. A ReactionSystem can be marked as complete using the `complete` function."
 
-# Used to, when required, display a warning about conservation law removal and remake.
-function check_cons_warning(remove_conserved, remove_conserved_warn)
-    (remove_conserved && remove_conserved_warn) || return
-    @warn "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 ($CONSERVED_CONSTANT_SYMBOL), 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`."
-end
 
 ### System Conversions ###
 
@@ -482,19 +467,16 @@ Keyword args and default values:
 - `remove_conserved=false`, if set to `true` will calculate conservation laws of the
   underlying set of reactions (ignoring constraint equations), and then apply them to reduce
   the number of equations.
-- `remove_conserved_warn = true`: If `true`, if also `remove_conserved = true`, there will be
-  a warning regarding limitations of modifying problems generated from the created system.
 """
 function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        include_zero_odes = true, remove_conserved = false, remove_conserved_warn = true,
+        include_zero_odes = true, remove_conserved = false,
         checks = false, default_u0 = Dict(), default_p = Dict(),
         defaults = _merge(Dict(default_u0), Dict(default_p)),
         kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, ODESystem)
-    check_cons_warning(remove_conserved, remove_conserved_warn)
 
     fullrs = Catalyst.flatten(rs)
     remove_conserved && conservationlaws(fullrs)
@@ -529,19 +511,15 @@ Keyword args and default values:
 - `remove_conserved=false`, if set to `true` will calculate conservation laws of the
   underlying set of reactions (ignoring constraint equations), and then apply them to reduce
   the number of equations.
-- `remove_conserved_warn = true`: If `true`, if also `remove_conserved = true`, there will be
-  a warning regarding limitations of modifying problems generated from the created system.
 """
 function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         include_zero_odes = true, remove_conserved = false, checks = false,
-        remove_conserved_warn = true, default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
         all_differentials_permitted = false, kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, NonlinearSystem)
-    check_cons_warning(remove_conserved, remove_conserved_warn)
     if !isautonomous(rs)
         error("Attempting to convert a non-autonomous `ReactionSystem` (e.g. where some rate depend on $(get_iv(rs))) to a `NonlinearSystem`. This is not possible. if you are intending to compute system steady states, consider creating and solving a `SteadyStateProblem.")
     end
@@ -612,19 +590,15 @@ Notes:
 - `remove_conserved=false`, if set to `true` will calculate conservation laws of the
   underlying set of reactions (ignoring constraint equations), and then apply them to reduce
   the number of equations.
-- `remove_conserved_warn = true`: If `true`, if also `remove_conserved = true`, there will be
-  a warning regarding limitations of modifying problems generated from the created system.
 """
 function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
         name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         include_zero_odes = true, checks = false, remove_conserved = false,
-        remove_conserved_warn = true, default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
         kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, SDESystem)
-    check_cons_warning(remove_conserved, remove_conserved_warn)
 
     flatrs = Catalyst.flatten(rs)
 
@@ -708,12 +682,12 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan,
         p = DiffEqBase.NullParameters(), args...;
         check_length = false, name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        include_zero_odes = true, remove_conserved = false, remove_conserved_warn = true,
-        checks = false, structural_simplify = false, kwargs...)
+        include_zero_odes = true, remove_conserved = false, checks = false, 
+        structural_simplify = false, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks,
-        remove_conserved, remove_conserved_warn)
+        remove_conserved)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
     if structural_simplify
@@ -732,12 +706,12 @@ function DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0,
         p = DiffEqBase.NullParameters(), args...;
         name = nameof(rs), include_zero_odes = true,
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        remove_conserved = false, remove_conserved_warn = true, checks = false,
-        check_length = false, all_differentials_permitted = false, kwargs...)
+        remove_conserved = false, checks = false, check_length = false, 
+        all_differentials_permitted = false, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     nlsys = convert(NonlinearSystem, rs; name, combinatoric_ratelaws, include_zero_odes,
-        checks, all_differentials_permitted, remove_conserved, remove_conserved_warn)
+        checks, all_differentials_permitted, remove_conserved)
     nlsys = complete(nlsys)
     return NonlinearProblem(nlsys, u0map, pmap, args...; check_length,
         kwargs...)
@@ -748,11 +722,11 @@ function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan,
         p = DiffEqBase.NullParameters(), args...;
         name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         include_zero_odes = true, checks = false, check_length = false, remove_conserved = false,
-        remove_conserved_warn = true, structural_simplify = false, kwargs...)
+        structural_simplify = false, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     sde_sys = convert(SDESystem, rs; name, combinatoric_ratelaws,
-        include_zero_odes, checks, remove_conserved, remove_conserved_warn)
+        include_zero_odes, checks, remove_conserved)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
     if structural_simplify
@@ -876,12 +850,12 @@ function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0,
         p = DiffEqBase.NullParameters(), args...;
         check_length = false, name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        remove_conserved = false, remove_conserved_warn = true, include_zero_odes = true,
-        checks = false, structural_simplify = false, kwargs...)
+        remove_conserved = false, include_zero_odes = true, checks = false, 
+        structural_simplify = false, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
     osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks,
-        remove_conserved, remove_conserved_warn)
+        remove_conserved)
 
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
     if structural_simplify
diff --git a/src/steady_state_stability.jl b/src/steady_state_stability.jl
index de9a229ef8..707b8a79e4 100644
--- a/src/steady_state_stability.jl
+++ b/src/steady_state_stability.jl
@@ -118,7 +118,7 @@ function steady_state_jac(rs::ReactionSystem; u0 = [sp => 0.0 for sp in unknowns
     # Creates an `ODEProblem` with a Jacobian. Dummy values for `u0` and `ps` must be provided.
     ps = [p => 0.0 for p in parameters(rs)]
     return ODEProblem(rs, u0, 0, ps; jac = true, combinatoric_ratelaws,
-        remove_conserved = true, remove_conserved_warn = false)
+        remove_conserved = true)
 end
 
 # Converts a `u` vector from a vector of values to a map.
diff --git a/test/network_analysis/conservation_laws.jl b/test/network_analysis/conservation_laws.jl
index 9793d6f5d9..14cffbce22 100644
--- a/test/network_analysis/conservation_laws.jl
+++ b/test/network_analysis/conservation_laws.jl
@@ -11,9 +11,6 @@ seed = rand(rng, 1:100)
 # Fetch test networks.
 include("../test_networks.jl")
 
-# Except where we test the warnings, we do not want to print this warning.
-remove_conserved_warn = false
-
 ### Basic Tests ###
 
 # Tests basic functionality on system with known conservation laws.
@@ -117,10 +114,10 @@ let
     tspan = (0.0, 20.0)
 
     # Simulates model using ODEs and checks that simulations are identical.
-    osys = complete(convert(ODESystem, rn; remove_conserved = true, remove_conserved_warn))
+    osys = complete(convert(ODESystem, rn; remove_conserved = true))
     oprob1 = ODEProblem(osys, u0, tspan, p)
     oprob2 = ODEProblem(rn, u0, tspan, p)
-    oprob3 = ODEProblem(rn, u0, tspan, p; remove_conserved = true, remove_conserved_warn)
+    oprob3 = ODEProblem(rn, u0, tspan, p; remove_conserved = true)
     osol1 = solve(oprob1, Tsit5(); abstol = 1e-8, reltol = 1e-8, saveat= 0.2)
     osol2 = solve(oprob2, Tsit5(); abstol = 1e-8, reltol = 1e-8, saveat= 0.2)
     osol3 = solve(oprob3, Tsit5(); abstol = 1e-8, reltol = 1e-8, saveat= 0.2)
@@ -128,13 +125,13 @@ let
 
     # Checks that steady states found using nonlinear solving and steady state simulations are identical.
     @test_broken begin # Conservation law removal currently not working for NonlinearSystems due to MTK depricating something. https://github.com/SciML/ModelingToolkit.jl/issues/3458, https://github.com/SciML/ModelingToolkit.jl/issues/3411
-        nsys = complete(convert(NonlinearSystem, rn; remove_conserved = true, remove_conserved_warn))
+        nsys = complete(convert(NonlinearSystem, rn; remove_conserved = true))
         nprob1 = NonlinearProblem{true}(nsys, u0, p; guesses = [nsys.Γ => [0.0]])
         nprob2 = NonlinearProblem(rn, u0, p)
-        nprob3 = NonlinearProblem(rn, u0, p; remove_conserved = true, remove_conserved_warn)
+        nprob3 = NonlinearProblem(rn, u0, p; remove_conserved = true)
         ssprob1 = SteadyStateProblem{true}(osys, u0, p)
         ssprob2 = SteadyStateProblem(rn, u0, p)
-        ssprob3 = SteadyStateProblem(rn, u0, p; remove_conserved = true, remove_conserved_warn)
+        ssprob3 = SteadyStateProblem(rn, u0, p; remove_conserved = true)
         nsol1 = solve(nprob1, NewtonRaphson(); abstol = 1e-8)
         # Nonlinear problems cannot find steady states properly without removing conserved species.
         nsol3 = solve(nprob3, NewtonRaphson(); abstol = 1e-8)
@@ -147,10 +144,10 @@ let
     # Creates SDEProblems using various approaches.
     u0_sde = [A => 100.0, B => 20.0, C => 5.0, D => 10.0, E => 3.0, F1 => 8.0, F2 => 2.0,
         F3 => 20.0]
-    ssys = complete(convert(SDESystem, rn; remove_conserved = true, remove_conserved_warn))
+    ssys = complete(convert(SDESystem, rn; remove_conserved = true))
     sprob1 = SDEProblem(ssys, u0_sde, tspan, p)
     sprob2 = SDEProblem(rn, u0_sde, tspan, p)
-    sprob3 = SDEProblem(rn, u0_sde, tspan, p; remove_conserved = true, remove_conserved_warn)
+    sprob3 = SDEProblem(rn, u0_sde, tspan, p; remove_conserved = true)
 
     # Checks that the SDEs f and g function evaluates to the same thing.
     ind_us = ModelingToolkit.get_unknowns(ssys)
@@ -191,9 +188,9 @@ let
     u0_2 = [rn.X => 2.0, rn.Y => 3.0, rn.XY => 4.0]
     u0_3 = [:X => 2.0, :Y => 3.0, :XY => 4.0]
     ps = (kB => 2, kD => 1.5)
-    oprob1 = ODEProblem(rn, u0_1, 10.0, ps; remove_conserved = true, remove_conserved_warn)
-    oprob2 = ODEProblem(rn, u0_2, 10.0, ps; remove_conserved = true, remove_conserved_warn)
-    oprob3 = ODEProblem(rn, u0_3, 10.0, ps; remove_conserved = true, remove_conserved_warn)
+    oprob1 = ODEProblem(rn, u0_1, 10.0, ps; remove_conserved = true)
+    oprob2 = ODEProblem(rn, u0_2, 10.0, ps; remove_conserved = true)
+    oprob3 = ODEProblem(rn, u0_3, 10.0, ps; remove_conserved = true)
     @test solve(oprob1)[sps] ≈ solve(oprob2)[sps] ≈ solve(oprob3)[sps]
 end
 
@@ -205,7 +202,7 @@ let
     end
     u0 = Dict([:X1 => 100.0, :X2 => 120.0])
     ps = [:k1 => 0.2, :k2 => 0.15]
-    sprob = SDEProblem(rn, u0, 10.0, ps; remove_conserved = true, remove_conserved_warn)
+    sprob = SDEProblem(rn, u0, 10.0, ps; remove_conserved = true)
 
     # Checks that conservation laws hold in all simulations.
     sol = solve(sprob, ImplicitEM(); seed)
@@ -218,7 +215,7 @@ let
     rn = @reaction_network begin
         (k1,k2), X1 <--> X2
     end
-    osys = complete(convert(ODESystem, rn; remove_conserved = true, remove_conserved_warn))
+    osys = complete(convert(ODESystem, rn; remove_conserved = true))
     u0 = [osys.X1 => 1.0, osys.X2 => 1.0]
     ps_1 = [osys.k1 => 2.0, osys.k2 => 3.0]
     ps_2 = [osys.k1 => 2.0, osys.k2 => 3.0, osys.Γ => [4.0]]
@@ -245,7 +242,7 @@ let
         Reaction(k2, [X2], [X1])
     ]
     @named rs = ReactionSystem(rxs, t)
-    osys = convert(ODESystem, complete(rs); remove_conserved = true, remove_conserved_warn = false)
+    osys = convert(ODESystem, complete(rs); remove_conserved = true)
     osys = complete(osys)
     @unpack Γ = osys
 
@@ -288,7 +285,7 @@ let
     rn = @reaction_network begin
         (k1,k2), X1 <--> X2
     end
-    @test_throws ArgumentError convert(JumpSystem, rn; remove_conserved = true, remove_conserved_warn)
+    @test_throws ArgumentError convert(JumpSystem, rn; remove_conserved = true)
 end
 
 # Checks that `conserved` metadata is added correctly to parameters.
@@ -299,7 +296,7 @@ let
         (k1,k2), X1 <--> X2
         (k1,k2), Y1 <--> Y2
     end
-    osys = convert(ODESystem, rs; remove_conserved = true, remove_conserved_warn)
+    osys = convert(ODESystem, rs; remove_conserved = true)
 
     # Checks that the correct parameters have the `conserved` metadata.
     @test Catalyst.isconserved(osys.Γ[1])
@@ -308,42 +305,6 @@ let
     @test !Catalyst.isconserved(osys.k2)
 end
 
-# Checks that conservation law elimination warnings are generated in the correct cases.
-let
-    # Prepare model.
-    rn = @reaction_network begin
-        (k1,k2), X1 <--> X2
-    end
-    u0 = [:X1 => 1.0, :X2 => 2.0]
-    tspan = (0.0, 1.0)
-    ps = [:k1 => 3.0, :k2 => 4.0]
-
-    # Check warnings in system conversion.
-    for XSystem in [ODESystem, SDESystem, NonlinearSystem]
-        @test_nowarn convert(XSystem, rn)
-        @test_logs (:warn, r"You are creating a system or problem while eliminating conserved quantities. Please *") convert(XSystem, rn; remove_conserved = true)
-        @test_nowarn convert(XSystem, rn; remove_conserved_warn = false)
-        @test_nowarn convert(XSystem, rn; remove_conserved = true, remove_conserved_warn = false)
-    end
-
-    # Checks during problem creation (separate depending on whether they have a time span or not).
-    for XProblem in [ODEProblem, SDEProblem]
-        @test_nowarn XProblem(rn, u0, tspan, ps)
-        @test_logs (:warn, r"You are creating a system or problem while eliminating conserved quantities. Please *") XProblem(rn, u0, tspan, ps; remove_conserved = true)
-        @test_nowarn XProblem(rn, u0, tspan, ps; remove_conserved_warn = false)
-        @test_nowarn XProblem(rn, u0, tspan, ps; remove_conserved = true, remove_conserved_warn = false)
-    end
-    # Conservation law removal currently not working for NonlinearSystems due to MTK depricating something. https://github.com/SciML/ModelingToolkit.jl/issues/3458, https://github.com/SciML/ModelingToolkit.jl/issues/3411
-    # Currently just checks whether a NonlinearProblem can be created.
-    @test_broken for XProblem in [NonlinearProblem, SteadyStateProblem]
-        XProblem(rn, u0, ps; remove_conserved = true)
-        # @test_nowarn XProblem(rn, u0, ps)
-        # @test_logs (:warn, r"You are creating a system or problem while eliminating conserved quantities. Please *") XProblem(rn, u0, ps; remove_conserved = true)
-        # @test_nowarn XProblem(rn, u0, ps; remove_conserved_warn = false)
-        # @test_nowarn XProblem(rn, u0, ps; remove_conserved = true, remove_conserved_warn = false)
-    end
-end
-
 # Conservation law simulations for vectorised species.
 let
     # Prepares the model.
diff --git a/test/simulation_and_solving/solve_nonlinear.jl b/test/simulation_and_solving/solve_nonlinear.jl
index 05ec4d327a..751ee75476 100644
--- a/test/simulation_and_solving/solve_nonlinear.jl
+++ b/test/simulation_and_solving/solve_nonlinear.jl
@@ -90,7 +90,7 @@ end
     # Creates NonlinearProblem.
     u0 = [steady_state_network_3.X => rand(), steady_state_network_3.Y => rand() + 1.0, steady_state_network_3.Y2 => rand() + 3.0, steady_state_network_3.XY2 => 0.0]
     p = [:p => rand()+1.0, :d => 0.5, :k1 => 1.0, :k2 => 2.0, :k3 => 3.0, :k4 => 4.0]
-    nl_prob_1 = NonlinearProblem(steady_state_network_3, u0, p; remove_conserved = true, remove_conserved_warn = false)
+    nl_prob_1 = NonlinearProblem(steady_state_network_3, u0, p; remove_conserved = true)
     nl_prob_2 = NonlinearProblem(steady_state_network_3, u0, p)
 
     # Solves it using standard algorithm and simulation based algorithm.

From b9e46fb00d6a2161cbdbbe550495d809410c4e71 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Sat, 29 Mar 2025 23:37:04 +0000
Subject: [PATCH 02/11] add some tests

---
 test/upstream/mtk_structure_indexing.jl | 20 +++++++++++++++++++-
 1 file changed, 19 insertions(+), 1 deletion(-)

diff --git a/test/upstream/mtk_structure_indexing.jl b/test/upstream/mtk_structure_indexing.jl
index 30e4eb865a..07d3893da1 100644
--- a/test/upstream/mtk_structure_indexing.jl
+++ b/test/upstream/mtk_structure_indexing.jl
@@ -339,10 +339,13 @@ let
 
     # Checks for both ODE and SDE problems.
     for (XProblem, solver) in zip([ODEProblem, SDEProblem], [Tsit5(), ImplicitEM()])
-        # Create the problems which we wish to check..
+        # Create the problem which we wish to check.
         u0 = [X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, W => 6.0]
         ps = [k1 => 0.1, k2 => 0.2, V0 => 3.0]
         prob1 = XProblem(rs, u0, 0.001, ps; remove_conserved = true)
+        Γ = prob1.f.sys.Γ
+        
+        # Creates various `remake` version of the problem.
         prob2 = remake(prob1, u0 = [X1 => 10.0])
         prob3 = remake(prob2, u0 = [X2 => 20.0])
         prob4 = remake(prob1, u0 = [X2 => 20.0, Y1 => 30.0])
@@ -350,6 +353,9 @@ let
         prob6 = remake(prob1, u0 = [Y2 => 40.0], p = [k1 => 0.4])
         prob7 = remake(prob1, u0 = [X1 => 10.0, X2 => 20.0], p = [V0 => 50.0])
         prob8 = remake(prob1, u0 = [W => 60.0])
+        prob9 = remake(prob2; p = [Γ => [10.0, 20.0]])        
+        prob10 = remake(prob1; u0 = [Y1 => 20.0], p = [Γ => [20.0, 30.0], k1 => 0.4])
+        prob11 = remake(prob10, u0 = [X1 => 10.0], p = [k2 => 0.5])
 
         # Creates a testing function.
         function test_vals(prob, us_correct::Dict, ps_correct::Dict)
@@ -393,6 +399,18 @@ let
         test_vals(prob8,
             Dict(X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 3.0, W => 60.0),
             Dict(k1 => 0.1, k2 => 0.2, V0 => 3.0, v => 3.0, w => 60.0, Γ[1] => 3.0, Γ[2] => 7.0))
+        test_vals(prob8,
+            Dict(X1 => 1.0, X2 => 2.0, Y1 => 3.0, Y2 => 4.0, V => 3.0, W => 60.0),
+            Dict(k1 => 0.1, k2 => 0.2, V0 => 3.0, v => 3.0, w => 60.0, Γ[1] => 3.0, Γ[2] => 7.0))
+        test_vals(prob9,
+            Dict(X1 => 10.0, X2 => 0.0, Y1 => 3.0, Y2 => 17.0, V => 3.0, W => 6.0),
+            Dict(k1 => 0.1, k2 => 0.2, V0 => 3.0, v => 3.0, w => 6.0, Γ[1] => 10.0, Γ[2] => 20.0))
+        test_vals(prob10,
+            Dict(X1 => 1.0, X2 => 19.0, Y1 => 20.0, Y2 => 10.0, V => 3.0, W => 6.0),
+            Dict(k1 => 0.4, k2 => 0.2, V0 => 3.0, v => 3.0, w => 6.0, Γ[1] => 20.0, Γ[2] => 30.0))
+        test_vals(prob11,
+            Dict(X1 => 10.0, X2 => 10.0, Y1 => 20.0, Y2 => 10.0, V => 3.0, W => 6.0),
+            Dict(k1 => 0.4, k2 => 0.5, V0 => 3.0, v => 3.0, w => 6.0, Γ[1] => 20.0, Γ[2] => 30.0))
     end
 end
 

From a340eade722380efa075cf08170f5bafa4f06fe0 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Sun, 30 Mar 2025 12:06:12 +0100
Subject: [PATCH 03/11] add new test

---
 test/network_analysis/conservation_laws.jl | 78 ++++++++++++++++++++++
 1 file changed, 78 insertions(+)

diff --git a/test/network_analysis/conservation_laws.jl b/test/network_analysis/conservation_laws.jl
index 14cffbce22..55c987b7a3 100644
--- a/test/network_analysis/conservation_laws.jl
+++ b/test/network_analysis/conservation_laws.jl
@@ -229,6 +229,8 @@ let
     @test all(sol2[osys.X1 + osys.X2] .== 4.0)
 end
 
+### Problem `remake`ing and Updating Tests ###
+
 # Tests system problem updating when conservation laws are eliminated.
 # Checks that the correct values are used after the conservation law species are updated.
 # Here is an issue related to the broken tests: https://github.com/SciML/Catalyst.jl/issues/952
@@ -278,6 +280,82 @@ let
     @test integrator.ps[k2] == 0.6
 end
 
+# Goes through a chain of updating of conservation law constants/species, checking that
+# the values of relevant quantitites are correct after each step.
+# Generally, if `Γ` has not been explicitly updated, it will be updated to accomodate new species
+# values. If it has been explicitly udpated, the corresponding elimianted quantity will have its
+# vaue updated to accomodate new Γ/species values (also, the elimianted species's value can not longer be changed).
+let
+    # Prepares an `ODEProblem` with conserved quantities eliminated. Computes the conservation equation.
+    rn = @reaction_network begin
+        (k1,k2), 2X1 <--> X2
+        (k3,k4), X1 + X2 <--> X3
+    end
+    @unpack X1, X2, X3 = rn
+    u0 = [X1 => 1.0, X2 => 1.0, X3 => 1.0]
+    ps = [:k1 => 0.1, :k2 => 0.2, :k3 => 0.3, :k4 => 0.4]
+    prob_old = ODEProblem(rn, u0, 1.0, ps; remove_conserved = true)
+    conserved_quantity = conservationlaw_constants(rn)[1].rhs
+
+    # For a couple of iterations, updates the problem, ensuring that when a species is updated:
+    # - Only that species and the conservation constant have their values updated.
+    # The `≈` is because sometimes the computed values will not be fully exact.
+    for _ = 1:3
+        # Updates X2, checks the values of all species and Γ, then sets which is the old problem.
+        X2_new = rand(rng, 1.0:10.0)
+        prob_new = remake(prob_old; u0 = [:X2 => X2_new])
+        @test prob_old[:X1] ≈ prob_new[:X1]
+        @test X2_new ≈ prob_new[:X2]
+        @test prob_old[:X3] ≈ prob_new[:X3]
+        @test substitute(conserved_quantity, Dict([X1 => prob_old[X1], X2 => X2_new, X3 => prob_old[X3]])) ≈ prob_new.ps[:Γ][1]
+        prob_old = prob_new
+
+        # Updates X3, checks the values of all species and Γ, then sets which is the old problem.
+        X3_new = rand(rng, 1.0:10.0)
+        prob_new = remake(prob_old; u0 = [:X3 => X3_new])
+        @test prob_old[:X1] ≈ prob_new[:X1]
+        @test prob_old[:X2] ≈ prob_new[:X2]
+        @test X3_new ≈ prob_new[:X3]
+        @test substitute(conserved_quantity, Dict([X1 => prob_old[X1], X2 => prob_old[X2], X3 => X3_new])) ≈ prob_new.ps[:Γ][1]
+        prob_old = prob_new
+    end
+
+    # Similarly, but now also updates the conservation constant. Here, once Γ has been updated:
+    # - The conservation law constant will be kept fixed, and secondary updates are made to the
+    # eliminated species.
+    # Assumes that X3 is the eliminated species.
+    # The random Γ is ensured to be large enough not to generate negative values in the eliminated species.
+    for _ in 1:3
+        # Updates Γ, checks the values of all species and Γ, then sets which is the old problem.
+        Γ_new = substitute(conserved_quantity, Dict([X1 => prob_old[X1], X2 => prob_old[X2], X3 => 0])) + rand(rng, 0.0:5.0)
+        prob_new = remake(prob_old; p = [:Γ => [Γ_new]], warn_initialize_determined = false)
+        @test prob_old[:X1] ≈ prob_new[:X1]
+        @test prob_old[:X2] ≈ prob_new[:X2]
+        @test Γ_new ≈ prob_new.ps[:Γ][1]
+        @test substitute(conserved_quantity, Dict([X1 => prob_old[X1], X2 => prob_old[X2], X3 => prob_new[X3]])) ≈ prob_new.ps[:Γ][1]
+        prob_old = prob_new
+
+        # Updates X1 (non-eliminated species), checks the values of all species and Γ, then sets which is the old problem.
+        # Note that now, `X3` will have its value modified (not and `Γ` remains unchanged).
+        X1_new = rand(rng, 1.0:10.0)
+        prob_new = remake(prob_old; u0 = [:X1 => X1_new])
+        @test X1_new ≈ prob_new[:X1]
+        @test prob_old[:X2] ≈ prob_new[:X2]
+        @test prob_old.ps[:Γ][1] ≈ prob_new.ps[:Γ][1]
+        @test substitute(conserved_quantity, Dict([X1 => X1_new, X2 => prob_old[X2], X3 => prob_new[X3]])) ≈ prob_new.ps[:Γ][1]
+        prob_old = prob_new
+
+        # Updates X3 (the eliminated species). Right now, this will have no effect on `X3` (or the system).
+        X3_new = rand(rng, 1.0:10.0)
+        prob_new = remake(prob_old; u0 = [:X3 => X3_new], warn_initialize_determined = false)
+        @test prob_old[:X1] ≈ prob_new[:X1]
+        @test prob_old[:X2] ≈ prob_new[:X2]
+        @test prob_old[:X3] ≈ prob_new[:X3]
+        @test prob_old.ps[:Γ][1] ≈ prob_new.ps[:Γ][1]
+        prob_old = prob_new
+    end
+end
+
 ### Other Tests ###
 
 # Checks that `JumpSystem`s with conservation laws cannot be generated.

From 81fbbe80d2a37a357e988cb70df88dc1d09b155b Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Sun, 30 Mar 2025 12:16:08 +0100
Subject: [PATCH 04/11] save doc progress

---
 docs/src/model_creation/conservation_laws.md | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/docs/src/model_creation/conservation_laws.md b/docs/src/model_creation/conservation_laws.md
index 6086952d25..1d0533010d 100644
--- a/docs/src/model_creation/conservation_laws.md
+++ b/docs/src/model_creation/conservation_laws.md
@@ -86,3 +86,14 @@ Finally, the `conservationlaws` function yields a $m \times n$ matrix, where $n$
 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))
+```
+Ad advantage of this is that we can set the conserved quantity#s value directly in simulations. E.g. 
\ No newline at end of file

From fdabb28a71fa303bc31e281c9b4a70a43040486b Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Sun, 30 Mar 2025 13:33:23 +0100
Subject: [PATCH 05/11] init

---
 src/reactionsystem_conversions.jl | 30 +++++++++---------------------
 1 file changed, 9 insertions(+), 21 deletions(-)

diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl
index 7211492094..e0b322b21a 100644
--- a/src/reactionsystem_conversions.jl
+++ b/src/reactionsystem_conversions.jl
@@ -710,8 +710,6 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan,
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         include_zero_odes = true, remove_conserved = false, remove_conserved_warn = true,
         checks = false, structural_simplify = false, kwargs...)
-    u0map = symmap_to_varmap(rs, u0)
-    pmap = symmap_to_varmap(rs, p)
     osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks,
         remove_conserved, remove_conserved_warn)
 
@@ -724,7 +722,7 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan,
         osys = complete(osys)
     end
 
-    return ODEProblem(osys, u0map, tspan, pmap, args...; check_length, kwargs...)
+    return ODEProblem(osys, u0, tspan, p, args...; check_length, kwargs...)
 end
 
 # NonlinearProblem from AbstractReactionNetwork
@@ -734,12 +732,10 @@ function DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0,
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         remove_conserved = false, remove_conserved_warn = true, checks = false,
         check_length = false, all_differentials_permitted = false, kwargs...)
-    u0map = symmap_to_varmap(rs, u0)
-    pmap = symmap_to_varmap(rs, p)
     nlsys = convert(NonlinearSystem, rs; name, combinatoric_ratelaws, include_zero_odes,
         checks, all_differentials_permitted, remove_conserved, remove_conserved_warn)
     nlsys = complete(nlsys)
-    return NonlinearProblem(nlsys, u0map, pmap, args...; check_length,
+    return NonlinearProblem(nlsys, u0, p, args...; check_length,
         kwargs...)
 end
 
@@ -749,8 +745,6 @@ function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan,
         name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         include_zero_odes = true, checks = false, check_length = false, remove_conserved = false,
         remove_conserved_warn = true, structural_simplify = false, kwargs...)
-    u0map = symmap_to_varmap(rs, u0)
-    pmap = symmap_to_varmap(rs, p)
     sde_sys = convert(SDESystem, rs; name, combinatoric_ratelaws,
         include_zero_odes, checks, remove_conserved, remove_conserved_warn)
 
@@ -764,7 +758,7 @@ function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan,
     end
 
     p_matrix = zeros(length(get_unknowns(sde_sys)), numreactions(rs))
-    return SDEProblem(sde_sys, u0map, tspan, pmap, args...; check_length,
+    return SDEProblem(sde_sys, u0, tspan, p, args...; check_length,
         noise_rate_prototype = p_matrix, kwargs...)
 end
 
@@ -784,8 +778,8 @@ struct JumpInputs{S <: MT.JumpSystem, T <: SciMLBase.AbstractODEProblem}
 end
 
 """
-    jumpinput = JumpInputs(rs::ReactionSystem, u0map, tspan,
-                    pmap = DiffEqBase.NullParameters;
+    jumpinput = JumpInputs(rs::ReactionSystem, u0, tspan,
+                    p = DiffEqBase.NullParameters;
                     name = nameof(rs),
                     combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
                     checks = false, kwargs...)
@@ -817,13 +811,11 @@ plot(sol, idxs = :A)
 function JumpInputs(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters();
         name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         checks = false, kwargs...)
-    u0map = symmap_to_varmap(rs, u0)
-    pmap = symmap_to_varmap(rs, p)
     jsys = complete(convert(JumpSystem, rs; name, combinatoric_ratelaws, checks))
     if MT.has_variableratejumps(jsys)
-        prob = ODEProblem(jsys, u0map, tspan, pmap; kwargs...)
+        prob = ODEProblem(jsys, u0, tspan, p; kwargs...)
     else
-        prob = DiscreteProblem(jsys, u0map, tspan, pmap; kwargs...)
+        prob = DiscreteProblem(jsys, u0, tspan, p; kwargs...)
     end
     JumpInputs(jsys, prob)
 end
@@ -847,11 +839,9 @@ function DiffEqBase.DiscreteProblem(rs::ReactionSystem, u0, tspan::Tuple,
         name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         checks = false, kwargs...)
-    u0map = symmap_to_varmap(rs, u0)
-    pmap = symmap_to_varmap(rs, p)
     jsys = convert(JumpSystem, rs; name, combinatoric_ratelaws, checks)
     jsys = complete(jsys)
-    return DiscreteProblem(jsys, u0map, tspan, pmap, args...; kwargs...)
+    return DiscreteProblem(jsys, u0, tspan, p, args...; kwargs...)
 end
 
 # JumpProblem from AbstractReactionNetwork
@@ -878,8 +868,6 @@ function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0,
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         remove_conserved = false, remove_conserved_warn = true, include_zero_odes = true,
         checks = false, structural_simplify = false, kwargs...)
-    u0map = symmap_to_varmap(rs, u0)
-    pmap = symmap_to_varmap(rs, p)
     osys = convert(ODESystem, rs; name, combinatoric_ratelaws, include_zero_odes, checks,
         remove_conserved, remove_conserved_warn)
 
@@ -892,7 +880,7 @@ function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0,
         osys = complete(osys)
     end
 
-    return SteadyStateProblem(osys, u0map, pmap, args...; check_length, kwargs...)
+    return SteadyStateProblem(osys, u0, p, args...; check_length, kwargs...)
 end
 
 ### Symbolic Variable/Symbol Conversions ###

From 0fcb51dca501addd56797e9339195be5f8f33d05 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Sun, 30 Mar 2025 15:44:39 +0100
Subject: [PATCH 06/11] Updated docs

---
 docs/src/model_creation/conservation_laws.md | 59 +++++++++++++++++---
 1 file changed, 51 insertions(+), 8 deletions(-)

diff --git a/docs/src/model_creation/conservation_laws.md b/docs/src/model_creation/conservation_laws.md
index 1d0533010d..59efb7b44b 100644
--- a/docs/src/model_creation/conservation_laws.md
+++ b/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. 
\ No newline at end of file
+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
+```
\ No newline at end of file

From f2a88838205573a1abce5e7be7fdd02274c1237e Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Sun, 30 Mar 2025 17:35:12 +0200
Subject: [PATCH 07/11] Update HISTORY.md

Co-authored-by: Sam Isaacson <isaacsas@users.noreply.github.com>
---
 HISTORY.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/HISTORY.md b/HISTORY.md
index a3a3a1cbd1..4f8e26dd14 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -14,7 +14,7 @@
   use `Tsit5` install `OrdinaryDiffEqTsit5`, etc. The possible sub-libraries, each containing different solvers,
   can be viewed [here](https://github.com/SciML/OrdinaryDiffEq.jl/tree/master/lib).
 - It should now be safe to use `remake` on problems which have had conservation laws removed.
-  The warning regarding this when eliminating conservation laws have been removed (in conjecture
+  The warning regarding this when eliminating conservation laws have been removed (in conjunction
   with the `remove_conserved_warn` argument for disabling this warning).
 - New formula for inferring variables from equations (declared using the `@equations` options) in the DSL. The order of inference of species/variables/parameters is now:
     (1) Every symbol explicitly declared using `@species`, `@variables`, and `@parameters` are assigned to the correct category.

From 53995928835928d3ae1e5d2aa5c3781cb499e4cd Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Sun, 30 Mar 2025 19:06:11 +0100
Subject: [PATCH 08/11] repeat test for SDEProblem as well

---
 test/network_analysis/conservation_laws.jl | 126 +++++++++++----------
 1 file changed, 66 insertions(+), 60 deletions(-)

diff --git a/test/network_analysis/conservation_laws.jl b/test/network_analysis/conservation_laws.jl
index 55c987b7a3..c404b1bc95 100644
--- a/test/network_analysis/conservation_laws.jl
+++ b/test/network_analysis/conservation_laws.jl
@@ -244,11 +244,10 @@ let
         Reaction(k2, [X2], [X1])
     ]
     @named rs = ReactionSystem(rxs, t)
-    osys = convert(ODESystem, complete(rs); remove_conserved = true)
-    osys = complete(osys)
+    osys = complete(convert(ODESystem, complete(rs); remove_conserved = true))
     @unpack Γ = osys
 
-    # Creates an `ODEProblem`.
+    # Creates the various problem types.
     u0 = [X1 => 1.0, X2 => 2.0]
     ps = [k1 => 0.1, k2 => 0.2]
     oprob = ODEProblem(osys, u0, (0.0, 1.0), ps)
@@ -294,65 +293,72 @@ let
     @unpack X1, X2, X3 = rn
     u0 = [X1 => 1.0, X2 => 1.0, X3 => 1.0]
     ps = [:k1 => 0.1, :k2 => 0.2, :k3 => 0.3, :k4 => 0.4]
-    prob_old = ODEProblem(rn, u0, 1.0, ps; remove_conserved = true)
     conserved_quantity = conservationlaw_constants(rn)[1].rhs
 
-    # For a couple of iterations, updates the problem, ensuring that when a species is updated:
-    # - Only that species and the conservation constant have their values updated.
-    # The `≈` is because sometimes the computed values will not be fully exact.
-    for _ = 1:3
-        # Updates X2, checks the values of all species and Γ, then sets which is the old problem.
-        X2_new = rand(rng, 1.0:10.0)
-        prob_new = remake(prob_old; u0 = [:X2 => X2_new])
-        @test prob_old[:X1] ≈ prob_new[:X1]
-        @test X2_new ≈ prob_new[:X2]
-        @test prob_old[:X3] ≈ prob_new[:X3]
-        @test substitute(conserved_quantity, Dict([X1 => prob_old[X1], X2 => X2_new, X3 => prob_old[X3]])) ≈ prob_new.ps[:Γ][1]
-        prob_old = prob_new
-
-        # Updates X3, checks the values of all species and Γ, then sets which is the old problem.
-        X3_new = rand(rng, 1.0:10.0)
-        prob_new = remake(prob_old; u0 = [:X3 => X3_new])
-        @test prob_old[:X1] ≈ prob_new[:X1]
-        @test prob_old[:X2] ≈ prob_new[:X2]
-        @test X3_new ≈ prob_new[:X3]
-        @test substitute(conserved_quantity, Dict([X1 => prob_old[X1], X2 => prob_old[X2], X3 => X3_new])) ≈ prob_new.ps[:Γ][1]
-        prob_old = prob_new
-    end
-
-    # Similarly, but now also updates the conservation constant. Here, once Γ has been updated:
-    # - The conservation law constant will be kept fixed, and secondary updates are made to the
-    # eliminated species.
-    # Assumes that X3 is the eliminated species.
-    # The random Γ is ensured to be large enough not to generate negative values in the eliminated species.
-    for _ in 1:3
-        # Updates Γ, checks the values of all species and Γ, then sets which is the old problem.
-        Γ_new = substitute(conserved_quantity, Dict([X1 => prob_old[X1], X2 => prob_old[X2], X3 => 0])) + rand(rng, 0.0:5.0)
-        prob_new = remake(prob_old; p = [:Γ => [Γ_new]], warn_initialize_determined = false)
-        @test prob_old[:X1] ≈ prob_new[:X1]
-        @test prob_old[:X2] ≈ prob_new[:X2]
-        @test Γ_new ≈ prob_new.ps[:Γ][1]
-        @test substitute(conserved_quantity, Dict([X1 => prob_old[X1], X2 => prob_old[X2], X3 => prob_new[X3]])) ≈ prob_new.ps[:Γ][1]
-        prob_old = prob_new
-
-        # Updates X1 (non-eliminated species), checks the values of all species and Γ, then sets which is the old problem.
-        # Note that now, `X3` will have its value modified (not and `Γ` remains unchanged).
-        X1_new = rand(rng, 1.0:10.0)
-        prob_new = remake(prob_old; u0 = [:X1 => X1_new])
-        @test X1_new ≈ prob_new[:X1]
-        @test prob_old[:X2] ≈ prob_new[:X2]
-        @test prob_old.ps[:Γ][1] ≈ prob_new.ps[:Γ][1]
-        @test substitute(conserved_quantity, Dict([X1 => X1_new, X2 => prob_old[X2], X3 => prob_new[X3]])) ≈ prob_new.ps[:Γ][1]
-        prob_old = prob_new
-
-        # Updates X3 (the eliminated species). Right now, this will have no effect on `X3` (or the system).
-        X3_new = rand(rng, 1.0:10.0)
-        prob_new = remake(prob_old; u0 = [:X3 => X3_new], warn_initialize_determined = false)
-        @test prob_old[:X1] ≈ prob_new[:X1]
-        @test prob_old[:X2] ≈ prob_new[:X2]
-        @test prob_old[:X3] ≈ prob_new[:X3]
-        @test prob_old.ps[:Γ][1] ≈ prob_new.ps[:Γ][1]
-        prob_old = prob_new
+    # Loops through the tests for different problem types.
+    oprob_old = ODEProblem(rn, u0, 1.0, ps; remove_conserved = true)
+    sprob_old = SDEProblem(rn, u0, 1.0, ps; remove_conserved = true)
+
+    for prob in [oprob_old, sprob_old]
+        prob_old = prob
+        
+        # For a couple of iterations, updates the problem, ensuring that when a species is updated:
+        # - Only that species and the conservation constant have their values updated.
+        # The `≈` is because sometimes the computed values will not be fully exact.
+        for _ = 1:3
+            # Updates X2, checks the values of all species and Γ, then sets which is the old problem.
+            X2_new = rand(rng, 1.0:10.0)
+            prob_new = remake(prob_old; u0 = [:X2 => X2_new])
+            @test prob_old[:X1] ≈ prob_new[:X1]
+            @test X2_new ≈ prob_new[:X2]
+            @test prob_old[:X3] ≈ prob_new[:X3]
+            @test substitute(conserved_quantity, Dict([X1 => prob_old[X1], X2 => X2_new, X3 => prob_old[X3]])) ≈ prob_new.ps[:Γ][1]
+            prob_old = prob_new
+
+            # Updates X3, checks the values of all species and Γ, then sets which is the old problem.
+            X3_new = rand(rng, 1.0:10.0)
+            prob_new = remake(prob_old; u0 = [:X3 => X3_new])
+            @test prob_old[:X1] ≈ prob_new[:X1]
+            @test prob_old[:X2] ≈ prob_new[:X2]
+            @test X3_new ≈ prob_new[:X3]
+            @test substitute(conserved_quantity, Dict([X1 => prob_old[X1], X2 => prob_old[X2], X3 => X3_new])) ≈ prob_new.ps[:Γ][1]
+            prob_old = prob_new
+        end
+
+        # Similarly, but now also updates the conservation constant. Here, once Γ has been updated:
+        # - The conservation law constant will be kept fixed, and secondary updates are made to the
+        # eliminated species.
+        # Assumes that X3 is the eliminated species.
+        # The random Γ is ensured to be large enough not to generate negative values in the eliminated species.
+        for _ in 1:3
+            # Updates Γ, checks the values of all species and Γ, then sets which is the old problem.
+            Γ_new = substitute(conserved_quantity, Dict([X1 => prob_old[X1], X2 => prob_old[X2], X3 => 0])) + rand(rng, 0.0:5.0)
+            prob_new = remake(prob_old; p = [:Γ => [Γ_new]], warn_initialize_determined = false)
+            @test prob_old[:X1] ≈ prob_new[:X1]
+            @test prob_old[:X2] ≈ prob_new[:X2]
+            @test Γ_new ≈ prob_new.ps[:Γ][1]
+            @test substitute(conserved_quantity, Dict([X1 => prob_old[X1], X2 => prob_old[X2], X3 => prob_new[X3]])) ≈ prob_new.ps[:Γ][1]
+            prob_old = prob_new
+
+            # Updates X1 (non-eliminated species), checks the values of all species and Γ, then sets which is the old problem.
+            # Note that now, `X3` will have its value modified (not and `Γ` remains unchanged).
+            X1_new = rand(rng, 1.0:10.0)
+            prob_new = remake(prob_old; u0 = [:X1 => X1_new])
+            @test X1_new ≈ prob_new[:X1]
+            @test prob_old[:X2] ≈ prob_new[:X2]
+            @test prob_old.ps[:Γ][1] ≈ prob_new.ps[:Γ][1]
+            @test substitute(conserved_quantity, Dict([X1 => X1_new, X2 => prob_old[X2], X3 => prob_new[X3]])) ≈ prob_new.ps[:Γ][1]
+            prob_old = prob_new
+
+            # Updates X3 (the eliminated species). Right now, this will have no effect on `X3` (or the system).
+            X3_new = rand(rng, 1.0:10.0)
+            prob_new = remake(prob_old; u0 = [:X3 => X3_new], warn_initialize_determined = false)
+            @test prob_old[:X1] ≈ prob_new[:X1]
+            @test prob_old[:X2] ≈ prob_new[:X2]
+            @test prob_old[:X3] ≈ prob_new[:X3]
+            @test prob_old.ps[:Γ][1] ≈ prob_new.ps[:Γ][1]
+            prob_old = prob_new
+        end
     end
 end
 

From 32b02041f75a9f5b8d3ec625f77a161fe4fd9fd6 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Sun, 30 Mar 2025 19:49:00 +0100
Subject: [PATCH 09/11] finish missing SDE test stuff

---
 test/network_analysis/conservation_laws.jl | 11 ++++-------
 1 file changed, 4 insertions(+), 7 deletions(-)

diff --git a/test/network_analysis/conservation_laws.jl b/test/network_analysis/conservation_laws.jl
index c404b1bc95..7a0177fe22 100644
--- a/test/network_analysis/conservation_laws.jl
+++ b/test/network_analysis/conservation_laws.jl
@@ -285,7 +285,7 @@ end
 # values. If it has been explicitly udpated, the corresponding elimianted quantity will have its
 # vaue updated to accomodate new Γ/species values (also, the elimianted species's value can not longer be changed).
 let
-    # Prepares an `ODEProblem` with conserved quantities eliminated. Computes the conservation equation.
+    # Prepares the problem inputs and computes the conservation equation.
     rn = @reaction_network begin
         (k1,k2), 2X1 <--> X2
         (k3,k4), X1 + X2 <--> X3
@@ -296,12 +296,9 @@ let
     conserved_quantity = conservationlaw_constants(rn)[1].rhs
 
     # Loops through the tests for different problem types.
-    oprob_old = ODEProblem(rn, u0, 1.0, ps; remove_conserved = true)
-    sprob_old = SDEProblem(rn, u0, 1.0, ps; remove_conserved = true)
-
-    for prob in [oprob_old, sprob_old]
-        prob_old = prob
-        
+    oprob = ODEProblem(rn, u0, 1.0, ps; remove_conserved = true)
+    sprob = SDEProblem(rn, u0, 1.0, ps; remove_conserved = true)
+    for prob_old in [oprob, sprob]
         # For a couple of iterations, updates the problem, ensuring that when a species is updated:
         # - Only that species and the conservation constant have their values updated.
         # The `≈` is because sometimes the computed values will not be fully exact.

From 3557f207905cad1fbaa29c72eef5bf348d5c42c0 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Tue, 1 Apr 2025 10:47:26 +0100
Subject: [PATCH 10/11] update out of date bit in previous conservation law
 tutorial

---
 docs/src/model_creation/conservation_laws.md | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/src/model_creation/conservation_laws.md b/docs/src/model_creation/conservation_laws.md
index 59efb7b44b..8eafb77e66 100644
--- a/docs/src/model_creation/conservation_laws.md
+++ b/docs/src/model_creation/conservation_laws.md
@@ -1,5 +1,5 @@
 # [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)).
+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](@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
@@ -39,7 +39,7 @@ Using the `unknowns` function we can confirm that the ODE only has a single unkn
 ```@example conservation_laws
 unknowns(osys)
 ```
-Next, using `parameters` we note that an additional parameter, `Γ[1]` has been added to the system:
+Next, using `parameters` we note that an additional parameter, `Γ` has been added to the system:
 ```@example conservation_laws
 parameters(osys)
 ```
@@ -53,7 +53,7 @@ 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:
+Here, while `Γ` 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)
@@ -66,7 +66,7 @@ While `X₂` is an observable (and not unknown) of the ODE, we can [access it](@
 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`.
+    Generally, `remove_conserved = true` should not change any modelling 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.

From 5e2362489c4650bf3e26442786740e48ab04a60b Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Thu, 3 Apr 2025 18:02:40 +0200
Subject: [PATCH 11/11] Update docs/src/model_creation/conservation_laws.md

Co-authored-by: Sam Isaacson <isaacsas@users.noreply.github.com>
---
 docs/src/model_creation/conservation_laws.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/model_creation/conservation_laws.md b/docs/src/model_creation/conservation_laws.md
index 8eafb77e66..9adebaa5df 100644
--- a/docs/src/model_creation/conservation_laws.md
+++ b/docs/src/model_creation/conservation_laws.md
@@ -39,7 +39,7 @@ Using the `unknowns` function we can confirm that the ODE only has a single unkn
 ```@example conservation_laws
 unknowns(osys)
 ```
-Next, using `parameters` we note that an additional parameter, `Γ` has been added to the system:
+Next, using `parameters` we note that an additional (vector) parameter, `Γ` has been added to the system:
 ```@example conservation_laws
 parameters(osys)
 ```