Fundamental Decompositions and Multistationarity of Power-Law Kinetic Systems

TL;DR

This paper investigates the properties of the fundamental decomposition of chemical reaction networks, particularly its independence and incidence-independence, and applies these findings to enhance the multistationarity analysis of power-law kinetic systems.

Contribution

It provides necessary and sufficient conditions for the properties of the $\\mathscr{F}$-decomposition and improves the multistationarity algorithm for power-law systems by removing the need for certain transformations.

Findings

Derived conditions for independence and incidence-independence of the $\\mathscr{F}$-decomposition.

Identified network classes where the $\\mathscr{F}$-decomposition aligns with other decompositions.

Improved the multistationarity algorithm by eliminating the transformation step for certain systems.

Abstract

The fundamental decomposition of a chemical reaction network (also called its "$\mathscr{F}$-decomposition") is the set of subnetworks generated by the partition of its set of reactions into the "fundamental classes" introduced by Ji and Feinberg in 2011 as the basis of their "higher deficiency algorithm" for mass action systems. The first part of this paper studies the properties of the $\mathscr{F}$-decomposition, in particular, its independence (i.e., the network's stoichiometric subspace is the direct sum of the subnetworks' stoichiometric subspaces) and its incidence-independence (i.e., the image of the network's incidence map is the direct sum of the incidence maps' images of the subnetworks). We derive necessary and sufficient conditions for these properties and identify network classes where the $\mathscr{F}$-decomposition coincides with other known decompositions. The second part of the paper applies the above-mentioned results to improve the Multistationarity Algorithm for power-law kinetic systems (MSA), a general computational approach that we introduced in previous work. We show that for systems with non-reactant determined interactions but with an independent $\mathscr{F}$-decomposition, the transformation to a dynamically equivalent system with reactant-determined interactions -- required in the original MSA -- is not necessary. We illustrate this improvement with the subnetwork of Schmitz's carbon cycle model recently analyzed by Fortun et al.

Equations

$${\left({{I_{a}}}\right)_{ij}}=\left\{\begin{array}[]{rl}-1&{\rm{if\ }}{C_{i}}{\rm{\ is\ in\ the\ reactant\ complex\ of\ reaction\ }}{R_{j}},\\ 1&{\rm{if\ }}{C_{i}}{\rm{\ is\ in\ the\ product\ complex\ of\ reaction\ }}{R_{j}},\\ 0&{\rm{otherwise}}.\end{array}\right.$$

$$\delta \le {\delta }_1+{\delta }_2...+{\delta }_k.$$

$$\text{Im }{I}_{a,i}=I_a({\mathbb{R}}^{{\mathcal{R}}_i}).$$

$$E_{+}\left({\mathcal{N}}_1,K_1\right)\cap E_{+}\left({\mathcal{N}}_2,K_2\right)\cap ...\cap E_{+}\left({\mathcal{N}}_k,K_k\right)\subseteq E_{+}\left(\mathcal{N},K\right).$$

$$L_{\mathcal{O}}(\alpha )=\sum _{y\to y^{\prime }\in \mathcal{O}}{\alpha }_{y\to y^{\prime }}\left(y^{\prime }-y\right).$$

$$\begin{array}[]{lll}R_{1}:M_{1}\to M_{5}&&R_{5}:M_{1}\to M_{3}\\ R_{2}:M_{5}\to M_{1}&&R_{6}:M_{3}\to M_{4}\\ R_{3}:M_{5}\to M_{6}&&R_{7}:M_{4}\to M_{2}\\ R_{4}:M_{6}\to M_{1}&&R_{8}:M_{2}\to M_{1}\\ \end{array}$$

$$\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left(\kern 0.0pt\kern-2.5pt\kern-6.0pt\vbox{\kern-0.86108pt\vbox{\vbox{ \halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&& \kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr 5.0pt\hfil\@arstrut$\displaystyle~{}$\hfil\kern 5.0pt&5.0pt\hfil$\displaystyle v_{1}$\kern 5.0pt&5.0pt\hfil$\displaystyle v_{2}$\kern 5.0pt\crcr\kern 0.86108pt\cr 5.0pt\hfil\@arstrut$\displaystyle R_{1}$\hfil\kern 5.0pt&5.0pt\hfil$0$\kern 5.0pt&5.0pt\hfil$1$\kern 5.0pt\cr 5.0pt\hfil\@arstrut$\displaystyle R_{3}$\hfil\kern 5.0pt&5.0pt\hfil$0$\kern 5.0pt&5.0pt\hfil$1$\kern 5.0pt\cr 5.0pt\hfil\@arstrut$\displaystyle R_{4}$\hfil\kern 5.0pt&5.0pt\hfil$0$\kern 5.0pt&5.0pt\hfil$1$\kern 5.0pt\cr 5.0pt\hfil\@arstrut$\displaystyle R_{5}$\hfil\kern 5.0pt&5.0pt\hfil$1$\kern 5.0pt&5.0pt\hfil$0$\kern 5.0pt\cr 5.0pt\hfil\@arstrut$\displaystyle R_{6}$\hfil\kern 5.0pt&5.0pt\hfil$1$\kern 5.0pt&5.0pt\hfil$0$\kern 5.0pt\cr 5.0pt\hfil\@arstrut$\displaystyle R_{7}$\hfil\kern 5.0pt&5.0pt\hfil$1$\kern 5.0pt&5.0pt\hfil$0$\kern 5.0pt\cr 5.0pt\hfil\@arstrut$\displaystyle R_{8}$\hfil\kern 5.0pt&5.0pt\hfil$1$\kern 5.0pt&5.0pt\hfil$0$\kern 5.0pt\cr}}}}\right)$}}.$$

$$\left(\begin{array}{ccccc}{F}_0& & & & \\ & F_1& & \text{0}& \\ & & \ddots & & \\ & \text{0}& & \ddots & \\ & & & & F_w\end{array}\right)$$

$$\begin{array}[]{c}{S_{0}}+K\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{0}}K\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{1}}K{\mkern 1.0mu}\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}...{\mkern 1.0mu}\to{\mkern 1.0mu}{S_{k-1}}K\to{S_{k}}+K\\ {S_{k}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{k}}F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{\mkern 1.0mu}...\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{2}}F{\mkern 1.0mu}\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{1}}F\to{S_{0}}+F\end{array}$$

$$\begin{array}[]{c}{\alpha_{k}}\left({{S_{k}}+K-{S_{k-1}}K}\right)+{\alpha_{k-1}}\left({{S_{k-1}}K-{S_{k-2}}K}\right)+...+{\alpha_{1}}\left({{S_{1}}K-{S_{0}}K}\right)\\ +{\alpha_{0}}\left({{S_{0}}K-\left({{S_{0}}+K}\right)}\right)+{\beta_{0}}\left({{S_{0}}+F-{S_{1}}F}\right)+{\beta_{1}}\left({{S_{1}}F-{S_{2}}F}\right)\\ +...+{\beta_{k}}\left({{S_{k}}F-\left({{S_{k}}+F}\right)}\right)=0.\end{array}$$

$$\begin{array}[]{c}{S_{k}}\left({{\alpha_{k}}-{\beta_{k}}}\right)+{S_{k-1}}K\left({-{\alpha_{k}}+{\alpha_{k-1}}}\right)+...+{S_{0}}K\left({-{\alpha_{1}}+{\alpha_{0}}}\right)+{S_{0}}\left({-{\alpha_{0}}+{\beta_{0}}}\right)\\ +{S_{1}}F\left({-{\beta_{0}}+{\beta_{1}}}\right)+...+{S_{k}}F\left({-{\beta_{k-1}}+{\beta_{k}}}\right)+K\left({{\alpha_{k}}-{\alpha_{0}}}\right)+F\left({{\beta_{0}}-{\beta_{k}}}\right)=0.\end{array}$$

$$\begin{array}[]{c}{S_{0}}+K\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{0}}K{\mkern 1.0mu}\to{S_{1}}+K{\mkern 1.0mu}\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{1}}K\to{S_{2}}+K{\mkern 1.0mu}\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}...\to{S_{k}}+K\\ {S_{k}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}...{\mkern 1.0mu}{\mkern 1.0mu}\to{\mkern 1.0mu}{S_{2}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{2}}F\to{S_{1}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{1}}F\to{S_{0}}+F\end{array}$$

$$\begin{array}[]{*{20}{c}}{{S_{i}}+K\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{i}}K\to{S_{i+1}}+K}\\ {{S_{i+1}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{i+1}}F\to{S_{i}}+F}\end{array}{\text{ \ \ for }}i=0,1,...,k-1,{\text{ and }}$$

$$\begin{array}[]{*{20}{c}}{{\alpha_{0}}\left({{S_{0}}K-\left({{S_{0}}+K}\right)}\right)+{\beta_{0}}\left({{S_{1}}+K-{S_{0}}K}\right)+{\alpha_{1}}\left({{S_{1}}K-\left({{S_{1}}+K}\right)}\right)+{\beta_{1}}\left({{S_{2}}+K-{S_{1}}K}\right)}\\ {+...+{\alpha_{k-1}}\left({{S_{k-1}}K-\left({{S_{k-1}}+K}\right)}\right)+{\beta_{k-1}}\left({{S_{k}}+K-{S_{k-1}}K}\right)+}\\ {{\lambda_{k-1}}\left({{S_{k}}F-\left({{S_{k}}+F}\right)}\right)+{\gamma_{k-1}}\left({{S_{k-1}}+F-{S_{k}}F}\right)+...+}\\ {{\lambda_{1}}\left({{S_{2}}F-\left({{S_{2}}+F}\right)}\right)+{\gamma_{1}}\left({{S_{1}}+F-{S_{2}}F}\right)+{\lambda_{0}}\left({{S_{1}}F-\left({{S_{1}}+F}\right)}\right)+{\gamma_{0}}\left({{S_{0}}+F-{S_{1}}F}\right)=0.}\end{array}$$

$$\begin{array}[]{c}{S_{0}}K\left({{\alpha_{0}}-{\beta_{0}}}\right)+{S_{0}}\left({-{\alpha_{0}}+{\gamma_{0}}}\right)+{S_{1}}K\left({{\alpha_{1}}-{\beta_{1}}}\right)+...+{S_{k-1}}K\left({{\alpha_{k-1}}-{\beta_{k-1}}}\right)+\\ {S_{k}}F\left({{\lambda_{k-1}}-{\gamma_{k-1}}}\right)+...+{S_{2}}F\left({{\lambda_{1}}-{\gamma_{1}}}\right)+{S_{1}}F\left({{\lambda_{0}}-{\gamma_{0}}}\right)+\\ K\left({-{\alpha_{0}}+{\beta_{0}}-{\alpha_{1}}+{\beta_{1}}-...-{\alpha_{k-1}}+{\beta_{k-1}}}\right)+F\left({-{\lambda_{k-1}}+{\gamma_{k-1}}-...-{\lambda_{1}}+{\gamma_{1}}-{\lambda_{0}}+{\gamma_{0}}}\right)\\ +{S_{1}}\left({{\beta_{0}}-{\alpha_{1}}+{\gamma_{1}}-{\lambda_{0}}}\right)+{S_{2}}\left({{\beta_{1}}-{\alpha_{2}}+{\gamma_{2}}-{\lambda_{1}}}\right)+...+{S_{k}}\left({{\beta_{k-1}}-{\lambda_{k-1}}}\right)=0.\end{array}$$

