Network Analysis in Catalyst

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.

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

We first load Catalyst and construct our model of the repressilator

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 α K n δ γ β μ

\[ \begin{align*} \require{mhchem} \ce{ \varnothing &->[$\frac{\alpha K^{n}}{K^{n} + P_3^{n}}$] m_1}\\ \ce{ \varnothing &->[$\frac{\alpha K^{n}}{K^{n} + P_1^{n}}$] m_2}\\ \ce{ \varnothing &->[$\frac{\alpha K^{n}}{K^{n} + P_2^{n}}$] m_3}\\ \ce{ m_1 &<=>[$\delta$][$\gamma$] \varnothing}\\ \ce{ m_2 &<=>[$\delta$][$\gamma$] \varnothing}\\ \ce{ m_3 &<=>[$\delta$][$\gamma$] \varnothing}\\ \ce{ m_1 &->[$\beta$] m_1 + P_1}\\ \ce{ m_2 &->[$\beta$] m_2 + P_2}\\ \ce{ m_3 &->[$\beta$] m_3 + P_3}\\ \ce{ P_1 &->[$\mu$] \varnothing}\\ \ce{ P_2 &->[$\mu$] \varnothing}\\ \ce{ P_3 &->[$\mu$] \varnothing} \end{align*} \]

In the Using 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.

g = Graph(repressilator)

We also showed in the Using Catalyst tutorial that the reaction rate equation ODE model for the repressilator is

\[\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

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

\[\mathbf{x}(t) = \begin{pmatrix} x_1(t) \\ \vdots \\ x_M(t)) \end{pmatrix}.\]

For the repressilator, $\mathbf{x}(t)$ is just

x = species(repressilator)

6-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
 m₁(t)
 m₂(t)
 m₃(t)
 P₁(t)
 P₂(t)
 P₃(t)

The RRE ODEs satisfied by $\mathbf{x}(t)$ are then

\[\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

\[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

N = netstoichmat(repressilator)

6×15 Matrix{Int64}:
 1  0  0  -1  1   0  0   0  0  0  0  0   0   0   0
 0  1  0   0  0  -1  1   0  0  0  0  0   0   0   0
 0  0  1   0  0   0  0  -1  1  0  0  0   0   0   0
 0  0  0   0  0   0  0   0  0  1  0  0  -1   0   0
 0  0  0   0  0   0  0   0  0  0  1  0   0  -1   0
 0  0  0   0  0   0  0   0  0  0  0  1   0   0  -1

and by using the oderatelaw function

rxs = reactions(repressilator)
ν = oderatelaw.(rxs)

15-element Vector{SymbolicUtils.Symbolic{Real}}:
 Catalyst.hillr(P₃(t), α, K, n)
 Catalyst.hillr(P₁(t), α, K, n)
 Catalyst.hillr(P₂(t), α, K, n)
 δ*m₁(t)
 γ
 δ*m₂(t)
 γ
 δ*m₃(t)
 γ
 β*m₁(t)
 β*m₂(t)
 β*m₃(t)
 μ*P₁(t)
 μ*P₂(t)
 μ*P₃(t)

Note, as oderatelaw takes just one reaction as input we use broadcasting to apply it to each element of rxs.

Let's check this really gives the same ODEs as Catalyst. Here is what Catalyst generates by converting to an ODESystem

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)]

6-element Vector{SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}}:
 γ + Catalyst.hillr(P₃(t), α, K, n) - δ*m₁(t)
 γ + Catalyst.hillr(P₁(t), α, K, n) - δ*m₂(t)
 γ + Catalyst.hillr(P₂(t), α, K, n) - δ*m₃(t)
 β*m₁(t) - μ*P₁(t)
 β*m₂(t) - μ*P₂(t)
 β*m₃(t) - μ*P₃(t)

whereas our matrix-vector representation gives

odes2 = N * ν

