fix docstrings and docs

Author: vyudu

docs/src/network_analysis/crnt.md renamed to docs/src/network_analysis/crn_theory.md

@@ -1,23 +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.
+(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
@@ -34,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,
@@ -69,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
@@ -127,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
@@ -148,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
@@ -174,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
@@ -248,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
@@ -268,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)
@@ -296,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
@@ -310,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

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. 

src/network_analysis.jl

@@ -198,6 +198,8 @@ end
     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.
 
+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.
 """
 function laplacianmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false)
@@ -209,10 +211,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 +288,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.
+
+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.
+**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)