From 626b73758725f55fecc7079fc091e8802b09ea8e Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Mon, 2 Dec 2024 13:57:47 -0500
Subject: [PATCH 01/24] init

---
 .../crnt.md}                                  | 229 +-----------
 docs/src/network_analysis/odes.md             | 346 ++++++++++++++++++
 2 files changed, 360 insertions(+), 215 deletions(-)
 rename docs/src/{model_creation/network_analysis.md => network_analysis/crnt.md} (60%)
 create mode 100644 docs/src/network_analysis/odes.md

diff --git a/docs/src/model_creation/network_analysis.md b/docs/src/network_analysis/crnt.md
similarity index 60%
rename from docs/src/model_creation/network_analysis.md
rename to docs/src/network_analysis/crnt.md
index 3f77a99f55..9a445da2a0 100644
--- a/docs/src/model_creation/network_analysis.md
+++ b/docs/src/network_analysis/crnt.md
@@ -1,214 +1,4 @@
-# [Network Analysis in Catalyst](@id network_analysis)
-
-In this tutorial we introduce several of the Catalyst API functions for network
-analysis. A complete summary of the exported functions is given in the API
-section
-[`Network-Analysis-and-Representations`](https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Network-Analysis-and-Representations).
-
-Note, currently API functions for network analysis and conservation law analysis
-do not work with constant species (currently only generated by SBMLToolkit).
-
-## [Network representation of the Repressilator `ReactionSystem`](@id network_analysis_repressilator_representation)
-We first load Catalyst and construct our model of the repressilator
-```@example s1
-using Catalyst
-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
-```
-In the [Introduction to Catalyst](@ref introduction_to_catalyst)
-tutorial we showed how the above network could be visualized as a
-species-reaction graph. There, species are represented by the nodes of the graph
-and edges show the reactions in which a given species is a substrate or product.
-```julia
-g = Graph(repressilator)
-```
-![Repressilator solution](../assets/repressilator.svg)
-
-We also showed in the [Introduction to Catalyst](@ref introduction_to_catalyst) tutorial that
-the reaction rate equation ODE model for the repressilator is
-```math
-\begin{aligned}
-\frac{dm_1(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_3}\left( t \right) \right)^{n}} - \delta {m_1}\left( t \right) + \gamma \\
-\frac{dm_2(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_1}\left( t \right) \right)^{n}} - \delta {m_2}\left( t \right) + \gamma \\
-\frac{dm_3(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_2}\left( t \right) \right)^{n}} - \delta {m_3}\left( t \right) + \gamma \\
-\frac{dP_1(t)}{dt} =& \beta {m_1}\left( t \right) - \mu {P_1}\left( t \right) \\
-\frac{dP_2(t)}{dt} =& \beta {m_2}\left( t \right) - \mu {P_2}\left( t \right) \\
-\frac{dP_3(t)}{dt} =& \beta {m_3}\left( t \right) - \mu {P_3}\left( t \right)
-\end{aligned}
-```
-
-## [Matrix-vector reaction rate equation representation](@id network_analysis_matrix_vector_representation)
-In general, reaction rate equation (RRE) ODE models for chemical reaction networks can
-be represented as a first-order system of ODEs in a compact matrix-vector notation. Suppose
-we have a reaction network with ``K`` reactions and ``M`` species, labelled by the state vector
-```math
-\mathbf{x}(t) = \begin{pmatrix} x_1(t) \\ \vdots \\ x_M(t)) \end{pmatrix}.
-```
-For the repressilator, ``\mathbf{x}(t)`` is just
-```@example s1
-x = species(repressilator)
-```
-The RRE ODEs satisfied by $\mathbf{x}(t)$ are then
-```math
-\frac{d\mathbf{x}}{dt} = N \mathbf{v}(\mathbf{x}(t),t),
-```
-where ``N`` is a constant ``M`` by ``K`` matrix with ``N_{m k}`` the net
-stoichiometric coefficient of species ``m`` in reaction ``k``.
-``\mathbf{v}(\mathbf{x}(t),t)`` is the rate law vector, with
-``v_k(\mathbf{x}(t),t)`` the rate law for the ``k``th reaction. For example,
-for the first reaction of the repressilator above, the rate law is
-```math
-v_1(\mathbf{x}(t),t) = \frac{\alpha K^{n}}{K^{n} + \left( P_3(t) \right)^{n}}.
-```
-We can calculate each of these in Catalyst via
-```@example s1
-N = netstoichmat(repressilator)
-```
-and by using the [`oderatelaw`](@ref) function
-```@example s1
-rxs = reactions(repressilator)
-ν = oderatelaw.(rxs)
-```
-Note, as [`oderatelaw`](@ref) takes just one reaction as input we use
-broadcasting to apply it to each element of `rxs`.
-
-Let's check that this really gives the same ODEs as Catalyst. Here is what Catalyst
-generates by converting to an `ODESystem`
-```@example s1
-osys = convert(ODESystem, repressilator)
-
-# for display purposes we just pull out the right side of the equations
-odes = [eq.rhs for eq in equations(osys)]
-```
-whereas our matrix-vector representation gives
-```@example s1
-odes2 = N * ν
-```
-Let's check these are equal symbolically
-```@example s1
-isequal(odes, odes2)
-```
-
-## [Reaction complex representation](@id network_analysis_reaction_complexes)
-We now introduce a further decomposition of the RRE ODEs, which has been used to
-facilitate analysis of a variety of reaction network properties. Consider a simple
-reaction system like
-```@example s1
-rn = @reaction_network begin
- k*A, 2*A + 3*B --> A + 2*C + D
- b, C + D --> 2*A + 3*B
-end
-```
-We can think of the first reaction as converting the *reaction complex*,
-``2A+3B`` to the complex ``A+2C+D`` with rate ``kA``. Suppose we order our
-species the same way as Catalyst does, i.e.
-```math
-\begin{pmatrix}
-x_1(t)\\
-x_2(t)\\
-x_3(t)\\
-x_4(t)
-\end{pmatrix} =
-\begin{pmatrix}
-A(t)\\
-B(t)\\
-C(t)\\
-D(t)
-\end{pmatrix},
-```
-which should be the same as
-```@example s1
-species(rn)
-```
-We can describe a given reaction complex by the stoichiometric coefficients of
-each species within the complex. For the reactions in `rn` these vectors would
-be
-```math
-\begin{align*}
-2A+3B = \begin{pmatrix}
-2\\
-3\\
-0\\
-0
-\end{pmatrix}, &&
-A+2C+D = \begin{pmatrix}
-1\\
-0\\
-2\\
-1
-\end{pmatrix},
- &&
-C+D = \begin{pmatrix}
-0\\
-0\\
-1\\
-1
-\end{pmatrix}
-\end{align*}
-```
-Catalyst can calculate these representations as the columns of the complex
-stoichiometry matrix,
-```@example s1
-Z = complexstoichmat(rn)
-```
-If we have ``C`` complexes, ``Z`` is a ``M`` by ``C`` matrix with ``Z_{m c}``
-giving the stoichiometric coefficient of species ``m`` within complex ``c``.
-
-We can use this representation to provide another representation of the RRE
-ODEs. The net stoichiometry matrix can be factored as ``N = Z B``, where ``B``
-is called the incidence matrix of the reaction network,
-```@example s1
-B = incidencemat(rn)
-```
-Here ``B`` is a ``C`` by ``K`` matrix with ``B_{c k} = 1`` if complex ``c``
-appears as a product of reaction ``k``, and ``B_{c k} = -1`` if complex ``c`` is a
-substrate of reaction ``k``.
-
-Using our decomposition of ``N``, the RRE ODEs become
-```math
-\frac{dx}{dt} = Z B \mathbf{v}(\mathbf{x}(t),t).
-```
-Let's verify that ``N = Z B``,
-```@example s1
-N = netstoichmat(rn)
-N == Z*B
-```
-
-Reaction complexes give an alternative way to visualize a reaction network
-graph. Catalyst's [`complexgraph`](@ref) command will calculate the complexes of
-a network and then show how they are related. For example,
-```julia
-complexgraph(rn)
-```
-gives
-
-![Simple example complex graph](../assets/simple_complexgraph.svg)
-
-while for the repressilator we find
-```julia
-complexgraph(repressilator)
-```
-
-![Repressilator complex](../assets/repressilator_complexgraph.svg)
-
-Here ∅ represents the empty complex, black arrows show reactions converting
-substrate complexes into product complexes where the rate is just a number or
-parameter, and red arrows indicate the conversion of substrate complexes into
-product complexes where the rate is an expression involving chemical species.
-
-## [Aspects of reaction network structure](@id network_analysis_structural_aspects)
+# [Chemical Reaction Network Theory](@id network_analysis_structural_aspects)
 The reaction complex representation can be exploited via [Chemical Reaction
 Network Theory](https://en.wikipedia.org/wiki/Chemical_reaction_network_theory)
 to provide insight into possible steady state and time-dependent properties of
@@ -273,7 +63,7 @@ and,
 
 ![subnetwork_2](../assets/complex_subnets2.svg)
 
-#### [Deficiency of the network](@id network_analysis_structural_aspects_deficiency)
+# [Deficiency of the network](@id network_analysis_structural_aspects_deficiency)
 A famous theorem in Chemical Reaction Network Theory, the Deficiency Zero
 Theorem [^1], allows us to use knowledge of the net stoichiometry matrix and the
 linkage classes of a *mass action* RRE ODE system to draw conclusions about the
@@ -327,7 +117,7 @@ Quoting Feinberg [^1]
 > stoichiometry vectors] are as independent as the partition of complexes into
 > linkage classes will allow.
 
-#### [Reversibility of the network](@id network_analysis_structural_aspects_reversibility)
+# [Reversibility of the network](@id network_analysis_structural_aspects_reversibility)
 A reaction network is *reversible* if the "arrows" of the reactions are
 symmetric so that every reaction is accompanied by its reverse reaction.
 Catalyst's API provides the [`isreversible`](@ref) function to determine whether
@@ -376,7 +166,7 @@ isweaklyreversible(rn, subnets)
 Every reversible network is also weakly reversible, but not every weakly
 reversible network is reversible.
 
-#### [Deficiency Zero Theorem](@id network_analysis_structural_aspects_deficiency_zero_theorem)
+# [Deficiency Zero Theorem](@id network_analysis_structural_aspects_deficiency_zero_theorem)
 Knowing the deficiency and weak reversibility of a mass action chemical reaction
 network ODE model allows us to make inferences about the corresponding
 steady state behavior. Before illustrating how this works for one example, we
@@ -466,7 +256,16 @@ which we see is mass action and has deficiency zero, but is not weakly
 reversible. As such, we can conclude that for any choice of rate constants the
 RRE ODEs cannot have a positive equilibrium solution.
 
