From 6b1fed7479e7d718cbc4d533b0894a5ba8b25395 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sun, 28 Apr 2024 15:44:15 -0400
Subject: [PATCH 1/4] save progress

---
 .../introduction_to_catalyst.md               | 430 ++++++------------
 1 file changed, 137 insertions(+), 293 deletions(-)

diff --git a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
index c0a36a5754..4261bd8ad7 100644
--- a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
+++ b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
@@ -1,301 +1,164 @@
 # [Introduction to Catalyst](@id introduction_to_catalyst)
-In this tutorial we provide an introduction to using Catalyst to specify
-chemical reaction networks, and then to solve ODE, jump, and SDE models
-generated from them. At the end we show what mathematical rate laws and
-transition rate functions (i.e. intensities or propensities) are generated by
-Catalyst for ODE, SDE and jump process models.
+This tutorial provides a basic introduction on how to create chemical reaction network (CRN) models using Catalyst, and how to perform ODE, SDE, and jump simulation using these. [An alternative introduction](@ref catalyst_for_new_julia_users) is available for users with little familiarity with Julia programming. 
 
-Let's start by using the Catalyst [`@reaction_network`](@ref) macro
-to specify a simple chemical reaction network: the well-known repressilator.
-
-We first import the basic packages we'll need:
-
-```@example tut1
-# If not already installed, first hit "]" within a Julia REPL. Then type:
-# add Catalyst DifferentialEquations Plots Latexify
-
-using Catalyst, DifferentialEquations, Plots, Latexify
-```
-
-We now construct the reaction network. The basic types of arrows and predefined
-rate laws one can use are discussed in detail within the tutorial, [The Reaction
-DSL](@ref dsl_description). Here, we use a mix of first order, zero order, and repressive Hill
-function rate laws. Note, $\varnothing$ corresponds to the empty state, and is
-used for zeroth order production and first order degradation reactions:
-```@example tut1
-repressilator = @reaction_network Repressilator begin
-    hillr(P₃,α,K,n), ∅ --> m₁
-    hillr(P₁,α,K,n), ∅ --> m₂
-    hillr(P₂,α,K,n), ∅ --> m₃
-    (δ,γ), m₁ <--> ∅
-    (δ,γ), m₂ <--> ∅
-    (δ,γ), m₃ <--> ∅
-    β, m₁ --> m₁ + P₁
-    β, m₂ --> m₂ + P₂
-    β, m₃ --> m₃ + P₃
-    μ, P₁ --> ∅
-    μ, P₂ --> ∅
-    μ, P₃ --> ∅
-end
-show(stdout, MIME"text/plain"(), repressilator) # hide
-```
-showing that we've created a new network model named `Repressilator` with the
-listed chemical species and unknowns. [`@reaction_network`](@ref) returns a
-[`ReactionSystem`](@ref), which we saved in the `repressilator` variable. It can
-be converted to a variety of other mathematical models represented as
-`ModelingToolkit.AbstractSystem`s, or analyzed in various ways using the
-[Catalyst.jl API](@ref). For example, to see the chemical species, parameters,
-and reactions we can use
-```@example tut1
-species(repressilator)
-```
-```@example tut1
-parameters(repressilator)
-```
-and
-```@example tut1
-reactions(repressilator)
-```
-We can also use Latexify to see the corresponding reactions in Latex, which shows what
-the `hillr` terms mathematically correspond to
+Before starting this tutorial, we must (unless this has already been done) [install and activate](@ref ref) the Catalyst package:
 ```julia
-latexify(repressilator)
+using Pkg
+Pkg.add("Catalyst")
+using Catalyst
+```
+
+## [Creating basic Catalyst models](@id introduction_to_catalyst_model_creation)
+Catalyst exports the [`@reaction_network`](@ref) [macro]([macro](https://docs.julialang.org/en/v1/manual/metaprogramming/#man-macros)), which provides the main way of creating CRN models (alternative approaches are described [later](@ref ref)). This is followed by a `begin ... end` block which encapsulates the reactions of the system. E.g here we create a simple [birth-death model](@ref ref) (where a single species $X$ is both produced and degraded):
+```@example intro_1
+using Catalyst # hide
+bd_model = @reaction_network begin
+ p, 0 --> X
+ d, X --> 0
+end
 ```
-```@example tut1
-repressilator #hide
+Next, we create a simple [Michaelis-Menten enzyme kinetics model](@ref ref) (where an enzyme, $E$ converts a substrate, $S$, to a product, $P$):
+```@example intro_1
+mm_model = @reaction_network begin
+ kB, S + E --> SE
+ kD, SE --> S + E
+ kP, SE --> P + E
+end
 ```
-Assuming [Graphviz](https://graphviz.org/) is installed and command line
-accessible, within a Jupyter notebook we can also graph the reaction network by
-```julia
-g = Graph(repressilator)
+And finally, a [SIR model of an infectious disease](@ref ref) (where susceptible individuals, $I$, get infected by infected individuals, $I$, and later recovers, $R$):
+```@example intro_1
+sir_model = @reaction_network begin
+    β, S + I --> 2I
+    γ, I --> R
+end
 ```
-giving
 
-![Repressilator solution](../assets/repressilator.svg)
+In all three cases, each reaction consists of:
+- A rate (at which the reaction occurs).
+- A set of substrates (species that are consumed by the reaction).
+- A set of products (species that are produced by the reaction).
 
-The network graph shows a variety of information, representing each species as a
-blue node, and each reaction as an orange dot. Black arrows from species to
-reactions indicate reactants, and are labelled with their input stoichiometry.
-Similarly, black arrows from reactions to species indicate products, and are
-labelled with their output stoichiometry. In contrast, red arrows from a species
-to reactions indicate the species is used within the reactions' rate
-expressions. For the repressilator, the reactions
-```julia
-hillr(P₃,α,K,n), ∅ --> m₁
-hillr(P₁,α,K,n), ∅ --> m₂
-hillr(P₂,α,K,n), ∅ --> m₃
-```
-have rates that depend on the proteins, and hence lead to red arrows from each
-`Pᵢ`.
+For reactions with multiple substrates/products (e.g. `S + E --> SE`), these are separated by `+`. Where substrates/products contain multiple units of the same species (e.g. `S + I --> 2I`), this is indicated by pre-appending that species with the corresponding [stoichiometric coefficient](https://en.wikipedia.org/wiki/Stoichiometry). The absences of substrates/products (that is, in pure production/degradation reactions like `0 --> X` and `X --> 0`) are denoted with `0`. A more throughout description of how to create models using `@reaction_network` (also called *the Catalyst DSL*) is provided [here](@ref ref).
 
-Note, from the REPL or scripts one can always use [`savegraph`](@ref) to save
-the graph (assuming `Graphviz` is installed).
+## [Simulating Catalyst models](@id introduction_to_catalyst_model_simulation)
+There exist three primary modes of simulation CRN models:
+1. Using ordinary differential equations (ODEs).
+2. Using stochastic differential equations (SDEs).
+3. Using jump process simulations.
 
-## Mass action ODE models
-Let's now use our `ReactionSystem` to generate and solve a corresponding mass
-action ODE model. We first convert the system to a `ModelingToolkit.ODESystem`
-by
-```@example tut1
-repressilator = complete(repressilator)
-odesys = convert(ODESystem, repressilator)
-```
-(Here Latexify is used automatically to display `odesys` in Latex within Markdown
-documents or notebook environments like Pluto.jl.)
+Catalyst models can be simulated using all three modes. Below we give a simple introduction to each, with [more detailed documentation provided elsewhere](@ref ref).
 
-Before we can solve the ODEs, we need to specify the values of the parameters in
-the model, the initial condition, and the time interval to solve the model on.
-To do this we need to build mappings from the symbolic parameters and the
-species to the corresponding numerical values for parameters and initial
-conditions. We can build such mappings in several ways. One is to use Julia
-`Symbols` to specify the values like
-```@example tut1
-pmap  = (:α => .5, :K => 40, :n => 2, :δ => log(2)/120,
-         :γ => 5e-3, :β => log(2)/6, :μ => log(2)/60)
-u₀map = [:m₁ => 0., :m₂ => 0., :m₃ => 0., :P₁ => 20., :P₂ => 0., :P₃ => 0.]
-nothing   # hide
+### [ODE simulations](@id introduction_to_catalyst_model_simulation_ode)
+Let us consider the infectious disease model declared in the previous section. To simulate it we need (in addition to the model):
+* The simulation's initial condition. That is, the concentration (or copy numbers) of each species at the start of the simulation.
+* The simulation's timespan. That is, the time frame over which we wish to run the simulation.
+* The simulation's parameter values. That is, the values of the model's parameters for this simulation.
+
+The initial conditions and parameter values are declared as vectors of pairs. Here, each entry pairs the symbol corresponding to a specific quantity's name to its value (e.g. `:S => 99` denotes that the initial concentration of `S` should be `99`). The timespan is simply a tuple with the simulations *starting* and *final* times. Once we have set these values, we collect them in a so-called `ODEProblem`.
+```@example intro_1
+u0 = [:S => 99, :I => 1, :R => 0]
+tspan = (0.0, 20.0)
+ps = [:β => 0.025, :γ => 0.2]
+oprob = ODEProblem(sir_model, u0, tspan, ps)
 ```
-Alternatively, we can use ModelingToolkit-based symbolic species variables to
-specify these mappings like
-```@example tut1
-t = default_t()
-@parameters  α K n δ γ β μ
-@species m₁(t) m₂(t) m₃(t) P₁(t) P₂(t) P₃(t)
-psymmap  = (α => .5, K => 40, n => 2, δ => log(2)/120,
-         γ => 5e-3, β => 20*log(2)/120, μ => log(2)/60)
-u₀symmap = [m₁ => 0., m₂ => 0., m₃ => 0., P₁ => 20., P₂ => 0., P₃ => 0.]
-nothing   # hide
+To simulate our `ODEProblem`, we simply input it to the `solve` command. Before we can do so, however, we also need to activate the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package (which is used for all ODE simulations in Julia).
+```@example intro_1
+using OrdinaryDiffEq
+sol = solve(oprob)
+```
+Finally, we can plot the solution using the [Plots.jl](https://github.com/JuliaPlots/Plots.jl) package's `plot` function:
+```@example intro_1
+using Plots
+plot(sol)
 ```
-Knowing these mappings we can set up the `ODEProblem` we want to solve:
-
-```@example tut1
-# time interval to solve on
-tspan = (0., 10000.)
 
-# create the ODEProblem we want to solve
-oprob = ODEProblem(repressilator, u₀map, tspan, pmap)
-nothing   # hide
+### [SDE simulations](@id introduction_to_catalyst_model_simulation_sde)
+SDE simulations can be carried out in a very similar manner to ODE ones, but by creating an `SDEProblem` instead of an `ODEProblem`. However, we can re-use the problem inputs from previously:
+```@example intro_1
+sprob = SDEProblem(sir_model, u0, tspan, ps)
 ```
-By passing `repressilator` directly to the `ODEProblem`, Catalyst has to
-(internally) call `convert(ODESystem, repressilator)` again to generate the
-symbolic ODEs. We could instead pass `odesys` directly like
-```@example tut1
-odesys = complete(odesys)
-oprob2 = ODEProblem(odesys, u₀symmap, tspan, psymmap)
-nothing   # hide
+Next, we import the [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl) package (which is used for all SDE simulations in Julia) and use it to simulate our SDE. 
+```@example intro_1
+using StochasticDiffEq
+sol = solve(sprob, STrapezoid())
+sol = solve(sprob, STrapezoid(); seed = 1234) # hide
+nothing # hide
 ```
-`oprob` and `oprob2` are functionally equivalent, each representing the same
-underlying problem.
-
-!!! note
-    When passing `odesys` to `ODEProblem` we needed to use the symbolic
-    variable-based parameter mappings, `u₀symmap` and `psymmap`, while when
-    directly passing `repressilator` we could use either those or the
-    `Symbol`-based mappings, `u₀map` and `pmap`. `Symbol`-based mappings can
-    always be converted to `symbolic` mappings using [`symmap_to_varmap`](@ref),
-    see the [Basic Syntax](@ref basic_examples) section for more details.
+Note that we here have to provide a second argument to the `solver` command. This is our choice of [simulation method](@ref ref). While we can provide this for ODE simulations as well, OrdinaryDiffEq.jl is able to [automatically select a suitable solver for any problem](@ref ref). This is currently not possible for SDEs. For now, `STrapezoid` is often a passable default choice, however, for important applications it can be good to [closer study the available SDE solvers](@ref ref).
 
-
-!!! note
-    Above we have used `repressilator = complete(repressilator)` and `odesys = complete(odesys)` to mark these systems as *complete*, indicating to Catalyst and ModelingToolkit that these models are finalized. This must be done before any system is given as input to a `convert` call or some problem type. `ReactionSystem` models created through the @reaction_network` DSL (which is introduced elsewhere, and primarily used throughout these documentation) are always marked as complete when generated. Hence `complete` does not need to be called on them. Symbolically generated `ReactionSystem`s, `ReactionSystem`s generated via the `@network_component` macro, and any ModelingToolkit system generated by `convert` always needs to be manually marked as `complete` as we do for `odesys` above. An expanded description on *completeness* can be found [here](@ref completeness_note).
-
-At this point we are all set to solve the ODEs. We can now use any ODE solver
-from within the
-[DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)
-package. We'll use the recommended default explicit solver, `Tsit5()`, and then
-plot the solutions:
-
-```@example tut1
-sol = solve(oprob, Tsit5(), saveat=10.)
+Next, we can again use `plot` to plot the solution.
+```@example intro_1
 plot(sol)
 ```
-We see the well-known oscillatory behavior of the repressilator! For more on the
-choices of ODE solvers, see the [DifferentialEquations.jl
-documentation](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/).
-
----
-
-## Stochastic simulation algorithms (SSAs) for stochastic chemical kinetics
-Let's now look at a stochastic chemical kinetics model of the repressilator,
-modeling it with jump processes. Here, we will construct a
-[JumpProcesses](https://docs.sciml.ai/JumpProcesses/stable/) `JumpProblem` that uses
-Gillespie's `Direct` method, and then solve it to generate one realization of
-the jump process:
-
-```@example tut1
-# redefine the initial condition to be integer valued
-u₀map = [:m₁ => 0, :m₂ => 0, :m₃ => 0, :P₁ => 20, :P₂ => 0, :P₃ => 0]
-
-# next we create a discrete problem to encode that our species are integer-valued:
-dprob = DiscreteProblem(repressilator, u₀map, tspan, pmap)
-
-# now, we create a JumpProblem, and specify Gillespie's Direct Method as the solver:
-jprob = JumpProblem(repressilator, dprob, Direct(), save_positions=(false,false))
-
-# now, let's solve and plot the jump process:
-sol = solve(jprob, SSAStepper(), saveat=10.)
+Since SDEs are *stochastic* (random), if we simulate it again we get another result:
+```@example intro_1
+sol = solve(sprob, ImplicitEM())
+sol = solve(sprob, STrapezoid(); seed = 5678) # hide
 plot(sol)
 ```
 
-We see that oscillations remain, but become much noisier. Note, in constructing
-the `JumpProblem` we could have used any of the SSAs that are part of JumpProcesses
-instead of the `Direct` method, see the list of SSAs (i.e., constant rate jump
-aggregators) in the
-[documentation](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#Jump-Aggregators-for-Exact-Simulation).
-
-Common questions that arise in using the JumpProcesses SSAs (i.e. Gillespie methods)
-are collated in the [JumpProcesses FAQ](https://docs.sciml.ai/JumpProcesses/stable/faq/).
-
----
-## Chemical Langevin equation (CLE) stochastic differential equation (SDE) models
-At an intermediate physical scale between macroscopic ODE models and microscopic
-stochastic chemical kinetics models lies the CLE, given by a system of SDEs that
-add to each ODE above a noise term. As the repressilator has species that get
-very close to zero in size, it is not a good candidate to model with the CLE
-(where solutions can then go negative and become unphysical). Let's create a
-simpler reaction network for a birth-death process that will stay non-negative:
-
-```@example tut1
-bdp = @reaction_network begin
-  c₁, X --> 2X
-  c₂, X --> 0
-  c₃, 0 --> X
-end
-p = (:c₁ => 1.0, :c₂ => 2.0, :c₃ => 50.)
-u₀ = [:X => 5.]
-tspan = (0.,4.)
-nothing   # hide
+### [Jump simulations](@id introduction_to_catalyst_model_simulation_jumps)
+Finally, Catalyst models can be simulated using jump simulations. These simulate the exact (and randomised) occurrence of system reactions, and their effect on the system's state. Jump simulations require the creation of `JumpProblem`. These are a bit different from `ODEProblem`s and `SDEProblem`s in that we first create a `DiscreteProblem` (to which we provide the model, initial condition, time span, and parameter values), and then use this (in addition to again providing our model) as input to `JumpProblem`.
+```@example intro_1
+dprob = DiscreteProblem(sir_model, u0, tspan, ps)
+jprob = JumpProblem(sir_model, dprob)
+nothing # hide
 ```
-
-The corresponding Chemical Langevin Equation SDE is then
-
-```math
-dX(t) = \left( c_1 X\left( t \right) - c_2 X\left( t \right) + c_3 \right) dt + \sqrt{c_1 X(t)} dW_1(t) - \sqrt{c_2 X(t)} dW_2(t) + \sqrt{c_3} dW_3(t)
+Finally, this problem can be simulated and plotted just like ODEs and SDEs (however, this time requiring the [JumpProcesses.jl](https://github.com/SciML/JumpProcesses.jl) simulation package):
+```@example intro_1
+using JumpProcesses
+sol = solve(jprob)
+sol = solve(jprob; seed = 12345) # hide
+plot(sol)
 ```
+Looking at the plot, we can actually distinguish the individual jumps of the simulation.
 
-where each $W_i(t)$ denotes an independent Brownian Motion. We can solve the CLE
-model by creating an `SDEProblem` and solving it similarly to what we did for ODEs
-above:
-
-```@example tut1
-# SDEProblem for CLE
-sprob = SDEProblem(bdp, u₀, tspan, p)
+### [Basic model and simulation querying](@id introduction_to_catalyst_querying)
+In this section, we will consider some basic functions for querying the content of Catalyst models. If we wish to list the species a model consists of, we can use the `species` function:
+```@example intro_1
+species(sir_model)
+```
+Next, to list the parameters we can use the `parameters` function:
+```@example intro_1
+parameters(sir_model)
+```
+!!! note
+ The `species` and `parameters` functions return a model's species/parameters as so-called *symbolic variables*. This is a special type (which can be used to form *symbolic expressions) that is described in more detail [here](@ref ref).
 
-# solve and plot, tstops is used to specify enough points
-# that the plot looks well-resolved
-sol = solve(sprob, LambaEM(), tstops = range(0., step = 4e-3, length = 1001))
-plot(sol)
+Finally, to list the reactions we can use
+```@example intro_1
+reactions(sir_model)
 ```
+A full list of similar and related functions can be found in [the API](@ref ref).
 
-We again have complete freedom to select any of the
-StochasticDiffEq.jl SDE solvers, see the
-[documentation](https://docs.sciml.ai/stable/modules/DiffEqDocs/solvers/sde_solve/).
+Simulation solutions can also be queried in various ways. To receive the value of species $P$ across the solution `sol` we simply call
+```@example intro_1
+sol[:P]
+```
 
----
-## Specifying a complete model via the DSL
-In the previous examples we specified initial conditions and parameter values
-via mappings that were constructed after building our [`ReactionSystem`](@ref).
-Catalyst also supports specifying default values for these during
-`ReactionSystem` construction. For example, for the last SDE example we
-could have also built and simulated the complete model using the DSL like
-```@example tut1
-bdp2 = @reaction_network begin
-    @parameters c₁ = 1.0 c₂ = 2.0 c₃ = 50.0
-    @species X(t) = 5.0
-    c₁, X --> 2X
-    c₂, X --> 0
-    c₃, 0 --> X
-end
-tspan = (0., 4.)
-sprob2 = SDEProblem(bdp2, [], tspan)
+Finally, it is possible to print a model in [LaTeX format](https://en.wikipedia.org/wiki/LaTeX) using the [Latexify.jl](https://github.com/korsbo/Latexify.jl) package. To print it formatted as chemical reactions we simply call `latexify` with our model as input:
+```@example intro_1
+using Latexify
+latexify(sir_model)
 ```
-Let's now simulate both models, starting from the same random number generator
-seed, and check we get the same solutions
-```@example tut1
-using Random
-Random.seed!(1)
-sol = solve(sprob, LambaEM(), tstops = range(0., step = 4e-3, length = 1001))
-p1 = plot(sol)
-Random.seed!(1)
-sol2 = solve(sprob2, LambaEM(), tstops = range(0., step = 4e-3, length = 1001))
-p2 = plot(sol2)
-plot(p1, p2, layout = (2,1))
+If we instead wish to print equations the model is converted to for ODE simulations, we have to add the `form = :ode` argument:
+```@example intro_1
+latexify(sir_model; form = :ode)
 ```
 
-For details on what information can be specified via the DSL see the [The
-Reaction DSL](@ref dsl_description) tutorial.
+!!! note
+    To have `latexify`'s output printed in LaTeX format you must use a programming environment that actually supports this (like a [JuPyteR](https://github.com/JuliaLang/IJulia.jl) or [Pluto](https://github.com/fonsp/Pluto.jl) notebook). Otherwise, you will simply print the LaTeX code which would generate the LaTeX print. 
 
----
-## Reaction rate laws used in simulations
-In generating mathematical models from a [`ReactionSystem`](@ref), reaction
-rates are treated as *microscopic* rates. That is, for a general mass action
-reaction of the form $n_1 S_1 + n_2 S_2 + \dots n_M S_M \to \dots$ with
-stoichiometric substrate coefficients $\{n_i\}_{i=1}^M$ and rate constant $k$,
-the corresponding ODE and SDE rate laws are taken to be
+## [Reaction rate laws used in ODE simulations](@id introduction_to_catalyst_rate_equations)
+When generating mathematical models from Catalyst-generated [`ReactionSystem`](@ref) models, reaction rates are treated as *microscopic* rates. That is, for a general mass action reaction of the form 
+```math
+n_1 S_1 + n_2 S_2 + \dots n_M S_M \to \dots
+```
+ with stoichiometric substrate coefficients $\{n_i\}_{i=1}^M$ and rate $k$, the corresponding ODE and SDE rate laws are taken to be
 ```math
 k \prod_{i=1}^M \frac{(S_i)^{n_i}}{n_i!},
 ```
-while the jump process transition rate (i.e., the propensity or intensity
-function) is
+while the jump process transition rate (i.e., the propensity or intensity function) is
 ```math
 k \prod_{i=1}^M \frac{S_i (S_i-1) \dots (S_i-n_i+1)}{n_i!}.
 ```
@@ -306,19 +169,12 @@ k \frac{X^2}{2!} \frac{Y^3}{3!} \\
 giving the ODE model
 ```math
 \begin{align*}
-\frac{dX}{dt} &=  -2 k \frac{X^2}{2!} \frac{Y^3}{3!}, &
-\frac{dY}{dt} &=  -3 k \frac{X^2}{2!} \frac{Y^3}{3!}, &
+\frac{dX}{dt} &= -2 k \frac{X^2}{2!} \frac{Y^3}{3!}, &
+\frac{dY}{dt} &= -3 k \frac{X^2}{2!} \frac{Y^3}{3!}, &
 \frac{dZ}{dt} &= k \frac{X^2}{2!} \frac{Y^3}{3!}.
 \end{align*}
 ```
-This implicit rescaling of rate constants can be disabled through explicit
-conversion of a [`ReactionSystem`](@ref) to another system via
-[`Base.convert`](@ref) using the `combinatoric_ratelaws=false` keyword
-argument, i.e.
-```julia
-rn = @reaction_network ...
-convert(ODESystem, rn; combinatoric_ratelaws=false)
-```
+This implicit rescaling of rate constants can be disabled through [explicit providing the `combinatoric_ratelaws=false` argument](@ref ref) when a set of equations is generated.
 
 For the previous example using this keyword argument would give the rate law
 ```math
@@ -327,32 +183,20 @@ k X^2 Y^3
 and the ODE model
 ```math
 \begin{align*}
-\frac{dX}{dt} &=  -2 k X^2 Y^3, &
-\frac{dY}{dt} &=  -3 k X^2 Y^3, &
+\frac{dX}{dt} &= -2 k X^2 Y^3, &
+\frac{dY}{dt} &= -3 k X^2 Y^3, &
 \frac{dZ}{dt} &= k X^2 Y^3.
 \end{align*}
 ```
 
----
-## Notes
-1. For each of the preceding models we converted the `ReactionSystem` to, i.e.,
-   ODEs, jumps, or SDEs, we had two paths for conversion:
-
-    a. Convert to the corresponding ModelingToolkit system and then use it in
-       creating the corresponding problem.
-
-    b. Directly create the desired problem type from the `ReactionSystem`.
+## [Next steps](@id introduction_to_catalyst_next_steps)
+The above tutorial gives enough background to use create and simulate basic CRN models using Catalyst. The remaining documentation goes through additional features of the package. In some places, it will also provide introductions to basic theory, and general advice on how to carry out various workflows. While it would be possible to read it from start to finish, you might also select specific sections depending on your intended use of Catalyst:
+- If you have not read it already, this documentation's [home](@ref ref) page provides some useful information.
+- The [basic](@ref ref) and [advanced](@ref ref) modelling tutorials describe a range of options for model creation, and are useful to everyone who plans on using Catalyst extensively.
+- The introduction to [simulations](@ref ref), [spatial modelling](@ref ref), and [inverse problem solving](@ref ref) are useful to anyone who plans on using Catalyst for any of these purposes.
+- The [Advanced introduction to Catalyst](@ref ref) is not required to utilise any Catalyst feature. While it is not intended to be read at an early stage, it is still useful to anyone who has, or plans to, use Catalyst extensively.
 
-   The latter is more convenient, however, the former will be more efficient if
-   one needs to repeatedly create the associated `Problem`.
-2. ModelingToolkit offers many options for optimizing the generated ODEs and
-   SDEs, including options to build functions for evaluating Jacobians and/or
-   multithreaded versions of derivative evaluation functions. See the options
-   for
-   [`ODEProblem`s](https://mtk.sciml.ai/dev/systems/ODESystem/#DiffEqBase.ODEProblem)
-   and
-   [`SDEProblem`s](https://mtk.sciml.ai/dev/systems/SDESystem/#DiffEqBase.SDEProblem).
 
 ---
-## References
+## [References](@id introduction_to_catalyst_references)
 [^1]: [Torkel E. Loman, Yingbo Ma, Vasily Ilin, Shashi Gowda, Niklas Korsbo, Nikhil Yewale, Chris Rackauckas, Samuel A. Isaacson, *Catalyst: Fast and flexible modeling of reaction networks*, PLOS Computational Biology (2023).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011530)
\ No newline at end of file

From 1b6ed62dad36131c8d116abf8cf11a4a148e4956 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sun, 28 Apr 2024 18:19:33 -0400
Subject: [PATCH 2/4] save progress (add relevant references)

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

diff --git a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
index 4261bd8ad7..8f3b97b3f6 100644
--- a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
+++ b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
@@ -135,6 +135,7 @@ Simulation solutions can also be queried in various ways. To receive the value o
 ```@example intro_1
 sol[:P]
 ```
+A more throughout tutorial on how to query solutions (and other relevant structures) for parameters and species values can be found [here](@ref simulation_structure_interfacing).
 
 Finally, it is possible to print a model in [LaTeX format](https://en.wikipedia.org/wiki/LaTeX) using the [Latexify.jl](https://github.com/korsbo/Latexify.jl) package. To print it formatted as chemical reactions we simply call `latexify` with our model as input:
 ```@example intro_1
@@ -145,6 +146,7 @@ If we instead wish to print equations the model is converted to for ODE simulati
 ```@example intro_1
 latexify(sir_model; form = :ode)
 ```
+Printing of models using Latexify is described in more detail [here](@ref ref).
 
 !!! note
     To have `latexify`'s output printed in LaTeX format you must use a programming environment that actually supports this (like a [JuPyteR](https://github.com/JuliaLang/IJulia.jl) or [Pluto](https://github.com/fonsp/Pluto.jl) notebook). Otherwise, you will simply print the LaTeX code which would generate the LaTeX print. 

From 7960c1e9efa069019877ce58b74022d8e70148fe Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Wed, 29 May 2024 18:46:18 -0400
Subject: [PATCH 3/4] exapnd symbolic example

---
 .../introduction_to_catalyst.md               | 79 +++++++++++++++----
 1 file changed, 63 insertions(+), 16 deletions(-)

diff --git a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
index 8f3b97b3f6..3c245a04d7 100644
--- a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
+++ b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
@@ -1,5 +1,5 @@
 # [Introduction to Catalyst](@id introduction_to_catalyst)
-This tutorial provides a basic introduction on how to create chemical reaction network (CRN) models using Catalyst, and how to perform ODE, SDE, and jump simulation using these. [An alternative introduction](@ref catalyst_for_new_julia_users) is available for users with little familiarity with Julia programming. 
+This tutorial provides a basic introduction on how to create chemical reaction network (CRN) models using Catalyst, and how to perform ODE, SDE, and jump simulation using these. [An alternative introduction](@ref catalyst_for_new_julia_users) is available for users who are unfamiliar with Julia programming. 
 
 Before starting this tutorial, we must (unless this has already been done) [install and activate](@ref ref) the Catalyst package:
 ```julia
@@ -9,7 +9,7 @@ using Catalyst
 ```
 
 ## [Creating basic Catalyst models](@id introduction_to_catalyst_model_creation)
-Catalyst exports the [`@reaction_network`](@ref) [macro]([macro](https://docs.julialang.org/en/v1/manual/metaprogramming/#man-macros)), which provides the main way of creating CRN models (alternative approaches are described [later](@ref ref)). This is followed by a `begin ... end` block which encapsulates the reactions of the system. E.g here we create a simple [birth-death model](@ref ref) (where a single species $X$ is both produced and degraded):
+Catalyst exports the [`@reaction_network`](@ref) [macro]([macro](https://docs.julialang.org/en/v1/manual/metaprogramming/#man-macros)), which provides the main way of creating CRN models (alternative approaches are described [later](@ref ref)). This is followed by a `begin ... end` block which encapsulates the model's reactions. E.g here we create a simple [birth-death model](@ref ref) (where a single species $X$ is both produced and degraded):
 ```@example intro_1
 using Catalyst # hide
 bd_model = @reaction_network begin
@@ -41,7 +41,7 @@ In all three cases, each reaction consists of:
 For reactions with multiple substrates/products (e.g. `S + E --> SE`), these are separated by `+`. Where substrates/products contain multiple units of the same species (e.g. `S + I --> 2I`), this is indicated by pre-appending that species with the corresponding [stoichiometric coefficient](https://en.wikipedia.org/wiki/Stoichiometry). The absences of substrates/products (that is, in pure production/degradation reactions like `0 --> X` and `X --> 0`) are denoted with `0`. A more throughout description of how to create models using `@reaction_network` (also called *the Catalyst DSL*) is provided [here](@ref ref).
 
 ## [Simulating Catalyst models](@id introduction_to_catalyst_model_simulation)
-There exist three primary modes of simulation CRN models:
+There exist three primary modes for simulating CRN models:
 1. Using ordinary differential equations (ODEs).
 2. Using stochastic differential equations (SDEs).
 3. Using jump process simulations.
@@ -54,7 +54,7 @@ Let us consider the infectious disease model declared in the previous section. T
 * The simulation's timespan. That is, the time frame over which we wish to run the simulation.
 * The simulation's parameter values. That is, the values of the model's parameters for this simulation.
 
-The initial conditions and parameter values are declared as vectors of pairs. Here, each entry pairs the symbol corresponding to a specific quantity's name to its value (e.g. `:S => 99` denotes that the initial concentration of `S` should be `99`). The timespan is simply a tuple with the simulations *starting* and *final* times. Once we have set these values, we collect them in a so-called `ODEProblem`.
+The initial conditions and parameter values are declared as vectors of pairs. Here, each entry pairs the symbol corresponding to a specific quantity's name to its value (e.g. `:S => 99` denotes that the initial value of `S` should be `99`). The timespan is simply a tuple with the simulations *starting* and *final* times. Once we have set these values, we collect them in a so-called `ODEProblem`.
 ```@example intro_1
 u0 = [:S => 99, :I => 1, :R => 0]
 tspan = (0.0, 20.0)
@@ -84,7 +84,8 @@ sol = solve(sprob, STrapezoid())
 sol = solve(sprob, STrapezoid(); seed = 1234) # hide
 nothing # hide
 ```
-Note that we here have to provide a second argument to the `solver` command. This is our choice of [simulation method](@ref ref). While we can provide this for ODE simulations as well, OrdinaryDiffEq.jl is able to [automatically select a suitable solver for any problem](@ref ref). This is currently not possible for SDEs. For now, `STrapezoid` is often a passable default choice, however, for important applications it can be good to [closer study the available SDE solvers](@ref ref).
+!!! note
+    Here, we have to provide a second argument to the `solve` command. This is our choice of [simulation method](@ref ref). While we can provide this for ODE simulations as well, OrdinaryDiffEq.jl is able to [automatically select a suitable solver for any problem](@ref ref). This is currently not possible for SDEs. For now, `STrapezoid` is often a passable default choice, however, for important applications it can be good to [closer study the available SDE solvers](@ref ref).
 
 Next, we can again use `plot` to plot the solution.
 ```@example intro_1
@@ -92,23 +93,23 @@ plot(sol)
 ```
 Since SDEs are *stochastic* (random), if we simulate it again we get another result:
 ```@example intro_1
-sol = solve(sprob, ImplicitEM())
+sol = solve(sprob, STrapezoid())
 sol = solve(sprob, STrapezoid(); seed = 5678) # hide
 plot(sol)
 ```
 
 ### [Jump simulations](@id introduction_to_catalyst_model_simulation_jumps)
-Finally, Catalyst models can be simulated using jump simulations. These simulate the exact (and randomised) occurrence of system reactions, and their effect on the system's state. Jump simulations require the creation of `JumpProblem`. These are a bit different from `ODEProblem`s and `SDEProblem`s in that we first create a `DiscreteProblem` (to which we provide the model, initial condition, time span, and parameter values), and then use this (in addition to again providing our model) as input to `JumpProblem`.
+Finally, Catalyst models can be simulated using jump simulations. These simulate the exact (and randomised) occurrence of system reactions, and their effect on the system's state. Jump simulations require the creation of `JumpProblem`. These are a bit different from `ODEProblem`s and `SDEProblem`s in that we first create a `DiscreteProblem` (to which we provide the model, initial condition, time span, and parameter values), and then use this (in addition to again providing our model) as input to `JumpProblem`. We most also load the [JumpProcesses.jl](https://github.com/SciML/JumpProcesses.jl) package.
 ```@example intro_1
+using JumpProcesses
 dprob = DiscreteProblem(sir_model, u0, tspan, ps)
-jprob = JumpProblem(sir_model, dprob)
+jprob = JumpProblem(sir_model, dprob, Direct())
 nothing # hide
 ```
-Finally, this problem can be simulated and plotted just like ODEs and SDEs (however, this time requiring the [JumpProcesses.jl](https://github.com/SciML/JumpProcesses.jl) simulation package):
+Finally, this problem can be simulated and plotted just like ODEs and SDEs:
 ```@example intro_1
-using JumpProcesses
-sol = solve(jprob)
-sol = solve(jprob; seed = 12345) # hide
+sol = solve(jprob, SSAStepper())
+sol = solve(jprob, SSAStepper(); seed = 12345) # hide
 plot(sol)
 ```
 Looking at the plot, we can actually distinguish the individual jumps of the simulation.
@@ -123,7 +124,7 @@ Next, to list the parameters we can use the `parameters` function:
 parameters(sir_model)
 ```
 !!! note
- The `species` and `parameters` functions return a model's species/parameters as so-called *symbolic variables*. This is a special type (which can be used to form *symbolic expressions) that is described in more detail [here](@ref ref).
+ The `species` and `parameters` functions return a model's species/parameters as so-called *symbolic variables*. This is a special type (which can be used to form *symbolic expressions*) that is described in more detail [here](@ref ref).
 
 Finally, to list the reactions we can use
 ```@example intro_1
@@ -131,9 +132,9 @@ reactions(sir_model)
 ```
 A full list of similar and related functions can be found in [the API](@ref ref).
 
-Simulation solutions can also be queried in various ways. To receive the value of species $P$ across the solution `sol` we simply call
+Simulation solutions can also be queried in various ways. To retrieve the value of species $I$ across the solution `sol` we simply call
 ```@example intro_1
-sol[:P]
+sol[:I]
 ```
 A more throughout tutorial on how to query solutions (and other relevant structures) for parameters and species values can be found [here](@ref simulation_structure_interfacing).
 
@@ -156,7 +157,7 @@ When generating mathematical models from Catalyst-generated [`ReactionSystem`](@
 ```math
 n_1 S_1 + n_2 S_2 + \dots n_M S_M \to \dots
 ```
- with stoichiometric substrate coefficients $\{n_i\}_{i=1}^M$ and rate $k$, the corresponding ODE and SDE rate laws are taken to be
+with stoichiometric substrate coefficients $\{n_i\}_{i=1}^M$ and rate $k$, the corresponding ODE and SDE rate laws are taken to be
 ```math
 k \prod_{i=1}^M \frac{(S_i)^{n_i}}{n_i!},
 ```
@@ -191,12 +192,58 @@ and the ODE model
 \end{align*}
 ```
 
+## [Catalyst represents its models symbolically, and these can be generated programmatically](@id introduction_to_catalyst_sym_prog)
+Above we have described he basics of model creation and simulation in Catalyst. Before we close this introduction, however, we would also like to describe that Catalyst represents its model equations symbolically (i.e. storing the actual algebraic expressions). Internally, Catalyst harnesses this in various ways (to e.g. boost simulation speed), however, an experienced user can also take advantage of it.
+
+Throughout most of Catalyst's documentation, models are declared using the `@reaction_network` macro. However, to gain better access to internal representations, it is typically better to create the model *programmatically*. Below we will create a model programmatically. In the context of this short tutorial, we will not go into all the details of what is going on (these can, however, be found [here](@ref ref)).
+
+A [Repressilator](@ref ref) is a system with three species ($X$, $Y$, and $Z$), each repressing the production of the subsequent species in the circuit. We begin with declaring the model species and parameters.
+```@example intro_prog_sym
+using Catalyst # hide
+t = default_t()
+@species X(t) Y(t) Z(t)
+@parameters K n d
+nothing # hide
+```
+Next, we create a production and a degradation reaction for each species. Normally, we'd have to list each reaction, but since we are building our model programmatically, we can utilise a `for`-loop. Finally, the model is created using the the `ReactionSystem` constructor:
+```@example intro_prog_sym
+rxs = []
+for (sp1, sp2) in zip([X, Y, Z], [Z, X, Y])
+    push!(rxs, Reaction(K/(K + sp2^3), [], [sp1]))
+    push!(rxs, Reaction(d, [sp1], []))
+end
+@named repressilator = ReactionSystem(rxs, t)
+repressilator = complete(repressilator)
+```
+We now have a model, `repressilator`, which can be analysed and simulated just like the models we have created previously:
+```@example intro_prog_sym
+using OrdinaryDiffEq, Plots # hide
+u0 = [X => 0.8, Y => 0.5, Z => 0.5]
+tspan = (0.0, 10.)
+ps = [K => 0.1, n => 3, d => 0.5]
+oprob = ODEProblem(repressilator, u0, tspan, ps)
+sol = solve(oprob, Tsit5())
+plot(sol)
+```
+
+If we consider the code above, we use the species and parameters (e.g. `K` and `X`) to designate values and create our reactions. This is possible because the `@species` and `@parameters` macros create these as [*symbolic variables*](@ref ref). This enables us to form algebraic expressions using them. E.g. here we create an expression describing $X$'s production rate:
+```@example intro_prog_sym
+X_prod = K/(K + Z^3)
+```
+Such expressions are useful as they can be used to designate e.g. what values to plot. E.g. here we plot $X$'s concentration and production rate across the simulation:
+```@example intro_prog_sym
+plot(sol; idxs = [X, X_prod])
+```
+
+While the above description was brief, we hope it provide some idea how to harness Catalyst's ability for symbolic and programmatic modelling. With it, we conclude our introduction with some advice of which part of the documentation to read next. 
+
 ## [Next steps](@id introduction_to_catalyst_next_steps)
 The above tutorial gives enough background to use create and simulate basic CRN models using Catalyst. The remaining documentation goes through additional features of the package. In some places, it will also provide introductions to basic theory, and general advice on how to carry out various workflows. While it would be possible to read it from start to finish, you might also select specific sections depending on your intended use of Catalyst:
 - If you have not read it already, this documentation's [home](@ref ref) page provides some useful information.
 - The [basic](@ref ref) and [advanced](@ref ref) modelling tutorials describe a range of options for model creation, and are useful to everyone who plans on using Catalyst extensively.
 - The introduction to [simulations](@ref ref), [spatial modelling](@ref ref), and [inverse problem solving](@ref ref) are useful to anyone who plans on using Catalyst for any of these purposes.
 - The [Advanced introduction to Catalyst](@ref ref) is not required to utilise any Catalyst feature. While it is not intended to be read at an early stage, it is still useful to anyone who has, or plans to, use Catalyst extensively.
+- If you are interested to learn more on how to create your models programmatically, an extensive description can be found [here](@ref ref).
 
 
 ---

From 6e730f93be94db261378bd3566baef6364fa0d31 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Thu, 30 May 2024 20:58:25 -0400
Subject: [PATCH 4/4] misc writing improvmenets

---
 .../introduction_to_catalyst.md               | 39 ++++++++++---------
 1 file changed, 20 insertions(+), 19 deletions(-)

diff --git a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
index 3c245a04d7..65ff8ab6c6 100644
--- a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
+++ b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
@@ -1,7 +1,7 @@
 # [Introduction to Catalyst](@id introduction_to_catalyst)
 This tutorial provides a basic introduction on how to create chemical reaction network (CRN) models using Catalyst, and how to perform ODE, SDE, and jump simulation using these. [An alternative introduction](@ref catalyst_for_new_julia_users) is available for users who are unfamiliar with Julia programming. 
 
-Before starting this tutorial, we must (unless this has already been done) [install and activate](@ref ref) the Catalyst package:
+Before starting this tutorial, we must (unless this has already been done) [install and activate](@ref catalyst_for_new_julia_users_packages) the Catalyst package:
 ```julia
 using Pkg
 Pkg.add("Catalyst")
@@ -9,20 +9,20 @@ using Catalyst
 ```
 
 ## [Creating basic Catalyst models](@id introduction_to_catalyst_model_creation)
-Catalyst exports the [`@reaction_network`](@ref) [macro]([macro](https://docs.julialang.org/en/v1/manual/metaprogramming/#man-macros)), which provides the main way of creating CRN models (alternative approaches are described [later](@ref ref)). This is followed by a `begin ... end` block which encapsulates the model's reactions. E.g here we create a simple [birth-death model](@ref ref) (where a single species $X$ is both produced and degraded):
+Catalyst exports the [`@reaction_network`](@ref) [macro]([macro](https://docs.julialang.org/en/v1/manual/metaprogramming/#man-macros)), which provides the main way of creating CRN models (an alternative, and in many ways more potent, approach is described [later](@ref programmatic_CRN_construction)). This is followed by a `begin ... end` block which encapsulates the model's reactions. E.g here we create a simple [birth-death model](@ref ref) (where a single species $X$ is both produced and degraded):
 ```@example intro_1
 using Catalyst # hide
 bd_model = @reaction_network begin
- p, 0 --> X
- d, X --> 0
+    p, 0 --> X
+    d, X --> 0
 end
 ```
 Next, we create a simple [Michaelis-Menten enzyme kinetics model](@ref ref) (where an enzyme, $E$ converts a substrate, $S$, to a product, $P$):
 ```@example intro_1
 mm_model = @reaction_network begin
- kB, S + E --> SE
- kD, SE --> S + E
- kP, SE --> P + E
+    kB, S + E --> SE
+    kD, SE --> S + E
+    kP, SE --> P + E
 end
 ```
 And finally, a [SIR model of an infectious disease](@ref ref) (where susceptible individuals, $I$, get infected by infected individuals, $I$, and later recovers, $R$):
@@ -38,7 +38,7 @@ In all three cases, each reaction consists of:
 - A set of substrates (species that are consumed by the reaction).
 - A set of products (species that are produced by the reaction).
 
-For reactions with multiple substrates/products (e.g. `S + E --> SE`), these are separated by `+`. Where substrates/products contain multiple units of the same species (e.g. `S + I --> 2I`), this is indicated by pre-appending that species with the corresponding [stoichiometric coefficient](https://en.wikipedia.org/wiki/Stoichiometry). The absences of substrates/products (that is, in pure production/degradation reactions like `0 --> X` and `X --> 0`) are denoted with `0`. A more throughout description of how to create models using `@reaction_network` (also called *the Catalyst DSL*) is provided [here](@ref ref).
+For reactions with multiple substrates/products (e.g. `S + E --> SE`), these are separated by `+`. Where substrates/products contain multiple units of the same species (e.g. `S + I --> 2I`), this is indicated by pre-appending that species with the corresponding [stoichiometric coefficient](https://en.wikipedia.org/wiki/Stoichiometry). The absences of substrates/products (that is, in pure production/degradation reactions like `0 --> X` and `X --> 0`) are denoted with `0`. A more throughout description of how to create models using `@reaction_network` (also called *the Catalyst DSL*) is provided [here](@ref dsl_description).
 
 ## [Simulating Catalyst models](@id introduction_to_catalyst_model_simulation)
 There exist three primary modes for simulating CRN models:
@@ -76,6 +76,7 @@ plot(sol)
 SDE simulations can be carried out in a very similar manner to ODE ones, but by creating an `SDEProblem` instead of an `ODEProblem`. However, we can re-use the problem inputs from previously:
 ```@example intro_1
 sprob = SDEProblem(sir_model, u0, tspan, ps)
+nothing # hide
 ```
 Next, we import the [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl) package (which is used for all SDE simulations in Julia) and use it to simulate our SDE. 
 ```@example intro_1
@@ -85,7 +86,7 @@ sol = solve(sprob, STrapezoid(); seed = 1234) # hide
 nothing # hide
 ```
 !!! note
-    Here, we have to provide a second argument to the `solve` command. This is our choice of [simulation method](@ref ref). While we can provide this for ODE simulations as well, OrdinaryDiffEq.jl is able to [automatically select a suitable solver for any problem](@ref ref). This is currently not possible for SDEs. For now, `STrapezoid` is often a passable default choice, however, for important applications it can be good to [closer study the available SDE solvers](@ref ref).
+    Here, we have to provide a second argument to the `solve` command. This is our choice of [simulation method](@ref ref). While we can provide this for ODE simulations as well, OrdinaryDiffEq.jl is able to automatically select a suitable solver for any problem. This is currently not possible for SDEs. For now, `STrapezoid` is often a passable default choice, however, for important applications it can be good to [closer study the available SDE solvers](@ref ref).
 
 Next, we can again use `plot` to plot the solution.
 ```@example intro_1
@@ -124,7 +125,7 @@ Next, to list the parameters we can use the `parameters` function:
 parameters(sir_model)
 ```
 !!! note
- The `species` and `parameters` functions return a model's species/parameters as so-called *symbolic variables*. This is a special type (which can be used to form *symbolic expressions*) that is described in more detail [here](@ref ref).
+    The `species` and `parameters` functions return a model's species/parameters as so-called *symbolic variables*. This is a special type (which can be used to [form *symbolic expressions*](@ref ref)) that is described in more detail [here](@ref ref).
 
 Finally, to list the reactions we can use
 ```@example intro_1
@@ -136,7 +137,7 @@ Simulation solutions can also be queried in various ways. To retrieve the value
 ```@example intro_1
 sol[:I]
 ```
-A more throughout tutorial on how to query solutions (and other relevant structures) for parameters and species values can be found [here](@ref simulation_structure_interfacing).
+A more throughout tutorial on how to query solutions (and other relevant structures) for species and parameters values can be found [here](@ref simulation_structure_interfacing).
 
 Finally, it is possible to print a model in [LaTeX format](https://en.wikipedia.org/wiki/LaTeX) using the [Latexify.jl](https://github.com/korsbo/Latexify.jl) package. To print it formatted as chemical reactions we simply call `latexify` with our model as input:
 ```@example intro_1
@@ -177,9 +178,9 @@ giving the ODE model
 \frac{dZ}{dt} &= k \frac{X^2}{2!} \frac{Y^3}{3!}.
 \end{align*}
 ```
-This implicit rescaling of rate constants can be disabled through [explicit providing the `combinatoric_ratelaws=false` argument](@ref ref) when a set of equations is generated.
+This implicit rescaling of rate constants can be disabled through [explicit providing the `combinatoric_ratelaws = false` argument](@ref ref) when an ODE is generated (e.g. through `ODEProblem`).
 
-For the previous example using this keyword argument would give the rate law
+For the previous example, using the `combinatoric_ratelaws = false` keyword argument would give the rate law
 ```math
 k X^2 Y^3
 ```
@@ -193,9 +194,9 @@ and the ODE model
 ```
 
 ## [Catalyst represents its models symbolically, and these can be generated programmatically](@id introduction_to_catalyst_sym_prog)
-Above we have described he basics of model creation and simulation in Catalyst. Before we close this introduction, however, we would also like to describe that Catalyst represents its model equations symbolically (i.e. storing the actual algebraic expressions). Internally, Catalyst harnesses this in various ways (to e.g. boost simulation speed), however, an experienced user can also take advantage of it.
+Above we have described the basics of model creation and simulation in Catalyst. Before we close this introduction, however, we would also like to describe that Catalyst represents its model equations symbolically (i.e. storing the actual algebraic expressions). Internally, Catalyst harnesses this in various ways (to e.g. boost simulation speed), however, an experienced user can also take direct advantage of it.
 
-Throughout most of Catalyst's documentation, models are declared using the `@reaction_network` macro. However, to gain better access to internal representations, it is typically better to create the model *programmatically*. Below we will create a model programmatically. In the context of this short tutorial, we will not go into all the details of what is going on (these can, however, be found [here](@ref ref)).
+Throughout most of Catalyst's documentation, models are declared using the `@reaction_network` macro. However, to gain better access to internal representations, it is typically better to create the model *programmatically*. Below we will create a model programmatically. In the context of this short tutorial, we will not go into all the details of what is going on (these can, however, be found [here](@ref programmatic_CRN_construction)).
 
 A [Repressilator](@ref ref) is a system with three species ($X$, $Y$, and $Z$), each repressing the production of the subsequent species in the circuit. We begin with declaring the model species and parameters.
 ```@example intro_prog_sym
@@ -205,7 +206,7 @@ t = default_t()
 @parameters K n d
 nothing # hide
 ```
-Next, we create a production and a degradation reaction for each species. Normally, we'd have to list each reaction, but since we are building our model programmatically, we can utilise a `for`-loop. Finally, the model is created using the the `ReactionSystem` constructor:
+Next, we create a production and a degradation reaction for each species. Normally, we have to list each reaction. However, since we are building our model programmatically, we can utilise a `for`-loop. Finally, the model is created using the the `ReactionSystem` constructor:
 ```@example intro_prog_sym
 rxs = []
 for (sp1, sp2) in zip([X, Y, Z], [Z, X, Y])
@@ -230,12 +231,12 @@ If we consider the code above, we use the species and parameters (e.g. `K` and `
 ```@example intro_prog_sym
 X_prod = K/(K + Z^3)
 ```
-Such expressions are useful as they can be used to designate e.g. what values to plot. E.g. here we plot $X$'s concentration and production rate across the simulation:
+Such expressions are useful as they e.g. can designate what values to plot. Here we use this to plot $X$'s concentration and production rate across the simulation:
 ```@example intro_prog_sym
 plot(sol; idxs = [X, X_prod])
 ```
 
-While the above description was brief, we hope it provide some idea how to harness Catalyst's ability for symbolic and programmatic modelling. With it, we conclude our introduction with some advice of which part of the documentation to read next. 
+While the above description was brief, we hope it provide some idea of how to harness Catalyst's ability for symbolic and programmatic modelling. With it, we conclude our introduction with some advice of which part of the documentation to read next. 
 
 ## [Next steps](@id introduction_to_catalyst_next_steps)
 The above tutorial gives enough background to use create and simulate basic CRN models using Catalyst. The remaining documentation goes through additional features of the package. In some places, it will also provide introductions to basic theory, and general advice on how to carry out various workflows. While it would be possible to read it from start to finish, you might also select specific sections depending on your intended use of Catalyst:
@@ -243,7 +244,7 @@ The above tutorial gives enough background to use create and simulate basic CRN
 - The [basic](@ref ref) and [advanced](@ref ref) modelling tutorials describe a range of options for model creation, and are useful to everyone who plans on using Catalyst extensively.
 - The introduction to [simulations](@ref ref), [spatial modelling](@ref ref), and [inverse problem solving](@ref ref) are useful to anyone who plans on using Catalyst for any of these purposes.
 - The [Advanced introduction to Catalyst](@ref ref) is not required to utilise any Catalyst feature. While it is not intended to be read at an early stage, it is still useful to anyone who has, or plans to, use Catalyst extensively.
-- If you are interested to learn more on how to create your models programmatically, an extensive description can be found [here](@ref ref).
+- If you are interested to learn more on how to create your models programmatically, an extensive description can be found [here](@ref programmatic_CRN_construction).
 
 
 ---