diff --git a/HISTORY.md b/HISTORY.md
index b3cff61329..a5db50e576 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -7,29 +7,47 @@
   (at the time the release is made). If you need a dependency version increased,
   please open an issue and we can update it and make a new Catalyst release once
   testing against the newer dependency version is complete.
-- It is no longer recommended to install and use the full OrdinaryDiffEq library to access specific ODE solvers.
-  Instead, only install the specific OrdinaryDiffEq sub-libraries that contain the desired
-  solver. This significantly reduces installation and package loading times. I.e. to use the default
-  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 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.
-    (2) Every symbol used as a reaction reactant is inferred as a species.
-    (3) Every symbol not declared in (1) or (2) that occurs in an expression provided after `@equations` is inferred as a variable.
-    (4) Every symbol not declared in (1), (2), or (3) that occurs either as a reaction rate or stoichiometric coefficient is inferred to be a parameter.
-  E.g. in
-  ```julia
-  @reaction_network begin
-      @equations V1 + S ~ V2^2
-      (p + S + V1), S --> 0
-  end
-  ```
-  `S` is inferred as a species, `V1` and `V2` as variables, and `p` as a parameter. The previous special cases for the `@observables`, `@compounds`, and `@differentials` options still hold. Finally, the `@require_declaration` options (described in more detail below) can now be used to require everything to be explicitly declared.
-- New formula for determining whether the default differentials have been used within an `@equations` option. Now, if any expression `D(...)` is encountered (where `...` can be anything), this is inferred as usage of the default differential D. E.g. in the following equations `D` is inferred as a differential with respect to the default independent variable:
+- It is no longer recommended to install and use the full OrdinaryDiffEq library
+  to access specific ODE solvers. Instead, only install the specific
+  OrdinaryDiffEq sub-libraries that contain the desired solver. This
+  significantly reduces installation and package loading times. I.e. to use the
+  default 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 with the exception of `NonlinearProblem`s or `NonlinearSystem`s.
+  For NonlinearProblems it is safe to use `remake` if only updating `u0` values,
+  but it is not safe to update the value of the conserved constant, `Γ`. See
+  [the FAQ](https://docs.sciml.ai/Catalyst/stable/faqs/#faq_remake_nonlinprob)
+  for details.
+- 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.
+    3. Every symbol not declared in (1) or (2) that occurs in an expression
+       provided after `@equations` is inferred as a variable.
+    4. Every symbol not declared in (1), (2), or (3) that occurs either as a
+       reaction rate or stoichiometric coefficient is inferred to be a
+       parameter. E.g. in
+       ```julia
+       @reaction_network begin
+           @equations V1 + S ~ V2^2
+           (p + S + V1), S --> 0
+       end
+       ```
+       `S` is inferred as a species, `V1` and `V2` as variables, and `p` as a
+       parameter. The previous special cases for the `@observables`, `@compounds`,
+       and `@differentials` options still hold. Finally, the `@require_declaration`
+       options (described in more detail below) can now be used to require everything
+       to be explicitly declared.
+- New formula for determining whether the default differentials have been used
+  within an `@equations` option. Now, if any expression `D(...)` is encountered
+  (where `...` can be anything), this is inferred as usage of the default
+  differential D. E.g. in the following equations `D` is inferred as a
+  differential with respect to the default independent variable:
   ```julia
   @reaction_network begin
     @equations D(V) + V ~ 1
@@ -38,7 +56,9 @@
     @equations D(D(V)) ~ 1
   end
   ```
-  Please note that this cannot be used at the same time as `D` is used to represent a species, variable, or parameter (including if these are implicitly designated as such by e.g. appearing as a reaction reactant). 
+  Please note that this cannot be used at the same time as `D` is used to
+  represent a species, variable, or parameter (including if these are implicitly
+  designated as such by e.g. appearing as a reaction reactant). 
 - Array symbolics support is more consistent with ModelingToolkit v9. Parameter
   arrays are no longer scalarized by Catalyst, while species and variables
   arrays still are (as in ModelingToolkit). As such, parameter arrays should now
@@ -51,7 +71,11 @@
   *not* to do this as it has signifcant performance costs with ModelingToolkit
   v9. Note, scalarized parameter arrays passed to the two-argument
   `ReactionSystem` constructor may become unscalarized.
-- Functional (e.g. time-dependent) parameters can now be used in Catalyst models. These can e.g. be used to incorporate arbitrary time-dependent functions (as a parameter) in a model. For more details on how to use these, please read: https://docs.sciml.ai/Catalyst/stable/model_creation/functional_parameters/.
+- Functional (e.g. time-dependent) parameters can now be used in Catalyst
+  models. These can e.g. be used to incorporate arbitrary time-dependent
+  functions (as a parameter) in a model. For more details on how to use these,
+  please read:
+  https://docs.sciml.ai/Catalyst/stable/model_creation/functional_parameters/.
 - Scoped species/variables/parameters are now treated similar to the latest MTK
   releases (≥ 9.49).
 - A tutorial on making interactive plot displays using Makie has been added.
diff --git a/docs/src/faqs.md b/docs/src/faqs.md
index 25f09ca330..a91294fa56 100644
--- a/docs/src/faqs.md
+++ b/docs/src/faqs.md
@@ -271,7 +271,7 @@ so the reaction
 ```@example faq7
 using Catalyst
 rn = @reaction_network begin
-  k, X --> ∅
+    k, X --> ∅
 end
 convert(ODESystem, rn)
 ```
@@ -280,7 +280,7 @@ using any of the following non-filled arrows when declaring the reaction: `<=`,
 `⇐`, `⟽`, `=>`, `⇒`, `⟾`, `⇔`, `⟺`. This means that the reaction
 ```@example faq7
 rn = @reaction_network begin
-  k, X => ∅
+    k, X => ∅
 end
 convert(ODESystem, rn)
 ```
@@ -290,7 +290,7 @@ will be degraded at a constant rate even when very small or equal to 0).
 Note, stoichiometric coefficients are still included, i.e. the reaction
 ```@example faq7
 rn = @reaction_network begin
-  k, 2*X ⇒ ∅
+    k, 2*X ⇒ ∅
 end
 convert(ODESystem, rn)
 ```
@@ -303,7 +303,7 @@ ModelingToolkit. e.g., this is should work
 using Catalyst
 myHill(x) = 2*x^3/(x^3+1.5^3)
 rn = @reaction_network begin
-  myHill(X), ∅ --> X
+    myHill(X), ∅ --> X
 end
 ```
 In some cases, it may be necessary or desirable to register functions with
@@ -322,3 +322,137 @@ rn = @reaction_network begin
     (k1, k2), A <--> B
 end
 ```
