From 5b97ea53cd5d0ba0038a2b0e2985cddd29d7e75f Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Tue, 16 Jan 2024 15:26:54 -0500
Subject: [PATCH 01/10] update file

---
 docs/src/model_creation/dsl_description.md | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/docs/src/model_creation/dsl_description.md b/docs/src/model_creation/dsl_description.md
index c888b4dda9..e32bba7673 100644
--- a/docs/src/model_creation/dsl_description.md
+++ b/docs/src/model_creation/dsl_description.md
@@ -2,7 +2,7 @@
 This tutorial describes the syntax for building chemical reaction
 network models using Catalyst's domain-specific language (DSL). Examples showing
 how to both construct and solve ODE, SDE, and jump models are provided in [Basic
-Chemical Reaction Network Examples](@ref basic_CRN_examples). To learn more about the symbolic
+Chemical Reaction Network Examples](@ref basic_CRN_library). To learn more about the symbolic
 [`ReactionSystem`](@ref)s generated by the DSL, and how to use them directly, see
 the tutorial on [Programmatic Construction of Symbolic Reaction Systems](@ref programmatic_CRN_construction).
 
@@ -80,8 +80,7 @@ u0 = [:X => 1.0, :Y => 1.0, :XY => 1.0, :Z1 => 1.0, :Z2 => 1.0]
 oprob = ODEProblem(rn, u0, tspan, [])
 ```
 
-For more detailed examples, see the [Basic Chemical Reaction Network
-Examples](@ref basic_CRN_examples).
+For more detailed examples, see the example implementations in the [Library of Basic Chemical Reaction Network Models](@ref basic_CRN_library) section.
 
 ## Defining parameters and species
 Numeric parameter values do not need to be set when the model is created, i.e.

From 9356ed155b4c442484ef3b36b305b78d689ae831 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Tue, 16 Jan 2024 15:27:07 -0500
Subject: [PATCH 02/10] up

---
 .../examples/basic_CRN_examples.md            | 325 +++++++++++++-----
 1 file changed, 239 insertions(+), 86 deletions(-)

diff --git a/docs/src/model_creation/examples/basic_CRN_examples.md b/docs/src/model_creation/examples/basic_CRN_examples.md
index 3036f3a686..762761cbbc 100644
--- a/docs/src/model_creation/examples/basic_CRN_examples.md
+++ b/docs/src/model_creation/examples/basic_CRN_examples.md
@@ -1,121 +1,274 @@
-# [Basic Chemical Reaction Network Examples](@id basic_CRN_examples)
+# [Library of Basic Chemical Reaction Network Models](@id basic_CRN_library)
+Below we will present various simple and established chemical reaction network (CRN) models. Each model is given some brief background, implemented using the `@reaction_network` DSL, and basic simulations are performed.
 
-## Example: Birth-death process
-
-```@example bcrn1
-using Catalyst, DifferentialEquations, Plots
+## Birth-death process
+The birth-death process is one of the simplest possible CRN models. It consists of a single component ($X$) which is both produced and degraded at linear rates:
+```@example crn_library_birth_death
+using Catalyst
+bd_process = @reaction_network begin
+    (p,d), ∅ <--> X
+end
+```
+Next we define simulation conditions. Note that the initial condition is integer-valued (required to perform jump simulations).
+```@example crn_library_birth_death
+u0 = [:X => 1]
+tspan = (0.0, 10.0)
+ps = [:p => 1, :d => 0.2]
+```
+We can now simulate our model using all three interpretations. First we perform a reaction rate equation-based ODE simulation:
+```@example crn_library_birth_death
+using OrdinaryDiffEq
+oprob = ODEProblem(bd_process, u0, tspan, ps)
+osol = solve(oprob, Tsit5())
+nothing # hide
+```
+Next, a chemical Langevin equation-based SDE simulation:
+```@example crn_library_birth_death
+using StochasticDiffEq
+sprob = SDEProblem(bd_process, u0, tspan, ps)
+ssol = solve(sprob, ImplicitEM())
+nothing # hide
+```
+Next, a stochastic chemical kinetics-based jump simulation:
+```@example crn_library_birth_death
+using JumpProcesses
+dprob = DiscreteProblem(bd_process, u0, tspan, ps)
+jprob = JumpProblem(bd_process, dprob, Direct())
+jsol = solve(jprob, SSAStepper())
+nothing # hide
+```
+Finally, we plot the results:
+```@example crn_library_birth_death
+using Plots
+oplt = plot(osol; title = "Reaction rate equation (ODE)")
+splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
+jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
+plot(oplt, splt, jplt; size=(800,700), layout = (3,1))
+```
 
-rs = @reaction_network begin
-  c1, X --> 2X
-  c2, X --> 0
-  c3, 0 --> X
+## Michaelis-Menten enzyme kinetics
+[Michaelis-Menten enzyme kinetics](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics) is a simple description of an enzyme ($E$) transforming a substrate ($S$) into a product ($P$). Under certain assumptions it can be simplified to a singe function (a Michaelis-Menten function) and used as a reaction rate. Here we instead present the full system model:
+```@example crn_library_michaelis_menten
+using Catalyst
+mm_system = @reaction_network begin
+  kB, S + E --> SE
+  kD, SE --> S + E
+  kP, SE --> P + E
 end
-p = (:c1 => 1.0, :c2 => 2.0, :c3 => 50.)
-tspan = (0.,4.)
-u0 = [:X => 5.]
+```
+Next, we perform ODE, SDE, and jump simulations of the model:
+```@example crn_library_michaelis_menten
+u0 = [:S => 301, :E => 100, :SE => 0, :P => 0]
+tspan = (0., 100.)
+ps = [:kB => 0.00166, :kD => 0.0001, :kP => 0.1]
 
-# solve ODEs
-oprob = ODEProblem(rs, u0, tspan, p)
+using OrdinaryDiffEq
+oprob = ODEProblem(mm_system, u0, tspan, ps)
 osol  = solve(oprob, Tsit5())
 
-# solve for Steady states
-ssprob = SteadyStateProblem(rs, u0, p)
-sssol  = solve(ssprob, SSRootfind())
+using StochasticDiffEq
+sprob = SDEProblem(mm_system, u0, tspan, ps)
+ssol = solve(sprob, ImplicitEM())
 
-# solve SDEs
-sprob = SDEProblem(rs, u0, tspan, p)
-ssol  = solve(sprob, EM(), dt=.01)
-
-# solve JumpProblem
-u0 = [:X => 5]
-dprob = DiscreteProblem(rs, u0, tspan, p)
-jprob = JumpProblem(rs, dprob, Direct())
+using JumpProcesses
+dprob = DiscreteProblem(mm_system, u0, tspan, ps)
+jprob = JumpProblem(mm_system, dprob, Direct())
 jsol = solve(jprob, SSAStepper())
 
-plot(plot(osol; title = "Reaction Rate Equation ODEs"),
-     plot(ssol; title = "Chemical Langevin SDEs"),
-     plot(jsol; title = "Stochastic Chemical Kinetics Jump Processes");
-     layout = (3, 1))
+using Plots
+oplt = plot(osol; title = "Reaction rate equation (ODE)")
+splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
+jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
+plot(oplt, splt, jplt; size=(800,700), layout = (3,1))
 ```
 
-## Example: Michaelis-Menten enzyme kinetics
-
-```@example bcrn2
-using Catalyst, DifferentialEquations, Plots
-
-rs = @reaction_network begin
-  c1, S + E --> SE
-  c2, SE --> S + E
-  c3, SE --> P + E
+## SIR infection model
+The [SIR model](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model) is the simplest model of the spread of an infectious disease. While the real system is very different from the chemical and cellular processes typically modelled with CRNs, it (and several other epidemiological systems) can be modelled using the same CRN formalism. The SIR model consists of three species: susceptible ($S$), infected ($I$), and removed ($R$) individuals, and two reaction events: infection and recovery.
+```@example crn_library_sir
+using Catalyst
+sir_model = @reaction_network begin
+    α, S + I --> 2I
+    β, I --> R
 end
-p = (:c1 => 0.00166, :c2 => 0.0001, :c3 => 0.1)
-tspan = (0., 100.)
-u0 = [:S => 301., :E => 100., :SE => 0., :P => 0.]
+```
+First we perform a deterministic ODE simulation:
+```@example crn_library_sir
+using OrdinaryDiffEq, Plots
+u0 = [:S => 99, :I => 1, :R => 0]
+tspan = (0.0, 500.0)
+ps = [:α => 0.001, :β => 0.01]
 
-# solve ODEs
-oprob = ODEProblem(rs, u0, tspan, p)
-osol  = solve(oprob, Tsit5())
+# Solve ODEs.
+oprob = ODEProblem(sir_model, u0, tspan, ps)
+osol = solve(oprob, Tsit5())
+plot(osol; title = "Reaction rate equation (ODE)")
+```
+Next we perform 3 different Jump simulations. Note that for the stochastic model, the occurrence of a outbreak is not certain. Rather, there is a possibility that it fizzles out without a noteworthy peak.
+```@example crn_library_sir
+using JumpProcesses
+dprob = DiscreteProblem(sir_model, u0, tspan, ps)
+jprob = JumpProblem(sir_model, dprob, Direct())
 