$$\begin{array}[]{c}{S_{00}}+K\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{00}}K\to{S_{01}}K\to{S_{11}}+K{\rm{\ \ \ \ \ }}{S_{11}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{11}}F\to{S_{10}}F\to{S_{00}}+F\\ {S_{01}}K\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{01}}+K{\rm{\ \ \ \ \ }}{S_{10}}F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{10}}+F\\ {S_{10}}+K\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{10}}K\to{S_{11}}+K{\rm{\ \ \ \ \ }}{S_{01}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{01}}F\to{S_{00}}+F\end{array}$$

$$\begin{array}[]{c}{S_{0}}+K\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{0}}K\to{S_{1}}K\to{S_{2}}+K\\ {S_{2}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{2}}F\to{S_{1}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{1}}F\to{S_{0}}+F\end{array}$$

$$\begin{array}[]{lll}R_{1}:A_{1}+2A_{2}\to 2A_{1}+A_{2}&&R_{6}:A_{1}+2A_{4}\to 2A_{1}+A_{4}\\ R_{2}:A_{1}+A_{2}\to 2A_{2}&&R_{7}:A_{1}+A_{4}\to 2A_{4}\\ R_{3}:A_{2}\to A_{3}&&R_{8}:A_{4}\to A_{3}\\ R_{4}:A_{3}\to A_{2}&&R_{9}:A_{3}\to A_{4}\\ R_{5}:A_{4}+A_{5}\to 2A_{4}&&R_{10}:A_{1}+A_{2}+A_{4}\to A_{5}+A_{2}+A_{4}\\ \end{array}$$

$$\begin{array}[]{c}{S_{0}}+K\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{0}}K\to{S_{1}}K\to{S_{2}}+K\\ {S_{2}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{2}}F\to{S_{1}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{1}}F\to{S_{0}}+F\end{array}$$

$$\begin{array}[]{c}{S_{0}}+K\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{0}}K{\mkern 1.0mu}\to{S_{1}}+K{\mkern 1.0mu}\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{1}}K\to{S_{2}}+K{\mkern 1.0mu}\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}...\to{S_{k}}+K\\ {S_{k}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}...{\mkern 1.0mu}{\mkern 1.0mu}\to{\mkern 1.0mu}{S_{2}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{2}}F\to{S_{1}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{1}}F\to{S_{0}}+F\end{array}$$

$$\begin{array}[]{*{20}{c}}{{S_{i}}+K\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{i}}K\to{S_{i+1}}+K}\\ {{S_{i+1}}+F\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{S_{i+1}}F\to{S_{i}}+F}\end{array}{\text{ for }}i=0,1,...,k-1.$$

$$\sum \left(n_i-l_i\right)=\sum \left(6-2\right)=4k=\left(4k+2\right)-2=n-l.$$

$$\begin{array}[]{c}X\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}XT\to{X_{p}}\\ {X_{p}}+Y\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}{X_{p}}Y\to X+{Y_{p}}\\ XT+{Y_{p}}\mathbin{\lower 1.1625pt\hbox{$\buildrel\textstyle\rightarrow\over{\smash{\leftarrow}\vphantom{}}$}}XT{Y_{p}}\to XT+Y\end{array}$$

$$\begin{array}[]{l}R_{i}:x_{i}\to 2x_{i}\\ R_{-i}:x_{i}\to 0\\ R_{-i}^{\prime}:2x_{i}\to x_{i}\end{array}$$

$${\alpha }_1\left(C_2-C_1\right)+{\alpha }_2\left(C_3-C_2\right)+...+{\alpha }_{{\ell }_1}\left(C_1-C_{{\ell }_1}\right)+$$

$${\alpha }_{{\ell }_1+1}\left(C_{{\ell }_1+2}-C_{{\ell }_1+1}\right)+{\alpha }_{{\ell }_1+2}\left(C_{{\ell }_1+3}-C_{{\ell }_1+2}\right)+...+{\alpha }_{{\ell }_2}\left(C_{{\ell }_1+1}-C_{{\ell }_2}\right)+...+$$

$${\alpha }_{{\ell }_{k-1}+1}\left(C_{{\ell }_{k-1}+2}-C_1\right)+{\alpha }_{{\ell }_{k-1}+2}\left(C_{{\ell }_{k-1}+3}-C_{{\ell }_{{}_{k-1}}+2}\right)+...+{\alpha }_{{\ell }_k}\left(C_{{\ell }_{{}_{k-1}}+1}-C_{{\ell }_k}\right)=0.$$

$$C_1\left({\alpha }_{{\ell }_1}-{\alpha }_1\right)+C_2\left({\alpha }_1-{\alpha }_2\right)+C_3\left({\alpha }_2-{\alpha }_3\right)+...+C_{{\ell }_1}\left({\alpha }_{{\ell }_1-1}-{\alpha }_{{\ell }_1}\right)+$$

$$C_{{\ell }_1+1}\left({\alpha }_{{\ell }_2}-{\alpha }_{{\ell }_1+1}\right)+C_{{\ell }_1+2}\left({\alpha }_{{\ell }_1+1}-{\alpha }_{{\ell }_1+2}\right)+...+C_{{\ell }_2}\left({\alpha }_{{\ell }_2-1}-{\alpha }_{{\ell }_2}\right)+...+$$

$$C_{{\ell }_{{}_{k-1}}+1}\left({\alpha }_{{\ell }_k}-{\alpha }_{{\ell }_{k-1}+1}\right)+C_{{\ell }_{k-1}+2}\left({\alpha }_{{\ell }_{k-1}+1}-{\alpha }_{{\ell }_{k-1}+2}\right)+...+C_{{\ell }_k}\left({\alpha }_{{\ell }_k-1}-{\alpha }_{{\ell }_k}\right)=0.$$

Full text

1 Introduction

The fundamental decomposition (also called “$\mathscr{F}$-decomposition”) of a chemical reaction network (CRN) is the set of subnetworks generated by the partition of its set of reactions into the “fundamental classes” introduced by Ji and Feinberg in 2011 as the basis of their Higher Deficiency Algorithm (HDA) for mass action systems. For a CRN with only irreversible reactions and stoichiometric matrix $N$, the characteristic functions $\omega_{r}$ and $\omega_{r^{\prime}}$ of reactions $r$ and $r^{\prime}$ are in the same non-zero fundamental class $C_{i}$ if they are non-zero and pairwise dependent in the factor space ${\mathbb{R}^{\mathscr{R}}/({\text{Ker\ }}N)^{\perp}}$. Any reaction with $\omega_{r}=0$ in the factor space is assigned to the zero fundamental class $C_{0}$. In the general case, a reversible reaction pair is also assigned to the same fundamental class.

The first part of this paper studies the properties of the $\mathscr{F}$-decomposition, in particular, its independence (i.e., the network’s stoichiometric subspace is the direct sum of the subnetworks’ stoichiometric subspaces) and its incidence-independence (i.e., the image of the network’s incidence map is the direct sum of the incidence maps’ images of the subnetworks). For deficiency zero networks, these two properties coincide. M. Feinberg established the essential relationship between independent decompositions and the set of positive equilibria of a network (recalled as Theorem 2.21) in 1987 [7]. A corresponding relationship between incidence-independent, weakly reversible decompositions and complex-balanced equilibria of a weakly reversible network (recalled as Theorem 2.22) was recently documented by Farinas et al. [6]. We derive necessary and sufficient conditions for the independence and incidence of $\mathscr{F}$-decompositions, and identify network classes where the $\mathscr{F}$-decomposition coincides with other known decompositions.

In our previous work [13], we showed that the HDA of Ji and Feinberg can be extended to any power-law kinetic system which has reactant-determined interactions (we denote the set by PL-RDK), i.e., the reactions branching from the same reactant complex have identical kinetic order vectors (or “interactions”). By combining this extension with a method (called CF-RM) to transform a power-law kinetic system with non-reactant-determined interactions (denoted by PL-NDK) to a dynamically equivalent PL-RDK system, we developed a computational approach called the “Multistationarity Algorithm” (MSA) to determine multistationarity in a stoichiometric class for any power-law kinetic system.

The second part of the paper applies the results of the first part to improve the MSA for power-law kinetic systems. We show that for PL-NDK systems with an independent $\mathscr{F}$-decomposition, the extended HDA, can be directly applied just as in the case of a PL-RDK system. This is because the new transformation, denoted by CF-RI+, preserves reversibility/ irreversibility of reactions, leading to identical multistationarity computations for such PL-NDK systems and their PL-RDK transforms via CF-RI+. The subnetwork of Schmitz’s global carbon cycle model, a PL-NDK system studied by Fortun et al. [11], is used as a running example to illustrate this improvement of the MSA.

The paper is organized as follows: Section 2 collects the fundamentals of chemical reaction networks and kinetic systems required for the later sections, including relevant results from decomposition theory. After introducing the $\mathscr{F}$-decomposition and related constructs, Section 3 derives its basic properties, including bounds for the number of subnetworks and related necessary conditions for independence and incidence-independence. Network classes with independent or incidence-independent $\mathscr{F}$-decompositions are identified. In Section 4, Ji and Feinberg’s characterization of $\mathscr{F}$-decomposition subnetworks is used to develop a classification of $\mathscr{F}$-decompositions. Bounds for the deficiency and other properties are then obtained for the three $\mathscr{F}$-decomposition types. The reaction reversibility/ irreversibility preserving transformation CF-RI+ is introduced in Section 5. This forms the basis for the improvement of the Multistationarity Algorithm (MSA) derived in Section 6. Conclusions and an outlook constitute Section 7. Tables of acronyms and frequently used symbols are provided in Appendix A.

2 Fundamentals of Chemical Reaction Networks and Kinetic Systems

In this section, we recall some fundamental notions about chemical reaction networks and chemical kinetic systems. These concepts are provided in [1, 8]. Moreover, we present some important preliminaries on the decomposition theory which was introduced by Feinberg in [7].

2.1 Fundamentals of Chemical Reaction Networks

Definition 2.1**.**

A chemical reaction network (CRN) $\mathscr{N}$ is a triple $\left(\mathscr{S},\mathscr{C},\mathscr{R}\right)$ of nonempty finite sets where $\mathscr{S}$, $\mathscr{C}$, and $\mathscr{R}$ are the sets of $m$ species, $n$ complexes, and $r$ reactions, respectively, such that $\left({{C_{i}},{C_{i}}}\right)\notin\mathscr{R}$ for each $C_{i}\in\mathscr{C}$; and for each $C_{i}\in\mathscr{C}$, there exists $C_{j}\in\mathscr{C}$ such that $\left({{C_{i}},{C_{j}}}\right)\in\mathscr{R}$ or $\left({{C_{j}},{C_{i}}}\right)\in\mathscr{R}$.