+
+## [What to be aware of when using `remake` with conservation law elimination and NonlinearProblems?](@id faq_remake_nonlinprob)
+
+When constructing `NonlinearSystem`s or `NonlinearProblem`s with `remove_conserved = true`, i.e.
+```julia
+# for rn a ReactionSystem
+nsys = convert(NonlinearSystem, rn; remove_conserved = true)
+
+# or 
+nprob = NonlinearProblem(rn, u0, p; remove_conserved = true)
+```
+`remake` is currently unable to correctly update all `u0` values when the
+conserved constant(s), `Γ`, are updated. As an example consider the following
+```@example faq_remake
+using Catalyst, NonlinearSolve
+rn = @reaction_network begin
+    (k₁,k₂), X₁ <--> X₂
+    (k₃,k₄), X₁ + X₂ --> 2X₃
+end
+u0 = [:X₁ => 1.0, :X₂ => 2.0, :X₃ => 3.0]
+ps = [:k₁ => 0.1, :k₂ => 0.2, :k₃ => 0.3, :k₄ => 0.4]
+nlsys = convert(NonlinearSystem, rn; remove_conserved = true, conseqs_remake_warn = false)
+nlsys = complete(nlsys)
+equations(nlsys)
+```
+If we generate a `NonlinearProblem` from this system the conservation constant,
+`Γ[1]`, is automatically set to `X₁ + X₂ + X₃ = 6` and the initial values are
+those in `u0`. i.e if
+```@example faq_remake
+nlprob1 = NonlinearProblem(nlsys, u0, ps)
+```
+then
+```@example faq_remake
+nlprob1[(:X₁, :X₂, :X₃)] == (1.0, 2.0, 3.0)
+```
+and
+```@example faq_remake
+nlprob1.ps[:Γ][1] == 6.0
+```
+If we now try to change a value of `X₁`, `X₂`, or `X₃` using `remake`, the
+conserved constant will be recalculated. i.e. if
+```@example faq_remake
+nlprob2 = remake(nlprob1; u0 = [:X₂ => 3.0])
+```
+compare  
+```@example faq_remake
+println("Correct u0 is: ", (1.0, 3.0, 3.0), "\n", "remade value is: ", nlprob2[(:X₁, :X₂, :X₃)]) 
+```
+and
+```@example faq_remake
+println("Correct Γ is: ", 7.0, "\n", "remade value is: ", nlprob2.ps[:Γ][1])
+```
+However, if we try to directly change the value of `Γ` it is not always the case
+that a `u0` value will correctly update so that the conservation law is
+conserved. Consider
+```@example faq_remake
+nlprob3 = remake(nlprob1; u0 = [:X₂ => nothing], p = [:Γ => [4.0]])
+```
+Setting `[:X₂ => nothing]` for other problem types communicates that the
+`u0` value for `X₂` should be solved for. However, if we examine the values we
+find
+```@example faq_remake
+println("Correct u0 is: ", (1.0, 0.0, 3.0), "\n", "remade value is: ", nlprob3[(:X₁, :X₂, :X₃)]) 
+```
+and
+```@example faq_remake
+println("Correct Γ is: ", 4.0, "\n", "remade value is: ", nlprob3.ps[:Γ][1])
+```
+As such, the `u0` value for `X₂` has not updated, and the conservation law is
+now violated by the `u0` values, i.e,
+```@example faq_remake
+(nlprob3[:X₁] + nlprob3[:X₂] + nlprob3[:X₃]) == nlprob3.ps[:Γ][1] 
+```
+Currently, the only way to avoid this issue is to manually specify updated
+values for the `u0` components, which will ensure that `Γ` updates appropriately
+as in the first example. i.e. we manually set `X₂` to the value it should be and
+`Γ` will be updated accordingly:
+```@example faq_remake
+nlprob4 = remake(nlprob1; u0 = [:X₂ => 0.0])
+```
+so that
+```@example faq_remake
+println("Correct u0 is: ", (1.0, 0.0, 3.0), "\n", "remade value is: ", nlprob4[(:X₁, :X₂, :X₃)]) 
+```
+and
+```@example faq_remake
+println("Correct Γ is: ", 4.0, "\n", "remade value is: ", nlprob4.ps[:Γ][1])
+```
+
+Finally, we note there is one extra consideration to take into account if using
+`structural_simplify`. In this case one of `X₁`, `X₂`, or `X₃` will be moved to
+being an observed. It will then always correspond to the updated value if one
+tries to manually change `Γ`. Let's see what happens here directly
+```@example faq_remake
+nlsys = convert(NonlinearSystem, rn; remove_conserved = true, conseqs_remake_warn = false)
+nlsys = structural_simplify(nlsys)
+nlprob1 = NonlinearProblem(nlsys, u0, ps)
+```
+We can now try to change just `Γ` and implicitly the observed variable that was
+removed will be assumed to have changed its initial value to compensate for it.
+Let's confirm this. First we find the observed variable that was elminated.
+```@example faq_remake
+obs_unknown = only(observed(nlsys)).lhs
+```
+We can figure out its index in `u0` via
+```@example faq_remake
+obs_symbol = ModelingToolkit.getname(obs_unknown)
+obsidx = findfirst(p -> p[1] == obs_symbol, u0)
+```
+Let's now remake 
+```@example faq_remake
+nlprob2 = remake(nlprob1; u0 = [obs_unknown => nothing], p = [:Γ => [8.0]])
+```
+Here we indicate that the observed variable should be treated as unspecified
+during initialization. Since the observed variable is not considered an unknown,
+everything now works, with the observed variable's assumed initial value
+adjusted to allow `Γ = 8`:
+```@example faq_remake
+correct_u0 = last.(u0)
+correct_u0[obsidx] = 8 - sum(correct_u0) + correct_u0[obsidx]
+println("Correct u0 is: ", (1.0, 2.0, 5.0), "\n", "remade value is: ", nlprob2[(:X₁, :X₂, :X₃)]) 
+```
+and `Γ` becomes
+```@example faq_remake
+println("Correct Γ is: ", 8.0, "\n", "remade value is: ", nlprob2.ps[:Γ][1])
+```
+Unfortunately, as with our first example, trying to enforce that a
+non-eliminated species should have its initial value updated instead of the
+observed species will not work.
+
+*Summary:* it is not recommended to directly update `Γ` via `remake`, but to
+instead update values of the initial guesses in `u0` to obtain a desired `Γ`. At
+this time the behavior when updating `Γ` can result in `u0` values that do not
+satisfy the conservation law defined by `Γ` as illustrated above. 
\ No newline at end of file
diff --git a/docs/src/model_creation/conservation_laws.md b/docs/src/model_creation/conservation_laws.md
index 9adebaa5df..5bb77eb649 100644
--- a/docs/src/model_creation/conservation_laws.md
+++ b/docs/src/model_creation/conservation_laws.md
@@ -133,7 +133,7 @@ Generally, for a conservation law where we have:
 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`).
+    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`). Moreover, care is needed when using `remake` with `NonlinearProblem`s for which conserved equations have been removed. See [the FAQ](https://docs.sciml.ai/Catalyst/stable/faqs/#faq_remake_nonlinprob) for details on what is supported and what may lead to `u0` values that do not satisfy the conservation law in the special case of `NonlinearProblem`s.
 
 ### [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:
diff --git a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
index 6528a93041..9b207bc2e8 100644
--- a/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
+++ b/ext/CatalystBifurcationKitExtension/bifurcation_kit_extension.jl
@@ -37,6 +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)
+    nsys = convert(NonlinearSystem, rs; defaults, remove_conserved = true, conseqs_remake_warn = false)
     return complete(nsys)
 end
diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
index 4e70b619f7..2b353d420f 100644
--- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
+++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl
@@ -1,16 +1,16 @@
 ### Homotopy Continuation Based Steady State Finding ###
 
 """
-    hc_steady_states(rs::ReactionSystem, ps; filter_negative=true, neg_thres=-1e-20, u0=typeof(ps)(), kwargs...)
+    hc_steady_states(rs::ReactionSystem, ps; filter_negative = true, neg_thres = -1e-16, u0 = typeof(ps)(), kwargs...)
 
 Uses homotopy continuation via HomotopyContinuation.jl to find the steady states of the ODE system corresponding to the provided reaction system.
 
 Arguments:
 - `rs::ReactionSystem`: The reaction system for which we want to find the steady states.
 - `ps`: The parameter values for which we want to find the steady states.
-- `filter_negative=true`: If set to true, solutions with any species concentration <neg_thres is removed from the output.
-- `neg_thres=-1e-20`: Determine the minimum values for which a species concentration is to be considered non-negative. Species concentrations ``> neg_thres`` but `< 0.0` are set to `0.0`.
-- `u0=nothing`: Initial conditions for which we want to find the steady states. For systems with conservation laws this are required to compute conserved quantities. Initial conditions are not required for all species, only those involved in conserved quantities (if this set is unknown, it is recommended to provide initial conditions for all species).
+- `filter_negative = true`: If set to true, solutions with any species concentration <neg_thres is removed from the output.
+- `neg_thres = -1e-15`: Determine the minimum values for which a species concentration is to be considered non-negative. Species concentrations ``> neg_thres`` but `< 0.0` are set to `0.0`.
+- `u0 = nothing`: Initial conditions for which we want to find the steady states. For systems with conservation laws this are required to compute conserved quantities. Initial conditions are not required for all species, only those involved in conserved quantities (if this set is unknown, it is recommended to provide initial conditions for all species).
 - `kwargs...`: any additional arguments (like `show_progress= true`) are passed into HomotopyContinuation.jl's `solve` call.
 
 Examples
@@ -36,7 +36,7 @@ Notes:
 ```
   """
 function Catalyst.hc_steady_states(rs::ReactionSystem, ps; filter_negative = true,
-        neg_thres = -1e-20, u0 = [], kwargs...)
+        neg_thres = -1e-15, u0 = [], kwargs...)
     if !isautonomous(rs)
         error("Attempting to compute steady state for a non-autonomous system (e.g. where some rate depend on $(get_iv(rs))). This is not possible.")
     end
@@ -51,7 +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))
+    ns = complete(convert(NonlinearSystem, rs; remove_conserved = true, conseqs_remake_warn = false))
     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)
@@ -164,7 +164,7 @@ function reorder_sols!(sols, ss_poly, rs::ReactionSystem)
 end
 
 # Filters away solutions with negative species concentrations (and for neg_thres < val < 0.0, sets val=0.0).
-function filter_negative_f(sols; neg_thres = -1e-20)
+function filter_negative_f(sols; neg_thres = -1e-15)
     for sol in sols, idx in 1:length(sol)
         (neg_thres < sol[idx] < 0) && (sol[idx] = 0)
     end
diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl
index 0b0d231850..32fc1710dd 100644
--- a/src/reactionsystem_conversions.jl
+++ b/src/reactionsystem_conversions.jl
@@ -381,12 +381,16 @@ end
 
 # merge constraint components with the ReactionSystem components
 # also handles removing BC and constant species
-function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved = false)
+function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved = false,
+        treat_conserved_as_eqs = false)
     # if there are BC species, put them after the independent species
     rssts = get_unknowns(rs)
     sts = any(isbc, rssts) ? vcat(ists, filter(isbc, rssts)) : ists
     ps = get_ps(rs)
-
+    initeqs = Equation[]
+    defs = MT.defaults(rs)
+    obs = MT.observed(rs)
+    
     # make dependent species observables and add conservation constants as parameters
     if remove_conserved
         nps = get_networkproperties(rs)
@@ -394,14 +398,23 @@ function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved
         # add the conservation constants as parameters and set their values
         ps = copy(ps)
         push!(ps, nps.conservedconst)
-        defs = copy(MT.defaults(rs))
 
-        # add the dependent species as observed
-        obs = copy(MT.observed(rs))
-        append!(obs, nps.conservedeqs)
-    else
-        defs = MT.defaults(rs)
-        obs = MT.observed(rs)
+        if treat_conserved_as_eqs
+            # add back previously removed dependent species
+            sts = union(sts, nps.depspecs)       
+
+            # treat conserved eqs as normal eqs
+            append!(eqs, conservedequations(rs))
+
+            # add initialization equations for conserved parameters
+            initialmap = Dict(u => Initial(u) for u in species(rs))
+            conseqs = conservationlaw_constants(rs)
+            initeqs = [Symbolics.substitute(conseq, initialmap) for conseq in conseqs]        
+        else
+            # add the dependent species as observed
+            obs = copy(obs)
+            append!(obs, conservedequations(rs))
+        end
     end
 
     ceqs = Equation[eq for eq in get_eqs(rs) if eq isa Equation]
@@ -419,9 +432,10 @@ function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved
         append!(eqs, ceqs)
     end
 
-    eqs, sts, ps, obs, defs
+    eqs, sts, ps, obs, defs, initeqs
 end
 
+
 # used by flattened systems that don't support constraint equations currently
 function error_if_constraints(::Type{T}, sys::ReactionSystem) where {T <: MT.AbstractSystem}
     any(eq -> eq isa Equation, get_eqs(sys)) &&
@@ -495,6 +509,23 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs)
         kwargs...)
 end
 
