diff --git a/HISTORY.md b/HISTORY.md
index ff92e9b5f5..739578bd84 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -39,9 +39,8 @@ briefly summarised in the following bullet points:
   new merged and/or composed `ReactionSystem`s from multiple component systems.
 
 #### General changes
-- The `default_t()` and `default_time_deriv()` functions are now the preferred
-  approaches for creating the default time independent variable and its
-  differential. i.e.
+- `default_t()` and `default_time_deriv()` functions should be used for creating
+  the default time independent variable and its differential. i.e.
   ```julia
   # do
   t = default_t()
diff --git a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
index d46d5a9397..cd04731536 100644
--- a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
+++ b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
@@ -3,7 +3,11 @@ 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 [1]. 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.
+Catalyst for ODE, SDE and jump process models. The [Mathematical Models Catalyst
+can Generate](@ref math_models_in_catalyst) documentation illustrates the
+abstract mathematical models Catalyst reaction models can be converted to,
+but please note it assumes one has already read this tutorial as a
+prerequisite.
 
 We begin by installing Catalyst and any needed packages into a new environment.
 This step can be skipped if you have already installed them in your current,
@@ -350,6 +354,8 @@ and the ODE model
 \end{align*}
 ```
 
+A more detailed summary of the precise mathematical equations Catalyst can generate is available in the [Mathematical Models Catalyst can Generate](@ref math_models_in_catalyst) documentation.
+
 ---
 ## Notes
 1. For each of the preceding models we converted the `ReactionSystem` to, i.e.,
diff --git a/docs/src/introduction_to_catalyst/math_models_intro.md b/docs/src/introduction_to_catalyst/math_models_intro.md
index afcd82e7b7..1357178f56 100644
--- a/docs/src/introduction_to_catalyst/math_models_intro.md
+++ b/docs/src/introduction_to_catalyst/math_models_intro.md
@@ -1,4 +1,4 @@
-# [Mathematical Models Catalyst Can Generate](@id math_models_in_catalyst)
+# [Mathematical Models Catalyst can Generate](@id math_models_in_catalyst)
 We now describe the types of mathematical models that Catalyst can generate from
 chemical reaction networks (CRNs), corresponding to reaction rate equation (RRE)
 ordinary differential equation (ODE) models, Chemical Langevin equation (CLE)
@@ -7,13 +7,16 @@ stochastic differential equation (SDE) models, and stochastic chemical kinetics
 models that Catalyst can support, along with concrete examples. Note that we
 restrict ourselves to models involving only chemical reactions, and do not
 consider more general models that Catalyst can support such as coupling in
-non-reaction ODEs, algebraic equations, or events. Please see [here]() for more
-details on how Catalyst supports such functionality.
+non-reaction ODEs, algebraic equations, or events. Please see the broader
+documentation for more details on how Catalyst supports such functionality.
+
+!!! note
+    This documentation assumes you have already read the [Introduction to Catalyst](@ref introduction_to_catalyst) tutorial.
 
 ## General Chemical Reaction Notation
 Suppose we have a reaction network with ``K`` reactions and ``M`` species, with the species labeled by $S_1$, $S_2$, $\dots$, $S_M$. We denote by
 ```math
-\mathbf{X}(t) = \begin{pmatrix} X_1(t) \\ \vdots \\ X_M(t)) \end{pmatrix}.
+\mathbf{X}(t) = \begin{pmatrix} X_1(t) \\ \vdots \\ X_M(t)) \end{pmatrix}
 ```
 the state vector for the amount of each species, i.e. $X_m(t)$ represents the amount of species $S_m$ at time $t$. This could be either a concentration or a count (i.e. "number of molecules" units), but for consistency between modeling representations we will assume the latter in the remainder of this introduction.
 
@@ -21,7 +24,7 @@ The $k$th chemical reaction is given by
 ```math
 \alpha_1^k S_1 + \alpha_2^k S_2 + \dots \alpha_M^k S_M \to \beta_1^k S_1 + \beta_2^k S_2 + \dots \beta_M^k S_M
 ```