Definition 2.2**.**

The molecularity matrix, denoted by $Y$, is an $m\times n$ matrix such that $Y_{ij}$ is the stoichiometric coefficient of species $X_{i}$ in complex $C_{j}$. The incidence matrix, denoted by $I_{a}$, is an $n\times r$ matrix such that

[TABLE]

The stoichiometric matrix, denoted by $N$, is the $m\times r$ matrix given by $N=YI_{a}$.

Let $\mathscr{I}=\mathscr{S},\mathscr{C}$ or $\mathscr{R}$. We denote the standard basis for $\mathbb{R}^{\mathscr{I}}$ by $\left\{\omega_{i}\in\mathbb{R}^{\mathscr{I}}\mid i\in\mathscr{I}\right\}$.

Definition 2.3**.**

Let $\mathscr{N}=(\mathscr{S,C,R})$ be a CRN. The incidence map $I_{a}:\mathbb{R}^{\mathscr{R}}\rightarrow\mathbb{R}^{\mathscr{C}}$ is the linear map such that for each reaction $r:C_{i}\rightarrow C_{j}\in\mathscr{R}$, the basis vector $\omega_{r}$ to the vector $\omega_{C_{j}}-\omega_{C_{i}}\in\mathscr{C}$.

Definition 2.4**.**

The reaction vectors for a given reaction network $\left(\mathscr{S},\mathscr{C},\mathscr{R}\right)$ are the elements of the set $\left\{{C_{j}}-{C_{i}}\in\mathbb{R}^{\mathscr{S}}|\left({{C_{i}},{C_{j}}}\right)\in\mathscr{R}\right\}.$

Definition 2.5**.**

The stoichiometric subspace of a reaction network $\left(\mathscr{S},\mathscr{C},\mathscr{R}\right)$, denoted by $S$, is the linear subspace of $\mathbb{R}^{\mathscr{S}}$ given by $S=span\left\{{{C_{j}}-{C_{i}}\in\mathbb{R}^{\mathscr{S}}|\left({{C_{i}},{C_{j}}}\right)\in\mathscr{R}}\right\}.$ The rank of the network, denoted by $s$, is given by $s=\dim S$. The set $\left({x+S}\right)\cap\mathbb{R}_{\geq 0}^{\mathscr{S}}$ is said to be a stoichiometric compatibility class of $x\in\mathbb{R}_{\geq 0}^{\mathscr{S}}$.

Definition 2.6**.**

Two vectors $x,x^{*}\in{\mathbb{R}^{\mathscr{S}}}$ are stoichiometrically compatible if $x-x^{*}$ is an element of the stoichiometric subspace $S$.

We can view complexes as vertices and reactions as edges. With this, CRNs can be seen as graphs. At this point, if we are talking about geometric properties, vertices are complexes and edges are reactions. 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 CRN. The maximal strongly connected subgraphs where there are no edges 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 and the number of strong linkage classes by $l$ and $sl$, respectively. A CRN is said to be weakly reversible if $sl=l$.

Definition 2.7**.**

For a CRN, the deficiency is given by $\delta=n-l-s$ where $n$ is the number of complexes, $l$ is the number of linkage classes, and $s$ is the dimension of the stoichiometric subspace $S$.

2.2 Fundamentals of Chemical Kinetic Systems

Definition 2.8**.**

A kinetics $K$ for a reaction network $\left(\mathscr{S},\mathscr{C},\mathscr{R}\right)$ is an assignment to each reaction $r:y\to y^{\prime}\in\mathscr{R}$ of a rate function ${K_{r}}:{\Omega_{K}}\to{\mathbb{R}_{\geq 0}}$ such that $\mathbb{R}_{>0}^{\mathscr{S}}\subseteq{\Omega_{K}}\subseteq\mathbb{R}_{\geq 0}^{\mathscr{S}}$, $c\wedge d\in{\Omega_{K}}$ if $c,d\in{\Omega_{K}}$, and ${K_{r}}\left(c\right)\geq 0$ for each $c\in{\Omega_{K}}$. Furthermore, it satisfies the positivity property: supp $y$ $\subset$ supp $c$ if and only if $K_{r}(c)>0$. The system $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ is called a chemical kinetic system.

Definition 2.9**.**

The species formation rate function (SFRF) of a chemical kinetic system is given by $f\left(x\right)=NK(x)=\displaystyle\sum\limits_{{C_{i}}\to{C_{j}}\in\mathscr{R}}{{K_{{C_{i}}\to{C_{j}}}}\left(x\right)\left({{C_{j}}-{C_{i}}}\right)}.$

The ordinary differential equation (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$.

Definition 2.10**.**

The set of positive equilibria of a chemical kinetic system $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ is given by ${E_{+}}\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)=\left\{{x\in\mathbb{R}^{\mathscr{S}}_{>0}|f\left(x\right)=0}\right\}.$

A CRN is said to admit multiple equilibria if there exist positive rate constants such that the ODE system admits more than one stoichiometrically compatible equilibria.

Definition 2.11**.**

A kinetics $K$ is complex factorizable if, for $K(x)=kI_{K}(x)$, the interaction map $I_{K}:\mathbb{R}^{\mathscr{S}}\to\mathbb{R}^{\mathscr{R}}$ factorizes via the space of complexes $\mathbb{R}^{\mathscr{C}}:I_{K}=I_{k}\circ\psi_{K}$ with $\psi_{K}:\mathbb{R}^{\mathscr{S}}\to\mathbb{R}^{\mathscr{C}}$ as factor map and $I_{k}=diag(k)\circ\rho^{\prime}$ with $\rho^{\prime}:\mathbb{R}^{\mathscr{C}}\to\mathbb{R}^{\mathscr{R}}$ assigning the value at a reactant complex to all its reactions.

Definition 2.12**.**

A kinetics $K$ is a power-law kinetics (PLK) if ${K_{i}}\left(x\right)={k_{i}}{{x^{{F_{i}}}}}$ for $i=1,...,r$ where ${k_{i}}\in{\mathbb{R}_{>0}}$ and ${F_{ij}}\in{\mathbb{R}}$. The power-law kinetics is defined by an $r\times m$ matrix $F$, called the kinetic order matrix and a vector $k\in\mathbb{R}^{\mathscr{R}}$, called the rate vector.

If the kinetic order matrix is the transpose of the molecularity matrix, then the system becomes the well-known mass action kinetics (MAK).

Definition 2.13**.**

A PLK system has reactant-determined kinetics (of type PL-RDK) if for any two reactions $i,j$ with identical reactant complexes, the corresponding rows of kinetic orders in $F$ are identical, i.e., ${f_{ik}}={f_{jk}}$ for $k=1,2,...,m$. A PLK system has non-reactant-determined kinetics (of type PL-NDK) if there exist two reactions with the same reactant complexes whose corresponding rows in $F$ are not identical.

2.3 Review of Decomposition Theory

In this subsection, we recall some definitions and earlier results from the decomposition theory of chemical reaction networks.

Definition 2.14**.**

A decomposition of $\mathscr{N}$ is a set of subnetworks $\{\mathscr{N}_{1},\mathscr{N}_{2},...,\mathscr{N}_{k}\}$ of $\mathscr{N}$ induced by a partition $\{\mathscr{R}_{1},\mathscr{R}_{2},...,\mathscr{R}_{k}\}$ of its reaction set $\mathscr{R}$.

We denote a decomposition with $\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 [12]. It also follows immediately that, for the corresponding stoichiometric subspaces, ${S}={S}_{1}+{S}_{2}+...+{S}_{k}$. It is also useful to consider refinements and coarsenings of decompositions.

Definition 2.15**.**

A network decomposition $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup...\cup\mathscr{N}_{k}$ is a refinement of $\mathscr{N}={\mathscr{N}^{\prime}}_{1}\cup{\mathscr{N}^{\prime}}_{2}\cup...\cup{\mathscr{N}^{\prime}}_{k^{\prime}}$ (and the latter a coarsening of the former) if it is induced by a refinement $\{\mathscr{R}_{1},\mathscr{R}_{2},...,\mathscr{R}_{k}\}$ of $\{{\mathscr{R}^{\prime}}_{1}\cup{\mathscr{R}^{\prime}}_{2}\cup...\cup{\mathscr{R}^{\prime}}_{k^{\prime}}\}$, i.e., each ${\mathscr{R}}_{i}$ is contained in an ${\mathscr{R}^{\prime}}_{j}$.

In [7], Feinberg introduced the important concept of independent decomposition.

Definition 2.16**.**

A network decomposition $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup...\cup\mathscr{N}_{k}$ is independent if its stoichiometric subspace is a direct sum of the subnetwork stoichiometric subspaces.

In Lemma 1 of [11], it was shown that for an independent decomposition,

[TABLE]

Definition 2.17**.**

A decomposition $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup...\cup\mathscr{N}_{k}$ with $\mathscr{N}_{i}=(\mathscr{S}_{i},\mathscr{C}_{i},\mathscr{R}_{i})$ is a $\mathscr{C}$-decomposition* if for each pair of distinct $i$ and $j$, $\mathscr{C}_{i}$ and $\mathscr{C}_{j}$ are disjoint.*

Example 2.18**.**

Linkage classes form the primary example of a decomposition of a CRN. They are special in the sense that the reaction set partition inducing them is also a complex set partition. In [6], the term “$\mathscr{C}$-decomposition” was introduced for such a decomposition and it was shown that any $\mathscr{C}$-decomposition is a coarsening of the linkage class decomposition.

For any decomposition, it also holds that $\text{Im }I_{a}=\text{Im }I_{a,1}+...+\text{Im }I_{a,k}$, where

[TABLE]

Definition 2.19**.**

A decomposition is incidence-independent if ${\rm{Im\ }}I_{a}$ is a direct sum of the ${\rm{Im\ }}I_{a,i}$. It is bi-independent if it is both independent and incidence-independent.

An equivalent formulation of showing incidence-independent is to satisfy $n-l=\sum{\left({{n_{i}}-{l_{i}}}\right)}$, where $n_{i}$ is the number of complexes and $l_{i}$ is the number of linkage classes, in each subnetwork $i$.

In [6], it was shown that for any incidence-independent decomposition, $\delta\geq\delta_{1}+\delta_{2}...+\delta_{k}$ and that $\mathscr{C}$-decompositions form a subset of incidence-independent decompositions. Hence, independent linkage classes form the primary example of a bi-independent decomposition.

The following proposition is easily verified:

Proposition 2.20**.**

A decomposition $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup...\cup\mathscr{N}_{k}$ is independent or incidence-independent and $\sum{{\delta_{i}}=\delta}$ iff $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup...\cup\mathscr{N}_{k}$ is bi-independent.

In particular, for a zero deficiency network, independence and incidence-independence are equivalent.

Feinberg established the following basic relation between an independent decomposition and the set of positive equilibria of a kinetics on the network:

Theorem 2.21**.**

*(Feinberg Decomposition Theorem [7])

