Skip to content

Commit

Permalink
init
Browse files Browse the repository at this point in the history
  • Loading branch information
@vyudu
vyudu committed Dec 2, 2024
1 parent 19d9652 commit 626b737
Show file tree
Hide file tree
Showing 2 changed files with 360 additions and 215 deletions.
229 changes: 14 additions & 215 deletions docs/src/model_creation/network_analysis.md → docs/src/network_analysis/crnt.md
View file
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 Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading

0 comments on commit 626b737

Please sign in to comment.