-## [Caching of Network Properties in `ReactionSystems`](@id network_analysis_caching_properties)
+# [Complex and Detailed Balance](@id complex_and_detailed_balance)
+The reason that the deficiency zero theorem puts such strong restrictions on the steady state properties of the reaction network is because it implies that the reaction network will be complex balanced for any set of rate constants and parameters. The fact that this holds from a purely structural property of the graph, regardless of kinetics, is what makes it so useful. In some cases it might be desirable to check complex balance (and a related property called detailed balance) on their own, though.
+
+A reaction network's steady state is complex-balanced if the total production of each *complex* is zero at the steady state. A reaction network's steady state is detailed balanced if every reaction is balanced by its reverse reaction at the steady-state (note that this requires every reaction be reversible). Detailed balance at a given steady state implies complex balance for that steady state. If a reaction network has at least one complex (detailed) balanced steady state, we say that it is complex (detailed) balanced. 
+
+Being complex (detailed) balanced is an incredibly powerful property. Having at least one positive steady state that is complex (detailed) balance implies that complex (detailed) balance obtains at *every* positive steady state, there will be exactly one steady state in every positive subspace, and this steady state is asymptotically stable. (For proofs of these results, please consult Martin Feinberg's *Foundations of Chemical Reaction Network Theory*[^1]). 
+
+Catalyst exposes API functions to determine whether a reaction network is complex or detailed balanced, `iscomplexbalanced` and `isdetailedbalanced`. Since complex and detailed balance are properties of a reaction network with a particular assignment of rate constants and parameters, both of these functions require a parameter map.
+
+# [Caching of Network Properties in `ReactionSystems`](@id network_analysis_caching_properties)
 When calling many of the network API functions, Catalyst calculates and caches
 in `rn` a variety of information. For example the first call to
 ```julia
diff --git a/docs/src/network_analysis/odes.md b/docs/src/network_analysis/odes.md
new file mode 100644
index 0000000000..de68822147
--- /dev/null
+++ b/docs/src/network_analysis/odes.md
@@ -0,0 +1,346 @@
+# [Network Analysis in Catalyst](@id network_analysis)
+
+In this tutorial we introduce several of the Catalyst API functions for network
+analysis. A complete summary of the exported functions is given in the API
+section
+[`Network-Analysis-and-Representations`](https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Network-Analysis-and-Representations).
+
+Note, currently API functions for network analysis and conservation law analysis
+do not work with constant species (currently only generated by SBMLToolkit).
+
+## [Network representation of the Repressilator `ReactionSystem`](@id network_analysis_repressilator_representation)
+We first load Catalyst and construct our model of the repressilator
+```@example s1
+using Catalyst
+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
+```
+In the [Introduction to Catalyst](@ref introduction_to_catalyst)
+tutorial we showed how the above network could be visualized as a
+species-reaction graph. There, species are represented by the nodes of the graph
+and edges show the reactions in which a given species is a substrate or product.
+```julia
+g = Graph(repressilator)
+```
+![Repressilator solution](../assets/repressilator.svg)
+
+We also showed in the [Introduction to Catalyst](@ref introduction_to_catalyst) tutorial that
+the reaction rate equation ODE model for the repressilator is
+```math
+\begin{aligned}
+\frac{dm_1(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_3}\left( t \right) \right)^{n}} - \delta {m_1}\left( t \right) + \gamma \\
+\frac{dm_2(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_1}\left( t \right) \right)^{n}} - \delta {m_2}\left( t \right) + \gamma \\
+\frac{dm_3(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_2}\left( t \right) \right)^{n}} - \delta {m_3}\left( t \right) + \gamma \\
+\frac{dP_1(t)}{dt} =& \beta {m_1}\left( t \right) - \mu {P_1}\left( t \right) \\
+\frac{dP_2(t)}{dt} =& \beta {m_2}\left( t \right) - \mu {P_2}\left( t \right) \\
+\frac{dP_3(t)}{dt} =& \beta {m_3}\left( t \right) - \mu {P_3}\left( t \right)
+\end{aligned}
+```
+
+## [Matrix-vector reaction rate equation representation](@id network_analysis_matrix_vector_representation)
+In general, reaction rate equation (RRE) ODE models for chemical reaction networks can
+be represented as a first-order system of ODEs in a compact matrix-vector notation. Suppose
+we have a reaction network with ``K`` reactions and ``M`` species, labelled by the state vector
+```math
+\mathbf{x}(t) = \begin{pmatrix} x_1(t) \\ \vdots \\ x_M(t)) \end{pmatrix}.
+```
+For the repressilator, ``\mathbf{x}(t)`` is just
+```@example s1
+x = species(repressilator)
+```
+The RRE ODEs satisfied by $\mathbf{x}(t)$ are then
+```math
+\frac{d\mathbf{x}}{dt} = N \mathbf{v}(\mathbf{x}(t),t),
+```
+where ``N`` is a constant ``M`` by ``K`` matrix with ``N_{m k}`` the net
+stoichiometric coefficient of species ``m`` in reaction ``k``.
+``\mathbf{v}(\mathbf{x}(t),t)`` is the rate law vector, with
+``v_k(\mathbf{x}(t),t)`` the rate law for the ``k``th reaction. For example,
+for the first reaction of the repressilator above, the rate law is
+```math
+v_1(\mathbf{x}(t),t) = \frac{\alpha K^{n}}{K^{n} + \left( P_3(t) \right)^{n}}.
+```
+We can calculate each of these in Catalyst via
+```@example s1
+N = netstoichmat(repressilator)
+```
+and by using the [`oderatelaw`](@ref) function
+```@example s1
+rxs = reactions(repressilator)
+ν = oderatelaw.(rxs)
+```
+Note, as [`oderatelaw`](@ref) takes just one reaction as input we use
+broadcasting to apply it to each element of `rxs`.
+
+Let's check that this really gives the same ODEs as Catalyst. Here is what Catalyst
+generates by converting to an `ODESystem`
+```@example s1
+osys = convert(ODESystem, repressilator)
+
+# for display purposes we just pull out the right side of the equations
+odes = [eq.rhs for eq in equations(osys)]
+```
+whereas our matrix-vector representation gives
+```@example s1
+odes2 = N * ν
+```
+Let's check these are equal symbolically
+```@example s1
+isequal(odes, odes2)
+```
+
+## [Reaction complex representation](@id network_analysis_reaction_complexes)
+We now introduce a further decomposition of the RRE ODEs, which has been used to
+facilitate analysis of a variety of reaction network properties. Consider a simple
+reaction system like
+```@example s1
+rn = @reaction_network begin
+ k*A, 2*A + 3*B --> A + 2*C + D
+ b, C + D --> 2*A + 3*B
+end
+```
+We can think of the first reaction as converting the *reaction complex*,
+``2A+3B`` to the complex ``A+2C+D`` with rate ``kA``. Suppose we order our
+species the same way as Catalyst does, i.e.
+```math
+\begin{pmatrix}
+x_1(t)\\
+x_2(t)\\
+x_3(t)\\
+x_4(t)
+\end{pmatrix} =
+\begin{pmatrix}
+A(t)\\
+B(t)\\
+C(t)\\
+D(t)
+\end{pmatrix},
+```
+which should be the same as
+```@example s1
+species(rn)
+```
+We can describe a given reaction complex by the stoichiometric coefficients of
+each species within the complex. For the reactions in `rn` these vectors would
+be
+```math
+\begin{align*}
+2A+3B = \begin{pmatrix}
+2\\
+3\\
+0\\
+0
+\end{pmatrix}, &&
+A+2C+D = \begin{pmatrix}
+1\\
+0\\
+2\\
+1
+\end{pmatrix},
+ &&
+C+D = \begin{pmatrix}
+0\\
+0\\
+1\\
+1
+\end{pmatrix}
+\end{align*}
+```
+Catalyst can calculate these representations as the columns of the complex
+stoichiometry matrix,
+```@example s1
+Z = complexstoichmat(rn)
+```
+If we have ``C`` complexes, ``Z`` is a ``M`` by ``C`` matrix with ``Z_{m c}``
+giving the stoichiometric coefficient of species ``m`` within complex ``c``.
+
+We can use this representation to provide another representation of the RRE
+ODEs. The net stoichiometry matrix can be factored as ``N = Z B``, where ``B``
+is called the incidence matrix of the reaction network,
+```@example s1
+B = incidencemat(rn)
+```
+Here ``B`` is a ``C`` by ``K`` matrix with ``B_{c k} = 1`` if complex ``c``
+appears as a product of reaction ``k``, and ``B_{c k} = -1`` if complex ``c`` is a
+substrate of reaction ``k``.
+
+Using our decomposition of ``N``, the RRE ODEs become
+```math
+\frac{dx}{dt} = Z B \mathbf{v}(\mathbf{x}(t),t).
+```
+Let's verify that ``N = Z B``,
+```@example s1
+N = netstoichmat(rn)
+N == Z*B
+```
+
+Reaction complexes give an alternative way to visualize a reaction network
+graph. Catalyst's [`complexgraph`](@ref) command will calculate the complexes of
+a network and then show how they are related. For example,
+```julia
+complexgraph(rn)
+```
+gives
+
+![Simple example complex graph](../assets/simple_complexgraph.svg)
+
+while for the repressilator we find
+```julia
+complexgraph(repressilator)
+```
+
+![Repressilator complex](../assets/repressilator_complexgraph.svg)
+
+Here ∅ represents the empty complex, black arrows show reactions converting
+substrate complexes into product complexes where the rate is just a number or
+parameter, and red arrows indicate the conversion of substrate complexes into
+product complexes where the rate is an expression involving chemical species.
+
+# Full decomposition of the reaction network ODEs (flux matrix and mass-action vector)
+So far we have covered two equivalent descriptions of the chemical reaction network ODEs: 
+```math
+\begin{align}
+\frac{d\mathbf{x}}{dt} &= N \mathbf{v}(\mathbf{x}(t),t) \\
+&= Z B \mathbf{v}(\mathbf{x}(t),t)
+\end{align}
+```
+
+In this section we discuss a further decomposition of the ODEs. Recall that the reaction rate vector $\mathbf{v}$, which is a vector of length $R$ whose elements are the rate expressions for each reaction. Its elements can be written as 
+```math
+\mathbf{v}_{y \rightarrow y'} = k_{y \rightarrow y'} \mathbf{x}^y,
+```
+where $\mathbf{x}^y = \prod_s x_s^{y_s}$, the mass-action product of the complex $y$, where $y$ is the substrate complex of the reaction $y \rightarrow y'$. We can define a new vector called the mass action vector $\Phi(\mathbf{x}(t))$, a vector of length $C$ whose elements are the mass action products of each complex: 
+```@example s1
+Φ = massactionvector(rn)
+```
+
+An important thing to note is this function assumes [combinatoric ratelaws](@introduction_to_catalyst_ratelaws), meaning that mass-action products will get rescaled by factorial factors. For instance, note that the mass-action product for the complex `2A + 3B` has a factor of 1/12, corresponding to 1/(2! 3!). This option can be turned off with `combinatoric_ratelaws = false`.
+
+```@example s1
+Φ_2 = massactionvector(rn; combinatoric_ratelaws = false)
+```
+
+Then the reaction rate vector $\mathbf{v}$ can be written as 
+```math
+\mathbf{v}(\mathbf{x}(t)) = K \Phi(\mathbf{x}(t))
+```
+where $K$ is an $R$-by-$C$ matrix called the flux matrix, where $K_{rc}$ is the rate constant of reaction $r$ if $c$ is the index of the substrate complex of reaction $r$, and 0 otherwise. In Catalyst, the API function for $K$ is `fluxmat`: 
+```@example s1
+K = fluxmat(rn)
+```
+
+Since we have that $\mathbf{v} = K\Phi$, we can rewrite the above decompositions as follows: 
+```math
+\frac{d\mathbf{x}}{dt} &= N \mathbf{v}(\mathbf{x}(t),t) \\
+&= N K \Phi(\mathbf{x}(t),t) \\
+&= Z B K \Phi(\mathbf{x}(t),t). 
+```
+
+The final matrix to discuss is the product of $A_k = BK$, which is a $C$-by-$C$ matrix that turns out to be exactly the negative of the [graph Laplacian](https://en.wikipedia.org/wiki/Laplacian_matrix) of the weighted graph whose nodes are reaction complexes and whose edges represent reactions, weighted by the rate constants. The API function for $A_k$ is the `laplacianmat`: 
+```@example s1
+A_k = incidencemat(rn)
+```
+We can check that 
+```@example s1
+A_k == B * K
+```
+
+In sum, we have that
+```math
+\frac{d\mathbf{x}}{dt} &= N \mathbf{v}(\mathbf{x}(t),t) \\
+&= N K \Phi(\mathbf{x}(t),t) \\
+&= Z B K \Phi(\mathbf{x}(t),t). 
+&= Z A_k \Phi(\mathbf{x}(t),t). 
+```
+
+All three of the objects introduced in this section (the flux matrix, mass-action vector, Laplacian matrix) will return symbolic outputs by default, but can be made to return numerical outputs if values are specified. 
+For example, `massactionvector` will return a numerical output if a set of species concentrations is supplied using a dictionary, tuple, or vector of Symbol-value pairs.
+```@example s1
+concmap = Dict([:A => 3., :B => 5., :C => 2.4, :D => 1.5])
+massactionvector(rn, concmap)
+```
+
+`fluxmat` and `laplacianmat` will return numeric matrices if a set of rate constants and other aprameters are supplied the same way.
+```@example s1
+parammap = Dict([:k => 12., b => 8.])
+fluxmat(rn, parammap)
+```
+
+```@example s1
+laplacianmat(rn, parammap)
+```
+
+# Symbolic ODE functions
+In some cases it might be useful to generate the function defining the system of ODEs as a symbolic Julia function that can be used for further analysis. This can be done using Symbolics' [`build_function`](https://docs.sciml.ai/Symbolics/stable/getting_started/#Building-Functions), which takes a symbolic expression and a set of desired arguments, and converts it into a Julia function taking those arguments.
+
+Let's build the full symbolic function corresponding to our ODE system. `build_function` will return two expressions, one for a function that outputs a new vector for the result, and one for a function that modifies the input in-place. Either expression can then be evaluated to return a Julia function.
+```@example s1
+parammap = Dict([:k => 12., b => 8.])
+K = fluxmat(rn, parammap)
+odes = N * K * Φ
+f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn))
+ode_func = eval(f_oop_expr)
+
+concvec = [3., 5., 2., 6.]
+ode_func(concvec)
+```
+The generated `ode_func` now corresponds to the $f(\mathbf{x}(t))$ on the right-hand side of $\frac{d\mathbf{x}(t)}{dt} = f(\mathbf{x}(t))$. Given a vector of species concentrations $c$, `ode_func` will return the rate of change of each species. Steady state concentration vectors `c_ss` will satisfy `ode_func(c_ss) = zeros(length(species(rn)))`.
+
+Above we have generated a numeric rate matrix to substitute the rate constants into the symbolic expressions. We could have used a symbolic rate matrix, but then we would need to define the parameters `k, b`, so that the function `ode_func` knows what `k` and `b` in its output refer to.
+```@example s1
+@parameters k b 
+K = fluxmat(rn)
+odes = N * K * Φ
+f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn))
+ode_func = eval(f_oop_expr)
+
+concvec = [3., 5., 2., 6.]
+ode_func(concvec)
+```
+
+Alternatively, if we use a symbolic rate matrix, we could define our function to take in both species concentrations and parameter values as arguments:
+```@example s1
+K = fluxmat(rn)
+odes = N * K * Φ
+f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn), parameters(rn))
+ode_func = eval(f_oop_expr)
+
+concvec = [3., 5., 2., 6]; rateconsts = [12., 4.]
+ode_func(concvec, rateconsts)
+```
+
+Note also that `ode_func` can take any vector with the right dimension (i.e. the number of species), not just a vector of `Number`, so it can be used to build, e.g. a vector of polynomials in Nemo for commutative algebraic methods.
+
+# API Section for matrices and vectors
+We have that: 
+- $N$ is the `netstoichmat`
+- $Z$ is the `complexstoichmat`
+- $B$ is the `incidencemat`
+- $K$ is the `fluxmat`
+- $A_k$ is the `laplacianmat`
+- $\Phi$ is the `massactionvector`
+
+```@docs
+netstoichmat
+complexstoichmat
+incidencemat
+fluxmat
+laplacianmat
+massactionvector
+```
+
+---
+## References
+[^1]: [Feinberg, M. *Foundations of Chemical Reaction Network Theory*, Applied Mathematical Sciences 202, Springer (2019).](https://link.springer.com/book/10.1007/978-3-030-03858-8?noAccess=true)

From bb070e9a4158aac2d53473cc6f956d790f503f75 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Sun, 8 Dec 2024 15:11:33 -0500
Subject: [PATCH 02/24] up

---
 docs/pages.jl                                 |  5 ++
 docs/src/network_analysis/crnt.md             | 80 ++++++++++---------
 .../network_analysis/network_properties.md    | 37 +++++++++
 docs/src/network_analysis/odes.md             | 21 +++--
 4 files changed, 100 insertions(+), 43 deletions(-)
 create mode 100644 docs/src/network_analysis/network_properties.md

diff --git a/docs/pages.jl b/docs/pages.jl
index f053a5b6ce..616b6550aa 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -25,6 +25,11 @@ pages = Any[
             "model_creation/examples/smoluchowski_coagulation_equation.md"
         ]
     ],
+    "Network Analysis" => Any[
+        "network_analysis/odes.md",
+        "network_analysis/crnt.md",
+        "network_analysis/",
+    ],
     "Model simulation and visualization" => Any[
         "model_simulation/simulation_introduction.md",
         "model_simulation/simulation_plotting.md",
diff --git a/docs/src/network_analysis/crnt.md b/docs/src/network_analysis/crnt.md
index 9a445da2a0..d3194189b5 100644
--- a/docs/src/network_analysis/crnt.md
+++ b/docs/src/network_analysis/crnt.md
@@ -24,7 +24,7 @@ complexgraph(rn)
 
 ![network_1](../assets/complex_rn.svg)
 
-#### [Linkage classes and sub-networks of the reaction network](@id network_analysis_structural_aspects_linkage)
+# [Linkage classes and sub-networks of the reaction network](@id network_analysis_structural_aspects_linkage)
 The preceding reaction complex graph shows that `rn` is composed of two
 disconnected sub-graphs, one containing the complexes ``A+B``, ``C``, ``D+E``, and
 ``F``, the other containing the complexes ``2A``, ``B + G``, and ``H``. These sets,
@@ -63,6 +63,7 @@ and,
 
 ![subnetwork_2](../assets/complex_subnets2.svg)
 
+
 # [Deficiency of the network](@id network_analysis_structural_aspects_deficiency)
 A famous theorem in Chemical Reaction Network Theory, the Deficiency Zero
 Theorem [^1], allows us to use knowledge of the net stoichiometry matrix and the
@@ -117,6 +118,7 @@ Quoting Feinberg [^1]
 > stoichiometry vectors] are as independent as the partition of complexes into
 > linkage classes will allow.
 
+
 # [Reversibility of the network](@id network_analysis_structural_aspects_reversibility)
 A reaction network is *reversible* if the "arrows" of the reactions are
 symmetric so that every reaction is accompanied by its reverse reaction.
@@ -256,51 +258,55 @@ which we see is mass action and has deficiency zero, but is not weakly
 reversible. As such, we can conclude that for any choice of rate constants the
 RRE ODEs cannot have a positive equilibrium solution.
 
-# [Complex and Detailed Balance](@id complex_and_detailed_balance)
-The reason that the deficiency zero theorem puts such strong restrictions on the steady state properties of the reaction network is because it implies that the reaction network will be complex balanced for any set of rate constants and parameters. The fact that this holds from a purely structural property of the graph, regardless of kinetics, is what makes it so useful. In some cases it might be desirable to check complex balance (and a related property called detailed balance) on their own, though.
+```@docs
+satisfiesdeficiencyzero
+```
 
-A reaction network's steady state is complex-balanced if the total production of each *complex* is zero at the steady state. A reaction network's steady state is detailed balanced if every reaction is balanced by its reverse reaction at the steady-state (note that this requires every reaction be reversible). Detailed balance at a given steady state implies complex balance for that steady state. If a reaction network has at least one complex (detailed) balanced steady state, we say that it is complex (detailed) balanced. 
+# Deficiency One Theorem
+Very analogous to the deficiency zero theorem is the deficiency one theorem. The deficiency one theorem applies to a network with the following properties: 
+1. The deficiency of each *linkage class* of the network is at most 1, 
+2. The sum of the linkage class deficiencies is the total deficiency of the network, and
+3. Each linkage class has at most one terminal linkage class. 
 
-Being complex (detailed) balanced is an incredibly powerful property. Having at least one positive steady state that is complex (detailed) balance implies that complex (detailed) balance obtains at *every* positive steady state, there will be exactly one steady state in every positive subspace, and this steady state is asymptotically stable. (For proofs of these results, please consult Martin Feinberg's *Foundations of Chemical Reaction Network Theory*[^1]). 
+If these conditions are met, then the network will have at most one steady state in each stoichiometric compatibility class for any choice of rate constants and parameters.
+Unlike the deficiency zero theorem, networks obeying the deficiency one theorem are not guaranteed to have stable solutions.
+```@docs
+satisfiesdeficiencyone
+```
 
-Catalyst exposes API functions to determine whether a reaction network is complex or detailed balanced, `iscomplexbalanced` and `isdetailedbalanced`. Since complex and detailed balance are properties of a reaction network with a particular assignment of rate constants and parameters, both of these functions require a parameter map.
+# [Complex and Detailed Balance](@id complex_and_detailed_balance)
+A reaction network's steady state is **complex-balanced** if the total production of each *complex* is zero at the steady state. A reaction network's steady state is **detailed balanced** if every reaction is balanced by its reverse reaction at the steady-state (this corresponds to the usual notion of chemical equilibrium; note that this requires every reaction be reversible). 
 
-# [Caching of Network Properties in `ReactionSystems`](@id network_analysis_caching_properties)
-When calling many of the network API functions, Catalyst calculates and caches
-in `rn` a variety of information. For example the first call to
-```julia
-rcs,B = reactioncomplexes(rn)
-```
-calculates, caches, and returns the reaction complexes, `rcs`, and the incidence
-matrix, `B`, of `rn`. Subsequent calls simply return `rcs` and `B` from the
-cache.
+Note that detailed balance at a given steady state implies complex balance for that steady state, i.e. detailed balance is a stronger property than complex balance.
 
-Similarly, the first call to
-```julia
-N = netstoichmat(rn)
+Remarkably, having just one positive steady state that is complex (detailed) balance implies that complex (detailed) balance obtains at *every* positive steady state, so we say that a network is complex (detailed) balanced if any one of its steady states are complex (detailed) balanced. Additionally, there will be exactly one steady state in every positive stoichiometric compatibility class, and this steady state is asymptotically stable. (For proofs of these results, please consult Martin Feinberg's *Foundations of Chemical Reaction Network Theory*[^1]). So knowing that a network is complex balanced is really quite powerful.
+
+Let's check that the reaction network defined above is complex balanced by providing a set of rates:
+```@example s1
+rates = Dict([:k1 => 2.4, :k2 => 4., :k3 => 10., :k4 => 5.5, :k5 => 0.4])
+iscomplexbalanced(rn, rates)
 ```
-calculates, caches and returns the net stoichiometry matrix. Subsequent calls
-then simply return the cached value of `N`. Caching such information means users
-do not need to manually know which subsets of network properties are needed for
-a given calculation (like the deficiency). Generally only
-```julia
-rcs,B = reactioncomplexes(rn)    # must be called once to cache rcs and B
-any_other_network_property(rn)
+
+Complex balance obtains for some sets of rates but not others: 
+```@example s1
 ```
-should work to calculate a desired network property, with the API doc strings
-indicating when `reactioncomplexes(rn)` must be called at least once before a
-given function is used.
-
-Because of the caching of network properties, subsequent calls to most API
-functions will be fast, simply returning the previously calculated and cached
-values. In some cases it may be desirable to reset the cache and recalculate
-these properties. This can be done by calling
-```julia
-Catalyst.reset_networkproperties!(rn)
+
+We can do a similar check for detailed balance. Let us make the reaction network 
+```@example s1
+rn1 = @reaction_network begin
+  (k1,k2),A <--> 2B
+  (k3,k4), A + C <--> D
+  (k5,k6), B + E --> C + D
+end
+isdetailedbalanced(rn, rates)
 ```
-Network property functions will then recalculate their associated properties and
-cache the new values the next time they are called.
 
+The reason that the deficiency zero theorem puts such strong restrictions on the steady state properties of the reaction network is because it implies that the reaction network will be complex balanced for any set of rate constants and parameters. The fact that this holds from a purely structural property of the graph, regardless of kinetics, is what makes it so useful. But in some cases it might be desirable to check complex balance on its own, as for higher deficiency networks.
+
+```@docs
+iscomplexbalanced
+isdetailedbalanced
+```
 ---
 ## References
 [^1]: [Feinberg, M. *Foundations of Chemical Reaction Network Theory*, Applied Mathematical Sciences 202, Springer (2019).](https://link.springer.com/book/10.1007/978-3-030-03858-8?noAccess=true)
diff --git a/docs/src/network_analysis/network_properties.md b/docs/src/network_analysis/network_properties.md
new file mode 100644
index 0000000000..bced9cb89c
--- /dev/null
+++ b/docs/src/network_analysis/network_properties.md
@@ -0,0 +1,37 @@
+# The `NetworkProperties` Object
+
+# [Caching of Network Properties in `ReactionSystems`](@id network_analysis_caching_properties)
+When calling many of the network API functions, Catalyst calculates and caches
+in `rn` a variety of information. For example the first call to
+```julia
+rcs,B = reactioncomplexes(rn)
+```
+calculates, caches, and returns the reaction complexes, `rcs`, and the incidence
+matrix, `B`, of `rn`. Subsequent calls simply return `rcs` and `B` from the
+cache.
+
+Similarly, the first call to
+```julia
+N = netstoichmat(rn)
+```
+calculates, caches and returns the net stoichiometry matrix. Subsequent calls
+then simply return the cached value of `N`. Caching such information means users
+do not need to manually know which subsets of network properties are needed for
+a given calculation (like the deficiency). Generally only
+```julia
+rcs,B = reactioncomplexes(rn)    # must be called once to cache rcs and B
+any_other_network_property(rn)
+```
+should work to calculate a desired network property, with the API doc strings
+indicating when `reactioncomplexes(rn)` must be called at least once before a
+given function is used.
+
+Because of the caching of network properties, subsequent calls to most API
+functions will be fast, simply returning the previously calculated and cached
+values. In some cases it may be desirable to reset the cache and recalculate
+these properties. This can be done by calling
+```julia
+Catalyst.reset_networkproperties!(rn)
+```
+Network property functions will then recalculate their associated properties and
+cache the new values the next time they are called.
diff --git a/docs/src/network_analysis/odes.md b/docs/src/network_analysis/odes.md
index de68822147..c796047122 100644
--- a/docs/src/network_analysis/odes.md
+++ b/docs/src/network_analysis/odes.md
@@ -323,14 +323,23 @@ ode_func(concvec, rateconsts)
 
 Note also that `ode_func` can take any vector with the right dimension (i.e. the number of species), not just a vector of `Number`, so it can be used to build, e.g. a vector of polynomials in Nemo for commutative algebraic methods.
 
+# Properties of matrix null spaces
+The null spaces of the matrices discussed in this section often have special meaning. Below we will discuss some of these properties.
+
+Recall that we may write the net stoichiometry matrix ``N = YB``.
+
+[Conservation laws](@ref conservation_laws) arise as left null eigenvectors of the net stoichiometry matrix ``N``, and cycles arise as right null eigenvectors of the stoichiometry matrix. A cycle may be understood as a sequence of reactions that leaves the overall species composition unchanged. These do not necessarily have to correspond to actual cycles in the graph.
+
+[Complex balance](@ref complex_and_detailed_balance) can be compactly formulated as the following: a set of steady state reaction fluxes is complex-balanced if it is in the nullspace of the incidence matrix ``B``.
+
 # API Section for matrices and vectors
 We have that: 
-- $N$ is the `netstoichmat`
-- $Z$ is the `complexstoichmat`
-- $B$ is the `incidencemat`
-- $K$ is the `fluxmat`
-- $A_k$ is the `laplacianmat`
-- $\Phi$ is the `massactionvector`
+- ``N`` is the `netstoichmat`
+- ``Z`` is the `complexstoichmat`
+- ``B`` is the `incidencemat`
+- ``K`` is the `fluxmat`
+- ``A_k`` is the `laplacianmat`
+- ``\Phi`` is the `massactionvector`
 
 ```@docs
 netstoichmat

From 9f5ca724a725ccff910b855f8f505b31ecbe8244 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 28 Jan 2025 11:28:25 -0500
Subject: [PATCH 03/24] up

---
 docs/pages.jl                                 |  3 +-
 .../src/network_analysis/advanced_analysis.md |  3 ++
 docs/src/network_analysis/crnt.md             | 40 +++++++++++++++----
 .../network_analysis/network_properties.md    |  4 ++
 docs/src/network_analysis/odes.md             | 24 +++++------
 src/Catalyst.jl                               |  1 +
 6 files changed, 54 insertions(+), 21 deletions(-)
 create mode 100644 docs/src/network_analysis/advanced_analysis.md

diff --git a/docs/pages.jl b/docs/pages.jl
index 4eac095897..18ffe2e347 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -28,7 +28,8 @@ pages = Any[
     "Network Analysis" => Any[
         "network_analysis/odes.md",
         "network_analysis/crnt.md",
-        "network_analysis/",
+        "network_analysis/network_properties.md",
+        "network_analysis/advanced_analysis.md"
     ],
     "Model simulation and visualization" => Any[
         "model_simulation/simulation_introduction.md",
diff --git a/docs/src/network_analysis/advanced_analysis.md b/docs/src/network_analysis/advanced_analysis.md
new file mode 100644
index 0000000000..3bf024ee1e
--- /dev/null
+++ b/docs/src/network_analysis/advanced_analysis.md
@@ -0,0 +1,3 @@
+# [Advanced Analysis using CatalystNetworkAnalysis](@ref CatalystNetworkAnalysis)
+
+For more advanced chemical reaction network analysis, there is a special CatalystNetworkAnalysis package. This package has functionality for a wider range of analysis for chemical reaction networks, including the following:
diff --git a/docs/src/network_analysis/crnt.md b/docs/src/network_analysis/crnt.md
index d3194189b5..a42b4f6e3e 100644
--- a/docs/src/network_analysis/crnt.md
+++ b/docs/src/network_analysis/crnt.md
@@ -1,10 +1,15 @@
 # [Chemical Reaction Network Theory](@id network_analysis_structural_aspects)
-The reaction complex representation can be exploited via [Chemical Reaction
-Network Theory](https://en.wikipedia.org/wiki/Chemical_reaction_network_theory)
-to provide insight into possible steady state and time-dependent properties of
-RRE ODE models and stochastic chemical kinetics models. We'll now illustrate
-some of the types of network properties that Catalyst can determine, using the
-reaction complex representation in these calculations.
+The systems of ODEs or stochastic chemical kinetics models that arise from chemical 
+reaction networks can often have their steady-state properties known in advance, 
+simply by analyzing the graph structure of the network. The subfield of chemistry 
+and math studying this relationship is called [Chemical Reaction Network Theory](https://en.wikipedia.org/wiki/Chemical_reaction_network_theory).
+
+Broadly, results from chemical reaction network theory relate a purely
+graph-structural property (e.g. deficiency) to dynamical properties of the reaction system 
+(e.g. complex balance).
+
+We'll now illustrate some of the types of network properties that Catalyst can determine, 
+using the [reaction complex representation](@ref network_analysis_reaction_complexes) in these calculations.
 
 Consider the following reaction network.
 ```@example s1
@@ -67,7 +72,7 @@ and,
 # [Deficiency of the network](@id network_analysis_structural_aspects_deficiency)
 A famous theorem in Chemical Reaction Network Theory, the Deficiency Zero
 Theorem [^1], allows us to use knowledge of the net stoichiometry matrix and the
-linkage classes of a *mass action* RRE ODE system to draw conclusions about the
+linkage classes of a *mass action* [RRE ODE system](@ref network_analysis_matrix_vector_representation) to draw conclusions about the
 system's possible steady states. In this section we'll see how Catalyst can
 calculate a network's deficiency.
 
@@ -274,7 +279,7 @@ Unlike the deficiency zero theorem, networks obeying the deficiency one theorem
 satisfiesdeficiencyone
 ```
 
-# [Complex and Detailed Balance](@id complex_and_detailed_balance)
+# [Complex and Detailed Balance](@id network_analysis_complex_and_detailed_balance)
 A reaction network's steady state is **complex-balanced** if the total production of each *complex* is zero at the steady state. A reaction network's steady state is **detailed balanced** if every reaction is balanced by its reverse reaction at the steady-state (this corresponds to the usual notion of chemical equilibrium; note that this requires every reaction be reversible). 
 
 Note that detailed balance at a given steady state implies complex balance for that steady state, i.e. detailed balance is a stronger property than complex balance.
@@ -307,6 +312,25 @@ The reason that the deficiency zero theorem puts such strong restrictions on the
 iscomplexbalanced
 isdetailedbalanced
 ```
+
+# [Concentration Robustness](@id network_analysis_concentration_robustness)
+Certain reaction networks have species that do not change their concentration,
+regardless of the system is perturbed to a different stoichiometric compatibility
+class. This is a very useful property to have in biological contexts, where it might
+be important to keep the concentration of a critical species relatively stable in
+the face of changes in its environment. 
+
+Determining every species with concentration-robustness in a network is in general very difficult. However, there are certain cases where there are sufficient conditions
+that can be checked relatively easily. One example is for deficiency one networks.
+
+**Theorem (a sufficient condition for concentration robustness for deficiency one networks)**: If there are two *non-terminal* reaction complexes that differ only in species ``s``, then the system is absolutely concentration robust with respect to ``s``. 
+
+This is the check provided by the API function `robustspecies(rn)`. More general concentration robustness analysis can be done using the [CatalystNetworkAnalysis](@ref) package.
+
+```@docs
+robustspecies
+```
+
 ---
 ## References
 [^1]: [Feinberg, M. *Foundations of Chemical Reaction Network Theory*, Applied Mathematical Sciences 202, Springer (2019).](https://link.springer.com/book/10.1007/978-3-030-03858-8?noAccess=true)
diff --git a/docs/src/network_analysis/network_properties.md b/docs/src/network_analysis/network_properties.md
index bced9cb89c..92b25c5055 100644
--- a/docs/src/network_analysis/network_properties.md
+++ b/docs/src/network_analysis/network_properties.md
@@ -35,3 +35,7 @@ Catalyst.reset_networkproperties!(rn)
 ```
 Network property functions will then recalculate their associated properties and
 cache the new values the next time they are called.
+
+```@docs
+NetworkProperties
+```
diff --git a/docs/src/network_analysis/odes.md b/docs/src/network_analysis/odes.md
index c796047122..1c7b718f10 100644
--- a/docs/src/network_analysis/odes.md
+++ b/docs/src/network_analysis/odes.md
@@ -291,23 +291,23 @@ parammap = Dict([:k => 12., b => 8.])
 K = fluxmat(rn, parammap)
 odes = N * K * Φ
 f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn))
-ode_func = eval(f_oop_expr)
+f = eval(f_oop_expr)
 
-concvec = [3., 5., 2., 6.]
-ode_func(concvec)
+c = [3., 5., 2., 6.]
+f(c)
 ```
-The generated `ode_func` now corresponds to the $f(\mathbf{x}(t))$ on the right-hand side of $\frac{d\mathbf{x}(t)}{dt} = f(\mathbf{x}(t))$. Given a vector of species concentrations $c$, `ode_func` will return the rate of change of each species. Steady state concentration vectors `c_ss` will satisfy `ode_func(c_ss) = zeros(length(species(rn)))`.
+The generated `f` now corresponds to the $f(\mathbf{x}(t))$ on the right-hand side of $\frac{d\mathbf{x}(t)}{dt} = f(\mathbf{x}(t))$. Given a vector of species concentrations $c$, `ode_func` will return the rate of change of each species. Steady state concentration vectors `c_ss` will satisfy `f(c_ss) = zeros(length(species(rn)))`.
 
-Above we have generated a numeric rate matrix to substitute the rate constants into the symbolic expressions. We could have used a symbolic rate matrix, but then we would need to define the parameters `k, b`, so that the function `ode_func` knows what `k` and `b` in its output refer to.
+Above we have generated a numeric rate matrix to substitute the rate constants into the symbolic expressions. We could have used a symbolic rate matrix, but then we would need to define the parameters `k, b`, so that the function `f` knows what `k` and `b` in its output refer to.
 ```@example s1
 @parameters k b 
 K = fluxmat(rn)
 odes = N * K * Φ
 f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn))
-ode_func = eval(f_oop_expr)
+f = eval(f_oop_expr)
 
-concvec = [3., 5., 2., 6.]
-ode_func(concvec)
+c = [3., 5., 2., 6.]
+f(c)
 ```
 
 Alternatively, if we use a symbolic rate matrix, we could define our function to take in both species concentrations and parameter values as arguments:
@@ -315,13 +315,13 @@ Alternatively, if we use a symbolic rate matrix, we could define our function to
 K = fluxmat(rn)
 odes = N * K * Φ
 f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn), parameters(rn))
-ode_func = eval(f_oop_expr)
+f = eval(f_oop_expr)
 
-concvec = [3., 5., 2., 6]; rateconsts = [12., 4.]
-ode_func(concvec, rateconsts)
+c = [3., 5., 2., 6]; ks = [12., 4.]
+f(c, ks)
 ```
 
-Note also that `ode_func` can take any vector with the right dimension (i.e. the number of species), not just a vector of `Number`, so it can be used to build, e.g. a vector of polynomials in Nemo for commutative algebraic methods.
+Note also that `f` can take any vector with the right dimension (i.e. the number of species), not just a vector of `Number`, so it can be used to build, e.g. a vector of polynomials in Nemo for commutative algebraic methods.
 
 # Properties of matrix null spaces
 The null spaces of the matrices discussed in this section often have special meaning. Below we will discuss some of these properties.
diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index 73fbc7586f..430c847b11 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -133,6 +133,7 @@ export complexoutgoingmat, incidencematgraph, linkageclasses, stronglinkageclass
        terminallinkageclasses, deficiency, subnetworks
 export linkagedeficiencies, isreversible, isweaklyreversible
 export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants
+export iscomplexbalanced, isdetailedbalanced, robustspecies
 
 # registers CRN specific functions using Symbolics.jl
 include("registered_functions.jl")

From 0b7c4de2364bf8bb53a533ab6b7fd011210f34b0 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 28 Jan 2025 11:28:40 -0500
Subject: [PATCH 04/24] up

---
 docs/src/network_analysis/crnt.md | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/src/network_analysis/crnt.md b/docs/src/network_analysis/crnt.md
index a42b4f6e3e..569cf8d277 100644
--- a/docs/src/network_analysis/crnt.md
+++ b/docs/src/network_analysis/crnt.md
@@ -1,14 +1,14 @@
 # [Chemical Reaction Network Theory](@id network_analysis_structural_aspects)
-The systems of ODEs or stochastic chemical kinetics models that arise from chemical 
-reaction networks can often have their steady-state properties known in advance, 
-simply by analyzing the graph structure of the network. The subfield of chemistry 
+The systems of ODEs or stochastic chemical kinetics models that arise from chemical
+reaction networks can often have their steady-state properties known in advance,
+simply by analyzing the graph structure of the network. The subfield of chemistry
 and math studying this relationship is called [Chemical Reaction Network Theory](https://en.wikipedia.org/wiki/Chemical_reaction_network_theory).
 
 Broadly, results from chemical reaction network theory relate a purely
-graph-structural property (e.g. deficiency) to dynamical properties of the reaction system 
+graph-structural property (e.g. deficiency) to dynamical properties of the reaction system
 (e.g. complex balance).
 
-We'll now illustrate some of the types of network properties that Catalyst can determine, 
+We'll now illustrate some of the types of network properties that Catalyst can determine,
 using the [reaction complex representation](@ref network_analysis_reaction_complexes) in these calculations.
 
 Consider the following reaction network.

From 569c36fcb9dd652e5964471b851708c358ec20b4 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 28 Jan 2025 11:32:20 -0500
Subject: [PATCH 05/24] Up

---
 src/Catalyst.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index 75630f0ab5..185bf7afa0 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -115,11 +115,12 @@ export @reaction_network, @network_component, @reaction, @species
 # Network analysis functionality.
 include("network_analysis.jl")
 export reactioncomplexmap, reactioncomplexes, incidencemat
-export complexstoichmat, laplacianmat, fluxmat, massactionvector, complexoutgoingmat
+export complexstoichmat, laplacianmat, fluxmat, massactionvector, complexoutgoingmat, adjacencymat
 export incidencematgraph, linkageclasses, stronglinkageclasses,
        terminallinkageclasses, deficiency, subnetworks
 export linkagedeficiencies, isreversible, isweaklyreversible
 export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants
+export satisfiesdeficiencyone, satisfiesdeficiencyzero
 export iscomplexbalanced, isdetailedbalanced, robustspecies
 
 # registers CRN specific functions using Symbolics.jl

From c0d309ddfb1767d3c813f084c3eba6d53d56db83 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 28 Jan 2025 14:54:27 -0500
Subject: [PATCH 06/24] add API functions

---
 docs/src/api.md                   | 16 ++++++++++++++--
 docs/src/index.md                 |  2 +-
 docs/src/network_analysis/crnt.md |  9 +++++----
 docs/src/network_analysis/odes.md | 11 ++++++-----
 4 files changed, 26 insertions(+), 12 deletions(-)

diff --git a/docs/src/api.md b/docs/src/api.md
index f78a5f70e4..a47538c044 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -242,6 +242,9 @@ Catalyst.flatten
 Note, currently API functions for network analysis and conservation law analysis
 do not work with constant species (currently only generated by SBMLToolkit).
 
+For more information about these functions, please see the sections of the docs on
+[network ODE representation](@ref network_analysis_odes) and [chemical reaction network theory](@ref network_analysis_structural_aspects).
+
 ```@docs
 conservationlaws
 conservedquantities
@@ -252,15 +255,24 @@ ReactionComplex
 reactioncomplexmap
 reactioncomplexes
 incidencemat
+incidencematgraph
 complexstoichmat
 complexoutgoingmat
-incidencematgraph
+fluxmat
+adjacencymat
+laplacianmat
+massactionvector
 linkageclasses
 deficiency
-subnetworks
 linkagedeficiencies
+satisfiesdeficiencyone
+satisfiesdeficiencyzero
+subnetworks
 isreversible
 isweaklyreversible
+iscomplexbalanced
+isdetailedbalanced
+robustspecies
 reset_networkproperties!
 ```
 
diff --git a/docs/src/index.md b/docs/src/index.md
index c30cd42ca3..1360231470 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -23,7 +23,7 @@ etc).
 - The [Catalyst.jl API](@ref api) provides functionality for building networks programmatically and for composing multiple networks together.
 - Leveraging ModelingToolkit, generated models can be converted to symbolic reaction rate equation ODE models, symbolic Chemical Langevin Equation models, and symbolic stochastic chemical kinetics (jump process) models. These can be simulated using any [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/) [ODE/SDE/jump solver](@ref simulation_intro), and can be used within `EnsembleProblem`s for carrying out [parallelized parameter sweeps and statistical sampling](@ref ensemble_simulations). Plot recipes are available for [visualization of all solutions](@ref simulation_plotting).
 - Non-integer (e.g. `Float64`) stoichiometric coefficients [are supported](@ref dsl_description_stoichiometries_decimal) for generating ODE models, and symbolic expressions for stoichiometric coefficients [are supported](@ref parametric_stoichiometry) for all system types.
-- A [network analysis suite](@ref network_analysis) permits the computation of linkage classes, deficiencies, reversibility, and other network properties.
+- A [network analysis suite](@ref network_analysis_structural_aspects) permits the computation of linkage classes, deficiencies, reversibility, and other network properties.
 - [Conservation laws can be detected and utilized](@ref conservation_laws) to reduce system sizes, and to generate non-singular Jacobians (e.g. during conversion to ODEs, SDEs, and steady state equations).
 - Catalyst reaction network models can be [coupled with differential and algebraic equations](@ref constraint_equations_coupling_constraints) (which are then incorporated during conversion to ODEs, SDEs, and steady state equations).
 - Models can be [coupled with events](@ref constraint_equations_events) that affect the system and its state during simulations.
diff --git a/docs/src/network_analysis/crnt.md b/docs/src/network_analysis/crnt.md
index c42b907e85..70500e504d 100644
--- a/docs/src/network_analysis/crnt.md
+++ b/docs/src/network_analysis/crnt.md
@@ -4,6 +4,11 @@ reaction networks can often have their steady-state properties known in advance,
 simply by analyzing the graph structure of the network. The subfield of chemistry
 and math studying this relationship is called [Chemical Reaction Network Theory](https://en.wikipedia.org/wiki/Chemical_reaction_network_theory).
 
+In this tutorial we introduce several of the Catalyst API functions for network
+analysis. A complete summary of the exported functions is given in the API
+section
+[`Network-Analysis-and-Representations`](https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Network-Analysis-and-Representations).
+
 Broadly, results from chemical reaction network theory relate a purely
 graph-structural property (e.g. deficiency) to dynamical properties of the reaction system
 (e.g. complex balance). We'll now illustrate some of the types of network properties that 
@@ -269,14 +274,10 @@ and,
  plot_complexes(subnets[2])
 ```
 
-<<<<<<< HEAD:docs/src/network_analysis/crnt.md
 ![subnetwork_2](../assets/complex_subnets2.svg)
 
 
-# [Deficiency of the network](@id network_analysis_structural_aspects_deficiency)
-=======
 #### [Deficiency of the network](@id network_analysis_structural_aspects_deficiency)
->>>>>>> origin:docs/src/model_creation/network_analysis.md
 A famous theorem in Chemical Reaction Network Theory, the Deficiency Zero
 Theorem [^1], allows us to use knowledge of the net stoichiometry matrix and the
 linkage classes of a *mass action* [RRE ODE system](@ref network_analysis_matrix_vector_representation) to draw conclusions about the
diff --git a/docs/src/network_analysis/odes.md b/docs/src/network_analysis/odes.md
index 1c7b718f10..5aa532c558 100644
--- a/docs/src/network_analysis/odes.md
+++ b/docs/src/network_analysis/odes.md
@@ -1,9 +1,10 @@
-# [Network Analysis in Catalyst](@id network_analysis)
+# [Decomposing the Reaction Network ODEs](@id network_analysis_odes)
 
-In this tutorial we introduce several of the Catalyst API functions for network
-analysis. A complete summary of the exported functions is given in the API
-section
-[`Network-Analysis-and-Representations`](https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Network-Analysis-and-Representations).
+In this tutorial we will discuss the specific mathematical 
+structure of the ODEs that arise from the mass-action dynamics
+of chemical reaction networks, and decompose them as a product
+of matrices that describe the network.  A complete summary of 
+the exported functions is given in the API section [`Network-Analysis-and-Representations`](https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Network-Analysis-and-Representations).
 
 Note, currently API functions for network analysis and conservation law analysis
 do not work with constant species (currently only generated by SBMLToolkit).

From ac5702bb9f642a930f8f09ca6cef318280f3d26c Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 28 Jan 2025 16:11:08 -0500
Subject: [PATCH 07/24] up

---
 docs/pages.jl                                   | 1 -
 docs/src/network_analysis/crnt.md               | 8 ++++----
 docs/src/network_analysis/network_properties.md | 2 +-
 docs/src/network_analysis/odes.md               | 2 +-
 src/network_analysis.jl                         | 4 ++--
 5 files changed, 8 insertions(+), 9 deletions(-)

diff --git a/docs/pages.jl b/docs/pages.jl
index 903b181b40..4fc7c8053b 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -16,7 +16,6 @@ pages = Any[
         "model_creation/model_file_loading_and_export.md",
         "model_creation/model_visualisation.md",
         "model_creation/reactionsystem_content_accessing.md",
-        "model_creation/network_analysis.md",
         "model_creation/chemistry_related_functionality.md",
         "Examples" => Any[
             "model_creation/examples/basic_CRN_library.md",
diff --git a/docs/src/network_analysis/crnt.md b/docs/src/network_analysis/crnt.md
index 70500e504d..2ed1e2c503 100644
--- a/docs/src/network_analysis/crnt.md
+++ b/docs/src/network_analysis/crnt.md
@@ -469,7 +469,7 @@ which we see is mass action and has deficiency zero, but is not weakly
 reversible. As such, we can conclude that for any choice of rate constants the
 RRE ODEs cannot have a positive equilibrium solution.
 
-```@docs
+```@docs canonical=false
 satisfiesdeficiencyzero
 ```
 
@@ -481,7 +481,7 @@ Very analogous to the deficiency zero theorem is the deficiency one theorem. The
 
 If these conditions are met, then the network will have at most one steady state in each stoichiometric compatibility class for any choice of rate constants and parameters.
 Unlike the deficiency zero theorem, networks obeying the deficiency one theorem are not guaranteed to have stable solutions.
-```@docs
+```@docs canonical=false
 satisfiesdeficiencyone
 ```
 
@@ -514,7 +514,7 @@ isdetailedbalanced(rn, rates)
 
 The reason that the deficiency zero theorem puts such strong restrictions on the steady state properties of the reaction network is because it implies that the reaction network will be complex balanced for any set of rate constants and parameters. The fact that this holds from a purely structural property of the graph, regardless of kinetics, is what makes it so useful. But in some cases it might be desirable to check complex balance on its own, as for higher deficiency networks.
 
-```@docs
+```@docs canonical=false
 iscomplexbalanced
 isdetailedbalanced
 ```
@@ -533,7 +533,7 @@ that can be checked relatively easily. One example is for deficiency one network
 
 This is the check provided by the API function `robustspecies(rn)`. More general concentration robustness analysis can be done using the [CatalystNetworkAnalysis](@ref) package.
 
-```@docs
+```@docs canonical=false
 robustspecies
 ```
 
diff --git a/docs/src/network_analysis/network_properties.md b/docs/src/network_analysis/network_properties.md
index 92b25c5055..cb787615c5 100644
--- a/docs/src/network_analysis/network_properties.md
+++ b/docs/src/network_analysis/network_properties.md
@@ -36,6 +36,6 @@ Catalyst.reset_networkproperties!(rn)
 Network property functions will then recalculate their associated properties and
 cache the new values the next time they are called.
 
-```@docs
+```@docs canonical=false
 NetworkProperties
 ```
diff --git a/docs/src/network_analysis/odes.md b/docs/src/network_analysis/odes.md
index 5aa532c558..3e0b12e2fb 100644
--- a/docs/src/network_analysis/odes.md
+++ b/docs/src/network_analysis/odes.md
@@ -342,7 +342,7 @@ We have that:
 - ``A_k`` is the `laplacianmat`
 - ``\Phi`` is the `massactionvector`
 
-```@docs
+```@docs canonical=false
 netstoichmat
 complexstoichmat
 incidencemat
diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index 2800f36ddc..d5aa9e0bb8 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -1222,7 +1222,7 @@ end
 """
     satisfiesdeficiencyone(rn::ReactionSystem)
 
-    Check if a reaction network obeys the conditions of the deficiency one theorem, which ensures that there is only one equilibrium for every positive stoichiometric compatibility class.
+Check if a reaction network obeys the conditions of the deficiency one theorem, which ensures that there is only one equilibrium for every positive stoichiometric compatibility class.
 """
 function satisfiesdeficiencyone(rn::ReactionSystem)
     all(r -> ismassaction(r, rn), reactions(rn)) ||
@@ -1244,7 +1244,7 @@ end
 """
     satisfiesdeficiencyzero(rn::ReactionSystem)
 
-    Check if a reaction network obeys the conditions of the deficiency zero theorem, which ensures that there is only one equilibrium for every positive stoichiometric compatibility class, this equilibrium is asymptotically stable, and this equilibrium is complex balanced.
+Check if a reaction network obeys the conditions of the deficiency zero theorem, which ensures that there is only one equilibrium for every positive stoichiometric compatibility class, this equilibrium is asymptotically stable, and this equilibrium is complex balanced.
 """
 function satisfiesdeficiencyzero(rn::ReactionSystem)
     all(r -> ismassaction(r, rn), reactions(rn)) ||

From a5b3d3fab00295757648043db6d71cf44294efc1 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 28 Jan 2025 16:15:38 -0500
Subject: [PATCH 08/24] up

---
 docs/src/network_analysis/crnt.md               | 8 ++++----
 docs/src/network_analysis/network_properties.md | 2 +-
 docs/src/network_analysis/odes.md               | 2 +-
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/docs/src/network_analysis/crnt.md b/docs/src/network_analysis/crnt.md
index 2ed1e2c503..8789f30750 100644
--- a/docs/src/network_analysis/crnt.md
+++ b/docs/src/network_analysis/crnt.md
@@ -469,7 +469,7 @@ which we see is mass action and has deficiency zero, but is not weakly
 reversible. As such, we can conclude that for any choice of rate constants the
 RRE ODEs cannot have a positive equilibrium solution.
 
-```@docs canonical=false
+```@docs; canonical=false
 satisfiesdeficiencyzero
 ```
 
@@ -481,7 +481,7 @@ Very analogous to the deficiency zero theorem is the deficiency one theorem. The
 
 If these conditions are met, then the network will have at most one steady state in each stoichiometric compatibility class for any choice of rate constants and parameters.
 Unlike the deficiency zero theorem, networks obeying the deficiency one theorem are not guaranteed to have stable solutions.
-```@docs canonical=false
+```@docs; canonical=false
 satisfiesdeficiencyone
 ```
 
@@ -514,7 +514,7 @@ isdetailedbalanced(rn, rates)
 
 The reason that the deficiency zero theorem puts such strong restrictions on the steady state properties of the reaction network is because it implies that the reaction network will be complex balanced for any set of rate constants and parameters. The fact that this holds from a purely structural property of the graph, regardless of kinetics, is what makes it so useful. But in some cases it might be desirable to check complex balance on its own, as for higher deficiency networks.
 
-```@docs canonical=false
+```@docs; canonical=false
 iscomplexbalanced
 isdetailedbalanced
 ```
@@ -533,7 +533,7 @@ that can be checked relatively easily. One example is for deficiency one network
 
 This is the check provided by the API function `robustspecies(rn)`. More general concentration robustness analysis can be done using the [CatalystNetworkAnalysis](@ref) package.
 
-```@docs canonical=false
+```@docs; canonical=false
 robustspecies
 ```
 
diff --git a/docs/src/network_analysis/network_properties.md b/docs/src/network_analysis/network_properties.md
index cb787615c5..31586e39da 100644
--- a/docs/src/network_analysis/network_properties.md
+++ b/docs/src/network_analysis/network_properties.md
@@ -36,6 +36,6 @@ Catalyst.reset_networkproperties!(rn)
 Network property functions will then recalculate their associated properties and
 cache the new values the next time they are called.
 
-```@docs canonical=false
+```@docs; canonical=false
 NetworkProperties
 ```
diff --git a/docs/src/network_analysis/odes.md b/docs/src/network_analysis/odes.md
index 3e0b12e2fb..28ecbd79fc 100644
--- a/docs/src/network_analysis/odes.md
+++ b/docs/src/network_analysis/odes.md
@@ -342,7 +342,7 @@ We have that:
 - ``A_k`` is the `laplacianmat`
 - ``\Phi`` is the `massactionvector`
 
-```@docs canonical=false
+```@docs; canonical=false
 netstoichmat
 complexstoichmat
 incidencemat

From bde1c35195762dfb3b7b8e702b342d6fe79ca6c4 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 4 Feb 2025 20:14:35 -0500
Subject: [PATCH 09/24] init multidocumenter

---
 docs/make.jl                                  | 28 +++++++++++++++++++
 .../src/network_analysis/advanced_analysis.md |  3 +-
 docs/src/network_analysis/crnt.md             | 17 +++++------
 3 files changed, 39 insertions(+), 9 deletions(-)

diff --git a/docs/make.jl b/docs/make.jl
index f732463793..5516a9756a 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -2,6 +2,7 @@ using Documenter
 using Catalyst, ModelingToolkit
 # Add packages for plotting
 using GraphMakie, CairoMakie
+using MultiDocumenter
 
 docpath = Base.source_dir()
 assetpath = joinpath(docpath, "src", "assets")
@@ -48,5 +49,32 @@ makedocs(sitename = "Catalyst.jl",
     pagesonly = true,
     warnonly = [:missing_docs])
 
+clonedir = mktempdir()
+
+docs = [MultiDocumenter.MultiDocRef(
+        upstream = joinpath(@__DIR__, "build"),
+        path = "docs",
+        name = "Catalyst",
+        fix_canonical_url = false,
+        ),
+        MultiDocumenter.MultiDocRef(
+        upstream = joinpath(clonedir, "build"),
+        path = "NetworkAnalysis",
+        name = "CatalystNetworkAnalysis",
+        giturl = "https://github.com/SciML/CatalystNetworkAnalysis.jl",
+        )]
+
 deploydocs(repo = "github.com/SciML/Catalyst.jl.git";
     push_preview = true)
+
+outpath = joinpath(@__DIR__, "build")
+MultiDocumenter.make(outpath, docs;
+                     assets_dir = "docs/src/assets",
+                     search_engine = MultiDocumenter.SearchConfig(index_versions = [
+                                                                      "stable",
+                                                                  ],
+                                                                  engine = MultiDocumenter.FlexSearch),
+                     brand_image = MultiDocumenter.BrandImage("https://sciml.ai",
+                                                              joinpath("assets",
+                                                                       "logo.png")))
+
diff --git a/docs/src/network_analysis/advanced_analysis.md b/docs/src/network_analysis/advanced_analysis.md
index 3bf024ee1e..897673a95d 100644
--- a/docs/src/network_analysis/advanced_analysis.md
+++ b/docs/src/network_analysis/advanced_analysis.md
@@ -1,3 +1,4 @@
 # [Advanced Analysis using CatalystNetworkAnalysis](@ref CatalystNetworkAnalysis)
 
-For more advanced chemical reaction network analysis, there is a special CatalystNetworkAnalysis package. This package has functionality for a wider range of analysis for chemical reaction networks, including the following:
+For more advanced chemical reaction network analysis, there is a special CatalystNetworkAnalysis package. This package has functionality for a wider range of analysis for chemical reaction networks. For more information, please see the [CatalystNetworkAnalysis]() section accessible from the header of this page.
+
diff --git a/docs/src/network_analysis/crnt.md b/docs/src/network_analysis/crnt.md
index 8789f30750..13fd4cc342 100644
--- a/docs/src/network_analysis/crnt.md
+++ b/docs/src/network_analysis/crnt.md
@@ -1,4 +1,5 @@
 # [Chemical Reaction Network Theory](@id network_analysis_structural_aspects)
+
 The systems of ODEs or stochastic chemical kinetics models that arise from chemical
 reaction networks can often have their steady-state properties known in advance,
 simply by analyzing the graph structure of the network. The subfield of chemistry
@@ -110,7 +111,7 @@ Let's check these are equal symbolically
 isequal(odes, odes2)
 ```
 
-## [Reaction complex representation](@id network_analysis_reaction_complexes)
+### [Reaction complex representation](@id network_analysis_reaction_complexes)
 We now introduce a further decomposition of the RRE ODEs, which has been used to
 facilitate analysis of a variety of reaction network properties. Consider a simple
 reaction system like
@@ -213,7 +214,7 @@ substrate complexes into product complexes where the rate is just a number or
 parameter, and red arrows indicate the conversion of substrate complexes into
 product complexes where the rate is an expression involving chemical species.
 
-## [Aspects of reaction network structure](@id network_analysis_structural_aspects)
+### [Aspects of reaction network structure](@id network_analysis_structural_aspects)
 The reaction complex representation can be exploited via [Chemical Reaction
 Network Theory](https://en.wikipedia.org/wiki/Chemical_reaction_network_theory)
 to provide insight into possible steady state and time-dependent properties of
@@ -277,7 +278,7 @@ and,
 ![subnetwork_2](../assets/complex_subnets2.svg)
 
 
-#### [Deficiency of the network](@id network_analysis_structural_aspects_deficiency)
+### [Deficiency of the network](@id network_analysis_structural_aspects_deficiency)
 A famous theorem in Chemical Reaction Network Theory, the Deficiency Zero
 Theorem [^1], allows us to use knowledge of the net stoichiometry matrix and the
 linkage classes of a *mass action* [RRE ODE system](@ref network_analysis_matrix_vector_representation) to draw conclusions about the
@@ -332,7 +333,7 @@ Quoting Feinberg [^1]
 > linkage classes will allow.
 
 
-# [Reversibility of the network](@id network_analysis_structural_aspects_reversibility)
+### [Reversibility of the network](@id network_analysis_structural_aspects_reversibility)
 A reaction network is *reversible* if the "arrows" of the reactions are
 symmetric so that every reaction is accompanied by its reverse reaction.
 Catalyst's API provides the [`isreversible`](@ref) function to determine whether
@@ -379,7 +380,7 @@ isweaklyreversible(rn, subnets)
 Every reversible network is also weakly reversible, but not every weakly
 reversible network is reversible.
 
-# [Deficiency Zero Theorem](@id network_analysis_structural_aspects_deficiency_zero_theorem)
+### [Deficiency Zero Theorem](@id network_analysis_structural_aspects_deficiency_zero_theorem)
 Knowing the deficiency and weak reversibility of a mass action chemical reaction
 network ODE model allows us to make inferences about the corresponding
 steady state behavior. Before illustrating how this works for one example, we
@@ -473,7 +474,7 @@ RRE ODEs cannot have a positive equilibrium solution.
 satisfiesdeficiencyzero
 ```
 
-# Deficiency One Theorem
+### Deficiency One Theorem
 Very analogous to the deficiency zero theorem is the deficiency one theorem. The deficiency one theorem applies to a network with the following properties: 
 1. The deficiency of each *linkage class* of the network is at most 1, 
 2. The sum of the linkage class deficiencies is the total deficiency of the network, and
@@ -485,7 +486,7 @@ Unlike the deficiency zero theorem, networks obeying the deficiency one theorem
 satisfiesdeficiencyone
 ```
 
-# [Complex and Detailed Balance](@id network_analysis_complex_and_detailed_balance)
+### [Complex and Detailed Balance](@id network_analysis_complex_and_detailed_balance)
 A reaction network's steady state is **complex-balanced** if the total production of each *complex* is zero at the steady state. A reaction network's steady state is **detailed balanced** if every reaction is balanced by its reverse reaction at the steady-state (this corresponds to the usual notion of chemical equilibrium; note that this requires every reaction be reversible). 
 
 Note that detailed balance at a given steady state implies complex balance for that steady state, i.e. detailed balance is a stronger property than complex balance.
@@ -519,7 +520,7 @@ iscomplexbalanced
 isdetailedbalanced
 ```
 
-# [Concentration Robustness](@id network_analysis_concentration_robustness)
+### [Concentration Robustness](@id network_analysis_concentration_robustness)
 Certain reaction networks have species that do not change their concentration,
 regardless of the system is perturbed to a different stoichiometric compatibility
 class. This is a very useful property to have in biological contexts, where it might

From 959e79ab589d605f9badc1fd6d86f8d7e4e40f74 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 4 Feb 2025 20:21:56 -0500
Subject: [PATCH 10/24] up

---
 docs/Project.toml                 | 1 +
 docs/src/network_analysis/crnt.md | 5 -----
 2 files changed, 1 insertion(+), 5 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index b4ef01b543..c3bdf41270 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -18,6 +18,7 @@ Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+MultiDocumenter = "87ed4bf0-c935-4a67-83c3-2a03bee4197c"
 NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
diff --git a/docs/src/network_analysis/crnt.md b/docs/src/network_analysis/crnt.md
index 13fd4cc342..c44a1c401a 100644
--- a/docs/src/network_analysis/crnt.md
+++ b/docs/src/network_analysis/crnt.md
@@ -221,7 +221,6 @@ to provide insight into possible steady state and time-dependent properties of
 RRE ODE models and stochastic chemical kinetics models. We'll now illustrate
 some of the types of network properties that Catalyst can determine, using the
 reaction complex representation in these calculations.
->>>>>>> origin:docs/src/model_creation/network_analysis.md
 
 Consider the following reaction network.
 ```@example s1
@@ -499,10 +498,6 @@ rates = Dict([:k1 => 2.4, :k2 => 4., :k3 => 10., :k4 => 5.5, :k5 => 0.4])
 iscomplexbalanced(rn, rates)
 ```
 
-Complex balance obtains for some sets of rates but not others: 
-```@example s1
-```
-
 We can do a similar check for detailed balance. Let us make the reaction network 
 ```@example s1
 rn1 = @reaction_network begin

From 739971ce61eae2a09a027b1562f64d8a890b918a Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 4 Feb 2025 21:39:06 -0500
Subject: [PATCH 11/24] up

---
 docs/src/assets/Project.toml      | 91 +++++++++++++++++++++++++++++++
 docs/src/network_analysis/odes.md |  3 +-
 2 files changed, 93 insertions(+), 1 deletion(-)
 create mode 100644 docs/src/assets/Project.toml

diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml
new file mode 100644
index 0000000000..c3bdf41270
--- /dev/null
+++ b/docs/src/assets/Project.toml
@@ -0,0 +1,91 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634"
+GlobalSensitivity = "af5da776-676b-467e-8baf-acd8249e4f0f"
+GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
+IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
+JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
+Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+MultiDocumenter = "87ed4bf0-c935-4a67-83c3-2a03bee4197c"
+NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
+OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b"
+OptimizationEvolutionary = "cb963754-43f6-435e-8d4b-99009ff27753"
+OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
+OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
+OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
+OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8"
+OrdinaryDiffEqDefault = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
+OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
+PEtab = "48d54b35-e43e-4a66-a5a1-dde6b987cf69"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
+StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+
+[compat]
+BenchmarkTools = "1.5"
+BifurcationKit = "0.4.4"
+CairoMakie = "0.12, 0.13"
+Catalyst = "14.4"
+DataFrames = "1.6"
+DiffEqBase = "6.159.0"
+Distributions = "0.25"
+Documenter = "1.4.1"
+DynamicalSystems = "3.3"
+GlobalSensitivity = "2.6"
+GraphMakie = "0.5"
+Graphs = "1.11.1"
+HomotopyContinuation = "2.9"
+IncompleteLU = "0.2"
+JumpProcesses = "9.13.2"
+Latexify = "0.16.5"
+LinearSolve = "2.30"
+ModelingToolkit = "< 9.60"
+NetworkLayout = "0.4"
+NonlinearSolve = "3.12, 4"
+Optim = "1.9"
+Optimization = "4"
+OptimizationBBO = "0.4"
+OptimizationEvolutionary = "0.4"
+OptimizationNLopt = "0.3"
+OptimizationOptimJL = "0.4"
+OptimizationOptimisers = "0.3"
+OrdinaryDiffEqBDF = "1"
+OrdinaryDiffEqDefault = "1"
+OrdinaryDiffEqRosenbrock = "1"
+OrdinaryDiffEqSDIRK = "1"
+OrdinaryDiffEqTsit5 = "1"
+OrdinaryDiffEqVerner = "1"
+PEtab = "3.5"
+Plots = "1.40"
+QuasiMonteCarlo = "0.3"
+SciMLBase = "2.46"
+SciMLSensitivity = "7.60"
+SpecialFunctions = "2.4"
+StaticArrays = "1.9"
+SteadyStateDiffEq = "2.2"
+StochasticDiffEq = "6.65"
+Symbolics = "6.22"
diff --git a/docs/src/network_analysis/odes.md b/docs/src/network_analysis/odes.md
index 28ecbd79fc..ed715635f1 100644
--- a/docs/src/network_analysis/odes.md
+++ b/docs/src/network_analysis/odes.md
@@ -4,7 +4,8 @@ In this tutorial we will discuss the specific mathematical
 structure of the ODEs that arise from the mass-action dynamics
 of chemical reaction networks, and decompose them as a product
 of matrices that describe the network.  A complete summary of 
-the exported functions is given in the API section [`Network-Analysis-and-Representations`](https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Network-Analysis-and-Representations).
+the exported functions is given in the API section [`Network-Analysis-and-Representations`](https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Network-Analysis-and-Representations). See Feinberg's *Foundations of Chemical Reaction
+Network Theory*[^1] for more discussion about the concepts on this page.
 
 Note, currently API functions for network analysis and conservation law analysis
 do not work with constant species (currently only generated by SBMLToolkit).

From fbf8c63720a4a1ec7f452505034cd026405a0464 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Tue, 4 Feb 2025 21:39:32 -0500
Subject: [PATCH 12/24] up

---
 docs/src/assets/Project.toml | 91 ------------------------------------
 1 file changed, 91 deletions(-)
 delete mode 100644 docs/src/assets/Project.toml

diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml
deleted file mode 100644
index c3bdf41270..0000000000
--- a/docs/src/assets/Project.toml
+++ /dev/null
@@ -1,91 +0,0 @@
-[deps]
-BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
-Catalyst = "479239e8-5488-4da2-87a7-35f2df7eef83"
-DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
-DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634"
-GlobalSensitivity = "af5da776-676b-467e-8baf-acd8249e4f0f"
-GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
-Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
-HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
-IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
-JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
-Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
-LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
-ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
-MultiDocumenter = "87ed4bf0-c935-4a67-83c3-2a03bee4197c"
-NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a"
-NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-Optim = "429524aa-4258-5aef-a3af-852621145aeb"
-Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
-OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b"
-OptimizationEvolutionary = "cb963754-43f6-435e-8d4b-99009ff27753"
-OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
-OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
-OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
-OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8"
-OrdinaryDiffEqDefault = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
-OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
-OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
-OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
-OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
-PEtab = "48d54b35-e43e-4a66-a5a1-dde6b987cf69"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
-SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
-SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
-StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
-Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-
-[compat]
-BenchmarkTools = "1.5"
-BifurcationKit = "0.4.4"
-CairoMakie = "0.12, 0.13"
-Catalyst = "14.4"
-DataFrames = "1.6"
-DiffEqBase = "6.159.0"
-Distributions = "0.25"
-Documenter = "1.4.1"
-DynamicalSystems = "3.3"
-GlobalSensitivity = "2.6"
-GraphMakie = "0.5"
-Graphs = "1.11.1"
-HomotopyContinuation = "2.9"
-IncompleteLU = "0.2"
-JumpProcesses = "9.13.2"
-Latexify = "0.16.5"
-LinearSolve = "2.30"
-ModelingToolkit = "< 9.60"
-NetworkLayout = "0.4"
-NonlinearSolve = "3.12, 4"
-Optim = "1.9"
-Optimization = "4"
-OptimizationBBO = "0.4"
-OptimizationEvolutionary = "0.4"
-OptimizationNLopt = "0.3"
-OptimizationOptimJL = "0.4"
-OptimizationOptimisers = "0.3"
-OrdinaryDiffEqBDF = "1"
-OrdinaryDiffEqDefault = "1"
-OrdinaryDiffEqRosenbrock = "1"
-OrdinaryDiffEqSDIRK = "1"
-OrdinaryDiffEqTsit5 = "1"
-OrdinaryDiffEqVerner = "1"
-PEtab = "3.5"
-Plots = "1.40"
-QuasiMonteCarlo = "0.3"
-SciMLBase = "2.46"
-SciMLSensitivity = "7.60"
-SpecialFunctions = "2.4"
-StaticArrays = "1.9"
-SteadyStateDiffEq = "2.2"
-StochasticDiffEq = "6.65"
-Symbolics = "6.22"

From 64459286640c8bbecaa894cc385ef4046811b540 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 5 Feb 2025 09:52:19 -0500
Subject: [PATCH 13/24] fix references

---
 .../src/network_analysis/advanced_analysis.md |   4 +-
 docs/src/network_analysis/crnt.md             | 211 +-----------------
 docs/src/network_analysis/odes.md             |  25 +--
 3 files changed, 15 insertions(+), 225 deletions(-)

diff --git a/docs/src/network_analysis/advanced_analysis.md b/docs/src/network_analysis/advanced_analysis.md
index 897673a95d..908136695e 100644
--- a/docs/src/network_analysis/advanced_analysis.md
+++ b/docs/src/network_analysis/advanced_analysis.md
@@ -1,4 +1,4 @@
-# [Advanced Analysis using CatalystNetworkAnalysis](@ref CatalystNetworkAnalysis)
+# [Advanced Analysis using CatalystNetworkAnalysis](@id catalyst_network_analysis)
 
-For more advanced chemical reaction network analysis, there is a special CatalystNetworkAnalysis package. This package has functionality for a wider range of analysis for chemical reaction networks. For more information, please see the [CatalystNetworkAnalysis]() section accessible from the header of this page.
+For more advanced chemical reaction network analysis, there is a special CatalystNetworkAnalysis package. This package has functionality for a wider range of analysis for chemical reaction networks. For more information, please see the CatalystNetworkAnalysis section accessible from the header of this page.
 
diff --git a/docs/src/network_analysis/crnt.md b/docs/src/network_analysis/crnt.md
index c44a1c401a..ef93c2bcb6 100644
--- a/docs/src/network_analysis/crnt.md
+++ b/docs/src/network_analysis/crnt.md
@@ -18,211 +18,7 @@ Catalyst can determine, using the [reaction complex representation](@ref network
 Note, currently API functions for network analysis and conservation law analysis
 do not work with constant species (currently only generated by SBMLToolkit).
 
-## [Network representation of the Repressilator `ReactionSystem`](@id network_analysis_repressilator_representation)
-We first load Catalyst and construct our model of the repressilator
-```@example s1
-using Catalyst
-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
-```
-In the [Introduction to Catalyst](@ref introduction_to_catalyst)
-tutorial we showed how the above network could be visualized as a
-species-reaction graph. There, species are represented by the nodes of the graph
-and edges show the reactions in which a given species is a substrate or product.
-```@example s1
-using Catalyst, GraphMakie, NetworkLayout, CairoMakie
-g = plot_network(repressilator)
-```
-
-We also showed in the [Introduction to Catalyst](@ref introduction_to_catalyst) tutorial that
-the reaction rate equation ODE model for the repressilator is
-```math
-\begin{aligned}
-\frac{dm_1(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_3}\left( t \right) \right)^{n}} - \delta {m_1}\left( t \right) + \gamma \\
-\frac{dm_2(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_1}\left( t \right) \right)^{n}} - \delta {m_2}\left( t \right) + \gamma \\
-\frac{dm_3(t)}{dt} =& \frac{\alpha K^{n}}{K^{n} + \left( {P_2}\left( t \right) \right)^{n}} - \delta {m_3}\left( t \right) + \gamma \\
-\frac{dP_1(t)}{dt} =& \beta {m_1}\left( t \right) - \mu {P_1}\left( t \right) \\
-\frac{dP_2(t)}{dt} =& \beta {m_2}\left( t \right) - \mu {P_2}\left( t \right) \\
-\frac{dP_3(t)}{dt} =& \beta {m_3}\left( t \right) - \mu {P_3}\left( t \right)
-\end{aligned}
-```
-
-## [Matrix-vector reaction rate equation representation](@id network_analysis_matrix_vector_representation)
-In general, reaction rate equation (RRE) ODE models for chemical reaction networks can
-be represented as a first-order system of ODEs in a compact matrix-vector notation. Suppose
-we have a reaction network with ``K`` reactions and ``M`` species, labelled by the state vector
-```math
-\mathbf{x}(t) = \begin{pmatrix} x_1(t) \\ \vdots \\ x_M(t)) \end{pmatrix}.
-```
-For the repressilator, ``\mathbf{x}(t)`` is just
-```@example s1
-x = species(repressilator)
-```
-The RRE ODEs satisfied by $\mathbf{x}(t)$ are then
-```math
-\frac{d\mathbf{x}}{dt} = N \mathbf{v}(\mathbf{x}(t),t),
-```
-where ``N`` is a constant ``M`` by ``K`` matrix with ``N_{m k}`` the net
-stoichiometric coefficient of species ``m`` in reaction ``k``.
-``\mathbf{v}(\mathbf{x}(t),t)`` is the rate law vector, with
-``v_k(\mathbf{x}(t),t)`` the rate law for the ``k``th reaction. For example,
-for the first reaction of the repressilator above, the rate law is
-```math
-v_1(\mathbf{x}(t),t) = \frac{\alpha K^{n}}{K^{n} + \left( P_3(t) \right)^{n}}.
-```
-We can calculate each of these in Catalyst via
-```@example s1
-N = netstoichmat(repressilator)
-```
-and by using the [`oderatelaw`](@ref) function
-```@example s1
-rxs = reactions(repressilator)
-ν = oderatelaw.(rxs)
-```
-Note, as [`oderatelaw`](@ref) takes just one reaction as input we use
-broadcasting to apply it to each element of `rxs`.
-
-Let's check that this really gives the same ODEs as Catalyst. Here is what Catalyst
-generates by converting to an `ODESystem`
-```@example s1
-osys = convert(ODESystem, repressilator)
-
-# for display purposes we just pull out the right side of the equations
-odes = [eq.rhs for eq in equations(osys)]
-```
-whereas our matrix-vector representation gives
-```@example s1
-odes2 = N * ν
-```
-Let's check these are equal symbolically
-```@example s1
-isequal(odes, odes2)
-```
-
-### [Reaction complex representation](@id network_analysis_reaction_complexes)
-We now introduce a further decomposition of the RRE ODEs, which has been used to
-facilitate analysis of a variety of reaction network properties. Consider a simple
-reaction system like
-```@example s1
-rn = @reaction_network begin
- k*A, 2*A + 3*B --> A + 2*C + D
- b, C + D --> 2*A + 3*B
-end
-```
-We can think of the first reaction as converting the *reaction complex*,
-``2A+3B`` to the complex ``A+2C+D`` with rate ``kA``. Suppose we order our
-species the same way as Catalyst does, i.e.
-```math
-\begin{pmatrix}
-x_1(t)\\
-x_2(t)\\
-x_3(t)\\
-x_4(t)
-\end{pmatrix} =
-\begin{pmatrix}
-A(t)\\
-B(t)\\
-C(t)\\
-D(t)
-\end{pmatrix},
-```
-which should be the same as
-```@example s1
-species(rn)
-```
-We can describe a given reaction complex by the stoichiometric coefficients of
-each species within the complex. For the reactions in `rn` these vectors would
-be
-```math
-\begin{align*}
-2A+3B = \begin{pmatrix}
-2\\
-3\\
-0\\
-0
-\end{pmatrix}, &&
-A+2C+D = \begin{pmatrix}
-1\\
-0\\
-2\\
-1
-\end{pmatrix},
- &&
-C+D = \begin{pmatrix}
-0\\
-0\\
-1\\
-1
-\end{pmatrix}
-\end{align*}
-```
-Catalyst can calculate these representations as the columns of the complex
-stoichiometry matrix,
-```@example s1
-Z = complexstoichmat(rn)
-```
-If we have ``C`` complexes, ``Z`` is a ``M`` by ``C`` matrix with ``Z_{m c}``
-giving the stoichiometric coefficient of species ``m`` within complex ``c``.
-
-We can use this representation to provide another representation of the RRE
-ODEs. The net stoichiometry matrix can be factored as ``N = Z B``, where ``B``
-is called the incidence matrix of the reaction network,
-```@example s1
-B = incidencemat(rn)
-```
-Here ``B`` is a ``C`` by ``K`` matrix with ``B_{c k} = 1`` if complex ``c``
-appears as a product of reaction ``k``, and ``B_{c k} = -1`` if complex ``c`` is a
-substrate of reaction ``k``.
-
-Using our decomposition of ``N``, the RRE ODEs become
-```math
-\frac{dx}{dt} = Z B \mathbf{v}(\mathbf{x}(t),t).
-```
-Let's verify that ``N = Z B``,
-```@example s1
-N = netstoichmat(rn)
-N == Z*B
-```
-
-Reaction complexes give an alternative way to visualize a reaction network
-graph. Catalyst's [`plot_complexes`](@ref) command will calculate the complexes of
-a network and then show how they are related. For example,
-```@example s1
-using Catalyst, GraphMakie, NetworkLayout, CairoMakie
-plot_complexes(rn)
-```
-
-while for the repressilator we find
-```@example s1
-plot_complexes(repressilator)
-```
-
-Here ∅ represents the empty complex, black arrows show reactions converting
-substrate complexes into product complexes where the rate is just a number or
-parameter, and red arrows indicate the conversion of substrate complexes into
-product complexes where the rate is an expression involving chemical species.
-
-### [Aspects of reaction network structure](@id network_analysis_structural_aspects)
-The reaction complex representation can be exploited via [Chemical Reaction
-Network Theory](https://en.wikipedia.org/wiki/Chemical_reaction_network_theory)
-to provide insight into possible steady state and time-dependent properties of
-RRE ODE models and stochastic chemical kinetics models. We'll now illustrate
-some of the types of network properties that Catalyst can determine, using the
-reaction complex representation in these calculations.
-
-Consider the following reaction network.
+Throughout this seciton, we will consider the [reaction complex representation](@ref network_analysis_reaction_complexes) of the following reaction network.
 ```@example s1
 rn = @reaction_network begin
     (k1,k2), A + B <--> C