Let $P(\mathscr{R})=\{\mathscr{R}_{1},\mathscr{R}_{2},...,\mathscr{R}_{k}\}$ be a partition of a CRN $\mathscr{N}$ and let $K$ be a kinetics on $\mathscr{N}$. If $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup...\cup\mathscr{N}_{k}$ is the network decomposition of $P(\mathscr{R})$ and ${E_{+}}\left(\mathscr{N}_{i},{K}_{i}\right)=\left\{{x\in\mathbb{R}^{\mathscr{S}}_{>0}|N_{i}K_{i}(x)=0}\right\}$ then*

[TABLE]

If the network decomposition is independent, then equality holds.

The analogue of Feinberg’s 1987 result for incidence-independent decompositions and complex-balanced equilibria is shown in [6]:

Theorem 2.22**.**

Let $\mathscr{N}$ be a network, $K$ any kinetics and $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup...\cup\mathscr{N}_{k}$ an incidence-independent decomposition of weakly reversible subnetworks. Then $\mathscr{N}$ is weakly reversible and

• i.

${Z_{+}}\left({\mathscr{N},K}\right)=\cap{Z_{+}}\left({{\mathscr{N}_{i}},K}\right)$* for each subnetwork ${\mathscr{N}_{i}}$.*

• ii.

If ${Z_{+}}\left({\mathscr{N},K}\right)\neq\emptyset$ then ${Z_{+}}\left({{\mathscr{N}_{i}},K}\right)\neq\emptyset$ for each subnetwork ${\mathscr{N}_{i}}$.

• iii.

If the decomposition is a $\mathscr{C}$-decomposition and $K$ a complex factorizable kinetics then ${Z_{+}}\left({{\mathscr{N}_{i}},K}\right)\neq\emptyset$ for each subnetwork ${\mathscr{N}_{i}}$ implies that ${Z_{+}}\left({{\mathscr{N}},K}\right)\neq\emptyset$.

Definition 2.23**.**

The network decomposition $\mathscr{N}=\mathscr{N^{\prime}}\cup\mathscr{N^{\prime\prime}}$ is said to be trivial if $\mathscr{N^{\prime}}$ is a subnetwork whose stoichiometric subspace coincides with that of $\mathscr{N}$.

3 Orientations, and $\mathscr{O}$-, $\mathscr{P}$-, and $\mathscr{F}$-decompositions of CRNs

In this section, we review the concepts and properties underlying HDA and its extension to PL-RDK systems in the context of decomposition theory.

3.1 Review of Orientations, $\mathscr{O}$-, $\mathscr{P}$-, and $\mathscr{F}$-decompositions

The main references for this subsection are [13] and [14].

Definition 3.1**.**

A subset $\mathscr{O}$ of $\mathscr{R}$ is said to be an orientation if for every reaction $y\to y^{\prime}\in\mathscr{R}$, either $y\to y^{\prime}\in\mathscr{O}$ or $y^{\prime}\to y\in\mathscr{O}$, but not both.

For an orientation $\mathscr{O}$, we define a linear map ${L_{\mathscr{O}}}:{\mathbb{R}^{\mathscr{O}}}\to S$ such that

[TABLE]

Let ${{r}_{\text{irr}}}$ and ${{r}_{\text{rev}}}$ be the number of irreversible reactions and reversible reaction pairs in $\mathscr{N}$, respectively. Clearly, $r={{r}_{\text{irr}}}+2{{r}_{\text{rev}}}$. In the succeeding disscussion, we will be using the notation $\mathscr{N}_{\mathscr{O}}$ to denote the subnetwork of $\mathscr{N}$ with respect to the orientation $\mathscr{O}$. The following proposition collects some basic properties of orientations:

Proposition 3.2**.**

Let $\mathscr{N}$ be a CRN and $\mathscr{O}_{\mathscr{N}}$ be the set of orientations of $\mathscr{N}$. Then

• i.

$\left|{\mathscr{O}_{\mathscr{N}}}\right|={2^{r^{\prime\prime}}}$* where $r^{\prime\prime}=r_{\rm{rev}}$, and*

• ii.

the map $P:{\mathscr{O}_{\mathscr{N}}}\to\{\mathscr{N}_{\mathscr{O}}\}$ is bijective. Hence, there are ${2^{r^{\prime\prime}}}$ $\mathscr{O}$-subnetworks in $\mathscr{N}$.

Each orientation $\mathscr{O}$ defines a partition of $\mathscr{N}$ into $\mathscr{O}$ and its complement $\mathscr{O}^{\prime}$, which generates the following decomposition:

Definition 3.3**.**

For an orientation $\mathscr{O}$ on $\mathscr{N}$, the $\mathscr{O}$-decomposition* of $\mathscr{N}$ consists of the subnetworks $\mathscr{N}_{\mathscr{O}}$ and $\mathscr{N}_{\mathscr{O}^{\prime}}$, i.e., $\mathscr{N}=\mathscr{N}_{\mathscr{O}}\cup\mathscr{N}_{\mathscr{O^{\prime}}}$.*

A basic property of an $\mathscr{O}$-decomposition is the following proposition.

Proposition 3.4**.**

Let $\mathscr{N}=\left(\mathscr{S},\mathscr{C},\mathscr{R}\right)$ be a CRN and $\mathscr{O}$ be an orientation. Then the $\mathscr{O}$-decomposition is a trivial decomposition.

Proof.

Let $\mathscr{O}$ be an orientation and $\mathscr{O}^{\prime}=\mathscr{R}\backslash\mathscr{O}$. Then $\mathscr{R}=\mathscr{O}\cup\mathscr{O}^{\prime}$ is a partition with corresponding network decomposition $\mathscr{N}={\mathscr{N}_{\mathscr{O}}}\cup{\mathscr{N}_{\mathscr{O^{\prime}}}}$. Note that $S=S_{\mathscr{O}}+S_{\mathscr{O^{\prime}}}=S_{\mathscr{O}}$. Hence, the stoichiometric subspace of ${\mathscr{N}_{\mathscr{O}}}$ is the same as that of the whole stoichiometric subspace of $\mathscr{N}$. ∎

Corollary 3.5**.**

For any CRN $\mathscr{N}$ and orientation $\mathscr{O}$, $\delta\left(\mathscr{N}\right)=\delta\left({{\mathscr{N}_{\mathscr{O}}}}\right)$.

We now review the important concept of “equivalence classes” from [14]. Let $\left\{{{v^{l}}}\right\}_{l=1}^{d}$ be a basis for $Ker{L_{\mathscr{O}}}$. If for $y\to y^{\prime}\in\mathscr{O}$, $v_{y\to y^{\prime}}^{l}=0$ for all $1\leq l\leq d$ then the reaction $y\to y^{\prime}$ belongs to the zeroth equivalence class $P_{0}$. For $y\to y^{\prime},{\overline{y}\to\overline{y}^{\prime}}\in\mathscr{O}\backslash P_{0}$, if there exists $\alpha\neq 0$ such that $v_{y\to y^{\prime}}^{l}=\alpha v_{\overline{y}\to\overline{y}^{\prime}}^{l}$ for all $1\leq l\leq d$, then the two reactions are in the same equivalence class denoted by $P_{i}$, $i\neq 0$.

Definition 3.6**.**

The $\mathscr{P}$-decomposition* of the $\mathscr{O}$-subnetwork $\mathscr{N}_{\mathscr{O}}$ is the decomposition induced by the partition of $\mathscr{O}$ into equivalence classes.*

Proposition 3.7**.**

Let $\mathscr{O}_{\mathscr{N}}$ be the set of orientations of a CRN.

• i.

The map $P_{d}:\mathscr{O}_{\mathscr{N}}\to\{\mathscr{P}{\text{-decompositions}}\}$ is bijective. For any two orientations $\mathscr{O}$ and $\mathscr{O}^{\prime}$, $\dim{\rm Ker\ }L_{\mathscr{O}}=\dim{\rm Ker\ }L_{\mathscr{O^{\prime}}}={{r}_{\rm irr}}$ + ${{r}_{\rm{rev}}}-s$. Hence, the $\mathscr{P}$-decomposition of $\mathscr{N}_{\mathscr{O}}$ has the same number of subnetworks as that of $\mathscr{N}_{\mathscr{O}^{\prime}}$ (denoted by $w$ if $P_{0}=\emptyset$ and $w+1$ if $P_{0}\neq\emptyset$).

• ii.

For the $i^{\text{t}h}$ subnetworks of the $\mathscr{P}$-decompositions of $\mathscr{N}_{\mathscr{O}}$ and $\mathscr{N}_{\mathscr{O}^{\prime}}$, respectively, the stoichiometric subspaces and the incidence map images coincide, i.e., $S_{i}=S_{i}^{{}^{\prime}}$, and ${\rm Im}\ I_{a,i}={\rm{Im}}\ I_{a,i}^{{}^{\prime}}$.

Proof.

For (i), we note that $\dim{\rm Ker\ }L_{\mathscr{O}}={{r}_{\text{irr}}}$ + ${{r}_{\text{rev}}}-s$ follows immediately from the Rank-Nullity Theorem. For (ii), in the stoichiometric subspaces $S_{i}$ and $S_{i}^{{}^{\prime}}$, the reaction vectors from irreversible reactions are identical, and those from reversible pairs are either identical or the negative of each other. Hence, the stoichiometric subspaces coincide. Similarly, an element of ${\rm Im}\ I_{a,i}$, say $z=\sum{{\alpha_{i}}\left({{\omega_{{y_{i}}}}^{\prime}-{\omega_{{y_{i}}}}}\right)}+\sum{{\beta_{i}}\left({{\omega_{{y_{i}}}}^{\prime}-{\omega_{{y_{i}}}}}\right)}$, where the first sum is over identical reactions (irreversible and identical pair reaction) and the second from the converse pair reaction. Rewriting the latter as $-\sum{{\beta_{i}}\left({{\omega_{{y_{i}}}}-{\omega_{{y_{i}}}}^{\prime}}\right)}$ shows that it belongs to ${\rm{Im}}\ I_{a,i}^{{}^{\prime}}$, and vice versa. ∎

We can now introduce the central concept of “fundamental classes”, which is the basis of the Higher Deficiency Algorithm of Ji and Feinberg. The reactions $y\to y^{\prime}$ and $\overline{y}\to\overline{y}^{\prime}$ in $\mathscr{R}$ belong to the same fundamental class if at least one of the following is satisfied [14].

• i.

$y\to y^{\prime}$ and $\overline{y}\to\overline{y}^{\prime}$ are the same reaction.

• ii.

$y\to y^{\prime}$ and $\overline{y}\to\overline{y}^{\prime}$ are reversible pair.

• iii.

Either $y\to y^{\prime}$ or $y^{\prime}\to y$, and either $\overline{y}\to\overline{y}^{\prime}$ or $\overline{y}^{\prime}\to\overline{y}$ are in the same equivalence class on $\mathscr{O}$.

We can easily see that the orientation $\mathscr{O}$ is partitioned into equivalence classes while the reaction set $\mathscr{R}$ is partitioned into fundamental classes.

Definition 3.8**.**