+
+const NONLIN_PROB_REMAKE_WARNING = """
+    Note, when constructing `NonlinearSystem`s with `remove_conserved = true`, possibly via \
+    calling `NonlinearProblem`, `remake(::NonlinearProblem)` has some \
+    limitations. If in `remake` the value of the conserved constant is \
+    explicitly updated, it is not possible to have one variable's `u0` value \
+    auto-update to preserve the conservation law. See the [FAQ \
+    entry](https://docs.sciml.ai/Catalyst/stable/faqs/#faq_remake_nonlinprob) for more \
+    details. This warning can be disabled by passing `conseqs_remake_warn = false`."""
+
+function is_autonomous_error(iv) 
+    return """
+    Attempting to convert a non-autonomous `ReactionSystem` (e.g. where some rate depends \
+    on $(iv)) to a `NonlinearSystem`. This is not possible. if you are intending to \
+    compute system steady states, consider creating and solving a `SteadyStateProblem."""
+end
+
 """
 ```julia
 Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
@@ -503,34 +534,40 @@ Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.NonlinearSystem`.
 
 Keyword args and default values:
-- `combinatoric_ratelaws=true` uses factorial scaling factors in calculating the rate law,
+- `combinatoric_ratelaws = true` uses factorial scaling factors in calculating the rate law,
   i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
   `combinatoric_ratelaws=false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
   ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
   itself defaults to true).
-- `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 = false`, if set to `true` will calculate conservation laws of the
+  underlying set of reactions (ignoring coupled ODE or algebraic equations). For each
+  conservation law one steady-state equation is eliminated, and replaced with the
+  conservation law. This ensures a non-singular Jacobian.
+- `conseqs_remake_warn = true`, set to false to disable warning about `remake` and
+  conservation laws. See the [FAQ
+  entry](https://docs.sciml.ai/Catalyst/stable/faqs/#faq_remake_nonlinprob) for more
+  details.
 """
 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,
-        default_u0 = Dict(), default_p = Dict(), defaults = _merge(Dict(default_u0), Dict(default_p)),
+        remove_conserved = false, conseqs_remake_warn = true, checks = false,
+        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)
-    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
+    remove_conserved && conseqs_remake_warn && (@warn NONLIN_PROB_REMAKE_WARNING)
+    isautonomous(rs) || error(is_autonomous_error(get_iv(rs)))
 
     # Generates system equations.
     fullrs = Catalyst.flatten(rs)
     remove_conserved && conservationlaws(fullrs)
     ists, ispcs = get_indep_sts(fullrs, remove_conserved)
     eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved,
-        as_odes = false, include_zero_odes)
-    eqs, us, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
+        as_odes = false, include_zero_odes = false)
+    eqs, us, ps, obs, defs, initeqs = addconstraints!(eqs, fullrs, ists, ispcs; 
+        remove_conserved, treat_conserved_as_eqs = true)
 
     # Throws a warning if there are differential equations in non-standard format.
     # Next, sets all differential terms to `0`.
@@ -539,7 +576,7 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
 
     NonlinearSystem(eqs, us, ps;
         name,
-        observed = obs,
+        observed = obs, initialization_eqs = initeqs,
         defaults = _merge(defaults, defs),
         checks,
         kwargs...)
@@ -699,16 +736,43 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan,
     return ODEProblem(osys, u0, tspan, p, args...; check_length, kwargs...)
 end
 
-# NonlinearProblem from AbstractReactionNetwork
-function DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0,
+"""
+```julia
+DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0,
         p = DiffEqBase.NullParameters(), args...;
-        name = nameof(rs), include_zero_odes = true,
-        combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
+        name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         remove_conserved = false, checks = false, check_length = false, 
-        all_differentials_permitted = false, kwargs...)
-    nlsys = convert(NonlinearSystem, rs; name, combinatoric_ratelaws, include_zero_odes,
-        checks, all_differentials_permitted, remove_conserved)
-    nlsys = complete(nlsys)
+        structural_simplify = remove_conserved, all_differentials_permitted = false, 
+        kwargs...)
+```
+
+Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.NonlinearSystem`.
+
+Keyword args and default values:
+- `combinatoric_ratelaws=true` uses factorial scaling factors in calculating the rate law,
+  i.e. for `2S -> 0` at rate `k` the ratelaw would be `k*S^2/2!`. Set
+  `combinatoric_ratelaws=false` for a ratelaw of `k*S^2`, i.e. the scaling factor is
+  ignored. Defaults to the value given when the `ReactionSystem` was constructed (which
+  itself defaults to true).
+- `remove_conserved=false`, if set to `true` will calculate conservation laws of the
+  underlying set of reactions (ignoring coupled ODE or algebraic equations). For each
+  conservation law one steady-state equation is eliminated, and replaced with the
+  conservation law. This ensures a non-singular Jacobian. When using this option, it is
+  recommended to call `ModelingToolkit.structural_simplify` on the converted system to then
+  eliminate the conservation laws from the system equations.
+- `conseqs_remake_warn = true`, set to false to disable warning about `remake` and
+  conservation laws. See the [FAQ
+  entry](https://docs.sciml.ai/Catalyst/stable/faqs/#faq_remake_nonlinprob) for more
+  details.
+"""
+function DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0,
+        p = DiffEqBase.NullParameters(), args...;
+        name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
+        remove_conserved = false, conseqs_remake_warn = true, checks = false, check_length = false, 
+        structural_simplify = false, all_differentials_permitted = false, kwargs...)
+    nlsys = convert(NonlinearSystem, rs; name, combinatoric_ratelaws, checks, 
+        all_differentials_permitted, remove_conserved, conseqs_remake_warn)
+    nlsys = structural_simplify ? MT.structural_simplify(nlsys) : complete(nlsys)
     return NonlinearProblem(nlsys, u0, p, args...; check_length,
         kwargs...)
 end
diff --git a/test/compositional_modelling/component_based_model_creation.jl b/test/compositional_modelling/component_based_model_creation.jl
index 50b1795321..bd96762f3e 100644
--- a/test/compositional_modelling/component_based_model_creation.jl
+++ b/test/compositional_modelling/component_based_model_creation.jl
@@ -99,7 +99,7 @@ let
     @named nsys = NonlinearSystem(connections, [], [])
     @named ssrepressilator = ReactionSystem(t; systems = [nsys, sys₁, sys₂, sys₃])
     ssrepressilator = complete(ssrepressilator)
-    @named nlrepressilator = convert(NonlinearSystem, ssrepressilator, include_zero_odes = false)
+    @named nlrepressilator = convert(NonlinearSystem, ssrepressilator)
     sys2 = structural_simplify(nlrepressilator)
     @test length(equations(sys2)) <= 6
     nlprob = NonlinearProblem(sys2, u₀_nl, pvals)