-# solve JumpProblem
-u0 = [:S => 301, :E => 100, :SE => 0, :P => 0]
-dprob = DiscreteProblem(rs, u0, tspan, p)
-jprob = JumpProblem(rs, dprob, Direct())
-jsol = solve(jprob, SSAStepper())
+jsol1 = solve(jprob, SSAStepper())
+jsol2 = solve(jprob, SSAStepper())
+jsol3 = solve(jprob, SSAStepper())
+jsol1 = solve(jprob, SSAStepper(); seed=1) # hide
+jsol2 = solve(jprob, SSAStepper(); seed=2) # hide
+jsol3 = solve(jprob, SSAStepper(); seed=3) # hide
 
-plot(plot(osol; title = "Reaction Rate Equation ODEs"),
-     plot(jsol; title = "Stochastic Chemical Kinetics Jump Processes");
-     layout = (2, 1))
+jplt1 = plot(jsol1; title = "Outbreak")
+jplt2 = plot(jsol2; title = "Outbreak")
+jplt3 = plot(jsol3; title = "No outbreak")
+plot(jplt1, jplt2, jplt3; size=(800,700), layout = (3,1))
 ```
 
-## Example: SIR infection model
-```@example bcrn3 
-using Catalyst, DifferentialEquations, Plots
+## The Wilhelm model
+The Wilhelm model was introduced in [*Wilhelm (2009)*](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) as the smallest CRN model (with constant rates) that exhibits bistability.
+```@example crn_library_wilhelm
+wilhelm_model = @reaction_network begin
+    k1, Y --> 2X
+    k2, 2X --> X + Y
+    k3, X + Y --> Y
+    k4, X --> 0
+end
+```
+We can simulate the model for two different initial conditions, demonstrating the existence of two different stable steady states. 
+```@example crn_library_wilhelm
+using OrdinaryDiffEq, Plots
+u0_1 = [:X => 1.5, :Y => 0.5]
+u0_2 = [:X => 2.5, :Y => 0.5]
+tspan = (0., 10.)
+ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5]
 
-rs = @reaction_network begin
-    α, S + I --> 2I
-    β, I --> R
+oprob1 = ODEProblem(wilhelm_model, u0_1, tspan, ps)
+oprob2 = ODEProblem(wilhelm_model, u0_2, tspan, ps)
+osol1 = solve(oprob1, Tsit5())
+osol2 = solve(oprob2, Tsit5())
+oplt1 = plot(osol1; idxs = :X, label = "X(0) = 1.5")
+oplt2 = plot!(osol2; idxs = :X, label = "X(0) = 2.5", yguide = "X", size = (800,700))
+```
+
+## Simple self-activation loop
+The simplest self-activation loop consist of a single species (here called $X$) which activates its own production. If its production rate is modelled with a hill function with $n>1$, the system may exhibit bistability.
+```@example crn_library_self_activation
+using Catalyst
+sa_loop = @reaction_network begin
+    v₀ + hill(X,v,K,n), ∅ --> X
+    d, X --> ∅
 end
-p = [:α => .1/100, :β => .01]
-tspan = (0.0,500.0)
-u0 = [:S => 99.0, :I => 1.0, :R => 0.0]
+```
+A simple example of such a loop is a transcription factor which activates its own gene. Here, $v₀$ represents a basic transcription rate (leakage) in the absence of the transcription factor. 
 
-# Solve ODEs.
-oprob = ODEProblem(rs, u0, tspan, p)
-osol = solve(oprob)
+We simulate the self-activation loop from a single initial condition using both deterministic (ODE) and stochastic (jump) simulations. We note that while the deterministic simulation reaches a single steady state, the stochastic one switches between two different states.
+```@example crn_library_self_activation
+using OrdinaryDiffEq, Plots
+u0 = [:X => 4]
+tspan = (0.0, 1000.0)
+ps = [:v₀ => 0.1, :v => 2.0, :K => 10.0, :n => 2, :d => 0.1]
 
-# Solve Jumps.
-dprob = DiscreteProblem(rs, u0, tspan, p)
-jprob = JumpProblem(rs, dprob, Direct())
+oprob = ODEProblem(sa_loop, u0, tspan, ps)
+osol = solve(oprob, Tsit5())
+
+dprob = DiscreteProblem(sa_loop, u0, tspan, ps)
+jprob = JumpProblem(sa_loop, dprob, Direct())
 jsol = solve(jprob, SSAStepper())
+jsol = solve(jprob, SSAStepper(); seed = 2091) # hide
 
-plot(plot(osol; title = "Reaction Rate Equation ODEs"),
-     plot(jsol; title = "Stochastic Chemical Kinetics Jump Processes");
-     layout = (2, 1))
+plot(osol; label = "Reaction rate equation (ODE)")
+plot!(jsol; label = "Stochastic chemical kinetics (Jump)", yguide = "X", size = (800,600))
 ```
 
-## Example: Brusselator chemical reaction network
-```@example bcrn4
-using Catalyst, DifferentialEquations, Plots
-
-rs = @reaction_network begin
-    @parameters A B
+## The Brusselator
+The [Brusselator](https://en.wikipedia.org/wiki/Brusselator) is a well known (theoretical) CRN model able to produce oscillations (its name is a portmanteau of "Brussels" and "oscillator").
+```@example crn_library_brusselator
+using Catalyst
+brusselator = @reaction_network begin
     A, ∅ --> X
     1, 2X + Y --> 3X
     B, X --> Y
     1, X --> ∅
 end
-tspan = (0.0,50.0)
+```
+It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this results is based on models generated when *combinatorial adjustment of rates is not performed*. Since Catalyst automatically perform these adjustments, and one reaction contain a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in he second case the system is able to generate oscillations.
+```@example crn_library_brusselator
+using OrdinaryDiffEq, Plots
 u0 = [:X => 1.0, :Y => 1.0]
+tspan = (0., 50.)
+ps1 = [:A => 1.0, :B => 1.0]
+ps2 = [:A => 1.0, :B => 1.8]
 
-# Non-oscillation parameter set
-oprob1 = ODEProblem(rs, u0, tspan, [:A => 1.0, :B => 1.0])
-osol1 = solve(oprob1)
+oprob1 = ODEProblem(brusselator, u0, tspan, ps1)
+oprob2 = ODEProblem(brusselator, u0, tspan, ps2)
+osol1  = solve(oprob1, Rodas5P())
+osol2  = solve(oprob2, Rodas5P())
+oplt1 = plot(osol1; title = "No Oscillation")
+oplt2 = plot(osol2; title = "Oscillation")
 