@@ -238,8 +34,7 @@ with graph
 plot_complexes(rn)
 ```
 
-
-# [Linkage classes and sub-networks of the reaction network](@id network_analysis_structural_aspects_linkage)
+### [Linkage classes and sub-networks of the reaction network](@id network_analysis_structural_aspects_linkage)
 The preceding reaction complex graph shows that `rn` is composed of two
 disconnected sub-graphs, one containing the complexes ``A+B``, ``C``, ``D+E``, and
 ``F``, the other containing the complexes ``2A``, ``B + G``, and ``H``. These sets,
@@ -527,7 +322,7 @@ that can be checked relatively easily. One example is for deficiency one network
 
 **Theorem (a sufficient condition for concentration robustness for deficiency one networks)**: If there are two *non-terminal* reaction complexes that differ only in species ``s``, then the system is absolutely concentration robust with respect to ``s``. 
 
-This is the check provided by the API function `robustspecies(rn)`. More general concentration robustness analysis can be done using the [CatalystNetworkAnalysis](@ref) package.
+This is the check provided by the API function `robustspecies(rn)`. More general concentration robustness analysis can be done using the [CatalystNetworkAnalysis](@ref catalyst_network_analysis) package.
 
 ```@docs; canonical=false
 robustspecies
diff --git a/docs/src/network_analysis/odes.md b/docs/src/network_analysis/odes.md
index ed715635f1..f3b2247e18 100644
--- a/docs/src/network_analysis/odes.md
+++ b/docs/src/network_analysis/odes.md
@@ -109,12 +109,12 @@ facilitate analysis of a variety of reaction network properties. Consider a simp
 reaction system like
 ```@example s1
 rn = @reaction_network begin
- k*A, 2*A + 3*B --> A + 2*C + D
+ k, 2*A + 3*B --> A + 2*C + D
  b, C + D --> 2*A + 3*B
 end
 ```
 We can think of the first reaction as converting the *reaction complex*,