@@ -112,7 +112,7 @@ let
     # Flattening.
     fsys = Catalyst.flatten(ssrepressilator)
     fsys = complete(fsys)
-    @named nlrepressilator = convert(NonlinearSystem, fsys, include_zero_odes = false)
+    @named nlrepressilator = convert(NonlinearSystem, fsys)
     sys2 = structural_simplify(nlrepressilator)
     @test length(equations(sys2)) <= 6
     nlprob = NonlinearProblem(sys2, u₀_nl, pvals)
@@ -130,7 +130,7 @@ let
                                 [])
     @named repressilator2 = ReactionSystem(connections, t; systems = [sys₁, sys₂, sys₃])
     repressilator2 = complete(repressilator2)
-    @named nlrepressilator = convert(NonlinearSystem, repressilator2, include_zero_odes = false)
+    @named nlrepressilator = convert(NonlinearSystem, repressilator2)
     sys2 = structural_simplify(nlrepressilator)
     @test length(equations(sys2)) <= 6
     nlprob = NonlinearProblem(sys2, u₀_nl, pvals)
@@ -249,7 +249,7 @@ let
     repressilator2 = Catalyst.flatten(repressilator2)
     repressilator2 = extend(csys, repressilator2)
     repressilator2 = complete(repressilator2)
-    @named nlrepressilator = convert(NonlinearSystem, repressilator2, include_zero_odes = false)
+    @named nlrepressilator = convert(NonlinearSystem, repressilator2)
     sys2 = structural_simplify(nlrepressilator)
     @test length(equations(sys2)) <= 6
     nlprob = NonlinearProblem(sys2, u₀_nl, pvals)
diff --git a/test/network_analysis/conservation_laws.jl b/test/network_analysis/conservation_laws.jl
index 7a0177fe22..78961ea84e 100644
--- a/test/network_analysis/conservation_laws.jl
+++ b/test/network_analysis/conservation_laws.jl
@@ -124,22 +124,32 @@ let
     @test osol1[sps] ≈ osol2[sps] ≈ osol3[sps]
 
     # 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))
-        nprob1 = NonlinearProblem{true}(nsys, u0, p; guesses = [nsys.Γ => [0.0]])
-        nprob2 = NonlinearProblem(rn, u0, p)
-        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)
-        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)
-        sssol1 = solve(ssprob1, DynamicSS(Tsit5()); abstol = 1e-8, reltol = 1e-8)
-        sssol2 = solve(ssprob2, DynamicSS(Tsit5()); abstol = 1e-8, reltol = 1e-8)
-        sssol3 = solve(ssprob3, DynamicSS(Tsit5()); abstol = 1e-8, reltol = 1e-8)
-        @test nsol1[sps] ≈ nsol3[sps] ≈ sssol1[sps] ≈ sssol2[sps] ≈ sssol3[sps]
-    end
+    nsys = convert(NonlinearSystem, rn; remove_conserved = true, conseqs_remake_warn = false)
+    nsys_ss = structural_simplify(nsys)
+    nsys = complete(nsys)
+    nprob1 = NonlinearProblem{true}(nsys, u0, p)
+    nprob1b = NonlinearProblem{true}(nsys_ss, u0, p)
+    nprob2 = NonlinearProblem(rn, u0, p; remove_conserved = true, conseqs_remake_warn = false)
+    nprob2b = NonlinearProblem(rn, u0, p; remove_conserved = true, conseqs_remake_warn = false, 
+        structural_simplify = true)
+    nsol1 = solve(nprob1, NewtonRaphson(); abstol = 1e-8)
+    nsol1b = solve(nprob1b, NewtonRaphson(); abstol = 1e-8)
+    nsol2 = solve(nprob2, NewtonRaphson(); abstol = 1e-8)
+    nsol2b = solve(nprob2b, NewtonRaphson(); abstol = 1e-8)
+    # Nonlinear problems cannot find steady states properly without removing conserved species.
+
+    ssprob1 = SteadyStateProblem{true}(osys, u0, p)
+    ssprob2 = SteadyStateProblem(rn, u0, p)
+    ssprob3 = SteadyStateProblem(rn, u0, p; remove_conserved = true)
+    sssol1 = solve(ssprob1, DynamicSS(Tsit5()); abstol = 1e-8, reltol = 1e-8)
+    sssol2 = solve(ssprob2, DynamicSS(Tsit5()); abstol = 1e-8, reltol = 1e-8)
+    sssol3 = solve(ssprob3, DynamicSS(Tsit5()); abstol = 1e-8, reltol = 1e-8)
+    @test nsol1[sps] ≈ nsol1b[sps] 
+    @test nsol1[sps] ≈ nsol2[sps] 
+    @test nsol1[sps] ≈ nsol2b[sps] 
+    @test nsol1[sps] ≈ sssol1[sps]
+    @test nsol1[sps] ≈ sssol2[sps] 
+    @test nsol1[sps] ≈ sssol3[sps]
 
     # 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,
@@ -216,11 +226,12 @@ let
         (k1,k2), X1 <--> X2
     end
     osys = complete(convert(ODESystem, rn; remove_conserved = true))
-    u0 = [osys.X1 => 1.0, osys.X2 => 1.0]
+    u0_1 = [osys.X1 => 1.0, osys.X2 => 1.0]
+    u0_2 = [osys.X1 => 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]]
-    oprob1 = ODEProblem(osys, u0, 10.0, ps_1)
-    oprob2 = ODEProblem(osys, u0, 10.0, ps_2)
+    oprob1 = ODEProblem(osys, u0_1, 10.0, ps_1)
+    oprob2 = ODEProblem(osys, u0_2, 10.0, ps_2)
 
     # Checks that the solutions generates the correct conserved quantities.
     sol1 = solve(oprob1; saveat = 1.0)
@@ -259,8 +270,6 @@ let
     @test oprob.ps[k2] == 0.2
     @test oprob.ps[Γ[1]] == 3.0
 
-    # Currently, any kind of updating of species or the conservation parameter(s) is not possible.
-
     # Update problem parameters using `remake`.
     oprob_new = remake(oprob; p = [k1 => 0.3, k2 => 0.4])
     @test oprob_new.ps[k1] == 0.3
@@ -279,11 +288,59 @@ let
     @test integrator.ps[k2] == 0.6
 end
 
