up

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",

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)

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.

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