-``2A+3B`` to the complex ``A+2C+D`` with rate ``kA``. Suppose we order our
+``2A+3B`` to the complex ``A+2C+D`` with rate ``k``. Suppose we order our
 species the same way as Catalyst does, i.e.
 ```math
 \begin{pmatrix}
@@ -189,22 +189,17 @@ N == Z*B
 ```
 
 Reaction complexes give an alternative way to visualize a reaction network
-graph. Catalyst's [`complexgraph`](@ref) command will calculate the complexes of
-a network and then show how they are related. For example,
-```julia
-complexgraph(rn)
+graph. Catalyst's [`plot_complexes`](@ref) command will calculate the complexes of
+a network and then show how they are related. For example, we can run
+```@example s1
+plot_complexes(rn)
 ```
-gives
-
-![Simple example complex graph](../assets/simple_complexgraph.svg)
 
 while for the repressilator we find
 ```julia
-complexgraph(repressilator)
+plot_complexes(repressilator)
 ```
 
-![Repressilator complex](../assets/repressilator_complexgraph.svg)
-
 Here ∅ represents the empty complex, black arrows show reactions converting
 substrate complexes into product complexes where the rate is just a number or
 parameter, and red arrows indicate the conversion of substrate complexes into
@@ -228,7 +223,7 @@ where $\mathbf{x}^y = \prod_s x_s^{y_s}$, the mass-action product of the complex
 Φ = massactionvector(rn)
 ```
 
