Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: beginning to refactor network analysis docs #1142

Draft
wants to merge 2 commits into
base: master
Choose a base branch
from
+447 −245
Draft

WIP: beginning to refactor network analysis docs #1142

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
626b737
init
vyudu Dec 2, 2024
bb070e9
up
vyudu Dec 8, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions docs/pages.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
295 changes: 50 additions & 245 deletions docs/src/model_creation/network_analysis.md → docs/src/network_analysis/crnt.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -234,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,
Expand Down Expand Up @@ -273,7 +63,8 @@ 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
Expand Down Expand Up @@ -327,7 +118,8 @@ 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
Expand Down Expand Up @@ -376,7 +168,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
Expand Down Expand Up @@ -466,42 +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.

## [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)
```@docs
satisfiesdeficiencyzero
```
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)
# 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.

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
```
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 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).

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:
```@example s1
rates = Dict([:k1 => 2.4, :k2 => 4., :k3 => 10., :k4 => 5.5, :k5 => 0.4])
iscomplexbalanced(rn, rates)
```
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)

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
(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)
37 changes: 37 additions & 0 deletions docs/src/network_analysis/network_properties.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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.
Loading
Loading