Decomposing the Reaction Network ODEs
In this tutorial we will discuss the specific mathematical structure of the ODEs that arise from the mass-action dynamics 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. 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 do not work with constant species (which are generated by SBML, and can be declared in Catalyst as well.
Network representation of the Repressilator ReactionSystem
We first load Catalyst and construct our model of the repressilator
using Catalyst, CairoMakie, GraphMakie, NetworkLayout
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
\[ \begin{align*} \varnothing &\xrightarrow{\frac{K^{n} \alpha}{K^{n} + \mathtt{P_3}^{n}}} \mathrm{\mathtt{m_1}} \\ \varnothing &\xrightarrow{\frac{K^{n} \alpha}{\mathtt{P_1}^{n} + K^{n}}} \mathrm{\mathtt{m_2}} \\ \varnothing &\xrightarrow{\frac{K^{n} \alpha}{K^{n} + \mathtt{P_2}^{n}}} \mathrm{\mathtt{m_3}} \\ \mathrm{\mathtt{m_1}} &\xrightleftharpoons[\gamma]{\delta} \varnothing \\ \mathrm{\mathtt{m_2}} &\xrightleftharpoons[\gamma]{\delta} \varnothing \\ \mathrm{\mathtt{m_3}} &\xrightleftharpoons[\gamma]{\delta} \varnothing \\ \mathrm{\mathtt{m_1}} &\xrightarrow{\beta} \mathrm{\mathtt{m_1}} + \mathrm{\mathtt{P_1}} \\ \mathrm{\mathtt{m_2}} &\xrightarrow{\beta} \mathrm{\mathtt{m_2}} + \mathrm{\mathtt{P_2}} \\ \mathrm{\mathtt{m_3}} &\xrightarrow{\beta} \mathrm{\mathtt{m_3}} + \mathrm{\mathtt{P_3}} \\ \mathrm{\mathtt{P_1}} &\xrightarrow{\mu} \varnothing \\ \mathrm{\mathtt{P_2}} &\xrightarrow{\mu} \varnothing \\ \mathrm{\mathtt{P_3}} &\xrightarrow{\mu} \varnothing \end{align*} \]
In the Model Visualization 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 = plot_network(repressilator)

We also showed in the Introduction to Catalyst tutorial that the reaction rate equation (RRE) 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{SymbolicUtils.BasicSymbolic{Real}}:
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.BasicSymbolic{Real}}:
((K^n)*α) / (P₃(t)^n + K^n)
((K^n)*α) / (K^n + P₁(t)^n)
((K^n)*α) / (K^n + P₂(t)^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 that 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.BasicSymbolic{Real}}:
((K^n)*α) / (P₃(t)^n + K^n) + γ - m₁(t)*δ
((K^n)*α) / (K^n + P₁(t)^n) + γ - m₂(t)*δ
((K^n)*α) / (K^n + P₂(t)^n) + γ - m₃(t)*δ
-P₁(t)*μ + m₁(t)*β
m₂(t)*β - P₂(t)*μ
-P₃(t)*μ + m₃(t)*β
whereas our matrix-vector representation gives
odes2 = N * ν
6-element Vector{Any}:
((K^n)*α) / (P₃(t)^n + K^n) + γ - m₁(t)*δ
((K^n)*α) / (K^n + P₁(t)^n) + γ - m₂(t)*δ
((K^n)*α) / (K^n + P₂(t)^n) + γ - m₃(t)*δ
-P₁(t)*μ + m₁(t)*β
m₂(t)*β - P₂(t)*μ
-P₃(t)*μ + m₃(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, 2A + 3B --> A + 2C + D
b, C + D --> 2A + 3B
end
\[ \begin{align*} 2 \mathrm{A} + 3 \mathrm{B} &\xrightarrow{k} \mathrm{A} + 2 \mathrm{C} + \mathrm{D} \\ \mathrm{C} + \mathrm{D} &\xrightarrow{b} 2 \mathrm{A} + 3 \mathrm{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 $k$. 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{SymbolicUtils.BasicSymbolic{Real}}:
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 plot_complexes
command will calculate the complexes of a network and then show how they are related. For example, we can run
plot_complexes(rn)

while for the repressilator we find
plot_complexes(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 the conversion of substrate complexes into product complexes where the rate is an expression involving chemical species.
Full decomposition of the reaction network ODEs (flux matrix and mass-action vector)
So far we have covered two equivalent descriptions of the chemical reaction network ODEs:
\[\begin{align} \frac{d\mathbf{x}}{dt} &= N \mathbf{v}(\mathbf{x}(t),t) \\ &= Z B \mathbf{v}(\mathbf{x}(t),t). \end{align}\]
In this section we discuss a further decomposition of the ODEs. Recall the reaction rate vector $\mathbf{v}$, which is a vector of length $R$ whose elements are the rate expressions for each reaction. Its elements can be written as
\[\mathbf{v}_{y \rightarrow y'} = k_{y \rightarrow y'} \mathbf{x}^y,\]
where $\mathbf{x}^y = \prod_s x_s^{y_s}$ denotes the mass-action product of the substrate complex $y$ from the $y \rightarrow y'$ reaction. We can define a new vector called the mass action vector $\Phi(\mathbf{x}(t))$, a vector of length $C$ whose elements are the mass action products of each complex:
Φ = massactionvector(rn)
\[ \begin{equation} \left[ \begin{array}{c} \frac{1}{12} \left( A\left( t \right) \right)^{2} \left( B\left( t \right) \right)^{3} \\ \frac{1}{2} \left( C\left( t \right) \right)^{2} D\left( t \right) A\left( t \right) \\ D\left( t \right) C\left( t \right) \\ \end{array} \right] \end{equation} \]
An important thing to note is this function assumes combinatoric ratelaws, meaning that mass-action products will get rescaled by factorial factors. For instance, note that the mass-action product for the complex 2A + 3B
has a factor of 1/12, corresponding to 1/(2! 3!). This option can be turned off with combinatoric_ratelaws = false
.
Φ_2 = massactionvector(rn; combinatoric_ratelaws = false)
\[ \begin{equation} \left[ \begin{array}{c} \left( A\left( t \right) \right)^{2} \left( B\left( t \right) \right)^{3} \\ \left( C\left( t \right) \right)^{2} D\left( t \right) A\left( t \right) \\ D\left( t \right) C\left( t \right) \\ \end{array} \right] \end{equation} \]
Then the reaction rate vector $\mathbf{v}$ can be written as
\[\mathbf{v}(\mathbf{x}(t)) = K \Phi(\mathbf{x}(t))\]
where $K$ is an $R$-by-$C$ matrix called the flux matrix, where $K_{rc}$ is the rate constant of reaction $r$ if $c$ is the index of the substrate complex of reaction $r$, and 0 otherwise. In Catalyst, the API function for $K$ is fluxmat
:
K = fluxmat(rn)
\[ \begin{equation} \left[ \begin{array}{ccc} k & 0 & 0 \\ 0 & 0 & b \\ \end{array} \right] \end{equation} \]
Since we have that $\mathbf{v} = K\Phi$, we can rewrite the above decompositions as follows:
\[\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 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
:
A_k = laplacianmat(rn)
\[ \begin{equation} \left[ \begin{array}{ccc} - k & 0 & b \\ k & 0 & 0 \\ 0 & 0 & - b \\ \end{array} \right] \end{equation} \]
We can check that
isequal(A_k, B * K)
true
Note that we have used isequal
instead of ==
here because laplacianmat
returns a Matrix{Num}
, since some of its entries are symbolic rate constants (symbolic variables and Num
s cannot be compared using ==
, since a == b
is interpreted as a symbolic expression).
In summary, we have that
\[\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, and Laplacian matrix) will return symbolic outputs by default, but can be made to return numerical outputs if values are specified. For example, massactionvector
will return a numerical output if a set of species concentrations is supplied using a dictionary, tuple, or vector of Symbol-value pairs.
concmap = Dict([:A => 3., :B => 5., :C => 2.4, :D => 1.5])
massactionvector(rn, concmap)
3-element Vector{Float64}:
93.75
12.96
3.5999999999999996
fluxmat
and laplacianmat
will return numeric matrices if a set of rate constants and other parameters are supplied the same way.
parammap = Dict([:k => 12., :b => 8.])
fluxmat(rn, parammap)
2×3 Matrix{Float64}:
12.0 0.0 0.0
0.0 0.0 8.0
laplacianmat(rn, parammap)
3×3 Matrix{Float64}:
-12.0 0.0 8.0
12.0 0.0 0.0
0.0 0.0 -8.0
Symbolic ODE functions
In some cases it might be useful to generate the function defining the system of ODEs as a symbolic Julia function that can be used for further analysis. This can be done using Symbolics' build_function
, which takes a symbolic expression and a set of desired arguments, and converts it into a Julia function taking those arguments.
Let's build the full symbolic function corresponding to our ODE system. build_function
will return two expressions, one for a function that outputs a new vector for the result, and one for a function that modifies the input in-place. Either expression can then be evaluated to return a Julia function.
parammap = Dict([:k => 12., :b => 8.])
K = fluxmat(rn, parammap)
odes = N * K * Φ
f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn))
f = eval(f_oop_expr)
c = [3., 5., 2., 6.]
f(c)
4-element Vector{Float64}:
-933.0
-3087.0
2154.0
1029.0
The generated f
now corresponds to the $f(\mathbf{x}(t))$ on the right-hand side of $\frac{d\mathbf{x}(t)}{dt} = f(\mathbf{x}(t))$. Given a vector of species concentrations $c$, f
will return the rate of change of each species. Steady state concentration vectors c_ss
will satisfy f(c_ss) = zeros(length(species(rn)))
.
Above we have generated a numeric rate matrix to substitute the rate constants into the symbolic expressions. We could have used a symbolic rate matrix, but then we would need to define the parameters k, b
, so that the function f
knows what k
and b
in its output refer to.
@parameters k b
K = fluxmat(rn)
odes = N * K * Φ
f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn))
f = eval(f_oop_expr)
c = [3., 5., 2., 6.]
f(c)
\[ \begin{equation} \left[ \begin{array}{c} 24 b - 93.75 k \\ 36 b - 281.25 k \\ - 12 b + 187.5 k \\ - 12 b + 93.75 k \\ \end{array} \right] \end{equation} \]
Alternatively, if we use a symbolic rate matrix, we could define our function to take in both species concentrations and parameter values as arguments:
K = fluxmat(rn)
odes = N * K * Φ
f_oop_expr, f_iip_expr = Symbolics.build_function(odes, species(rn), parameters(rn))
f = eval(f_oop_expr)
c = [3., 5., 2., 6]; ks = [12., 4.]
f(c, ks)
4-element Vector{Float64}:
-1029.0
-3231.0
2202.0
1077.0
Note also that f
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 = ZB$, where Z
is the complex stoichiometry matrix and B
is the incidence matrix of the graph.
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 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
Catalyst.netstoichmat
— Functionnetstoichmat(rn, sparse=false)
Returns the net stoichiometry matrix, $N$, with $N_{i j}$ the net stoichiometric coefficient of the ith species within the jth reaction.
Notes:
- Set sparse=true for a sparse matrix representation
- Caches the matrix internally within
rn
so subsequent calls are fast. - Note that constant species are not treated as reactants, but just components that modify the associated rate law. As such they do not contribute to the net stoichiometry matrix.
Catalyst.complexstoichmat
— Functioncomplexstoichmat(network::ReactionSystem; sparse=false)
Given a ReactionSystem
and vector of reaction complexes, return a matrix with positive entries of size number of species by number of complexes, where the non-zero positive entries in the kth column denote stoichiometric coefficients of the species participating in the kth reaction complex.
Notes:
- Set sparse=true for a sparse matrix representation
Catalyst.incidencemat
— Functionincidencemat(rn::ReactionSystem; sparse=false)
Calculate the incidence matrix of rn
, see reactioncomplexes
.
Notes:
- Is cached in
rn
so that future calls, assuming the same sparsity, will also be fast.
Catalyst.fluxmat
— Functionfluxmat(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
. Has the useful property that $\frac{dx}{dt} = S*K*Φ(x)$, where S is the netstoichmat
or net stoichiometry matrix and $Φ(x)$ is the massactionvector
. 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.
Catalyst.laplacianmat
— Functionlaplacianmat(rn::ReactionSystem, pmap=Dict(); sparse=false)
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
and $A_k$ is the negative of the graph Laplacian, and $Φ$ is the massactionvector
. $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.
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.
Catalyst.massactionvector
— Functionmassactionvector(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 $rac{dx}{dt} = Y A_k Φ(x)$, where $Y$ is the complexstoichmat
and $A_k$ is the negative of the laplacianmat
. 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. 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.