-# Oscillation parameter set
-oprob2 = ODEProblem(rs, u0, tspan, [:A => 1.0, :B => 3.0])
-osol2 = solve(oprob2)
+plot(oplt1, oplt2; layout = (1,2), size(800,700))
+```
+
+## The repressilator
+The repressilator was introduced in [*Elowitz & Leibler (2000)*](https://www.nature.com/articles/35002125) as a simple system that is able to generate oscillations (most notably, they demonstrated this both in a model and in a synthetic in vivo implementation in *Escherichia col*). It consists of three genes, repressing each other in a cycle. Here, we will implement it using three species ($X$, $Y$, and $Z$) which production rates are (repressing) [Hill functions](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)).
+```@example crn_library_brusselator
+using Catalyst
+repressilator = @reaction_network begin
+    hillr(Z,v,K,n), ∅ --> X
+    hillr(X,v,K,n), ∅ --> Y
+    hillr(Y,v,K,n), ∅ --> Z
+    d, (X, Y, Z) --> ∅
+end
+```
+Whether it oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model is able to sustain oscillations in both cases. This is an example of the phenomena of *noise-induced oscillation*.
+```@example crn_library_brusselator
+using OrdinaryDiffEq, StochasticDiffEq, Plots
+u0 = [:X => 50.0, :Y => 15.0, :Z => 15.0]
+tspan = (0., 200.)
+ps1 = [:v => 10.0, :K => 20.0, :n => 3, :d => 0.1]
+ps2 = [:v => 10.0, :K => 50.0, :n => 3, :d => 0.1]
+
+oprob1 = ODEProblem(repressilator, u0, tspan, ps1)
+oprob2 = ODEProblem(repressilator, u0, tspan, ps2)
+osol1  = solve(oprob1, Tsit5())
+osol2  = solve(oprob2, Tsit5())
+oplt1 = plot(osol1; title = "Oscillation (ODE, K = 20)")
+oplt2 = plot(osol2; title = "No oscillation (ODE, K = 50)")
+
+sprob1 = SDEProblem(repressilator, u0, tspan, ps1)
+sprob2 = SDEProblem(repressilator, u0, tspan, ps2)
+ssol1  = solve(sprob1, ImplicitEM())
+ssol2  = solve(sprob2, ImplicitEM())
+ssol1  = solve(sprob1, ImplicitEM(); seed = 1) # hide
+ssol2  = solve(sprob2, ImplicitEM(); seed = 100) # hide
+splt1 = plot(ssol1; title = "Oscillation (SDE, K = 20)")
+splt2 = plot(ssol2; title = "Oscillation (SDE, K = 50)")
+
+plot(oplt1, oplt2, splt1, splt2; layout = (2,2), size = (800,600))
+```
+
+## The Willamowski–Rössler model
+The Willamowski–Rössler model was introduced in [*Willamowski & Rössler (1979)*](https://www.degruyter.com/document/doi/10.1515/zna-1980-0308/html?lang=en) as an example of a simple CRN model which exhibits [*chaotic behaviours*](https://en.wikipedia.org/wiki/Chaos_theory). This means that small changes in initial conditions can produce relatively large changes in the system's trajectory.
+```@example crn_library_chaos
+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
+```
+Here we first simulate the model for a single initial conditions, showing in both time-state space and phase space how how it reaches a [*strange attractor*](https://www.dynamicmath.xyz/strange-attractors/).
+```@example crn_library_chaos
+using OrdinaryDiffEq, Plots
+u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5]
+tspan = (0.0, 50.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(plot(osol1; title = "No oscillation (B < 1 + A^2)"),
-     plot(osol2; title = "Oscillation (B > 1 + A^2)");
-     layout = (2, 1))
+plt1 = plot(sol; title = "Time-state space")
+plt2 = plot(sol; idxs = (:X, :Y, :Z), title = "Phase space")
+plot(plt1, plt2; layout = (1,2), size = (800,400))
 ```
\ No newline at end of file

From fe29787dd950414cc73cd321309b3dec69875167 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Tue, 16 Jan 2024 15:27:23 -0500
Subject: [PATCH 03/10] better file name

---
 .../example_networks/basic_CRN_library.md     | 274 ++++++++++++++++++
 1 file changed, 274 insertions(+)
 create mode 100644 docs/src/catalyst_functionality/example_networks/basic_CRN_library.md

diff --git a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md b/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
new file mode 100644
index 0000000000..762761cbbc
--- /dev/null
+++ b/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
@@ -0,0 +1,274 @@
+# [Library of Basic Chemical Reaction Network Models](@id basic_CRN_library)
+Below we will present various simple and established chemical reaction network (CRN) models. Each model is given some brief background, implemented using the `@reaction_network` DSL, and basic simulations are performed.
+
+## Birth-death process
+The birth-death process is one of the simplest possible CRN models. It consists of a single component ($X$) which is both produced and degraded at linear rates:
+```@example crn_library_birth_death
+using Catalyst
+bd_process = @reaction_network begin
+    (p,d), ∅ <--> X
+end
+```
+Next we define simulation conditions. Note that the initial condition is integer-valued (required to perform jump simulations).
+```@example crn_library_birth_death
+u0 = [:X => 1]
+tspan = (0.0, 10.0)
+ps = [:p => 1, :d => 0.2]
+```
+We can now simulate our model using all three interpretations. First we perform a reaction rate equation-based ODE simulation:
+```@example crn_library_birth_death
+using OrdinaryDiffEq
+oprob = ODEProblem(bd_process, u0, tspan, ps)
+osol = solve(oprob, Tsit5())
+nothing # hide
+```
+Next, a chemical Langevin equation-based SDE simulation:
+```@example crn_library_birth_death
+using StochasticDiffEq
+sprob = SDEProblem(bd_process, u0, tspan, ps)
+ssol = solve(sprob, ImplicitEM())
+nothing # hide
+```
+Next, a stochastic chemical kinetics-based jump simulation:
+```@example crn_library_birth_death
+using JumpProcesses
+dprob = DiscreteProblem(bd_process, u0, tspan, ps)
+jprob = JumpProblem(bd_process, dprob, Direct())
+jsol = solve(jprob, SSAStepper())
+nothing # hide
+```
+Finally, we plot the results:
+```@example crn_library_birth_death
+using Plots
+oplt = plot(osol; title = "Reaction rate equation (ODE)")
+splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
+jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
+plot(oplt, splt, jplt; size=(800,700), layout = (3,1))
+```
+
+## Michaelis-Menten enzyme kinetics
+[Michaelis-Menten enzyme kinetics](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics) is a simple description of an enzyme ($E$) transforming a substrate ($S$) into a product ($P$). Under certain assumptions it can be simplified to a singe function (a Michaelis-Menten function) and used as a reaction rate. Here we instead present the full system model:
+```@example crn_library_michaelis_menten
+using Catalyst
+mm_system = @reaction_network begin
+  kB, S + E --> SE
+  kD, SE --> S + E
+  kP, SE --> P + E
+end
+```
+Next, we perform ODE, SDE, and jump simulations of the model:
+```@example crn_library_michaelis_menten
+u0 = [:S => 301, :E => 100, :SE => 0, :P => 0]
+tspan = (0., 100.)
+ps = [:kB => 0.00166, :kD => 0.0001, :kP => 0.1]
+
+using OrdinaryDiffEq
+oprob = ODEProblem(mm_system, u0, tspan, ps)
+osol  = solve(oprob, Tsit5())
+
+using StochasticDiffEq
+sprob = SDEProblem(mm_system, u0, tspan, ps)
+ssol = solve(sprob, ImplicitEM())
+
+using JumpProcesses
+dprob = DiscreteProblem(mm_system, u0, tspan, ps)
+jprob = JumpProblem(mm_system, dprob, Direct())
+jsol = solve(jprob, SSAStepper())
+
+using Plots
+oplt = plot(osol; title = "Reaction rate equation (ODE)")
+splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
+jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
+plot(oplt, splt, jplt; size=(800,700), layout = (3,1))
+```
+
+## SIR infection model
+The [SIR model](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model) is the simplest model of the spread of an infectious disease. While the real system is very different from the chemical and cellular processes typically modelled with CRNs, it (and several other epidemiological systems) can be modelled using the same CRN formalism. The SIR model consists of three species: susceptible ($S$), infected ($I$), and removed ($R$) individuals, and two reaction events: infection and recovery.
+```@example crn_library_sir
+using Catalyst
+sir_model = @reaction_network begin
+    α, S + I --> 2I
+    β, I --> R
+end
+```
+First we perform a deterministic ODE simulation:
+```@example crn_library_sir
+using OrdinaryDiffEq, Plots
+u0 = [:S => 99, :I => 1, :R => 0]
+tspan = (0.0, 500.0)
+ps = [:α => 0.001, :β => 0.01]
+
+# Solve ODEs.
+oprob = ODEProblem(sir_model, u0, tspan, ps)
+osol = solve(oprob, Tsit5())
+plot(osol; title = "Reaction rate equation (ODE)")
+```
+Next we perform 3 different Jump simulations. Note that for the stochastic model, the occurrence of a outbreak is not certain. Rather, there is a possibility that it fizzles out without a noteworthy peak.
+```@example crn_library_sir
+using JumpProcesses
+dprob = DiscreteProblem(sir_model, u0, tspan, ps)
+jprob = JumpProblem(sir_model, dprob, Direct())
+
+jsol1 = solve(jprob, SSAStepper())
+jsol2 = solve(jprob, SSAStepper())
+jsol3 = solve(jprob, SSAStepper())
+jsol1 = solve(jprob, SSAStepper(); seed=1) # hide
+jsol2 = solve(jprob, SSAStepper(); seed=2) # hide
+jsol3 = solve(jprob, SSAStepper(); seed=3) # hide
+
+jplt1 = plot(jsol1; title = "Outbreak")
+jplt2 = plot(jsol2; title = "Outbreak")
+jplt3 = plot(jsol3; title = "No outbreak")
+plot(jplt1, jplt2, jplt3; size=(800,700), layout = (3,1))
+```
+
+## The Wilhelm model
+The Wilhelm model was introduced in [*Wilhelm (2009)*](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) as the smallest CRN model (with constant rates) that exhibits bistability.
+```@example crn_library_wilhelm
+wilhelm_model = @reaction_network begin
+    k1, Y --> 2X
+    k2, 2X --> X + Y
+    k3, X + Y --> Y
+    k4, X --> 0
+end
+```
+We can simulate the model for two different initial conditions, demonstrating the existence of two different stable steady states. 
+```@example crn_library_wilhelm
+using OrdinaryDiffEq, Plots
+u0_1 = [:X => 1.5, :Y => 0.5]
+u0_2 = [:X => 2.5, :Y => 0.5]
+tspan = (0., 10.)
+ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5]
+
+oprob1 = ODEProblem(wilhelm_model, u0_1, tspan, ps)
+oprob2 = ODEProblem(wilhelm_model, u0_2, tspan, ps)
+osol1 = solve(oprob1, Tsit5())
+osol2 = solve(oprob2, Tsit5())
+oplt1 = plot(osol1; idxs = :X, label = "X(0) = 1.5")
+oplt2 = plot!(osol2; idxs = :X, label = "X(0) = 2.5", yguide = "X", size = (800,700))
+```
+
+## Simple self-activation loop
+The simplest self-activation loop consist of a single species (here called $X$) which activates its own production. If its production rate is modelled with a hill function with $n>1$, the system may exhibit bistability.
+```@example crn_library_self_activation
+using Catalyst
+sa_loop = @reaction_network begin
+    v₀ + hill(X,v,K,n), ∅ --> X
+    d, X --> ∅
+end
+```
+A simple example of such a loop is a transcription factor which activates its own gene. Here, $v₀$ represents a basic transcription rate (leakage) in the absence of the transcription factor. 
+
+We simulate the self-activation loop from a single initial condition using both deterministic (ODE) and stochastic (jump) simulations. We note that while the deterministic simulation reaches a single steady state, the stochastic one switches between two different states.
+```@example crn_library_self_activation
+using OrdinaryDiffEq, Plots
+u0 = [:X => 4]
+tspan = (0.0, 1000.0)
+ps = [:v₀ => 0.1, :v => 2.0, :K => 10.0, :n => 2, :d => 0.1]
+
+oprob = ODEProblem(sa_loop, u0, tspan, ps)
+osol = solve(oprob, Tsit5())
+
+dprob = DiscreteProblem(sa_loop, u0, tspan, ps)
+jprob = JumpProblem(sa_loop, dprob, Direct())
+jsol = solve(jprob, SSAStepper())
+jsol = solve(jprob, SSAStepper(); seed = 2091) # hide
+
+plot(osol; label = "Reaction rate equation (ODE)")
+plot!(jsol; label = "Stochastic chemical kinetics (Jump)", yguide = "X", size = (800,600))
+```
+
+## The Brusselator
+The [Brusselator](https://en.wikipedia.org/wiki/Brusselator) is a well known (theoretical) CRN model able to produce oscillations (its name is a portmanteau of "Brussels" and "oscillator").
+```@example crn_library_brusselator
+using Catalyst
+brusselator = @reaction_network begin
+    A, ∅ --> X
+    1, 2X + Y --> 3X
+    B, X --> Y
+    1, X --> ∅
+end
+```
+It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this results is based on models generated when *combinatorial adjustment of rates is not performed*. Since Catalyst automatically perform these adjustments, and one reaction contain a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in he second case the system is able to generate oscillations.
+```@example crn_library_brusselator
+using OrdinaryDiffEq, Plots
+u0 = [:X => 1.0, :Y => 1.0]
+tspan = (0., 50.)
+ps1 = [:A => 1.0, :B => 1.0]
+ps2 = [:A => 1.0, :B => 1.8]
+
+oprob1 = ODEProblem(brusselator, u0, tspan, ps1)
+oprob2 = ODEProblem(brusselator, u0, tspan, ps2)
+osol1  = solve(oprob1, Rodas5P())
+osol2  = solve(oprob2, Rodas5P())
+oplt1 = plot(osol1; title = "No Oscillation")
+oplt2 = plot(osol2; title = "Oscillation")
+
+plot(oplt1, oplt2; layout = (1,2), size(800,700))
+```
+
+## The repressilator
+The repressilator was introduced in [*Elowitz & Leibler (2000)*](https://www.nature.com/articles/35002125) as a simple system that is able to generate oscillations (most notably, they demonstrated this both in a model and in a synthetic in vivo implementation in *Escherichia col*). It consists of three genes, repressing each other in a cycle. Here, we will implement it using three species ($X$, $Y$, and $Z$) which production rates are (repressing) [Hill functions](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)).
+```@example crn_library_brusselator
+using Catalyst
+repressilator = @reaction_network begin
+    hillr(Z,v,K,n), ∅ --> X
+    hillr(X,v,K,n), ∅ --> Y
+    hillr(Y,v,K,n), ∅ --> Z
+    d, (X, Y, Z) --> ∅
+end
+```
+Whether it oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model is able to sustain oscillations in both cases. This is an example of the phenomena of *noise-induced oscillation*.
+```@example crn_library_brusselator
+using OrdinaryDiffEq, StochasticDiffEq, Plots
+u0 = [:X => 50.0, :Y => 15.0, :Z => 15.0]
+tspan = (0., 200.)
+ps1 = [:v => 10.0, :K => 20.0, :n => 3, :d => 0.1]
+ps2 = [:v => 10.0, :K => 50.0, :n => 3, :d => 0.1]
+
+oprob1 = ODEProblem(repressilator, u0, tspan, ps1)
+oprob2 = ODEProblem(repressilator, u0, tspan, ps2)
+osol1  = solve(oprob1, Tsit5())
+osol2  = solve(oprob2, Tsit5())
+oplt1 = plot(osol1; title = "Oscillation (ODE, K = 20)")
+oplt2 = plot(osol2; title = "No oscillation (ODE, K = 50)")
+
+sprob1 = SDEProblem(repressilator, u0, tspan, ps1)
+sprob2 = SDEProblem(repressilator, u0, tspan, ps2)
+ssol1  = solve(sprob1, ImplicitEM())
+ssol2  = solve(sprob2, ImplicitEM())
+ssol1  = solve(sprob1, ImplicitEM(); seed = 1) # hide
+ssol2  = solve(sprob2, ImplicitEM(); seed = 100) # hide
+splt1 = plot(ssol1; title = "Oscillation (SDE, K = 20)")
+splt2 = plot(ssol2; title = "Oscillation (SDE, K = 50)")
+
+plot(oplt1, oplt2, splt1, splt2; layout = (2,2), size = (800,600))
+```
+
+## The Willamowski–Rössler model
+The Willamowski–Rössler model was introduced in [*Willamowski & Rössler (1979)*](https://www.degruyter.com/document/doi/10.1515/zna-1980-0308/html?lang=en) as an example of a simple CRN model which exhibits [*chaotic behaviours*](https://en.wikipedia.org/wiki/Chaos_theory). This means that small changes in initial conditions can produce relatively large changes in the system's trajectory.
+```@example crn_library_chaos
+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
+```
+Here we first simulate the model for a single initial conditions, showing in both time-state space and phase space how how it reaches a [*strange attractor*](https://www.dynamicmath.xyz/strange-attractors/).
+```@example crn_library_chaos
+using OrdinaryDiffEq, Plots
+u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5]
+tspan = (0.0, 50.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())
+
+plt1 = plot(sol; title = "Time-state space")
+plt2 = plot(sol; idxs = (:X, :Y, :Z), title = "Phase space")
+plot(plt1, plt2; layout = (1,2), size = (800,400))
+```
\ No newline at end of file

From c07dbd9983cdc19e54f41d6c15f6930c39f90ea0 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Tue, 16 Jan 2024 16:25:28 -0500
Subject: [PATCH 04/10] up

---
 docs/Project.toml                             |  1 +
 .../example_networks/basic_CRN_library.md     | 50 +++++++++++--------
 2 files changed, 29 insertions(+), 22 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index e7454091d2..64dc3323e6 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -7,6 +7,7 @@ DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
+JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
 Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
diff --git a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md b/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
index 762761cbbc..02a147321c 100644
--- a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
+++ b/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
@@ -9,13 +9,14 @@ bd_process = @reaction_network begin
     (p,d), ∅ <--> X
 end
 ```
