Weakly Reversible CF-Decompositions of Chemical Kinetic Systems

The reactant map $\rho:\mathscr{R}\to\mathscr{C}$ associates a reaction $y\to y^{\prime}$ to its reactant $y$. Denote by $n_{r}$ the number of reactant complexes. The map is surjective if and only if $n_{r}=n$. In this case, the network is cycle terminal. It is injective if and only if $n_{r}=r$. In this case, the network is non-branching. The reactant map $\rho$ induces a further useful mapping called reactions map. The reactions map $\rho^{\prime}$ associates the coordinate of reactant complexes to the coordinates of all reactions of which the complex is a reactant, i.e., $\rho^{\prime}:\mathbb{R}^{n}\to\mathbb{R}^{r}$ such that $f:\mathscr{C}\to\mathbb{R}$ mapped to $f\circ\rho$.

We associate the following important linear maps to a positive element $k\in\mathbb{R}^{r}$: the $k$-diagonal map ${\rm{diag}}(k)$ that maps $\omega_{i}$ to $k_{i}\omega_{i}$ for $i\in\mathscr{R}$, the $k$-incidence map $I_{k}$, which is defined as the composition ${\rm{diag}}(k)\circ\rho^{\prime}$, and the Laplacian map $A_{k}:\mathbb{R}^{n}\to\mathbb{R}^{n}$, which is given by $A_{k}=I_{a}\circ I_{k}$.

From the classic graph theory, a directed graph or simply digraph, denoted by $D$, consists of a non-empty finite set of vertices and a finite set ordered pairs called arcs. For an arc $(u,v)$, the first vertex $u$ is its tail and the second vertex $v$ is its head [4].

A trail is a walk such that all arcs are distinct. If the vertices of walk $W$ are distinct, $W$ is a called path. If the vertices $x_{1},x_{2},\ldots,x_{k-1}$ are distinct with $k\geq 3$ and $x_{1}=x_{k}$, the walk is called a cycle [4].

For this paper, we consider “cycles” to be simple cycles, i.e., no repeated vertices and no repeated arcs.

Chemical reaction networks are digraphs where complexes are vertices and reactions are arcs. If there is a path between two vertices $C_{i}$ and $C_{j}$, then they are said to be connected. If there is a directed path from vertex $C_{i}$ to vertex $C_{j}$ and vice versa, then they are said to be strongly connected. If any two vertices of a subgraph are (strongly) connected, then the subgraph is said to be a (strongly) connected component. The (strong) connected components are precisely the (strong) linkage classes of a chemical reaction network. The maximal strongly connected subgraphs where there are no arcs from a complex in the subgraph to a complex outside the subgraph is said to be the terminal strong linkage classes. We denote the number of linkage classes, the number of strong linkage classes, and the number of terminal strong linkage classes by $l,sl,$ and $t$, respectively. A chemical reaction network is said to be weakly reversible if $sl=l$, i.e., each linkage class is strongly connected, and it is said to be $t$-minimal if $t=l$.

The ODE or dynamical system of a chemical kinetics system is $\dfrac{{dx}}{{dt}}=f\left(x\right)$. An equilibrium or steady state is a zero of $f$.

A chemical reaction network is said to admit multiple equilibria if there exist positive rate constants such that the ODE system admits more than one stoichiometrically compatible equilibria. On the other hand, the set of complex balanced equilibria [21] is given by

A positive vector $c\in\mathbb{R}^{m}$ is complex balanced if $K\left(c\right)$ is contained in $\ker{I_{a}}$. A chemical kinetic system is complex balanced if it has a complex balanced equilibrium.

If the kinetic order values are the corresponding stoichiometric coefficients, then the system has the mass action kinetics (MAK).

We now consider the complex factorizable kinetics (CFK) and its complement, the non-complex factorizable kinetics (NFK) [1].

We have $K={\rm{diag}}\left(k\right)\circ\rho^{\prime}\circ{\Psi_{K}}$ via the definition of the $k$-incidence map $I_{k}$. This implies that a complex factorizable kinetics $K$ has a decomposition into diagonal matrix of rate constants and an interaction map $\rho^{\prime}\circ{\Psi_{K}}:\mathbb{R}_{\geq 0}^{m}\to\mathbb{R}_{\geq 0}^{r}$. Note that the values of the interaction map are “reaction-determined”. It basically means that they are determined by the values on the reactant complexes. CF kinetics (CFK) generalizes the key structural property of MAK. The complex formation rate function decomposes as $g={A_{k}}\circ{\Psi_{K}}$ while the species formation rate function factorizes as $f=Y\circ{A_{k}}\circ{\Psi_{K}}$. The CFK systems in the set of PLK systems are precisely the PL-RDK systems [1].

