diff --git a/docs/Project.toml b/docs/Project.toml
index cf248f2507..e7454091d2 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -35,7 +35,7 @@ DataFrames = "1"
 DiffEqParamEstim = "2.1"
 DifferentialEquations = "7.7"
 Distributions = "0.25"
-Documenter = "0.27"
+Documenter = "1.4.1"
 HomotopyContinuation = "2.6"
 Latexify = "0.15, 0.16"
 ModelingToolkit = "9.5"
diff --git a/docs/src/inverse_problems/petab_ode_param_fitting.md b/docs/old_files/petab_ode_param_fitting.md
similarity index 100%
rename from docs/src/inverse_problems/petab_ode_param_fitting.md
rename to docs/old_files/petab_ode_param_fitting.md
diff --git a/docs/pages.jl b/docs/pages.jl
index 04f255eede..faf19c5cac 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -9,10 +9,10 @@ pages = Any[
         "model_creation/dsl_description.md",
         # DSL - Advanced
         "model_creation/programmatic_CRN_construction.md",
-        "model_creation/constraint_equations.md",
         "model_creation/compositional_modeling.md",
+        "model_creation/constraint_equations.md",
         # Events.
-        # Distributed parameters, rates, and initial conditions.
+        "model_creation/parametric_stoichiometry.md",# Distributed parameters, rates, and initial conditions.
         # Loading and writing models to files.
         # Model visualisation.
         "model_creation/network_analysis.md",
@@ -27,7 +27,7 @@ pages = Any[
         # Simulation introduction.
         # Simulation plotting.
         "model_simulation/simulation_structure_interfacing.md",
-        # Monte Carlo/Ensemble simulations.
+        "model_simulation/TOBEREMOVED_advanced_simulations.md",# Monte Carlo/Ensemble simulations.        
         # Stochastic simulation statistical analysis.
         # ODE Performance considerations/advice.
         # SDE Performance considerations/advice.
@@ -48,7 +48,7 @@ pages = Any[
         # ODE parameter fitting using Turing.
         # SDE/Jump fitting.
         # Non-parameter fitting optimisation.
-        "inverse_problems/07_structural_identifiability.md",
+        "inverse_problems/structural_identifiability.md",
         # Practical identifiability.
         # GLobal and local sensitivity analysis.
         "Inverse problem examples" => Any[
diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml
new file mode 100644
index 0000000000..e7454091d2
--- /dev/null
+++ b/docs/src/assets/Project.toml
@@ -0,0 +1,57 @@
+[deps]
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
+Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DiffEqParamEstim = "1130ab10-4a5a-5621-a13d-e4788d82bd4c"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
+Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
+OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
+OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
+Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+
+[compat]
+BifurcationKit = "0.3"
+Catalyst = "13"
+DataFrames = "1"
+DiffEqParamEstim = "2.1"
+DifferentialEquations = "7.7"
+Distributions = "0.25"
+Documenter = "1.4.1"
+HomotopyContinuation = "2.6"
+Latexify = "0.15, 0.16"
+ModelingToolkit = "9.5"
+NonlinearSolve = "3.4.0"
+Optim = "1"
+Optimization = "3.19"
+OptimizationNLopt = "0.1.8"
+OptimizationOptimJL = "0.1.14"
+OptimizationOptimisers = "0.1.1"
+OrdinaryDiffEq = "6"
+Plots = "1.36"
+SciMLBase = "2.13"
+SciMLSensitivity = "7.19"
+Setfield = "1.1"
+SpecialFunctions = "2.1"
+SteadyStateDiffEq = "2.0.1"
+StochasticDiffEq = "6"
+StructuralIdentifiability = "0.5.1"
+Symbolics = "5.14"
diff --git a/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md b/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md
index e1e8af412c..5e23f8235f 100644
--- a/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md
+++ b/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md
@@ -67,7 +67,7 @@ using Plots
 ```
 Here, if we restart Julia, these commands *do need to be rerun.
 
-A more comprehensive (but still short) introduction to package management in Julia is provided at [the end of this documentation page](catalyst_for_new_julia_users_packages). It contains some useful information and is hence highly recommended reading. For a more detailed introduction to Julia package management, please read [the Pkg documentation](https://docs.julialang.org/en/v1/stdlib/Pkg/).
+A more comprehensive (but still short) introduction to package management in Julia is provided at [the end of this documentation page](@ref catalyst_for_new_julia_users_packages). It contains some useful information and is hence highly recommended reading. For a more detailed introduction to Julia package management, please read [the Pkg documentation](https://docs.julialang.org/en/v1/stdlib/Pkg/).
 
 ## Simulating a basic Catalyst model
 Now that we have some basic familiarity with Julia, and have installed and imported the required packages, we will create and simulate a basic chemical reaction network model using Catalyst.
diff --git a/docs/src/model_creation/dsl_description.md b/docs/src/model_creation/dsl_description.md
index 555a107803..c888b4dda9 100644
--- a/docs/src/model_creation/dsl_description.md
+++ b/docs/src/model_creation/dsl_description.md
@@ -41,10 +41,11 @@ variables in Julia.
 The chemical reaction model is generated by the `@reaction_network` macro and
 stored in the `rn` variable (a normal Julia variable, which does not need to be
 called `rn`). It corresponds to a [`ReactionSystem`](@ref), a symbolic
-representation of the chemical network.  The generated `ReactionSystem` can be
+representation of the chemical network. The generated `ReactionSystem` can be
 converted to a symbolic differential equation model via
 ```@example tut2
-osys  = convert(ODESystem, rn)
+osys = convert(ODESystem, rn)
+osys = complete(osys)
 ```
 
 We can then convert the symbolic ODE model into a compiled, optimized
@@ -670,7 +671,8 @@ sol[:Xtot]
 ```
 similarly, we can plot the values of $Xtot$ and $Ytot$ using
 ```@example obs1
-plot(sol; idxs=[:Xtot, :Ytot], label=["Total X" "Total Y"])
+using Plots
+plot(sol; idxs = [rn.Xtot, rn.Ytot], label = ["Total X" "Total Y"])
 ```
 
 If we only wish to provide a single observable, the `begin ... end` block is note required. E.g., to record only the total amount of $X$ we can use:
diff --git a/docs/src/model_creation/example/basic_CRN_examples.md b/docs/src/model_creation/examples/basic_CRN_examples.md
similarity index 100%
rename from docs/src/model_creation/example/basic_CRN_examples.md
rename to docs/src/model_creation/examples/basic_CRN_examples.md
diff --git a/docs/src/model_creation/example/hodgkin_huxley_equation.md b/docs/src/model_creation/examples/hodgkin_huxley_equation.md
similarity index 100%
rename from docs/src/model_creation/example/hodgkin_huxley_equation.md
rename to docs/src/model_creation/examples/hodgkin_huxley_equation.md
diff --git a/docs/src/model_creation/example/smoluchowski_coagulation_equation.md b/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md
similarity index 100%
rename from docs/src/model_creation/example/smoluchowski_coagulation_equation.md
rename to docs/src/model_creation/examples/smoluchowski_coagulation_equation.md
diff --git a/docs/src/model_creation/TOBEREMOVED_parametric_stoichiometry.md b/docs/src/model_creation/parametric_stoichiometry.md
similarity index 100%
rename from docs/src/model_creation/TOBEREMOVED_parametric_stoichiometry.md
rename to docs/src/model_creation/parametric_stoichiometry.md
diff --git a/docs/src/model_simulation/TOBEREMOVED_advanced_simulations.md b/docs/src/model_simulation/TOBEREMOVED_advanced_simulations.md
index 2df209f660..55f4fd8758 100644
--- a/docs/src/model_simulation/TOBEREMOVED_advanced_simulations.md
+++ b/docs/src/model_simulation/TOBEREMOVED_advanced_simulations.md
@@ -177,7 +177,7 @@ Here, we access the system's unknown `X` through the `integrator`, and add
 `5.0` to its amount. We can now simulate our system using the
 callback:
 ```@example ex2
-sol = solve(oprob; callback = ps_cb)
+sol = solve(oprob, Tsit5(); callback = ps_cb)
 plot(sol)
 ```
 
@@ -194,16 +194,16 @@ p = [:k => 1.0]
 oprob = ODEProblem(rn, u0, tspan, p)
 
 condition = [5.0]
-affect!(integrator) = integrator[:k] = 5.0
+affect!(integrator) = integrator.ps[:k] = 5.0
 ps_cb = PresetTimeCallback(condition, affect!)
 
-sol = solve(oprob; callback = ps_cb)
+sol = solve(oprob, Tsit5(); callback = ps_cb)
 plot(sol)
 ```
 The result looks as expected. However, what happens if we attempt to run the
 simulation again?
 ```@example ex2
-sol = solve(oprob; callback = ps_cb)
+sol = solve(oprob, Tsit5(); callback = ps_cb)
 plot(sol)
 ```
 The plot looks different, even though we simulate the same problem. Furthermore,
@@ -226,15 +226,15 @@ p = [:k => 1.0]
 oprob = ODEProblem(rn, u0, tspan, p)
 
 condition = [5.0]
-affect!(integrator) = integrator[:k] = 5.0
+affect!(integrator) = integrator.ps[:k] = 5.0
 ps_cb = PresetTimeCallback(condition, affect!)
 
-sol = solve(deepcopy(oprob); callback = ps_cb)
+sol = solve(deepcopy(oprob), Tsit5(); callback = ps_cb)
 plot(sol)
 ```
 where we parse a copy of our `ODEProblem` to the solver (using `deepcopy`). We can now run
 ```@example ex2
-sol = solve(deepcopy(oprob); callback = ps_cb)
+sol = solve(deepcopy(oprob), Tsit5(); callback = ps_cb)
 plot(sol)
 ```
 and get the expected result.
@@ -245,15 +245,15 @@ has to bundle them together in a `CallbackSet`, here follows one example:
 rn = @reaction_network begin
     (k,1), X1 <--> X2
 end
-u0 = [:X1 => 10.0,:X2 => 0.0]
+u0 = [:X1 => 10.0, :X2 => 0.0]
 tspan = (0.0, 20.0)
 p = [:k => 1.0]
 oprob = ODEProblem(rn, u0, tspan, p)
 
 ps_cb_1 = PresetTimeCallback([3.0, 7.0], integ -> integ[:X1] += 5.0)
-ps_cb_2 = PresetTimeCallback([5.0], integ -> integ[:k] = 5.0)
+ps_cb_2 = PresetTimeCallback([5.0], integ -> integ.ps[:k] = 5.0)
 
-sol = solve(deepcopy(oprob); callback=CallbackSet(ps_cb_1, ps_cb_2))
+sol = solve(deepcopy(oprob), Tsit5(); callback=CallbackSet(ps_cb_1, ps_cb_2))
 plot(sol)
 ```
 
@@ -279,7 +279,7 @@ jprob = JumpProblem(rn, dprob, Direct())
 condition = [5.0]
 function affect!(integrator)
     integrator[:X1] += 5.0
-    integrator[:k] += 2.0
+    integrator.ps[:k] += 2.0
     reset_aggregated_jumps!(integrator)
     nothing
 end
@@ -350,12 +350,12 @@ It is possible to use a different noise scaling expression for each reaction. He
 ```@example ex3
 rn_4 = @reaction_network begin
     (p, d), 0 <--> X
-    (p, d), 0 <--> Y, ([noise_scaling=0.0, noise_scaling=0.0])
+    (p, d), 0 <--> Y, ([noise_scaling=0.0], [noise_scaling=0.0])
 end
 
 u0_4 = [:X => 10.0, :Y => 10.0]
 tspan = (0.0, 10.0)
-p_4 = [p => 10.0, d => 1.]
+p_4 = [:p => 10.0, :d => 1.]
 
 sprob_4 = SDEProblem(rn_4, u0_4, tspan, p_4)
 sol_4 = solve(sprob_4, ImplicitEM())
diff --git a/docs/src/model_simulation/simulation_structure_interfacing.md b/docs/src/model_simulation/simulation_structure_interfacing.md
index 6675d1ef1b..672b14b3a6 100644
--- a/docs/src/model_simulation/simulation_structure_interfacing.md
+++ b/docs/src/model_simulation/simulation_structure_interfacing.md
@@ -3,6 +3,9 @@ When simulating a model, one begins with creating a [problem](https://docs.sciml
 
 Generally, when we have a structure `simulation_struct` and want to interface with the unknown (or parameter) `G`, we use `simulation_struct[:G]` to access the value, and `simulation_struct[:G] = 5.0` to set it to a new value. However, see the following examples for full details.
 
+!!! note
+    The following tutorial will describe how to interface with problems, integrators, and solutions using `[]` notation. An alternative is to use the `ModelingToolkit.getu`, `ModelingToolkit.getp`, `ModelingToolkit.setu`, and `ModelingToolkit.setp` functions. These requires an additional step to use, however, they can also improve performance when a very large number interfaces are carried out.
+
 ## Interfacing problem objects
 
 We begin by demonstrating how we can interface with problem objects. We will demonstrate using a `ODEProblem`, however, it works similarly for other problem types.
@@ -24,7 +27,7 @@ oprob[:X1]
 ```
 with the notation being identical for parameters:
 ```@example ex1
-oprob[:k1]
+oprob.ps[:k1]
 ```
 
 If we want to change an unknown's initial condition value, we use the following notation
@@ -33,6 +36,9 @@ oprob[:X1] = 10.0
 ```
 with parameters using the same notation.
 
+!!! note
+    When interfacing with a parameter, `.ps` must be appended to the structure uses (e.g. `oprob`). This is not done when species are interfaced with.
+
 #### [Remaking problems using the `remake` function](@id simulation_structure_interfacing_remake)
 Typically, when modifying problems, it is recommended to use the `remake` function. Unlike when we do `oprob[:X1] = 10.0` (which modifies the problem in question), `remake` creates a new problem object. The `remake` function takes a problem as input, and any fields you wish to modify (and their new values) as optional inputs. Thus, we can do:
 ```@example ex1
@@ -74,7 +80,7 @@ integrator[:X1]
 ```
 or a parameter:
 ```@example ex1
-integrator[:k1]
+integrator.ps[:k1]
 ```
 Similarly, we can update their values using:
 ```@example ex1
@@ -95,7 +101,7 @@ sol[:X1]
 ```
 while when we access a parameter we only get a single value:
 ```@example ex1
-sol[:k1]
+sol.ps[:k1]
 ```
 Finally, we note that we cannot change the values of solution unknowns or parameters (i.e. both `sol[:X1] = 0.0` and `sol[:k1] = 0.0` generate errors).
 
@@ -103,6 +109,7 @@ Finally, we note that we cannot change the values of solution unknowns or parame
 
 Catalyst is built on an *intermediary representation* implemented by (ModelingToolkit.jl)[https://github.com/SciML/ModelingToolkit.jl]. ModelingToolkit is a modelling framework where one first declares a set of symbolic variables and parameters using e.g.
 ```@example ex2
+using Catalyst # hide
 using ModelingToolkit
 t = default_t()
 @parameters σ ρ β
@@ -124,7 +131,7 @@ u0 = [X1 => 1.0, X2 => 5.0]
 p = [:k1 => 5.0, :k2 => 2.0]
 oprob = ODEProblem(rn, u0, (0.0,10.0), p)
 
-oprob[k1]
+oprob.ps[k1]
 ```
 To interface with integrators and solutions we use a similar syntax.
 
diff --git a/docs/src/steady_state_functionality/nonlinear_solve.md b/docs/src/steady_state_functionality/nonlinear_solve.md
index 0d4a903714..acf20d56c9 100644
--- a/docs/src/steady_state_functionality/nonlinear_solve.md
+++ b/docs/src/steady_state_functionality/nonlinear_solve.md
@@ -66,7 +66,7 @@ Like previously, the found steady state is unaffected by the initial `u` guess.
 
 
 ## [Systems with conservation laws](@id nonlinear_solve_conservation_laws)
-As described in the section on homotopy continuation, when finding the steady states of systems with conservation laws, [additional considerations have to be taken](homotopy_continuation_conservation_laws). E.g. consider the following two-state system:
+As described in the section on homotopy continuation, when finding the steady states of systems with conservation laws, [additional considerations have to be taken](@ref homotopy_continuation_conservation_laws). E.g. consider the following two-state system:
 ```@example nonlinear_solve2
 using Catalyst, NonlinearSolve # hide
 two_state_model = @reaction_network begin