The $\mathscr{F}$-decomposition* of $\mathscr{N}$ is the decomposition generated by the partition of $\mathscr{R}$ into fundamental classes.*

Proposition 3.9**.**

Any $\mathscr{P}$-decomposition generates the (unique) $\mathscr{F}$-decomposition of $\mathscr{N}$.

Proof.

Note that each partition set (“equivalence class”) of any $\mathscr{P}$-decomposition is contained in a unique partition set (“fundamental class”) of the $\mathscr{F}$-decomposition. Hence, every subnetwork of any $\mathscr{P}$-decomposition is contained in a unique subnetwork of the $\mathscr{F}$-decomposition. Conversely, each $\mathscr{F}$-subnetwork contains a unique $\mathscr{P}$-subnetwork. ∎

An upper bound for the number of fundamental classes is clearly given by $(r_{\text{irr}}+r_{\text{rev}})$, which is the number of reactions in any orientation $\mathscr{O}$. It follows that $w\leq r_{\text{irr}}+r_{\text{rev}}$. This upper bound is sharp in the sense that there are CRNs for which $w=r_{\text{irr}}+r_{\text{rev}}$. The following proposition provides a lower bound for $w$:

Proposition 3.10**.**

For any CRN, $w\geq r_{\text{irr}}+r_{\text{rev}}-s$.

Proof.

$\omega_{r}$ and $\omega_{r^{\prime}}$ are equivalent iff $\omega_{r}=\alpha\omega_{r^{\prime}}$ in ${\mathbb{R}^{\mathscr{O}}/({\text{Ker\ }}L_{\mathscr{O}})^{\perp}}$. The latter means that their cosets are pairwise linearly dependent in the factor space. Hence, the number of equivalence classes is at least $\text{dim }{\mathbb{R}^{\mathscr{O}}/({\text{Ker\ }}L_{\mathscr{O}})^{\perp}}=r_{\text{irr}}+r_{\text{rev}}-s$. Since the fundamental class $C_{0}$ is mapped to 0 in the factor space, it follows that $w\geq r_{\text{irr}}+r_{\text{rev}}-s$. ∎

Corollary 3.11**.**

$w\geq r_{\text{irr}}+r_{\text{rev}}-(n-l)$**

Running Example 3.12**.**

Consider the subnetwork $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}$ of the Schmitz’s carbon cycle model in [18] that shows how the movement of carbon among different pools which represent major parts of the Earth [11]. We label the reactions of the subnetwork, together with kinetic orders, depicted in Figure 1.

The following are the reactions of the subnetwork.

[TABLE]

We choose the orientation $\mathscr{O}=\{R_{1},R_{3},R_{4},R_{5},R_{6},R_{7},R_{8}\}$. Hence, a basis for $Ker{L_{\mathscr{O}}}$ is

[TABLE]

This shows that the $\mathscr{P}$-decomposition induces the following fundamental classes: $\{R_{1},R_{3},R_{4}\}$ and $\{R_{5},R_{6},R_{7},R_{8}\}$. Hence, it follows that the $\mathscr{F}$-decomposition has these fundamental classes: $\{R_{1},R_{2},R_{3},R_{4}\}$ and $\{R_{5},R_{6},R_{7},R_{8}\}$, which are precisely the subnetworks $\mathscr{N}_{1}$ and $\mathscr{N}_{2}$ of $\mathscr{N}$, respectively.

3.2 Independence and Incidence-Independence of the $\mathscr{P}$- and $\mathscr{F}$-decompositions

This subsection shows that the essential properties of any $\mathscr{P}$-decomposition of any $\mathscr{O}$-subnetwork $\mathscr{N}_{\mathscr{O}}$ are fully reflected in the unique $\mathscr{F}$-decomposition of $\mathscr{N}$. We record some basic properties of independent $\mathscr{F}$-decompositions as well as provide an example of a class of CRNs with independent $\mathscr{F}$-decompositions.

Theorem 3.13**.**

Let $\mathscr{N}_{\mathscr{O}}$ be the subnetwork of $\mathscr{N}$ defined by the orientation $\mathscr{O}$ being a subset of $\mathscr{R}$. Then the following holds:

• i.

The $\mathscr{P}$-decomposition of $\mathscr{N}_{\mathscr{O}}$ is independent if and only if the $\mathscr{F}$-decomposition of $\mathscr{N}$ is independent.

• ii.

The $\mathscr{P}$-decomposition of $\mathscr{N}_{\mathscr{O}}$ is incidence-independent if and only if the $\mathscr{F}$-decomposition of $\mathscr{N}$ is incidence-independent.

• iii.

The $\mathscr{P}$-decomposition of $\mathscr{N}_{\mathscr{O}}$ is bi-independent if and only if the $\mathscr{F}$-decomposition of $\mathscr{N}$ is bi-independent.

Proof.

For (i), suppose the $\mathscr{F}$-decomposition is independent. Then, $S$ is the direct sum of the stoichiometric subspaces of the subnetworks corresponding to the fundamental classes. Note, ${P_{i}}\subseteq{C_{i}}$ for all $i$. By the definition of orientation, if $A\to B,B\to A\in\mathscr{R}$ then $A\to B\in P_{i}$ or $B\to A\in P_{i}$, but not both, for some $i$. Thus, $S$ is still the direct sum of the stoichiometric subspaces of the subnetworks corresponding to the equivalence classes. Therefore, the $\mathscr{P}$-decomposition is independent. On the other hand, suppose the $\mathscr{P}$-decomposition is independent. For the remaining reactions in $\mathscr{R}\backslash\mathscr{O}$, by definition, the reversible pairs must belong to the same fundamental class. Hence, $S$ is the direct sum of the stoichiometric subspaces of the subnetworks corresponding to the fundamental classes. Therefore, the $\mathscr{F}$-decomposition is independent. On the other hand, for (ii), suppose the $\mathscr{F}$-decomposition is incidence-independent. By definition, the incidence matrix of the the network is the direct sum of the incidence matrices of the fundamental classes as depicted below where each $F_{i}$ is indexed by complexes (rows) and reactions (columns).

[TABLE]

Note ${P_{i}}\subseteq{C_{i}}$ for all $i$. Also, by definition of orientation which is partitioned by equivalence classes $P_{i}$’s, if a reaction is reversible, one must belong to the orientation and the other must not. Hence, if we remove one of these two reactions corresponding to two columns in the incidence matrix, the dimension of the resulting matrix will not change. This implies the incidence-independence of the $\mathscr{P}$-decomposition. Conversely, suppose the $\mathscr{P}$-decomposition is incidence-independent. Then adding the reversible pair in $P_{i}$ for any $i$ will not change the dimension of the incidence matrix. Thus, $\mathscr{F}$-decomposition is incidence-independent. Statement (iii) follows from (i) and (ii). ∎

Running Example 3.14**.**

We again consider the subnetwork $\mathscr{N}$ of the Schmitz’s carbon cycle model. Since the dimension of the stoichiometric subspaces of the fundamental classes under the $\mathscr{F}$-decomposition is equal to the dimension of the stoichiometric subspaces of $\mathscr{N}$, the $\mathscr{F}$-decomposition is independent. On the other hand, $n-l=6-1=5$ and $\sum{\left({{n_{i}}-{l_{i}}}\right)}=\left({3-1}\right)+\left({4-1}\right)=5$, which proves the incidence-independence of the $\mathscr{F}$-decomposition. Therefore, the said decomposition is bi-independent.

3.3 Independent $\mathscr{F}$-decompositions

In this section, we present a useful necessary condition for an independent $\mathscr{F}$-decomposition and two CRN classes whose $\mathscr{F}$-decompositions are always independent. In Section 4, independent $\mathscr{F}$-decompositions is further analyzed by using the classification introduced by H. Ji into subnetwork types.

3.3.1 A necessary condition for independent $\mathscr{F}$-decompositions

We begin with a general property of independent decompositions.

Proposition 3.15**.**

Let $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup...\cup\mathscr{N}_{k}$ be a CRN decomposition. If the decomposition is independent, then $k\leq s$. Consequently, $k\leq n-l$.

Proof.

Since $s_{i}\geq 1$, we have $\sum s_{i}\geq k$. If $k>s$, then $\sum s_{i}>s$, i.e., the decomposition is dependent, which shows the first claim. The second follows from $0\leq\delta$ or $s\leq n-l$. ∎

Corollary 3.16**.**

For an independent $\mathscr{F}$-decomposition, $w\leq s$. Consequently, $w\leq n-l$.

Proof.

For the $\mathscr{F}$-decomposition, $w=k-1$ (if $P_{0}$ is non-empty) or $w=k$ (otherwise), and the claims follow. ∎

Example 3.17**.**

In [13], we presented the CRNs of a popular model of anaerobic yeast fermentation (Section 3.1) and a model of terrestrial carbon recovery (Section 3.2). For the first CRN, we have $s=7$ and $w=11$ , and for the second one, $s=4$ and $w=6$ . It follows that the $\mathscr{F}$-decompositions of both CRNs are not independent.

3.3.2 S-system CRNs have independent $\mathscr{F}$-decompositions

First, we show that for S-system CRNs, the $\mathscr{F}$-decomposition is a familiar construct. We recall a definition and a result from [6]:

Definition 3.18**.**

Let $R_{j}$ and $P_{j}$ be the set of variables regulating the inflow and outflow reactions of the species $X_{j}$ (i.e., dependent variable) of an S-system, respectively. The species is called reversible if $R_{j}=P_{j}$. Otherwise, it is called irreversible. An S-system is called reversible (irreversible) if all its species are reversible (irreversible).

Our claim is simply the following proposition.

Proposition 3.19**.**

For any S-system embedded CRN, the $\mathscr{F}$-decomposition is the species decomposition.

Proof.

We denote the inflow reaction in $\mathscr{R}_{i}$ with $r_{i}$, the outflow with $r_{-i}$, and the corresponding basis vectors with $\omega_{i}$ and $\omega_{-i}$, respectively. We set $m^{\prime}:=m-m_{\rm rev}$, and index the irreversible species $X_{1},X_{2},...,X_{m^{\prime}}$. Since $|\mathscr{R}|=2m$ and $s=m$, for any orientation, $|O|=2m-m_{\rm{rev}}$ and $\dim{\rm Ker}L_{\mathscr{O}}=m-m_{\text{rev}}$. The $m-m_{\text{rev}}$ vectors $\omega_{i}+\omega_{-i}$, $i=1,2,...,m^{\prime}$ in ${\text{K}er}L_{\mathscr{O}}$ are linearly independent, hence form a basis. On the other hand, the $m$ vectors $\omega_{i}+\omega_{-i},\chi_{j}$ with $i=1,2,...,m^{\prime}$, and $j=1,2,...,m_{\text{rev}}$ and $\chi_{j}$ the reaction from a reversible pair included in the orientation, form a basis for ${\text{K}er}^{\perp}L_{\mathscr{O}}$ . From the $\mathscr{F}$-decomposition definition, the reactions $\omega_{i}$ and $\omega_{-i}$ are equivalent, $i=1,2,...,m$. If $k\neq i$, $\left\langle{\omega_{k}-\alpha\omega_{i},\omega_{i}+\omega_{-i}}\right\rangle=-\alpha$ , so that if $\alpha$ is nonzero, then the $k$-th inflow reaction is not equivalent. Similarly, the $k$-th outflow reaction is not equivalent. Hence, the $\mathscr{F}$-equivalence classes are precisely the $\mathscr{R}_{i}$’s. ∎