-An important thing to note is this function assumes [combinatoric ratelaws](@introduction_to_catalyst_ratelaws), meaning that mass-action products will get rescaled by factorial factors. For instance, note that the mass-action product for the complex `2A + 3B` has a factor of 1/12, corresponding to 1/(2! 3!). This option can be turned off with `combinatoric_ratelaws = false`.
+An important thing to note is this function assumes [combinatoric ratelaws](@ref introduction_to_catalyst_ratelaws), meaning that mass-action products will get rescaled by factorial factors. For instance, note that the mass-action product for the complex `2A + 3B` has a factor of 1/12, corresponding to 1/(2! 3!). This option can be turned off with `combinatoric_ratelaws = false`.
 
 ```@example s1
 Φ_2 = massactionvector(rn; combinatoric_ratelaws = false)
@@ -276,7 +271,7 @@ massactionvector(rn, concmap)
 
 `fluxmat` and `laplacianmat` will return numeric matrices if a set of rate constants and other aprameters are supplied the same way.
 ```@example s1
-parammap = Dict([:k => 12., b => 8.])
+parammap = Dict([:k => 12., :b => 8.])
 fluxmat(rn, parammap)
 ```
 
@@ -332,7 +327,7 @@ Recall that we may write the net stoichiometry matrix ``N = YB``.
 
 [Conservation laws](@ref conservation_laws) arise as left null eigenvectors of the net stoichiometry matrix ``N``, and cycles arise as right null eigenvectors of the stoichiometry matrix. A cycle may be understood as a sequence of reactions that leaves the overall species composition unchanged. These do not necessarily have to correspond to actual cycles in the graph.
 
-[Complex balance](@ref complex_and_detailed_balance) can be compactly formulated as the following: a set of steady state reaction fluxes is complex-balanced if it is in the nullspace of the incidence matrix ``B``.
+[Complex balance](@ref network_analysis_complex_and_detailed_balance) can be compactly formulated as the following: a set of steady state reaction fluxes is complex-balanced if it is in the nullspace of the incidence matrix ``B``.
 
 # API Section for matrices and vectors
 We have that: 

From e9edce1ffca30eff486996110d41ac57e045ad62 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 5 Feb 2025 11:03:13 -0500
Subject: [PATCH 14/24] load Catalyst

---
 docs/src/network_analysis/crnt.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/src/network_analysis/crnt.md b/docs/src/network_analysis/crnt.md
index ef93c2bcb6..4b3715fe55 100644
--- a/docs/src/network_analysis/crnt.md
+++ b/docs/src/network_analysis/crnt.md
@@ -20,6 +20,7 @@ do not work with constant species (currently only generated by SBMLToolkit).
 
 Throughout this seciton, we will consider the [reaction complex representation](@ref network_analysis_reaction_complexes) of the following reaction network.
 ```@example s1
+using Catalyst
 rn = @reaction_network begin
     (k1,k2), A + B <--> C
     k3, C --> D+E

From 0dcc7a5bd6ae1d3a151f392cafeb39540b6d72e2 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 5 Feb 2025 12:29:19 -0500
Subject: [PATCH 15/24] fix

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