-Next we define simulation conditions. Note that the initial condition is integer-valued (required to perform jump simulations).
+Next, we define simulation conditions. Note that the initial condition is integer-valued (required to perform jump simulations).
 ```@example crn_library_birth_death
 u0 = [:X => 1]
 tspan = (0.0, 10.0)
 ps = [:p => 1, :d => 0.2]
+nothing # hide
 ```
-We can now simulate our model using all three interpretations. First we perform a reaction rate equation-based ODE simulation:
+We can now simulate our model using all three interpretations. First, we perform a reaction rate equation-based ODE simulation:
 ```@example crn_library_birth_death
 using OrdinaryDiffEq
 oprob = ODEProblem(bd_process, u0, tspan, ps)
@@ -43,11 +44,11 @@ using Plots
 oplt = plot(osol; title = "Reaction rate equation (ODE)")
 splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
 jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
-plot(oplt, splt, jplt; size=(800,700), layout = (3,1))
+plot(oplt, splt, jplt; lw = 3, size=(800,700), layout = (3,1))
 ```
 
 ## Michaelis-Menten enzyme kinetics
-[Michaelis-Menten enzyme kinetics](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics) is a simple description of an enzyme ($E$) transforming a substrate ($S$) into a product ($P$). Under certain assumptions it can be simplified to a singe function (a Michaelis-Menten function) and used as a reaction rate. Here we instead present the full system model:
+[Michaelis-Menten enzyme kinetics](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics) is a simple description of an enzyme ($E$) transforming a substrate ($S$) into a product ($P$). Under certain assumptions, it can be simplified to a single function (a Michaelis-Menten function) and used as a reaction rate. Here we instead present the full system model:
 ```@example crn_library_michaelis_menten
 using Catalyst
 mm_system = @reaction_network begin
@@ -79,7 +80,8 @@ using Plots
 oplt = plot(osol; title = "Reaction rate equation (ODE)")
 splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
 jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
-plot(oplt, splt, jplt; size=(800,700), layout = (3,1))
+plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1)) 
+plot!(bottom_margin = 3Plots.Measures.mm) # hide
 ```
 
 ## SIR infection model
@@ -91,7 +93,7 @@ sir_model = @reaction_network begin
     β, I --> R
 end
 ```
-First we perform a deterministic ODE simulation:
+First, we perform a deterministic ODE simulation:
 ```@example crn_library_sir
 using OrdinaryDiffEq, Plots
 u0 = [:S => 99, :I => 1, :R => 0]
@@ -101,9 +103,9 @@ ps = [:α => 0.001, :β => 0.01]
 # Solve ODEs.
 oprob = ODEProblem(sir_model, u0, tspan, ps)
 osol = solve(oprob, Tsit5())
-plot(osol; title = "Reaction rate equation (ODE)")
+plot(osol; title = "Reaction rate equation (ODE)", size=(800,350))
 ```
-Next we perform 3 different Jump simulations. Note that for the stochastic model, the occurrence of a outbreak is not certain. Rather, there is a possibility that it fizzles out without a noteworthy peak.
+Next, we perform 3 different Jump simulations. Note that for the stochastic model, the occurrence of an outbreak is not certain. Rather, there is a possibility that it fizzles out without a noteworthy peak.
 ```@example crn_library_sir
 using JumpProcesses
 dprob = DiscreteProblem(sir_model, u0, tspan, ps)
@@ -119,12 +121,13 @@ jsol3 = solve(jprob, SSAStepper(); seed=3) # hide
 jplt1 = plot(jsol1; title = "Outbreak")
 jplt2 = plot(jsol2; title = "Outbreak")
 jplt3 = plot(jsol3; title = "No outbreak")