+### Jacobian Tests ###
+
+# Checks that conservation law elimination generates a system with a non-singular Jacobian.
+# Does this for different system types (ODE, SDE, and Nonlinear).
+# Checks singularity by checking whether the Jacobian have high enough a condition number 
+# (when evaluated at the steady state).
+let
+    # Creates the model (contains both conserved and non-conserved species).
+    rn = @reaction_network begin
+        (p,d), 0 <--> X
+        d, XY --> Y
+        (k1,k2), X + Y <--> XY
+        (k3,k4), 2Y <--> Y2
+    end
+
+    # Finds a steady state at which we will compute the Jacobian.
+    u0_ss_init = [:X => 0.0, :Y => 4.0, :XY => 0.0, :Y2 => 0.0]
+    ps = [:p => 1.0, :d => 0.5, :k1 => 1.0, :k2 => 2.0, :k3 => 3.0, :k4 => 4.0]
+    ssprob = SteadyStateProblem(rn, u0_ss_init, ps)
+    ss = solve(ssprob, DynamicSS(Vern7()))
+    ss = [sp => ss[sp] for sp in unknowns(rn)]
+
+    # Creates a function which evaluates whether the Jacobian is singular.
+    # Singularity means infinite condition number (here it is about 1e17).
+    function is_singular(prob; infthres = 1e12)
+        J = zeros(length(prob.u0), length(prob.u0))
+        ModelingToolkit.is_time_dependent(prob) ? prob.f.jac(J, prob.u0, prob.p, 0.0) : prob.f.jac(J, prob.u0, prob.p)
+        return cond(J) > infthres
+    end
+
+    # Creates problems (with Jacobian and W/O conservation law elimination) which singularity we wish to test.
+    oprob = ODEProblem(rn, ss, 1.0, ps; jac = true)
+    oprob_rc = ODEProblem(rn, ss, 1.0, ps; jac = true, remove_conserved = true)
+    sprob = SDEProblem(rn, ss, 1.0, ps; jac = true)
+    sprob_rc = SDEProblem(rn, ss, 1.0, ps; jac = true, remove_conserved = true)
+    nlprob = NonlinearProblem(rn, ss, ps; jac = true)
+    nlprob_rc = NonlinearProblem(rn, ss, ps; jac = true, remove_conserved = true, conseqs_remake_warn = false)
+
+    # Checks that removing conservation laws generates non-singular Jacobian (and else that it is singular).
+    @test is_singular(oprob) == true
+    @test !is_singular(oprob_rc) == true
+    @test is_singular(sprob) == true
+    @test !is_singular(sprob_rc) == true
+    @test is_singular(nlprob) == true
+    @test !is_singular(nlprob_rc) == true
+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).
+# the values of relevant quantities are correct after each step.
+# Generally, if `Γ` has not been explicitly updated, it will be updated to acomodate new species
+# values. If it has been explicitly updated, the corresponding eliminated quantity will have its
+# value updated to acomodate new Γ/species values (also, the eliminated species's value can not longer be changed).
+# Also checks that quantities are correctly updated in integrators and solutions derived from problems.
 let
     # Prepares the problem inputs and computes the conservation equation.
     rn = @reaction_network begin
@@ -298,7 +355,9 @@ let
     # Loops through the tests for different problem types.
     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]
+    # nlprob = NonlinearProblem(rn, u0, ps; remove_conserved = true, conseqs_remake_warn = false)
+    @test_broken false # Cases of `remake` does not work for `NonlinearSystem`s with `remove_conserved = true`.
+    for (prob_old, solver) in zip([oprob, sprob], [Tsit5(), ImplicitEM()])
         # 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.
@@ -320,6 +379,20 @@ let
             @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
+
+            # Checks that integrator and solutions have identical content to problem.
+            integrator = init(prob_new, solver)
+            sol = solve(prob_new, solver; maxiters = 1, verbose = false)
+            @test prob_new.ps[:Γ][1] == integrator.ps[:Γ][1] == sol.ps[:Γ][1]
+            if ModelingToolkit.is_time_dependent(prob_new)
+                @test prob_new[:X1] == integrator[:X1] == sol[:X1][1]
+                @test prob_new[:X2] == integrator[:X2] == sol[:X2][1]
+                @test prob_new[:X3] == integrator[:X3] == sol[:X3][1]
+            else
+                @test prob_new[:X1] == integrator[:X1]
+                @test prob_new[:X2] == integrator[:X2]
+                @test prob_new[:X3] == integrator[:X3]
+            end
         end
 
         # Similarly, but now also updates the conservation constant. Here, once Γ has been updated:
@@ -330,7 +403,7 @@ let
         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)
+            prob_new = remake(prob_old; p = [:Γ => [Γ_new]])
             @test prob_old[:X1] ≈ prob_new[:X1]
             @test prob_old[:X2] ≈ prob_new[:X2]
             @test Γ_new ≈ prob_new.ps[:Γ][1]
@@ -349,12 +422,26 @@ let
 
             # 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)
+            prob_new = remake(prob_old; u0 = [:X3 => X3_new])
             @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
+
+            # Checks that integrator and solutions have identical content to problem.
+            integrator = init(prob_new, solver)
+            sol = solve(prob_new, solver; maxiters = 1, verbose = false)
+            @test prob_new.ps[:Γ][1] == integrator.ps[:Γ][1] == sol.ps[:Γ][1]
+            if ModelingToolkit.is_time_dependent(prob_new)
+                @test prob_new[:X1] == integrator[:X1] == sol[:X1][1]
+                @test prob_new[:X2] == integrator[:X2] == sol[:X2][1]
+                @test prob_new[:X3] == integrator[:X3] == sol[:X3][1]
+            else
+                @test prob_new[:X1] == integrator[:X1]
+                @test prob_new[:X2] == integrator[:X2]
+                @test prob_new[:X3] == integrator[:X3]
+            end
         end
     end
 end
@@ -408,3 +495,107 @@ let
     # @test sol[X[1]][end] ≈ 8.0
     # @test sol[X[2]][end] ≈ 4.0
 end