-with $\alpha^k = (\alpha_1^k,\dots,\alpha_M^k)$ its substrate stoichiometry vector, $\beta^k = (\beta_1^k,\dots,\beta_M^k)$ its product stoichiometry vector, and $\nu^k = \beta^k - \alpha^k$ its net stoichiometry vector. $\nu^k$ corresponds to the change in $\mathbf{X}(t)$ when reaction $k$ occurs, i.e. $\mathbf{X}(t) \to \mathbf{X}(t) + \nu^k$. Along with the stoichiometry vectors, we assume each reaction has a reaction rate law (ODEs/SDEs) or propensity (jump process) function, $a_k(\mathbf{X}(t))$.
+with $\alpha^k = (\alpha_1^k,\dots,\alpha_M^k)$ its substrate stoichiometry vector, $\beta^k = (\beta_1^k,\dots,\beta_M^k)$ its product stoichiometry vector, and $\nu^k = \beta^k - \alpha^k$ its net stoichiometry vector. $\nu^k$ corresponds to the change in $\mathbf{X}(t)$ when reaction $k$ occurs, i.e. $\mathbf{X}(t) \to \mathbf{X}(t) + \nu^k$. Along with the stoichiometry vectors, we assume each reaction has a reaction rate law (ODEs/SDEs) or propensity (jump process) function, $a_k(\mathbf{X}(t),t)$.
 
 As explained in [the Catalyst introduction](@ref introduction_to_catalyst), for a mass action reaction where the preceding reaction has a fixed rate constant, $k$, this function would be the rate law
 ```math
@@ -33,7 +36,8 @@ a_k(\mathbf{X}(t)) = k \prod_{m=1}^M \frac{X_m(t) (X_m(t)-1) \dots (X_m(t)-\sigm
 ```
 for stochastic chemical kinetics jump process models.
 
-For example, for the reaction $2A + B \overset{k}{\to} 3 C$ we would have
+### Rate Law vs. Propensity Example:
+For the reaction $2A + B \overset{k}{\to} 3 C$ we would have
 ```math
 \mathbf{X}(t) = (A(t), B(t), C(t))
 ```
@@ -69,17 +73,19 @@ while the jump process propensity function is
 a(\mathbf{X}(t)) = k A (A-1) B.
 ```
 
-### Encoding General Rate Laws in Catalyst
-TODO, show how to factor for optimal representations.
-
 ## Reaction Rate Equation (RRE) ODE Models
 The RRE ODE models Catalyst creates for a general system correspond to the coupled system of ODEs given by
 ```math
-\frac{d X_m}{dt} =\sum_{k=1}^K \nu_m^k a_k(\mathbf{X}(t)), \quad m = 1,\dots,M
+\frac{d X_m}{dt} =\sum_{k=1}^K \nu_m^k a_k(\mathbf{X}(t),t), \quad m = 1,\dots,M.
 ```
-These models can be generated by creating `ODEProblem`s from Catalyst `ReactionSystem`s, and solved using the solvers in [OrdinaryDiffEq.jl](/~https://github.com/SciML/OrdinaryDiffEq.jl).
+These models can be generated by creating `ODEProblem`s from Catalyst `ReactionSystem`s, and solved using the solvers in [OrdinaryDiffEq.jl](/~https://github.com/SciML/OrdinaryDiffEq.jl). Similarly, creating `NonlinearProblem`s or `SteadyStateProblem`s will generate the coupled algebraic system of steady-state equations associated with a RRE ODE model, i.e.
+```math
+0 =\sum_{k=1}^K \nu_m^k a_k(\bar{\mathbf{X}}), \quad m = 1,\dots,M
+```
+for a steady-state $\bar{\mathbf{X}}$. Note, here we have assumed the rate laws are [autonomous](https://en.wikipedia.org/wiki/Autonomous_system_(mathematics)) so that the equations are well-defined.
 
-For example, see the generated ODEs for the following network
+### RRE ODE Example
+Let's see the generated ODEs for the following network
 ```@example math_examples
 using Catalyst, ModelingToolkit, Latexify
 rn = @reaction_network begin