diff --git a/docs/src/network_analysis/odes.md b/docs/src/network_analysis/odes.md
index f3b2247e18..a93f9316ef 100644
--- a/docs/src/network_analysis/odes.md
+++ b/docs/src/network_analysis/odes.md
@@ -284,7 +284,7 @@ In some cases it might be useful to generate the function defining the system of
 
 Let's build the full symbolic function corresponding to our ODE system. `build_function` will return two expressions, one for a function that outputs a new vector for the result, and one for a function that modifies the input in-place. Either expression can then be evaluated to return a Julia function.
 ```@example s1
-parammap = Dict([:k => 12., b => 8.])
+parammap = Dict([:k => 12., :b => 8.])
 K = fluxmat(rn, parammap)
 odes = N * K * Φ
 f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn))

From 217570b5cbdf57dd7848a59e5d1bcb1042a41eb6 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 5 Feb 2025 16:45:53 -0500
Subject: [PATCH 16/24] up

---
 docs/src/network_analysis/network_properties.md | 6 ------
 1 file changed, 6 deletions(-)

diff --git a/docs/src/network_analysis/network_properties.md b/docs/src/network_analysis/network_properties.md
index 31586e39da..b3c9cc8f2f 100644
--- a/docs/src/network_analysis/network_properties.md
+++ b/docs/src/network_analysis/network_properties.md
@@ -1,5 +1,3 @@
-# The `NetworkProperties` Object
-
 # [Caching of Network Properties in `ReactionSystems`](@id network_analysis_caching_properties)
 When calling many of the network API functions, Catalyst calculates and caches
 in `rn` a variety of information. For example the first call to
@@ -35,7 +33,3 @@ Catalyst.reset_networkproperties!(rn)
 ```
 Network property functions will then recalculate their associated properties and
 cache the new values the next time they are called.
-
-```@docs; canonical=false
-NetworkProperties
-```

From fea6563e24266afb981aa18caca14d7a85ecdc79 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 5 Feb 2025 18:48:04 -0500
Subject: [PATCH 17/24] up

---
 docs/src/api/core_api.md             | 367 +++++++++++++++++++++++++++
 docs/src/api/network_analysis_api.md |  37 +++
 2 files changed, 404 insertions(+)
 create mode 100644 docs/src/api/core_api.md
 create mode 100644 docs/src/api/network_analysis_api.md