+
+# Check conservation law elimination warnings (and the disabling of these) for `NonlinearSystem`s
+# and `NonlinearProblem`s.
+let
+    # Create models.
+    rn = @reaction_network begin
+        (k1,k2), X1 <--> X2
+        (k3,k4), X1 + X2 <--> 2X3
+    end
+    u0 = [:X1 => 1.0, :X2 => 2.0, :X3 => 3.0]
+    ps = [:k1 => 0.1, :k2 => 0.2, :k3 => 0.3, :k4 => 0.4]
+    
+    # Checks that the warning si given and can be supressed for the variosu cases.
+    @test_nowarn convert(NonlinearSystem, rn; remove_conserved = true, conseqs_remake_warn = false)
+    @test_logs (:warn, r"Note, when constructing*") convert(NonlinearSystem, rn; remove_conserved = true, conseqs_remake_warn = true)
+    @test_nowarn NonlinearProblem(rn, u0, ps; remove_conserved = true, conseqs_remake_warn = false)
+    @test_logs (:warn, Catalyst.NONLIN_PROB_REMAKE_WARNING) NonlinearProblem(rn, u0, ps; remove_conserved = true, conseqs_remake_warn = true)
+
+    # @variables X1 X2 X3
+    # @parameters k1 k2 k3 k4 Γ[1:1] = missing [guess = ones(1)]
+    # eqs = [
+    #     0 ~ -k1*X1 + k2*X2 - k3*X1*X2 + (1//2)*k4*((-X1 - X2 + Γ[1])^2),
+    #     0 ~ k1*X1 - k2*X2 - k3*X1*X2 + (1//2)*k4*((-X1 - X2 + Γ[1])^2),
+    #     0 ~ -X1 - X2 - X3 + Γ[1]
+    # ]
+    # initeqs = [Γ[1] ~ Initial(X1) + Initial(X3) + Initial(X2)]
+    # @named nlsys = NonlinearSystem(eqs, [X1, X2, X3], [k1, k2, k3, k4, Γ]; 
+    #     initialization_eqs = initeqs)
+
+    
+    # WITHOUT structural_simplify
+    nlsys = convert(NonlinearSystem, rn; remove_conserved = true, 
+        conseqs_remake_warn = false)
+    nlsys1 = complete(nlsys)
+    nlprob1 = NonlinearProblem(nlsys1, u0, ps)
+
+    @test nlprob1[:X1] == 1.0
+    @test nlprob1[:X2] == 2.0
+    @test nlprob1[:X3] == 3.0
+    @test nlprob1.ps[:Γ][1] == 6.0
+    integ1 = init(nlprob1, NewtonRaphson())
+    @test integ1[:X1] == 1.0
+    @test integ1[:X2] == 2.0
+    @test integ1[:X3] == 3.0
+    @test integ1.ps[:Γ][1] == 6.0
+
+    nlprob1b = remake(nlprob1; u0 = [:X3 => nothing], p = [:Γ => [4.0]])
+    @test nlprob1b[:X1] == 1.0   
+    @test nlprob1b[:X2] == 2.0   
+    @test_broken nlprob1b[:X3] == 1.0
+    @test nlprob1b.ps[:Γ][1] == 4.0
+    integ1 = init(nlprob1b, NewtonRaphson())
+    @test integ1[:X1] == 1.0
+    @test integ1[:X2] == 2.0
+    @test_broken integ1[:X3] == 1.0
+    @test integ1.ps[:Γ][1] == 4.0
+
+    nlprob1c = remake(nlprob1; u0 = [:X2 => nothing], p = [:Γ => [4.0]])
+    @test nlprob1c[:X1] == 1.0
+    @test_broken nlprob1c[:X2] == 0.0
+    @test nlprob1c[:X3] == 3.0
+    @test nlprob1c.ps[:Γ][1] == 4.0
+    integ1 = init(nlprob1b, NewtonRaphson())
+    @test integ1[:X1] == 1.0
+    @test_broken integ1[:X2] == 0.0
+    @test integ1[:X3] == 3.0
+    @test integ1.ps[:Γ][1] == 4.0
+    
+    # WITH structural_simplify
+    nlsys2 = structural_simplify(nlsys)
+    nlprob2 = NonlinearProblem(nlsys2, u0, ps)
+
+    @test nlprob2[:X1] == 1.0
+    @test nlprob2[:X2] == 2.0
+    @test nlprob2[:X3] == 3.0
+    @test nlprob2.ps[:Γ][1] == 6.0
+    integ2 = init(nlprob2, NewtonRaphson())
+    @test integ2[:X1] == 1.0
+    @test integ2[:X2] == 2.0
+    @test integ2[:X3] == 3.0
+    @test integ2.ps[:Γ][1] == 6.0
+
+    nlprob2b = remake(nlprob2; u0 = [:X3 => nothing], p = [:Γ => [4.0]])
+    @test nlprob2b[:X1] == 1.0   
+    @test nlprob2b[:X2] == 2.0   
+    @test nlprob2b[:X3] == 1.0   
+    @test nlprob2b.ps[:Γ][1] == 4.0 
+    integ2 = init(nlprob2b, NewtonRaphson())
+    @test integ2[:X1] == 1.0  
+    @test integ2[:X2] == 2.0  
+    @test integ2[:X3] == 1.0  
+    @test integ2.ps[:Γ][1] == 4.0 
+
+    nlprob2c = remake(nlprob2; u0 = [:X2 => nothing], p = [:Γ => [4.0]])
+    @test nlprob2c[:X1] == 1.0   
+    @test_broken nlprob2c[:X2] == 0.0 
+    @test_broken nlprob2c[:X3] == 3.0
+    @test nlprob2c.ps[:Γ][1] == 4.0 
+    integ2 = init(nlprob2c, NewtonRaphson())
+    @test integ2[:X1] == 1.0  
+    @test_broken integ2[:X2] == 0.0 
+    @test_broken integ2[:X3] == 3.0
+    @test integ2.ps[:Γ][1] == 4.0 
+end
\ No newline at end of file
diff --git a/test/reactionsystem_core/coupled_equation_crn_systems.jl b/test/reactionsystem_core/coupled_equation_crn_systems.jl
index ea6b11ac3f..88140bac2e 100644
--- a/test/reactionsystem_core/coupled_equation_crn_systems.jl
+++ b/test/reactionsystem_core/coupled_equation_crn_systems.jl
@@ -167,7 +167,7 @@ let
     u0 = [X => 0.1, A => 16.1/2]
     nlprob = NonlinearProblem(coupled_rs, u0, ps, structural_simplify = true)
     nlsol = solve(nlprob)
-    @test nlsol ≈ [2.0, 3.0]
+    @test nlsol[[X,A]] ≈ [2.0, 3.0]
 
     # Checks that the correct steady state is found through SteadyStateProblem.
     u0 = [X => 0.1]
@@ -380,7 +380,7 @@ end
     ps = [A1 => 0.1, B1 => 1.0, C1 => 10.0]
 
     # Create ODE structures.
-    oprob = ODEProblem(coupled_rs, u0, tspan, ps; structural_simplify = true)
+    oprob = ODEProblem(coupled_rs, u0, tspan, ps; structural_simplify = true, warn_initialize_determined = false)
     oint = init(oprob, Tsit5())
     osol = solve(oprob, Tsit5())
 
@@ -485,7 +485,7 @@ let # SDEs are currently broken with structural simplify (https://github.com/Sci
     @test_throws Exception SDEProblem(coupled_rs, u0, tspan, ps)
 
     # Checks the algebraic equation holds.