We now formally define the RID (rate constant-interaction map decomposable) kinetics.

Decomposition theory was initiated by M. Feinberg in 1987 in his review [8]. He introduced the following definition of a decomposition:

The use of the term “decomposition” is not consistent in the CRNT literature as some authors require only that $\mathscr{R}=\cup\mathscr{R}_{i}$, but not that the subsets are pairwise disjoint. We denote this looser concept as a covering of the network.

We denote a decomposition by $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup...\cup\mathscr{N}_{k}$ since $\mathscr{N}$ is a union of the subnetworks in the sense of [15]. It also follows that,

for the corresponding stoichiometric subspaces.

M. Feinberg also introduced the important concept of independence of decompositions:

A basic property of an independent decomposition was first recorded in [14]:

Decompositions can have the following useful relation:

Farinas et al. [7] noted the following relationship:

M. Feinberg in [8] demonstrated the importance of the independence concept through the following Theorem:

Note however that independence does not ensure that the intersection is non-empty though there is a positive equilibrium for each subnetwork. Additional properties of the decomposition and the kinetics are necessary to achieve this. To our knowledge, there are only a few such results.

The corresponding concepts and results for complex balanced equilibria were established only recently in [7]:

Clearly, incidence independence is equivalent to

Incidence independence has the following analogous properties to independence:

An important result in [7] is the analogue of Feinberg´s Theorem:

In contrast to the case of positive equilibria, there is a class of decompositions ensuring the converse of statement iii. To state this, we introduce the concept of the set $\mathscr{C}_{\mathscr{D}}$ of common complexes of a decomposition, first introduced in [11]:

The $\mathscr{C}$-decompositions are those which also partition the set of complexes. In [7], it is shown that any [7] decomposition is a coarsening of the linkage class decomposition. More importantly, the following converse to Statement iii. holds:

It is shown in [7] too, that any $\mathscr{C}^{*}$-decomposition is incidence independent. As observed in [11], for a decomposition with $|\mathscr{C}_{\mathscr{D}}|>1$, this is no longer true in general. For example, the simple network $R_{1}:X\to Y$, $R_{2}:Y\to X$ has the decomposition $\{R_{1}\}\cup\{R_{2}\}$ with $|\mathscr{C}_{\mathscr{D}}|=2$, $n-l=1<1+1=2$, and hence, not incidence independent.

The combination of the independence and incidence independence in a decomposition leads to several interesting properties:

A zero deficiency decomposition (ZDD) is a decomposition whose subnetworks have zero deficiency. In a deficiency zero network, we obtain the following interesting equivalences:

Despite its early founding by M. Feinberg in 1987, decomposition theory has received little attention in the CRNT community. Recently however, its usefulness beyond results on equilibria existence and parametrization, e.g., in the relatively new field of concentration robustness. Interesting results for larger and high deficiency systems have been derived for more general kinetic systems such as power law systems [12, 13] and Hill-type systems [17].

We refer to Figure 2. R. Schmitz’s model [27] of the earth’s carbon cycle in the pre-industrial state consists of six carbon pools – $M_{1}$(atmosphere), $M_{2}$(warm ocean surface water), $M_{3}$(cool ocean surface water), $M_{4}$ (deep ocean waters), $M_{5}$(terrestrial biota) and $M_{6}$(soil and detritus) – and 13 (directed) mass transfers among them. A power law representation of the model was first investigated by Fortun et al. in [14]: the weakly reversible system $(\mathscr{N},K)$ has zero deficiency since its complexes are monospecies and a PL-NDK system with 3 NF nodes $M_{1}$, $M_{2}$ and $M_{3}$ (s. Figure and kinetic order matrix). Table 3.1 lists the reaction sets and the CF-subsets of the NF nodes.

The system has a weakly reversible PL-RDK decomposition $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup\mathscr{N}_{3}$, with the subnetworks induced by the following subsets of the reaction set:

The decomposition is also clearly a zero deficiency decomposition (ZDD). Its set of $\mathscr{C}_{\mathscr{D}}$ of common complexes = $\{M_{1},M_{2},M_{3},M_{4}\}$.