diff --git a/docs/src/api/core_api.md b/docs/src/api/core_api.md
new file mode 100644
index 0000000000..176717244e
--- /dev/null
+++ b/docs/src/api/core_api.md
@@ -0,0 +1,367 @@
+# [Catalyst.jl API](@id api)
+```@meta
+CurrentModule = Catalyst
+```
+
+## Reaction network generation and representation
+Catalyst provides the [`@reaction_network`](@ref) macro for generating a
+complete network, stored as a [`ReactionSystem`](@ref), which in turn is
+composed of [`Reaction`](@ref)s. `ReactionSystem`s can be converted to other
+`ModelingToolkit.AbstractSystem`s, including a `ModelingToolkit.ODESystem`,
+`ModelingToolkit.SDESystem`, or `ModelingToolkit.JumpSystem`.
+
+When using the [`@reaction_network`](@ref) macro, Catalyst will automatically
+attempt to detect what is a species and what is a parameter. Everything that
+appear as a substrate or product in some reaction will be treated as a species,
+while all remaining symbols will be considered parameters (corresponding to
+those symbols that only appear within rate expressions and/or as stoichiometric
+coefficients). I.e. in
+```julia
+rn = @reaction_network begin
+    k*X, Y --> W
+end
+```
+`Y` and `W` will all be classified as chemical species, while `k` and `X` will
+be classified as parameters.
+
+The [`ReactionSystem`](@ref) generated by the [`@reaction_network`](@ref) macro
+is a `ModelingToolkit.AbstractSystem` that symbolically represents a system of
+chemical reactions. In some cases it can be convenient to bypass the macro and
+directly generate a collection of symbolic [`Reaction`](@ref)s and a
+corresponding [`ReactionSystem`](@ref) encapsulating them. Below we illustrate
+with a simple SIR example how a system can be directly constructed, and
+demonstrate how to then generate from the [`ReactionSystem`](@ref) and solve
+corresponding chemical reaction ODE models, chemical Langevin equation SDE
+models, and stochastic chemical kinetics jump process models.
+
+```@example ex1
+using Catalyst, OrdinaryDiffEqTsit5, StochasticDiffEq, JumpProcesses, Plots
+t = default_t()
+@parameters β γ
+@species S(t) I(t) R(t)
+
+rxs = [Reaction(β, [S,I], [I], [1,1], [2])
+       Reaction(γ, [I], [R])]
+@named rs = ReactionSystem(rxs, t)
+rs = complete(rs)
+
+u₀map    = [S => 999.0, I => 1.0, R => 0.0]
+parammap = [β => 1/10000, γ => 0.01]
+tspan    = (0.0, 250.0)
+
+# solve as ODEs
+odesys = convert(ODESystem, rs)
+odesys = complete(odesys)
+oprob = ODEProblem(odesys, u₀map, tspan, parammap)
+sol = solve(oprob, Tsit5())
+p1 = plot(sol, title = "ODE")
+
+# solve as SDEs
+sdesys = convert(SDESystem, rs)
+sdesys = complete(sdesys)
+sprob = SDEProblem(sdesys, u₀map, tspan, parammap)
+sol = solve(sprob, EM(), dt=.01, saveat = 2.0)
+p2 = plot(sol, title = "SDE")
+
+# solve as jump process
+jumpsys = convert(JumpSystem, rs)
+jumpsys = complete(jumpsys)
+u₀map    = [S => 999, I => 1, R => 0]
+dprob = DiscreteProblem(jumpsys, u₀map, tspan, parammap)
+jprob = JumpProblem(jumpsys, dprob, Direct())
+sol = solve(jprob)
+sol = solve(jprob; seed = 1234) # hide
+p3 = plot(sol, title = "jump")
+plot(p1, p2, p3; layout = (3,1))
+Catalyst.PNG(plot(p1, p2, p3; layout = (3,1), fmt = :png, dpi = 200)) # hide
+```
+
+```@docs
+@reaction_network
+@network_component
+make_empty_network
+@reaction
+Reaction
+ReactionSystem
+```
+
+## [Options for the `@reaction_network` DSL](@id api_dsl_options)
+We have [previously described](@ref dsl_advanced_options) how options permits the user to supply non-reaction information to [`ReactionSystem`](@ref) created through the DSL. Here follows a list
+of all options currently available.
+- [`parameters`]:(@ref dsl_advanced_options_declaring_species_and_parameters) Allows the designation of a set of symbols as system parameter.
+- [`species`](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system species.
+- [`variables`](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system non-species variables.
+- [`ivs`](@ref dsl_advanced_options_ivs): Allows the designation of a set of symbols as system independent variables.
+- [`compounds`](@ref chemistry_functionality_compounds): Allows the designation of compound species.
+- [`observables`](@ref dsl_advanced_options_observables): Allows the designation of compound observables.
+- [`default_noise_scaling`](@ref simulation_intro_SDEs_noise_saling): Enables the setting of a default noise scaling expression.
+- [`differentials`](@ref constraint_equations_coupling_constraints): Allows the designation of differentials.
+- [`equations`](@ref constraint_equations_coupling_constraints): Allows the creation of algebraic and/or differential equations.
+- [`continuous_events`](@ref constraint_equations_events): Allows the creation of continuous events.
+- [`discrete_events`](@ref constraint_equations_events): Allows the creation of discrete events.
+- [`combinatoric_ratelaws`](@ref faq_combinatoric_ratelaws): Takes a single option (`true` or `false`), which sets whether to use combinatorial rate laws.
+
+## [ModelingToolkit and Catalyst accessor functions](@id api_accessor_functions)
+A [`ReactionSystem`](@ref) is an instance of a
+`ModelingToolkit.AbstractTimeDependentSystem`, and has a number of fields that
+can be accessed using the Catalyst API and the [ModelingToolkit.jl Abstract
+System
+Interface](https://docs.sciml.ai/ModelingToolkit/stable/basics/AbstractSystem/).
+Below we overview these components.
+
+There are three basic sets of convenience accessors that will return information
+either from a top-level system, the top-level system and all sub-systems that
+are also `ReactionSystem`s (i.e. the full reaction-network), or the top-level
+system, all subs-systems, and all constraint systems (i.e. the full model). To
+retrieve info from just a base [`ReactionSystem`](@ref) `rn`, ignoring
+sub-systems of `rn`, one can use the ModelingToolkit accessors (these provide
+direct access to the corresponding internal fields of the `ReactionSystem`)
+
+* `ModelingToolkit.get_unknowns(rn)` is a vector that collects all the species
+  defined within `rn`, ordered by species and then non-species variables.
+* `Catalyst.get_species(rn)` is a vector of all the species variables in the system. The
+  entries in `get_species(rn)` correspond to the first `length(get_species(rn))`
+  components in `get_unknowns(rn)`.
+* `ModelingToolkit.get_ps(rn)` is a vector that collects all the parameters
+  defined *within* reactions in `rn`.
+* `ModelingToolkit.get_eqs(rn)` is a vector that collects all the
+  [`Reaction`](@ref)s and `Symbolics.Equation` defined within `rn`, ordering all
+  `Reaction`s before `Equation`s.
+* `Catalyst.get_rxs(rn)` is a vector of all the [`Reaction`](@ref)s in `rn`, and
+  corresponds to the first `length(get_rxs(rn))` entries in `get_eqs(rn)`.
+* `ModelingToolkit.get_iv(rn)` is the independent variable used in the system
+  (usually `t` to represent time).
+* `ModelingToolkit.get_systems(rn)` is a vector of all sub-systems of `rn`.
+* `ModelingToolkit.get_defaults(rn)` is a dictionary of all the default values
+  for parameters and species in `rn`.
+
+The preceding accessors do not allocate, directly accessing internal fields of
+the `ReactionSystem`.
+
+To retrieve information from the full reaction network represented by a system
+`rn`, which corresponds to information within both `rn` and all sub-systems, one
+can call:
+
+* `ModelingToolkit.unknowns(rn)` returns all species *and variables* across the
+  system, *all sub-systems*, and all constraint systems. Species are ordered
+  before non-species variables in `unknowns(rn)`, with the first `numspecies(rn)`
+  entries in `unknowns(rn)` being the same as `species(rn)`.
+* [`species(rn)`](@ref) is a vector collecting all the chemical species within
+  the system and any sub-systems that are also `ReactionSystems`.
+* `ModelingToolkit.parameters(rn)` returns all parameters across the
+  system, *all sub-systems*, and all constraint systems.
+* `ModelingToolkit.equations(rn)` returns all [`Reaction`](@ref)s and all
+  `Symbolics.Equations` defined across the system, *all sub-systems*, and all
+  constraint systems. `Reaction`s are ordered ahead of `Equation`s with the
+  first `numreactions(rn)` entries in `equations(rn)` being the same as
+  `reactions(rn)`.
+* [`reactions(rn)`](@ref) is a vector of all the `Reaction`s within the system
+  and any sub-systems that are also `ReactionSystem`s.
+
+These accessors will generally allocate new arrays to store their output unless
+there are no subsystems. In the latter case the usually return the same vector
+as the corresponding `get_*` function.
+
+Below we list the remainder of the Catalyst API accessor functions mentioned
+above.
+
+## Basic system properties
+See [Programmatic Construction of Symbolic Reaction Systems](@ref
+programmatic_CRN_construction) for examples and [ModelingToolkit and Catalyst
+Accessor Functions](@ref api_accessor_functions) for more details on the basic
+accessor functions.
+
+```@docs
+species
+Catalyst.get_species
+nonspecies
+reactions
+Catalyst.get_rxs
+nonreactions
+numspecies
+numparams
+numreactions
+speciesmap
+paramsmap
+isautonomous
+```
+
+## Coupled reaction/equation system properties
+The following system property accessor functions are primarily relevant to reaction system [coupled
+to differential and/or algebraic equations](@ref constraint_equations).
+```@docs
+ModelingToolkit.has_alg_equations
+ModelingToolkit.alg_equations
+ModelingToolkit.has_diff_equations
+ModelingToolkit.diff_equations
+```
+
+## Basic species properties
+The following functions permits the querying of species properties.
+```@docs
+isspecies
+Catalyst.isconstant
+Catalyst.isbc
+Catalyst.isvalidreactant
+```
+
+## Basic reaction properties
+```@docs
+ismassaction
+dependents
+dependants
+substoichmat
+prodstoichmat
+netstoichmat
+reactionrates
+```
+
+## [Reaction metadata](@id api_rx_metadata)
+The following functions permits the retrieval of [reaction metadata](@ref dsl_advanced_options_reaction_metadata).
+```@docs
+Catalyst.hasnoisescaling
+Catalyst.getnoisescaling
+Catalyst.hasdescription
+Catalyst.getdescription
+Catalyst.hasmisc
+Catalyst.getmisc
+```
+
+## [Functions to extend or modify a network](@id api_network_extension_and_modification)
+`ReactionSystem`s can be programmatically extended using [`ModelingToolkit.extend`](@ref) and
+[`ModelingToolkit.compose`](@ref).
+
+```@docs
+setdefaults!
+ModelingToolkit.extend
+ModelingToolkit.compose
+Catalyst.flatten
+```
+
+## Network comparison
+```@docs
+==(rn1::Reaction, rn2::Reaction)
+isequivalent
+==(rn1::ReactionSystem, rn2::ReactionSystem)
+```
+
+## [Network visualization](@id network_visualization)
+[Latexify](https://korsbo.github.io/Latexify.jl/stable/) can be used to convert
+networks to LaTeX equations by
+```julia
+using Latexify
+latexify(rn)
+```
+An optional argument, `form` allows using `latexify` to display a reaction
+network's ODE (as generated by the reaction rate equation) or SDE (as generated
+by the chemical Langevin equation) form:
+```julia
+latexify(rn; form=:ode)
+```
+```julia
+latexify(rn; form=:sde)
+```
+(As of writing this, an upstream bug causes the SDE form to be erroneously
+displayed as the ODE form)
+
+Finally, another optional argument (`expand_functions=true`) automatically expands functions defined by Catalyst (such as `mm`). To disable this, set `expand_functions=false`.
+
+Reaction networks can be plotted using the `GraphMakie` extension, which is loaded whenever all of `Catalyst`, `GraphMakie`, and `NetworkLayout` are loaded (note that a Makie backend, like `CairoMakie`, must be loaded as well). The two functions for plotting networks are `plot_network` and `plot_complexes`, which are two distinct representations. 
+```@docs
+plot_network(::ReactionSystem)
+plot_complexes(::ReactionSystem)
+```
+
+## [Rate laws](@id api_rate_laws)
+As the underlying [`ReactionSystem`](@ref) is comprised of `ModelingToolkit`
+expressions, one can directly access the generated rate laws, and using
+`ModelingToolkit` tooling generate functions or Julia `Expr`s from them.
+```@docs
+oderatelaw
+jumpratelaw
+mm
+mmr
+hill
+hillr
+hillar
+```
+
+## Transformations
+```@docs
+Base.convert
+JumpInputs
+ModelingToolkit.structural_simplify
+set_default_noise_scaling
+```
+
+## Chemistry-related functionalities
+Various functionalities primarily relevant to modelling of chemical systems (but potentially also in biology).
+```@docs
+@compound
+@compounds
+iscompound
+components
+coefficients
+component_coefficients
+```
+
+## Unit validation
+```@docs
+validate(rx::Reaction; info::String = "")
+validate(rs::ReactionSystem, info::String="")
+```
+
+## Utility functions
+```@docs
+symmap_to_varmap
+```
+
+## [Spatial modelling](@id api_lattice_simulations)
+The first step of spatial modelling is to create a so-called `LatticeReactionSystem`:
+```@docs
+LatticeReactionSystem
+```
+
+The following functions can be used to querying the properties of `LatticeReactionSystem`s:
+```@docs
+reactionsystem
+Catalyst.spatial_reactions
+Catalyst.lattice
+Catalyst.num_verts
+Catalyst.num_edges
+Catalyst.num_species
+Catalyst.spatial_species
+Catalyst.vertex_parameters
+Catalyst.edge_parameters
+Catalyst.edge_iterator
+Catalyst.is_transport_system
+has_cartesian_lattice
+has_masked_lattice
+has_grid_lattice
+has_graph_lattice
+grid_size
+grid_dims
+```
+In addition, most accessor functions for normal `ReactionSystem`s (such as `species` and `parameters`) works when applied to `LatticeReactionSystem`s as well.
+
+The following two helper functions can be used to create non-uniform parameter values.
+```@docs
+make_edge_p_values
+make_directed_edge_values
+```
+
+The following functions can be used to access, or change, species or parameter values stored in problems, integrators, and solutions that are based on `LatticeReactionSystem`s.
+```@docs
+lat_getu
+lat_setu!
+lat_getp
+lat_setp!
+rebuild_lat_internals!
+```
+
+Finally, we provide the following helper functions to plot and animate spatial lattice simulations.
+```@docs
+lattice_plot
+lattice_animation
+lattice_kymograph
+```
diff --git a/docs/src/api/network_analysis_api.md b/docs/src/api/network_analysis_api.md
new file mode 100644
index 0000000000..dad88e5ca6
--- /dev/null
+++ b/docs/src/api/network_analysis_api.md
@@ -0,0 +1,37 @@
+## [Network analysis and representations](@id api_network_analysis)
+Note, currently API functions for network analysis and conservation law analysis
+do not work with constant species (currently only generated by SBMLToolkit).
+
+For more information about these functions, please see the sections of the docs on
+[network ODE representation](@ref network_analysis_odes) and [chemical reaction network theory](@ref network_analysis_structural_aspects).
+
+```@docs
+conservationlaws
+conservedquantities
+conservedequations
+conservationlaw_constants
+ReactionComplexElement
+ReactionComplex
+reactioncomplexmap
+reactioncomplexes
+incidencemat
+incidencematgraph
+complexstoichmat
+complexoutgoingmat
+fluxmat
+adjacencymat
+laplacianmat
+massactionvector
+linkageclasses
+deficiency
+linkagedeficiencies
+satisfiesdeficiencyone
+satisfiesdeficiencyzero
+subnetworks
+isreversible
+isweaklyreversible
+iscomplexbalanced
+isdetailedbalanced
+robustspecies
+reset_networkproperties!
+```

From efd85622d987585fe07b85f2e022aa927d90f3f1 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 5 Feb 2025 18:48:16 -0500
Subject: [PATCH 18/24] up

---
 docs/src/api.md | 405 ------------------------------------------------
 1 file changed, 405 deletions(-)
 delete mode 100644 docs/src/api.md

diff --git a/docs/src/api.md b/docs/src/api.md
deleted file mode 100644
index a47538c044..0000000000
--- a/docs/src/api.md
+++ /dev/null
@@ -1,405 +0,0 @@
-# [Catalyst.jl API](@id api)
-```@meta
-CurrentModule = Catalyst
-```
-
-## Reaction network generation and representation
-Catalyst provides the [`@reaction_network`](@ref) macro for generating a
-complete network, stored as a [`ReactionSystem`](@ref), which in turn is
-composed of [`Reaction`](@ref)s. `ReactionSystem`s can be converted to other
-`ModelingToolkit.AbstractSystem`s, including a `ModelingToolkit.ODESystem`,
-`ModelingToolkit.SDESystem`, or `ModelingToolkit.JumpSystem`.
-
-When using the [`@reaction_network`](@ref) macro, Catalyst will automatically
-attempt to detect what is a species and what is a parameter. Everything that
-appear as a substrate or product in some reaction will be treated as a species,
-while all remaining symbols will be considered parameters (corresponding to
-those symbols that only appear within rate expressions and/or as stoichiometric
-coefficients). I.e. in
-```julia
-rn = @reaction_network begin
-    k*X, Y --> W
-end
-```
-`Y` and `W` will all be classified as chemical species, while `k` and `X` will
-be classified as parameters.
-
-The [`ReactionSystem`](@ref) generated by the [`@reaction_network`](@ref) macro
-is a `ModelingToolkit.AbstractSystem` that symbolically represents a system of
-chemical reactions. In some cases it can be convenient to bypass the macro and
-directly generate a collection of symbolic [`Reaction`](@ref)s and a
-corresponding [`ReactionSystem`](@ref) encapsulating them. Below we illustrate
-with a simple SIR example how a system can be directly constructed, and
-demonstrate how to then generate from the [`ReactionSystem`](@ref) and solve
-corresponding chemical reaction ODE models, chemical Langevin equation SDE
-models, and stochastic chemical kinetics jump process models.
-
-```@example ex1
-using Catalyst, OrdinaryDiffEqTsit5, StochasticDiffEq, JumpProcesses, Plots
-t = default_t()
-@parameters β γ
-@species S(t) I(t) R(t)
-
-rxs = [Reaction(β, [S,I], [I], [1,1], [2])
-       Reaction(γ, [I], [R])]
-@named rs = ReactionSystem(rxs, t)
-rs = complete(rs)
-
-u₀map    = [S => 999.0, I => 1.0, R => 0.0]
-parammap = [β => 1/10000, γ => 0.01]
-tspan    = (0.0, 250.0)
-
-# solve as ODEs
-odesys = convert(ODESystem, rs)
-odesys = complete(odesys)
-oprob = ODEProblem(odesys, u₀map, tspan, parammap)
-sol = solve(oprob, Tsit5())
-p1 = plot(sol, title = "ODE")
-
-# solve as SDEs
-sdesys = convert(SDESystem, rs)
-sdesys = complete(sdesys)
-sprob = SDEProblem(sdesys, u₀map, tspan, parammap)
-sol = solve(sprob, EM(), dt=.01, saveat = 2.0)
-p2 = plot(sol, title = "SDE")
-
-# solve as jump process
-jumpsys = convert(JumpSystem, rs)
-jumpsys = complete(jumpsys)
-u₀map    = [S => 999, I => 1, R => 0]
-dprob = DiscreteProblem(jumpsys, u₀map, tspan, parammap)
-jprob = JumpProblem(jumpsys, dprob, Direct())
-sol = solve(jprob)
-sol = solve(jprob; seed = 1234) # hide
-p3 = plot(sol, title = "jump")
-plot(p1, p2, p3; layout = (3,1))
-Catalyst.PNG(plot(p1, p2, p3; layout = (3,1), fmt = :png, dpi = 200)) # hide
-```
-
-```@docs
-@reaction_network
-@network_component
-make_empty_network
-@reaction
-Reaction
-ReactionSystem
-```
-
-## [Options for the `@reaction_network` DSL](@id api_dsl_options)
-We have [previously described](@ref dsl_advanced_options) how options permits the user to supply non-reaction information to [`ReactionSystem`](@ref) created through the DSL. Here follows a list
-of all options currently available.
-- [`parameters`]:(@ref dsl_advanced_options_declaring_species_and_parameters) Allows the designation of a set of symbols as system parameter.
-- [`species`](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system species.
-- [`variables`](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system non-species variables.
-- [`ivs`](@ref dsl_advanced_options_ivs): Allows the designation of a set of symbols as system independent variables.
-- [`compounds`](@ref chemistry_functionality_compounds): Allows the designation of compound species.
-- [`observables`](@ref dsl_advanced_options_observables): Allows the designation of compound observables.
-- [`default_noise_scaling`](@ref simulation_intro_SDEs_noise_saling): Enables the setting of a default noise scaling expression.
-- [`differentials`](@ref constraint_equations_coupling_constraints): Allows the designation of differentials.
-- [`equations`](@ref constraint_equations_coupling_constraints): Allows the creation of algebraic and/or differential equations.
-- [`continuous_events`](@ref constraint_equations_events): Allows the creation of continuous events.
-- [`discrete_events`](@ref constraint_equations_events): Allows the creation of discrete events.
-- [`combinatoric_ratelaws`](@ref faq_combinatoric_ratelaws): Takes a single option (`true` or `false`), which sets whether to use combinatorial rate laws.
-
-## [ModelingToolkit and Catalyst accessor functions](@id api_accessor_functions)
-A [`ReactionSystem`](@ref) is an instance of a
-`ModelingToolkit.AbstractTimeDependentSystem`, and has a number of fields that
-can be accessed using the Catalyst API and the [ModelingToolkit.jl Abstract
-System
-Interface](https://docs.sciml.ai/ModelingToolkit/stable/basics/AbstractSystem/).
-Below we overview these components.
-
-There are three basic sets of convenience accessors that will return information
-either from a top-level system, the top-level system and all sub-systems that
-are also `ReactionSystem`s (i.e. the full reaction-network), or the top-level
-system, all subs-systems, and all constraint systems (i.e. the full model). To
-retrieve info from just a base [`ReactionSystem`](@ref) `rn`, ignoring
-sub-systems of `rn`, one can use the ModelingToolkit accessors (these provide
-direct access to the corresponding internal fields of the `ReactionSystem`)
-
-* `ModelingToolkit.get_unknowns(rn)` is a vector that collects all the species
-  defined within `rn`, ordered by species and then non-species variables.
-* `Catalyst.get_species(rn)` is a vector of all the species variables in the system. The
-  entries in `get_species(rn)` correspond to the first `length(get_species(rn))`
-  components in `get_unknowns(rn)`.
-* `ModelingToolkit.get_ps(rn)` is a vector that collects all the parameters
-  defined *within* reactions in `rn`.
-* `ModelingToolkit.get_eqs(rn)` is a vector that collects all the
-  [`Reaction`](@ref)s and `Symbolics.Equation` defined within `rn`, ordering all
-  `Reaction`s before `Equation`s.
-* `Catalyst.get_rxs(rn)` is a vector of all the [`Reaction`](@ref)s in `rn`, and
-  corresponds to the first `length(get_rxs(rn))` entries in `get_eqs(rn)`.
-* `ModelingToolkit.get_iv(rn)` is the independent variable used in the system
-  (usually `t` to represent time).
-* `ModelingToolkit.get_systems(rn)` is a vector of all sub-systems of `rn`.
-* `ModelingToolkit.get_defaults(rn)` is a dictionary of all the default values
-  for parameters and species in `rn`.
-
-The preceding accessors do not allocate, directly accessing internal fields of
-the `ReactionSystem`.
-
-To retrieve information from the full reaction network represented by a system
-`rn`, which corresponds to information within both `rn` and all sub-systems, one
-can call:
-
-* `ModelingToolkit.unknowns(rn)` returns all species *and variables* across the
-  system, *all sub-systems*, and all constraint systems. Species are ordered
-  before non-species variables in `unknowns(rn)`, with the first `numspecies(rn)`
-  entries in `unknowns(rn)` being the same as `species(rn)`.
-* [`species(rn)`](@ref) is a vector collecting all the chemical species within
-  the system and any sub-systems that are also `ReactionSystems`.
-* `ModelingToolkit.parameters(rn)` returns all parameters across the
-  system, *all sub-systems*, and all constraint systems.
-* `ModelingToolkit.equations(rn)` returns all [`Reaction`](@ref)s and all
-  `Symbolics.Equations` defined across the system, *all sub-systems*, and all
-  constraint systems. `Reaction`s are ordered ahead of `Equation`s with the
-  first `numreactions(rn)` entries in `equations(rn)` being the same as
-  `reactions(rn)`.
-* [`reactions(rn)`](@ref) is a vector of all the `Reaction`s within the system
-  and any sub-systems that are also `ReactionSystem`s.
-
-These accessors will generally allocate new arrays to store their output unless
-there are no subsystems. In the latter case the usually return the same vector
-as the corresponding `get_*` function.
-
-Below we list the remainder of the Catalyst API accessor functions mentioned
-above.
-
-## Basic system properties
-See [Programmatic Construction of Symbolic Reaction Systems](@ref
-programmatic_CRN_construction) for examples and [ModelingToolkit and Catalyst
-Accessor Functions](@ref api_accessor_functions) for more details on the basic
-accessor functions.
-
-```@docs
-species
-Catalyst.get_species
-nonspecies
-reactions
-Catalyst.get_rxs
-nonreactions
-numspecies
-numparams
-numreactions
-speciesmap
-paramsmap
-isautonomous
-```
-
-## Coupled reaction/equation system properties
-The following system property accessor functions are primarily relevant to reaction system [coupled
-to differential and/or algebraic equations](@ref constraint_equations).
-```@docs
-ModelingToolkit.has_alg_equations
-ModelingToolkit.alg_equations
-ModelingToolkit.has_diff_equations
-ModelingToolkit.diff_equations
-```
-
-## Basic species properties
-The following functions permits the querying of species properties.
-```@docs
-isspecies
-Catalyst.isconstant
-Catalyst.isbc
-Catalyst.isvalidreactant
-```
-
-## Basic reaction properties
-```@docs
-ismassaction
-dependents
-dependants
-substoichmat
-prodstoichmat
-netstoichmat
-reactionrates
-```
-
-## [Reaction metadata](@id api_rx_metadata)
-The following functions permits the retrieval of [reaction metadata](@ref dsl_advanced_options_reaction_metadata).
-```@docs
-Catalyst.hasnoisescaling
-Catalyst.getnoisescaling
-Catalyst.hasdescription
-Catalyst.getdescription
-Catalyst.hasmisc
-Catalyst.getmisc
-```
-
-## [Functions to extend or modify a network](@id api_network_extension_and_modification)
-`ReactionSystem`s can be programmatically extended using [`ModelingToolkit.extend`](@ref) and
-[`ModelingToolkit.compose`](@ref).
-
-```@docs
-setdefaults!
-ModelingToolkit.extend
-ModelingToolkit.compose
-Catalyst.flatten
-```
-
-## Network analysis and representations
-Note, currently API functions for network analysis and conservation law analysis
-do not work with constant species (currently only generated by SBMLToolkit).
-
-For more information about these functions, please see the sections of the docs on
-[network ODE representation](@ref network_analysis_odes) and [chemical reaction network theory](@ref network_analysis_structural_aspects).
-
-```@docs
-conservationlaws
-conservedquantities
-conservedequations
-conservationlaw_constants
-ReactionComplexElement
-ReactionComplex
-reactioncomplexmap
-reactioncomplexes
-incidencemat
-incidencematgraph
-complexstoichmat
-complexoutgoingmat
-fluxmat
-adjacencymat
-laplacianmat
-massactionvector
-linkageclasses
-deficiency
-linkagedeficiencies
-satisfiesdeficiencyone
-satisfiesdeficiencyzero
-subnetworks
-isreversible
-isweaklyreversible
-iscomplexbalanced
-isdetailedbalanced
-robustspecies
-reset_networkproperties!
-```
-
-## Network comparison
-```@docs
-==(rn1::Reaction, rn2::Reaction)
-isequivalent
-==(rn1::ReactionSystem, rn2::ReactionSystem)
-```
-
-## [Network visualization](@id network_visualization)
-[Latexify](https://korsbo.github.io/Latexify.jl/stable/) can be used to convert
-networks to LaTeX equations by
-```julia
-using Latexify
-latexify(rn)
-```
-An optional argument, `form` allows using `latexify` to display a reaction
-network's ODE (as generated by the reaction rate equation) or SDE (as generated
-by the chemical Langevin equation) form:
-```julia
-latexify(rn; form=:ode)
-```
-```julia
-latexify(rn; form=:sde)
-```
-(As of writing this, an upstream bug causes the SDE form to be erroneously
-displayed as the ODE form)
-
-Finally, another optional argument (`expand_functions=true`) automatically expands functions defined by Catalyst (such as `mm`). To disable this, set `expand_functions=false`.
-
-Reaction networks can be plotted using the `GraphMakie` extension, which is loaded whenever all of `Catalyst`, `GraphMakie`, and `NetworkLayout` are loaded (note that a Makie backend, like `CairoMakie`, must be loaded as well). The two functions for plotting networks are `plot_network` and `plot_complexes`, which are two distinct representations. 
-```@docs
-plot_network(::ReactionSystem)
-plot_complexes(::ReactionSystem)
-```
-
-## [Rate laws](@id api_rate_laws)
-As the underlying [`ReactionSystem`](@ref) is comprised of `ModelingToolkit`
-expressions, one can directly access the generated rate laws, and using
-`ModelingToolkit` tooling generate functions or Julia `Expr`s from them.
-```@docs
-oderatelaw
-jumpratelaw
-mm
-mmr
-hill
-hillr
-hillar
-```
-
-## Transformations
-```@docs
-Base.convert
-JumpInputs
-ModelingToolkit.structural_simplify
-set_default_noise_scaling
-```
-
-## Chemistry-related functionalities
-Various functionalities primarily relevant to modelling of chemical systems (but potentially also in biology).
-```@docs
-@compound
-@compounds
-iscompound
-components
-coefficients
-component_coefficients
-```
-
-## Unit validation
-```@docs
-validate(rx::Reaction; info::String = "")
-validate(rs::ReactionSystem, info::String="")
-```
-
-## Utility functions
-```@docs
-symmap_to_varmap
-```
-
-## [Spatial modelling](@id api_lattice_simulations)
-The first step of spatial modelling is to create a so-called `LatticeReactionSystem`:
-```@docs
-LatticeReactionSystem
-```
-
-The following functions can be used to querying the properties of `LatticeReactionSystem`s:
-```@docs
-reactionsystem
-Catalyst.spatial_reactions
-Catalyst.lattice
-Catalyst.num_verts
-Catalyst.num_edges
-Catalyst.num_species
-Catalyst.spatial_species
-Catalyst.vertex_parameters
-Catalyst.edge_parameters
-Catalyst.edge_iterator
-Catalyst.is_transport_system
-has_cartesian_lattice
-has_masked_lattice
-has_grid_lattice
-has_graph_lattice
-grid_size
-grid_dims
-```
-In addition, most accessor functions for normal `ReactionSystem`s (such as `species` and `parameters`) works when applied to `LatticeReactionSystem`s as well.
-
-The following two helper functions can be used to create non-uniform parameter values.
-```@docs
-make_edge_p_values
-make_directed_edge_values
-```
-
-The following functions can be used to access, or change, species or parameter values stored in problems, integrators, and solutions that are based on `LatticeReactionSystem`s.
-```@docs
-lat_getu
-lat_setu!
-lat_getp
-lat_setp!
-rebuild_lat_internals!
-```
-
-Finally, we provide the following helper functions to plot and animate spatial lattice simulations.
-```@docs
-lattice_plot
-lattice_animation
-lattice_kymograph
-```

From 1d8c4cdcdd9a379baae68d36307795f20f7c24ad Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 5 Feb 2025 18:49:06 -0500
Subject: [PATCH 19/24] change pages

---
 docs/pages.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/docs/pages.jl b/docs/pages.jl
index 4fc7c8053b..61008601bd 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -67,6 +67,9 @@ pages = Any[
         "spatial_modelling/spatial_jump_simulations.md"
     ],
     "FAQs" => "faqs.md",
-    "API" => "api.md",
+    "API" => Any[
+                 "core_api.md",
+                 "network_analysis_api.md"
+                ],
     "Developer Documentation" => "devdocs/dev_guide.md"
 ]

From 25f138a4a8a6105cb495da2e37cd9a944bcea939 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 5 Feb 2025 18:49:50 -0500
Subject: [PATCH 20/24] up

---
 docs/pages.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/pages.jl b/docs/pages.jl
index 61008601bd..7c5e2d2b64 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -68,8 +68,8 @@ pages = Any[
     ],
     "FAQs" => "faqs.md",
     "API" => Any[
-                 "core_api.md",
-                 "network_analysis_api.md"
+                 "api/core_api.md",
+                 "api/network_analysis_api.md"
                 ],
     "Developer Documentation" => "devdocs/dev_guide.md"
 ]

From d4cb6a71eca7bb4d4227e469df67ed59d6cfba48 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Wed, 5 Feb 2025 18:51:46 -0500
Subject: [PATCH 21/24] up

---
 docs/src/api/network_analysis_api.md | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/docs/src/api/network_analysis_api.md b/docs/src/api/network_analysis_api.md
index dad88e5ca6..f435478c11 100644
--- a/docs/src/api/network_analysis_api.md
+++ b/docs/src/api/network_analysis_api.md
@@ -1,4 +1,8 @@
-## [Network analysis and representations](@id api_network_analysis)
+# [Network analysis and representations](@id api_network_analysis)
+```@meta
+CurrentModule = Catalyst
+```
+
 Note, currently API functions for network analysis and conservation law analysis
 do not work with constant species (currently only generated by SBMLToolkit).
 

From 9f61b34e15d4961e771e34f9245e4c1f861d54f9 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Sun, 9 Feb 2025 19:15:04 -0500
Subject: [PATCH 22/24] fix docstrings and docs

---
 .../{crnt.md => crn_theory.md}                | 61 ++++++++++---------
 docs/src/network_analysis/odes.md             | 10 ++-
 src/network_analysis.jl                       | 47 +++++++-------
 3 files changed, 64 insertions(+), 54 deletions(-)
 rename docs/src/network_analysis/{crnt.md => crn_theory.md} (85%)

diff --git a/docs/src/network_analysis/crnt.md b/docs/src/network_analysis/crn_theory.md
similarity index 85%
rename from docs/src/network_analysis/crnt.md
rename to docs/src/network_analysis/crn_theory.md
index 4b3715fe55..48e93ef898 100644
--- a/docs/src/network_analysis/crnt.md
+++ b/docs/src/network_analysis/crn_theory.md
@@ -1,24 +1,27 @@
+> Note, currently API functions for network analysis and conservation law analysis
+> do not work with constant species (currently only generated by SBMLToolkit).
+
 # [Chemical Reaction Network Theory](@id network_analysis_structural_aspects)
 
 The systems of ODEs or stochastic chemical kinetics models that arise from chemical
-reaction networks can often have their steady-state properties known in advance,
+reaction networks can often have their steady-state properties (number of steady states, etc.) known in advance,
 simply by analyzing the graph structure of the network. The subfield of chemistry
 and math studying this relationship is called [Chemical Reaction Network Theory](https://en.wikipedia.org/wiki/Chemical_reaction_network_theory).
 
-In this tutorial we introduce several of the Catalyst API functions for network
+In this tutorial we give a broad overview of chemical reaction network theory, building on the discussion of the structure of 
+the reaction network ODEs in the previous section. We also introduce several of the Catalyst API functions for network
 analysis. A complete summary of the exported functions is given in the API
 section
-[`Network-Analysis-and-Representations`](https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Network-Analysis-and-Representations).
+[Network Analysis and Representations](@ref api_network_analysis).
 
 Broadly, results from chemical reaction network theory relate a purely
 graph-structural property (e.g. deficiency) to dynamical properties of the reaction system
-(e.g. complex balance). We'll now illustrate some of the types of network properties that 
-Catalyst can determine, using the [reaction complex representation](@ref network_analysis_reaction_complexes) in these calculations.
-
-Note, currently API functions for network analysis and conservation law analysis
-do not work with constant species (currently only generated by SBMLToolkit).
+(e.g. complex balance). The relevant graph here is the one corresponding to the [reaction complex representation](@ref network_analysis_reaction_complexes)
+of the network, where the nodes represent the reaction complexes and the edges represent reactions.
+Let us illustrate some of the types of network properties that 
+Catalyst can determine.
 
-Throughout this seciton, we will consider the [reaction complex representation](@ref network_analysis_reaction_complexes) of the following reaction network.
+To begin, consider the following reaction network:
 ```@example s1
 using Catalyst
 rn = @reaction_network begin
@@ -35,7 +38,7 @@ with graph
 plot_complexes(rn)
 ```
 
-### [Linkage classes and sub-networks of the reaction network](@id network_analysis_structural_aspects_linkage)
+## [Linkage classes and sub-networks of the reaction network](@id network_analysis_structural_aspects_linkage)
 The preceding reaction complex graph shows that `rn` is composed of two
 disconnected sub-graphs, one containing the complexes ``A+B``, ``C``, ``D+E``, and
 ``F``, the other containing the complexes ``2A``, ``B + G``, and ``H``. These sets,
@@ -70,10 +73,7 @@ and,
  plot_complexes(subnets[2])
 ```
 
-![subnetwork_2](../assets/complex_subnets2.svg)
-
-
-### [Deficiency of the network](@id network_analysis_structural_aspects_deficiency)
+## [Deficiency of the network](@id network_analysis_structural_aspects_deficiency)
 A famous theorem in Chemical Reaction Network Theory, the Deficiency Zero
 Theorem [^1], allows us to use knowledge of the net stoichiometry matrix and the
 linkage classes of a *mass action* [RRE ODE system](@ref network_analysis_matrix_vector_representation) to draw conclusions about the
@@ -128,7 +128,7 @@ Quoting Feinberg [^1]
 > linkage classes will allow.
 
 
-### [Reversibility of the network](@id network_analysis_structural_aspects_reversibility)
+## [Reversibility of the network](@id network_analysis_structural_aspects_reversibility)
 A reaction network is *reversible* if the "arrows" of the reactions are
 symmetric so that every reaction is accompanied by its reverse reaction.
 Catalyst's API provides the [`isreversible`](@ref) function to determine whether
@@ -149,7 +149,7 @@ isreversible(rn)
 ```
 Consider another example,
 ```@example s1
-rn = @reaction_network begin
+rn2 = @reaction_network begin
   (k1,k2),A <--> B
   k3, A + C --> D
   k4, D --> B+E
@@ -175,7 +175,7 @@ isweaklyreversible(rn, subnets)
 Every reversible network is also weakly reversible, but not every weakly
 reversible network is reversible.
 
-### [Deficiency Zero Theorem](@id network_analysis_structural_aspects_deficiency_zero_theorem)
+## [Deficiency Zero Theorem](@id network_analysis_structural_aspects_deficiency_zero_theorem)
 Knowing the deficiency and weak reversibility of a mass action chemical reaction
 network ODE model allows us to make inferences about the corresponding
 steady state behavior. Before illustrating how this works for one example, we
@@ -249,16 +249,16 @@ asymptotically stable.
 
 As a final example, consider the following network from [^1]
 ```@example s1
-rn = @reaction_network begin
+def0_rn = @reaction_network begin
   (k1,k2),A <--> 2B
   (k3,k4), A + C <--> D
   k5, B+E --> C + D
 end
-reactioncomplexes(rn)
-subnets = subnetworks(rn)
-isma = all(rx -> ismassaction(rx,rn), reactions(rn))
-def = deficiency(rn)
-iswr = isweaklyreversible(rn, subnets)
+reactioncomplexes(def0_rn)
+subnets = subnetworks(def0_rn)
+isma = all(rx -> ismassaction(rx,def0_rn), reactions(def0_rn))
+def = deficiency(def0_rn)
+iswr = isweaklyreversible(def0_rn, subnets)
 isma,def,iswr
 ```
 which we see is mass action and has deficiency zero, but is not weakly
@@ -269,26 +269,29 @@ RRE ODEs cannot have a positive equilibrium solution.
 satisfiesdeficiencyzero
 ```
 
-### Deficiency One Theorem
+## Deficiency One Theorem
 Very analogous to the deficiency zero theorem is the deficiency one theorem. The deficiency one theorem applies to a network with the following properties: 
 1. The deficiency of each *linkage class* of the network is at most 1, 
 2. The sum of the linkage class deficiencies is the total deficiency of the network, and
-3. Each linkage class has at most one terminal linkage class. 
+3. Each linkage class has at most one terminal linkage class, which is a linkage class that is 1) strongly connected, and 2) has no outgoing reactions.
+
+For the set of reactions $A \to B, B \to A, B \to C$, $\{A, B, C\}$ is a linkage class, and $\{A, B\}$ is a strong linkage class (since A is reachable from B and vice versa). However, $\{A, B\}$ is not a terminal linkage class, because the reaction $B \to C$ goes to a complex outside the linkage class. 
 
 If these conditions are met, then the network will have at most one steady state in each stoichiometric compatibility class for any choice of rate constants and parameters.
 Unlike the deficiency zero theorem, networks obeying the deficiency one theorem are not guaranteed to have stable solutions.
+
 ```@docs; canonical=false
 satisfiesdeficiencyone
 ```
 
-### [Complex and Detailed Balance](@id network_analysis_complex_and_detailed_balance)
+## [Complex and Detailed Balance](@id network_analysis_complex_and_detailed_balance)
 A reaction network's steady state is **complex-balanced** if the total production of each *complex* is zero at the steady state. A reaction network's steady state is **detailed balanced** if every reaction is balanced by its reverse reaction at the steady-state (this corresponds to the usual notion of chemical equilibrium; note that this requires every reaction be reversible). 
 
 Note that detailed balance at a given steady state implies complex balance for that steady state, i.e. detailed balance is a stronger property than complex balance.
 
 Remarkably, having just one positive steady state that is complex (detailed) balance implies that complex (detailed) balance obtains at *every* positive steady state, so we say that a network is complex (detailed) balanced if any one of its steady states are complex (detailed) balanced. Additionally, there will be exactly one steady state in every positive stoichiometric compatibility class, and this steady state is asymptotically stable. (For proofs of these results, please consult Martin Feinberg's *Foundations of Chemical Reaction Network Theory*[^1]). So knowing that a network is complex balanced is really quite powerful.
 
-Let's check that the reaction network defined above is complex balanced by providing a set of rates:
+Let's check whether the reaction network defined above is complex balanced by providing a set of rates:
 ```@example s1
 rates = Dict([:k1 => 2.4, :k2 => 4., :k3 => 10., :k4 => 5.5, :k5 => 0.4])
 iscomplexbalanced(rn, rates)
@@ -297,7 +300,7 @@ iscomplexbalanced(rn, rates)
 We can do a similar check for detailed balance. Let us make the reaction network 
 ```@example s1
 rn1 = @reaction_network begin
-  (k1,k2),A <--> 2B
+  (k1,k2), A <--> 2B
   (k3,k4), A + C <--> D
   (k5,k6), B + E --> C + D
 end
@@ -311,7 +314,7 @@ iscomplexbalanced
 isdetailedbalanced
 ```
 
-### [Concentration Robustness](@id network_analysis_concentration_robustness)
+## [Concentration Robustness](@id network_analysis_concentration_robustness)
 Certain reaction networks have species that do not change their concentration,
 regardless of the system is perturbed to a different stoichiometric compatibility
 class. This is a very useful property to have in biological contexts, where it might
diff --git a/docs/src/network_analysis/odes.md b/docs/src/network_analysis/odes.md
index a93f9316ef..af3ebec43f 100644
--- a/docs/src/network_analysis/odes.md
+++ b/docs/src/network_analysis/odes.md
@@ -1,10 +1,10 @@
 # [Decomposing the Reaction Network ODEs](@id network_analysis_odes)
 
 In this tutorial we will discuss the specific mathematical 
-structure of the ODEs that arise from the mass-action dynamics
+structure of the [ODEs that arise from the mass-action dynamics](@ref math_models_in_catalyst_rre_odes)
 of chemical reaction networks, and decompose them as a product
 of matrices that describe the network.  A complete summary of 
-the exported functions is given in the API section [`Network-Analysis-and-Representations`](https://docs.sciml.ai/Catalyst/stable/api/catalyst_api/#Network-Analysis-and-Representations). See Feinberg's *Foundations of Chemical Reaction
+the exported functions is given in the API section [Network Analysis and Representations](@ref api_network_analysis). Please consult Feinberg's *Foundations of Chemical Reaction
 Network Theory*[^1] for more discussion about the concepts on this page.
 
 Note, currently API functions for network analysis and conservation law analysis
@@ -240,14 +240,16 @@ K = fluxmat(rn)
 
 Since we have that $\mathbf{v} = K\Phi$, we can rewrite the above decompositions as follows: 
 ```math
+\begin{align}
 \frac{d\mathbf{x}}{dt} &= N \mathbf{v}(\mathbf{x}(t),t) \\
 &= N K \Phi(\mathbf{x}(t),t) \\
 &= Z B K \Phi(\mathbf{x}(t),t). 
+\end{align}
 ```
 
 The final matrix to discuss is the product of $A_k = BK$, which is a $C$-by-$C$ matrix that turns out to be exactly the negative of the [graph Laplacian](https://en.wikipedia.org/wiki/Laplacian_matrix) of the weighted graph whose nodes are reaction complexes and whose edges represent reactions, weighted by the rate constants. The API function for $A_k$ is the `laplacianmat`: 
 ```@example s1
-A_k = incidencemat(rn)
+A_k = laplacianmat(rn)
 ```
 We can check that 
 ```@example s1
@@ -256,10 +258,12 @@ A_k == B * K
 
 In sum, we have that
 ```math
+\begin{align}
 \frac{d\mathbf{x}}{dt} &= N \mathbf{v}(\mathbf{x}(t),t) \\
 &= N K \Phi(\mathbf{x}(t),t) \\
 &= Z B K \Phi(\mathbf{x}(t),t). 
 &= Z A_k \Phi(\mathbf{x}(t),t). 
+\end{align}
 ```
 
 All three of the objects introduced in this section (the flux matrix, mass-action vector, Laplacian matrix) will return symbolic outputs by default, but can be made to return numerical outputs if values are specified. 
diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index d5aa9e0bb8..978a3d7652 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -195,8 +195,9 @@ end
 @doc raw"""
     laplacianmat(rn::ReactionSystem, pmap=Dict(); sparse=false)
 
-    Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the graph Laplacian, and ``Φ`` is the [`massactionvector`](@ref). ``A_k`` is an n-by-n matrix, where n is the number of complexes, where ``A_{ij} = k_{ij}`` if a reaction exists between the two complexes and 0 otherwise. 
-    Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. 
+Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the graph Laplacian, and ``Φ`` is the [`massactionvector`](@ref). ``A_k`` is an n-by-n matrix, where n is the number of complexes, where ``A_{ij} = k_{ij}`` if a reaction exists between the two complexes and 0 otherwise. 
+
+Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. 
 
     **Warning**: Unlike other Catalyst functions, the `laplacianmat` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs, and to ensure that multiplying the sparse form of the matrix will work.
 """
@@ -209,10 +210,10 @@ end
 @doc raw"""
     fluxmat(rn::ReactionSystem, pmap = Dict(); sparse=false)
 
-    Return an r×c matrix ``K`` such that, if complex ``j`` is the substrate complex of reaction ``i``, then ``K_{ij} = k``, the rate constant for this reaction. Mostly a helper function for the network Laplacian, [`laplacianmat`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix and ``Φ(x)`` is the [`massactionvector`](@ref).
+Return an r×c matrix ``K`` such that, if complex ``j`` is the substrate complex of reaction ``i``, then ``K_{ij} = k``, the rate constant for this reaction. Mostly a helper function for the network Laplacian, [`laplacianmat`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix and ``Φ(x)`` is the [`massactionvector`](@ref).
     Returns a symbolic matrix by default, but will return a numerical matrix if rate constants are specified as a `Tuple`, `Vector`, or `Dict` of symbol-value pairs via `pmap`.
 
-    **Warning**: Unlike other Catalyst functions, the `fluxmat` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs, and to ensure that multiplying the sparse form of the matrix will work.
+**Warning**: Unlike other Catalyst functions, the `fluxmat` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs, and to ensure that multiplying the sparse form of the matrix will work.
 """
 function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false)
     deps = Set()
@@ -286,11 +287,13 @@ end
 """
     massactionvector(rn::ReactionSystem, scmap = Dict(); combinatoric_ratelaws = true)
 
-    Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be ``A*B``, and for the complex 2X + Y, the corresponding entry would be ``X^2*Y``. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the [`laplacianmat`](@ref). This utility returns ``Φ(x)``.
-    Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap.
-    If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws). Will default to the default for the system.
+Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be ``A*B``, and for the complex 2X + Y, the corresponding entry would be ``X^2*Y``. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the [`laplacianmat`](@ref). This utility returns ``Φ(x)``.
+
+Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap.
 
-    **Warning**: Unlike other Catalyst functions, the `massactionvector` function will return a `Vector{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs.
+If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws). Will default to the default for the system.
+
+**Warning**: Unlike other Catalyst functions, the `massactionvector` function will return a `Vector{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs.
 """
 function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn))
     r = numreactions(rn)
@@ -1050,17 +1053,16 @@ end
 """
     adjacencymat(rs::ReactionSystem, pmap = Dict(); sparse = false)
 
-    Given a reaction system with n complexes, outputs an n-by-n matrix where R_{ij} is the rate 
-    constant of the reaction between complex i and complex j. Accepts a dictionary, vector, or tuple 
-    of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
+Given a reaction system with n complexes, outputs an n-by-n matrix where R_{ij} is the rate 
+constant of the reaction between complex i and complex j. Accepts a dictionary, vector, or tuple 
+of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
 
-    Equivalent to the adjacency matrix of the weighted graph whose weights are the 
-    reaction rates. 
-    Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. 
+Equivalent to the adjacency matrix of the weighted graph whose weights are the reaction rates. 
+Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. 
 
-    The difference between this matrix and [`fluxmat`](@ref) is quite subtle. The non-zero entries to both matrices are the rate constants. However:
-        - In `fluxmat`, the rows represent complexes and columns represent reactions, and an entry (c, r) is non-zero if reaction r's source complex is c. 
-        - In `adjacencymat`, the rows and columns both represent complexes, and an entry (c1, c2) is non-zero if there is a reaction c1 --> c2.
+The difference between this matrix and [`fluxmat`](@ref) is quite subtle. The non-zero entries to both matrices are the rate constants. However:
+    - In `fluxmat`, the rows represent complexes and columns represent reactions, and an entry (c, r) is non-zero if reaction r's source complex is c. 
+    - In `adjacencymat`, the rows and columns both represent complexes, and an entry (c1, c2) is non-zero if there is a reaction c1 --> c2.
 """
 function adjacencymat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false)
     deps = Set()