-    sprob = SDEProblem(coupled_rs, u0, tspan, ps; structural_simplify = true)
+    sprob = SDEProblem(coupled_rs, u0, tspan, ps; structural_simplify = true, warn_initialize_determined = false)
     ssol = solve(sprob, ImplicitEM())
     @test (2 .+ ps[k1] * ssol[:A]) ≈ (3 .+ ps[k2] * ssol[:X])
 end
@@ -610,7 +610,7 @@ let
         coupled_sir_ordered = complete(coupled_sir_ordered)
 
         # Checks that ODE an simulation of the system achieves the correct steady state.
-        oprob_ordered = ODEProblem(coupled_sir_ordered, u0, tspan, ps; structural_simplify = true)
+        oprob_ordered = ODEProblem(coupled_sir_ordered, u0, tspan, ps; structural_simplify = true, warn_initialize_determined = false)
         solve(oprob_ordered, Vern7(); abstol = 1e-8, reltol = 1e-8, saveat = 1.0)
     end
 
@@ -630,7 +630,7 @@ let
         coupled_sir_messy = complete(coupled_sir_messy)
 
         # Checks that ODE an simulation of the system achieves the correct steady state.
-        oprob_messy = ODEProblem(coupled_sir_messy, u0, tspan, ps; structural_simplify = true)
+        oprob_messy = ODEProblem(coupled_sir_messy, u0, tspan, ps; structural_simplify = true, warn_initialize_determined = false)
         solve(oprob_messy, Vern7(); abstol = 1e-8, reltol = 1e-8, saveat = 1.0)
     end
 
@@ -698,8 +698,8 @@ let
     # Test most likely redundant, but seem useful to have one test like this to be sure.
     u0 = [X1 => 0.1, X2 => 0.2, X3 => 0.2, X_tot => 0.6, N => 10.0, X_conc => 10.0]
     ps = [p => 1.0, k1 => 1.2, k2 => 1.5, d => 2.0, v => 0.2, n => 0.5, x_scale => 2.0]
-    oprob_prog = ODEProblem(rs_prog, u0, (0.0, 10.0), ps; structural_simplify = true)
-    oprob_dsl = ODEProblem(rs_dsl, u0, (0.0, 10.0), ps; structural_simplify = true)
+    oprob_prog = ODEProblem(rs_prog, u0, (0.0, 10.0), ps; structural_simplify = true, warn_initialize_determined = false)
+    oprob_dsl = ODEProblem(rs_dsl, u0, (0.0, 10.0), ps; structural_simplify = true, warn_initialize_determined = false)
     @test solve(oprob_prog, Rosenbrock23()) == solve(oprob_dsl, Rosenbrock23())
 end
 
@@ -998,7 +998,7 @@ let
     rs = complete(rs)
     u0 = [S1 => 1.0, S2 => 2.0]
     ps = [p1 => 2.0, p2 => 3.0]
-    @test_throws Exception ODEProblem(rs, u0, (0.0, 1.0), ps; structural_simplify = true)
+    @test_throws Exception ODEProblem(rs, u0, (0.0, 1.0), ps; structural_simplify = true, warn_initialize_determined = false)
 
     # Coupled system overconstrained due to additional algebraic equations (with variables).
     eqs = [
diff --git a/test/simulation_and_solving/simulate_ODEs.jl b/test/simulation_and_solving/simulate_ODEs.jl
index 491c45c953..e12cac10a2 100644
--- a/test/simulation_and_solving/simulate_ODEs.jl
+++ b/test/simulation_and_solving/simulate_ODEs.jl
@@ -1,7 +1,7 @@
 ### Prepares Tests ###
 
 # Fetch packages.
-using Catalyst, OrdinaryDiffEqVerner, OrdinaryDiffEqRosenbrock, Test
+using Catalyst, OrdinaryDiffEqRosenbrock, OrdinaryDiffEqTsit5, OrdinaryDiffEqVerner, Test
 
 # Sets stable rng number.
 using StableRNGs
@@ -176,7 +176,7 @@ let
     u0 = [:X1 => 1.0, :X2 => 3.0]
     ps = [:k1 => 2.0, :k2 => 3.0]
     oprob = ODEProblem(rn, u0, 1.0, ps)
-    osol = solve(oprob)
+    osol = solve(oprob, Tsit5())
     @test eltype(osol[:X1]) == eltype(osol[:X2]) == typeof(oprob[:X1]) == typeof(oprob[:X2]) == Float64
     @test eltype(osol.t) == typeof(oprob.tspan[1]) == typeof(oprob.tspan[2]) == Float64
 
@@ -184,7 +184,7 @@ let
     u0 = [:X1 => 1, :X2 => 3]
     ps = [:k1 => 2, :k2 => 3]
     oprob = ODEProblem(rn, u0, 1, ps)
-    osol = solve(oprob)
+    osol = solve(oprob, Tsit5())
     @test eltype(osol[:X1]) == eltype(osol[:X2]) == typeof(oprob[:X1]) == typeof(oprob[:X2]) == Float64
     @test eltype(osol.t) == Float64
 
@@ -192,7 +192,7 @@ let
     u0 = [:X1 => 1.0f0, :X2 => 3.0f0]
     ps = [:k1 => 2.0f0, :k2 => 3.0f0]
     oprob = ODEProblem(rn, u0, 1.0f0, ps)
-    osol = solve(oprob)
+    osol = solve(oprob, Tsit5())
     @test eltype(osol[:X1]) == eltype(osol[:X2]) == typeof(oprob[:X1]) == typeof(oprob[:X2]) == Float32
     @test eltype(osol.t) == typeof(oprob.tspan[1]) == typeof(oprob.tspan[2]) == Float32
 end
diff --git a/test/simulation_and_solving/solve_nonlinear.jl b/test/simulation_and_solving/solve_nonlinear.jl
index 9a6ef06533..11b4e111c8 100644
--- a/test/simulation_and_solving/solve_nonlinear.jl
+++ b/test/simulation_and_solving/solve_nonlinear.jl
@@ -78,7 +78,7 @@ end
 
 # Checks for system with conservation laws.
 # Checks using interfacing with output solution.
-@test_broken let  # 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
+let
     # Creates steady state network, unpack the parameter values.
     steady_state_network_3 = @reaction_network begin
         (p,d), 0 <--> X
@@ -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)
+    nl_prob_1 = NonlinearProblem(steady_state_network_3, u0, p; remove_conserved = true, conseqs_remake_warn = false)
     nl_prob_2 = NonlinearProblem(steady_state_network_3, u0, p)
 
     # Solves it using standard algorithm and simulation based algorithm.