Corollary 3.20**.**

The $\mathscr{F}$-decomposition of the embedded CRN of an S-system is independent.

Proof.

It follows from [6] Theorem 1 that the species decomposition is independent, which according to the previous proposition is identical with the $\mathscr{F}$-decomposition. ∎

3.3.3 CRNs of phosphorylation/dephosphorylation systems have independent $\mathscr{F}$-decompositions

In view of their ubiquitous occurrence in cellular signaling networks, phosphorylation/ dephosphorylation (PD) systems have been extensively studied in the CRNT literature. A recent review by Conradi and Shiu [3] lists several classes of multi-site PD processes which have been modeled with mass action systems – in the following proposition, we show that the CRNs of multisite processive and distributive PD processes have distinctively different but both independent $\mathscr{F}$-decompositions.

Proposition 3.21**.**

Let $\mathscr{N}$ be the following CRN for $k$-site processive phosphorylation/ dephosphorylation:

[TABLE]

then the fundamental classes generating the $\mathscr{F}$-decomposition is the full reaction set $\mathscr{R}$. Hence, the $\mathscr{F}$-decomposition consists only of $\mathscr{N}$ and is (trivially) independent.

Proof.

We choose the orientation consisting of the forward reactions and obtain:

[TABLE]

Thus,

[TABLE]

It follows that ${\alpha_{0}}={\alpha_{1}}=...={\alpha_{n}}={\beta_{0}}={\beta_{1}}=...={\beta_{n}}$ and the conclusion follows. ∎

Proposition 3.22**.**

Let $\mathscr{N}$ be the following CRN for $k$-site distributive phosphorylation/ dephosphorylation:

[TABLE]

then

• i.

the fundamental classes generating the $\mathscr{F}$-decomposition are of the form

[TABLE]

• ii.

the $\mathscr{F}$-decomposition is independent.

Proof.

We choose the orientation consisting of the forward reactions and obtain:

[TABLE]

Hence, we have

[TABLE]

Further manipulation gives ${\alpha_{i}}={\beta_{i}}={\lambda_{i}}={\gamma_{i}}$ for each $i=0,...,k-1$, which yields (i). In addition, each of the stoichiometric subspaces of the $k$ subnetworks has dimension 3. Since the dimension of the stoichiometric subspace of the whole network is $3k$, (ii) holds. ∎

Remark 3.23**.**

It is interesting to note that the linkage classes of the distributive CRN are not independent, in contrast to the $\mathscr{F}$-decomposition. The subnetworks of the latter are potentially useful for determining positive equilibria for any kinetics. The review [3] also presents models for dual-site PD with ERK mechanism and mixed processive-distributive mechanism. To complete the picture, we also determine their $\mathscr{F}$-decompositions in the following examples.

Example 3.24**.**

The CRN of dual-site PD with the ERK mechanism is given by:

[TABLE]

in which the $\mathscr{F}$-decomposition consists only of the whole network and hence independent.

Example 3.25**.**

The CRN of dual-site PD with the mixed-mode mechanism is given by:

[TABLE]

in which the $\mathscr{F}$-decomposition consists only of the whole network and hence independent.

3.4 Incidence-independent $\mathscr{F}$-decompositions

This section begins with a necessary condition for an incidence-independent $\mathscr{F}$-decomposition analogous to the proposition in Section 3.3.1 for an independent $\mathscr{F}$-decomposition. We then present a sufficient condition for incidence-independence, namely when the $\mathscr{F}$-decomposition is a $\mathscr{C}$-decomposition, and show that various subsets of S-system CRNs fulfill the condition. Finally, in Section 3.4.3, we show that the CRNs of PD processes are also incidence-independent (and hence bi-independent).

3.4.1 A necessary condition for incidence-independent $\mathscr{F}$-decompositions

The following proposition is the analogue of Proposition 3.15 in Section 3.3.1:

Proposition 3.26**.**

Let $\mathscr{N}=\mathscr{N}_{1}\cup\mathscr{N}_{2}\cup...\cup\mathscr{N}_{k}$ be a CRN decomposition. If the decomposition is incidence-independent, then $k\leq n-l$. If $\mathscr{N}$ has zero deficiency, then $k\leq s$.

Proof.

Since $n_{i}-l_{i}\geq 1$, we have $\sum\left(n_{i}-l_{i}\right)\geq k$. If $k>n-l$, then $\sum\left(n_{i}-l_{i}\right)>n-l$, i.e., the decomposition is incidence-dependent. If the deficiency is zero, $n-l=s$. ∎

Corollary 3.27**.**

For an incidence-independent $\mathscr{F}$-decomposition, $w\leq n-l$. If $\mathscr{N}$ has zero deficiency, then $w\leq s$.

Proof.

For the $\mathscr{F}$-decomposition, $w=k-1$ (if $P_{0}$ is non-empty) or $w=k$ (otherwise), and the claims follow. ∎

Unfortunately, the condition does not seem to be as useful as its analogue as the following example shows.

Example 3.28**.**

The CRN of a popular Generalized Mass Action (GMA) model of anaerobic yeast fermentation (denoted by ERM0-G) was analyzed in [13] for its capacity for multistationarity. The computation of its $\mathscr{F}$-decomposition showed that $w=11<12=n-l$, so that the stated condition is fulfilled. However, the same computation showed that the sum of $n_{i}-l_{i}$ (over the 11 subnetworks) = 13, so that the $\mathscr{F}$-decomposition is not incidence-independent.

3.4.2 A sufficient condition: when the $\mathscr{F}$-decomposition is a $\mathscr{C}$-decomposition

The $\mathscr{F}$-decompositions of S-system CRNs, though always independent, are in general not incidence-independent. To show this, we recall an example from [6].

Example 3.29**.**

The (embedded) S-system CRN of a model of the gene regulatory system of Mycobacterium tuberculosis (Mtb) in the non-replicating phase (NRP) of its life cycle was shown to have $m=40$ species, $n=98$ complexes, $r=80$ irreversible reactions and $l=19$ linkage classes. Since for any (embedded) S-system CRN, the network rank $s=m$, the deficiency $\delta=98-19-40=39<40$. Theorem 1 in [6] implies that the $\mathscr{F}$-decomposition, i.e., species decomposition, is incidence-independent.

However, various subsets of S-system CRNs always have incidence-independent $\mathscr{F}$-decompositions. Before discussing these examples, we show that in many cases, the incidence-independence is due to the fact that the $\mathscr{F}$-decompositions are $\mathscr{C}$-decompositions. We recall the following proposition and proof from [6]:

Proposition 3.30**.**

Any $\mathscr{C}$-decomposition is incidence-independent.

Proof.

Theorem 5 in [6] characterizes $\mathscr{C}$-decompositions as follows: a decomposition is a $\mathscr{C}$-decomposition if and only if it is a coarsening of the linkage class decomposition. Basically, this means that any subnetwork of a $\mathscr{C}$-decomposition is the disjoint union of (some) linkage classes. Since it is also shown (Proposition 3 of [6]) that any coarsening of an incidence-independent decomposition is incidence-independent, it follows that any $\mathscr{C}$-decomposition is incidence-independent. ∎

Corollary 3.31**.**

If the $\mathscr{F}$-decomposition of a CRN is a $\mathscr{C}$-decomposition, then $2w\leq n$.

Proof.

For any $\mathscr{C}$-decomposition, the number of subnetworks $k\leq l$, hence $w\leq l$. Since a $\mathscr{C}$-decomposition is incidence-independent, we have $w+l\leq n$. Combining the two inequalities shows the claim. ∎

Example 3.32**.**

The $\mathscr{F}$-decomposition of any (embedded) S-system CRN with m irreversible species and distinct complexes consists of subnetworks with 2 linkage classes containing the inflow and outflow reaction for each species respectively, and hence is a $\mathscr{C}$-decomposition. Since each CRN has $4m$ complexes, we will denote this subset by Ssys4m.

Example 3.33**.**

For any (embedded) S-system CRN with m reversible species, the $\mathscr{F}$-decomposition coincides with the linkage class decomposition. Each such CRN has $2m$ complexes and the subset will be denoted by Ssys2m.

A further class of CRNs whose $\mathscr{F}$-decompositions are $\mathscr{C}$-decompositions is provided by a special case of Theorem 4.14 in Section 4.

Example 3.34**.**

Let the CRN $\mathscr{N}=\{\mathscr{N}_{i}|\mathscr{N}_{i}=(\mathscr{C}_{i},\mathscr{R}_{i})\}$ with sequence of long monomolecular directed cycles, i.e., of length $\geq$ 3, and $|\mathscr{C}_{i}\cap\mathscr{C}_{j}|=0$ for distinct $i,j=1,2,...,k$.

It is shown in Theorem 4.14, that the $\mathscr{F}$-decomposition consists precisely of the $C_{i}$, which are simultaneously the linkage classes of the network.

We now discuss the CRN of the Heck et al. model of “terrestrial carbon recovery” from [13].

Example 3.35**.**

The CRN is given by:

[TABLE]

The $\mathscr{F}$-decomposition is a curiosity: it is almost a $\mathscr{C}$-decomposition, with 4 of 6 fundamental classes coinciding with the linkage class reaction sets $\mathscr{C}_{1}=\{R_{1}\}$, $\mathscr{C}_{2}=\{R_{2}\}$, $\mathscr{C}_{3}=\{R_{3},R_{4},R_{8},R_{9}\}$, and $\mathscr{C}_{4}=\{R_{6}\}$. For the remaining two, $\mathscr{C}_{5}=\{R_{5},R_{10}\}\neq\{R_{5},R_{7}\}$ and $\mathscr{C}_{6}=\{R_{7}\}\neq\{R_{10}\}$, but fortuitously, the images of their incidence maps are respectively isomorphic. Hence the $\mathscr{F}$-decomposition is also incidence-independent.

As the final example in this section, we present a subset of S-system CRNs whose $\mathscr{F}$-decompositions are always incidence-independent but which are not $\mathscr{C}$-decompositions.

Example 3.36**.**

The set of (embedded) S-system CRNs with self-regulating species and non-regulated outflows have $\mathscr{F}$-subnetworks of the form $\left\{{0\leftarrow{x_{i}}\to 2{x_{i}}}\right\}$, i.e., there are $2m+1$ complexes and one linkage class, inferring $n-l=2m$. Clearly, it is not a $\mathscr{C}$-decomposition. Since each $\mathscr{F}$-subnetwork has 3 complexes and a linkage class, the sum of $n_{i}-l_{i}=2m$, too, showing incidence-independence. We denote this subset as Ssys2m+1.

