diff --git a/HISTORY.md b/HISTORY.md
index 753403e7af..ff92e9b5f5 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -4,18 +4,44 @@
 
 ## Catalyst 14.0
 
-#### 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 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 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 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.
+#### 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 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 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]`).
+-  ModelingToolkit currently does not support the safe application of the
+   `remake` function, or safe direct mutation, for 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 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
-- 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. 
+- 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()
@@ -25,110 +51,141 @@ Catalyst v14 was prompted by the (breaking) release of ModelingToolkit v9, which
   @variables t
   @species A(t)
 - 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).
+  ```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)
+- Updated 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 two `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]
+- 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 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.
-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
-goodwind_oscillator = @reaction_network begin
-    (mmr(P,pₘ,1), dₘ), 0 <--> M
-    (pₑ*M,dₑ), 0 <--> E
-    (pₚ*E,dₚ), 0 <--> P
-end
-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).
-- A more detailed of how this extension works can be found [here](https://docs.sciml.ai/Catalyst/stable/inverse_problems/structural_identifiability/).
+  # 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 keyword 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.
+  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
+  ```
+  couples the ODE $dV/dt = 1 - V$ to the reaction system.
+- Coupled reaction networks and differential equation (or algebraic differential
+  equation) systems can now be converted to `SDESystem`s and `NonlinearSystem`s.
 
-#### Bifurcation analysis extension 
-- Add a CatalystBifurcationKitExtension, permitting BifurcationKit's `BifurcationProblem`s to be created from Catalyst reaction networks. Example usage:
-```julia
-using Catalyst
-wilhelm_2009_model = @reaction_network begin
-    k1, Y --> 2X
-    k2, 2X --> X + Y
-    k3, X + Y --> Y
-    k4, X --> 0
-    k5, 0 --> X
-end
+#### Structural identifiability extension
+- Added CatalystStructuralIdentifiabilityExtension, which permits
+  StructuralIdentifiability.jl to be applied directly to Catalyst systems. E.g.
+  use
+  ```julia
+  using Catalyst, StructuralIdentifiability
+  goodwind_oscillator = @reaction_network begin
+      (mmr(P,pₘ,1), dₘ), 0 <--> M
+      (pₑ*M,dₑ), 0 <--> E
+      (pₚ*E,dₚ), 0 <--> P
+  end
+  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).
+- A more detailed of how this extension works can be found
+  [here](https://docs.sciml.ai/Catalyst/stable/inverse_problems/structural_identifiability/).
+
+#### Bifurcation analysis extension
+- Add a CatalystBifurcationKitExtension, permitting BifurcationKit's
+  `BifurcationProblem`s to be created from Catalyst reaction networks. Example
+  usage:
+  ```julia
+  using Catalyst
+  wilhelm_2009_model = @reaction_network begin
+      k1, Y --> 2X
+      k2, 2X --> X + Y
+      k3, X + Y --> Y
+      k4, X --> 0
+      k5, 0 --> X
+  end
 
-using BifurcationKit
-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)
+  using BifurcationKit
+  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)
 