@@ -97,11 +103,12 @@ osys = convert(ODESystem, rn; combinatoric_ratelaws = false)
 ## Chemical Langevin Equation (CLE) SDE Models
 The CLE SDE models Catalyst creates for a general system correspond to the coupled system of SDEs given by
 ```math
-d X_m = \sum_{k=1}^K \nu_m^k a_k(\mathbf{X}(t)) dt + \sum_{k=1}^K \nu_m^k \sqrt{a_k(\mathbf{X}(t))} dW_k(t), \quad m = 1,\dots,M,
+d X_m = \sum_{k=1}^K \nu_m^k a_k(\mathbf{X}(t),t) dt + \sum_{k=1}^K \nu_m^k \sqrt{a_k(\mathbf{X}(t),t)} dW_k(t), \quad m = 1,\dots,M,
 ```
 where each $W_k(t)$ represents an independent, standard Brownian Motion. Realizations of these processes can be generated by creating `SDEProblem`s from Catalyst `ReactionSystem`s, and sampling the processes using the solvers in [StochasticDiffEq.jl](/~https://github.com/SciML/StochasticDiffEq.jl).
 
-For example, for the same network as above,
+### CLE SDE Example
+Consider the same network as above,
 ```julia
 rn = @reaction_network begin
     k₁, 2A + B --> 3C
@@ -109,7 +116,7 @@ rn = @reaction_network begin
     k₃, 0 --> A
 end
 ```
-we would obtain the CLE SDEs
+We obtain the CLE SDEs
 ```math
 \begin{align}
 dA(t) &= \left(- k_1 A^{2} B - k_2 A + k_3 \right) dt
@@ -123,17 +130,18 @@ dC(t) &= \frac{3}{2} k_1 A^{2} B \, dt + 3 \sqrt{\frac{k_1}{2} A^{2} B} \, dW_1(
 ## Stochastic Chemical Kinetics Jump Process Models
 The stochastic chemical kinetics jump process models Catalyst creates for a general system correspond to the coupled system of jump processes, in the time change representation, given by
 ```math
-X_m(t) = X_m(0) + \sum_{k=1}^k \nu_m^k Y_k\left( \int_{0}^t a_k(\mathbf{X}(s^-) \, ds \right), \quad m = 1,\dots,M.
+X_m(t) = X_m(0) + \sum_{k=1}^K \nu_m^k Y_k\left( \int_{0}^t a_k(\mathbf{X}(s^-),s) \, ds \right), \quad m = 1,\dots,M.
 ```
 Here each $Y_k(t)$ denotes an independent unit rate Poisson counting process with $Y_k(0) = 0$, which counts the number of times the $k$th reaction has occurred up to time $t$. Realizations of these processes can be generated by creating `JumpProblem`s from Catalyst `ReactionSystem`s, and sampling the processes using the solvers in [JumpProcesses.jl](/~https://github.com/SciML/JumpProcesses.jl).
 
 Let $P(\mathbf{x},t) = \operatorname{Prob}[\mathbf{X}(t) = \mathbf{x}]$ represent the probability the state of the system, $\mathbf{X}(t)$, has the concrete value $\mathbf{x}$ at time $t$. The forward equation, i.e. Chemical Master Equation (CME), associated with $\mathbf{X}(t)$ is then the coupled system of ODEs over all possible values for $\mathbf{x}$ given by
 ```math
-\frac{dP}{dt}(\mathbf{x},t) = \sum_{k=1}^k \left[ a_k(\mathbf{x} - \nu^k) P(\mathbf{x} - \nu^k,t) - a_k(\mathbf{x}) P(\mathbf{x},t) \right].
+\frac{dP}{dt}(\mathbf{x},t) = \sum_{k=1}^k \left[ a_k(\mathbf{x} - \nu^k,t) P(\mathbf{x} - \nu^k,t) - a_k(\mathbf{x},t) P(\mathbf{x},t) \right].
 ```
 While Catalyst does not currently support generating and solving for $P(\mathbf{x},t)$, for sufficiently small models the [FiniteStateProjection.jl](/~https://github.com/SciML/FiniteStateProjection.jl) package can be used to generate and solve such models directly from Catalyst [`ReactionSystem`](@ref)s.
 
-For example, for the same network as above,
+### Stochastic Chemical Kinetics Jump Process Example
+Consider the same network as above,
 ```julia
 rn = @reaction_network begin
     k₁, 2A + B --> 3C
@@ -141,7 +149,7 @@ rn = @reaction_network begin
     k₃, 0 --> A
 end
 ```
-the time change process representation would be
+The time change process representation would be
 ```math
 \begin{align*}
 A(t) &= A(0) - 2 Y_1\left( \frac{k_1}{2} \int_0^t A(s^-)(A(s^-)-1) B(s^-) \, ds \right) - Y_2 \left( k_2 \int_0^t A(s^-) \, ds \right) + Y_3 \left( k_3 t \right) \\
@@ -149,11 +157,12 @@ B(t) &= B(0) - Y_1\left( \frac{k_1}{2} \int_0^t A(s^-)(A(s^-)-1) B(s^-) \, ds \r
 C(t) &= C(0) + 3 Y_1\left( \frac{k_1}{2} \int_0^t A(s^-)(A(s^-)-1) B(s^-) \, ds \right),
 \end{align*}
 ```
-while the CME would be
+while the CME would be the coupled (infinite) system of ODEs over all realizable values of the non-negative integers $a$, $b$, and $c$ given by
 ```math
 \begin{align*}
 \frac{dP}{dt}(a,b,c,t) &= \left[\tfrac{k_1}{2} (a+2) (a+1) (b+1) P(a+2,b+1,c-3,t) - \tfrac{k_1}{2} a (a-1) b P(a,b,c,t)\right] \\
 &\phantom{=} + \left[k_2 (a+1) P(a+1,b,c,t) - k_2 a P(a,b,c,t)\right] \\
 &\phantom{=} + \left[k_3 P(a-1,b,c,t) - k_3 P(a,b,c,t)\right].
 \end{align*}
-```
\ No newline at end of file
+```
+If we initially have $A(0) = a_0$, $B(0) = b_0$, and $C(0) = c_0$ then we would have one ODE for each of possible state $(a,b,c)$ where $a \in \{0,1,\dots\}$ (i.e. $a$ can be any non-negative integer), $b \in \{0,1,\dots,b_0\}$, and $c \in \{c_0, c_0 + 1,\dots, c_0 + 3 b_0\}$. Other initial conditions would lead to different possible ranges for $a$, $b$, and $c$.
\ No newline at end of file