@@ -1176,8 +1178,9 @@ end
 """
     cycles(rs::ReactionSystem)
 
-    Returns the matrix of a basis of cycles (or flux vectors), or a basis for reaction fluxes for which the system is at steady state. 
-    These correspond to right eigenvectors of the stoichiometric matrix. Equivalent to [`fluxmodebasis`](@ref). 
+Returns the matrix of a basis of cycles (or flux vectors), or a basis for reaction fluxes for which the system is at steady state. 
+
+These correspond to right eigenvectors of the stoichiometric matrix. Equivalent to [`fluxmodebasis`](@ref). 
 """
 function cycles(rs::ReactionSystem)
     nps = get_networkproperties(rs)
@@ -1211,7 +1214,7 @@ end
 """
     fluxvectors(rs::ReactionSystem)
 
-    See documentation for [`cycles`](@ref). 
+See documentation for [`cycles`](@ref). 
 """
 function fluxvectors(rs::ReactionSystem)
     cycles(rs)
@@ -1256,9 +1259,9 @@ end
 """
     robustspecies(rn::ReactionSystem)
 
-    Return a vector of indices corresponding to species that are concentration robust, i.e. for every positive equilbrium, the concentration of species s will be the same. 
+Return a vector of indices corresponding to species that are concentration robust, i.e. for every positive equilbrium, the concentration of species s will be the same. 
 
-    Note: This function currently only works for networks of deficiency one, and is not currently guaranteed to return *all* the concentration-robust species in the network. Any species returned by the function will be robust, but this may not include all of them. Use with caution. Support for higher deficiency networks and necessary conditions for robustness will be coming in the future.  
+Note: This function currently only works for networks of deficiency one, and is not currently guaranteed to return *all* the concentration-robust species in the network. Any species returned by the function will be robust, but this may not include all of them. Use with caution. More extensive robustness analysis can be performed using the CatalystNetworkAnalysis package.
 """
 function robustspecies(rn::ReactionSystem)
     complexes, D = reactioncomplexes(rn)

From d8140b606ccaab37bcc7df9461a362ce41417fa8 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Mon, 10 Feb 2025 09:26:38 -0500
Subject: [PATCH 23/24] fix file name

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

diff --git a/docs/pages.jl b/docs/pages.jl
index 7c5e2d2b64..235b7380ec 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -26,7 +26,7 @@ pages = Any[
     ],
     "Network Analysis" => Any[
         "network_analysis/odes.md",
-        "network_analysis/crnt.md",
+        "network_analysis/crn_theory.md",
         "network_analysis/network_properties.md",
         "network_analysis/advanced_analysis.md"
     ],

From af5737c75de0dc73d47b83a14d101e2dcf1cf2e8 Mon Sep 17 00:00:00 2001
From: vyudu <vincent.duyuan@gmail.com>
Date: Mon, 10 Feb 2025 15:26:26 -0500
Subject: [PATCH 24/24] fix blocks

---
 docs/src/network_analysis/crn_theory.md | 11 +++--------
 docs/src/network_analysis/odes.md       |  5 ++++-
 2 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/docs/src/network_analysis/crn_theory.md b/docs/src/network_analysis/crn_theory.md
index 48e93ef898..6f607d8fa8 100644
--- a/docs/src/network_analysis/crn_theory.md
+++ b/docs/src/network_analysis/crn_theory.md
@@ -149,7 +149,7 @@ isreversible(rn)
 ```
 Consider another example,
 ```@example s1
-rn2 = @reaction_network begin
+rn = @reaction_network begin
   (k1,k2),A <--> B
   k3, A + C --> D
   k4, D --> B+E
@@ -291,19 +291,14 @@ Note that detailed balance at a given steady state implies complex balance for t
 
 Remarkably, having just one positive steady state that is complex (detailed) balance implies that complex (detailed) balance obtains at *every* positive steady state, so we say that a network is complex (detailed) balanced if any one of its steady states are complex (detailed) balanced. Additionally, there will be exactly one steady state in every positive stoichiometric compatibility class, and this steady state is asymptotically stable. (For proofs of these results, please consult Martin Feinberg's *Foundations of Chemical Reaction Network Theory*[^1]). So knowing that a network is complex balanced is really quite powerful.
 
-Let's check whether the reaction network defined above is complex balanced by providing a set of rates:
+Let's check whether the deficiency 0 reaction network that we defined above is complex balanced by providing a set of rates:
 ```@example s1
 rates = Dict([:k1 => 2.4, :k2 => 4., :k3 => 10., :k4 => 5.5, :k5 => 0.4])
 iscomplexbalanced(rn, rates)
 ```
 
-We can do a similar check for detailed balance. Let us make the reaction network 
+We can do a similar check for detailed balance.
 ```@example s1
-rn1 = @reaction_network begin
-  (k1,k2), A <--> 2B
-  (k3,k4), A + C <--> D
-  (k5,k6), B + E --> C + D
-end
 isdetailedbalanced(rn, rates)
 ```
 
diff --git a/docs/src/network_analysis/odes.md b/docs/src/network_analysis/odes.md
index af3ebec43f..f0a9a1e9c4 100644
--- a/docs/src/network_analysis/odes.md
+++ b/docs/src/network_analysis/odes.md
@@ -253,9 +253,12 @@ A_k = laplacianmat(rn)
 ```
 We can check that 
 ```@example s1
-A_k == B * K
+isequal(A_k, B * K)
 ```
 
+Note that we have used `isequal` instead of `==` here because `laplacianmat`
+returns a `Matrix{Num}`, since some of its entries are symbolic rate constants.
+
 In sum, we have that
 ```math
 \begin{align}