-p_span = (2.0, 20.0)
-opts_br = ContinuationPar(p_min = p_span[1], p_max = p_span[2], max_steps = 1000)
+  p_span = (2.0, 20.0)
+  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", guide = "X")
-```
-- Automatically handles elimination of conservation laws for computing bifurcation diagrams.
+  using Plots
+  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.
 
 ## Catalyst 13.5
diff --git a/Project.toml b/Project.toml
index 0a302b5529..5d2088785a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "14.0.0-DEV"
+version = "14.0.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
diff --git a/README.md b/README.md
index bc431716ca..18a432c9a1 100644
--- a/README.md
+++ b/README.md
@@ -69,9 +69,9 @@ be found in its corresponding research paper, [Catalyst: Fast and flexible model
 - Support for [parallelization of all simulations](https://docs.sciml.ai/Catalyst/stable/model_simulation/ode_simulation_performance/#ode_simulation_performance_parallelisation), including parallelization of [ODE simulations on GPUs](https://docs.sciml.ai/Catalyst/stable/model_simulation/ode_simulation_performance/#ode_simulation_performance_parallelisation_GPU) using [DiffEqGPU.jl](https://github.com/SciML/DiffEqGPU.jl).
 - [Latexify](https://korsbo.github.io/Latexify.jl/stable/) can be used to [generate LaTeX expressions](https://docs.sciml.ai/Catalyst/stable/model_creation/model_visualisation/#visualisation_latex) corresponding to generated mathematical models or the underlying set of reactions.
 - [Graphviz](https://graphviz.org/) can be used to generate and [visualize reaction network graphs](https://docs.sciml.ai/Catalyst/stable/model_creation/model_visualisation/#visualisation_graphs) (reusing the Graphviz interface created in [Catlab.jl](https://algebraicjulia.github.io/Catlab.jl/stable/)).
-- Model steady states can be computed through homotopy continuation using [HomotopyContinuation.jl](https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl) (which can find *all* steady states of systems with multiple ones), by forward ODE simulations using [SteadyStateDiffEq.jl](https://github.com/SciML/SteadyStateDiffEq.jl), or by numerically solving steady-state nonlinear equations using [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl).
+- Model steady states can be [computed through homotopy continuation](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/homotopy_continuation/) using [HomotopyContinuation.jl](https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl) (which can find *all* steady states of systems with multiple ones), by [forward ODE simulations](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/nonlinear_solve/#steady_state_solving_simulation) using [SteadyStateDiffEq.jl](https://github.com/SciML/SteadyStateDiffEq.jl), or by [numerically solving steady-state nonlinear equations](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/nonlinear_solve/#steady_state_solving_nonlinear) using [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl).
 - [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) can be used to [compute bifurcation diagrams](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/bifurcation_diagrams/) of model steady states (including finding periodic orbits).
-- [DynamicalSystems.jl](https://github.com/JuliaDynamics/DynamicalSystems.jl) can be used to compute model [basins of attraction](@ref dynamical_systems_basins_of_attraction), [Lyapunov spectrums](@ref dynamical_systems_lyapunov_exponents), and other dynamical system properties.
+- [DynamicalSystems.jl](https://github.com/JuliaDynamics/DynamicalSystems.jl) can be used to compute model [basins of attraction]([@ref dynamical_systems_basins_of_attraction](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/dynamical_systems/#dynamical_systems_basins_of_attraction)), [Lyapunov spectrums](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/dynamical_systems/#dynamical_systems_lyapunov_exponents), and other dynamical system properties.
 - [StructuralIdentifiability.jl](https://github.com/SciML/StructuralIdentifiability.jl) can be used to [perform structural identifiability analysis](https://docs.sciml.ai/Catalyst/stable/inverse_problems/structural_identifiability/).
 - [Optimization.jl](https://github.com/SciML/Optimization.jl), [DiffEqParamEstim.jl](https://github.com/SciML/DiffEqParamEstim.jl), and [PEtab.jl](https://github.com/sebapersson/PEtab.jl) can all be used to [fit model parameters to data](https://sebapersson.github.io/PEtab.jl/stable/Define_in_julia/).
 - [GlobalSensitivity.jl](https://github.com/SciML/GlobalSensitivity.jl) can be used to perform [global sensitivity analysis](https://docs.sciml.ai/Catalyst/stable/inverse_problems/global_sensitivity_analysis/) of model behaviors.
diff --git a/docs/src/api.md b/docs/src/api.md
index 8a5e543973..459dbd1c1b 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -77,6 +77,7 @@ plot(p1, p2, p3; layout = (3,1))
 
 ```@docs
 @reaction_network
+@network_component
 make_empty_network
 @reaction
 Reaction