The weakly reversible subnetwork $\mathscr{N}^{\prime}:=\mathscr{N}_{1}\cup\mathscr{N}_{2}$ is of special interest since the PL-RDK decomposition for it is also independent. Being also ZDD in a deficiency zero network, the decomposition is bi-independent. It is also a $\mathscr{C}^{*}$-decomposition. The kinetic order matrix of the subsystem consists of the first 8 rows of the matrix in Figure 2, showing that the subsystem is PL-NDK with one NF-node $M_{1}$ with 2 CF-subsets $\{r_{2}\}$ and $\{r_{5}\}$.

We introduce some concepts that are essential to this paper about digraphs provided in [4, 5].

Theorem 4.5 was taken from [4], which was attributed to Veblen, while Theorem 4.6 was taken from the version in [5], which was called the Euler’s Theorem. These results are useful to determine whether a digraph can be decomposed into cycles.

Since we are studying cycles originating in a reactant complex, it suffices to consider a connected digraph, i.e., a linkage class. We will call a linkage class an NF-linkage class if it contains an NF-reactant complex, otherwise a CF-linkage class.

According to Theorems 4.5 and 4.6, a digraph is Eulerian if and only if it has a decomposition into directed cycles. In the lucky scenario that a linkage class is Eulerian, we obtain a solution to our problem:

We will define concepts in general for a kinetic system, but clearly the interesting ones that we will focus on are the weakly reversible systems. Let $y$ be any complex in a weakly reversible linkage class $\mathscr{N}$. If $K$ is an RIDK kinetics, let $\{{\rm{CF}}_{i}(y)\}$ be the CF-subsets of $y$, which partition the set $\mathscr{R}(y)$ of $y$’s branching reactions, $i=1,\ldots,N_{R}(y)$. Note that $N_{R}(y)=1$ iff $y$ is a CF-node. Also, a simple cycle of $y$ contains only 1 branching reaction.

The covering is not necessarily a decomposition because an $i$-cycle and a $j$-cycle where $i\neq j$ can have reactions in common. An example of a non-trivial decomposition is given in Example 4.10.

We now introduce the following important necessary condition for a weakly reversible CF-decomposition of a linkage class with a single NF node. This result is general whether Eulerian or non-Eulerian.

The necessary condition also suggests the following procedure for generating a weakly reversible CF-decomposition in a single linkage class and single NF node RIDK system:

1. 1.

   Choose a CF-subset of the NF-node $y$, denote it with CF${}_{1}(y)$. For each branching reaction in CF₁, choose a simple cycle of $y$. Let $\mathscr{N}_{0,1}$ be the subnetwork generated by the set $\mathscr{R}_{0,1}$ of reactions of these cycles.

2. 2.

   If $\mathscr{R}\backslash\mathscr{R}_{0,1}\neq\varnothing$ and the subnetwork generated by $\mathscr{R}\backslash\mathscr{R}_{0,1}$ is not weakly reversible, then the system has no weakly reversible CF-decomposition.

3. 3.

   Otherwise, for another CF subset, CF${}_{2}(y)$, choose simple cycles which do not contain any reactions from $\mathscr{R}_{0,1}$. If there is a branching reaction in CF${}_{2}(y)$ for which this is not possible, then the system does not have a weakly reversible CF-decomposition. Otherwise, let $\mathscr{N}_{0,2}$ be the subnetwork generated by the set $\mathscr{R}_{0,2}$ of reactions of these cycles.

4. 4.

   If $\mathscr{R}\backslash(\mathscr{R}_{0,1}\cup\mathscr{R}_{0,2})\neq\varnothing$ and the subnetwork generated by $\mathscr{R}\backslash(\mathscr{R}_{0,1}\cup\mathscr{R}_{0,2})$ is not weakly reversible, then the system does not have a weakly reversible CF-decomposition.

5. 5.

   Continue accordingly until all CF subsets for $1,\ldots,N_{R}(y)$ are covered. If the complement of $\cup\mathscr{R}_{0,i}$ is not empty, find for each reaction a cycle of $y$ with the property that if the cycle’s branching reaction belongs to CF_j, no reactions from $\cup\mathscr{R}_{0,i}$, $i\neq j$ are contained in the cycle. If there is an instance for which this is not possible, then the system does not have a weakly reversible CF-decomposition. Otherwise, add the cycles to the subnetworks with the same CF_j.

The last step of the procedure defines an appropriate weakly reversible CF-decomposition with $N_{R}(y)$ subnetworks, $y$ in each subnetwork and by construction, the reactions of each subnetwork as those of the $i$-cycles contained in the subnetwork. Because of the many choices made in the procedure, this decomposition is not unique. We provide Algorithm 1 for this procedure.