From 4529be4d981f2e054c3cb7c2dd13079a223c7d77 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sun, 9 Jun 2024 15:45:14 -0400
Subject: [PATCH 01/33] init

---
 README.md                                     |   3 +-
 docs/make.jl                                  |   1 -
 docs/src/v14_migration_guide.md               | 188 ++++++++++++++++++
 .../petab_ode_param_fitting.md                |   0
 4 files changed, 190 insertions(+), 2 deletions(-)
 create mode 100644 docs/src/v14_migration_guide.md
 rename docs/{src/inverse_problems => unpublished}/petab_ode_param_fitting.md (100%)

diff --git a/README.md b/README.md
index 0f8d03e474..24f172a91c 100644
--- a/README.md
+++ b/README.md
@@ -34,7 +34,8 @@ etc).
 **NOTE:** version 14 is a breaking release, prompted by the release of ModelingToolkit.jl version 9. This caused several breaking changes in how Catalyst models are represented and interfaced with.
 
 Breaking changes and new functionality are summarized in the
-[HISTORY.md](HISTORY.md) file.
+[HISTORY.md](HISTORY.md) file. Furthermore, a migration guide on how to adapt your workflows to the 
+the new v14 update can be found [here](https://docs.sciml.ai/Catalyst/stable/v14_migration_guide/).
 
 ## Tutorials and documentation
 
diff --git a/docs/make.jl b/docs/make.jl
index 8595cd6ffa..af7e2ffd4a 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -40,7 +40,6 @@ makedocs(sitename = "Catalyst.jl",
     doctest = false,
     clean = true,
     pages = pages,
-    pagesonly = true,
     warnonly = [:missing_docs])
 
 deploydocs(repo = "github.com/SciML/Catalyst.jl.git";
diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
new file mode 100644
index 0000000000..20f1c18be6
--- /dev/null
+++ b/docs/src/v14_migration_guide.md
@@ -0,0 +1,188 @@
+# Version 14 Migration Guide
+
+Catalyst is built on the [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl) modelling language. A recent update of ModelingToolkit from version 8 to version 9 have required a corresponding update to Catalyst (from version 13 to 14). This update have introduced a couple of breaking changes, all of which will be detailed below.
+
+!!! note
+    Catalyst version 14 also introduces a large number of new features. These will not be discussed here, however, they are described in Catalyst's [history file]([@ref ](https://github.com/SciML/Catalyst.jl/blob/master/HISTORY.md)).
+
+## System completeness
+In ModelingToolkit v19 (and thus also Catalyst v14) all systems (e.g. `ReactionSystem`s and `ODESystem`s) are either *complete* or *incomplete* (completeness was a thing previously, however, recent updates mean that the user now have to be aware of this). Complete and incomplete systems differ in that
+- Only complete models can be used as inputs to simulations or certain tools for model analysis.
+- Only incomplete models can be [composed with other models to form hierarchical models](@ref compositional_modeling).
+
+A model's completeness depends on how it was created:
+- Models created programmatically (using the `ReactionSystem` constructor) are *incomplete*.
+- Models created using the `@reaction_network` DSL are *complete*.
+- To create *incomplete models using the DSL*, use the `@network_component` macro.
+- Models generated through the `compose` and `extend` functions are *incomplete*.
+
+Furthermore, any systems generated through e.g. `convert(ODESystem, rs)` are also complete. 
+
+Complete models can be generated from incomplete models through the `complete` function. Here is a workflow where we take completeness into account in the simulation of a simple birth-death process.
+```@example v14_migration_1
+using Catalyst
+t = default_t()
+@species X(t)
+@parameters p d
+rxs = [
+    Reaction(p, [], [X]),
+    Reaction(d, [X], [])
+]
+@named rs = ReactionSystem(rxs, t)
+```
+Here we have created an incomplete model. If we think our model is ready (i.e. we do not wish to compose it with additional models) we mark it as complete:
+```@example v14_migration_1
+rs = complete(rs)
+```
+Here, `complete` does not change the input model, but simply creates a new complete model. We hence overwrite our model variable (`rs`) with `complete`'s output. We can confirm that our model is complete using the `iscomplete` function
+```@example v14_migration_1
+Catalyst.iscomplete(rs)
+```
+We can now go on and use our model for e.g. simulations:
+```@example v14_migration_1
+using OrdinaryDiffEq, Plots
+u0 = [X => 0.1]
+tspan = (0.0, 10.0)
+ps = [d => 1.0, d => 0.2]
+oprob = ODEProblem(rs, u0, tspan, ps)
+sol = solve(oprob)
+plot(sol)
+```
+
+if we wish to first create a `ODESystem` from our `ReactionSystem`, the created `ODESystem` will be incomplete:
+```@example v14_migration_1
+osys = convert(ODESystem, rs)
+Catalyst.iscomplete(osys)
+```
+(note that `rs` must be complete before it can be converted to a `ODESystem`)
+
+If we now wish to create an `ODEProblem` from our `ODESystem`, we must first mark it as complete (using similar syntax as for our `ReactionSystem`):
+```@example v14_migration_1
+osys = complete(osys)
+oprob = ODEProblem(osys, u0, tspan, ps)
+sol = solve(oprob)
+plot(sol)
+```
+
+## Unknowns instead of states
+Previously, "states" was used as a term for a systems variables (both species and non-species variables). Now, the term "unknowns" is preferred instead. This means that there have been a number of changes to function names (e.g. `states` => `unknowns`, `get_states` => `get_unknowns`). 
+
+E.g. here we declare a `ReactionSystem` model containing both species and non-species unknowns:
+```@example v14_migration_2
+using Catalyst
+t = default_t()
+D = default_time_deriv()
+@species X(t)
+@variables V(t)
+@parameters p d Vmax
+
+eqs = [
+    Reaction(p, [], [X]),
+    Reaction(d, [X], []),
+    D(V) ~ Vmax - V*X*d/p
+]
+@named rs = ReactionSystem(eqs, t)
+```
+We can now use `unknowns` to retrieve all unknowns
+```@example v14_migration_2
+unknowns(rs)
+```
+Meanwhile, `species` and `nonspecies` (like previously) returns all species or non-species unknowns, respectively:
+```@example v14_migration_2
+species(rs)
+```
+```@example v14_migration_2
+nonspecies(rs)
+```
+
+## Lost support for most units
+As part of its v9 update, ModelingToolkit changed how units was handled. This includes using the package [DynamicQuantities.jl](https://github.com/SymbolicML/DynamicQuantities.jl) to manage units (instead of [Uniful.jl](https://github.com/PainterQubits/Unitful.jl), like previously).
+
+While this should lead to long-term improvements, unfortunately, as part of the process support for most units where removed. Currently, only the main SI units are supported (`s`, `m`, `kg`, `A`, `K`, `mol`, and `cd`). Composite units (e.g. `N = kg/(m^2)`) are no longer supported. Furthermore, prefix units (e.g. `mm = m/1000`) are not supported either. This means that most units relevant to Catalyst (such as `µM`) cannot be used directly. While composite units can still be written out in full and used (e.f. `mol/(m^3)`) this is hardly user friendly.
+
+The maintainers of ModelingTOolkit have been notified of this issue. We are unsure when this will be fixed, however, we do not think it will be a permanent change.
+
+## Removed support for system-mutating functions
+According to ModelingToolkit policy, created systems should not be modified. In accordance with this, the following functions have been deprecated: `addparam!`, `addreaction!`, `addspecies!`, `@add_reactions`, and `merge!`. Please use `ModelingToolkit.extend` and `ModelingToolkit.compose` to generate new merged and/or composed `ReactionSystems` from multiple component systems.
+
+It is still possible to add default values to a created `ReactionSystem`, i.e. the `setdefaults!` function is still supported.
+
+## New interface for creating time variable (`t`) and its differential (`D`)
+
+Previously, the time independent variable (typically called `t`) was declared using
+```@example v14_migration_3
+using Catalyst
+@variables t
+nothing # hide
+```
+A new, preferred, interface have now been introduced:
+```@example v14_migration_3
+t = default_t()
+nothing # hide
+```
+
+Similarly, the time differential (primarily relevant when creating combined reaction-equation models) used to be declared through
+```@example v14_migration_3
+D = Differential(t)
+nothing # hide
+```
+where the preferred method now being
+```@example v14_migration_3
+Dt = default_time_deriv()
+nothing # hide
+```
+
+## New interface for accessing problem/integrator/solution parameter (and species) value
+
+Previously, it was possible to index problems to query them for their parameter values. e.g.
+```@example v14_migration_4
+using Catalyst
+rn = @reaction_network begin
+    (p,d), 0 <--> X
+end
+oprob = ODEProblem(rn, [:X => 1.0], (0.0, 1.0), [:p => 1.0, :d => 0.2])
+nothing # hide
+```julia
+oprob[:p]
+```
+This *is no longer supported*. When you wish to query a problem (or integrator or solution) for a parameter value (or to update a parameter value), you must append `.ps` to the problem:
+```@example v14_migration_4
+oprob.ps[:p]
+```
+
+Furthermore, a few new functions (`getp`, `getu`, `setp`, `setu`) have been introduced. These can improve performance when querying a structure for a value multiple times (especially for very large models). These are describe in more detail [here](@ref simulation_structure_interfacing_functions).
+
+For more details on how to query various structures for parameter and species values, please read [this documentation page](@ref simulation_structure_interfacing).
+
+## Other notes
+Finally, here are some additional, minor, notes regarding the new update.
+
+### Modification of problems with conservation laws broken
+While it is possible to update e.g. `ODEProblem`s using the [`remake`](@ref simulation_structure_interfacing_problems_remake) function, this is not possible if the `remove_conserved = true` option was used. E.g. while
+```@example v14_migration_5
+using Catalyst, OrdinaryDiffEq
+rn = @reaction_network begin
+    (k1,k2), X1 <--> X2
+end
+u0 = [:X1 => 1.0, :X2 => 2.0]
+ps = [:k1 => 0.5, :k2 => 3.0]
+oprob = ODEProblem(rn, u0, (0.0, 10.0), ps; remove_conserved = true)
+sol(oprob)
+# hide
+```
+is perfectly fine, attempting to then modify `oprob` (in any manner!) is not possible (without introducing a bug):
+```@example v14_migration_5
+oprob_remade = remake(oprob; u0 = [:X1 => 5.0]) # NEVER do this.
+sol(oprob)
+# hide
+```
+
+This bug was likely present on earlier versions as well, but was only recently discovered. While we hope it will be fixed soon, the bug is in ModelingToolkit, and will not be fixed until its maintainers find the time to do so.
+
+### Depending on parameter order even more dangerous than before
+In early versions of Catalyst, parameters and species were provided as vectors (e.g. `[1.0, 2.0]`) rather than maps (e.g. `[p => 1.0, d => 2.0]`). While it has been *strongly* recommended to use the map form for a while, the vector form is still (technically) possible (however, we would not guarantee that this would produce expected behaviours). However, due to recent internal ModelingToolkit updates, the risk of unexpected behaviour when providing parameter values as vector is even larger than before.
+
+*Users should never use vector-forms to represent parameter and species values*
+
+### Additional deprecated functions
+The `reactionparams`, `numreactionparams`, and `reactionparamsmap` functions have been deprecated.
\ No newline at end of file
diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/unpublished/petab_ode_param_fitting.md
similarity index 100%
rename from docs/src/inverse_problems/petab_ode_param_fitting.md
rename to docs/unpublished/petab_ode_param_fitting.md

From 332fb318cccd1ea4cb4aff32e85ab243bb27daeb Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sun, 9 Jun 2024 16:40:21 -0400
Subject: [PATCH 02/33] update history file

---
 HISTORY.md                      | 132 +++++++++++++++++++-------------
 docs/make.jl                    |   1 +
 docs/src/v14_migration_guide.md |  67 ++++++++--------
 3 files changed, 114 insertions(+), 86 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index 270d64b177..d1558c122d 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -3,8 +3,82 @@
 ## Catalyst unreleased (master branch)
 
 ## Catalyst 14.0
-- The `reactionparams`, `numreactionparams`, and `reactionparamsmap`  functions have been removed.
+
+#### Breaking changes 
+Catalyst v14 was prompted by the release of ModelingToolkit v9. This introduced several breaking changes to the package. A summary of these (and how to handle them) can be found [here](https://docs.sciml.ai/Catalyst/stable/v14_migration_guide/). These are briefly summarised in the following bullet points:
+- `ReactionSystem`s must now be *complete* before they are exposed to most forms of simulation and analysis. With the exception of `ReactionSystem`s created through the `@reaction_network` macro, all `ReactionSystem`s are created incomplete. The `complete` function can be used to generate complete `ReactionSystem`s from incomplete ones.
+- The `states` function has been replaced with `unknowns`. The `get_states` function has been replaced with `get_unknowns`.
+- Support for most units (with the exception of `s`, `m`, `kg`, `A`, `K`, `mol`, and `cd`) has been dropped until further notice.
+- Problem parameter values are now accessed through `prob.ps[p]` (rather than `prob[p]`).
+- A significant bug prevents the safe application of the `remake` function on problems for which `remove_conserved = true` was used.
+- The `reactionparams`, `numreactionparams`, and `reactionparamsmap` functions have been removed.
 - To be more consistent with ModelingToolkit's immutability requirement for systems, we have removed API functions that mutate `ReactionSystem`s such as `addparam!`, `addreaction!`, `addspecies`, `@add_reactions`, and `merge!`. Please use `ModelingToolkit.extend` and `ModelingToolkit.compose` to generate new merged and/or composed `ReactionSystem`s from multiple component systems.
+
+#### General changes
+- The `default_t()` and `default_time_deriv()` functions are now the preferred approaches for creating the default time independent variable and its differential.
+- It is now possible to add metadata to individual reactions, e.g. using:
+```julia
+rn = @reaction_network begin
+    @parameters η
+    k, 2X --> X2, [description="Dimerisation"]
+end
+getdescription(rn)
+```
+a more detailed description can be found [here](https://docs.sciml.ai/Catalyst/dev/model_creation/dsl_advanced/#dsl_advanced_options_reaction_metadata).
+- `SDEProblem` no longer takes the `noise_scaling` argument. Noise scaling is now handled through the `noise_scaling` metadata (described in more detail [here](https://docs.sciml.ai/Catalyst/stable/model_simulation/simulation_introduction/#simulation_intro_SDEs_noise_saling))
+- Fields of the internal `Reaction` structure have been changed. `ReactionSystems`s saved using `serialize` on previous Catalyst versions cannot be loaded using this (or later) versions.
+- A new function, `save_reactionsystem`, which permits the writing of `ReactionSystem` models to files, has been created. A thorough description of this function can be found [here](https://docs.sciml.ai/Catalyst/stable/model_creation/model_file_loading_and_export/#Saving-Catalyst-models-to,-and-loading-them-from,-Julia-files)
+- Update how compounds are created. E.g. use
+```julia
+@variables t C(t) O(t)
+@compound CO2 ~ C + 2O
+```
+to create a compound species `CO2` that consists of `C` and 2 `O`.
+- Added documentation for chemistry-related functionality (compound creation and reaction balancing).
+- Added function `isautonomous` to check if a `ReactionSystem` is autonomous.
+- Added function `steady_state_stability` to compute stability for steady states. Example:
+```julia
+# Creates model.
+rn = @reaction_network begin
+    (p,d), 0 <--> X
+end
+p = [:p => 1.0, :d => 0.5]
+
+# Finds (the trivial) steady state, and computes stability.
+steady_state = [2.0]
+steady_state_stability(steady_state, rn, p)
+```
+Here, `steady_state_stability` takes an optional argument `tol = 10*sqrt(eps())`, which is used to determine whether a eigenvalue real part is reliably less than 0.
+- Added a DSL option, `@combinatoric_ratelaws`, which can be used to toggle whether to use combinatorial rate laws within the DSL. Example:
+```julia
+# Creates model.
+rn = @reaction_network begin
+    @combinatoric_ratelaws false
+    (kB,kD), 2X <--> X2
+end
+```
+- Added a DSL option, `@observables` for [creating observables](https://docs.sciml.ai/Catalyst/stable/model_creation/dsl_advanced/#dsl_advanced_options_observables) (this feature was already supported for programmatic modelling).
+- Added DSL options `@continuous_events` and `@discrete_events` to add events to a model as part of its creation (this feature was already supported for programmatic modelling). Example:
+```julia
+rn = @reaction_network begin
+    @continuous_events begin
+        [X ~ 1.0] => [X ~ X + 1.0]
+    end
+    d, X --> 0
+end
+```
+- Added DSL option `@equations` to add (algebraic or differential) equations to a model as part of its creation (this feature was already supported for programmatic modelling). Example:
+```julia
+rn = @reaction_network begin
+    @equations begin
+        D(V) ~ 1 - V
+    end
+    (p/V,d/V), 0 <--> X
+end
+```
+- Coupled reaction network + differential equation (or algebraic differential equation) systems can now be converted to `SDESystem`s and `NonlinearSystem`s.
+  
+#### Structural identifiability extension 
 - Added CatalystStructuralIdentifiabilityExtension, which permits StructuralIdentifiability.jl function to be applied directly to Catalyst systems. E.g. use
 ```julia
 using Catalyst, StructuralIdentifiability
@@ -17,47 +91,9 @@ assess_identifiability(goodwind_oscillator; measured_quantities=[:M])
 ```
 to assess (global) structural identifiability for all parameters and variables of the `goodwind_oscillator` model (under the presumption that we can measure `M` only).
 - Automatically handles conservation laws for structural identifiability problems (eliminates these internally to speed up computations).
-- Adds a tutorial to illustrate the use of the extension.
-- Enable adding metadata to individual reactions, e.g:
-```julia
-rn = @reaction_network begin
-    @parameters η
-    k, 2X --> X2, [noise_scaling=η]
-end
-getnoisescaling(rn)
-```
-- `SDEProblem` no longer takes the `noise_scaling` argument (see above for new approach to handle noise scaling).
-- Changed fields of internal `Reaction` structure. `ReactionSystems`s saved using `serialize` on previous Catalyst versions cannot be loaded using this (or later) versions.
-- Simulation of spatial ODEs now supported. For full details, please see https://github.com/SciML/Catalyst.jl/pull/644 and upcoming documentation. Note that these methods are currently considered alpha, with the interface and approach changing even in non-breaking Catalyst releases.
-- LatticeReactionSystem structure represents a spatial reaction network:
-  ```julia
-  rn = @reaction_network begin
-      (p,d), 0 <--> X
-  end
-  tr = @transport_reaction D X
-  lattice = Graphs.grid([5, 5])
-  lrs = LatticeReactionSystem(rn, [tr], lattice)
-```
-- Here, if a `u0` or `p` vector is given with scalar values:
-  ```julia
-  u0 = [:X => 1.0]
-  p = [:p => 1.0, :d => 0.5, :D => 0.1]
-  ```
-  this value will be used across the entire system. If their values are instead vectors, different values are used across the spatial system. Here
-  ```julia
-  X0 = zeros(25)
-  X0[1] = 1.0
-  u0 = [:X => X0]
-  ```
-  X's value will be `1.0` in the first vertex, but `0.0` in the remaining one (the system have 25 vertexes in total). SInce th parameters `p` and `d` are part of the non-spatial reaction network, their values are tied to vertexes. However, if the `D` parameter (which governs diffusion between vertexes) is given several values, these will instead correspond to the specific edges (and transportation along those edges.)
+- A more detailed of how this extension works can be found [here](https://docs.sciml.ai/Catalyst/stable/inverse_problems/structural_identifiability/).
 
-- Update how compounds are created. E.g. use
-```julia
-@variables t C(t) O(t)
-@compound CO2 ~ C + 2O
-```
-to create a compound species `CO2` that consists of `C` and 2 `O`.
-- Added documentation for chemistry related functionality (compound creation and reaction balancing).
+#### Bifurcation identifiability extension 
 - Add a CatalystBifurcationKitExtension, permitting BifurcationKit's `BifurcationProblem`s to be created from Catalyst reaction networks. Example usage:
 ```julia
 using Catalyst
@@ -86,20 +122,6 @@ plot(bif_dia; xguide="k1", yguide="X")
 ```
 - Automatically handles elimination of conservation laws for computing bifurcation diagrams.
 - Updated Bifurcation documentation with respect to this new feature.
-- Added function `isautonomous` to check if a `ReactionSystem` is autonomous.
-- Added function `steady_state_stability` to compute stability for steady states. Example:
-```julia
-# Creates model.
-rn = @reaction_network begin
-    (p,d), 0 <--> X
-end
-p = [:p => 1.0, :d => 0.5]
-
-# Finds (the trivial) steady state, and computes stability.
-steady_state = [2.0]
-steady_state_stability(steady_state, rn, p)
-```
-Here, `steady_state_stability` take an optional argument `tol = 10*sqrt(eps())`, which is used to determine whether a eigenvalue real part is reliably less that 0.
 
 ## Catalyst 13.5
 - Added a CatalystHomotopyContinuationExtension extension, which exports the `hc_steady_state` function if HomotopyContinuation is exported. `hc_steady_state` finds the steady states of a reaction system using the homotopy continuation method. This feature is only available for julia versions 1.9+. Example:
diff --git a/docs/make.jl b/docs/make.jl
index af7e2ffd4a..0d354c641a 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -40,6 +40,7 @@ makedocs(sitename = "Catalyst.jl",
     doctest = false,
     clean = true,
     pages = pages,
+    pagesonly = false,
     warnonly = [:missing_docs])
 
 deploydocs(repo = "github.com/SciML/Catalyst.jl.git";
diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 20f1c18be6..88d1ff6d00 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -1,19 +1,19 @@
 # Version 14 Migration Guide
 
-Catalyst is built on the [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl) modelling language. A recent update of ModelingToolkit from version 8 to version 9 have required a corresponding update to Catalyst (from version 13 to 14). This update have introduced a couple of breaking changes, all of which will be detailed below.
+Catalyst is built on the [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl) modelling language. A recent update of ModelingToolkit from version 8 to version 9 has required a corresponding update to Catalyst (from version 13 to 14). This update has introduced a couple of breaking changes, all of which will be detailed below.
 
 !!! note
-    Catalyst version 14 also introduces a large number of new features. These will not be discussed here, however, they are described in Catalyst's [history file]([@ref ](https://github.com/SciML/Catalyst.jl/blob/master/HISTORY.md)).
+ Catalyst version 14 also introduces several new features. These will not be discussed here, however, they are described in Catalyst's [history file](https://github.com/SciML/Catalyst.jl/blob/master/HISTORY.md).
 
 ## System completeness
-In ModelingToolkit v19 (and thus also Catalyst v14) all systems (e.g. `ReactionSystem`s and `ODESystem`s) are either *complete* or *incomplete* (completeness was a thing previously, however, recent updates mean that the user now have to be aware of this). Complete and incomplete systems differ in that
-- Only complete models can be used as inputs to simulations or certain tools for model analysis.
-- Only incomplete models can be [composed with other models to form hierarchical models](@ref compositional_modeling).
+In ModelingToolkit v9 (and thus also Catalyst v14) all systems (e.g. `ReactionSystem`s and `ODESystem`s) are either *complete* or *incomplete* (completeness was already a thing, however, recent updates mean that the user now has to be aware of this). Complete and incomplete systems differ in that
+- Only complete systems can be used as inputs to simulations or most tools for model analysis.
+- Only incomplete systems can be [composed with other systems to form hierarchical models](@ref compositional_modeling).
 
 A model's completeness depends on how it was created:
 - Models created programmatically (using the `ReactionSystem` constructor) are *incomplete*.
 - Models created using the `@reaction_network` DSL are *complete*.
-- To create *incomplete models using the DSL*, use the `@network_component` macro.
+- To create *incomplete models using the DSL*, use the `@network_component` macro (in all other aspects identical to `@reaction_network`).
 - Models generated through the `compose` and `extend` functions are *incomplete*.
 
 Furthermore, any systems generated through e.g. `convert(ODESystem, rs)` are also complete. 
@@ -30,11 +30,11 @@ rxs = [
 ]
 @named rs = ReactionSystem(rxs, t)
 ```
-Here we have created an incomplete model. If we think our model is ready (i.e. we do not wish to compose it with additional models) we mark it as complete:
+Here we have created an incomplete model. If our model is ready (i.e. we do not wish to compose it with additional models) we mark it as complete:
 ```@example v14_migration_1
 rs = complete(rs)
 ```
-Here, `complete` does not change the input model, but simply creates a new complete model. We hence overwrite our model variable (`rs`) with `complete`'s output. We can confirm that our model is complete using the `iscomplete` function
+Here, `complete` does not change the input model, but simply creates a new (complete) model. We hence overwrite our model variable (`rs`) with `complete`'s output. We can confirm that our model is complete using the `Catalyst.iscomplete` function:
 ```@example v14_migration_1
 Catalyst.iscomplete(rs)
 ```
@@ -49,12 +49,12 @@ sol = solve(oprob)
 plot(sol)
 ```
 
-if we wish to first create a `ODESystem` from our `ReactionSystem`, the created `ODESystem` will be incomplete:
+If we wish to first convert out `ReactionSystem` to an `ODESystem`, the `ODESystem` will be incomplete:
 ```@example v14_migration_1
 osys = convert(ODESystem, rs)
 Catalyst.iscomplete(osys)
 ```
-(note that `rs` must be complete before it can be converted to a `ODESystem`)
+(note that `rs` must be complete before it can be converted to an `ODESystem` or any other system type)
 
 If we now wish to create an `ODEProblem` from our `ODESystem`, we must first mark it as complete (using similar syntax as for our `ReactionSystem`):
 ```@example v14_migration_1
@@ -65,7 +65,7 @@ plot(sol)
 ```
 
 ## Unknowns instead of states
-Previously, "states" was used as a term for a systems variables (both species and non-species variables). Now, the term "unknowns" is preferred instead. This means that there have been a number of changes to function names (e.g. `states` => `unknowns`, `get_states` => `get_unknowns`). 
+Previously, "states" was used as a term for system variables (both species and non-species variables). Now, the term "unknowns" will be used instead. This means that there have been a number of changes to function names (e.g. `states` => `unknowns` and `get_states` => `get_unknowns`). 
 
 E.g. here we declare a `ReactionSystem` model containing both species and non-species unknowns:
 ```@example v14_migration_2
@@ -96,9 +96,9 @@ nonspecies(rs)
 ```
 
 ## Lost support for most units
-As part of its v9 update, ModelingToolkit changed how units was handled. This includes using the package [DynamicQuantities.jl](https://github.com/SymbolicML/DynamicQuantities.jl) to manage units (instead of [Uniful.jl](https://github.com/PainterQubits/Unitful.jl), like previously).
+As part of its v9 update, ModelingToolkit changed how units were handled. This includes using the package [DynamicQuantities.jl](https://github.com/SymbolicML/DynamicQuantities.jl) to manage units (instead of [Uniful.jl](https://github.com/PainterQubits/Unitful.jl), like previously).
 
-While this should lead to long-term improvements, unfortunately, as part of the process support for most units where removed. Currently, only the main SI units are supported (`s`, `m`, `kg`, `A`, `K`, `mol`, and `cd`). Composite units (e.g. `N = kg/(m^2)`) are no longer supported. Furthermore, prefix units (e.g. `mm = m/1000`) are not supported either. This means that most units relevant to Catalyst (such as `µM`) cannot be used directly. While composite units can still be written out in full and used (e.f. `mol/(m^3)`) this is hardly user friendly.
+While this should lead to long-term improvements, unfortunately, as part of the process support for most units was removed. Currently, only the main SI units are supported (`s`, `m`, `kg`, `A`, `K`, `mol`, and `cd`). Composite units (e.g. `N = kg/(m^2)`) are no longer supported. Furthermore, prefix units (e.g. `mm = m/1000`) are not supported either. This means that most units relevant to Catalyst (such as `µM`) cannot be used directly. While composite units can still be written out in full and used (e.g. `kg/(m^2)`) this is hardly user-friendly.
 
 The maintainers of ModelingTOolkit have been notified of this issue. We are unsure when this will be fixed, however, we do not think it will be a permanent change.
 
@@ -108,14 +108,13 @@ According to ModelingToolkit policy, created systems should not be modified. In
 It is still possible to add default values to a created `ReactionSystem`, i.e. the `setdefaults!` function is still supported.
 
 ## New interface for creating time variable (`t`) and its differential (`D`)
-
-Previously, the time independent variable (typically called `t`) was declared using
+Previously, the time-independent variable (typically called `t`) was declared using
 ```@example v14_migration_3
 using Catalyst
 @variables t
 nothing # hide
 ```
-A new, preferred, interface have now been introduced:
+A new, preferred, interface has now been introduced:
 ```@example v14_migration_3
 t = default_t()
 nothing # hide
@@ -126,43 +125,48 @@ Similarly, the time differential (primarily relevant when creating combined reac
 D = Differential(t)
 nothing # hide
 ```
-where the preferred method now being
+where the preferred method is now
 ```@example v14_migration_3
-Dt = default_time_deriv()
+D = default_time_deriv()
 nothing # hide
 ```
 
-## New interface for accessing problem/integrator/solution parameter (and species) value
+!!! note
+    If you look at ModelingToolkit documentation, these defaults are instead retrieved using `using ModelingToolkit: t_nounits as t, D_nounits as D`. This will also work, however, in Catalyst we have opted to instead use `default_t` and `default_time_deriv` as our main approach.
 
-Previously, it was possible to index problems to query them for their parameter values. e.g.
+## New interface for accessing problem/integrator/solution parameter (and species) value
+Previously, it was possible to directly index problems to query them for their parameter values. e.g.
 ```@example v14_migration_4
 using Catalyst
 rn = @reaction_network begin
-    (p,d), 0 <--> X
+ (p,d), 0 <--> X
 end
-oprob = ODEProblem(rn, [:X => 1.0], (0.0, 1.0), [:p => 1.0, :d => 0.2])
+u0 = [:X => 1.0]
+ps = [:p => 1.0, :d => 0.2]
+oprob = ODEProblem(rn, u0, (0.0, 1.0), ps)
 nothing # hide
+```
 ```julia
 oprob[:p]
 ```
-This *is no longer supported*. When you wish to query a problem (or integrator or solution) for a parameter value (or to update a parameter value), you must append `.ps` to the problem:
+This is *no longer supported*. When you wish to query a problem (or integrator or solution) for a parameter value (or to update a parameter value), you must append `.ps` to the problem variable name:
 ```@example v14_migration_4
 oprob.ps[:p]
 ```
 
-Furthermore, a few new functions (`getp`, `getu`, `setp`, `setu`) have been introduced. These can improve performance when querying a structure for a value multiple times (especially for very large models). These are describe in more detail [here](@ref simulation_structure_interfacing_functions).
+Furthermore, a few new functions (`getp`, `getu`, `setp`, `setu`) have been introduced. These can improve performance when querying a structure for a value multiple times (especially for very large models). These are described in more detail [here](@ref simulation_structure_interfacing_functions).
 
 For more details on how to query various structures for parameter and species values, please read [this documentation page](@ref simulation_structure_interfacing).
 
 ## Other notes
 Finally, here are some additional, minor, notes regarding the new update.
 
-### Modification of problems with conservation laws broken
-While it is possible to update e.g. `ODEProblem`s using the [`remake`](@ref simulation_structure_interfacing_problems_remake) function, this is not possible if the `remove_conserved = true` option was used. E.g. while
+#### Modification of problems with conservation laws broken
+While it is possible to update e.g. `ODEProblem`s using the [`remake`](@ref simulation_structure_interfacing_problems_remake) function, this is currently not possible if the `remove_conserved = true` option was used. E.g. while
 ```@example v14_migration_5
 using Catalyst, OrdinaryDiffEq
 rn = @reaction_network begin
-    (k1,k2), X1 <--> X2
+ (k1,k2), X1 <--> X2
 end
 u0 = [:X1 => 1.0, :X2 => 2.0]
 ps = [:k1 => 0.5, :k2 => 3.0]
@@ -170,19 +174,20 @@ oprob = ODEProblem(rn, u0, (0.0, 10.0), ps; remove_conserved = true)
 sol(oprob)
 # hide
 ```
-is perfectly fine, attempting to then modify `oprob` (in any manner!) is not possible (without introducing a bug):
+is perfectly fine, attempting to then modify `oprob` (in any manner!) is not possible:
 ```@example v14_migration_5
 oprob_remade = remake(oprob; u0 = [:X1 => 5.0]) # NEVER do this.
 sol(oprob)
 # hide
 ```
+This might generate a silent error, where the remade problem is different from the intended one. 
 
 This bug was likely present on earlier versions as well, but was only recently discovered. While we hope it will be fixed soon, the bug is in ModelingToolkit, and will not be fixed until its maintainers find the time to do so.
 
-### Depending on parameter order even more dangerous than before
-In early versions of Catalyst, parameters and species were provided as vectors (e.g. `[1.0, 2.0]`) rather than maps (e.g. `[p => 1.0, d => 2.0]`). While it has been *strongly* recommended to use the map form for a while, the vector form is still (technically) possible (however, we would not guarantee that this would produce expected behaviours). However, due to recent internal ModelingToolkit updates, the risk of unexpected behaviour when providing parameter values as vector is even larger than before.
+#### Depending on parameter order is even more dangerous than before
+In early versions of Catalyst, parameters and species were provided as vectors (e.g. `[1.0, 2.0]`) rather than maps (e.g. `[p => 1.0, d => 2.0]`). While we have already *strongly* recommended users to use the map form (or they might produce unintended results), the vector form is still (technically). Due to recent internal ModelingToolkit updates, the risk of unexpected behaviour when providing parameter values as vectors is even larger than before.
 
 *Users should never use vector-forms to represent parameter and species values*
 
-### Additional deprecated functions
+#### Additional deprecated functions
 The `reactionparams`, `numreactionparams`, and `reactionparamsmap` functions have been deprecated.
\ No newline at end of file

From 36682681b3f5b4154bad2a8745ce8b8c1858f320 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sun, 9 Jun 2024 16:41:24 -0400
Subject: [PATCH 03/33] minor history file updates

---
 HISTORY.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/HISTORY.md b/HISTORY.md
index d1558c122d..d3d88ff999 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -49,7 +49,7 @@ steady_state = [2.0]
 steady_state_stability(steady_state, rn, p)
 ```
 Here, `steady_state_stability` takes an optional argument `tol = 10*sqrt(eps())`, which is used to determine whether a eigenvalue real part is reliably less than 0.
-- Added a DSL option, `@combinatoric_ratelaws`, which can be used to toggle whether to use combinatorial rate laws within the DSL. Example:
+- Added a DSL option, `@combinatoric_ratelaws`, which can be used to toggle whether to use combinatorial rate laws within the DSL (this feature was already supported for programmatic modelling). Example:
 ```julia
 # Creates model.
 rn = @reaction_network begin

From 884cc7f1e930cdf30224212b6096fdce3119f4ac Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sun, 9 Jun 2024 16:44:09 -0400
Subject: [PATCH 04/33] migration guide code fixes

---
 docs/src/v14_migration_guide.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 88d1ff6d00..5d426a5e60 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -3,7 +3,7 @@
 Catalyst is built on the [ModelingToolkit.jl](https://github.com/SciML/ModelingToolkit.jl) modelling language. A recent update of ModelingToolkit from version 8 to version 9 has required a corresponding update to Catalyst (from version 13 to 14). This update has introduced a couple of breaking changes, all of which will be detailed below.
 
 !!! note
- Catalyst version 14 also introduces several new features. These will not be discussed here, however, they are described in Catalyst's [history file](https://github.com/SciML/Catalyst.jl/blob/master/HISTORY.md).
+    Catalyst version 14 also introduces several new features. These will not be discussed here, however, they are described in Catalyst's [history file](https://github.com/SciML/Catalyst.jl/blob/master/HISTORY.md).
 
 ## System completeness
 In ModelingToolkit v9 (and thus also Catalyst v14) all systems (e.g. `ReactionSystem`s and `ODESystem`s) are either *complete* or *incomplete* (completeness was already a thing, however, recent updates mean that the user now has to be aware of this). Complete and incomplete systems differ in that
@@ -43,7 +43,7 @@ We can now go on and use our model for e.g. simulations:
 using OrdinaryDiffEq, Plots
 u0 = [X => 0.1]
 tspan = (0.0, 10.0)
-ps = [d => 1.0, d => 0.2]
+ps = [p => 1.0, d => 0.2]
 oprob = ODEProblem(rs, u0, tspan, ps)
 sol = solve(oprob)
 plot(sol)

From 34260b2bd2b73851722efd9b57b7908815b0a4bb Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Mon, 10 Jun 2024 05:35:14 +0200
Subject: [PATCH 05/33] Update HISTORY.md

---
 HISTORY.md | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index d3d88ff999..435dd22ff7 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -93,7 +93,7 @@ to assess (global) structural identifiability for all parameters and variables o
 - Automatically handles conservation laws for structural identifiability problems (eliminates these internally to speed up computations).
 - A more detailed of how this extension works can be found [here](https://docs.sciml.ai/Catalyst/stable/inverse_problems/structural_identifiability/).
 
-#### Bifurcation identifiability extension 
+#### Bifurcation analysis extension 
 - Add a CatalystBifurcationKitExtension, permitting BifurcationKit's `BifurcationProblem`s to be created from Catalyst reaction networks. Example usage:
 ```julia
 using Catalyst
@@ -110,15 +110,15 @@ bif_par = :k1
 u_guess = [:X => 5.0, :Y => 2.0]
 p_start = [:k1 => 4.0, :k2 => 1.0, :k3 => 1.0, :k4 => 1.5, :k5 => 1.25]
 plot_var = :X
-bprob = BifurcationProblem(wilhelm_2009_model, u_guess, p_start, bif_par; plot_var=plot_var)
+bprob = BifurcationProblem(wilhelm_2009_model, u_guess, p_start, bif_par; plot_var = plot_var)
 
 p_span = (2.0, 20.0)
-opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps=1000)
+opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000)
 
-bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside=true)
+bif_dia = bifurcationdiagram(bprob, PALC(), 2, (args...) -> opts_br; bothside = true)
 
 using Plots
-plot(bif_dia; xguide="k1", yguide="X")
+plot(bif_dia; xguide = "k1", guide = "X")
 ```
 - Automatically handles elimination of conservation laws for computing bifurcation diagrams.
 - Updated Bifurcation documentation with respect to this new feature.

From a13d440f4db86bfde707dc6d746680b6e30c5fb8 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:51:13 -0400
Subject: [PATCH 06/33] 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 435dd22ff7..f7a0b2ef55 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -5,7 +5,7 @@
 ## Catalyst 14.0
 
 #### Breaking changes 
-Catalyst v14 was prompted by the release of ModelingToolkit v9. This introduced several breaking changes to the package. A summary of these (and how to handle them) can be found [here](https://docs.sciml.ai/Catalyst/stable/v14_migration_guide/). These are briefly summarised in the following bullet points:
+Catalyst v14 was prompted by the (breaking) release of ModelingToolkit v9, which introduced several breaking changes to Catalyst. A summary of these (and how to handle them) can be found [here](https://docs.sciml.ai/Catalyst/stable/v14_migration_guide/). These are briefly summarised in the following bullet points:
 - `ReactionSystem`s must now be *complete* before they are exposed to most forms of simulation and analysis. With the exception of `ReactionSystem`s created through the `@reaction_network` macro, all `ReactionSystem`s are created incomplete. The `complete` function can be used to generate complete `ReactionSystem`s from incomplete ones.
 - The `states` function has been replaced with `unknowns`. The `get_states` function has been replaced with `get_unknowns`.
 - Support for most units (with the exception of `s`, `m`, `kg`, `A`, `K`, `mol`, and `cd`) has been dropped until further notice.

From c3d848e55c598b6c0f3a04d3f7553e32b2d508e0 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:51:42 -0400
Subject: [PATCH 07/33] 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 f7a0b2ef55..f9ed15e7b1 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -6,7 +6,7 @@
 
 #### Breaking changes 
 Catalyst v14 was prompted by the (breaking) release of ModelingToolkit v9, which introduced several breaking changes to Catalyst. A summary of these (and how to handle them) can be found [here](https://docs.sciml.ai/Catalyst/stable/v14_migration_guide/). These are briefly summarised in the following bullet points:
-- `ReactionSystem`s must now be *complete* before they are exposed to most forms of simulation and analysis. With the exception of `ReactionSystem`s created through the `@reaction_network` macro, all `ReactionSystem`s are created incomplete. The `complete` function can be used to generate complete `ReactionSystem`s from incomplete ones.
+- `ReactionSystem`s must now be marked *complete* before they are exposed to most forms of simulation and analysis. With the exception of `ReactionSystem`s created through the `@reaction_network` macro, all `ReactionSystem`s are *not* marked complete upon construction. The `complete` function can be used to mark  `ReactionSystem`s as complete. To construct a `ReactionSystem` that is not marked complete via the DSL the new `@network_component` macro can be used.
 - The `states` function has been replaced with `unknowns`. The `get_states` function has been replaced with `get_unknowns`.
 - Support for most units (with the exception of `s`, `m`, `kg`, `A`, `K`, `mol`, and `cd`) has been dropped until further notice.
 - Problem parameter values are now accessed through `prob.ps[p]` (rather than `prob[p]`).

From 11030967813bda4bfdd15dbf9035533bb3029255 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:52:03 -0400
Subject: [PATCH 08/33] 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 f9ed15e7b1..e534bd6ee2 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -8,7 +8,7 @@
 Catalyst v14 was prompted by the (breaking) release of ModelingToolkit v9, which introduced several breaking changes to Catalyst. A summary of these (and how to handle them) can be found [here](https://docs.sciml.ai/Catalyst/stable/v14_migration_guide/). These are briefly summarised in the following bullet points:
 - `ReactionSystem`s must now be marked *complete* before they are exposed to most forms of simulation and analysis. With the exception of `ReactionSystem`s created through the `@reaction_network` macro, all `ReactionSystem`s are *not* marked complete upon construction. The `complete` function can be used to mark  `ReactionSystem`s as complete. To construct a `ReactionSystem` that is not marked complete via the DSL the new `@network_component` macro can be used.
 - The `states` function has been replaced with `unknowns`. The `get_states` function has been replaced with `get_unknowns`.
-- Support for most units (with the exception of `s`, `m`, `kg`, `A`, `K`, `mol`, and `cd`) has been dropped until further notice.
+- Support for most units (with the exception of `s`, `m`, `kg`, `A`, `K`, `mol`, and `cd`) has currently been dropped by ModelingToolkit, and hence they are unavailable via Catalyst too. Its is expected that eventually support for relevant chemical units such as molar will return to ModelingToolkit (and should then immediately work in Catalyst too).
 - Problem parameter values are now accessed through `prob.ps[p]` (rather than `prob[p]`).
 - A significant bug prevents the safe application of the `remake` function on problems for which `remove_conserved = true` was used.
 - The `reactionparams`, `numreactionparams`, and `reactionparamsmap` functions have been removed.

From 66cd51adf678934301477ab628532390200bdb47 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:52:20 -0400
Subject: [PATCH 09/33] 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 e534bd6ee2..ca8d97a9cf 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -10,7 +10,7 @@ Catalyst v14 was prompted by the (breaking) release of ModelingToolkit v9, which
 - The `states` function has been replaced with `unknowns`. The `get_states` function has been replaced with `get_unknowns`.
 - Support for most units (with the exception of `s`, `m`, `kg`, `A`, `K`, `mol`, and `cd`) has currently been dropped by ModelingToolkit, and hence they are unavailable via Catalyst too. Its is expected that eventually support for relevant chemical units such as molar will return to ModelingToolkit (and should then immediately work in Catalyst too).
 - Problem parameter values are now accessed through `prob.ps[p]` (rather than `prob[p]`).
-- A significant bug prevents the safe application of the `remake` function on problems for which `remove_conserved = true` was used.
+- A significant bug prevents the safe application of the `remake` function on problems for which `remove_conserved = true` was used when updating the values of initial conditions. Instead, the values of each conserved constant must be directly specified.
 - The `reactionparams`, `numreactionparams`, and `reactionparamsmap` functions have been removed.
 - To be more consistent with ModelingToolkit's immutability requirement for systems, we have removed API functions that mutate `ReactionSystem`s such as `addparam!`, `addreaction!`, `addspecies`, `@add_reactions`, and `merge!`. Please use `ModelingToolkit.extend` and `ModelingToolkit.compose` to generate new merged and/or composed `ReactionSystem`s from multiple component systems.
 

From d1897c59195f42462b986911f98cadcb43e074a0 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:52:30 -0400
Subject: [PATCH 10/33] 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 ca8d97a9cf..3713a9d47e 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -11,7 +11,7 @@ Catalyst v14 was prompted by the (breaking) release of ModelingToolkit v9, which
 - Support for most units (with the exception of `s`, `m`, `kg`, `A`, `K`, `mol`, and `cd`) has currently been dropped by ModelingToolkit, and hence they are unavailable via Catalyst too. Its is expected that eventually support for relevant chemical units such as molar will return to ModelingToolkit (and should then immediately work in Catalyst too).
 - Problem parameter values are now accessed through `prob.ps[p]` (rather than `prob[p]`).
 - A significant bug prevents the safe application of the `remake` function on problems for which `remove_conserved = true` was used when updating the values of initial conditions. Instead, the values of each conserved constant must be directly specified.
-- The `reactionparams`, `numreactionparams`, and `reactionparamsmap` functions have been removed.
+- The `reactionparams`, `numreactionparams`, and `reactionparamsmap` functions have been deprecated and removed.
 - To be more consistent with ModelingToolkit's immutability requirement for systems, we have removed API functions that mutate `ReactionSystem`s such as `addparam!`, `addreaction!`, `addspecies`, `@add_reactions`, and `merge!`. Please use `ModelingToolkit.extend` and `ModelingToolkit.compose` to generate new merged and/or composed `ReactionSystem`s from multiple component systems.
 
 #### General changes

From 3d55c45b25f66ff2e7bf69d7e09c40ade29e29cc Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:52:47 -0400
Subject: [PATCH 11/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 5d426a5e60..9ab3297b89 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -120,7 +120,7 @@ t = default_t()
 nothing # hide
 ```
 
-Similarly, the time differential (primarily relevant when creating combined reaction-equation models) used to be declared through
+Similarly, the time differential (primarily relevant when creating combined reaction-ODE models) used to be declared through
 ```@example v14_migration_3
 D = Differential(t)
 nothing # hide

From ab2ad80b59235a5cfad4a5cbed630f986dc738d7 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:52:57 -0400
Subject: [PATCH 12/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 9ab3297b89..187491edfe 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -174,7 +174,7 @@ oprob = ODEProblem(rn, u0, (0.0, 10.0), ps; remove_conserved = true)
 sol(oprob)
 # hide
 ```
-is perfectly fine, attempting to then modify `oprob` (in any manner!) is not possible:
+is perfectly fine, attempting to then modify any initial conditions or the value of the conservation constant in `oprob` will silently fail:
 ```@example v14_migration_5
 oprob_remade = remake(oprob; u0 = [:X1 => 5.0]) # NEVER do this.
 sol(oprob)

From 0f4930e0c294aabde2841412a806dfb11e9d0cda Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:53:12 -0400
Subject: [PATCH 13/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 187491edfe..6061460dee 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -182,7 +182,7 @@ sol(oprob)
 ```
 This might generate a silent error, where the remade problem is different from the intended one. 
 
-This bug was likely present on earlier versions as well, but was only recently discovered. While we hope it will be fixed soon, the bug is in ModelingToolkit, and will not be fixed until its maintainers find the time to do so.
+This bug was likely present on earlier versions as well, but was only recently discovered. While we hope it will be fixed soon, the issue is in ModelingToolkit, and will not be fixed until its maintainers find the time to do so.
 
 #### Depending on parameter order is even more dangerous than before
 In early versions of Catalyst, parameters and species were provided as vectors (e.g. `[1.0, 2.0]`) rather than maps (e.g. `[p => 1.0, d => 2.0]`). While we have already *strongly* recommended users to use the map form (or they might produce unintended results), the vector form is still (technically). Due to recent internal ModelingToolkit updates, the risk of unexpected behaviour when providing parameter values as vectors is even larger than before.

From c486367a1355d4b3d9ae0223956dbb4f75f56e75 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:53:20 -0400
Subject: [PATCH 14/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 6061460dee..d4d3c2d64e 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -180,7 +180,7 @@ oprob_remade = remake(oprob; u0 = [:X1 => 5.0]) # NEVER do this.
 sol(oprob)
 # hide
 ```
-This might generate a silent error, where the remade problem is different from the intended one. 
+This might generate a silent error, where the remade problem is different from the intended one (the value of the conserved constant will not be updated correctly).
 
 This bug was likely present on earlier versions as well, but was only recently discovered. While we hope it will be fixed soon, the issue is in ModelingToolkit, and will not be fixed until its maintainers find the time to do so.
 

From 13f3fe0083dbd93d2b3ab81cb19e4639283f5746 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:53:49 -0400
Subject: [PATCH 15/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index d4d3c2d64e..a6aecfa198 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -185,7 +185,7 @@ This might generate a silent error, where the remade problem is different from t
 This bug was likely present on earlier versions as well, but was only recently discovered. While we hope it will be fixed soon, the issue is in ModelingToolkit, and will not be fixed until its maintainers find the time to do so.
 
 #### Depending on parameter order is even more dangerous than before
-In early versions of Catalyst, parameters and species were provided as vectors (e.g. `[1.0, 2.0]`) rather than maps (e.g. `[p => 1.0, d => 2.0]`). While we have already *strongly* recommended users to use the map form (or they might produce unintended results), the vector form is still (technically). Due to recent internal ModelingToolkit updates, the risk of unexpected behaviour when providing parameter values as vectors is even larger than before.
+In early versions of Catalyst, parameters and species were provided as vectors (e.g. `[1.0, 2.0]`) rather than maps (e.g. `[p => 1.0, d => 2.0]`). While we previously *strongly* recommended users to use the map form (or they might produce unintended results), the vector form was still supported (technically). Due to recent internal ModelingToolkit updates, the purely numeric form is no longer supported and should never be used -- it will potentially lead to incorrect values for parameters and/or initial conditions. Note that if `rn` is a complete `ReactionSystem` you can now specify such mappings via `[rn.p => 1.0, rn.d => 2.0]`.
 
 *Users should never use vector-forms to represent parameter and species values*
 

From 0c38ab8a6eb46bc6edf39d10d8075ae886de8a0f Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:54:35 -0400
Subject: [PATCH 16/33] Update HISTORY.md

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

diff --git a/HISTORY.md b/HISTORY.md
index 3713a9d47e..6318a804c2 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -15,7 +15,15 @@ Catalyst v14 was prompted by the (breaking) release of ModelingToolkit v9, which
 - To be more consistent with ModelingToolkit's immutability requirement for systems, we have removed API functions that mutate `ReactionSystem`s such as `addparam!`, `addreaction!`, `addspecies`, `@add_reactions`, and `merge!`. Please use `ModelingToolkit.extend` and `ModelingToolkit.compose` to generate new merged and/or composed `ReactionSystem`s from multiple component systems.
 
 #### General changes
-- The `default_t()` and `default_time_deriv()` functions are now the preferred approaches for creating the default time independent variable and its differential.
+- The `default_t()` and `default_time_deriv()` functions are now the preferred approaches for creating the default time independent variable and its differential. i.e. 
+  ```julia
+  # do
+  t = default_t()
+  @species A(t)
+
+  # avoid
+  @variables t
+  @species A(t)
 - It is now possible to add metadata to individual reactions, e.g. using:
 ```julia
 rn = @reaction_network begin

From 236c8c076fc108a19106375a14ebcdd23669a799 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:57:03 -0400
Subject: [PATCH 17/33] 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 6318a804c2..2ee70fc0b2 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -56,7 +56,7 @@ p = [:p => 1.0, :d => 0.5]
 steady_state = [2.0]
 steady_state_stability(steady_state, rn, p)
 ```
-Here, `steady_state_stability` takes an optional argument `tol = 10*sqrt(eps())`, which is used to determine whether a eigenvalue real part is reliably less than 0.
+Here, `steady_state_stability` takes an optional argument `tol = 10*sqrt(eps())`, which is used to check that the real part of all eigenvalues are at least `tol` away from zero. Eigenvalues within `tol` of zero indicate that stability may not be reliably calculated.
 - Added a DSL option, `@combinatoric_ratelaws`, which can be used to toggle whether to use combinatorial rate laws within the DSL (this feature was already supported for programmatic modelling). Example:
 ```julia
 # Creates model.

From aa031fca707616c9b5c9c73e6c5d0e6fbef326c2 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:57:22 -0400
Subject: [PATCH 18/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index a6aecfa198..dc07ffdf63 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -11,7 +11,7 @@ In ModelingToolkit v9 (and thus also Catalyst v14) all systems (e.g. `ReactionSy
 - Only incomplete systems can be [composed with other systems to form hierarchical models](@ref compositional_modeling).
 
 A model's completeness depends on how it was created:
-- Models created programmatically (using the `ReactionSystem` constructor) are *incomplete*.
+- Models created programmatically (using the `ReactionSystem` constructor) are *not marked as complete* by default.
 - Models created using the `@reaction_network` DSL are *complete*.
 - To create *incomplete models using the DSL*, use the `@network_component` macro (in all other aspects identical to `@reaction_network`).
 - Models generated through the `compose` and `extend` functions are *incomplete*.

From 06ff56bb8893d3c50f871acd1a20498d887dacf0 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:57:42 -0400
Subject: [PATCH 19/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index dc07ffdf63..0d31eb13a5 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -6,7 +6,7 @@ Catalyst is built on the [ModelingToolkit.jl](https://github.com/SciML/ModelingT
     Catalyst version 14 also introduces several new features. These will not be discussed here, however, they are described in Catalyst's [history file](https://github.com/SciML/Catalyst.jl/blob/master/HISTORY.md).
 
 ## System completeness
-In ModelingToolkit v9 (and thus also Catalyst v14) all systems (e.g. `ReactionSystem`s and `ODESystem`s) are either *complete* or *incomplete* (completeness was already a thing, however, recent updates mean that the user now has to be aware of this). Complete and incomplete systems differ in that
+In ModelingToolkit v9 (and thus also Catalyst v14) all systems (e.g. `ReactionSystem`s and `ODESystem`s) are either *complete* or *incomplete*. Complete and incomplete systems differ in that
 - Only complete systems can be used as inputs to simulations or most tools for model analysis.
 - Only incomplete systems can be [composed with other systems to form hierarchical models](@ref compositional_modeling).
 

From af04a85780ffeee06e8a0e3a9ef138f6feb30eff Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:57:52 -0400
Subject: [PATCH 20/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 0d31eb13a5..e5a9e1e06b 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -12,7 +12,7 @@ In ModelingToolkit v9 (and thus also Catalyst v14) all systems (e.g. `ReactionSy
 
 A model's completeness depends on how it was created:
 - Models created programmatically (using the `ReactionSystem` constructor) are *not marked as complete* by default.
-- Models created using the `@reaction_network` DSL are *complete*.
+- Models created using the `@reaction_network` DSL are *automatically marked as complete*.
 - To create *incomplete models using the DSL*, use the `@network_component` macro (in all other aspects identical to `@reaction_network`).
 - Models generated through the `compose` and `extend` functions are *incomplete*.
 

From 589812f1a7422bebe5f3e5df9e9d2509365411b4 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:58:22 -0400
Subject: [PATCH 21/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index e5a9e1e06b..b4c78c1ec5 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -13,7 +13,7 @@ In ModelingToolkit v9 (and thus also Catalyst v14) all systems (e.g. `ReactionSy
 A model's completeness depends on how it was created:
 - Models created programmatically (using the `ReactionSystem` constructor) are *not marked as complete* by default.
 - Models created using the `@reaction_network` DSL are *automatically marked as complete*.
-- To create *incomplete models using the DSL*, use the `@network_component` macro (in all other aspects identical to `@reaction_network`).
+- To *use the DSL to create models that are not marked as complete*, use the `@network_component` macro (which in all other aspects is identical to `@reaction_network`).
 - Models generated through the `compose` and `extend` functions are *incomplete*.
 
 Furthermore, any systems generated through e.g. `convert(ODESystem, rs)` are also complete. 

From 2d1a383072b3037017cf233a162f7b5083c25cf6 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 17:58:32 -0400
Subject: [PATCH 22/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index b4c78c1ec5..4ca82952b6 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -14,7 +14,7 @@ A model's completeness depends on how it was created:
 - Models created programmatically (using the `ReactionSystem` constructor) are *not marked as complete* by default.
 - Models created using the `@reaction_network` DSL are *automatically marked as complete*.
 - To *use the DSL to create models that are not marked as complete*, use the `@network_component` macro (which in all other aspects is identical to `@reaction_network`).
-- Models generated through the `compose` and `extend` functions are *incomplete*.
+- Models generated through the `compose` and `extend` functions are *not marked as complete*.
 
 Furthermore, any systems generated through e.g. `convert(ODESystem, rs)` are also complete. 
 

From a30f105fe8d35c5db0c40a447e94ad9af5c58cf9 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 18:00:37 -0400
Subject: [PATCH 23/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 4ca82952b6..878e9bbb92 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -30,7 +30,7 @@ rxs = [
 ]
 @named rs = ReactionSystem(rxs, t)
 ```
-Here we have created an incomplete model. If our model is ready (i.e. we do not wish to compose it with additional models) we mark it as complete:
+Here we have created a model that is not marked as complete. If our model is ready (i.e. we do not wish to compose it with additional models) we mark it as complete:
 ```@example v14_migration_1
 rs = complete(rs)
 ```

From 52fca236d206f4904209eb9ae239bbff0e0e4a5a Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 18:00:56 -0400
Subject: [PATCH 24/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 878e9bbb92..02fcdb9b21 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -34,7 +34,7 @@ Here we have created a model that is not marked as complete. If our model is rea
 ```@example v14_migration_1
 rs = complete(rs)
 ```
-Here, `complete` does not change the input model, but simply creates a new (complete) model. We hence overwrite our model variable (`rs`) with `complete`'s output. We can confirm that our model is complete using the `Catalyst.iscomplete` function:
+Here, `complete` does not change the input model, but simply creates a new model that is tagged as complete. We hence overwrite our model variable (`rs`) with `complete`'s output. We can confirm that our model is complete using the `Catalyst.iscomplete` function:
 ```@example v14_migration_1
 Catalyst.iscomplete(rs)
 ```

From 4ee1126df5b9a286fd5ce6c72c19469de063b740 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 18:01:12 -0400
Subject: [PATCH 25/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 02fcdb9b21..ad200439e6 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -114,7 +114,7 @@ using Catalyst
 @variables t
 nothing # hide
 ```
-A new, preferred, interface has now been introduced:
+MTKv9 has introduced a standard global time variable, and as such a new, preferred, interface has been developed:
 ```@example v14_migration_3
 t = default_t()
 nothing # hide

From 2e6511ac26ed1725172e95cdbf46a7052512bf29 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 18:01:29 -0400
Subject: [PATCH 26/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index ad200439e6..c6cf8fd531 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -49,7 +49,7 @@ sol = solve(oprob)
 plot(sol)
 ```
 
-If we wish to first convert out `ReactionSystem` to an `ODESystem`, the `ODESystem` will be incomplete:
+If we wish to first manually convert our `ReactionSystem` to an `ODESystem`, the generated `ODESystem` will *not* be marked as complete
 ```@example v14_migration_1
 osys = convert(ODESystem, rs)
 Catalyst.iscomplete(osys)

From e83139e424a58a6aea057a0b17c9c30161dd7c16 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 18:01:46 -0400
Subject: [PATCH 27/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index c6cf8fd531..ce6adf9d6e 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -65,7 +65,7 @@ plot(sol)
 ```
 
 ## Unknowns instead of states
-Previously, "states" was used as a term for system variables (both species and non-species variables). Now, the term "unknowns" will be used instead. This means that there have been a number of changes to function names (e.g. `states` => `unknowns` and `get_states` => `get_unknowns`). 
+Previously, "states" was used as a term for system variables (both species and non-species variables). MTKv9 has switched to using the term "unknowns" instead. This means that there have been a number of changes to function names (e.g. `states` => `unknowns` and `get_states` => `get_unknowns`). 
 
 E.g. here we declare a `ReactionSystem` model containing both species and non-species unknowns:
 ```@example v14_migration_2

From 253bd3fc3bf190878dcea9a70761b3b575b5e9d7 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 18:02:03 -0400
Subject: [PATCH 28/33] Update docs/src/v14_migration_guide.md

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index ce6adf9d6e..8634945793 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -103,7 +103,7 @@ While this should lead to long-term improvements, unfortunately, as part of the
 The maintainers of ModelingTOolkit have been notified of this issue. We are unsure when this will be fixed, however, we do not think it will be a permanent change.
 
 ## Removed support for system-mutating functions
-According to ModelingToolkit policy, created systems should not be modified. In accordance with this, the following functions have been deprecated: `addparam!`, `addreaction!`, `addspecies!`, `@add_reactions`, and `merge!`. Please use `ModelingToolkit.extend` and `ModelingToolkit.compose` to generate new merged and/or composed `ReactionSystems` from multiple component systems.
+According to the ModelingToolkit system API, systems should not be mutable. In accordance with this, the following functions have been deprecated: `addparam!`, `addreaction!`, `addspecies!`, `@add_reactions`, and `merge!`. Please use `ModelingToolkit.extend` and `ModelingToolkit.compose` to generate new merged and/or composed `ReactionSystems` from multiple component systems.
 
 It is still possible to add default values to a created `ReactionSystem`, i.e. the `setdefaults!` function is still supported.
 

From e3e17cee910d151679e2a5f93dbc7a89341e8209 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Fri, 28 Jun 2024 18:04:12 -0400
Subject: [PATCH 29/33] change wrong statement on convert compelteness

---
 docs/src/v14_migration_guide.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 8634945793..8e1e9737c5 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -16,7 +16,7 @@ A model's completeness depends on how it was created:
 - To *use the DSL to create models that are not marked as complete*, use the `@network_component` macro (which in all other aspects is identical to `@reaction_network`).
 - Models generated through the `compose` and `extend` functions are *not marked as complete*.
 
-Furthermore, any systems generated through e.g. `convert(ODESystem, rs)` are also complete. 
+Furthermore, any systems generated through e.g. `convert(ODESystem, rs)`  are *not marked as complete*. 
 
 Complete models can be generated from incomplete models through the `complete` function. Here is a workflow where we take completeness into account in the simulation of a simple birth-death process.
 ```@example v14_migration_1

From ea51a88e7b613294291faeafbf8d48cc2f476ef8 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Mon, 1 Jul 2024 11:40:40 -0400
Subject: [PATCH 30/33] doc code fix

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

diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 8e1e9737c5..176ec89dd1 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -139,7 +139,7 @@ Previously, it was possible to directly index problems to query them for their p
 ```@example v14_migration_4
 using Catalyst
 rn = @reaction_network begin
- (p,d), 0 <--> X
+    (p,d), 0 <--> X
 end
 u0 = [:X => 1.0]
 ps = [:p => 1.0, :d => 0.2]
@@ -166,18 +166,18 @@ While it is possible to update e.g. `ODEProblem`s using the [`remake`](@ref simu
 ```@example v14_migration_5
 using Catalyst, OrdinaryDiffEq
 rn = @reaction_network begin
- (k1,k2), X1 <--> X2
+    (k1,k2), X1 <--> X2
 end
 u0 = [:X1 => 1.0, :X2 => 2.0]
 ps = [:k1 => 0.5, :k2 => 3.0]
 oprob = ODEProblem(rn, u0, (0.0, 10.0), ps; remove_conserved = true)
-sol(oprob)
+solve(oprob)
 # hide
 ```
 is perfectly fine, attempting to then modify any initial conditions or the value of the conservation constant in `oprob` will silently fail:
 ```@example v14_migration_5
 oprob_remade = remake(oprob; u0 = [:X1 => 5.0]) # NEVER do this.
-sol(oprob)
+solve(oprob_remade)
 # hide
 ```
 This might generate a silent error, where the remade problem is different from the intended one (the value of the conserved constant will not be updated correctly).

From 5dc4654b9601dc20ab434a61d5923acdf04e0424 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Mon, 1 Jul 2024 19:45:02 -0400
Subject: [PATCH 31/33] update tests

---
 test/upstream/mtk_problem_inputs.jl     |  8 ++++----
 test/upstream/mtk_structure_indexing.jl | 27 ++++++++++++++-----------
 2 files changed, 19 insertions(+), 16 deletions(-)

diff --git a/test/upstream/mtk_problem_inputs.jl b/test/upstream/mtk_problem_inputs.jl
index c1f86ca900..8712c6b504 100644
--- a/test/upstream/mtk_problem_inputs.jl
+++ b/test/upstream/mtk_problem_inputs.jl
@@ -218,10 +218,10 @@ let
     ]
 
     # Loops through all potential parameter sets, checking their inputs yield errors.
-    for ps in [ps_valid; ps_invalid], u0 in [u0_valid; u0s_invalid]
+    for ps in [[ps_valid]; ps_invalid], u0 in [[u0_valid]; u0s_invalid]
         # Handles problems with/without tspan separately. Special check ensuring that valid inputs passes.
         for XProblem in [ODEProblem, SDEProblem, DiscreteProblem]
-            if (ps == ps_valid) && (u0 == u0_valid)
+            if isequal(ps, ps_valid) && isequal(u0, u0_valid)
                 XProblem(rn, u0, (0.0, 1.0), ps); @test true;
             else
                 # Several of these cases do not throw errors (https://github.com/SciML/ModelingToolkit.jl/issues/2624).
@@ -231,7 +231,7 @@ let
             end
         end
         for XProblem in [NonlinearProblem, SteadyStateProblem]
-            if (ps == ps_valid) && (u0 == u0_valid)
+            if isequal(ps, ps_valid) && isequal(u0, u0_valid)
                 XProblem(rn, u0, ps); @test true;
             else
                 @test_broken false
@@ -240,4 +240,4 @@ let
             end
         end
     end
-end
\ No newline at end of file
+end
diff --git a/test/upstream/mtk_structure_indexing.jl b/test/upstream/mtk_structure_indexing.jl
index 78b70eb279..e45c8ba413 100644
--- a/test/upstream/mtk_structure_indexing.jl
+++ b/test/upstream/mtk_structure_indexing.jl
@@ -46,11 +46,11 @@ begin
     eproblems = [eoprob, esprob, edprob, ejprob, enprob, essprob]
 
     # Creates integrators.
-    oint = init(oprob, Tsit5(); save_everystep=false)
-    sint = init(sprob, ImplicitEM(); save_everystep=false)
+    oint = init(oprob, Tsit5(); save_everystep = false)
+    sint = init(sprob, ImplicitEM(); save_everystep = false)
     jint = init(jprob, SSAStepper())
-    nint = init(nprob, NewtonRaphson(); save_everystep=false)
-    @test_broken ssint = init(ssprob, DynamicSS(Tsit5()); save_everystep=false) # https://github.com/SciML/SciMLBase.jl/issues/660
+    nint = init(nprob, NewtonRaphson(); save_everystep = false)
+    @test_broken ssint = init(ssprob, DynamicSS(Tsit5()); save_everystep = false) # https://github.com/SciML/SciMLBase.jl/issues/660
     integrators = [oint, sint, jint, nint]
     
     # Creates solutions.
@@ -64,14 +64,17 @@ end
 
 # Tests problem indexing and updating.
 let 
+
     @test_broken false # A few cases fails for SteadyStateProblem: https://github.com/SciML/SciMLBase.jl/issues/660
     @test_broken false # Most cases broken for Ensemble problems: https://github.com/SciML/SciMLBase.jl/issues/661
-    for prob in deepcopy(problems[1:end-1])
+    @test_broken false # Cannot run deepcopy of SteadyStateProblem and Nonlinear Problems.
+    @test_broken false # A few cases fails for JumpProblems.
+    for prob in deepcopy([oprob, sprob, dprob])
         # Get u values (including observables).
         @test prob[X] == prob[model.X] == prob[:X] == 4
         @test prob[XY] == prob[model.XY] == prob[:XY] == 9
         @test prob[[XY,Y]] == prob[[model.XY,model.Y]] == prob[[:XY,:Y]] == [9, 5]
-        @test_broken prob[(XY,Y)] == prob[(model.XY,model.Y)] == prob[(:XY,:Y)] == (9, 5)
+        @test prob[(XY,Y)] == prob[(model.XY,model.Y)] == prob[(:XY,:Y)] == (9, 5)
         @test getu(prob, X)(prob) == getu(prob, model.X)(prob) == getu(prob, :X)(prob) == 4
         @test getu(prob, XY)(prob) == getu(prob, model.XY)(prob) == getu(prob, :XY)(prob) == 9 
         @test getu(prob, [XY,Y])(prob) == getu(prob, [model.XY,model.Y])(prob) == getu(prob, [:XY,:Y])(prob) == [9, 5]  
@@ -117,8 +120,10 @@ end
 
 # Test remake function.
 let 
+    @test_broken false # A few cases fails for JumpProblems.
+    @test_broken false # Cannot run deepcopy of SteadyStateProblem and Nonlinear Problems.
     @test_broken false # Currently cannot be run for Ensemble problems: https://github.com/SciML/SciMLBase.jl/issues/661 (as indexing cannot be used to check values).
-    for prob in deepcopy(problems)
+    for prob in deepcopy([oprob, sprob, dprob])
         # Remake for all u0s.
         rp = remake(prob; u0 = [X => 1, Y => 2])
         @test rp[[X, Y]] == [1, 2]
@@ -155,10 +160,8 @@ end
 
 # Test integrator indexing.
 let 
-    @test_broken false # NOTE: Multiple problems for `nint` (https://github.com/SciML/SciMLBase.jl/issues/662).
-    @test_broken false # NOTE: Multiple problems for `jint` (https://github.com/SciML/SciMLBase.jl/issues/654).
     @test_broken false # NOTE: Cannot even create a `ssint` (https://github.com/SciML/SciMLBase.jl/issues/660).
-    for int in deepcopy([oint, sint])
+    for int in deepcopy([oint, sint, jint, nint])
         # Get u values.
         @test int[X] == int[model.X] == int[:X] == 4
         @test int[XY] == int[model.XY] == int[:XY] == 9
@@ -262,7 +265,7 @@ let
             @test sol[X] == sol[model.X] == sol[:X]
             @test sol[XY] == sol[model.XY][1] == sol[:XY]
             @test sol[[XY,Y]] == sol[[model.XY,model.Y]] == sol[[:XY,:Y]]
-            @test_broken sol[(XY,Y)] == sol[(model.XY,model.Y)] == sol[(:XY,:Y)]
+            @test sol[(XY,Y)] == sol[(model.XY,model.Y)] == sol[(:XY,:Y)]
             @test getu(sol, X)(sol) == getu(sol, model.X)(sol)[1] == getu(sol, :X)(sol)
             @test getu(sol, XY)(sol) == getu(sol, model.XY)(sol)[1] == getu(sol, :XY)(sol)
             @test getu(sol, [XY,Y])(sol) == getu(sol, [model.XY,model.Y])(sol) == getu(sol, [:XY,:Y])(sol)
@@ -282,7 +285,7 @@ end
 # Tests plotting.
 let 
     @test_broken false # Currently broken for `ssol` (https://github.com/SciML/SciMLBase.jl/issues/580) 
-    for sol in deepcopy([osol, jsol])
+    for sol in deepcopy([osol, ssol, jsol])
         # Single variable.
         @test length(plot(sol; idxs = X).series_list) == 1
         @test length(plot(sol; idxs = XY).series_list) == 1

From 74277db0ea5e5ac4a3560670b335898bdc6c1e55 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Mon, 1 Jul 2024 20:40:38 -0400
Subject: [PATCH 32/33] test updates

---
 test/upstream/mtk_problem_inputs.jl | 713 +++++++++++++---------------
 1 file changed, 340 insertions(+), 373 deletions(-)

diff --git a/test/upstream/mtk_problem_inputs.jl b/test/upstream/mtk_problem_inputs.jl
index 788645b009..61a09f0b9d 100644
--- a/test/upstream/mtk_problem_inputs.jl
+++ b/test/upstream/mtk_problem_inputs.jl
@@ -1,9 +1,8 @@
-#! format: off
-
 ### Prepares Tests ###
 
 # Fetch packages
-using Catalyst, JumpProcesses, NonlinearSolve, OrdinaryDiffEq, SteadyStateDiffEq, StochasticDiffEq, Test
+using Catalyst, JumpProcesses, NonlinearSolve, OrdinaryDiffEq, Plots, SteadyStateDiffEq, StochasticDiffEq, Test
+import ModelingToolkit: getp, getu, setp, setu
 
 # Sets rnd number.
 using StableRNGs
@@ -13,411 +12,379 @@ seed = rand(rng, 1:100)
 
 ### Basic Tests ###
 
-# Prepares a models and initial conditions/parameters (of different forms) to be used as problem inputs.
-begin 
+# Prepares a model and its problems, integrators, and solutions.
+begin
+    # Creates the model and unpacks its context.
     model = @reaction_network begin
-        @species Z(t) = Z0
-        @parameters k2=0.5 Z0
+        @observables XY ~ X + Y
         (kp,kd), 0 <--> X
         (k1,k2), X <--> Y
-        (k1,k2), Y <--> Z
     end
-    @unpack X, Y, Z, kp, kd, k1, k2, Z0 = model
-
-    u0_alts = [
-        # Vectors not providing default values.
-        [X => 4, Y => 5],
-        [model.X => 4, model.Y => 5],
-        [:X => 4, :Y => 5],
-        # Vectors providing default values.
-        [X => 4, Y => 5, Z => 10],
-        [model.X => 4, model.Y => 5, model.Z => 10],
-        [:X => 4, :Y => 5, :Z => 10],
-        # Dicts not providing default values.
-        Dict([X => 4, Y => 5]),
-        Dict([model.X => 4, model.Y => 5]),
-        Dict([:X => 4, :Y => 5]),
-        # Dicts providing default values.
-        Dict([X => 4, Y => 5, Z => 10]),
-        Dict([model.X => 4, model.Y => 5, model.Z => 10]),
-        Dict([:X => 4, :Y => 5, :Z => 10]),
-        # Tuples not providing default values.
-        (X => 4, Y => 5),
-        (model.X => 4, model.Y => 5),
-        (:X => 4, :Y => 5),
-        # Tuples providing default values.
-        (X => 4, Y => 5, Z => 10),
-        (model.X => 4, model.Y => 5, model.Z => 10),
-        (:X => 4, :Y => 5, :Z => 10)
-    ]
+    @unpack XY, X, Y, kp, kd, k1, k2 = model
+
+    # Sets problem inputs (to be used for all problem creations).
+    u0_vals = [X => 4, Y => 5]
     tspan = (0.0, 10.0)
-    p_alts = [
-        # Vectors not providing default values.
-        [kp => 1.0, kd => 0.1, k1 => 0.25, Z0 => 10],
-        [model.kp => 1.0, model.kd => 0.1, model.k1 => 0.25, model.Z0 => 10],
-        [:kp => 1.0, :kd => 0.1, :k1 => 0.25, :Z0 => 10],
-        # Vectors providing default values.
-        [kp => 1.0, kd => 0.1, k1 => 0.25, k2 => 0.5, Z0 => 10],
-        [model.kp => 1.0, model.kd => 0.1, model.k1 => 0.25, model.k2 => 0.5, model.Z0 => 10],
-        [:kp => 1.0, :kd => 0.1, :k1 => 0.25, :k2 => 0.5, :Z0 => 10],
-        # Dicts not providing default values.
-        Dict([kp => 1.0, kd => 0.1, k1 => 0.25, Z0 => 10]),
-        Dict([model.kp => 1.0, model.kd => 0.1, model.k1 => 0.25, model.Z0 => 10]),
-        Dict([:kp => 1.0, :kd => 0.1, :k1 => 0.25, :Z0 => 10]),
-        # Dicts providing default values.
-        Dict([kp => 1.0, kd => 0.1, k1 => 0.25, k2 => 0.5, Z0 => 10]),
-        Dict([model.kp => 1.0, model.kd => 0.1, model.k1 => 0.25, model.k2 => 0.5, model.Z0 => 10]),
-        Dict([:kp => 1.0, :kd => 0.1, :k1 => 0.25, :k2 => 0.5, :Z0 => 10]),
-        # Tuples not providing default values.
-        (kp => 1.0, kd => 0.1, k1 => 0.25, Z0 => 10),
-        (model.kp => 1.0, model.kd => 0.1, model.k1 => 0.25, model.Z0 => 10),
-        (:kp => 1.0, :kd => 0.1, :k1 => 0.25, :Z0 => 10),
-        # Tuples providing default values.
-        (kp => 1.0, kd => 0.1, k1 => 0.25, k2 => 0.5, Z0 => 10),
-        (model.kp => 1.0, model.kd => 0.1, model.k1 => 0.25, model.k2 => 0.5, model.Z0 => 10),
-        (:kp => 1.0, :kd => 0.1, :k1 => 0.25, :k2 => 0.5, :Z0 => 10),
-    ]
+    p_vals = [kp => 1.0, kd => 0.1, k1 => 0.25, k2 => 0.5]
+
+    # Creates problems.
+    oprob = ODEProblem(model, u0_vals, tspan, p_vals)
+    sprob = SDEProblem(model,u0_vals, tspan, p_vals)
+    dprob = DiscreteProblem(model, u0_vals, tspan, p_vals)
+    jprob = JumpProblem(model, deepcopy(dprob), Direct(); rng)
+    nprob = NonlinearProblem(model, u0_vals, p_vals)
+    ssprob = SteadyStateProblem(model, u0_vals, p_vals)
+    problems = [oprob, sprob, dprob, jprob, nprob, ssprob]
+
+    # Creates an `EnsembleProblem` for each problem.
+    eoprob = EnsembleProblem(oprob)
+    esprob = EnsembleProblem(sprob)
+    edprob = EnsembleProblem(dprob)
+    ejprob = EnsembleProblem(jprob)
+    enprob = EnsembleProblem(nprob)
+    essprob = EnsembleProblem(ssprob)
+    eproblems = [eoprob, esprob, edprob, ejprob, enprob, essprob]
+
+    # Creates integrators.
+    oint = init(oprob, Tsit5(); save_everystep = false)
+    sint = init(sprob, ImplicitEM(); save_everystep = false)
+    jint = init(jprob, SSAStepper())
+    nint = init(nprob, NewtonRaphson(); save_everystep = false)
+    @test_broken ssint = init(ssprob, DynamicSS(Tsit5()); save_everystep = false) # https://github.com/SciML/SciMLBase.jl/issues/660
+    integrators = [oint, sint, jint, nint]
+    
+    # Creates solutions.
+    osol = solve(oprob, Tsit5())
+    ssol = solve(sprob, ImplicitEM(); seed)
+    jsol = solve(jprob, SSAStepper(); seed)
+    nsol = solve(nprob, NewtonRaphson())
+    sssol = solve(nprob, DynamicSS(Tsit5()))
+    sols = [osol, ssol, jsol, nsol, sssol]
 end
 
-# Perform ODE simulations (singular and ensemble).
+# Tests problem indexing and updating.
 let 
-    # Creates normal and ensemble problems.
-    base_oprob = ODEProblem(model, u0_alts[1], tspan, p_alts[1])
-    base_sol = solve(base_oprob, Tsit5(); saveat = 1.0)
-    base_eprob = EnsembleProblem(base_oprob)
-    base_esol = solve(base_eprob, Tsit5(); trajectories = 2, saveat = 1.0)
-
-    # Simulates problems for all input types, checking that identical solutions are found.
-    for u0 in u0_alts, p in p_alts
-        oprob = remake(base_oprob; u0, p)
-        @test base_sol == solve(oprob, Tsit5(); saveat = 1.0)
-        eprob = remake(base_eprob; u0, p)
-        @test base_esol == solve(eprob, Tsit5(); trajectories = 2, saveat = 1.0)
+    @test_broken false # A few cases fails for JumpProblem: https://github.com/SciML/ModelingToolkit.jl/issues/2838
+    @test_broken false # A few cases fails for SteadyStateProblem: https://github.com/SciML/SciMLBase.jl/issues/660
+    @test_broken false # Most cases broken for Ensemble problems: https://github.com/SciML/SciMLBase.jl/issues/661
+    for prob in deepcopy([oprob, sprob, dprob, nprob])
+        # Get u values (including observables).
+        @test prob[X] == prob[model.X] == prob[:X] == 4
+        @test prob[XY] == prob[model.XY] == prob[:XY] == 9
+        @test prob[[XY,Y]] == prob[[model.XY,model.Y]] == prob[[:XY,:Y]] == [9, 5]
+        @test prob[(XY,Y)] == prob[(model.XY,model.Y)] == prob[(:XY,:Y)] == (9, 5)
+        @test getu(prob, X)(prob) == getu(prob, model.X)(prob) == getu(prob, :X)(prob) == 4
+        @test getu(prob, XY)(prob) == getu(prob, model.XY)(prob) == getu(prob, :XY)(prob) == 9 
+        @test getu(prob, [XY,Y])(prob) == getu(prob, [model.XY,model.Y])(prob) == getu(prob, [:XY,:Y])(prob) == [9, 5]  
+        @test getu(prob, (XY,Y))(prob) == getu(prob, (model.XY,model.Y))(prob) == getu(prob, (:XY,:Y))(prob) == (9, 5)
+
+        # Set u values.
+        prob[X] = 20
+        @test prob[X] == 20
+        prob[model.X] = 30
+        @test prob[X] == 30
+        prob[:X] = 40
+        @test prob[X] == 40
+        setu(prob, X)(prob, 50)
+        @test prob[X] == 50
+        setu(prob, model.X)(prob, 60)
+        @test prob[X] == 60
+        setu(prob, :X)(prob, 70)
+        @test prob[X] == 70
+
+        # Get p values.
+        @test prob.ps[kp] == prob.ps[model.kp] == prob.ps[:kp] == 1.0    
+        @test prob.ps[[k1,k2]] == prob.ps[[model.k1,model.k2]] == prob.ps[[:k1,:k2]] == [0.25, 0.5]
+        @test prob.ps[(k1,k2)] == prob.ps[(model.k1,model.k2)] == prob.ps[(:k1,:k2)] == (0.25, 0.5)
+        @test getp(prob, kp)(prob) == getp(prob, model.kp)(prob) == getp(prob, :kp)(prob) == 1.0
+        @test getp(prob, [k1,k2])(prob) == getp(prob, [model.k1,model.k2])(prob) == getp(prob, [:k1,:k2])(prob) == [0.25, 0.5]
+        @test getp(prob, (k1,k2))(prob) == getp(prob, (model.k1,model.k2))(prob) == getp(prob, (:k1,:k2))(prob) == (0.25, 0.5)
+        
+        # Set p values.
+        prob.ps[kp] = 2.0
+        @test prob.ps[kp] == 2.0
+        prob.ps[model.kp] = 3.0
+        @test prob.ps[kp] == 3.0
+        prob.ps[:kp] = 4.0
+        @test prob.ps[kp] == 4.0
+        setp(prob, kp)(prob, 5.0)
+        @test prob.ps[kp] == 5.0
+        setp(prob, model.kp)(prob, 6.0)
+        @test prob.ps[kp] == 6.0
+        setp(prob, :kp)(prob, 7.0)
+        @test prob.ps[kp] == 7.0
     end
 end
 
-# Perform SDE simulations (singular and ensemble).
+# Test remake function.
 let 
-    # Creates normal and ensemble problems.
-    base_sprob = SDEProblem(model, u0_alts[1], tspan, p_alts[1])
-    base_sol = solve(base_sprob, ImplicitEM(); seed, saveat = 1.0)
-    base_eprob = EnsembleProblem(base_sprob)
-    base_esol = solve(base_eprob, ImplicitEM(); seed, trajectories = 2, saveat = 1.0)
-
-    # Simulates problems for all input types, checking that identical solutions are found.
-    for u0 in u0_alts, p in p_alts
-        sprob = remake(base_sprob; u0, p)
-        @test base_sol == solve(sprob, ImplicitEM(); seed, saveat = 1.0)
-        eprob = remake(base_eprob; u0, p)
-        @test base_esol == solve(eprob, ImplicitEM(); seed, trajectories = 2, saveat = 1.0)
+    @test_broken false # Cannot check result for JumpProblem: https://github.com/SciML/ModelingToolkit.jl/issues/2838
+    @test_broken false # Cannot deepcopy SteadyStateProblem :https://github.com/SciML/ModelingToolkit.jl/issues/2837
+    @test_broken false # Currently cannot be run for Ensemble problems: https://github.com/SciML/SciMLBase.jl/issues/661 (as indexing cannot be used to check values).
+    for prob in deepcopy([oprob, sprob, dprob, nprob])
+        # Remake for all u0s.
+        rp = remake(prob; u0 = [X => 1, Y => 2])
+        @test rp[[X, Y]] == [1, 2]
+        rp = remake(prob; u0 = [model.X => 3, model.Y => 4])
+        @test rp[[X, Y]] == [3, 4]
+        rp = remake(prob; u0 = [:X => 5, :Y => 6])
+        @test rp[[X, Y]] == [5, 6]
+
+        # Remake for a single u0.
+        rp = remake(prob; u0 = [Y => 7])
+        @test rp[[X, Y]] == [4, 7]
+        rp = remake(prob; u0 = [model.Y => 8])
+        @test rp[[X, Y]] == [4, 8]
+        rp = remake(prob; u0 = [:Y => 9])
+        @test rp[[X, Y]] == [4, 9]
+
+        # Remake for all ps.
+        rp = remake(prob; p = [kp => 1.0, kd => 2.0, k1 => 3.0, k2 => 4.0])
+        @test rp.ps[[kp, kd, k1, k2]] == [1.0, 2.0, 3.0, 4.0]
+        rp = remake(prob; p = [model.kp => 5.0, model.kd => 6.0, model.k1 => 7.0, model.k2 => 8.0])
+        @test rp.ps[[kp, kd, k1, k2]] == [5.0, 6.0, 7.0, 8.0]
+        rp = remake(prob; p = [:kp => 9.0, :kd => 10.0, :k1 => 11.0, :k2 => 12.0])
+        @test rp.ps[[kp, kd, k1, k2]] == [9.0, 10.0, 11.0, 12.0]
+
+        # Remake for a single p.
+        rp = remake(prob; p = [k2 => 13.0])
+        @test rp.ps[[kp, kd, k1, k2]] == [1.0, 0.1, 0.25, 13.0]
+        rp = remake(prob; p = [model.k2 => 14.0])
+        @test rp.ps[[kp, kd, k1, k2]] == [1.0, 0.1, 0.25, 14.0]
+        rp = remake(prob; p = [:k2 => 15.0])
+        @test rp.ps[[kp, kd, k1, k2]] == [1.0, 0.1, 0.25, 15.0]
     end
 end
 
-# Perform jump simulations (singular and ensemble).
+# Test integrator indexing.
 let 
-    # Creates normal and ensemble problems.
-    base_dprob = DiscreteProblem(model, u0_alts[1], tspan, p_alts[1])
-    base_jprob = JumpProblem(model, base_dprob, Direct(); rng)
-    base_sol = solve(base_jprob, SSAStepper(); seed, saveat = 1.0)
-    base_eprob = EnsembleProblem(base_jprob)
-    base_esol = solve(base_eprob, SSAStepper(); seed, trajectories = 2, saveat = 1.0)
-
-    # Simulates problems for all input types, checking that identical solutions are found.
-    for u0 in u0_alts, p in p_alts
-        jprob = remake(base_jprob; u0, p)
-        @test base_sol == solve(base_jprob, SSAStepper(); seed, saveat = 1.0)
-        eprob = remake(base_eprob; u0, p)
-        @test base_esol == solve(eprob, SSAStepper(); seed, trajectories = 2, saveat = 1.0)
+    @test_broken false # NOTE: Cannot even create a `ssint` (https://github.com/SciML/SciMLBase.jl/issues/660).
+    for int in deepcopy([oint, sint, jint, nint])
+        # Get u values.
+        @test int[X] == int[model.X] == int[:X] == 4
+        @test int[XY] == int[model.XY] == int[:XY] == 9
+        @test int[[XY,Y]] == int[[model.XY,model.Y]] == int[[:XY,:Y]] == [9, 5]
+        @test int[(XY,Y)] == int[(model.XY,model.Y)] == int[(:XY,:Y)] == (9, 5)
+        @test getu(int, X)(int) == getu(int, model.X)(int) == getu(int, :X)(int) == 4
+        @test getu(int, XY)(int) == getu(int, model.XY)(int) == getu(int, :XY)(int) == 9 
+        @test getu(int, [XY,Y])(int) == getu(int, [model.XY,model.Y])(int) == getu(int, [:XY,:Y])(int) == [9, 5]  
+        @test getu(int, (XY,Y))(int) == getu(int, (model.XY,model.Y))(int) == getu(int, (:XY,:Y))(int) == (9, 5)
+
+        # Set u values.
+        int[X] = 20
+        @test int[X] == 20
+        int[model.X] = 30
+        @test int[X] == 30
+        int[:X] = 40
+        @test int[X] == 40
+        setu(int, X)(int, 50)
+        @test int[X] == 50
+        setu(int, model.X)(int, 60)
+        @test int[X] == 60
+        setu(int, :X)(int, 70)
+        @test int[X] == 70
+
+        # Get p values.
+        @test int.ps[kp] == int.ps[model.kp] == int.ps[:kp] == 1.0    
+        @test int.ps[[k1,k2]] == int.ps[[model.k1,model.k2]] == int.ps[[:k1,:k2]] == [0.25, 0.5]
+        @test int.ps[(k1,k2)] == int.ps[(model.k1,model.k2)] == int.ps[(:k1,:k2)] == (0.25, 0.5)
+        @test getp(int, kp)(int) == getp(int, model.kp)(int) == getp(int, :kp)(int) == 1.0
+        @test getp(int, [k1,k2])(int) == getp(int, [model.k1,model.k2])(int) == getp(int, [:k1,:k2])(int) == [0.25, 0.5]
+        @test getp(int, (k1,k2))(int) == getp(int, (model.k1,model.k2))(int) == getp(int, (:k1,:k2))(int) == (0.25, 0.5)
+        
+        # Set p values.
+        int.ps[kp] = 2.0
+        @test int.ps[kp] == 2.0
+        int.ps[model.kp] = 3.0
+        @test int.ps[kp] == 3.0
+        int.ps[:kp] = 4.0
+        @test int.ps[kp] == 4.0
+        setp(int, kp)(int, 5.0)
+        @test int.ps[kp] == 5.0
+        setp(int, model.kp)(int, 6.0)
+        @test int.ps[kp] == 6.0
+        setp(int, :kp)(int, 7.0)
+        @test int.ps[kp] == 7.0
     end
 end
 
-# Solves a nonlinear problem (EnsembleProblems are not possible for these).
-let
-    base_nlprob = NonlinearProblem(model, u0_alts[1], p_alts[1])
-    base_sol = solve(base_nlprob, NewtonRaphson())
-    for u0 in u0_alts, p in p_alts
-        nlprob = remake(base_nlprob; u0, p)
-        @test base_sol == solve(nlprob, NewtonRaphson())
+# Test solve's save_idxs argument.
+# Currently, `save_idxs` is broken with symbolic stuff (https://github.com/SciML/ModelingToolkit.jl/issues/1761).
+let 
+    for (prob, solver) in zip(deepcopy([oprob, sprob, jprob]), [Tsit5(), ImplicitEM(), SSAStepper()])
+        # Save single variable
+        @test_broken solve(prob, solver; seed, save_idxs=X)[X][1] == 4
+        @test_broken solve(prob, solver; seed, save_idxs=model.X)[X][1] == 4
+        @test_broken solve(prob, solver; seed, save_idxs=:X)[X][1] == 4
+
+        # Save observable.
+        @test_broken solve(prob, solver; seed, save_idxs=XY)[XY][1] == 9
+        @test_broken solve(prob, solver; seed, save_idxs=model.XY)[XY][1] == 9
+        @test_broken solve(prob, solver; seed, save_idxs=:XY)[XY][1] == 9
+
+        # Save vector of stuff.
+        @test_broken solve(prob, solver; seed, save_idxs=[XY,Y])[[XY,Y]][1] == [9, 5]
+        @test_broken solve(prob, solver; seed, save_idxs=[model.XY,model.Y])[[model.XY,model.Y]][1] == [9, 5]
+        @test_broken solve(prob, solver; seed, save_idxs=[:XY,:Y])[[:XY,:Y]][1] == [9, 5]
     end
 end
 
-# Perform steady state simulations (singular and ensemble).
+# Tests solution indexing.
 let 
-    # Creates normal and ensemble problems.
-    base_ssprob = SteadyStateProblem(model, u0_alts[1], p_alts[1])
-    base_sol = solve(base_ssprob, DynamicSS(Tsit5()))
-    base_eprob = EnsembleProblem(base_ssprob)
-    base_esol = solve(base_eprob, DynamicSS(Tsit5()); trajectories = 2)
-
-    # Simulates problems for all input types, checking that identical solutions are found.
-    for u0 in u0_alts, p in p_alts
-        ssprob = remake(base_ssprob; u0, p)
-        @test base_sol == solve(ssprob, DynamicSS(Tsit5()))
-        eprob = remake(base_eprob; u0, p)
-        @test base_esol == solve(eprob, DynamicSS(Tsit5()); trajectories = 2)
+    for sol in deepcopy([osol, ssol, jsol])
+        # Get u values.
+        @test sol[X][1] == sol[model.X][1] == sol[:X][1] == 4
+        @test sol[XY][1] == sol[model.XY][1] == sol[:XY][1] == 9
+        @test sol[[XY,Y]][1] == sol[[model.XY,model.Y]][1] == sol[[:XY,:Y]][1] == [9, 5]
+        @test sol[(XY,Y)][1] == sol[(model.XY,model.Y)][1] == sol[(:XY,:Y)][1] == (9, 5)
+        @test getu(sol, X)(sol)[1] == getu(sol, model.X)(sol)[1] == getu(sol, :X)(sol)[1] == 4
+        @test getu(sol, XY)(sol)[1] == getu(sol, model.XY)(sol)[1] == getu(sol, :XY)(sol)[1] == 9 
+        @test getu(sol, [XY,Y])(sol)[1] == getu(sol, [model.XY,model.Y])(sol)[1] == getu(sol, [:XY,:Y])(sol)[1] == [9, 5]  
+        @test getu(sol, (XY,Y))(sol)[1] == getu(sol, (model.XY,model.Y))(sol)[1] == getu(sol, (:XY,:Y))(sol)[1] == (9, 5)       
+
+        # Get u values via idxs and functional call.
+        @test osol(0.0; idxs=X) == osol(0.0; idxs=model.X) == osol(0.0; idxs=:X) == 4
+        @test osol(0.0; idxs=XY) == osol(0.0; idxs=model.XY) == osol(0.0; idxs=:XY) == 9
+        @test osol(0.0; idxs = [XY,Y]) == osol(0.0; idxs = [model.XY,model.Y]) == osol(0.0; idxs = [:XY,:Y]) == [9, 5]
+        @test_broken osol(0.0; idxs = (XY,Y)) == osol(0.0; idxs = (model.XY,model.Y)) == osol(0.0; idxs = (:XY,:Y)) == (9, 5) # https://github.com/SciML/SciMLBase.jl/issues/711
+
+        # Get p values.
+        @test sol.ps[kp] == sol.ps[model.kp] == sol.ps[:kp] == 1.0    
+        @test sol.ps[[k1,k2]] == sol.ps[[model.k1,model.k2]] == sol.ps[[:k1,:k2]] == [0.25, 0.5]
+        @test sol.ps[(k1,k2)] == sol.ps[(model.k1,model.k2)] == sol.ps[(:k1,:k2)] == (0.25, 0.5)
+        @test getp(sol, kp)(sol) == getp(sol, model.kp)(sol) == getp(sol, :kp)(sol) == 1.0
+        @test getp(sol, [k1,k2])(sol) == getp(sol, [model.k1,model.k2])(sol) == getp(sol, [:k1,:k2])(sol) == [0.25, 0.5]
+        @test getp(sol, (k1,k2))(sol) == getp(sol, (model.k1,model.k2))(sol) == getp(sol, (:k1,:k2))(sol) == (0.25, 0.5)
     end
-end
 
-### Vector Species/Parameters Tests ###
-
-begin 
-    # Declares the model (with vector species/parameters, with/without default values, and observables).
-    t = default_t()
-    @species X(t)[1:2] Y(t)[1:2] = [10.0, 20.0] XY(t)[1:2]
-    @parameters p[1:2] d[1:2] = [0.2, 0.5]
-    rxs = [
-        Reaction(p[1], [], [X[1]]),
-        Reaction(p[2], [], [X[2]]),
-        Reaction(d[1], [X[1]], []),
-        Reaction(d[2], [X[2]], []),
-        Reaction(p[1], [], [Y[1]]),
-        Reaction(p[2], [], [Y[2]]),
-        Reaction(d[1], [Y[1]], []),
-        Reaction(d[2], [Y[2]], [])
-    ]
-    observed = [XY[1] ~ X[1] + Y[1], XY[2] ~ X[2] + Y[2]]
-    @named model_vec = ReactionSystem(rxs, t; observed)
-    model_vec = complete(model_vec)
-
-    # Declares various u0 versions (scalarised and vector forms).
-    u0_alts_vec = [
-        # Vectors not providing default values.
-        [X => [1.0, 2.0]],
-        [X[1] => 1.0, X[2] => 2.0],
-        [model_vec.X => [1.0, 2.0]],
-        [model_vec.X[1] => 1.0, model_vec.X[2] => 2.0],
-        [:X => [1.0, 2.0]],
-        # Vectors providing default values.
-        [X => [1.0, 2.0], Y => [10.0, 20.0]],
-        [X[1] => 1.0, X[2] => 2.0, Y[1] => 10.0, Y[2] => 20.0],
-        [model_vec.X => [1.0, 2.0], model_vec.Y => [10.0, 20.0]],
-        [model_vec.X[1] => 1.0, model_vec.X[2] => 2.0, model_vec.Y[1] => 10.0, model_vec.Y[2] => 20.0],
-        [:X => [1.0, 2.0], :Y => [10.0, 20.0]],
-        # Dicts not providing default values.
-        Dict([X => [1.0, 2.0]]),
-        Dict([X[1] => 1.0, X[2] => 2.0]),
-        Dict([model_vec.X => [1.0, 2.0]]),
-        Dict([model_vec.X[1] => 1.0, model_vec.X[2] => 2.0]),
-        Dict([:X => [1.0, 2.0]]),
-        # Dicts providing default values.
-        Dict([X => [1.0, 2.0], Y => [10.0, 20.0]]),
-        Dict([X[1] => 1.0, X[2] => 2.0, Y[1] => 10.0, Y[2] => 20.0]),
-        Dict([model_vec.X => [1.0, 2.0], model_vec.Y => [10.0, 20.0]]),
-        Dict([model_vec.X[1] => 1.0, model_vec.X[2] => 2.0, model_vec.Y[1] => 10.0, model_vec.Y[2] => 20.0]),
-        Dict([:X => [1.0, 2.0], :Y => [10.0, 20.0]]),
-        # Tuples not providing default values.
-        (X => [1.0, 2.0]),
-        (X[1] => 1.0, X[2] => 2.0),
-        (model_vec.X => [1.0, 2.0]),
-        (model_vec.X[1] => 1.0, model_vec.X[2] => 2.0),
-        (:X => [1.0, 2.0]),
-        # Tuples providing default values.
-        (X => [1.0, 2.0], Y => [10.0, 20.0]),
-        (X[1] => 1.0, X[2] => 2.0, Y[1] => 10.0, Y[2] => 20.0),
-        (model_vec.X => [1.0, 2.0], model_vec.Y => [10.0, 20.0]),
-        (model_vec.X[1] => 1.0, model_vec.X[2] => 2.0, model_vec.Y[1] => 10.0, model_vec.Y[2] => 20.0),
-        (:X => [1.0, 2.0], :Y => [10.0, 20.0]),
-    ]
-
-    # Declares various ps versions (vector forms only).
-    p_alts_vec = [
-        # Vectors not providing default values.
-        [p => [1.0, 2.0]],
-        [model_vec.p => [1.0, 2.0]],
-        [:p => [1.0, 2.0]],
-        # Vectors providing default values.
-        [p => [4.0, 5.0], d => [0.2, 0.5]],
-        [model_vec.p => [4.0, 5.0], model_vec.d => [0.2, 0.5]],
-        [:p => [4.0, 5.0], :d => [0.2, 0.5]],
-        # Dicts not providing default values.
-        Dict([p => [1.0, 2.0]]),
-        Dict([model_vec.p => [1.0, 2.0]]),
-        Dict([:p => [1.0, 2.0]]),
-        # Dicts providing default values.
-        Dict([p => [4.0, 5.0], d => [0.2, 0.5]]),
-        Dict([model_vec.p => [4.0, 5.0], model_vec.d => [0.2, 0.5]]),
-        Dict([:p => [4.0, 5.0], :d => [0.2, 0.5]]),
-        # Tuples not providing default values.
-        (p => [1.0, 2.0]),
-        (model_vec.p => [1.0, 2.0]),
-        (:p => [1.0, 2.0]),
-        # Tuples providing default values.
-        (p => [4.0, 5.0], d => [0.2, 0.5]),
-        (model_vec.p => [4.0, 5.0], model_vec.d => [0.2, 0.5]),
-        (:p => [4.0, 5.0], :d => [0.2, 0.5]),
-    ]
-
-    # Declares a timespan.
-    tspan = (0.0, 10.0)
-end
-
-# Perform ODE simulations (singular and ensemble).
-let 
-    # Creates normal and ensemble problems.
-    base_oprob = ODEProblem(model_vec, u0_alts_vec[1], tspan, p_alts_vec[1])
-    base_sol = solve(base_oprob, Tsit5(); saveat = 1.0)
-    base_eprob = EnsembleProblem(base_oprob)
-    base_esol = solve(base_eprob, Tsit5(); trajectories = 2, saveat = 1.0)
-
-    # Simulates problems for all input types, checking that identical solutions are found.
-    @test_broken false # Cannot remake problem (https://github.com/SciML/ModelingToolkit.jl/issues/2804).
-    # for u0 in u0_alts_vec, p in p_alts_vec
-    #     oprob = remake(base_oprob; u0, p)
-    #     @test base_sol == solve(oprob, Tsit5(); saveat = 1.0)
-    #     eprob = remake(base_eprob; u0, p)
-    #     @test base_esol == solve(eprob, Tsit5(); trajectories = 2, saveat = 1.0)
-    # end
+    # Handles nonlinear and steady state solutions differently.
+    let
+        for sol in deepcopy([nsol, sssol])
+            # Get u values.
+            @test sol[X] == sol[model.X] == sol[:X]
+            @test sol[XY] == sol[model.XY][1] == sol[:XY]
+            @test sol[[XY,Y]] == sol[[model.XY,model.Y]] == sol[[:XY,:Y]]
+            @test sol[(XY,Y)] == sol[(model.XY,model.Y)] == sol[(:XY,:Y)]
+            @test getu(sol, X)(sol) == getu(sol, model.X)(sol)[1] == getu(sol, :X)(sol)
+            @test getu(sol, XY)(sol) == getu(sol, model.XY)(sol)[1] == getu(sol, :XY)(sol)
+            @test getu(sol, [XY,Y])(sol) == getu(sol, [model.XY,model.Y])(sol) == getu(sol, [:XY,:Y])(sol)
+            @test_broken getu(sol, (XY,Y))(sol) == getu(sol, (model.XY,model.Y))(sol) == getu(sol, (:XY,:Y))(sol)[1] # https://github.com/SciML/SciMLBase.jl/issues/710
+
+            # Get p values.
+            @test sol.ps[kp] == sol.ps[model.kp] == sol.ps[:kp]
+            @test sol.ps[[k1,k2]] == sol.ps[[model.k1,model.k2]] == sol.ps[[:k1,:k2]]
+            @test sol.ps[(k1,k2)] == sol.ps[(model.k1,model.k2)] == sol.ps[(:k1,:k2)]
+            @test getp(sol, kp)(sol) == getp(sol, model.kp)(sol) == getp(sol, :kp)(sol)
+            @test getp(sol, [k1,k2])(sol) == getp(sol, [model.k1,model.k2])(sol) == getp(sol, [:k1,:k2])(sol)
+            @test getp(sol, (k1,k2))(sol) == getp(sol, (model.k1,model.k2))(sol) == getp(sol, (:k1,:k2))(sol)
+        end
+    end
 end
 
-# Perform SDE simulations (singular and ensemble).
+# Tests plotting.
 let 
-    # Creates normal and ensemble problems.
-    base_sprob = SDEProblem(model_vec, u0_alts_vec[1], tspan, p_alts_vec[1])
-    base_sol = solve(base_sprob, ImplicitEM(); seed, saveat = 1.0)
-    base_eprob = EnsembleProblem(base_sprob)
-    base_esol = solve(base_eprob, ImplicitEM(); seed, trajectories = 2, saveat = 1.0)
-
-    # Simulates problems for all input types, checking that identical solutions are found.
-    @test_broken false # Cannot remake problem (https://github.com/SciML/ModelingToolkit.jl/issues/2804).
-    # for u0 in u0_alts_vec, p in p_alts_vec
-    #     sprob = remake(base_sprob; u0, p)
-    #     @test base_sol == solve(sprob, ImplicitEM(); seed, saveat = 1.0)
-    #     eprob = remake(base_eprob; u0, p)
-    #     @test base_esol == solve(eprob, ImplicitEM(); seed, trajectories = 2, saveat = 1.0)
-    # end
+    for sol in deepcopy([osol, jsol, ssol])
+        # Single variable.
+        @test length(plot(sol; idxs = X).series_list) == 1
+        @test length(plot(sol; idxs = XY).series_list) == 1
+        @test length(plot(sol; idxs = model.X).series_list) == 1
+        @test length(plot(sol; idxs = model.XY).series_list) == 1
+        @test length(plot(sol; idxs = :X).series_list) == 1
+        @test length(plot(sol; idxs = :XY).series_list) == 1
+
+        # As vector.
+        @test length(plot(sol; idxs = [X,Y]).series_list) == 2
+        @test length(plot(sol; idxs = [XY,Y]).series_list) == 2
+        @test length(plot(sol; idxs = [model.X,model.Y]).series_list) == 2
+        @test length(plot(sol; idxs = [model.XY,model.Y]).series_list) == 2
+        @test length(plot(sol; idxs = [:X,:Y]).series_list) == 2
+        @test length(plot(sol; idxs = [:XY,:Y]).series_list) == 2
+
+        # As tuple.
+        @test length(plot(sol; idxs = (X, Y)).series_list) == 1
+        @test length(plot(sol; idxs = (XY, Y)).series_list) == 1
+        @test length(plot(sol; idxs = (model.X, model.Y)).series_list) == 1
+        @test length(plot(sol; idxs = (model.XY, model.Y)).series_list) == 1
+        @test length(plot(sol; idxs = (:X, :Y)).series_list) == 1
+        @test length(plot(sol; idxs = (:XY, :Y)).series_list) == 1
+    end     
 end
 
-# Perform jump simulations (singular and ensemble).
-let 
-    # Creates normal and ensemble problems.
-    base_dprob = DiscreteProblem(model_vec, u0_alts_vec[1], tspan, p_alts_vec[1])
-    base_jprob = JumpProblem(model_vec, base_dprob, Direct(); rng)
-    base_sol = solve(base_jprob, SSAStepper(); seed, saveat = 1.0)
-    base_eprob = EnsembleProblem(base_jprob)
-    base_esol = solve(base_eprob, SSAStepper(); seed, trajectories = 2, saveat = 1.0)
-
-    # Simulates problems for all input types, checking that identical solutions are found.
-    @test_broken false # Cannot remake problem (https://github.com/SciML/ModelingToolkit.jl/issues/2804).
-    # for u0 in u0_alts_vec, p in p_alts_vec
-    #     jprob = remake(base_jprob; u0, p)
-    #     @test base_sol == solve(base_jprob, SSAStepper(); seed, saveat = 1.0)
-    #     eprob = remake(base_eprob; u0, p)
-    #     @test base_esol == solve(eprob, SSAStepper(); seed, trajectories = 2, saveat = 1.0)
-    # end
-end
 
-# Solves a nonlinear problem (EnsembleProblems are not possible for these).
-let
-    base_nlprob = NonlinearProblem(model_vec, u0_alts_vec[1], p_alts_vec[1])
-    base_sol = solve(base_nlprob, NewtonRaphson())
-    @test_broken false # Cannot remake problem (https://github.com/SciML/ModelingToolkit.jl/issues/2804).
-    # for u0 in u0_alts_vec, p in p_alts_vec
-    #     nlprob = remake(base_nlprob; u0, p)
-    #     @test base_sol == solve(nlprob, NewtonRaphson())
-    # end
-end
+### Tests For Hierarchical System ###
 
-# Perform steady state simulations (singular and ensemble).
-let 
-    # Creates normal and ensemble problems.
-    base_ssprob = SteadyStateProblem(model_vec, u0_alts_vec[1], p_alts_vec[1])
-    base_sol = solve(base_ssprob, DynamicSS(Tsit5()))
-    base_eprob = EnsembleProblem(base_ssprob)
-    base_esol = solve(base_eprob, DynamicSS(Tsit5()); trajectories = 2)
-
-    # Simulates problems for all input types, checking that identical solutions are found.
-    @test_broken false # Cannot remake problem (https://github.com/SciML/ModelingToolkit.jl/issues/2804).
-    # for u0 in u0_alts_vec, p in p_alts_vec
-    #     ssprob = remake(base_ssprob; u0, p)
-    #     @test base_sol == solve(ssprob, DynamicSS(Tsit5()))
-    #     eprob = remake(base_eprob; u0, p)
-    #     @test base_esol == solve(eprob, DynamicSS(Tsit5()); trajectories = 2)
-    # end
-end
+# TODO
 
-### Checks Errors On Faulty Inputs ###
+### Mass Action Jump Rate Updating Correctness ###
 
-# Checks various erroneous problem inputs, ensuring that these throw errors.
+# Checks that the rates of mass action jumps are correctly updated after parameter values are changed.
 let
-    # Declares the model.
+    # Creates the model.
     rn = @reaction_network begin
-        (k1,k2), X1 <--> X2
-    end
-    @unpack k1, k2, X1, X2 = rn
-    t = default_t()
-    @species X3(t)
-    @parameters k3
-
-    # Creates systems (so these are not recreated in each problem call).
-    osys = convert(ODESystem, rn)
-    ssys = convert(SDESystem, rn)
-    nsys = convert(NonlinearSystem, rn)
-
-    # Declares valid initial conditions and parameter values
-    u0_valid = [X1 => 1, X2 => 2]
-    ps_valid = [k1 => 0.5, k2 => 0.1]
-
-    # Declares invalid initial conditions and parameters. This includes both cases where values are
-    # missing, or additional ones are given. Includes vector/Tuple/Dict forms.
-    u0s_invalid = [
-        # Missing a value.
-        [X1 => 1],
-        [rn.X1 => 1],
-        [:X1 => 1],
-        Dict([X1 => 1]),
-        Dict([rn.X1 => 1]),
-        Dict([:X1 => 1]),
-        (X1 => 1),
-        (rn.X1 => 1),
-        (:X1 => 1),
-        # Contain an additional value.
-        [X1 => 1, X2 => 2, X3 => 3],
-        [:X1 => 1, :X2 => 2, :X3 => 3],
-        Dict([X1 => 1, X2 => 2, X3 => 3]),
-        Dict([:X1 => 1, :X2 => 2, :X3 => 3]),
-        (X1 => 1, X2 => 2, X3 => 3),
-        (:X1 => 1, :X2 => 2, :X3 => 3)
-    ]
-    ps_invalid = [
-        # Missing a value.
-        [k1 => 1.0],
-        [rn.k1 => 1.0],
-        [:k1 => 1.0],
-        Dict([k1 => 1.0]),
-        Dict([rn.k1 => 1.0]),
-        Dict([:k1 => 1.0]),
-        (k1 => 1.0),
-        (rn.k1 => 1.0),
-        (:k1 => 1.0),
-        # Contain an additional value.
-        [k1 => 1.0, k2 => 2.0, k3 => 3.0],
-        [:k1 => 1.0, :k2 => 2.0, :k3 => 3.0],
-        Dict([k1 => 1.0, k2 => 2.0, k3 => 3.0]),
-        Dict([:k1 => 1.0, :k2 => 2.0, :k3 => 3.0]),
-        (k1 => 1.0, k2 => 2.0, k3 => 3.0),
-        (:k1 => 1.0, :k2 => 2.0, :k3 => 3.0)
-    ]
-
-    # Loops through all potential parameter sets, checking their inputs yield errors.
-    for ps in [[ps_valid]; ps_invalid], u0 in [[u0_valid]; u0s_invalid]
-        # Handles problems with/without tspan separately. Special check ensuring that valid inputs passes.
-        for XProblem in [ODEProblem, SDEProblem, DiscreteProblem]
-            if isequal(ps, ps_valid) && isequal(u0, u0_valid)
-                XProblem(rn, u0, (0.0, 1.0), ps); @test true;
-            else
-                @test_broken false
-                continue
-                @test_throws Exception XProblem(xsys, u0, (0.0, 1.0), ps)
-            end
-        end
-        for XProblem in [NonlinearProblem, SteadyStateProblem]
-            if isequal(ps, ps_valid) && isequal(u0, u0_valid)
-                XProblem(rn, u0, ps); @test true;
-            else
-                @test_broken false
-                continue
-                @test_throws Exception XProblem(xsys, u0, ps)
-            end
-        end
+        p1*p2, A + B --> C
     end
+    @unpack p1, p2 = rn
+
+    # Creates a JumpProblem and integrator. Checks that the initial mass action rate is correct.
+    u0 = [:A => 1, :B => 2, :C => 3]
+    ps = [:p1 => 3.0, :p2 => 2.0]
+    dprob = DiscreteProblem(rn, u0, (0.0, 1.0), ps)
+    jprob = JumpProblem(rn, dprob, Direct())
+    jint = init(jprob, SSAStepper())
+    @test jprob.massaction_jump.scaled_rates[1] == 6.0
+
+    # Checks that the mass action rate is correctly updated after normal indexing.
+    jprob.ps[p1] = 4.0
+    @test jprob.massaction_jump.scaled_rates[1] == 8.0
+    jprob.ps[rn.p1] = 5.0
+    @test jprob.massaction_jump.scaled_rates[1] == 10.0
+    jprob.ps[:p1] = 6.0
+    @test jprob.massaction_jump.scaled_rates[1] == 12.0
+    setp(jprob, p1)(jprob, 7.0)
+    @test jprob.massaction_jump.scaled_rates[1] == 14.0
+    setp(jprob, rn.p1)(jprob, 8.0)
+    @test jprob.massaction_jump.scaled_rates[1] == 16.0
+    setp(jprob, :p1)(jprob, 3.0)
+    @test jprob.massaction_jump.scaled_rates[1] == 6.0
+
+    # Check that the mass action rate is correctly updated when `remake` is used.
+    # Checks both when partial and full parameter vectors are provided to `remake`.
+    @test remake(jprob; p = [p1 => 4.0]).massaction_jump.scaled_rates[1] == 8.0
+    @test remake(jprob; p = [rn.p1 => 5.0]).massaction_jump.scaled_rates[1] == 10.0
+    @test remake(jprob; p = [:p1 => 6.0]).massaction_jump.scaled_rates[1] == 12.0
+    @test remake(jprob; p = [p1 => 4.0, p2 => 3.0]).massaction_jump.scaled_rates[1] == 12.0
+    @test remake(jprob; p = [rn.p1 => 5.0, rn.p2 => 4.0]).massaction_jump.scaled_rates[1] == 20.0
+    @test remake(jprob; p = [:p1 => 6.0, :p2 => 5.0]).massaction_jump.scaled_rates[1] == 30.0
+
+    # Checks that updating an integrators parameter values does not affect mass action rate until after
+    # `reset_aggregated_jumps!` have been applied as well (wt which point the correct rate is achieved).
+    jint.ps[p1] = 4.0
+    @test jint.cb.condition.ma_jumps.scaled_rates[1] == 30.0
+    reset_aggregated_jumps!(jint)
+    @test jint.cb.condition.ma_jumps.scaled_rates[1] == 8.0
+    
+    jint.ps[rn.p1] = 5.0
+    @test jint.cb.condition.ma_jumps.scaled_rates[1] == 8.0
+    reset_aggregated_jumps!(jint)
+    @test jint.cb.condition.ma_jumps.scaled_rates[1] == 10.0
+    
+    jint.ps[:p1] = 6.0
+    @test jint.cb.condition.ma_jumps.scaled_rates[1] == 10.0
+    reset_aggregated_jumps!(jint)
+    @test jint.cb.condition.ma_jumps.scaled_rates[1] == 12.0
+    
+    setp(jint, p1)(jint, 7.0)
+    @test jint.cb.condition.ma_jumps.scaled_rates[1] == 12.0
+    reset_aggregated_jumps!(jint)
+    @test jint.cb.condition.ma_jumps.scaled_rates[1] == 14.0
+    
+    setp(jint, rn.p1)(jint, 8.0)
+    @test jint.cb.condition.ma_jumps.scaled_rates[1] == 14.0
+    reset_aggregated_jumps!(jint)
+    @test jint.cb.condition.ma_jumps.scaled_rates[1] == 16.0
+    
+    setp(jint, :p1)(jint, 3.0)
+    @test jint.cb.condition.ma_jumps.scaled_rates[1] == 16.0
+    reset_aggregated_jumps!(jint)
+    @test jint.cb.condition.ma_jumps.scaled_rates[1] == 6.0
 end

From d6f1588504f77633e7bab735a975601ea69bb6f5 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Mon, 1 Jul 2024 23:33:12 -0400
Subject: [PATCH 33/33] fix

---
 test/upstream/mtk_structure_indexing.jl | 22 ++++++++++------------
 1 file changed, 10 insertions(+), 12 deletions(-)

diff --git a/test/upstream/mtk_structure_indexing.jl b/test/upstream/mtk_structure_indexing.jl
index ae65360128..61a09f0b9d 100644
--- a/test/upstream/mtk_structure_indexing.jl
+++ b/test/upstream/mtk_structure_indexing.jl
@@ -64,14 +64,15 @@ end
 
 # Tests problem indexing and updating.
 let 
+    @test_broken false # A few cases fails for JumpProblem: https://github.com/SciML/ModelingToolkit.jl/issues/2838
     @test_broken false # A few cases fails for SteadyStateProblem: https://github.com/SciML/SciMLBase.jl/issues/660
     @test_broken false # Most cases broken for Ensemble problems: https://github.com/SciML/SciMLBase.jl/issues/661
-    for prob in deepcopy(problems[1:end-1])
+    for prob in deepcopy([oprob, sprob, dprob, nprob])
         # Get u values (including observables).
         @test prob[X] == prob[model.X] == prob[:X] == 4
         @test prob[XY] == prob[model.XY] == prob[:XY] == 9
         @test prob[[XY,Y]] == prob[[model.XY,model.Y]] == prob[[:XY,:Y]] == [9, 5]
-        @test_broken prob[(XY,Y)] == prob[(model.XY,model.Y)] == prob[(:XY,:Y)] == (9, 5)
+        @test prob[(XY,Y)] == prob[(model.XY,model.Y)] == prob[(:XY,:Y)] == (9, 5)
         @test getu(prob, X)(prob) == getu(prob, model.X)(prob) == getu(prob, :X)(prob) == 4
         @test getu(prob, XY)(prob) == getu(prob, model.XY)(prob) == getu(prob, :XY)(prob) == 9 
         @test getu(prob, [XY,Y])(prob) == getu(prob, [model.XY,model.Y])(prob) == getu(prob, [:XY,:Y])(prob) == [9, 5]  
@@ -117,8 +118,10 @@ end
 
 # Test remake function.
 let 
+    @test_broken false # Cannot check result for JumpProblem: https://github.com/SciML/ModelingToolkit.jl/issues/2838
+    @test_broken false # Cannot deepcopy SteadyStateProblem :https://github.com/SciML/ModelingToolkit.jl/issues/2837
     @test_broken false # Currently cannot be run for Ensemble problems: https://github.com/SciML/SciMLBase.jl/issues/661 (as indexing cannot be used to check values).
-    for prob in deepcopy(problems)
+    for prob in deepcopy([oprob, sprob, dprob, nprob])
         # Remake for all u0s.
         rp = remake(prob; u0 = [X => 1, Y => 2])
         @test rp[[X, Y]] == [1, 2]
@@ -155,10 +158,8 @@ end
 
 # Test integrator indexing.
 let 
-    @test_broken false # NOTE: Multiple problems for `nint` (https://github.com/SciML/SciMLBase.jl/issues/662).
-    @test_broken false # NOTE: Multiple problems for `jint` (https://github.com/SciML/SciMLBase.jl/issues/654).
     @test_broken false # NOTE: Cannot even create a `ssint` (https://github.com/SciML/SciMLBase.jl/issues/660).
-    for int in deepcopy([oint, sint])
+    for int in deepcopy([oint, sint, jint, nint])
         # Get u values.
         @test int[X] == int[model.X] == int[:X] == 4
         @test int[XY] == int[model.XY] == int[:XY] == 9
@@ -258,13 +259,12 @@ let
 
     # Handles nonlinear and steady state solutions differently.
     let
-        @test_broken false # Currently a problem for nonlinear solutions and steady state solutions (https://github.com/SciML/SciMLBase.jl/issues/720).
-        for sol in deepcopy([])
+        for sol in deepcopy([nsol, sssol])
             # Get u values.
             @test sol[X] == sol[model.X] == sol[:X]
             @test sol[XY] == sol[model.XY][1] == sol[:XY]
             @test sol[[XY,Y]] == sol[[model.XY,model.Y]] == sol[[:XY,:Y]]
-            @test_broken sol[(XY,Y)] == sol[(model.XY,model.Y)] == sol[(:XY,:Y)]
+            @test sol[(XY,Y)] == sol[(model.XY,model.Y)] == sol[(:XY,:Y)]
             @test getu(sol, X)(sol) == getu(sol, model.X)(sol)[1] == getu(sol, :X)(sol)
             @test getu(sol, XY)(sol) == getu(sol, model.XY)(sol)[1] == getu(sol, :XY)(sol)
             @test getu(sol, [XY,Y])(sol) == getu(sol, [model.XY,model.Y])(sol) == getu(sol, [:XY,:Y])(sol)
@@ -283,8 +283,7 @@ end
 
 # Tests plotting.
 let 
-    @test_broken false # Currently broken for `ssol` (https://github.com/SciML/SciMLBase.jl/issues/580) 
-    for sol in deepcopy([osol, jsol])
+    for sol in deepcopy([osol, jsol, ssol])
         # Single variable.
         @test length(plot(sol; idxs = X).series_list) == 1
         @test length(plot(sol; idxs = XY).series_list) == 1
@@ -389,4 +388,3 @@ let
     reset_aggregated_jumps!(jint)
     @test jint.cb.condition.ma_jumps.scaled_rates[1] == 6.0
 end
-