diff --git a/docs/pages.jl b/docs/pages.jl
index dbb52df4e6..4320e7162d 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -3,6 +3,7 @@ pages = Any[
     "Introduction to Catalyst" => Any[
         "introduction_to_catalyst/catalyst_for_new_julia_users.md",
         "introduction_to_catalyst/introduction_to_catalyst.md"
+        "introduction_to_catalyst/math_models_intro.md"
     ],
     "Model Creation and Properties" => Any[
         "model_creation/dsl_basics.md",
diff --git a/docs/src/introduction_to_catalyst/math_models_intro.md b/docs/src/introduction_to_catalyst/math_models_intro.md
new file mode 100644
index 0000000000..afcd82e7b7
--- /dev/null
+++ b/docs/src/introduction_to_catalyst/math_models_intro.md
@@ -0,0 +1,159 @@
+# [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)
+stochastic differential equation (SDE) models, and stochastic chemical kinetics
+(jump process) models. For each we show the abstract representations for the
+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.
+
+## 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}.
+```
+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.
+
+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))$.
+
+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
+a_k(\mathbf{X}(t)) = k \prod_{m=1}^M \frac{(X_m(t))^{\sigma_m^k}}{\sigma_m^k!},
+```
+for RRE ODE and CLE SDE models, and the propensity function
+```math
+a_k(\mathbf{X}(t)) = k \prod_{m=1}^M \frac{X_m(t) (X_m(t)-1) \dots (X_m(t)-\sigma_m^k+1)}{\sigma_m^k!},
+```
+for stochastic chemical kinetics jump process models.
+
+For example, for the reaction $2A + B \overset{k}{\to} 3 C$ we would have
+```math
+\mathbf{X}(t) = (A(t), B(t), C(t))
+```
+with $\sigma_1 = 2$, $\sigma_2 = 1$, $\sigma_3 = 0$, $\beta_1 = 0$, $\beta_2 =
+0$, $\beta_3 = 3$, $\nu_1 = -2$, $\nu_2 = -1$, and $\nu_3 = 3$. For an ODE/SDE
+model we would have the rate law
+```math
+a(\mathbf{X}(t)) = \frac{k}{2} A^2 B
+```
+while for a jump process model we would have the propensity function
+```math
+a(\mathbf{X}(t)) = \frac{k}{2} A (A-1) B.
+```
+
+Note, if the combinatoric factors are already included in one's rate constants,
+the implicit rescaling of rate constants can be disabled through use of the
+`combinatoric_ratelaws = false` argument to [`Base.convert`](@ref) or whatever
+Problem is being generated, i.e.
+```julia
+rn = @reaction_network ...
+osys = convert(ODESystem, rn; combinatoric_ratelaws = false)
+oprob = ODEProblem(osys, ...)
+
+# or
+oprob = ODEProblem(rn, ...; combinatoric_ratelaws = false)
+```
+In this case our ODE/SDE rate law would be
+```math
+a(\mathbf{X}(t)) = k A^2 B
+```
+while the jump process propensity function is
+```math
+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
+```
+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).
+
+For example, see the generated ODEs for the following network
+```@example math_examples
+using Catalyst, ModelingToolkit, Latexify
+rn = @reaction_network begin
+    k₁, 2A + B --> 3C
+    k₂, A --> 0
+    k₃, 0 --> A
+end
+osys = convert(ODESystem, rn)
+```
+Likewise, the following drops the combinatoric scaling factors, giving unscaled ODEs
+```@example math_examples
+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,
+```
+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,
+```julia
+rn = @reaction_network begin
+    k₁, 2A + B --> 3C
+    k₂, A --> 0
+    k₃, 0 --> A
+end
+```
+we would obtain the CLE SDEs
+```math
+\begin{align}
+dA(t) &= \left(- k_1 A^{2} B - k_2 A + k_3 \right) dt
+ - 2 \sqrt{\tfrac{k_1}{2} A^{2} B} \, dW_1(t)  - \sqrt{k_2 A} \, dW_2(t) + \sqrt{k_3} \, dW_3(t)
+\\
+dB(t) &= - \frac{k_1}{2} A^{2} B \, dt - \sqrt{\frac{k_1}{2} A^{2} B} \, dW_1(t) \\
+dC(t) &= \frac{3}{2} k_1 A^{2} B \, dt + 3 \sqrt{\frac{k_1}{2} A^{2} B} \, dW_1(t).
+\end{align}
+```
+
+## 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.
+```
+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].
+```
+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,
+```julia
+rn = @reaction_network begin
+    k₁, 2A + B --> 3C
+    k₂, A --> 0
+    k₃, 0 --> A
+end
+```
+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) \\
+B(t) &= B(0) - Y_1\left( \frac{k_1}{2} \int_0^t A(s^-)(A(s^-)-1) B(s^-) \, ds \right) \\
+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
+```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