6-element Vector{Any}:
 γ + Catalyst.hillr(P₃(t), α, K, n) - δ*m₁(t)
 γ + Catalyst.hillr(P₁(t), α, K, n) - δ*m₂(t)
 γ + Catalyst.hillr(P₂(t), α, K, n) - δ*m₃(t)
 β*m₁(t) - μ*P₁(t)
 β*m₂(t) - μ*P₂(t)
 β*m₃(t) - μ*P₃(t)

Let's check these are equal symbolically

isequal(odes, odes2)

true

Reaction Complex Representation

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

rn = @reaction_network begin
 k*A, 2*A + 3*B --> A + 2*C + D
 b, C + D --> 2*A + 3*B
end k b

\[ \begin{align*} \require{mhchem} \ce{ 2 A + 3 B &->[$A k$] A + 2 C + D}\\ \ce{ C + D &->[$b$] 2 A + 3 B} \end{align*} \]

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.

\[\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

species(rn)

4-element Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}:
 A(t)
 B(t)
 C(t)
 D(t)

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

\[\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,

Z = complexstoichmat(rn)

4×3 Matrix{Int64}:
 2  1  0
 3  0  0
 0  2  1
 0  1  1

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,

B = incidencemat(rn)

3×2 Matrix{Int64}:
 -1   1
  1   0
  0  -1

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

\[\frac{dx}{dt} = Z B \mathbf{v}(\mathbf{x}(t),t).\]

Let's verify that $N = Z B$,

N = netstoichmat(rn)
N == Z*B

true

Reaction complexes give an alternative way to visualize a reaction network graph. Catalyst's complexgraph command will calculate the complexes of a network and then show how they are related. For example,

complexgraph(rn)

gives

while for the repressilator we find

complexgraph(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 conversion of substrate complexes into product complexes where the rate is an expression involving chemical species.

Aspects of Reaction Network Structure

The reaction complex representation can be exploited via 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.

rn = @reaction_network begin
     (k1,k2), A + B <--> C
     k3, C --> D+E
     (k4,k5), D+E <--> F
     (k6,k7), 2A <--> B+G
     k8, B+G --> H
     k9, H --> 2A
end k1 k2 k3 k4 k5 k6 k7 k8 k9

\[ \begin{align*} \require{mhchem} \ce{ A + B &<=>[$k1$][$k2$] C}\\ \ce{ C &->[$k3$] D + E}\\ \ce{ D + E &<=>[$k4$][$k5$] F}\\ \ce{ 2 A &<=>[$k6$][$k7$] B + G}\\ \ce{ B + G &->[$k8$] H}\\ \ce{ H &->[$k9$] 2 A} \end{align*} \]

with graph

complexgraph(rn)

Linkage classes and sub-networks of the reaction network

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, $\{A+B, C, D+E, F\}$ and $\{2A, B + G,H\}$ are called the "linkage classes" of the reaction network. The function linkageclasses will calculate these for a given network, returning a vector of the integer indices of reaction complexes participating in each set of linkage-classes. Note, indices of reaction complexes can be determined from the ordering returned by reactioncomplexes.

# we must first calculate the reaction complexes -- they are cached in rn
reactioncomplexes(rn)

# now we can calculate the linkage classes
lcs = linkageclasses(rn)

2-element Vector{Vector{Int64}}:
 [1, 2, 3, 4]
 [5, 6, 7]

It can often be convenient to obtain the disconnected sub-networks as distinct ReactionSystems, which are returned by the subnetworks function:

subnets = subnetworks(rn)

# check the reactions in each subnetwork
reactions.(subnets)

2-element Vector{Vector{Reaction}}:
 [k1, A + B --> C, k2, C --> A + B, k3, C --> D + E, k4, D + E --> F, k5, F --> D + E]
 [k6, 2*A --> B + G, k7, B + G --> 2*A, k8, B + G --> H, k9, H --> 2*A]

The graphs of the reaction complexes in the two sub-networks are then

  complexgraph(subnets[1])

and,

 complexgraph(subnets[2])

Deficiency of the network

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 system's possible steady-states. In this section we'll see how Catalyst can calculate a network's deficiency.

The rank, $r$, of a reaction network is defined as the dimension of the subspace spanned by the net stoichiometry vectors of the reaction-network ^[1], i.e. the span of the columns of the net stoichiometry matrix N. It corresponds to the number of independent species in a chemical reaction network. That is, if we calculate the linear conservation laws of a network, and use them to eliminate the dependent species of the network, we will have $r$ independent species remaining. For our current example the conservation laws are given by

# first we calculate the conservation laws -- they are cached in rn
conservationlaws(rn)

# then we display them as equations for the dependent variables
conservedequations(rn)

3-element Vector{Equation}:
 E(t) ~ D(t) + _ConLaw[1]
 G(t) ~ B(t) + C(t) + D(t) + F(t) + _ConLaw[2]
 H(t) ~ _ConLaw[3] - 0.5A(t) - B(t) - 1.5C(t) - 1.5D(t) - 1.5F(t)

Here the parameters _ConLaw[i] represent the constants of the three conservation laws, and we see that there are three dependent species that could be eliminated. As

numspecies(rn)

8

we find that there are five independent species. Let's check this is correct:

using LinearAlgebra
rank(netstoichmat(rn)) == 5

true

So we know that the rank of our reaction network is five.

The deficiency, $\delta$, of a reaction network is defined as

\[\delta = \textrm{(number of complexes)} - \textrm{(number of linkage classes)} - \textrm{(rank)}.\]

For our network this is $7 - 2 - 5 = 0$, which we can calculate in Catalyst as

# first we calculate the reaction complexes of rn and cache them in rn
reactioncomplexes(rn)

# then we can calculate the deficiency
δ = deficiency(rn)

0

Quoting Feinberg ^[1]

Deficiency zero networks are ones for which the reaction vectors [i.e. net stoichiometry vectors] are as independent as the partition of complexes into linkage classes will allow.

Reversibility of the network

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 function to determine whether a reaction network is reversible. As an example, consider

rn = @reaction_network begin
  (k1,k2),A <--> B
  (k3,k4),A + C <--> D
  (k5,k6),D <--> B+E
  (k7,k8),B+E <--> A+C
end k1 k2 k3 k4 k5 k6 k7 k8

# calculate the set of reaction complexes
reactioncomplexes(rn)

# test if the system is reversible
isreversible(rn)

true

Consider another example,

rn = @reaction_network begin
  (k1,k2),A <--> B
  k3, A + C --> D
  k4, D --> B+E
  k5, B+E --> A+C
end k1 k2 k3 k4 k5
reactioncomplexes(rn)
isreversible(rn)

false

complexgraph(rn)

It is evident from the preceding graph that the network is not reversible. However, it satisfies a weaker property in that there is a path from each reaction complex back to itself within its associated subgraph. This is known as weak reversiblity. One can test a network for weak reversibility by using the isweaklyreversible function:

# need subnetworks from the reaction network first
subnets = subnetworks(rn)
isweaklyreversible(rn, subnets)

true

Every reversible network is also weakly reversible, but not every weakly reversible network is reversible.

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 need one last definition.

Recall that in the matrix-vector representation for the RRE ODEs, the entries, $N_{m k}$, of the stoichiometry matrix, $N$, give the net change in species $m$ due to reaction $k$. If we let $\mathbf{N}_k$ denote the $k$th column of this matrix, this vector corresponds to the change in the species state vector, $\mathbf{x}(t)$, due to reaction $k$, i.e. when reaction $k$ occurs $\mathbf{x}(t) \to \mathbf{x}(t) + \mathbf{N}_k$. Moreover, by integrating the ODE

\[\frac{d\mathbf{x}}{dt} = N \mathbf{v}(\mathbf{x}(t)) = \sum_{k=1}^{K} v_k(\mathbf{x}(t)) \, \mathbf{N}_k\]

we find

\[\mathbf{x}(t) = \mathbf{x}(0) + \sum_{k=1}^K \left(\int_0^t v_k(\mathbf{x})(s) \, ds\right) \mathbf{N}_k,\]

which demonstrates that $\mathbf{x}(t) - \mathbf{x}(0)$ is always given by a linear combination of the stochiometry vectors, i.e.

\[\DeclareMathOperator{\span}{span} \mathbf{x}(t) - \mathbf{x}(0) \in \span\{\mathbf{N}_k \}.\]

In particular, this says that $\mathbf{x}(t)$ lives in the translation of the $\span\{\mathbf{N}_k \}$ by $\mathbf{x}(0)$ which we write as $(\mathbf{x}(0) + \span\{\mathbf{N}_k\})$. In fact, since the solution should stay non-negative, if we let $\bar{\mathbb{R}}_+^{M}$ denote the subset of vectors in $\mathbb{R}^{M}$ with non-negative components, the possible physical values for the solution, $\mathbf{x}(t)$, must be in the set

\[(\mathbf{x}(0) + \span\{\mathbf{N}_k\}) \cap \bar{\mathbb{R}}_+^{M}.\]

This set is called the stoichiometric compatibility class of $\mathbf{x}(t)$. The key property of stoichiometric compatibility classes is that they are invariant under the RRE ODE's dynamics, i.e. a solution will always remain within the subspace given by the stoichiometric compatibility class. Finally, we note that the positive stoichiometric compatibility class generated by $\mathbf{x}(0)$ is just $(\mathbf{x}(0) + \span\{\mathbf{N}_k\}) \cap \mathbb{R}_+^{M}$, where $\mathbb{R}_+^{M}$ denotes the vectors in $\mathbb{R}^M$ with strictly positive components.

With these definitions we can now see how knowing the deficiency and weak reversibility of the network can tell us about its steady-state behavior. Consider the previous example, which we know is weakly reversible. Its deficiency is

deficiency(rn)

0

We also verify that the system is purely mass action (though it is apparent from the network's definition):

all(rx -> ismassaction(rx, rn), reactions(rn))

true

We can therefore apply the Deficiency Zero Theorem to draw conclusions about the system's steady-state behavior. The Deficiency Zero Theorem (roughly) says that a mass action network with deficiency zero satisfies

1. If the network is weakly reversible, then independent of the reaction rate constants the RRE ODEs have exactly one equilibrium solution within each positive stoichiometric compatibility class. That equilibrium is locally asymptotically stable.
2. If the network is not weakly reversible, then the RRE ODEs cannot admit a positive equilibrium solution.

See ^[1] for a more precise statement, proof, and additional examples.

We can therefore conclude that for any initial condition that is positive, and hence in some positive stoichiometric compatibility class, rn will have exactly one equilibrium solution which will be positive and locally asymptotically stable.

As a final example, consider the following network from ^[1]

rn = @reaction_network begin
  (k1,k2),A <--> 2B
  (k3,k4), A + C <--> D
  k5, B+E --> C + D
end k1 k2 k3 k4 k5
reactioncomplexes(rn)
subnets = subnetworks(rn)
isma = all(rx -> ismassaction(rx,rn), reactions(rn))
def = deficiency(rn)
iswr = isweaklyreversible(rn, subnets)
isma,def,iswr

(true, 0, false)

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

When calling many of the network API functions, Catalyst calculates and caches in rn a variety of information. For example the first call to

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

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

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, for example after modifying a network (see addspecies!, addparam!, or addreaction!). This can be done by calling

Catalyst.reset_networkproperties!(rn)

Network property functions will then recalculate their associated properties and cache the new values the next time they are called.

References