-plot(jplt1, jplt2, jplt3; size=(800,700), layout = (3,1))
+plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1))
 ```
 
 ## The Wilhelm model
 The Wilhelm model was introduced in [*Wilhelm (2009)*](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) as the smallest CRN model (with constant rates) that exhibits bistability.
 ```@example crn_library_wilhelm
+using Catalyst
 wilhelm_model = @reaction_network begin
     k1, Y --> 2X
     k2, 2X --> X + Y
@@ -144,12 +147,13 @@ oprob1 = ODEProblem(wilhelm_model, u0_1, tspan, ps)
 oprob2 = ODEProblem(wilhelm_model, u0_2, tspan, ps)
 osol1 = solve(oprob1, Tsit5())
 osol2 = solve(oprob2, Tsit5())
-oplt1 = plot(osol1; idxs = :X, label = "X(0) = 1.5")
-oplt2 = plot!(osol2; idxs = :X, label = "X(0) = 2.5", yguide = "X", size = (800,700))
+plot(osol1; lw = 4, idxs = :X, label = "X(0) = 1.5")
+plot!(osol2; lw = 4, idxs = :X, label = "X(0) = 2.5", yguide = "X", size = (800,350))
+plot!(bottom_margin = 3Plots.Measures.mm) # hide
 ```
 
 ## Simple self-activation loop
-The simplest self-activation loop consist of a single species (here called $X$) which activates its own production. If its production rate is modelled with a hill function with $n>1$, the system may exhibit bistability.
+The simplest self-activation loop consists of a single species (here called $X$) which activates its own production. If its production rate is modelled with a hill function with $n>1$, the system may exhibit bistability.
 ```@example crn_library_self_activation
 using Catalyst
 sa_loop = @reaction_network begin
@@ -161,7 +165,7 @@ A simple example of such a loop is a transcription factor which activates its ow
 
 We simulate the self-activation loop from a single initial condition using both deterministic (ODE) and stochastic (jump) simulations. We note that while the deterministic simulation reaches a single steady state, the stochastic one switches between two different states.
 ```@example crn_library_self_activation
-using OrdinaryDiffEq, Plots
+using JumpProcesses, OrdinaryDiffEq, Plots
 u0 = [:X => 4]
 tspan = (0.0, 1000.0)
 ps = [:v₀ => 0.1, :v => 2.0, :K => 10.0, :n => 2, :d => 0.1]
@@ -174,12 +178,13 @@ jprob = JumpProblem(sa_loop, dprob, Direct())
 jsol = solve(jprob, SSAStepper())
 jsol = solve(jprob, SSAStepper(); seed = 2091) # hide
 
-plot(osol; label = "Reaction rate equation (ODE)")
-plot!(jsol; label = "Stochastic chemical kinetics (Jump)", yguide = "X", size = (800,600))
+plot(osol; lw = 3, label = "Reaction rate equation (ODE)")
+plot!(jsol; lw = 3, label = "Stochastic chemical kinetics (Jump)", yguide="X", size = (800,350))
+plot!(bottom_margin = 3Plots.Measures.mm, left_margin=3Plots.Measures.mm) # hide
 ```
 
 ## The Brusselator
