From 979a166af5bc3434b571907c47fb6f32b9fce96d Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Tue, 16 Jan 2024 10:34:57 -0500
Subject: [PATCH 01/17] init

---
 docs/Project.toml                             |   2 +
 .../dynamical_systems.md                      | 149 ++++++++++++++++++
 .../homotopy_continuation.md                  |   2 +-
 3 files changed, 152 insertions(+), 1 deletion(-)
 create mode 100644 docs/src/catalyst_applications/dynamical_systems.md

diff --git a/docs/Project.toml b/docs/Project.toml
index e7454091d2..372921dbde 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,11 +1,13 @@
 [deps]
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 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"
+DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
new file mode 100644
index 0000000000..e118e39893
--- /dev/null
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -0,0 +1,149 @@
+# [Analysing model steady state properties with DynamicalSystems.jl](@id dynamical_systems)
+The [DynamicalSystems.jl package](https://github.com/JuliaDynamics/DynamicalSystems.jl) implements a wide range of methods for analysing dynamical systems. This includes both continuous-time systems (i.e. ODEs) and discrete-times ones (difference equations, these are, however, not relevant to chemical reaction network modelling). For ODEs, methods primarily include those for analysing a system's steady state and stability properties. Here we give two examples of how DynamicalSystems.jl can be used, with the package's [documentation describing many more features](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/tutorial/). Finally, it should also be noted that DynamicalSystems.jl contain several tools for [analysing data measured from dynamical systems](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/contents/#Exported-submodules).
+
+## [Finding basins of attraction](@id dynamical_systems_basins_of_attraction)
+Given enough time, an ODE will eventually reach a so-called [*attractor*](https://en.wikipedia.org/wiki/Attractor). For chemical reaction networks (CRNs), this will (almost always) be either a *steady state* or a *limit cycle*. Since ODEs are deterministic, which attractor a simulation will reach is uniquely determined by the initial condition (assuming parameter values are fixed). Conversely, each attractor is associated with a set of initial conditions such that model simulations originating in these will tend to that attractor. These sets are called *basins of attraction*. Here, phase space (the space of all possible states of the system) can be divided into a number of basins of attraction equal to the number of attractors. 
+
+DynamicalSystems.jl provides a simple interface for finding an ODE's basins of attraction across any given subspace of phase space. In this example we will use the bistable [Wilhelm model](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) (which steady states we have previous [computed using homotopy continuation](@ref homotopy_continuation_basic_example)). As a first step, we create a `ODEProblem` corresponding to the model which basins of attraction we wish to compute. For this application, `u0` is unused, and its value is irrelevant. The value of `tspan` is not very important, however, it must provide long enough a time range (else basins will not be successfully found). 
+```@example dynamical_systems_basins
+using Catalyst
+wilhelm_model = @reaction_network begin
+    k1, Y --> 2X
+    k2, 2X --> X + Y
+    k3, X + Y --> Y
+    k4, X --> 0
+end
+
+u0 = [:X => 0.0, :Y => .0]
+tspan = (0.0, 10.0)
+ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5]
+oprob = ODEProblem(wilhelm_model, u0, tspan, ps)
+nothing # hide
+```
+Next, for any application of DynamicalSystems.jl, our `ODEProblem` must be converted into a so-called `CoupledODEs` structure. This is done by combining the ODE with the solver (and potential solver options) with which we wish to simulate it (just like when it is simulated using `solve`). Here, we will simply designate the `Tsit5` numeric solver (but provide no other options).
+```@example dynamical_systems_basins
+using DynamicalSystems, OrdinaryDiffEq
+ds = CoupledODEs(oprob, (alg = Tsit5(),))
+nothing # hide
+```
+We can now compute the basins of attraction. This is done by first designating a grid designate which subspace of phase-space we wish to investigate (here, the corresponding basin of attraction is found for every point on the grid). Next, we first create a `AttractorsViaRecurrences` struct and then use that as input to the `basins_of_attraction` function.
+```@example dynamical_systems_basins
+# We provide one grid of values for each species. These are then bundled into a tuple.
+x_grid = 0.0:0.01:5.0
+y_grid = 0.0:0.01:10.0
+grid = (x_grid, y_grid)
+avr = AttractorsViaRecurrences(ds, grid)
+basins, attractors = basins_of_attraction(avr, grid)
+nothing # hide
+```
+Here, `attractors` is a vector containing the attractors (in this case two fixed points, one at $(0.0,0.0)$ and one at $(4.5,6.0)$). Next, `basins` is a matrix of equal size to `grid`, where each value is an integer describing to which attractor's basin that state belong. 
+
+DynamicalSystems.jl also provide a simple interface for plotting the resulting basins. This uses [Makie.jl](https://docs.makie.org/stable/), and alternative plotting package to [Plots.jl](https://github.com/JuliaPlots/Plots.jl) (which is typically the preferred plotting package within the context of Catalyst). Generally, Makie is good at creating animations or interactive graphics (however, it is also a popular competitor to Plots.jl for normal plotting). 
+```@example dynamical_systems_basins
+using CairoMakie
+heatmap_basins_attractors(grid, basins, attractors)
+```
+Here, in addition to the basins of attraction, the system's three steady states are marked (the one at the intersection of the two basins is unstable).
+
+!!! warning
+    Both Makie and Plots.jl exports a function called `plot`. Hence, if both these packages are imported into the same session, calls to `plot` must be prepended with the package one wishes to use (e.g. `Plots.plot(sol)`).
+
+More information on how to compute basins of attractions for ODEs using DynamicalSystems.jl can be found [here](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/attractors/stable/basins/).
+
+## [Computing Lyapunov exponents](@id dynamical_systems_lyapunov_exponents)
+[*Lyapunov exponents*](https://en.wikipedia.org/wiki/Lyapunov_exponent) are scalar values that can be computed for any state of an ODE. For an ODE with $n$ variables, for each state, a total of $n$ Lyapunov exponents can be computed (and these are collectively called the *Lyapunov spectrum*). Large Lyapunov exponents indicates that trajectories infinitesimally close to the given state diverge from each other. Conversely, small Lyapunov exponents suggests that trajectories converge to each others. 
+
+While Lyapunov exponents can be used for other purposes, they are primarily used to characterise [*chaotic behaviours*](https://en.wikipedia.org/wiki/Chaos_theory) (where small changes in initial conditions has large effect on the resulting trajectories). Generally, an ODE exhibit chaotic behaviours in a point in phase space if that point have *at least one* positive Lyapunov exponent. Practically, Lyapunov exponents can be computed using DynamicalSystems.jl's `lyapunovspectrum` function. Here we will use it to investigate two models, one which exhibits chaos and one which do not.
+
+First, let us consider the [Willamowski–Rössler](https://www.degruyter.com/document/doi/10.1515/zna-1980-0308/html?lang=en) model, which is known to exhibit chaotic behaviours. 
+```@example dynamical_systems_lyapunov
+using Catalyst
+wr_model = @reaction_network begin
+    k1, 2X --> 3X
+    k2, X --> 2X
+    k3, Z + 2X --> 2Z
+    k4, Y + X --> 2Y
+    k5, Y --> ∅
+    k6, 2Z --> ∅
+    k7, Z --> ∅
+end
+```
+We can simulate the model, noting that its behaviour seems chaotic.
+```@example dynamical_systems_lyapunov
+using OrdinaryDiffEq, Plots
+
+u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5]
+tspan = (0.0, 100.0)
+p = [:k1 => 2.1, :k2 => 0.7, :k3 => 2.9, :k4 => 1.1, :k5 => 1.0, :k6 => 0.5, :k7 => 2.7]
+
+oprob = ODEProblem(wr_model, u0, tspan, p)
+sol = solve(oprob, Rodas5P())
+plot(sol; idxs=(:X, :Y, :Z))
+```
+Next, like when we [computed basins of attraction](@ref dynamical_systems_basins_of_attraction), we create a `CoupledODEs` corresponding to the model and state for which we wish to compute our Lyapunov spectrum (here, the `tspan` of the used `ODEProblem` have no impact on the computed Lyapunov spectrum).
+```@example dynamical_systems_lyapunov
+using DynamicalSystems
+ds = CoupledODEs(oprob, (alg = Rodas5P(autodiff=false),))
+nothing # hide
+```
+Here, the `autodiff=false` argument is required when Lyapunov spectrums are computed. We can now provide our `CoupledODEs` (`ds`) to `lyapunovspectrum` to compute the lyapunov spectrum. This function requires a second argument (here set to `100`). Generally setting this to a higher value will increase accuracy, but also increase runtime (since `lyapunovspectrum` is fast for most systems, setting this to a large value is recommended).
+```@example dynamical_systems_lyapunov
+lyapunovspectrum(ds, 100)
+```
+Here, the largest exponent is positive, suggesting that the model is chaotic in the point $(1.5, 1.5, 1.5)$.
+
+Next, we consider the [Brusselator] model. First we simulate the model for two similar initial conditions, confirming that they converge to the same limit cycle:
+```@example dynamical_systems_lyapunov
+brusselator = @reaction_network begin
+    A, ∅ --> X
+    1, 2X + Y --> 3X
+    B, X --> Y
+    1, X --> ∅
+end
+
+u0_1 = [:X => 1.0, :Y => 1.0]
+u0_2 = [:X => 1.2, :Y => 1.0]
+tspan = (0., 25.)
+ps = [:A => 1.0, :B => 4.0]
+
+oprob1 = ODEProblem(brusselator, u0_1, tspan, ps)
+oprob2 = ODEProblem(brusselator, u0_2, tspan, ps)
+osol1  = solve(oprob1, Tsit5())
+osol2  = solve(oprob2, Tsit5())
+plot(osol1; idxs = (:X, :Y))
+plot!(osol2; idxs = (:X, :Y))
+```
+Next, we compute the Lyapunov spectrum at one of the initial conditions:
+```@example dynamical_systems_lyapunov
+ds = CoupledODEs(oprob1, (alg = Rodas5P(autodiff=false),))
+lyapunovspectrum(ds, 100)
+```
+Here, all Lyapunov exponents are negative, confirming that the brusselator (at this point) is non-chaotic.
+
+More details on how to compute Lyapunov exponents using DynamicalSystems.jl can be found [here](https://juliadynamics.github.io/ChaosTools.jl/stable/lyapunovs/). A full overview of tools for analysing chaotic behaviours (using the "ChaosTools.jl subpackage) can be found [here](https://juliadynamics.github.io/ChaosTools.jl/stable/).
+
+
+---
+## [Citations](@id dynamical_systems_citations)
+If you use this functionality in your research, [in addition to Catalyst](@ref catalyst_citation), please cite the following paper to support the author of the DynamicalSystems.jl package:
+```
+@article{DynamicalSystems.jl-2018,
+  doi = {10.21105/joss.00598},
+  url = {https://doi.org/10.21105/joss.00598},
+  year  = {2018},
+  month = {mar},
+  volume = {3},
+  number = {23},
+  pages = {598},
+  author = {George Datseris},
+  title = {DynamicalSystems.jl: A Julia software library for chaos and nonlinear dynamics},
+  journal = {Journal of Open Source Software}
+}
+```
+
+
+---
+## References
+[^1]: [G. Datseris, U. Parlitz, *Nonlinear dynamics: A concise introduction interlaced with code*, Springer (2022).](https://link.springer.com/book/10.1007/978-3-030-91032-7)
+[^2]: [S. H. Strogatz, *Nonlinear Dynamics and Chaos*, Westview Press (1994).](http://users.uoa.gr/~pjioannou/nonlin/Strogatz,%20S.%20H.%20-%20Nonlinear%20Dynamics%20And%20Chaos.pdf)
+[^3]: [A. M. Lyapunov, *The general problem of the stability of motion*, International Journal of Control (1992).](https://www.tandfonline.com/doi/abs/10.1080/00207179208934253)
\ No newline at end of file
diff --git a/docs/src/steady_state_functionality/homotopy_continuation.md b/docs/src/steady_state_functionality/homotopy_continuation.md
index 1fb751bc15..b5fbe4f850 100644
--- a/docs/src/steady_state_functionality/homotopy_continuation.md
+++ b/docs/src/steady_state_functionality/homotopy_continuation.md
@@ -11,7 +11,7 @@ integer Hill exponents). The roots of these can reliably be found through a
 
 Catalyst contains a special homotopy continuation extension that is loaded whenever HomotopyContinuation.jl is. This exports a single function, `hc_steady_states`, that can be used to find the steady states of any `ReactionSystem` structure.
 
-## Basic example
+## [Basic example](@id homotopy_continuation_basic_example)
 For this tutorial, we will use a model from Wilhelm (2009)[^1] (which
 demonstrates bistability in a small chemical reaction network). We declare the
 model and the parameter set for which we want to find the steady states:

From 2fb5f3215e8045594a4906fd5249f6a19b646217 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Tue, 16 Jan 2024 15:24:20 -0500
Subject: [PATCH 02/17] remove printout

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

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index e118e39893..71c0c77173 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -33,7 +33,7 @@ x_grid = 0.0:0.01:5.0
 y_grid = 0.0:0.01:10.0
 grid = (x_grid, y_grid)
 avr = AttractorsViaRecurrences(ds, grid)
-basins, attractors = basins_of_attraction(avr, grid)
+basins, attractors = basins_of_attraction(avr, grid; show_progress = false)
 nothing # hide
 ```
 Here, `attractors` is a vector containing the attractors (in this case two fixed points, one at $(0.0,0.0)$ and one at $(4.5,6.0)$). Next, `basins` is a matrix of equal size to `grid`, where each value is an integer describing to which attractor's basin that state belong. 

From 4a15e96b56fc4a3bdfd2f3f044c7876624c0d425 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 09:11:36 -0500
Subject: [PATCH 03/17] Update
 docs/src/catalyst_applications/dynamical_systems.md

Co-authored-by: George Datseris <datseris.george@gmail.com>
---
 docs/src/catalyst_applications/dynamical_systems.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index 71c0c77173..8b7ccbb775 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -1,5 +1,5 @@
 # [Analysing model steady state properties with DynamicalSystems.jl](@id dynamical_systems)
-The [DynamicalSystems.jl package](https://github.com/JuliaDynamics/DynamicalSystems.jl) implements a wide range of methods for analysing dynamical systems. This includes both continuous-time systems (i.e. ODEs) and discrete-times ones (difference equations, these are, however, not relevant to chemical reaction network modelling). For ODEs, methods primarily include those for analysing a system's steady state and stability properties. Here we give two examples of how DynamicalSystems.jl can be used, with the package's [documentation describing many more features](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/tutorial/). Finally, it should also be noted that DynamicalSystems.jl contain several tools for [analysing data measured from dynamical systems](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/contents/#Exported-submodules).
+The [DynamicalSystems.jl package](https://github.com/JuliaDynamics/DynamicalSystems.jl) implements a wide range of methods for analysing dynamical systems. This includes both continuous-time systems (i.e. ODEs) and discrete-times ones (difference equations, these are, however, not relevant to chemical reaction network modelling). Here we give two examples of how DynamicalSystems.jl can be used, with the package's [documentation describing many more features](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/tutorial/). Finally, it should also be noted that DynamicalSystems.jl contain several tools for [analysing data measured from dynamical systems](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/contents/#Exported-submodules).
 
 ## [Finding basins of attraction](@id dynamical_systems_basins_of_attraction)
 Given enough time, an ODE will eventually reach a so-called [*attractor*](https://en.wikipedia.org/wiki/Attractor). For chemical reaction networks (CRNs), this will (almost always) be either a *steady state* or a *limit cycle*. Since ODEs are deterministic, which attractor a simulation will reach is uniquely determined by the initial condition (assuming parameter values are fixed). Conversely, each attractor is associated with a set of initial conditions such that model simulations originating in these will tend to that attractor. These sets are called *basins of attraction*. Here, phase space (the space of all possible states of the system) can be divided into a number of basins of attraction equal to the number of attractors. 

From 6beb7e1d469bfef5715afcc66d9eea828680cb5c Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 09:13:51 -0500
Subject: [PATCH 04/17] Update
 docs/src/catalyst_applications/dynamical_systems.md

Co-authored-by: George Datseris <datseris.george@gmail.com>
---
 docs/src/catalyst_applications/dynamical_systems.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index 8b7ccbb775..53dddda6f2 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -26,7 +26,7 @@ using DynamicalSystems, OrdinaryDiffEq
 ds = CoupledODEs(oprob, (alg = Tsit5(),))
 nothing # hide
 ```
-We can now compute the basins of attraction. This is done by first designating a grid designate which subspace of phase-space we wish to investigate (here, the corresponding basin of attraction is found for every point on the grid). Next, we first create a `AttractorsViaRecurrences` struct and then use that as input to the `basins_of_attraction` function.
+We can now compute the basins of attraction. This is done by first creating a grid that designates which subspace of phase-space we wish to investigate (here, the corresponding basin of attraction is found for every point on the grid). Next, we create a `AttractorsViaRecurrences` struct, that maps initial conditions to attractors, and then use that as input to the `basins_of_attraction` function.
 ```@example dynamical_systems_basins
 # We provide one grid of values for each species. These are then bundled into a tuple.
 x_grid = 0.0:0.01:5.0

From 0cccddb2301ee74ef71777923feda4d3665c9f6f Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 09:14:51 -0500
Subject: [PATCH 05/17] Update
 docs/src/catalyst_applications/dynamical_systems.md

Co-authored-by: George Datseris <datseris.george@gmail.com>
---
 docs/src/catalyst_applications/dynamical_systems.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index 53dddda6f2..4193239db4 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -36,7 +36,7 @@ avr = AttractorsViaRecurrences(ds, grid)
 basins, attractors = basins_of_attraction(avr, grid; show_progress = false)
 nothing # hide
 ```
-Here, `attractors` is a vector containing the attractors (in this case two fixed points, one at $(0.0,0.0)$ and one at $(4.5,6.0)$). Next, `basins` is a matrix of equal size to `grid`, where each value is an integer describing to which attractor's basin that state belong. 
+Here, `attractors` is a dictionary that maps attractor labels (the integers) to attractors. In this case we have two fixed points, one at $(0.0,0.0)$ and one at $(4.5,6.0)$. Next, `basins` is a matrix of equal size to `grid`, where each value is an integer describing to which attractor's basin that state belongs.
 
 DynamicalSystems.jl also provide a simple interface for plotting the resulting basins. This uses [Makie.jl](https://docs.makie.org/stable/), and alternative plotting package to [Plots.jl](https://github.com/JuliaPlots/Plots.jl) (which is typically the preferred plotting package within the context of Catalyst). Generally, Makie is good at creating animations or interactive graphics (however, it is also a popular competitor to Plots.jl for normal plotting). 
 ```@example dynamical_systems_basins

From eb32001a7a5a5e90d7ec3213afb33d44e07ed3f2 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 09:23:23 -0500
Subject: [PATCH 06/17] Update
 docs/src/catalyst_applications/dynamical_systems.md

Co-authored-by: George Datseris <datseris.george@gmail.com>
---
 docs/src/catalyst_applications/dynamical_systems.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index 4193239db4..bfcde97ee1 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -90,7 +90,7 @@ Here, the `autodiff=false` argument is required when Lyapunov spectrums are comp
 ```@example dynamical_systems_lyapunov
 lyapunovspectrum(ds, 100)
 ```
-Here, the largest exponent is positive, suggesting that the model is chaotic in the point $(1.5, 1.5, 1.5)$.
+Here, the largest exponent is positive, suggesting that the model is chaotic (or, more accurately, it has at least one  chaotic attractor, to which we go to from the initial condition $(1.5, 1.5, 1.5)).
 
 Next, we consider the [Brusselator] model. First we simulate the model for two similar initial conditions, confirming that they converge to the same limit cycle:
 ```@example dynamical_systems_lyapunov

From 3c60e91cb826142ed2350f59c96d625becd6ab5e Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 09:23:36 -0500
Subject: [PATCH 07/17] Update
 docs/src/catalyst_applications/dynamical_systems.md

Co-authored-by: George Datseris <datseris.george@gmail.com>
---
 docs/src/catalyst_applications/dynamical_systems.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index bfcde97ee1..0e7fb53f45 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -118,7 +118,7 @@ Next, we compute the Lyapunov spectrum at one of the initial conditions:
 ds = CoupledODEs(oprob1, (alg = Rodas5P(autodiff=false),))
 lyapunovspectrum(ds, 100)
 ```
-Here, all Lyapunov exponents are negative, confirming that the brusselator (at this point) is non-chaotic.
+Here, all Lyapunov exponents are negative, confirming that the brusselator is non-chaotic.
 
 More details on how to compute Lyapunov exponents using DynamicalSystems.jl can be found [here](https://juliadynamics.github.io/ChaosTools.jl/stable/lyapunovs/). A full overview of tools for analysing chaotic behaviours (using the "ChaosTools.jl subpackage) can be found [here](https://juliadynamics.github.io/ChaosTools.jl/stable/).
 

From 2b0f9a11f821170e8ce601c1d32b56091bab87d1 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 09:24:32 -0500
Subject: [PATCH 08/17] Update
 docs/src/catalyst_applications/dynamical_systems.md

Co-authored-by: George Datseris <datseris.george@gmail.com>
---
 docs/src/catalyst_applications/dynamical_systems.md | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index 0e7fb53f45..e990f12a48 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -140,7 +140,9 @@ If you use this functionality in your research, [in addition to Catalyst](@ref c
   journal = {Journal of Open Source Software}
 }
 ```
+## Learning more
 
+If you want to learn more about analyzing dynamical systems, including chaotic behavior, you can have a look at the textbook [Nonlinear Dynamics](https://link.springer.com/book/10.1007/978-3-030-91032-7). It utilizes DynamicalSystems.jl and provides a concise, hands-on approach to learning nonlinear dynamics and analyzing dynamical systems.
 
 ---
 ## References

From a821b382dbdff901cf20756004efcd1ea4066d05 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 09:26:55 -0500
Subject: [PATCH 09/17] Update
 docs/src/catalyst_applications/dynamical_systems.md

Co-authored-by: George Datseris <datseris.george@gmail.com>
---
 docs/src/catalyst_applications/dynamical_systems.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index e990f12a48..6cbff9fa6d 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -51,7 +51,7 @@ Here, in addition to the basins of attraction, the system's three steady states
 More information on how to compute basins of attractions for ODEs using DynamicalSystems.jl can be found [here](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/attractors/stable/basins/).
 
 ## [Computing Lyapunov exponents](@id dynamical_systems_lyapunov_exponents)
-[*Lyapunov exponents*](https://en.wikipedia.org/wiki/Lyapunov_exponent) are scalar values that can be computed for any state of an ODE. For an ODE with $n$ variables, for each state, a total of $n$ Lyapunov exponents can be computed (and these are collectively called the *Lyapunov spectrum*). Large Lyapunov exponents indicates that trajectories infinitesimally close to the given state diverge from each other. Conversely, small Lyapunov exponents suggests that trajectories converge to each others. 
+[*Lyapunov exponents*](https://en.wikipedia.org/wiki/Lyapunov_exponent) are scalar values that can be computed for any attractor of an ODE. For an ODE with $n$ variables, for each state, a total of $n$ Lyapunov exponents can be computed (and these are collectively called the *Lyapunov spectrum*). Positive Lyapunov exponents indicate that trajectories initially infinitesimally close diverge from each other. Conversely, negative Lyapunov exponents suggests that trajectories converge to each others. 
 
 While Lyapunov exponents can be used for other purposes, they are primarily used to characterise [*chaotic behaviours*](https://en.wikipedia.org/wiki/Chaos_theory) (where small changes in initial conditions has large effect on the resulting trajectories). Generally, an ODE exhibit chaotic behaviours in a point in phase space if that point have *at least one* positive Lyapunov exponent. Practically, Lyapunov exponents can be computed using DynamicalSystems.jl's `lyapunovspectrum` function. Here we will use it to investigate two models, one which exhibits chaos and one which do not.
 

From f6ac0ef42aae857f07fdb0fed878240dbec5ebc9 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 09:28:44 -0500
Subject: [PATCH 10/17] Update
 docs/src/catalyst_applications/dynamical_systems.md

Co-authored-by: George Datseris <datseris.george@gmail.com>
---
 docs/src/catalyst_applications/dynamical_systems.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index 6cbff9fa6d..189e5c49a0 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -53,7 +53,7 @@ More information on how to compute basins of attractions for ODEs using Dynamica
 ## [Computing Lyapunov exponents](@id dynamical_systems_lyapunov_exponents)
 [*Lyapunov exponents*](https://en.wikipedia.org/wiki/Lyapunov_exponent) are scalar values that can be computed for any attractor of an ODE. For an ODE with $n$ variables, for each state, a total of $n$ Lyapunov exponents can be computed (and these are collectively called the *Lyapunov spectrum*). Positive Lyapunov exponents indicate that trajectories initially infinitesimally close diverge from each other. Conversely, negative Lyapunov exponents suggests that trajectories converge to each others. 
 
-While Lyapunov exponents can be used for other purposes, they are primarily used to characterise [*chaotic behaviours*](https://en.wikipedia.org/wiki/Chaos_theory) (where small changes in initial conditions has large effect on the resulting trajectories). Generally, an ODE exhibit chaotic behaviours in a point in phase space if that point have *at least one* positive Lyapunov exponent. Practically, Lyapunov exponents can be computed using DynamicalSystems.jl's `lyapunovspectrum` function. Here we will use it to investigate two models, one which exhibits chaos and one which do not.
+While Lyapunov exponents can be used for other purposes, they are primarily used to characterise [*chaotic behaviours*](https://en.wikipedia.org/wiki/Chaos_theory) (where small changes in initial conditions has large effect on the resulting trajectories). Generally, an ODE exhibit chaotic behaviour if its attractor(s) have *at least one* positive Lyapunov exponent. Practically, Lyapunov exponents can be computed using DynamicalSystems.jl's `lyapunovspectrum` function. Here we will use it to investigate two models, one which exhibits chaos and one which do not.
 
 First, let us consider the [Willamowski–Rössler](https://www.degruyter.com/document/doi/10.1515/zna-1980-0308/html?lang=en) model, which is known to exhibit chaotic behaviours. 
 ```@example dynamical_systems_lyapunov

From aa523efb8ce31a3f365388f4e36f592123b0f6cd Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 09:29:06 -0500
Subject: [PATCH 11/17] Update
 docs/src/catalyst_applications/dynamical_systems.md

Co-authored-by: George Datseris <datseris.george@gmail.com>
---
 docs/src/catalyst_applications/dynamical_systems.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index 189e5c49a0..3ae6edeeef 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -55,7 +55,7 @@ More information on how to compute basins of attractions for ODEs using Dynamica
 
 While Lyapunov exponents can be used for other purposes, they are primarily used to characterise [*chaotic behaviours*](https://en.wikipedia.org/wiki/Chaos_theory) (where small changes in initial conditions has large effect on the resulting trajectories). Generally, an ODE exhibit chaotic behaviour if its attractor(s) have *at least one* positive Lyapunov exponent. Practically, Lyapunov exponents can be computed using DynamicalSystems.jl's `lyapunovspectrum` function. Here we will use it to investigate two models, one which exhibits chaos and one which do not.
 
-First, let us consider the [Willamowski–Rössler](https://www.degruyter.com/document/doi/10.1515/zna-1980-0308/html?lang=en) model, which is known to exhibit chaotic behaviours. 
+First, let us consider the [Willamowski–Rössler](https://www.degruyter.com/document/doi/10.1515/zna-1980-0308/html?lang=en) model, which is known to exhibit chaotic behaviour.
 ```@example dynamical_systems_lyapunov
 using Catalyst
 wr_model = @reaction_network begin

From d28250ff5f7bae94d39b05606d4b61d947af2da7 Mon Sep 17 00:00:00 2001
From: Torkel Loman <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 09:36:52 -0500
Subject: [PATCH 12/17] Update
 docs/src/catalyst_applications/dynamical_systems.md

Co-authored-by: George Datseris <datseris.george@gmail.com>
---
 docs/src/catalyst_applications/dynamical_systems.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index 3ae6edeeef..d1b7b9a32b 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -34,7 +34,7 @@ y_grid = 0.0:0.01:10.0
 grid = (x_grid, y_grid)
 avr = AttractorsViaRecurrences(ds, grid)
 basins, attractors = basins_of_attraction(avr, grid; show_progress = false)
-nothing # hide
+attractors
 ```
 Here, `attractors` is a dictionary that maps attractor labels (the integers) to attractors. In this case we have two fixed points, one at $(0.0,0.0)$ and one at $(4.5,6.0)$. Next, `basins` is a matrix of equal size to `grid`, where each value is an integer describing to which attractor's basin that state belongs.
 

From f0894531d765521cde52d0e527d26651187871aa Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 09:37:10 -0500
Subject: [PATCH 13/17] up

---
 docs/src/catalyst_applications/dynamical_systems.md | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index d1b7b9a32b..b11c949bad 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -24,13 +24,12 @@ Next, for any application of DynamicalSystems.jl, our `ODEProblem` must be conve
 ```@example dynamical_systems_basins
 using DynamicalSystems, OrdinaryDiffEq
 ds = CoupledODEs(oprob, (alg = Tsit5(),))
-nothing # hide
 ```
 We can now compute the basins of attraction. This is done by first creating a grid that designates which subspace of phase-space we wish to investigate (here, the corresponding basin of attraction is found for every point on the grid). Next, we create a `AttractorsViaRecurrences` struct, that maps initial conditions to attractors, and then use that as input to the `basins_of_attraction` function.
 ```@example dynamical_systems_basins
 # We provide one grid of values for each species. These are then bundled into a tuple.
-x_grid = 0.0:0.01:5.0
-y_grid = 0.0:0.01:10.0
+x_grid = 0.0:0.03:6.0
+y_grid = 0.0:0.03:9.0
 grid = (x_grid, y_grid)
 avr = AttractorsViaRecurrences(ds, grid)
 basins, attractors = basins_of_attraction(avr, grid; show_progress = false)
@@ -38,7 +37,7 @@ attractors
 ```
 Here, `attractors` is a dictionary that maps attractor labels (the integers) to attractors. In this case we have two fixed points, one at $(0.0,0.0)$ and one at $(4.5,6.0)$. Next, `basins` is a matrix of equal size to `grid`, where each value is an integer describing to which attractor's basin that state belongs.
 
-DynamicalSystems.jl also provide a simple interface for plotting the resulting basins. This uses [Makie.jl](https://docs.makie.org/stable/), and alternative plotting package to [Plots.jl](https://github.com/JuliaPlots/Plots.jl) (which is typically the preferred plotting package within the context of Catalyst). Generally, Makie is good at creating animations or interactive graphics (however, it is also a popular competitor to Plots.jl for normal plotting). 
+DynamicalSystems.jl also provides a simple interface for plotting the resulting basins. This uses [Makie.jl](https://docs.makie.org/stable/), an alternative plotting package to [Plots.jl](https://github.com/JuliaPlots/Plots.jl) (which is typically the preferred plotting package within the context of Catalyst). Generally, Makie is good at creating animations or interactive graphics (however, it is also a [popular competitor to Plots.jl for general-purpose plotting](https://juliapackagecomparisons.github.io/pages/plotting/)). 
 ```@example dynamical_systems_basins
 using CairoMakie
 heatmap_basins_attractors(grid, basins, attractors)

From cde0d45f7afe1d1fe7ced37006a7b1a776ae9229 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 14:34:38 -0500
Subject: [PATCH 14/17] up

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

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index b11c949bad..a870bf932e 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -2,9 +2,9 @@
 The [DynamicalSystems.jl package](https://github.com/JuliaDynamics/DynamicalSystems.jl) implements a wide range of methods for analysing dynamical systems. This includes both continuous-time systems (i.e. ODEs) and discrete-times ones (difference equations, these are, however, not relevant to chemical reaction network modelling). Here we give two examples of how DynamicalSystems.jl can be used, with the package's [documentation describing many more features](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/tutorial/). Finally, it should also be noted that DynamicalSystems.jl contain several tools for [analysing data measured from dynamical systems](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/contents/#Exported-submodules).
 
 ## [Finding basins of attraction](@id dynamical_systems_basins_of_attraction)
-Given enough time, an ODE will eventually reach a so-called [*attractor*](https://en.wikipedia.org/wiki/Attractor). For chemical reaction networks (CRNs), this will (almost always) be either a *steady state* or a *limit cycle*. Since ODEs are deterministic, which attractor a simulation will reach is uniquely determined by the initial condition (assuming parameter values are fixed). Conversely, each attractor is associated with a set of initial conditions such that model simulations originating in these will tend to that attractor. These sets are called *basins of attraction*. Here, phase space (the space of all possible states of the system) can be divided into a number of basins of attraction equal to the number of attractors. 
+Given enough time, an ODE will eventually reach a so-called [*attractor*](https://en.wikipedia.org/wiki/Attractor). For chemical reaction networks (CRNs), this will typically be either a *steady state* or a *limit cycle*. Since ODEs are deterministic, which attractor a simulation will reach is uniquely determined by the initial condition (assuming parameter values are fixed). Conversely, each attractor is associated with a set of initial conditions such that model simulations originating in these will tend to that attractor. These sets are called *basins of attraction*. Here, phase space (the space of all possible states of the system) can be divided into a number of basins of attraction equal to the number of attractors. 
 
-DynamicalSystems.jl provides a simple interface for finding an ODE's basins of attraction across any given subspace of phase space. In this example we will use the bistable [Wilhelm model](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) (which steady states we have previous [computed using homotopy continuation](@ref homotopy_continuation_basic_example)). As a first step, we create a `ODEProblem` corresponding to the model which basins of attraction we wish to compute. For this application, `u0` is unused, and its value is irrelevant. The value of `tspan` is not very important, however, it must provide long enough a time range (else basins will not be successfully found). 
+DynamicalSystems.jl provides a simple interface for finding an ODE's basins of attraction across any given subspace of phase space. In this example we will use the bistable [Wilhelm model](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) (which steady states we have previous [computed using homotopy continuation](@ref homotopy_continuation_basic_example)). As a first step, we create a `ODEProblem` corresponding to the model which basins of attraction we wish to compute. For this application, `u0` and `tspan` is unused, and their values are of little importance (the only exception is than `tspan`, for implementation reason, must provide a not to small interval, we recommend minimum `(0.0, 1.0)`). 
 ```@example dynamical_systems_basins
 using Catalyst
 wilhelm_model = @reaction_network begin
@@ -79,7 +79,7 @@ oprob = ODEProblem(wr_model, u0, tspan, p)
 sol = solve(oprob, Rodas5P())
 plot(sol; idxs=(:X, :Y, :Z))
 ```
-Next, like when we [computed basins of attraction](@ref dynamical_systems_basins_of_attraction), we create a `CoupledODEs` corresponding to the model and state for which we wish to compute our Lyapunov spectrum (here, the `tspan` of the used `ODEProblem` have no impact on the computed Lyapunov spectrum).
+Next, like when we [computed basins of attraction](@ref dynamical_systems_basins_of_attraction), we create a `CoupledODEs` corresponding to the model and state for which we wish to compute our Lyapunov spectrum. Lke previously, `tspan` must provide some small interval (at least `(0.0, 1.0)` is recommended), but else have no impact on the computed Lyapunov spectrum.
 ```@example dynamical_systems_lyapunov
 using DynamicalSystems
 ds = CoupledODEs(oprob, (alg = Rodas5P(autodiff=false),))
@@ -139,9 +139,13 @@ If you use this functionality in your research, [in addition to Catalyst](@ref c
   journal = {Journal of Open Source Software}
 }
 ```
+
+
+---
 ## Learning more
 
-If you want to learn more about analyzing dynamical systems, including chaotic behavior, you can have a look at the textbook [Nonlinear Dynamics](https://link.springer.com/book/10.1007/978-3-030-91032-7). It utilizes DynamicalSystems.jl and provides a concise, hands-on approach to learning nonlinear dynamics and analyzing dynamical systems.
+If you want to learn more about analysing dynamical systems, including chaotic behaviour, you can have a look at the textbook [Nonlinear Dynamics](https://link.springer.com/book/10.1007/978-3-030-91032-7). It utilizes DynamicalSystems.jl and provides a concise, hands-on approach to learning nonlinear dynamics and analysing dynamical systems [^1].
+
 
 ---
 ## References

From 12eb6c95aeb095d982365be97a0c3e6bb21889b3 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Fri, 2 Feb 2024 08:53:07 -0500
Subject: [PATCH 15/17] merge fixes

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

diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/catalyst_applications/dynamical_systems.md
index a870bf932e..41a89862c7 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/catalyst_applications/dynamical_systems.md
@@ -35,7 +35,7 @@ avr = AttractorsViaRecurrences(ds, grid)
 basins, attractors = basins_of_attraction(avr, grid; show_progress = false)
 attractors
 ```
-Here, `attractors` is a dictionary that maps attractor labels (the integers) to attractors. In this case we have two fixed points, one at $(0.0,0.0)$ and one at $(4.5,6.0)$. Next, `basins` is a matrix of equal size to `grid`, where each value is an integer describing to which attractor's basin that state belongs.
+Here, `attractors` is a dictionary that maps attractor labels (integers) to attractors. In this case we have two fixed points, one at $(0.0,0.0)$ and one at $(4.5,6.0)$. Next, `basins` is a matrix of equal size to `grid`, where each value is an integer describing to which attractor's basin that state belongs.
 
 DynamicalSystems.jl also provides a simple interface for plotting the resulting basins. This uses [Makie.jl](https://docs.makie.org/stable/), an alternative plotting package to [Plots.jl](https://github.com/JuliaPlots/Plots.jl) (which is typically the preferred plotting package within the context of Catalyst). Generally, Makie is good at creating animations or interactive graphics (however, it is also a [popular competitor to Plots.jl for general-purpose plotting](https://juliapackagecomparisons.github.io/pages/plotting/)). 
 ```@example dynamical_systems_basins
@@ -89,7 +89,7 @@ Here, the `autodiff=false` argument is required when Lyapunov spectrums are comp
 ```@example dynamical_systems_lyapunov
 lyapunovspectrum(ds, 100)
 ```
-Here, the largest exponent is positive, suggesting that the model is chaotic (or, more accurately, it has at least one  chaotic attractor, to which we go to from the initial condition $(1.5, 1.5, 1.5)).
+Here, the largest exponent is positive, suggesting that the model is chaotic (or, more accurately, it has at least one chaotic attractor, to which it go to from the initial condition $(1.5, 1.5, 1.5)).
 
 Next, we consider the [Brusselator] model. First we simulate the model for two similar initial conditions, confirming that they converge to the same limit cycle:
 ```@example dynamical_systems_lyapunov

From 90ff79607a35197e93a265717874601f4eb58e37 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Thu, 16 May 2024 15:31:10 -0400
Subject: [PATCH 16/17] rebase

---
 docs/pages.jl                                      |  2 +-
 .../dynamical_systems.md                           | 14 +++++++-------
 2 files changed, 8 insertions(+), 8 deletions(-)
 rename docs/src/{catalyst_applications => steady_state_functionality}/dynamical_systems.md (84%)

diff --git a/docs/pages.jl b/docs/pages.jl
index 0f6b95bfac..b54e11d80e 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -39,7 +39,7 @@ pages = Any[
         "steady_state_functionality/nonlinear_solve.md",
         "steady_state_functionality/steady_state_stability_computation.md",        
         "steady_state_functionality/bifurcation_diagrams.md"
-        # Dynamical systems analysis.
+        "steady_state_functionality/dynamical_systems.md"
     ],
     "Inverse Problems" => Any[
         # Inverse problems introduction.
diff --git a/docs/src/catalyst_applications/dynamical_systems.md b/docs/src/steady_state_functionality/dynamical_systems.md
similarity index 84%
rename from docs/src/catalyst_applications/dynamical_systems.md
rename to docs/src/steady_state_functionality/dynamical_systems.md
index 41a89862c7..fcd56129e0 100644
--- a/docs/src/catalyst_applications/dynamical_systems.md
+++ b/docs/src/steady_state_functionality/dynamical_systems.md
@@ -1,10 +1,10 @@
 # [Analysing model steady state properties with DynamicalSystems.jl](@id dynamical_systems)
-The [DynamicalSystems.jl package](https://github.com/JuliaDynamics/DynamicalSystems.jl) implements a wide range of methods for analysing dynamical systems. This includes both continuous-time systems (i.e. ODEs) and discrete-times ones (difference equations, these are, however, not relevant to chemical reaction network modelling). Here we give two examples of how DynamicalSystems.jl can be used, with the package's [documentation describing many more features](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/tutorial/). Finally, it should also be noted that DynamicalSystems.jl contain several tools for [analysing data measured from dynamical systems](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/contents/#Exported-submodules).
+The [DynamicalSystems.jl package](https://github.com/JuliaDynamics/DynamicalSystems.jl) implements a wide range of methods for analysing dynamical systems. This includes both continuous-time systems (i.e. ODEs) and discrete-times ones (difference equations, however, these are not relevant to chemical reaction network modelling). Here we give two examples of how DynamicalSystems.jl can be used, with the package's [documentation describing many more features](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/tutorial/). Finally, it should also be noted that DynamicalSystems.jl contain several tools for [analysing data measured from dynamical systems](https://juliadynamics.github.io/DynamicalSystemsDocs.jl/dynamicalsystems/dev/contents/#Exported-submodules).
 
 ## [Finding basins of attraction](@id dynamical_systems_basins_of_attraction)
 Given enough time, an ODE will eventually reach a so-called [*attractor*](https://en.wikipedia.org/wiki/Attractor). For chemical reaction networks (CRNs), this will typically be either a *steady state* or a *limit cycle*. Since ODEs are deterministic, which attractor a simulation will reach is uniquely determined by the initial condition (assuming parameter values are fixed). Conversely, each attractor is associated with a set of initial conditions such that model simulations originating in these will tend to that attractor. These sets are called *basins of attraction*. Here, phase space (the space of all possible states of the system) can be divided into a number of basins of attraction equal to the number of attractors. 
 
-DynamicalSystems.jl provides a simple interface for finding an ODE's basins of attraction across any given subspace of phase space. In this example we will use the bistable [Wilhelm model](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) (which steady states we have previous [computed using homotopy continuation](@ref homotopy_continuation_basic_example)). As a first step, we create a `ODEProblem` corresponding to the model which basins of attraction we wish to compute. For this application, `u0` and `tspan` is unused, and their values are of little importance (the only exception is than `tspan`, for implementation reason, must provide a not to small interval, we recommend minimum `(0.0, 1.0)`). 
+DynamicalSystems.jl provides a simple interface for finding an ODE's basins of attraction across any given subspace of phase space. In this example we will use the bistable [Wilhelm model](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) (which steady states we have previous [computed using homotopy continuation](@ref homotopy_continuation_basic_example)). As a first step, we create an `ODEProblem` corresponding to the model which basins of attraction we wish to compute. For this application, `u0` and `tspan` is unused, and their values are of little importance (the only exception is than `tspan`, for implementation reason, must provide a not too small interval, we recommend minimum `(0.0, 1.0)`). 
 ```@example dynamical_systems_basins
 using Catalyst
 wilhelm_model = @reaction_network begin
@@ -54,7 +54,7 @@ More information on how to compute basins of attractions for ODEs using Dynamica
 
 While Lyapunov exponents can be used for other purposes, they are primarily used to characterise [*chaotic behaviours*](https://en.wikipedia.org/wiki/Chaos_theory) (where small changes in initial conditions has large effect on the resulting trajectories). Generally, an ODE exhibit chaotic behaviour if its attractor(s) have *at least one* positive Lyapunov exponent. Practically, Lyapunov exponents can be computed using DynamicalSystems.jl's `lyapunovspectrum` function. Here we will use it to investigate two models, one which exhibits chaos and one which do not.
 
-First, let us consider the [Willamowski–Rössler](https://www.degruyter.com/document/doi/10.1515/zna-1980-0308/html?lang=en) model, which is known to exhibit chaotic behaviour.
+First, let us consider the [Willamowski–Rössler model](@ref ref), which is known to exhibit chaotic behaviour.
 ```@example dynamical_systems_lyapunov
 using Catalyst
 wr_model = @reaction_network begin
@@ -82,14 +82,14 @@ plot(sol; idxs=(:X, :Y, :Z))
 Next, like when we [computed basins of attraction](@ref dynamical_systems_basins_of_attraction), we create a `CoupledODEs` corresponding to the model and state for which we wish to compute our Lyapunov spectrum. Lke previously, `tspan` must provide some small interval (at least `(0.0, 1.0)` is recommended), but else have no impact on the computed Lyapunov spectrum.
 ```@example dynamical_systems_lyapunov
 using DynamicalSystems
-ds = CoupledODEs(oprob, (alg = Rodas5P(autodiff=false),))
+ds = CoupledODEs(oprob, (alg = Rodas5P(autodiff = false),))
 nothing # hide
 ```
-Here, the `autodiff=false` argument is required when Lyapunov spectrums are computed. We can now provide our `CoupledODEs` (`ds`) to `lyapunovspectrum` to compute the lyapunov spectrum. This function requires a second argument (here set to `100`). Generally setting this to a higher value will increase accuracy, but also increase runtime (since `lyapunovspectrum` is fast for most systems, setting this to a large value is recommended).
+Here, the `autodiff = false` argument is required when Lyapunov spectrums are computed. We can now provide our `CoupledODEs` (`ds`) to `lyapunovspectrum` to compute the lyapunov spectrum. This function requires a second argument (here set to `100`). Generally setting this to a higher value will increase accuracy, but also increase runtime (since `lyapunovspectrum` is fast for most systems, setting this to a large value is recommended).
 ```@example dynamical_systems_lyapunov
 lyapunovspectrum(ds, 100)
 ```
-Here, the largest exponent is positive, suggesting that the model is chaotic (or, more accurately, it has at least one chaotic attractor, to which it go to from the initial condition $(1.5, 1.5, 1.5)).
+Here, the largest exponent is positive, suggesting that the model is chaotic (or, more accurately, it has at least one chaotic attractor, to which is approached from the initial condition $(1.5,1.5,1.5)).
 
 Next, we consider the [Brusselator] model. First we simulate the model for two similar initial conditions, confirming that they converge to the same limit cycle:
 ```@example dynamical_systems_lyapunov
@@ -114,7 +114,7 @@ plot!(osol2; idxs = (:X, :Y))
 ```
 Next, we compute the Lyapunov spectrum at one of the initial conditions:
 ```@example dynamical_systems_lyapunov
-ds = CoupledODEs(oprob1, (alg = Rodas5P(autodiff=false),))
+ds = CoupledODEs(oprob1, (alg = Rodas5P(autodiff = false),))
 lyapunovspectrum(ds, 100)
 ```
 Here, all Lyapunov exponents are negative, confirming that the brusselator is non-chaotic.

From e4fe83e47095ce3851cb5428faaeb9f991c1ee61 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Fri, 24 May 2024 17:58:56 -0400
Subject: [PATCH 17/17] add missing comma

---
 docs/pages.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/pages.jl b/docs/pages.jl
index b54e11d80e..c39c87a771 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -38,7 +38,7 @@ pages = Any[
         "steady_state_functionality/homotopy_continuation.md",
         "steady_state_functionality/nonlinear_solve.md",
         "steady_state_functionality/steady_state_stability_computation.md",        
-        "steady_state_functionality/bifurcation_diagrams.md"
+        "steady_state_functionality/bifurcation_diagrams.md",
         "steady_state_functionality/dynamical_systems.md"
     ],
     "Inverse Problems" => Any[