diff --git a/docs/src/index.md b/docs/src/index.md
index a8ea0531da..7ba02ee07f 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -38,7 +38,7 @@ etc).
 - Support for [parallelization of all simulations](@ref ode_simulation_performance_parallelisation), including parallelization of [ODE simulations on GPUs](@ref ode_simulation_performance_parallelisation_GPU) using [DiffEqGPU.jl](https://github.com/SciML/DiffEqGPU.jl).
 - [Latexify](https://korsbo.github.io/Latexify.jl/stable/) can be used to [generate LaTeX expressions](@ref visualisation_latex) corresponding to generated mathematical models or the underlying set of reactions.
 - [Graphviz](https://graphviz.org/) can be used to generate and [visualize reaction network graphs](@ref visualisation_graphs) (reusing the Graphviz interface created in [Catlab.jl](https://algebraicjulia.github.io/Catlab.jl/stable/)).
-- Model steady states can be [computed through homotopy continuation](@ref homotopy_continuation) using [HomotopyContinuation.jl](https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl) (which can find *all* steady states of systems with multiple ones), by [forward ODE simulations](@ref steady_state_solving_simulation) using [SteadyStateDiffEq.jl)](https://github.com/SciML/SteadyStateDiffEq.jl, or by [numerically solving steady-state nonlinear equations](@ref steady_state_solving_nonlinear) using [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl).
+- Model steady states can be [computed through homotopy continuation](@ref homotopy_continuation) using [HomotopyContinuation.jl](https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl) (which can find *all* steady states of systems with multiple ones), by [forward ODE simulations](@ref steady_state_solving_simulation) using [SteadyStateDiffEq.jl)](https://github.com/SciML/SteadyStateDiffEq.jl), or by [numerically solving steady-state nonlinear equations](@ref steady_state_solving_nonlinear) using [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl).
 - [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) can be used to [compute bifurcation diagrams](@ref bifurcation_diagrams) of model steady states (including finding periodic orbits).
 - [DynamicalSystems.jl](https://github.com/JuliaDynamics/DynamicalSystems.jl) can be used to compute model [basins of attraction](@ref dynamical_systems_basins_of_attraction), [Lyapunov spectrums](@ref dynamical_systems_lyapunov_exponents), and other dynamical system properties.
 - [StructuralIdentifiability.jl](https://github.com/SciML/StructuralIdentifiability.jl) can be used to [perform structural identifiability analysis](@ref structural_identifiability).
diff --git a/docs/src/model_simulation/simulation_structure_interfacing.md b/docs/src/model_simulation/simulation_structure_interfacing.md
index 4f44863036..1c37df3672 100644
--- a/docs/src/model_simulation/simulation_structure_interfacing.md
+++ b/docs/src/model_simulation/simulation_structure_interfacing.md
@@ -1,7 +1,7 @@
 # [Interfacing problems, integrators, and solutions](@id simulation_structure_interfacing)
 When simulating a model, one begins with creating a [problem](https://docs.sciml.ai/DiffEqDocs/stable/basics/problem/). Next, a simulation is performed on the problem, during which the simulation's state is recorded through an [integrator](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/). Finally, the simulation output is returned as a [solution](https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/). This tutorial describes how to access (or modify) the state (or parameter) values of problem, integrator, and solution structures.
 
-Generally, when we have a structure `simulation_struct` and want to interface with the unknown (or parameter) `x`, we use `simulation_struct[:x]` to access the value, and `simulation_struct[:x] = 5.0` to set it to a new value. For situations where a value is accessed (or changed) a large number of times, it can *improve performance* to first create a [specialised getter/setter function](@ref simulation_structure_interfacing_functions).
+Generally, when we have a structure `simulation_struct` and want to interface with the unknown (or parameter) `x`, we use `simulation_struct[:x]` to access the value. For situations where a value is accessed (or changed) a large number of times, it can *improve performance* to first create a [specialised getter/setter function](@ref simulation_structure_interfacing_functions).
 
 ## [Interfacing problem objects](@id simulation_structure_interfacing_problems)
 
@@ -35,26 +35,6 @@ To retrieve several species initial condition (or parameter) values, simply give
 oprob[[:S₁, :S₂]]
 ```
 
-We can change a species's initial condition value using a similar notation. Here we increase the initial concentration of $C$ (and also confirm that the new value is stored in an updated `oprob`):
-```@example structure_indexing
-oprob[:C] = 0.1
-oprob[:C]
-```
-Again, parameter values can be changed using a similar notation, however, again requiring `oprob.ps` notation:
-```@example structure_indexing
-oprob.ps[:k₁] = 10.0
-oprob.ps[:k₁]
-```
-Finally, vectors can be used to update multiple quantities simultaneously
-```@example structure_indexing
-oprob[[:S₁, :S₂]] = [0.5, 0.3]
-oprob[[:S₁, :S₂]]
-```
-Generally, when updating problems, it is often better to use the [`remake` function](@ref simulation_structure_interfacing_problems_remake) (especially when several values are updated).
-
-!!! warning
-    Indexing *should not* be used not modify `JumpProblem`s. Here, [remake](@ref simulation_structure_interfacing_problems_remake) should be used exclusively.
-
 A problem's time span can be accessed through the `tspan` field:
 ```@example structure_indexing
 oprob.tspan
@@ -64,8 +44,8 @@ oprob.tspan
     Here we have used an `ODEProblem`to demonstrate all interfacing functionality. However, identical workflows work for the other problem types.
 
 ### [Remaking problems using the `remake` function](@id simulation_structure_interfacing_problems_remake)
-The `remake` function offers an (to indexing) alternative approach for updating problems. Unlike indexing, `remake` creates a new problem (rather than updating the old one). Furthermore, it permits the updating of several values simultaneously. The `remake` function takes the following inputs:
-- The problem that is remakes.
+To modify a problem, the `remake` function should be used. It takes an already created problem, and returns a new, updated, one (the input problem is unchanged). The `remake` function takes the following inputs:
+- The problem that it remakes.
 - (optionally) `u0`: A vector with initial conditions that should be updated. The vector takes the same form as normal initial condition vectors, but does not need to be complete (in which case only a subset of the initial conditions are updated).
 - (optionally) `tspan`: An updated time span (using the same format as time spans normally are given in).
 - (optionally) `p`: A vector with parameters that should be updated. The vector takes the same form as normal parameter vectors, but does not need to be complete (in which case only a subset of the parameters are updated).
@@ -74,22 +54,28 @@ Here we modify our problem to increase the initial condition concentrations of t
 ```@example structure_indexing
 using OrdinaryDiffEq
 oprob_new = remake(oprob; u0 = [:S₁ => 5.0, :S₂ => 2.5])
-oprob_new == oprob
+oprob_new != oprob
 ```
-Here, we instead use `remake` to simultaneously update a all three fields:
+Here, we instead use `remake` to simultaneously update all three fields:
 ```@example structure_indexing
 oprob_new_2 = remake(oprob; u0 = [:C => 0.2], tspan = (0.0, 20.0), p = [:k₁ => 2.0, :k₂ => 2.0])
 nothing # hide
 ```
 
+Typically, when using `remake` to update a problem, the common workflow is to overwrite the old one with the output. E.g. to set the value of `k₁` to `5.0` in `oprob`, you would do:
+```@example structure_indexing
+oprob = remake(oprob; p = [:k₁ => 5.0])
+nothing # hide
+```
+
 ## [Interfacing integrator objects](@id simulation_structure_interfacing_integrators)
 
-During a simulation, the solution is stored in an integrator object. Here, we will describe how to interface with these. The almost exclusive circumstance when integrator-interfacing is relevant is when simulation events are implemented through callbacks. However, to demonstrate integrator indexing in this tutorial, we will create one through the `init` function (while circumstances where one might [want to use `init` function exist](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/#Initialization-and-Stepping), since integrators are automatically created during simulations, these are rare).
+During a simulation, the solution is stored in an integrator object. Here, we will describe how to interface with these. The almost exclusive circumstance when integrator-interfacing is relevant is when simulation events are implemented through [callbacks](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/). However, to demonstrate integrator indexing in this tutorial, we will create one through the `init` function (while circumstances where one might [want to use `init` function exist](https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/#Initialization-and-Stepping), since integrators are automatically created during simulations, these are rare).
 ```@example structure_indexing
 integrator = init(oprob)
 nothing # hide
 ```
-We can interface with our integrator using an identical syntax as [was used for problems](@ref simulation_structure_interfacing_problems) (with the exception that `remake` is not available). Here we update, and then check the values of, first the species $C$ and then the parameter $k₁$:
+We can interface with our integrator using an identical syntax as [was used for problems](@ref simulation_structure_interfacing_problems). The primary exception is that there is no `remake` function for integrators. Instead, we can update species and parameter values using normal indexing. Here we update, and then check the values of, first the species $C$ and then the parameter $k₁$:
 ```@example structure_indexing
 integrator[:C] = 0.0
 integrator[:C]
@@ -164,7 +150,7 @@ get_S(oprob)
 ## [Interfacing using symbolic representations](@id simulation_structure_interfacing_symbolic_representation)
 When e.g. [programmatic modelling is used](@ref programmatic_CRN_construction), species and parameters can be represented as *symbolic variables*. These can be used to index a problem, just like symbol-based representations can. Here we create a simple [two-state model](@ref basic_CRN_library_two_states) programmatically, and use its symbolic variables to check, and update, an initial condition:
 ```@example structure_indexing_symbolic_variables
-using Catalyst
+using Catalyst, OrdinaryDiffEq
 t = default_t()
 @species X1(t) X2(t)
 @parameters k1 k2
@@ -180,7 +166,7 @@ tspan = (0.0, 1.0)
 ps = [k1 => 1.0, k2 => 2.0]
 oprob = ODEProblem(two_state_model, u0, tspan, ps)
 
-oprob[X1] = 5.0
+oprob = remake(oprob; u0 = [X1 => 5.0])
 oprob[X1]
 ```
 Symbolic variables can be used to access or update species or parameters for all the cases when `Symbol`s can (including when using `remake` or e.g. `getu`).
@@ -197,4 +183,4 @@ oprob[two_state_model.X1 + two_state_model.X2]
 This can be used to form symbolic expressions using model quantities when a model has been created using the DSL (as an alternative to @unpack). Alternatively, [creating an observable](@ref dsl_advanced_options_observables), and then interface using its `Symbol` representation, is also possible.
 
 !!! warning
-    With interfacing with a simulating structure using symbolic variables stored in a `ReactionSystem` model, ensure that the model is complete.
\ No newline at end of file
+    When accessing a simulation structure using symbolic variables from a `ReactionSystem` model, such as `rn.A` for `rn` a `ReactionSystem` and `A` a species within it, ensure that the model is complete.
diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index ed47a2a439..ed5074839b 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 *not marked as 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
@@ -64,8 +64,30 @@ sol = solve(oprob)
 plot(sol)
 ```
 
+Note, if we had instead used the [`@reaction_network`](@ref) DSL macro to build
+our model, i.e.
+```@example v14_migration_1
+rs2 = @reaction_network rs begin
+    p, ∅ --> X
+    d, X --> ∅
+end
+```
+then the model is automatically marked as complete
+```@example v14_migration_1
+Catalyst.iscomplete(rs2)
+```
+In contrast, if we used the [`@network_component`](@ref) DSL macro to build our
+model it is not marked as complete, and is equivalent to our original definition of `rs`
+```@example v14_migration_1
+rs3 = @network_component rs begin
+    p, ∅ --> X
+    d, X --> ∅
+end
+Catalyst.iscomplete(rs3)
+```
+
 ## Unknowns instead of states
-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`). 
+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
@@ -100,10 +122,10 @@ As part of its v9 update, ModelingToolkit changed how units were handled. This i
 
 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.
+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 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.
+According to the ModelingToolkit system API, systems should not be mutable. In accordance with this, the following functions have been deprecated and removed: `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.
 
@@ -132,9 +154,9 @@ nothing # hide
 ```
 
 !!! 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.
+    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 the functions `default_t()` and `default_time_deriv()` as our main approach.
 
-## New interface for accessing problem/integrator/solution parameter (and species) value
+## New interface for accessing problem/integrator/solution parameter (and species) values
 Previously, it was possible to directly index problems to query them for their parameter values. e.g.
 ```@example v14_migration_4
 using Catalyst
@@ -154,12 +176,11 @@ This is *no longer supported*. When you wish to query a problem (or integrator o
 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 described in more detail [here](@ref simulation_structure_interfacing_functions).
+Furthermore, a few new functions (`getp`, `getu`, `setp`, `setu`) have been introduced from [SymbolicIndexingInterface](https://github.com/SciML/SymbolicIndexingInterface.jl) to support efficient and systematic querying and/or updating of symbolic unknown/parameter values. Using these can *significantly* improve performance when querying or updating a value multiple times, for example within a callback. 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.
+## Other changes
 
 #### 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
@@ -174,7 +195,7 @@ oprob = ODEProblem(rn, u0, (0.0, 10.0), ps; remove_conserved = true)
 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:
+is perfectly fine, attempting to then modify any initial conditions or the value of the conservation constant in `oprob` will likely silently fail:
 ```@example v14_migration_5
 oprob_remade = remake(oprob; u0 = [:X1 => 5.0]) # NEVER do this.
 solve(oprob_remade)
diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index 7d40cf695c..584ca07db8 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -127,7 +127,8 @@ export @reaction_network, @network_component, @reaction, @species
 include("network_analysis.jl")
 export reactioncomplexmap, reactioncomplexes, incidencemat
 export complexstoichmat
-export complexoutgoingmat, incidencematgraph, linkageclasses, deficiency, subnetworks
+export complexoutgoingmat, incidencematgraph, linkageclasses, stronglinkageclasses,
+       terminallinkageclasses, deficiency, subnetworks
 export linkagedeficiencies, isreversible, isweaklyreversible
 export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants
 
diff --git a/src/dsl.jl b/src/dsl.jl
index 633ea7f64b..6148eb4814 100644
--- a/src/dsl.jl
+++ b/src/dsl.jl
@@ -175,7 +175,8 @@ end
 """
     @network_component
 
-As @reaction_network, but the output system is not complete.
+Equivalent to `@reaction_network` except the generated `ReactionSystem` is not marked as
+complete.
 """
 macro network_component(name::Symbol, ex::Expr)
     make_reaction_system(MacroTools.striplines(ex); name = :($(QuoteNode(name))))
diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index 81208e1bf7..1b0a18c1bf 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -260,25 +260,19 @@ end
 Construct a directed simple graph where nodes correspond to reaction complexes and directed
 edges to reactions converting between two complexes.
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-complexes,incidencemat = reactioncomplexes(sir)
 incidencematgraph(sir)
 ```
 """
 function incidencematgraph(rn::ReactionSystem)
     nps = get_networkproperties(rn)
     if Graphs.nv(nps.incidencegraph) == 0
-        isempty(nps.incidencemat) &&
-            error("Please call reactioncomplexes(rn) first to construct the incidence matrix.")
+        isempty(nps.incidencemat) && reactioncomplexes(rn)
         nps.incidencegraph = incidencematgraph(nps.incidencemat)
     end
     nps.incidencegraph
@@ -329,17 +323,12 @@ Given the incidence graph of a reaction network, return a vector of the
 connected components of the graph (i.e. sub-groups of reaction complexes that
 are connected in the incidence graph).
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-complexes,incidencemat = reactioncomplexes(sir)
 linkageclasses(sir)
 ```
 gives
@@ -359,6 +348,56 @@ end
 
 linkageclasses(incidencegraph) = Graphs.connected_components(incidencegraph)
 
+"""
+    stronglinkageclasses(rn::ReactionSystem)
+
+    Return the strongly connected components of a reaction network's incidence graph (i.e. sub-groups of reaction complexes such that every complex is reachable from every other one in the sub-group).
+"""
+
+function stronglinkageclasses(rn::ReactionSystem)
+    nps = get_networkproperties(rn)
+    if isempty(nps.stronglinkageclasses)
+        nps.stronglinkageclasses = stronglinkageclasses(incidencematgraph(rn))
+    end
+    nps.stronglinkageclasses
+end
+
+stronglinkageclasses(incidencegraph) = Graphs.strongly_connected_components(incidencegraph)
+
+"""
+    terminallinkageclasses(rn::ReactionSystem)
+
+    Return the terminal strongly connected components of a reaction network's incidence graph (i.e. sub-groups of reaction complexes that are 1) strongly connected and 2) every outgoing reaction from a complex in the component produces a complex also in the component).
+"""
+
+function terminallinkageclasses(rn::ReactionSystem)
+    nps = get_networkproperties(rn)
+    if isempty(nps.terminallinkageclasses)
+        slcs = stronglinkageclasses(rn)
+        tslcs = filter(lc -> isterminal(lc, rn), slcs)
+        nps.terminallinkageclasses = tslcs
+    end
+    nps.terminallinkageclasses
+end
+
+# Helper function for terminallinkageclasses. Given a linkage class and a reaction network, say whether the linkage class is terminal, 
+# i.e. all outgoing reactions from complexes in the linkage class produce a complex also in the linkage class
+function isterminal(lc::Vector, rn::ReactionSystem)
+    imat = incidencemat(rn)
+
+    for r in 1:size(imat, 2)
+        # Find the index of the reactant complex for a given reaction
+        s = findfirst(==(-1), @view imat[:, r])
+
+        # If the reactant complex is in the linkage class, check whether the product complex is also in the linkage class. If any of them are not, return false. 
+        if s in Set(lc)
+            p = findfirst(==(1), @view imat[:, r])
+            p in Set(lc) ? continue : return false
+        end
+    end
+    true
+end
+
 @doc raw"""
     deficiency(rn::ReactionSystem)
 
@@ -371,17 +410,12 @@ Here the deficiency, ``\delta``, of a network with ``n`` reaction complexes,
 \delta = n - \ell - s
 ```
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-rcs,incidencemat = reactioncomplexes(sir)
 δ = deficiency(sir)
 ```
 """
@@ -419,17 +453,12 @@ end
 
 Find subnetworks corresponding to each linkage class of the reaction network.
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-complexes,incidencemat = reactioncomplexes(sir)
 subnetworks(sir)
 ```
 """
@@ -456,17 +485,12 @@ end
 
 Calculates the deficiency of each sub-reaction network within `network`.
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-rcs,incidencemat = reactioncomplexes(sir)
 linkage_deficiencies = linkagedeficiencies(sir)
 ```
 """
@@ -487,17 +511,12 @@ end
 
 Given a reaction network, returns if the network is reversible or not.
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-rcs,incidencemat = reactioncomplexes(sir)
 isreversible(sir)
 ```
 """
@@ -511,30 +530,27 @@ end
 
 Determine if the reaction network with the given subnetworks is weakly reversible or not.
 
-Notes:
-- Requires the `incidencemat` to already be cached in `rn` by a previous call to
-  `reactioncomplexes`.
-
 For example,
 ```julia
 sir = @reaction_network SIR begin
     β, S + I --> 2I
     ν, I --> R
 end
-rcs,incidencemat = reactioncomplexes(sir)
 subnets = subnetworks(rn)
 isweaklyreversible(rn, subnets)
 ```
 """
 function isweaklyreversible(rn::ReactionSystem, subnets)
-    im = get_networkproperties(rn).incidencemat
-    isempty(im) &&
-        error("Error, please call reactioncomplexes(rn::ReactionSystem) to ensure the incidence matrix has been cached.")
-    sparseig = issparse(im)
+    nps = get_networkproperties(rn)
+    isempty(nps.incidencemat) && reactioncomplexes(rn)
+    sparseig = issparse(nps.incidencemat)
+
     for subnet in subnets
-        nps = get_networkproperties(subnet)
-        isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig)
+        subnps = get_networkproperties(subnet)
+        isempty(subnps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig)
     end
+
+    # A network is weakly reversible if all of its subnetworks are strongly connected
     all(Graphs.is_strongly_connected ∘ incidencematgraph, subnets)
 end
 
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index b504710eee..7bd09f555b 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -78,6 +78,7 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     isempty::Bool = true
     netstoichmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
     conservationmat::Matrix{I} = Matrix{I}(undef, 0, 0)
+    cyclemat::Matrix{I} = Matrix{I}(undef, 0, 0)
     col_order::Vector{Int} = Int[]
     rank::Int = 0
     nullity::Int = 0
@@ -93,6 +94,8 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     complexoutgoingmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
     incidencegraph::Graphs.SimpleDiGraph{Int} = Graphs.DiGraph()
     linkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
+    stronglinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
+    terminallinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
     deficiency::Int = 0
 end
 #! format: on
@@ -116,6 +119,7 @@ function reset!(nps::NetworkProperties{I, V}) where {I, V}
     nps.isempty && return
     nps.netstoichmat = Matrix{Int}(undef, 0, 0)
     nps.conservationmat = Matrix{I}(undef, 0, 0)
+    nps.cyclemat = Matrix{Int}(undef, 0, 0)
     empty!(nps.col_order)
     nps.rank = 0
     nps.nullity = 0
@@ -131,6 +135,8 @@ function reset!(nps::NetworkProperties{I, V}) where {I, V}
     nps.complexoutgoingmat = Matrix{Int}(undef, 0, 0)
     nps.incidencegraph = Graphs.DiGraph()
     empty!(nps.linkageclasses)
+    empty!(nps.stronglinkageclasses)
+    empty!(nps.terminallinkageclasses)
     nps.deficiency = 0
 
     # this needs to be last due to setproperty! setting it to false
diff --git a/test/network_analysis/network_properties.jl b/test/network_analysis/network_properties.jl
index 5507a9f655..811d6418cc 100644
--- a/test/network_analysis/network_properties.jl
+++ b/test/network_analysis/network_properties.jl
@@ -326,3 +326,86 @@ let
     @test Catalyst.iscomplexbalanced(rn, rates) == true 
 end
 
+### STRONG LINKAGE CLASS TESTS
+
+
+# a) Checks that strong/terminal linkage classes are correctly found. Should identify the (A, B+C) linkage class as non-terminal, since B + C produces D
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B + C
+        k3, B + C --> D
+        k4, D --> E
+        (k5, k6), E <--> 2F
+        k7, 2F --> D
+        (k8, k9), D + E <--> G
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    slcs = stronglinkageclasses(rn)
+    tslcs = terminallinkageclasses(rn)
+    @test length(slcs) == 3
+    @test length(tslcs) == 2
+    @test issubset([[1,2], [3,4,5], [6,7]], slcs)
+    @test issubset([[3,4,5], [6,7]], tslcs) 
+end
+
+# b) Makes the D + E --> G reaction irreversible. Thus, (D+E) becomes a non-terminal linkage class. Checks whether correctly identifies both (A, B+C) and (D+E) as non-terminal 
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B + C
+        k3, B + C --> D
+        k4, D --> E
+        (k5, k6), E <--> 2F
+        k7, 2F --> D
+        (k8, k9), D + E --> G
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    slcs = stronglinkageclasses(rn)
+    tslcs = terminallinkageclasses(rn)
+    @test length(slcs) == 4
+    @test length(tslcs) == 2
+    @test issubset([[1,2], [3,4,5], [6], [7]], slcs)
+    @test issubset([[3,4,5], [7]], tslcs) 
+end
+
+# From a), makes the B + C <--> D reaction reversible. Thus, the non-terminal (A, B+C) linkage class gets absorbed into the terminal (A, B+C, D, E, 2F) linkage class, and the terminal linkage classes and strong linkage classes coincide. 
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> B + C
+        (k3, k4), B + C <--> D
+        k5, D --> E
+        (k6, k7), E <--> 2F
+        k8, 2F --> D
+        (k9, k10), D + E <--> G
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    slcs = stronglinkageclasses(rn)
+    tslcs = terminallinkageclasses(rn)
+    @test length(slcs) == 2
+    @test length(tslcs) == 2
+    @test issubset([[1,2,3,4,5], [6,7]], slcs)
+    @test issubset([[1,2,3,4,5], [6,7]], tslcs) 
+end
+
+# Simple test for strong and terminal linkage classes
+let
+    rn = @reaction_network begin
+        (k1, k2), A <--> 2B 
+        k3, A --> C + D
+        (k4, k5), C + D <--> E
+        k6, 2B --> F
+        (k7, k8), F <--> 2G
+        (k9, k10), 2G <--> H
+        k11, H --> F
+    end
+
+    rcs, D = reactioncomplexes(rn)
+    slcs = stronglinkageclasses(rn)
+    tslcs = terminallinkageclasses(rn)
+    @test length(slcs) == 3
+    @test length(tslcs) == 2
+    @test issubset([[1,2], [3,4], [5,6,7]], slcs)
+    @test issubset([[3,4], [5,6,7]], tslcs) 
+end
diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl
index 514dbc136f..538ce0f174 100644
--- a/test/reactionsystem_core/reactionsystem.jl
+++ b/test/reactionsystem_core/reactionsystem.jl
@@ -816,3 +816,42 @@ end
 let
     @test_nowarn Catalyst.reactionsystem_uptodate_check()
 end
+
+# Test that functions using the incidence matrix properly cache it
+let
+    rn = @reaction_network begin
+        k1, A --> B
+        k2, B --> C
+        k3, C --> A
+    end
+
+    nps = Catalyst.get_networkproperties(rn)
+    @test isempty(nps.incidencemat) == true
+
+    img = incidencematgraph(rn)
+    @test size(nps.incidencemat) == (3,3)
+
+    Catalyst.reset!(nps)
+    lcs = linkageclasses(rn)
+    @test size(nps.incidencemat) == (3,3)
+
+    Catalyst.reset!(nps)
+    sns = subnetworks(rn)
+    @test size(nps.incidencemat) == (3,3)
+
+    Catalyst.reset!(nps)
+    δ = deficiency(rn)
+    @test size(nps.incidencemat) == (3,3)
+    
+    Catalyst.reset!(nps)
+    δ_l = linkagedeficiencies(rn)
+    @test size(nps.incidencemat) == (3,3)
+
+    Catalyst.reset!(nps)
+    rev = isreversible(rn)
+    @test size(nps.incidencemat) == (3,3)
+
+    Catalyst.reset!(nps)
+    weakrev = isweaklyreversible(rn, sns)
+    @test size(nps.incidencemat) == (3,3)
+end