3.4.3 CRNs of phosphorylation/dephosphorylation systems have incidence-independent $\mathscr{F}$-decompositions

In this section, we use the computations for $\mathscr{F}$-decompositions of the CRNs of PD systems in Section 3.3.3 to show that they are also incidence-independent (hence bi-independent). For multisite processive PD systems, this is trivial since the $\mathscr{F}$-decomposition has only one subnetwork. The cases of dual-site ERK mechanism and mixed-mechanism are discussed in the following example.

Example 3.37**.**

The CRN of dual-site PD with the mixed-mode mechanism is given by:

[TABLE]

in which the $\mathscr{F}$-decomposition consists only of the whole network and hence incidence-independent.

The following proposition completes the picture by providing the proof in the multisite distributive PD CRN to be incidence-independent.

Proposition 3.38**.**

The $k$-site distributive phosphorylation/dephosphorylation:

[TABLE]

has incidence-independent $\mathscr{F}$-decomposition.

Proof.

In Proposition 3.22, it was shown that the fundamental classes generating the $\mathscr{F}$-decomposition are of the form

[TABLE]

Since the network has $2(2k)=4k$ reactions, and the reactants and products are distinct, it has $4k+2$ complexes. Now,

[TABLE]

∎

Clearly, the image of $I_{a,\mathscr{O}}$, the restriction of the incidence map to $R^{\mathscr{O}}$, is equal to the image of $I_{a}$, since the image of a reverse reaction is just the negative of the forward reaction and the linkage classes are not changed. Hence, $r_{\text{irr}}+r_{\text{rev}}-(n-l)\geq 0$, being the $\text{dim }\text{Ker}I_{a,\mathscr{O}}$, or $n-l\leq r_{\text{irr}}+r_{\text{rev}}$. As shown in propositions above, for independent and incidence-independent $\mathscr{F}$-decompositions, $w\leq n-l$. To date, we have not yet found a CRN with a dependent and incidence-dependent $\mathscr{F}$-decomposition with $w>n-l$, so one can ask the question whether $n-l$ is an upper bound for $w$ in general.

4 Types of $\mathscr{F}$-decomposition and Network Properties

Ji classified the subnetworks occurring in a $\mathscr{P}$-decomposition into 3 types and summarized their properties as follows (Proposition 2.5.4 in [14]).

Lemma 4.1**.**

[14]* Let $\mathscr{N}=\left(\mathscr{S},\mathscr{C},\mathscr{R}\right)$ be a CRN and $\mathscr{O}$ be an orientation. Let $\mathscr{{N}}_{\mathscr{O},i}$ for $i=0,1,2,...,w$ be defined as the subnetwork generated by all reactions in $P_{i}$. Then one of the following holds:*

• i.

The reaction vectors for $\mathscr{{N}}_{\mathscr{O},i}$ are linearly independent, and the subnetwork $\mathscr{{N}}_{\mathscr{O},i}$ based on $P_{i}$ forms a forest (i.e., a graph with no cycle) with deficiency 0.

• ii.

The reaction vectors are minimally dependent, and the subnetwork $\mathscr{{N}}_{\mathscr{O},i}$ based on $P_{i}$ forms a forest with deficiency 1.

• iii.

The reaction vectors are minimally dependent, and the subnetwork $\mathscr{{N}}_{\mathscr{O},i}$ based on $P_{i}$ forms a big cycle (with at least three vertices) with deficiency 0.

We will denote the subnetwork classes in i, ii, and iii of Lemma 4.1 as Type I, Type II and Type III subnetworks respectively. Since our focus is on the $\mathscr{F}$-decomposition, we extend this classification as follows: an $\mathscr{F}$-subnetwork is of type I, II or III if it contains a $\mathscr{P}$-subnetwork of type I, II or III, respectively. Note that while the characterization of Type I and II $\mathscr{P}$-subnetworks as forests is lost, that Type III subnetwork as a big cycle is retained in the Type III $\mathscr{F}$-subnetwork. More importantly, the deficiency of each subnetwork type remains the same. We assign the numbers of fundamental classes for Types I, II and III with the symbols $w_{I}$, $w_{II}$ and $w_{III}$, respectively.

Definition 4.2**.**

An $\mathscr{F}$-decomposition is said to be

• i.

Type I* if it contains Type I subnetwork only.*

• ii.

Type II* if it contains Type II subnetwork only.*

• iii.

Type III* if it contains Type III subnetwork only.*

The first important consequence of the above classification is the following:

Proposition 4.3**.**

If a CRN has an independent $\mathscr{F}$-decomposition, then its deficiency $\delta\leq w_{II}$.

Proof.

Note that the $\mathscr{F}$-decomposition is independent, so $\delta\leq\delta_{1}+\delta_{2}...+\delta_{w}$. Since any Type I or Type III subnetworks have zero deficiency, we obtain the claim. ∎

Corollary 4.4**.**

A CRN whose $\mathscr{F}$-decomposition has no Type II subnetworks and is independent has zero deficiency. In particular, this is true for CRNs with Type I and Type III independent $\mathscr{F}$-decompositions.

The special case of the corollary above is significant because the properties of the positive equilibria sets of deficiency zero networks are well-known. In particular:

• i.

any positive equilibrium is complex-balanced [9],

• ii.

no positive equilibrium exists if the network is not weakly reversible (Deficiency Zero Theorem), and

• iii.

for some subsets of the power-law kinetics, existence and parametrization results on equilibria are proven, i.e., Deficiency Zero Theorem for MAK (Feinberg-Horn-Jackson), PL-RDK with zero kinetic deficiency (M${\rm\ddot{u}}$ller-Regensburger [16]), PL-TIK (Talabis et al. [20]), and PL-NDK with special independent decompositions (Fortun et al. [11]).

4.1 Independent Type I $\mathscr{F}$-decompositions

Proposition 4.5**.**

Let $s$ be the rank of a network $\mathscr{N}$. If for an orientation $\mathscr{O}$, $s=|\mathscr{O}|$, then $\delta\left(\mathscr{N}\right)=0$.

Proof.

Let $s$ be the rank of a network $\mathscr{N}$ with an orientation $\mathscr{O}$. Note that $s=\dim S$. Each reaction in $\mathscr{O}$ has a corresponding reaction vector. With the assumption that $s=|\mathscr{O}|$, the reaction vectors are linearly independent. From Lemma 4.1, there are Type I subnetworks only in the $\mathscr{P}$-decomposition. By definition, the $\mathscr{P}$-decomposition is of Type I. From Corollary 4.4, where the number of Type II subnetworks in the $\mathscr{P}$-decomposition is zero, it follows that $\delta\left(\mathscr{N}\right)=0$. ∎

Example 4.6**.**

If the network has only irreversible reactions, then the only orientation is the whole set of reactions. In this case, a network with independent Type I $\mathscr{P}$-decomposition is a trivial nullspace network.

Example 4.7**.**

The CRNs of the set Ssys2m introduced in Section 3.4.2 all have Type I $\mathscr{F}$-decompositions.

Proposition 4.8**.**

Let $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ be a PL-RDK system such that there is at least one irreversible reaction in $\mathscr{R}$. Let $\mathscr{{N}}_{\mathscr{O},i}$ for $i=0,1,2,...,w$ be defined as the subnetwork generated by all reactions in $P_{i}$. If the $\mathscr{P}$-decomposition is independent, and the reaction vectors in $\mathscr{{N}}_{\mathscr{O},i}$ are linearly independent for each $i$, then the system does not have the capacity to admit multiple equilibria.

Proof.

Suppose the $\mathscr{P}$-decomposition is independent and the reaction vectors in $\mathscr{{N}}_{\mathscr{O},i}$ are linearly independent for each $i$. Thus, $Ker{L_{\mathscr{O}}}$ is trivial (containing the zero vector only). So every reaction must be placed in the zeroth equivalence class $P_{0}$. But there is one irreversible reaction contradicting the rule in the higher deficiency algorithm [13, 14] that each reaction in $P_{0}$ must be reversible (with respect to $\mathscr{R}$). Therefore, the system does not have the capacity to admit multiple equilibria. ∎

4.2 Independent Type II $\mathscr{F}$-decompositions

The following proposition expresses an apparently rare relationship between deficiency and rank of a CRN.

Proposition 4.9**.**

For a CRN with an independent Type II $\mathscr{F}$-decomposition, $\delta\leq s$ or equivalently, $s\leq n-l\leq 2s$.

Proof.

We have for such a CRN, $\delta\leq w_{II}=w\leq s$. ∎

Example 4.10**.**

The embedded CRN of an S-system with only irreversible species is an example of an independent Type II $\mathscr{F}$-decomposition network. This shows that the formula in [6] can be seen in this case as an instance of Proposition 4.9. If the $\mathscr{F}$-decomposition is also incidence-independent, then the network deficiency is the number of species $m$. The sets Ssys4m and Ssys2m+1 discussed in Section 3.4.2 are subsets of this set of S-system CRNs.

Example 4.11**.**

The following CRN is the well-known model of the EnvZ-OmpR system of E. coli studied by Shinar and Feinberg in [19]:

[TABLE]

The $\mathscr{F}$-decomposition, like that of the multisite processive PD model, has only one subnetwork (i.e., the whole CRN) and is of Type II (its $\mathscr{F}$-subnetwork is clearly a forest of deficiency 1).

Example 4.12**.**

In [21] and [22], it is shown that evolutionary games with replicator dynamics can be represented as chemical kinetic systems. The following example is a symmetric population game with 2 pure strategies and non-linear, continuous payoff functions, also called “playing the field” games [4]. Let $x=(x_{1},x_{2})$ be the vector of pure strategies and $F$ a $2\times 2$ matrix of nonnegative real numbers with $F_{i}$ as its $i$-th row. The payoff function $f_{i}$ is defined as $f_{i}(x)=x^{F_{i}}$. A PLK representation is then given by the CRN:

[TABLE]

The rate functions for the reactions are given by $K_{i}(x)=x_{i}f_{i}(x)$, $K_{-i}^{\prime}(x)=x_{i}x_{1}f_{1}(x)$ and $K_{-i}(x)=x_{i}x_{2}f_{2}(x)$. The rate constants are set to 1 to ensure dynamic equivalence with the replicator equation.

To determine the $\mathscr{F}$-decomposition, we consider the orientation given by $\{R_{i},R_{-i}\}$. The subnetwork generated by this orientation is an S-system in 2 irreversible variables, and its $\mathscr{F}$-decomposition is independent. Since a $\mathscr{P}$-decomposition of the game’s CRN is independent, then its $\mathscr{F}$-decomposition is independent too.

4.3 Independent Type III $\mathscr{F}$-decompositions

Running Example 4.13**.**

The subnetwork of the Schmitz’s carbon cycle model is clearly an instance of an independent Type III $\mathscr{F}$-decomposition. We formulate a generalization in the following theorem:

Theorem 4.14**.**