-The [Brusselator](https://en.wikipedia.org/wiki/Brusselator) is a well known (theoretical) CRN model able to produce oscillations (its name is a portmanteau of "Brussels" and "oscillator").
+The [Brusselator](https://en.wikipedia.org/wiki/Brusselator) is a well-known (theoretical) CRN model able to produce oscillations (its name is a portmanteau of "Brussels" and "oscillator").
 ```@example crn_library_brusselator
 using Catalyst
 brusselator = @reaction_network begin
@@ -189,7 +194,7 @@ brusselator = @reaction_network begin
     1, X --> ∅
 end
 ```
-It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this results is based on models generated when *combinatorial adjustment of rates is not performed*. Since Catalyst automatically perform these adjustments, and one reaction contain a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in he second case the system is able to generate oscillations.
+It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this result is based on models generated when *combinatorial adjustment of rates is not performed*. Since Catalyst automatically perform these adjustments, and one reaction contains a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in the second case the system can generate oscillations.
 ```@example crn_library_brusselator
 using OrdinaryDiffEq, Plots
 u0 = [:X => 1.0, :Y => 1.0]
@@ -204,11 +209,11 @@ osol2  = solve(oprob2, Rodas5P())
 oplt1 = plot(osol1; title = "No Oscillation")
 oplt2 = plot(osol2; title = "Oscillation")
 
-plot(oplt1, oplt2; layout = (1,2), size(800,700))
+plot(oplt1, oplt2; lw = 3, layout = (2,1), size = (800,600))
 ```
 
 ## The repressilator
-The repressilator was introduced in [*Elowitz & Leibler (2000)*](https://www.nature.com/articles/35002125) as a simple system that is able to generate oscillations (most notably, they demonstrated this both in a model and in a synthetic in vivo implementation in *Escherichia col*). It consists of three genes, repressing each other in a cycle. Here, we will implement it using three species ($X$, $Y$, and $Z$) which production rates are (repressing) [Hill functions](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)).
+The repressilator was introduced in [*Elowitz & Leibler (2000)*](https://www.nature.com/articles/35002125) as a simple system that can generate oscillations (most notably, they demonstrated this both in a model and in a synthetic in vivo implementation in *Escherichia col*). It consists of three genes, repressing each other in a cycle. Here, we will implement it using three species ($X$, $Y$, and $Z$) whose production rates are (repressing) [Hill functions](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)).
 ```@example crn_library_brusselator
 using Catalyst
 repressilator = @reaction_network begin
@@ -218,7 +223,7 @@ repressilator = @reaction_network begin
     d, (X, Y, Z) --> ∅
 end
 ```
-Whether it oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model is able to sustain oscillations in both cases. This is an example of the phenomena of *noise-induced oscillation*.
+Whether it oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model can sustain oscillations in both cases. This is an example of the phenomena of *noise-induced oscillation*.
 ```@example crn_library_brusselator
 using OrdinaryDiffEq, StochasticDiffEq, Plots
 u0 = [:X => 50.0, :Y => 15.0, :Z => 15.0]
@@ -242,7 +247,7 @@ ssol2  = solve(sprob2, ImplicitEM(); seed = 100) # hide
 splt1 = plot(ssol1; title = "Oscillation (SDE, K = 20)")
 splt2 = plot(ssol2; title = "Oscillation (SDE, K = 50)")
 
-plot(oplt1, oplt2, splt1, splt2; layout = (2,2), size = (800,600))
+plot(oplt1, oplt2, splt1, splt2; lw = 2, layout = (2,2), size = (800,600))
 ```
 
 ## The Willamowski–Rössler model
@@ -259,7 +264,7 @@ wr_model = @reaction_network begin
     k7, Z --> ∅
 end
 ```
-Here we first simulate the model for a single initial conditions, showing in both time-state space and phase space how how it reaches a [*strange attractor*](https://www.dynamicmath.xyz/strange-attractors/).
+Here we first simulate the model for a single initial condition, showing in both time-state space and phase space how it reaches a [*strange attractor*](https://www.dynamicmath.xyz/strange-attractors/).
 ```@example crn_library_chaos
 using OrdinaryDiffEq, Plots
 u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5]
@@ -271,4 +276,5 @@ sol = solve(oprob, Rodas5P())
 plt1 = plot(sol; title = "Time-state space")
 plt2 = plot(sol; idxs = (:X, :Y, :Z), title = "Phase space")
 plot(plt1, plt2; layout = (1,2), size = (800,400))
+plot!(bottom_margin = 3Plots.Measures.mm) # hide
 ```
\ No newline at end of file

From fe62c5463e6886e651fcc296b76a3407ce8214cd Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Tue, 16 Jan 2024 16:40:56 -0500
Subject: [PATCH 05/10] format improvement

---
 .../example_networks/basic_CRN_library.md                       | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md b/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
index 02a147321c..7939a31543 100644
--- a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
+++ b/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
@@ -179,7 +179,7 @@ jsol = solve(jprob, SSAStepper())
 jsol = solve(jprob, SSAStepper(); seed = 2091) # hide
 
 plot(osol; lw = 3, label = "Reaction rate equation (ODE)")
-plot!(jsol; lw = 3, label = "Stochastic chemical kinetics (Jump)", yguide="X", size = (800,350))
+plot!(jsol; lw = 3, label = "Stochastic chemical kinetics (Jump)", yguide = "X", size = (800,350))
 plot!(bottom_margin = 3Plots.Measures.mm, left_margin=3Plots.Measures.mm) # hide
 ```
 

From 03de344c22d4d2670ad42ee390c254db53688a3d Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Wed, 17 Jan 2024 17:33:23 -0500
Subject: [PATCH 06/10] add two state model

---
 .../example_networks/basic_CRN_library.md     | 48 +++++++++++++++----
 1 file changed, 39 insertions(+), 9 deletions(-)

diff --git a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md b/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
index 7939a31543..0b0d247a65 100644
--- a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
+++ b/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
@@ -1,7 +1,7 @@
 # [Library of Basic Chemical Reaction Network Models](@id basic_CRN_library)
 Below we will present various simple and established chemical reaction network (CRN) models. Each model is given some brief background, implemented using the `@reaction_network` DSL, and basic simulations are performed.
 
-## Birth-death process
+## [Birth-death process](@id basic_CRN_library_bd)
 The birth-death process is one of the simplest possible CRN models. It consists of a single component ($X$) which is both produced and degraded at linear rates:
 ```@example crn_library_birth_death
 using Catalyst
@@ -47,7 +47,37 @@ jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
 plot(oplt, splt, jplt; lw = 3, size=(800,700), layout = (3,1))
 ```
 
-## Michaelis-Menten enzyme kinetics
+## [Two-state model](@id basic_CRN_library_two_states)
+The two-state model describes a component (here called $X$) which can exist in two different forms (here called $X₁$ and $X₂$). It switches between these forms at linear rates. First, we simulate the model using both ODEs and SDEs:
+```@example crn_library_two_states
+using Catalyst, OrdinaryDiffEq, StochasticDiffEq, Plots
+two_state_model = @reaction_network begin
+    (k1,k2), X₁ <--> X₂
+end
+
+u0 = [:X₁ => 50.0, :X₂ => 50.0]
+tspan = (0.0, 1.0)
+ps = [:k1 => 2.0, :k2 => 3.0]
+
+oprob = ODEProblem(two_state_model, u0, tspan, ps)
+osol = solve(oprob, Tsit5())
+oplt = plot(osol; title = "Reaction rate equation (ODE)")
+
+sprob = SDEProblem(two_state_model, u0, tspan, ps)
+ssol = solve(sprob, ImplicitEM())
+splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
+
+plot(oplt, splt; lw = 3, size = (800,550), layout = (2,1))
+```
+What is interesting about this model is that it has a *conserved quantity*, where $X₁ + X₂$ remains constant throughout the simulation (both in deterministic and stochastic cases). We can show this by instead plotting this conserved quantity.
+```@example crn_library_two_states
+@unpack X₁, X₂ = two_state_model
+oplt = plot(osol; idxs = X₁ + X₂, title = "Reaction rate equation (ODE)")
+splt = plot(ssol; idxs = X₁ + X₂, title = "Chemical Langevin equation (SDE)")
+plot(oplt, splt; lw = 3, size = (800,550), layout = (2,1))
+```
+
+## [Michaelis-Menten enzyme kinetics](@id basic_CRN_library_mm)
 [Michaelis-Menten enzyme kinetics](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics) is a simple description of an enzyme ($E$) transforming a substrate ($S$) into a product ($P$). Under certain assumptions, it can be simplified to a single function (a Michaelis-Menten function) and used as a reaction rate. Here we instead present the full system model:
 ```@example crn_library_michaelis_menten
 using Catalyst
@@ -84,7 +114,7 @@ plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1))
 plot!(bottom_margin = 3Plots.Measures.mm) # hide
 ```
 
-## SIR infection model
+## [SIR infection model](@id basic_CRN_library_sir)
 The [SIR model](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model) is the simplest model of the spread of an infectious disease. While the real system is very different from the chemical and cellular processes typically modelled with CRNs, it (and several other epidemiological systems) can be modelled using the same CRN formalism. The SIR model consists of three species: susceptible ($S$), infected ($I$), and removed ($R$) individuals, and two reaction events: infection and recovery.
 ```@example crn_library_sir
 using Catalyst
@@ -124,7 +154,7 @@ jplt3 = plot(jsol3; title = "No outbreak")
 plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1))
 ```
 
-## The Wilhelm model
+## [The Wilhelm model](@id basic_CRN_library_wilhelm)
 The Wilhelm model was introduced in [*Wilhelm (2009)*](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) as the smallest CRN model (with constant rates) that exhibits bistability.
 ```@example crn_library_wilhelm
 using Catalyst
@@ -152,7 +182,7 @@ plot!(osol2; lw = 4, idxs = :X, label = "X(0) = 2.5", yguide = "X", size = (800,
 plot!(bottom_margin = 3Plots.Measures.mm) # hide
 ```
 
-## Simple self-activation loop
+## [Simple self-activation loop](@id basic_CRN_library_self_activation)
 The simplest self-activation loop consists of a single species (here called $X$) which activates its own production. If its production rate is modelled with a hill function with $n>1$, the system may exhibit bistability.
 ```@example crn_library_self_activation
 using Catalyst
@@ -183,7 +213,7 @@ plot!(jsol; lw = 3, label = "Stochastic chemical kinetics (Jump)", yguide = "X",
 plot!(bottom_margin = 3Plots.Measures.mm, left_margin=3Plots.Measures.mm) # hide
 ```
 
-## The Brusselator
+## [The Brusselator](@id basic_CRN_library_brusselator)
 The [Brusselator](https://en.wikipedia.org/wiki/Brusselator) is a well-known (theoretical) CRN model able to produce oscillations (its name is a portmanteau of "Brussels" and "oscillator").
 ```@example crn_library_brusselator
 using Catalyst
@@ -209,10 +239,10 @@ osol2  = solve(oprob2, Rodas5P())
 oplt1 = plot(osol1; title = "No Oscillation")
 oplt2 = plot(osol2; title = "Oscillation")
 
-plot(oplt1, oplt2; lw = 3, layout = (2,1), size = (800,600))
+plot(oplt1, oplt2; lw = 3, size = (800,600), layout = (2,1))
 ```
 
-## The repressilator
+## [The repressilator](@id basic_CRN_library_)
 The repressilator was introduced in [*Elowitz & Leibler (2000)*](https://www.nature.com/articles/35002125) as a simple system that can generate oscillations (most notably, they demonstrated this both in a model and in a synthetic in vivo implementation in *Escherichia col*). It consists of three genes, repressing each other in a cycle. Here, we will implement it using three species ($X$, $Y$, and $Z$) whose production rates are (repressing) [Hill functions](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)).
 ```@example crn_library_brusselator
 using Catalyst
@@ -250,7 +280,7 @@ splt2 = plot(ssol2; title = "Oscillation (SDE, K = 50)")
 plot(oplt1, oplt2, splt1, splt2; lw = 2, layout = (2,2), size = (800,600))
 ```
 
-## The Willamowski–Rössler model
+## [The Willamowski–Rössler model](@id basic_CRN_library_wr)
 The Willamowski–Rössler model was introduced in [*Willamowski & Rössler (1979)*](https://www.degruyter.com/document/doi/10.1515/zna-1980-0308/html?lang=en) as an example of a simple CRN model which exhibits [*chaotic behaviours*](https://en.wikipedia.org/wiki/Chaos_theory). This means that small changes in initial conditions can produce relatively large changes in the system's trajectory.
 ```@example crn_library_chaos
 using Catalyst

From 2f7ef2c431d421b1cb9b7227c378b9893b2203f3 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Fri, 26 Jan 2024 10:02:06 -0500
Subject: [PATCH 07/10] up

---
 .../example_networks/basic_CRN_library.md     | 30 +++++++++++++++++++
 1 file changed, 30 insertions(+)

diff --git a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md b/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
index 0b0d247a65..5af2f482ac 100644
--- a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
+++ b/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
@@ -154,6 +154,36 @@ jplt3 = plot(jsol3; title = "No outbreak")
 plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1))
 ```
 
+## [Chemical cross coupling](@id basic_CRN_library_cc)
+In chemistry, [cross-coupling](https://en.wikipedia.org/wiki/Cross-coupling_reaction) is when a catalyst combines two substrates to form a product. In this example, the catalyst ($C$) first binds one substrate ($S₁$) to form an intermediary complex ($S₁Cat$). Next, the complex binds the second substrate ($S₂$) to form another complex ($CP$). Finally, the catalyst releases the now-formed product ($P$). This system is an extended version of the M[Michaelis-Menten system presented earlier](@ref basic_CRN_library_mm).
+```@example crn_library_cc
+using Catalyst
+cc_system = @reaction_network begin
+    k₁, S₁ + C --> S₁C
+    k₂, S₁C + S₂ --> CP
+    k₃, CP --> C + P
+end
+```
+Below, we perform a simple deterministic ODE simulation of teh system. Next, we plot both:
+- The concentration of the catalyst and the intermediaries.
+- The concentration of the substrates and the product.
+
+In two separate plots.
+```@example crn_library_cc
+using OrdinaryDiffEq, Plots
+u0 = [:S₁ => 1.0, :C => 0.05, :S₂ => 1.2, :S₁C => 0.0, :CP => 0.0, :P => 0.0]
+tspan = (0., 15.)
+ps = [:k₁ => 10.0, :k₂ => 5.0, :k₃ => 100.0] 
+
+# solve ODEs
+oprob = ODEProblem(cc_system, u0, tspan, ps)
+osol  = solve(oprob, Tsit5())
+
+plt1 = plot(osol; idxs = [:S₁, :S₂, :P], title = "Substrate and product dynamics")
+plt2 = plot(osol; idxs = [:C, :S₁C, :CP], title = "Catalyst and intermediaries dynamics")
+plot(plt1, plt2; lw = 3, size = (800,600), layout = (2,1))
+```
+
 ## [The Wilhelm model](@id basic_CRN_library_wilhelm)
 The Wilhelm model was introduced in [*Wilhelm (2009)*](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-3-90) as the smallest CRN model (with constant rates) that exhibits bistability.
 ```@example crn_library_wilhelm

From 00da3983af928f91dc0892b6647099efe918def9 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Fri, 26 Jan 2024 10:23:26 -0500
Subject: [PATCH 08/10] add cross-coupling example

---
 .../example_networks/basic_CRN_library.md              | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md b/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
index 5af2f482ac..ba3a4d4dd8 100644
--- a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
+++ b/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
@@ -154,8 +154,8 @@ jplt3 = plot(jsol3; title = "No outbreak")
 plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1))
 ```
 
-## [Chemical cross coupling](@id basic_CRN_library_cc)
-In chemistry, [cross-coupling](https://en.wikipedia.org/wiki/Cross-coupling_reaction) is when a catalyst combines two substrates to form a product. In this example, the catalyst ($C$) first binds one substrate ($S₁$) to form an intermediary complex ($S₁Cat$). Next, the complex binds the second substrate ($S₂$) to form another complex ($CP$). Finally, the catalyst releases the now-formed product ($P$). This system is an extended version of the M[Michaelis-Menten system presented earlier](@ref basic_CRN_library_mm).
+## [Chemical cross-coupling](@id basic_CRN_library_cc)
+In chemistry, [cross-coupling](https://en.wikipedia.org/wiki/Cross-coupling_reaction) is when a catalyst combines two substrates to form a product. In this example, the catalyst ($C$) first binds one substrate ($S₁$) to form an intermediary complex ($S₁C$). Next, the complex binds the second substrate ($S₂$) to form another complex ($CP$). Finally, the catalyst releases the now-formed product ($P$). This system is an extended version of the[Michaelis-Menten system presented earlier](@ref basic_CRN_library_mm).
 ```@example crn_library_cc
 using Catalyst
 cc_system = @reaction_network begin
@@ -164,16 +164,16 @@ cc_system = @reaction_network begin
     k₃, CP --> C + P
 end
 ```
-Below, we perform a simple deterministic ODE simulation of teh system. Next, we plot both:
-- The concentration of the catalyst and the intermediaries.
+Below, we perform a simple deterministic ODE simulation of the system. Next, we plot both:
 - The concentration of the substrates and the product.
+- The concentration of the catalyst and the intermediaries.
 
 In two separate plots.
 ```@example crn_library_cc
 using OrdinaryDiffEq, Plots
 u0 = [:S₁ => 1.0, :C => 0.05, :S₂ => 1.2, :S₁C => 0.0, :CP => 0.0, :P => 0.0]
 tspan = (0., 15.)
-ps = [:k₁ => 10.0, :k₂ => 5.0, :k₃ => 100.0] 
+ps = [:k₁ => 5.0, :k₂ => 5.0, :k₃ => 100.0] 
 
 # solve ODEs
 oprob = ODEProblem(cc_system, u0, tspan, ps)

From b84d285a40fc5ac9f6bd9075638e8270dee0b54e Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Thu, 16 May 2024 12:34:40 -0400
Subject: [PATCH 09/10] rebase

---
 .../example}/basic_CRN_library.md             | 67 ++++++++++---------
 1 file changed, 37 insertions(+), 30 deletions(-)
 rename docs/src/{catalyst_functionality/example_networks => model_creation/example}/basic_CRN_library.md (84%)

diff --git a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md b/docs/src/model_creation/example/basic_CRN_library.md
similarity index 84%
rename from docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
rename to docs/src/model_creation/example/basic_CRN_library.md
index ba3a4d4dd8..2ab3e52db8 100644
--- a/docs/src/catalyst_functionality/example_networks/basic_CRN_library.md
+++ b/docs/src/model_creation/example/basic_CRN_library.md
@@ -9,25 +9,26 @@ bd_process = @reaction_network begin
     (p,d), ∅ <--> X
 end
 ```
-Next, we define simulation conditions. Note that the initial condition is integer-valued (required to perform jump simulations).
+Next, we define simulation conditions. Note that the initial condition is integer-valued (more natural than decimal numbers for jump simulations).
 ```@example crn_library_birth_death
 u0 = [:X => 1]
 tspan = (0.0, 10.0)
-ps = [:p => 1, :d => 0.2]
+ps = [:p => 1.0, :d => 0.2]
 nothing # hide
 ```
 We can now simulate our model using all three interpretations. First, we perform a reaction rate equation-based ODE simulation:
 ```@example crn_library_birth_death
 using OrdinaryDiffEq
 oprob = ODEProblem(bd_process, u0, tspan, ps)
-osol = solve(oprob, Tsit5())
+osol = solve(oprob)
 nothing # hide
 ```
 Next, a chemical Langevin equation-based SDE simulation:
 ```@example crn_library_birth_death
 using StochasticDiffEq
 sprob = SDEProblem(bd_process, u0, tspan, ps)
-ssol = solve(sprob, ImplicitEM())
+ssol = solve(sprob, STrapezoid())
+ssol = solve(sprob, STrapezoid(); seed = 12) # hide
 nothing # hide
 ```
 Next, a stochastic chemical kinetics-based jump simulation:
@@ -36,6 +37,7 @@ using JumpProcesses
 dprob = DiscreteProblem(bd_process, u0, tspan, ps)
 jprob = JumpProblem(bd_process, dprob, Direct())
 jsol = solve(jprob, SSAStepper())
+jsol = solve(jprob, SSAStepper(); seed = 12) # hide
 nothing # hide
 ```
 Finally, we plot the results:
@@ -60,11 +62,12 @@ tspan = (0.0, 1.0)
 ps = [:k1 => 2.0, :k2 => 3.0]
 
 oprob = ODEProblem(two_state_model, u0, tspan, ps)
-osol = solve(oprob, Tsit5())
+osol = solve(oprob)
 oplt = plot(osol; title = "Reaction rate equation (ODE)")
 
 sprob = SDEProblem(two_state_model, u0, tspan, ps)
-ssol = solve(sprob, ImplicitEM())
+ssol = solve(sprob, STrapezoid())
+ssol = solve(sprob, STrapezoid(); seed = 12) # hide
 splt = plot(ssol; title = "Chemical Langevin equation (SDE)")
 
 plot(oplt, splt; lw = 3, size = (800,550), layout = (2,1))
@@ -74,8 +77,9 @@ What is interesting about this model is that it has a *conserved quantity*, wher
 @unpack X₁, X₂ = two_state_model
 oplt = plot(osol; idxs = X₁ + X₂, title = "Reaction rate equation (ODE)")
 splt = plot(ssol; idxs = X₁ + X₂, title = "Chemical Langevin equation (SDE)")
-plot(oplt, splt; lw = 3, size = (800,550), layout = (2,1))
+plot(oplt, splt; lw = 3, ylimit = (99,101), size = (800,450), layout = (2,1))
 ```
+Catalyst has special methods for working with conserved quantities, which are described [here](@ref ref).
 
 ## [Michaelis-Menten enzyme kinetics](@id basic_CRN_library_mm)
 [Michaelis-Menten enzyme kinetics](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics) is a simple description of an enzyme ($E$) transforming a substrate ($S$) into a product ($P$). Under certain assumptions, it can be simplified to a single function (a Michaelis-Menten function) and used as a reaction rate. Here we instead present the full system model:
@@ -95,16 +99,18 @@ ps = [:kB => 0.00166, :kD => 0.0001, :kP => 0.1]
 
 using OrdinaryDiffEq
 oprob = ODEProblem(mm_system, u0, tspan, ps)
-osol  = solve(oprob, Tsit5())
+osol  = solve(oprob)
 
 using StochasticDiffEq
 sprob = SDEProblem(mm_system, u0, tspan, ps)
-ssol = solve(sprob, ImplicitEM())
+ssol = solve(sprob, STrapezoid())
+ssol = solve(sprob, STrapezoid(); seed = 12) # hide
 
 using JumpProcesses
 dprob = DiscreteProblem(mm_system, u0, tspan, ps)
 jprob = JumpProblem(mm_system, dprob, Direct())
 jsol = solve(jprob, SSAStepper())
+jsol = solve(jprob, SSAStepper(); seed = 12) # hide
 
 using Plots
 oplt = plot(osol; title = "Reaction rate equation (ODE)")
@@ -113,6 +119,7 @@ jplt = plot(jsol; title = "Stochastic chemical kinetics (Jump)")
 plot(oplt, splt, jplt; lw = 2, size=(800,800), layout = (3,1)) 
 plot!(bottom_margin = 3Plots.Measures.mm) # hide
 ```
+Note that, due to the large amounts of the species involved, teh stochastic trajectories are very similar to the deterministic one.
 
 ## [SIR infection model](@id basic_CRN_library_sir)
 The [SIR model](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SIR_model) is the simplest model of the spread of an infectious disease. While the real system is very different from the chemical and cellular processes typically modelled with CRNs, it (and several other epidemiological systems) can be modelled using the same CRN formalism. The SIR model consists of three species: susceptible ($S$), infected ($I$), and removed ($R$) individuals, and two reaction events: infection and recovery.
@@ -132,7 +139,7 @@ ps = [:α => 0.001, :β => 0.01]
 
 # Solve ODEs.
 oprob = ODEProblem(sir_model, u0, tspan, ps)
-osol = solve(oprob, Tsit5())
+osol = solve(oprob)
 plot(osol; title = "Reaction rate equation (ODE)", size=(800,350))
 ```
 Next, we perform 3 different Jump simulations. Note that for the stochastic model, the occurrence of an outbreak is not certain. Rather, there is a possibility that it fizzles out without a noteworthy peak.
@@ -177,7 +184,7 @@ ps = [:k₁ => 5.0, :k₂ => 5.0, :k₃ => 100.0]
 
 # solve ODEs
 oprob = ODEProblem(cc_system, u0, tspan, ps)
-osol  = solve(oprob, Tsit5())
+osol  = solve(oprob)
 
 plt1 = plot(osol; idxs = [:S₁, :S₂, :P], title = "Substrate and product dynamics")
 plt2 = plot(osol; idxs = [:C, :S₁C, :CP], title = "Catalyst and intermediaries dynamics")
@@ -205,15 +212,15 @@ ps = [:k1 => 8.0, :k2 => 2.0, :k3 => 1.0, :k4 => 1.5]
 
 oprob1 = ODEProblem(wilhelm_model, u0_1, tspan, ps)
 oprob2 = ODEProblem(wilhelm_model, u0_2, tspan, ps)
-osol1 = solve(oprob1, Tsit5())
-osol2 = solve(oprob2, Tsit5())
+osol1 = solve(oprob1)
+osol2 = solve(oprob2)
 plot(osol1; lw = 4, idxs = :X, label = "X(0) = 1.5")
 plot!(osol2; lw = 4, idxs = :X, label = "X(0) = 2.5", yguide = "X", size = (800,350))
 plot!(bottom_margin = 3Plots.Measures.mm) # hide
 ```
 
 ## [Simple self-activation loop](@id basic_CRN_library_self_activation)
-The simplest self-activation loop consists of a single species (here called $X$) which activates its own production. If its production rate is modelled with a hill function with $n>1$, the system may exhibit bistability.
+The simplest self-activation loop consists of a single species (here called $X$) which activates its own production. If its production rate is modelled as a [Hill function](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)) with $n>1$, the system may exhibit bistability.
 ```@example crn_library_self_activation
 using Catalyst
 sa_loop = @reaction_network begin
@@ -231,12 +238,12 @@ tspan = (0.0, 1000.0)
 ps = [:v₀ => 0.1, :v => 2.0, :K => 10.0, :n => 2, :d => 0.1]
 
 oprob = ODEProblem(sa_loop, u0, tspan, ps)
-osol = solve(oprob, Tsit5())
+osol = solve(oprob)
 
 dprob = DiscreteProblem(sa_loop, u0, tspan, ps)
 jprob = JumpProblem(sa_loop, dprob, Direct())
 jsol = solve(jprob, SSAStepper())
-jsol = solve(jprob, SSAStepper(); seed = 2091) # hide
+jsol = solve(jprob, SSAStepper(); seed = 12) # hide
 
 plot(osol; lw = 3, label = "Reaction rate equation (ODE)")
 plot!(jsol; lw = 3, label = "Stochastic chemical kinetics (Jump)", yguide = "X", size = (800,350))
@@ -254,7 +261,7 @@ brusselator = @reaction_network begin
     1, X --> ∅
 end
 ```
-It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this result is based on models generated when *combinatorial adjustment of rates is not performed*. Since Catalyst automatically perform these adjustments, and one reaction contains a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in the second case the system can generate oscillations.
+It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this result is based on models generated when *combinatorial adjustment of rates is not performed*. Since Catalyst [automatically perform these adjustments](@ref ref), and one reaction contains a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in the second case the system can generate oscillations.
 ```@example crn_library_brusselator
 using OrdinaryDiffEq, Plots
 u0 = [:X => 1.0, :Y => 1.0]
@@ -264,16 +271,16 @@ ps2 = [:A => 1.0, :B => 1.8]
 
 oprob1 = ODEProblem(brusselator, u0, tspan, ps1)
 oprob2 = ODEProblem(brusselator, u0, tspan, ps2)
-osol1  = solve(oprob1, Rodas5P())
-osol2  = solve(oprob2, Rodas5P())
+osol1  = solve(oprob1)
+osol2  = solve(oprob2)
 oplt1 = plot(osol1; title = "No Oscillation")
 oplt2 = plot(osol2; title = "Oscillation")
 
 plot(oplt1, oplt2; lw = 3, size = (800,600), layout = (2,1))
 ```
 
-## [The repressilator](@id basic_CRN_library_)
-The repressilator was introduced in [*Elowitz & Leibler (2000)*](https://www.nature.com/articles/35002125) as a simple system that can generate oscillations (most notably, they demonstrated this both in a model and in a synthetic in vivo implementation in *Escherichia col*). It consists of three genes, repressing each other in a cycle. Here, we will implement it using three species ($X$, $Y$, and $Z$) whose production rates are (repressing) [Hill functions](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)).
+## [The Repressilator](@id basic_CRN_library_)
+The Repressilator was introduced in [*Elowitz & Leibler (2000)*](https://www.nature.com/articles/35002125) as a simple system that can generate oscillations (most notably, they demonstrated this both in a model and in a synthetic in vivo implementation in *Escherichia col*). It consists of three genes, repressing each other in a cycle. Here, we will implement it using three species ($X$, $Y$, and $Z$) whose production rates are (repressing) [Hill functions](https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)).
 ```@example crn_library_brusselator
 using Catalyst
 repressilator = @reaction_network begin
@@ -283,7 +290,7 @@ repressilator = @reaction_network begin
     d, (X, Y, Z) --> ∅
 end
 ```
-Whether it oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model can sustain oscillations in both cases. This is an example of the phenomena of *noise-induced oscillation*.
+Whether the Repressilator oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model can sustain oscillations in both cases. This is an example of the phenomena of *noise-induced oscillation*.
 ```@example crn_library_brusselator
 using OrdinaryDiffEq, StochasticDiffEq, Plots
 u0 = [:X => 50.0, :Y => 15.0, :Z => 15.0]
@@ -293,17 +300,17 @@ ps2 = [:v => 10.0, :K => 50.0, :n => 3, :d => 0.1]
 
 oprob1 = ODEProblem(repressilator, u0, tspan, ps1)
 oprob2 = ODEProblem(repressilator, u0, tspan, ps2)
-osol1  = solve(oprob1, Tsit5())
-osol2  = solve(oprob2, Tsit5())
+osol1  = solve(oprob1)
+osol2  = solve(oprob2)
 oplt1 = plot(osol1; title = "Oscillation (ODE, K = 20)")
 oplt2 = plot(osol2; title = "No oscillation (ODE, K = 50)")
 
 sprob1 = SDEProblem(repressilator, u0, tspan, ps1)
 sprob2 = SDEProblem(repressilator, u0, tspan, ps2)
-ssol1  = solve(sprob1, ImplicitEM())
-ssol2  = solve(sprob2, ImplicitEM())
-ssol1  = solve(sprob1, ImplicitEM(); seed = 1) # hide
-ssol2  = solve(sprob2, ImplicitEM(); seed = 100) # hide
+ssol1  = solve(sprob1, STrapezoid())
+ssol2  = solve(sprob2, STrapezoid())
+ssol1  = solve(sprob1, STrapezoid(); seed = 1) # hide
+ssol2  = solve(sprob2, STrapezoid(); seed = 100) # hide
 splt1 = plot(ssol1; title = "Oscillation (SDE, K = 20)")
 splt2 = plot(ssol2; title = "Oscillation (SDE, K = 50)")
 
@@ -324,14 +331,14 @@ wr_model = @reaction_network begin
     k7, Z --> ∅
 end
 ```
-Here we first simulate the model for a single initial condition, showing in both time-state space and phase space how it reaches a [*strange attractor*](https://www.dynamicmath.xyz/strange-attractors/).
+Here we simulate the model for a single initial condition, showing both time-state space and phase space how it reaches a [*strange attractor*](https://www.dynamicmath.xyz/strange-attractors/).
 ```@example crn_library_chaos
 using OrdinaryDiffEq, Plots
 u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5]
 tspan = (0.0, 50.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())
+sol = solve(oprob)
 
 plt1 = plot(sol; title = "Time-state space")
 plt2 = plot(sol; idxs = (:X, :Y, :Z), title = "Phase space")

From ed091cd3434b23521b0fa8646973e561fba8b73b Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Thu, 16 May 2024 16:55:09 -0400
Subject: [PATCH 10/10] up

---
 .../model_creation/{example => examples}/basic_CRN_library.md   | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)
 rename docs/src/model_creation/{example => examples}/basic_CRN_library.md (99%)

diff --git a/docs/src/model_creation/example/basic_CRN_library.md b/docs/src/model_creation/examples/basic_CRN_library.md
similarity index 99%
rename from docs/src/model_creation/example/basic_CRN_library.md
rename to docs/src/model_creation/examples/basic_CRN_library.md
index 2ab3e52db8..4955aa522d 100644
--- a/docs/src/model_creation/example/basic_CRN_library.md
+++ b/docs/src/model_creation/examples/basic_CRN_library.md
@@ -162,7 +162,7 @@ plot(jplt1, jplt2, jplt3; lw = 3, size=(800,700), layout = (3,1))
 ```
 
 ## [Chemical cross-coupling](@id basic_CRN_library_cc)
-In chemistry, [cross-coupling](https://en.wikipedia.org/wiki/Cross-coupling_reaction) is when a catalyst combines two substrates to form a product. In this example, the catalyst ($C$) first binds one substrate ($S₁$) to form an intermediary complex ($S₁C$). Next, the complex binds the second substrate ($S₂$) to form another complex ($CP$). Finally, the catalyst releases the now-formed product ($P$). This system is an extended version of the[Michaelis-Menten system presented earlier](@ref basic_CRN_library_mm).
+In chemistry, [cross-coupling](https://en.wikipedia.org/wiki/Cross-coupling_reaction) is when a catalyst combines two substrates to form a product. In this example, the catalyst ($C$) first binds one substrate ($S₁$) to form an intermediary complex ($S₁C$). Next, the complex binds the second substrate ($S₂$) to form another complex ($CP$). Finally, the catalyst releases the now-formed product ($P$). This system is an extended version of the [Michaelis-Menten system presented earlier](@ref basic_CRN_library_mm).
 ```@example crn_library_cc
 using Catalyst
 cc_system = @reaction_network begin