The following family of CRNs has bi-independent Type III $\mathscr{F}$-decomposition such that the $\mathscr{N}_{i}$’s are precisely the fundamental classes under the decomposition: $\mathscr{N}=\{\mathscr{N}_{i}|\mathscr{N}_{i}=(\mathscr{C}_{i},\mathscr{R}_{i})\}$ with a (possibly broken) chain of long monomolecular directed cycles, i.e. of length $\geq$ 3, and $|\mathscr{C}_{i}\cap\mathscr{C}_{j}|\leq 1$ if $j=i+1$ for $i=0,1,...,k-1$.

Proof.

To simplify the proof, we will only show the case with no break, since the one with break in the graph is rather obvious. Without loss of generality, we assume the orientation given by the graph in Figure 2. Note that there are exactly $\ell=\ell_{1}+\ell_{2}+...+\ell_{k}-(k-1)$ complexes in the network. In solving a basis for ${\rm Ker\ }L_{\mathscr{O}}$, we have the following equation:

[TABLE]

[TABLE]

[TABLE]

Hence, we have:

[TABLE]

[TABLE]

[TABLE]

For the first subnetwork, ${\alpha_{{\ell_{1}}}}={\alpha_{1}}={\alpha_{2}}=...={\alpha_{{q_{1}}-1}}$ and ${\alpha_{{q_{1}}}}={\alpha_{{q_{1}}+1}}=...={\alpha_{{\ell_{1}}}}$, which gives ${\alpha_{1}}={\alpha_{2}}=...={\alpha_{{\ell_{1}}}}$. A similar proof can be provided for the last subnetwork. For the remaining subnetworks, with $i=2,3,...,k-1$, we obtain the summand:

[TABLE]

which yields

[TABLE]

Note that the “apostrophe” symbol is used to differentiate the positions of the complexes from two consecutive subnetworks. Now, the term $\left({{C_{{q_{i}}-1}}}\right)^{\prime}\left({\left({{\alpha_{{q_{i}}-2}}}\right)^{\prime}-\left({{\alpha_{{q_{i}}-1}}}\right)^{\prime}}\right)$ yields $\left({{\alpha_{{q_{i}}-2}}}\right)^{\prime}=\left({{\alpha_{{q_{i}}-1}}}\right)^{\prime}$ which implies the equality of the $\alpha_{r}$’s with position $r<q_{i}$. Similarly, the term $\left({{C_{{q_{i}}+1}}}\right)^{\prime}\left({\left({{\alpha_{{q_{i}}}}}\right)^{\prime}-\left({{\alpha_{{q_{i}}+1}}}\right)^{\prime}}\right)$ yields $\left({{\alpha_{{q_{i}}}}}\right)^{\prime}-\left({{\alpha_{{q_{i}}+1}}}\right)^{\prime}$ which implies the equality of $\alpha_{r}$’s with position $r>q_{i}$. We also obtain

[TABLE]

Since the complex $C_{q_{i}}$ is also present in the $i+1$-st subnetwork, we get

[TABLE]

It follows that $\left({{\alpha_{{q_{i}}+1}}}\right)^{\prime}={\alpha_{{q_{i}}}}$. But ${\alpha_{{\ell_{i}}+1}}={\alpha_{{\ell_{i+1}}}}$, which proves the equality of the $\alpha$’s in a subnetwork. Thus, the subnetworks are precisely the fundamental classes which are independent. Indeed, the following is a basis for ${\rm Ker\ }L_{\mathscr{O}}$:

[TABLE]

where each $F_{i}$ is an $\ell_{i}\times 1$ matrix with entries all equal to 1.

Since $\mathscr{N}$ has zero deficiency, this also proves the incidence-independence of the $\mathscr{F}$-decomposition. ∎

5 The CF-RI+ Transformation Method

In this section, we present a transformation method whose key property is that it maps an irreversible reaction (a reversible pair of reactions) of the original system to an irreversible reaction (a reversible pair of reactions) of the target system. In other words, it is reversibility and irreversibility (RI) preserving. This method was based on the generic CF-RM method (transformation of complex factorizable kinetics by reactant multiples) which converts a PL-NDK to a PL-RDK system. We add in the notation CF-RI a sub-index “+” for two reasons: to indicate the “positive” (or preserving) relation and to highlight its partial coincidence with the CF-RM+ variant of CF-RM. However, in most cases, CF-RI+ adds new reactants which are not reactant multiples, so it is not a CF-RM variant.

5.1 Review of the CF-RM+ Method

We present the CF-RM transformation method in [17]. One can construct a PL-RDK system from a given PL-NDK system using this method. A CF-subset contains reactions having the same kinetic order vectors. At each reactant complex, the branching reactions are partitioned into CF-subsets. An NF-reactant complex has more than one CF-subset which makes the system NDK. For each subset, a complex is added to both the reactant and the product complexes of a reaction leaving the reaction vectors unchanged. The kinetic order matrix does not change as well.

The CF-RM method is given by the following steps.

Determine the set of reactant complexes $\rho\left(\mathscr{R}\right)$.

• 2.

Leave each CF-reactant complex unchanged.

• 3.

At an NF-reactant complex, select a CF-subset containing the highest number of reactions and leave this CF-subset unchanged. For each of the remaining $N_{R}(y)-1$ CF-subsets, choose successively a multiple of $y$ which is not among the current set of reactants. Different procedures are possible for the selection of a new reactant as long as it is different from those in the current reactant set. After each choice, the current set is updated.

CF-RM+ is a variant of CF-RM. All the steps are identical with the generic CF-RM method except that it uses additional criteria in the selection of the new reactant multiples. CF-RM+ chooses the reactant multiple so that the new reactant differs from all existing complexes and all the new product complexes in the CF-subset also differ from all existing complexes [17].

5.2 Details of the CF-RI+ Method

Note that the CF-RM+ method given in [17] updates the set of current complexes and complexes in the transform after each CF-subset of an NF-node is processed. The CF-RI+ method proceeds as follows:

Determine the reactant set $\rho({\mathscr{R})}$ and identify the subset $\rho({\mathscr{R})}_{CF}$ of CF-nodes.

• 2.

If the reaction set ${\mathscr{R}_{y}:=\rho^{-1}(y)}$ of a CF-node $y$ has no reversible reaction with an NF-node, then it is left unchanged.

• 3.

At an NF-node without reversible reactions, carry out the steps of CF-RM+.

• 4.

At an NF-node with a reversible reaction, among the CF-subsets without a reversible reaction (if there are any), select one with the highest number of reactions and leave this unchanged.

• 5.

For the remaining CF-subsets without a reversible reaction, carry out CF-RM+.

• 6.

For a CF-subset with a reversible reaction, carry out CF-RM+, but in addition, for each reversible reaction, also for the CF-subset of the reverse reaction (with the same “catalytic” complex). If the reactant complex of the reverse reaction is an NF-node, this additional step removes the original CF-subset from the reaction set of that NF-node. If this removal transforms the NF-node to a CF-node, then remove the node from the list of NF-nodes (to be processed).

Remark 5.1**.**

It is in the last step that the resulting new reactant may be a non-multiple of the original reactant, since the “catalytic” complex added is determined by the reactant of the other reaction in the reversible pair.

Two basic properties of the CF-RI+ transformation are collected in the following proposition:

Proposition 5.2**.**

Let $\mathscr{N}_{RI}$ be the CF-RI+ transform of $\mathscr{N}$.

• i.

If the CRN has no reversible reactions, then CF-RI+* = CF-RM*+.

• ii.

The stochiometric subspaces are equal, i.e., $S_{RI}=S$.

6 The ${\mathscr{F}}$-decomposition under the CF-RI+ Transformation

The following theorem implies that with or without the application of the CF-RM transformation, the computation on determining whether a PL-NDK system has the capacity to admit multiple equilibria using the Multistationarity Algorithm are the same with the assumption of the independence of the $\mathscr{F}$-decomposition.

Theorem 6.1**.**

Let $({\mathscr{N}},K)$ be a PL-NDK system and $({\mathscr{N}}_{RI},K_{RI})$ a CF-RI+ transform. Then

• i.

for any orientation ${\mathscr{O}}$ of ${\mathscr{N}}$, $|{\mathscr{O}}|=|{\mathscr{O}_{RI}}|$, and

• ii.

the ${\mathscr{F}}$-decomposition of ${\mathscr{N}}$ is independent if and only if the ${\mathscr{F}}$-decomposition of ${\mathscr{N}}_{RI}$ is independent.

Proof.

Since the transformation preserves the reversibility and irreversibility of the reactions, $|{\mathscr{O}}|=|{\mathscr{O}_{RI}}|$. Suppose the $\mathscr{F}$-decomposition of ${\mathscr{N}}$ is independent. By Theorem 3.13, $\mathscr{P}$-decomposition of ${\mathscr{N}}$ is independent. Note that the reaction vectors remain the same after the application of the transformation. Moreover, we assume that the reversibility and irreversibility of the reactions are retained and we have

[TABLE]

Hence, we can choose the same basis for $KerL_{\mathscr{O}}$ and $KerL_{\mathscr{O}_{RI}}$ such that the order of the rows of the reactions corresponding to the basis remains the same. Thus, the equivalence classes are retained under the transformation. Therefore, the $\mathscr{P}$-decomposition of ${\mathscr{N}}_{RI}$ is independent. It follows that the $\mathscr{F}$-decomposition of ${\mathscr{N}}_{RI}$ is independent. The same proof applies for the converse. ∎

Running Example 6.2**.**

We apply the CF-RI+ transform to the reaction network of the Schmitz’s carbon cycle model. We modify $R_{5}$ and obtain the following dynamically equivalent PL-RDK which also has an independent $\mathscr{F}$-decomposition.

[TABLE]

7 Conclusion and Outlook

We summarize our results and provide some direction for future research.

We introduced the ${\mathscr{O}}$-, ${\mathscr{P}}$-, and ${\mathscr{F}}$-decompositions underlying the HDA which is the basis of the multistationarity algorithm MSA for power-law kinetics. We derive properties of these decompositions such as independence and incidence-independence, and identify network classes where the $\mathscr{F}$-decomposition coincides with other known decompositions.

• 2.

We classified the $\mathscr{F}$-decomposition into three types according to the types of subnetworks induced by the decomposition. We explored the network properties of each of these types of decompositions.

• 3.

As our major examples, we have determined that the CRNs of phosphorylation/ dephosphorylation systems have bi-independent $\mathscr{F}$-decompositions. We also used a subnetwork of the Schmitz’s carbon cycle model as a running example. We have shown that the ${\mathscr{F}}$-decomposition is bi-independent. We generalized this type of subnetworks with a chain of long monomolecular directed cycles which can be possibly broken.

• 4.

We have shown that for independent $\mathscr{F}$-decomposition, the additional CF-RM transformation is not needed, and hence the MSA can be applied directly to the system.

• 5.

One can prove the results for a larger class containing the set of independent $\mathscr{F}$-decompositions.

Acknowledgement

BSH acknowledges the support of DOST-SEI (Department of Science and Technology-Science Education Institute), Philippines for the ASTHRDP Scholarship grant.

Appendix A Nomenclature

A.1 List of abbreviations

[TABLE]

A.2 List of important symbols

[TABLE]