A Computational Approach to Multistationarity of Power-Law Kinetic Systems

TL;DR

This paper introduces a computational algorithm to determine multistationarity in power-law kinetic systems, extending existing methods and applying them to biological and ecological models, revealing conditions for multiple steady states.

Contribution

The paper develops the Multistationarity Algorithm for PLK systems, extending deficiency-based methods to power-law kinetics and demonstrating its effectiveness on biological and ecological models.

Findings

Proved monostationarity of a yeast fermentation model

Established multistationarity in a global carbon cycle model

Extended deficiency one algorithm to PLK systems

Abstract

This paper presents a computational solution to determine if a chemical reaction network endowed with power-law kinetics (PLK system) has the capacity for multistationarity, i.e., whether there exist positive rate constants such that the corresponding differential equations admit multiple positive steady states within a stoichiometric class. The approach, which is called the "Multistationarity Algorithm for PLK systems" (MSA), combines (i) the extension of the "higher deficiency algorithm" of Ji and Feinberg for mass action to PLK systems with reactant-determined interactions, and (ii) a method that transforms any PLK system to a dynamically equivalent one with reactant-determined interactions. Using this algorithm, we obtain two new results: the monostationarity of a popular model of anaerobic yeast fermentation pathway, and the multistationarity of a global carbon cycle model with climate engineering, both in the generalized mass action format of biochemical systems theory. We also provide examples of the broader scope of our approach for deficiency one PLK systems in comparison to the extension of Feinberg's "deficiency one algorithm" to such systems.

Tables

Table 1. Table 1 : Nonterminal and terminal strong linkage classes of ERM0-G

Equivalence	Fundamental	Nonterminal Strong	Terminal Strong
Class	Class	Linkage Class	Linkage Class
$P_1=\{R_1,R_2\}$	$C_1=\{R_1,R_2\}$	$\left\{X_2\right\}$	$\left\{X_1+X_2\right\}$
		$\left\{X_1+X_5\right\}$	$\left\{X_2+X_5\right\}$
$P_2=\left\{R_3\right\}$	$C_2=\left\{R_3\right\}$	$\left\{2X_5+X_1\right\}$	$\left\{X_5+X_1\right\}$
$P_3=\left\{R_4\right\}$	$C_3=\left\{R_4\right\}$	$\left\{X_2+X_5\right\}$	$\left\{X_3+X_5\right\}$
$P_4=\left\{R_5\right\}$	$C_4=\left\{R_5\right\}$	$\left\{2X_5+X_2\right\}$	$\left\{X_5+X_2\right\}$
$P_5=\left\{R_6\right\}$	$C_5=\left\{R_6\right\}$	$\left\{X_2+X_5\right\}$	$\left\{X_5\right\}$
$P_6=\left\{R_7\right\}$	$C_6=\left\{R_7\right\}$	$\left\{X_2+X_5\right\}$	$\left\{X_2\right\}$
$P_7=\{R_8,R_{11}\}$	$C_7=\{R_8,R_{11}\}$	$\left\{X_3+X_5\right\}$	$\left\{X_4+X_5\right\}$
		$\left\{X_3+X_4+X_5\right\}$	$\left\{X_3+X_5\right\}$
$P_8=\left\{R_9\right\}$	$C_8=\left\{R_9\right\}$	$\left\{X_3+X_5\right\}$	$\left\{X_3+2X_5\right\}$
$P_9=\left\{R_{10}\right\}$	$C_9=\left\{R_{10}\right\}$	$\left\{X_3+X_4+X_5\right\}$	$\left\{X_4+X_5\right\}$
$P_{10}=\left\{R_{12}\right\}$	$C_{10}=\left\{R_{12}\right\}$	$\left\{X_3+X_4+X_5\right\}$	$\left\{X_3+X_4+2X_5\right\}$
$P_{11}=\left\{R_{13}\right\}$	$C_{11}=\left\{R_{13}\right\}$	$\left\{2X_5\right\}$	$\left\{X_5\right\}$

Table 2. Table 2: Network Numbers

	ERM0-G	Global Carbon
		Cycle Model
species	5	5
complexes	13	14
reactant complexes	8	9
reactions	13	10
irreversible reactions	13	6
linkage classes	1	6
strong linkage classes	13	12
terminal sl classes	5	6
rank of network	5	4
deficiency	7	4

Table 3. Table 3: Equilibria in Heck et al.’s Global Carbon Cycle Model

	$\mu$	$\sigma \in S$	$c$**	$c$*
$A_1$	$-0.062895375$	$-1$	16.40466123	15.40466123
$A_2$	$-0.262273058$	$-1$	4.334651228	3.334651228
$A_3$	$0.1$	1	9.508331945	10.50833194
$A_4$	0.291558477	$1$	2.954105867	3.954105867
$A_5$	$1.194680148$	1	0.434310288	1.434310288

Table 4. Table 4 : Equilibria in Example 4.1

	$\mu$	$\sigma \in S$	$c$**	$c$*
$A_1$	1	3	1.745930121	4.745930121
$A_2$	$0$	$0$	1	1
$A_3$	$1$	2	1.163953414	3.163953414

Table 5. Table 5 : Summary of values for Example 4.1

$y\to y^{\prime }$	${\kappa }_{y\to y^{\prime }}$	$c^{\ast \ast T_{.y}}$	$c^{\ast T_{.y}}$	$k_{y\to y^{\prime }}=\frac{{\kappa }_{y\to y^{\prime }}}{c^{\ast \ast T_{.y}}}$	${\kappa }_{y\to y^{\prime }}{e}^{T_{.y}\cdot \mu }$
$R_1$	1	1.425545974	3.875035717	0.701485619	2.718281828…$=1e$
$R_2$	2	1.745930121	4.745930121	1.145521219	5.436563657…$=2e$
$R_3$	2	1.163953414	3.163953414	1.718281828	5.436563657…$=2e$
$R_4$	1	1.163953414	3.163953414	0.859140914	2.718281828…$=1e$
$R_5$	1	1.163953414	3.163953414	0.859140914	2.718281828…$=1e$

Table 6. Table 6: Summary of values for for the Heck et al.’s Terrestrial Carbon Recovery Model

$y\to y^{\prime }$	${\kappa }_{y\to y^{\prime }}$	$k_{y\to y^{\prime }}$	${\kappa }_{y\to y^{\prime }}{e}^{T_{.y}\cdot \mu }$
$R_1$	$1.09\times {10}^{-5}$	$1.4014\times {10}^{-193}$	0.241022556
$R_2$	$1.09\times {10}^{-5}$	$1.4814\times {10}^{-159}$	0.00768157
$R_3$	0.694736046	0.160274959	0.534461087
$R_4$	0.694736046	0.073066028	0.767802074
$R_5$	1	1.222817833	8.426368578
$R_6$	1.58138827	$2.1781\times {10}^{43}$	12814.12188
$R_7$	0.58138827	$3.32525\times {10}^{57}$	12805.92885
$R_8$	1	0.338511903	1.338511903
$R_9$	1	0.105170918	1.105170918
$R_{10}$	1	$3.44864\times {10}^{-10}$	8.426368578

Equations

$$\sum _{y\to y^{\prime }\in \mathcal{R}}{\kappa }_{y\to y^{\prime }}\left(y^{\prime }-y\right)=0\mathrm{and}\sum _{y\to y^{\prime }\in \mathcal{R}}{\kappa }_{y\to y^{\prime }}{e}^{T_{.y}\cdot \mu }\left(y^{\prime }-y\right)=0$$

$$e^{T_{.y}\cdot \mu }=\prod _{s\in \mathcal{S}}{e}^{{\mu }_s{\left(T{.}_y\right)}_s}=\prod _{s\in \mathcal{S}}{\left(\frac{{c_s}^{\ast }}{{c_s}^{\ast \ast }}\right)}^{{\left(T{.}_y\right)}_s}=\frac{c{{}^{\ast }}^{T_{.y}}}{c{{}^{\ast \ast }}^{T_{.y}}}.$$

$$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}:X_{2}\to X_{1}+X_{2}&&R_{8}:X_{3}+X_{5}\to X_{4}+X_{5}\\ R_{2}:X_{1}+X_{5}\to X_{2}+X_{5}&&R_{9}:X_{3}+X_{5}\to X_{3}+2X_{5}\\ R_{3}:2X_{5}+X_{1}\to X_{5}+X_{1}&&R_{10}:X_{3}+X_{4}+X_{5}\to X_{4}+X_{5}\\ R_{4}:X_{2}+X_{5}\to X_{3}+X_{5}&&R_{11}:X_{3}+X_{4}+X_{5}\to X_{3}+X_{5}\\ R_{5}:2X_{5}+X_{2}\to X_{5}+X_{2}&&R_{12}:X_{3}+X_{4}+X_{5}\to X_{3}+X_{4}+2X_{5}\\ R_{6}:X_{2}+X_{5}\to X_{5}&&R_{13}:2X_{5}\to X_{5}\\ R_{7}:X_{2}+X_{5}\to X_{2}&\\ \end{array}$$

$$\bordermatrix{&v_{1}&v_{2}&v_{3}&v_{4}&v_{5}&v_{6}&v_{7}&v_{8}\cr R_{1}&0&1&0&0&1&1&0&0\cr R_{2}&0&1&0&0&1&1&0&0\cr R_{3}&-1&0&-1&1&0&0&1&-1\cr R_{4}&0&0&0&0&1&1&0&0\cr R_{5}&1&0&0&0&0&0&0&0\cr R_{6}&0&1&0&0&0&0&0&0\cr R_{7}&0&0&1&0&0&0&0&0\cr R_{8}&0&0&0&0&0&1&0&0\cr R_{9}&0&0&0&1&0&0&0&0\cr R_{10}&0&0&0&0&1&0&0&0\cr R_{11}&0&0&0&0&0&1&0&0\cr R_{12}&0&0&0&0&0&0&1&0\cr R_{13}&0&0&0&0&0&0&0&1\cr}$$

$${\Gamma }_W=\left\{x\in {\mathbb{R}}^{\mathcal{O}}|x\mathrm{has}\mathrm{support}\mathrm{in}W\right\}$$

$$D=\left\{y_i\to {y_i}^{\prime }\in W|g_W\left(y_i\to {y_i}^{\prime }\right)=0\right\}$$

$$ND=\left\{y_i\to {y_i}^{\prime }\in W|g_W\left(y_i\to {y_i}^{\prime }\right)\ne 0\right\}.$$

$$R_{+}^j=$$

$$R_{-}^j=$$

$$Q_{+}^j=$$

$$Q_{-}^j=$$

$$Q_1^j=$$

$$Q_2^j=$$

$$\sum _{y_i\to {y_i}^{\prime }\in ND}{b}_{y_i\to {y_i}^{\prime }}^j{g}_W\left(y_i\to {y_i}^{\prime }\right)=0,j=1,2,...,q$$

$$\sum _{y_i\to {y_i}^{\prime }\in ND}{\rho }_W\left(y_i\to {y_i}^{\prime }\right)b_{y_i\to {y_i}^{\prime }}^j{g}_W\left(y_i\to {y_i}^{\prime }\right)$$

$$+\sum _{y_i\to {y_i}^{\prime }\in D}{b}_{y_i\to {y_i}^{\prime }}^j{h}_W\left(y_i\to {y_i}^{\prime }\right)=0,j=1,2,...,q$$

$$h_W\left(y_i\to {y_i}^{\prime }\right)={\rho }_W\left(y_i\to {y_i}^{\prime }\right)g_W\left(y_i\to {y_i}^{\prime }\right),y_i\to {y_i}^{\prime }\in ND.$$

$$\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}$$

$$\bordermatrix{&v^{1}&v^{2}&v^{3}&v^{4}\cr R_{1}&1&0&1&0\cr R_{2}&1&0&0&0\cr R_{4}&0&0&1&0\cr R_{5}&0&0&0&1\cr R_{6}&0&1&-1&1\cr R_{7}&0&1&0&0\cr R_{8}&0&0&1&0\cr R_{10}&0&0&0&1\cr}$$

$$\bordermatrix{&v_{1}^{\prime}&v_{2}^{\prime}&v_{3}^{\prime}&v_{4}^{\prime}\cr w=1&1&0&1&0\cr w=2&1&0&0&0\cr w=3&0&0&1&0\cr w=4&0&0&0&1\cr w=5&0&1&-1&1\cr w=6&0&1&0&0\cr}\bordermatrix{&b_{1}&b_{2}\cr w=1&-1&0\cr w=2&1&0\cr w=3&1&1\cr w=4&0&-1\cr w=5&0&1\cr w=6&0&-1\cr}$$

$$\begin{array}[]{c}199.75{\mu_{{A_{1}}}}-86.03{\mu_{{A_{2}}}}={M_{1}}\\ 159.84{\mu_{{A_{1}}}}-63.32{\mu_{{A_{2}}}}={M_{2}}\\ {\mu_{{A_{4}}}}+1.54{\mu_{{A_{5}}}}={M_{4}}={\mu_{{A_{1}}}}+4.44{\mu_{{A_{2}}}}+11.52{\mu_{{A_{4}}}}\\ -43.80{\mu_{{A_{1}}}}+21.42{\mu_{{A_{4}}}}={M_{5}}\\ -56.13{\mu_{{A_{1}}}}+22.19{\mu_{{A_{4}}}}={M_{6}}\\ {\mu_{{A_{2}}}}<{\mu_{{A_{3}}}}<{\mu_{{A_{4}}}}<{M_{3}}\\ M_{2}<M_{1}<M_{3}\\ M_{5}<M_{6}<M_{3}\end{array}$$

$$\bordermatrix{&A_{1}&A_{2}&A_{3}\cr R_{1}:A_{1}+A_{3}\to A_{1}+A_{2}&0.5&0&0.5\cr R_{2}:A_{1}+A_{2}\to 2A_{3}&1&1&0\cr R_{3}:2A_{3}\to A_{1}+A_{3}&0&0&1\cr R_{4}:A_{3}\to 0&0&0&1\cr R_{5}:A_{3}\to A_{3}+A_{2}&0&0&1\cr}$$

$$\left({\begin{array}[]{*{20}{r}}1&{-1}\\ 1&0\\ 1&0\\ 0&1\\ 0&1\end{array}}\right)$$

$$\begin{array}[]{llll}R_{1}:0\to A_{1}&&&R_{3}:A_{1}\to 2A_{1}\\ R_{2}:A_{1}\to 0&&&R_{4}:2A_{1}\to 0\\ \end{array}$$

$$\left({\begin{array}[]{*{20}{r}}{-1}&2\\ 1&0\\ 0&1\end{array}}\right){\text{ and }}\left({\begin{array}[]{*{20}{r}}1&{-1}&2\\ 1&0&0\\ 0&1&0\\ 0&0&1\\ \end{array}}\right),{\text{ respectively. }}$$

$$f\left({{c^{**}}}\right)=\left({\begin{array}[]{*{20}{c}}1&{-1}&1&-2\\ \end{array}}\right)\left({\begin{array}[]{*{20}{c}}{2{{\left({1}\right)}^{0}}}\\ {1{{\left({1}\right)}^{0.5}}}\\ {1{{\left({1}\right)}^{1}}}\\ {1{{\left({1}\right)}^{0.5}}}\end{array}}\right){\rm\ and\ }f\left({{c^{*}}}\right)=\left({\begin{array}[]{*{20}{c}}1&{-1}&1&-2\\ \end{array}}\right)\left({\begin{array}[]{*{20}{c}}{2{{\left({4}\right)}^{0}}}\\ {1{{\left({4}\right)}^{0.5}}}\\ {1{{\left({4}\right)}^{1}}}\\ {1{{\left({4}\right)}^{0.5}}}\end{array}}\right)$$

$${\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.$$

$$S=span\left\{C_j-C_i\in {\mathbb{R}}^{\mathcal{S}}|\left(C_i,C_j\right)\in \mathcal{R}\right\}.$$

$$f\left(x\right)=NK(x)=\sum _{C_i\to C_j\in \mathcal{R}}{K}_{C_i\to C_j}\left(x\right)\left(C_j-C_i\right).$$

$$E_{+}\left(\mathcal{S},\mathcal{C},\mathcal{R},K\right)=\left\{x\in {\mathbb{R}}_{>0}^{\mathcal{S}}|f\left(x\right)=0\right\}.$$

$${\left({\widetilde{Y}}\right)_{ij}}=\left\{\begin{array}[]{cl}{\left(F\right)_{ki}}&{{\rm{if}}\ j\ {\rm{is\ a\ reactant\ complex\ of\ reaction\ }}k},\\ 0&{\rm otherwise}.\end{array}\right.$$

$$\sum _{y\to y^{\prime }\in \mathcal{R}}{\kappa }_{y\to y^{\prime }}\left(y^{\prime }-y\right)=0\mathrm{and}\sum _{y\to y^{\prime }\in \mathcal{R}}{\kappa }_{y\to y^{\prime }}{e}^{T_{.y}\cdot \mu }\left(y^{\prime }-y\right)=0.$$

Full text

††journal: Mathematical Biosciences

1 Introduction

Biochemical maps are widely used in biochemistry to visualize complex, non-linear interactions between molecules and other chemical entities depicted as nodes. Two types of directed interactions are captured in a biochemical map: transfer of mass (usually represented by solid arrows) from one node to another and regulation (transfer of information, often denoted by dotted line arrows) from a node to a mass transfer arrow. Biochemical maps are akin to the “box models” used in chemical engineering and related disciplines. Figure 1 shows the biochemical map for the anaerobic yeast fermentation pathway model of Galazzo and Bailey [12], and Curto et al. [4], which is provided in [20].

To describe the dynamics of a biochemical system, an interaction is typically assigned a rate function which captures some of its properties such as speed, type (activating or inhibiting) and strength. In Biochemical Systems Theory (BST), power-law functions, which are products of variables with real exponents, are the rate functions of choice. In the half century, since M. Savageau introduced the BST approach, numerous models in many areas of the life sciences have been developed. There are two main types of BST models: GMA (Generalized Mass Action) systems and S-systems. In a GMA system, a power-law rate function is assigned to each arrow in a biochemical map. In an S-system, which formally is a special GMA system, the incoming (outgoing) arrows of a node are “lumped” together to a single incoming (outgoing) arrow and assigned a power-law function each (often viewed as a “weighted average” of the individual arrows in a corresponding GMA system).

S-systems have many advantages for analysis of system dynamics, including a very simple solution to the question of multistationarity. Determining the positive equilibria is equivalent to setting the two power-law functions equal and taking the logarithm results in a system of linear equations in the logarithms of original variables (i.e., in logarithmic space). The entries of the matrix of this system, which is called the “A-matrix” in BST, are the differences of the corresponding exponents (called “kinetic orders”), the nonzero values for which a solution (i.e., an equilibrium) are sought are the ratios of the rate constants. Linear algebra then tells us when there is a unique, infinitely many or no equilibria. For the structurally more accurate GMA systems however, there are very few general results regarding equilibria addressing the question of multistationarity. One typically could compute results for systems with particular numerical parameters.

The potential of Chemical Reaction Network Theory (CRNT) to provide an approach to the multistationarity problem of GMA systems was a motivating factor for Arceo et al. [3, 2] to introduce representations of BST systems as power-law kinetic systems. This involved translating a biochemical map to a full-fledged chemical reaction network and using the concept of “embedded network” to obtain a chemical kinetic system (CKS) dynamically equivalent to the BST system. The analysis of embedded CKS of 15 BST systems (7 of which were GMA systems) revealed that their underlying networks were all non-weakly reversible and have deficiency greater than one (called “higher deficiency” in CRNT). This implied that they did not fall under the coverage of the available results of CRNT on equilibria from Müller-Regensburger [17] and Talabis et al. [19]. In addition, most of the embedded CKS were power-law systems with non-reactant-determined interactions (PL-NDK systems) and the only result for equilibria of PL-NDK systems is a recent Deficiency Zero Theorem from Fortun et al. [10]. The “Multistationarity Algorithm for PLK systems” (MSA) presented in this paper provides a computational solution for the multistationarity problem of GMA Systems.

M. Feinberg pioneered the use of computational approaches to multistationarity in chemical reaction networks in 1995 with the novel Deficiency One Algorithm (DOA) for mass action kinetics (MAK) [7]. It is a significant complement to the Deficiency One Theorem (DOT). In other words, if a reaction network with deficiency one that does not satisfy the conditions of the DOT, one may use the DOA. On the other hand, the higher deficiency algorithm (HDA) for MAK is given in the Ph.D. thesis of H. Ji [15]. It is a reformulation of the works of Feinberg’s DOA [7, 8] and P. Ellison’s Advanced Deficiency Algorithm (ADA) which already handles MAK with higher deficiency [5]. Their results were already implemented in CRNToolbox [16]. When you input a specific reaction network, the software outputs its basic properties. It also reports whether or not the MAK system has the capacity to admit multiple equilibria.

In the MSA, we consider a subset of the reaction set called an orientation, which induces a partition of the subnetwork into equivalence classes and correspondingly, the whole network is partitioned into fundamental classes. A system of equations and/or inequalities will be obtained and checked whether a solution called signature exists. If such a signature exists, then the kinetic system has the capacity for multistationarity for particular rate constants. If no signature exists for all possible construction of the system of equations and/or inequalities, then the kinetic system does not have the capacity for multistationarity, no matter what positive values the rate constants assume.

The multistationarity algorithm MSA consists of two new results: the extension of the HDA of Ji and Feinberg for MAK systems to power-law systems with reactant-determined interactions (PL-RDK systems) and its combination with a transformation method in [18] to address the multistationarity of PL-NDK systems. In particular, we present the following new results:

the monostationarity (for any set of rate constants) of the anaerobic yeast fermentation pathway of Curto et al. (Section 3.1), and 2. 2.

the multistationarity (for certain sets of rate constants) of the global carbon cycle model of Heck et al. (Section 3.2).

In a further new application, we apply the MSA to non-regular deficiency one PL-RDK systems, which are not covered by the extension of Feinberg’s DOA to PL-RDK systems by Fortun et al. [9]. A deficiency one regular network is required to be $t$-minimal and satisfy a “cut pair” condition besides being positive dependent. We also provide an example of a deficiency one PL-NDK system whose transform has deficiency 2, which can still be handled by HDA, but obviously not by DOA.

Remark 1.1**.**

In CRNT, “higher deficiency” refers to networks with deficiency greater than one. We desisted from further use of the term “Higher Deficiency Algorithm” because the algorithm is valid for any deficiency, and in fact, as remarked above, quite useful for deficiency one networks.

2 The Multistationarity Algorithm for Power-Law Systems with Reactant-Determined Interactions

For power-law kinetic systems with reactant-determined interactions, the multistationarity algorithm MSA is the extension of the HDA of Ji and Feinberg to these systems. Using the discussion of the higher deficiency theory and algorithm for MAK in [15], we develop a version for PL-RDK. The algorithm basically follows that of MAK where each reaction $y\to y^{\prime}$ corresponds to a column in the molecularity matrix $Y$. In PL-RDK, the reaction $y\to y^{\prime}$ corresponds to a row in the kinetic order matrix and correspondingly, to a column in the $T$-matrix (see Definition B.13 in Appendix B). With this, we modified the algorithm by considering the $T$-matrix of the kinetic system, for solving system of equations, instead of the molecularity matrix of the reaction network. For a background on CRNT, the reader may refer to Appendix B.

2.1 Parameters for Multistationarity

Let $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ be a PL-RDK system. If the system admits two positive and distinct equilibria, say $c^{*}$ and $c^{**}$, then $\sum\limits_{y\to y^{\prime}\in\mathscr{R}}{{k_{y\to y^{\prime}}}{{\left({{c^{*}}}\right)}^{T_{.y}}}\left({y^{\prime}-y}\right)}=0$ and $\sum\limits_{y\to y^{\prime}\in\mathscr{R}}{{k_{y\to y^{\prime}}}{{\left({{c^{**}}}\right)}^{T_{.y}}}\left({y^{\prime}-y}\right)=0}$ where $T_{.y}$ is the column of the $T$-matrix associated with the reactant complex $y$. Lemmas 2.1 and 2.2, which provide necessary and sufficient conditions for the existence of multiple equilibria, show the importance of the introduction of two new terms $\kappa$ and $\mu$. We modified Lemma 2.2.1 in [15] to come up with these lemmas.

Lemma 2.1**.**

Let $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ be a PL-RDK system. Suppose there exist a set of positive rate constants $\left\{{{k_{y\to y^{\prime}}}|y^{\prime}\to y\in\mathscr{R}}\right\}$ and two distinct, positive, stoichiometrically compatible ${c^{*}},{c^{**}}\in\mathbb{R}_{>0}^{\mathscr{S}}$. Then there exist $\kappa\in\mathbb{R}_{>0}^{\mathscr{R}}$ and $\mu\in{\mathbb{R}^{\mathscr{S}}}$ such that ${\mu_{s}}=\ln\left({\dfrac{{{c_{s}}^{*}}}{{{c_{s}}^{**}}}}\right)\ \forall s\in\mathscr{S}$ and ${\kappa_{y\to y^{\prime}}}={k_{y\to y^{\prime}}}{c}{{}^{**T_{.y}}}$ $\forall$ $y\to y^{\prime}\in\mathscr{R}$ where $\mu$ is both nonzero and stoichiometrically compatible with (the stoichiometric subspace) $S$.

Lemma 2.2**.**

Let $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ be a PL-RDK system. Suppose there exist $\kappa\in\mathbb{R}_{>0}^{\mathscr{R}}$ and $0\neq\mu\in{\mathbb{R}^{\mathscr{S}}}$ which is stoichiometrically compatible with $S$. Then there exist a set of positive rate constants $\left\{{{k_{y\to y^{\prime}}}|y^{\prime}\to y\in\mathscr{R}}\right\}$ and two distinct, positive, stoichiometrically compatible ${c^{*}},{c^{**}}\in\mathbb{R}_{>0}^{\mathscr{S}}$.

The first two equations in this section can be expressed as

[TABLE]

where $\mu=\ln\left({\dfrac{{{c^{*}}}}{{{c^{**}}}}}\right)$. The equivalence can be established by letting ${\mu_{s}}=\ln\left({\dfrac{{{c_{s}}^{*}}}{{{c_{s}}^{**}}}}\right)\ \forall s\ \in\mathscr{S}$. Hence, ${e^{{\mu_{s}}}}=\dfrac{{{c_{s}}^{*}}}{{{c_{s}}^{**}}}\;\forall s\;\in\mathscr{S}$. Then,

[TABLE]

Thus, ${c}{{}^{{*{{T_{.y}}}}}}={e^{{T_{.y}}\cdot\mu}}{c}{{}^{{{**}{{T_{.y}}}}}}$ and ${k_{y\to y^{\prime}}}{c}{{}^{*{{{T_{.y}}}}}}={k_{y\to y^{\prime}}}{e^{{T_{.y}}\cdot\mu}}{c}{{}^{{**}{{T_{.y}}}}}={\kappa_{y\to y^{\prime}}}{e^{{T_{.y}}\cdot\mu}}$ since we set ${\kappa_{y\to y^{\prime}}}={k_{y\to y^{\prime}}}{c}{{}^{**T_{.y}}}$ $\forall$ $y\to y^{\prime}\in\mathscr{R}$.

2.2 Fundamental Classes of Reactions

In this subsection, we introduce a subset of the reaction set $\mathscr{R}$ known as an “orientation”, denoted by $\mathscr{O}$, which will lead to the “unique” partition of $\mathscr{R}$.

Definition 2.3**.**

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]

Remark 2.4**.**

If the network $\mathscr{R}$ has no reversible reactions, then there is only one orientation $\mathscr{O}=\mathscr{R}$. In this case, $L_{\mathscr{R}}=YI_{a}=N$, the stoichiometric matrix.

Illustration 2.5**.**

Consider the GMA model of anaerobic fermentation pathway of Saccharomyces cerevisiae in Chapter 8 of [20] depicted in Figure 1. The model has the following reactions.

[TABLE]

In [2], the system was referred to as ERM0-G and it was shown that the system is PL-RDK. Each reaction is irreversible. Thus, ${\mathscr{O}}=\mathscr{R}.$

Let $\left\{{{\omega_{y\to y^{\prime}}}|y\to y^{\prime}\in{\mathbb{R}^{\mathscr{O}}}}\right\}$ be the standard basis for ${{\mathbb{R}^{\mathscr{O}}}}$. Reactions ${y\to y^{\prime}}$, $\overline{y}\to\overline{y}^{\prime}\in\mathscr{O}$ are related by $\sim$ if there exists a nonzero $\alpha$ such that ${\omega_{y\to y^{\prime}}}-\alpha{\omega_{\overline{y}\to\overline{y}^{\prime}}}\in Ker^{\perp}L_{\mathscr{O}}$. It is easy to see that the relation $\sim$ is indeed an equivalence relation. We define the equivalence class of a reaction $y\to y^{\prime}$ under the relation $\sim$ by $\left[y\to y^{\prime}\right]=\{{\widetilde{y}\to\widetilde{y}^{\prime}}\in\mathscr{O}|{\widetilde{y}\to\widetilde{y}^{\prime}}\sim{y\to y^{\prime}}\}$. We call the equivalence classes as $P_{i}$ where $i$ runs through the number of classes induced by the relation. In particular, if there is a reaction such that ${\omega_{y\to y^{\prime}}}\in Ker^{\perp}L_{\mathscr{O}}$, we label the equivalence class containing $y\to y^{\prime}$ as $P_{0}$. We set $P_{0}=\emptyset$ if there is no such equivalence class exists.

Lemma 2.6**.**

[15]** Let $\left(\mathscr{S},\mathscr{C},\mathscr{R}\right)$ be reaction network and $\mathscr{O}$ be an orientation. Let $\alpha\neq 0$. Then the following holds:

• i.

$y\to y^{\prime}\in{P_{0}}$* if and only if ${x_{y\to y^{\prime}}}=0$ for each $x\in KerL_{\mathscr{O}}$, and*

• ii.

${\omega_{y\to y^{\prime}}}-\alpha{\omega_{\overline{y}\to\overline{y}^{\prime}}}\in Ker^{\perp}L_{\mathscr{O}}$* if and only if ${x_{y\to y^{\prime}}}=\alpha{x_{\overline{y}\to\overline{y}^{\prime}}}\ \forall x\in KerL_{\mathscr{O}}$.*

Lemma 2.6 leads to partitioning the orientation $\mathscr{O}$ into the zeroth equivalence class $P_{0}$, if it exists, and the nonzeroth equivalence classes $P_{i}$ where $i\geq 1$.

Remark 2.7**.**

[15*]**

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 belong to the same equivalence class.*

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 statements holds [15].

• 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 label the fundamental class containing $P_{i}$ as $C_{i}$. A fundamental (equivalence) class is said to be reversible if each reaction in the fundamental (equivalence) class is reversible (with respect to $\mathscr{R}$). If at least one of the reactions is irreversible, then the class is said to be nonreversible. If $P_{i}$ is nonreversible, we pick any irreversible reaction to be the representative of $P_{i}$. Otherwise, we pick any irreversible reaction as the representative. We identify the $i$th reaction in $W$ with ${y_{i}}\to{y_{i}}^{\prime}$. We let $W$ be the collection of all such representatives from $P_{i}$ where $i=1,2,...,w$.

Illustration 2.8**.**

We consider ERM0-G in Illustration 2.5. Below is basis for $Ker{L_{\mathscr{O}}}$ obtained by solving $\sum\limits_{y\to y^{\prime}\in\mathscr{O}}{{\alpha_{y\to y^{\prime}}}\left({y^{\prime}-y}\right)}=0$.

[TABLE]

Notice that there is no row with zero entry so $P_{0}$ is empty. We partition $\mathscr{O}$ into its equivalence classes. The equivalence classes and the fundamental classes based on the given basis for $Ker{L_{\mathscr{O}}}$ are provided in the first and second columns of Table 1, respectively.

A fundamental class $C_{i}$ with $0\leq i\leq w$ is said to be degenerate if $g_{y_{i}\to y_{i}^{\prime}}=0$ while a fundamental class $C_{i}$ with $1\leq i\leq w$ is said to be nondegenerate if $g_{y_{i}\to y_{i}^{\prime}}\neq 0$. Let $i\geq 1$. For each nondegenerate fundamental class $C_{i}$, we assume a 3-shelf bookcase to store all reactions in $C_{i}$. Let $y\to y^{\prime}$ be a reaction in a nondegenerate fundamental class $C_{i}$, and ${\rho_{{y_{i}}\to{y_{i}}^{\prime}}}=\dfrac{{{h_{{y_{i}}\to{y_{i}}^{\prime}}}}}{{{g_{{y_{i}}\to{y_{i}}^{\prime}}}}}$ where ${g_{{y_{i}}\to{y_{i}}^{\prime}}}\neq 0$ and $i=1,...,w$. Then we define the shelving of $y\to y^{\prime}$ in the following manner.

• i.

$y\to y^{\prime}$ is on the upper shelf if ${e^{{T_{.y}}\cdot\mu}}>\rho_{y_{i}\to y_{i}^{\prime}}$.

• ii.

$y\to y^{\prime}$ is on the lower shelf if ${e^{{T_{.y}}\cdot\mu}}<\rho_{y_{i}\to y_{i}^{\prime}}$.

• iii.

$y\to y^{\prime}$ is on the middle shelf if ${e^{{T_{.y}}\cdot\mu}}=\rho_{y_{i}\to y_{i}^{\prime}}$.

Let ${M_{{y_{i}}\to{y_{i}}^{\prime}}}=\ln{\rho_{{y_{i}}\to{y_{i}}^{\prime}}}$ if ${\rho_{{y_{i}}\to{y_{i}}^{\prime}}}>0$. Otherwise, we take ${M_{{y_{i}}\to{y_{i}}^{\prime}}}$ to be an arbitrarily large and negative number. We restrict the problem by considering a representative of each of the classes instead of all the elements of the whole class. Lemmas 2.9 and 2.10 are extensions of Lemma 2.8.7 and Proposition 2.8.1 in [15] from MAK to PL-RDK.

Lemma 2.9**.**

*Let $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ be a PL-RDK system and $\mathscr{O}$ be an orientation. Let $\kappa\in\mathbb{R}_{>0}^{\mathscr{R}}$, and $\mu\in{\mathbb{R}^{\mathscr{S}}}$. Let $g,h\in\mathbb{R}^{\mathscr{O}}$ such that

{g_{y\to y^{\prime}}}=\left\{\begin{array}[]{ll}{\kappa_{y\to y^{\prime}}}-{\kappa_{y^{\prime}\to y}}&{\rm{if}}\ y\to y^{\prime}\in\mathscr{O}{\rm{\ is\ reversible}}\\ {\kappa_{y\to y^{\prime}}}&{\rm{if}}\ y\to y^{\prime}\in\mathscr{O}{\rm{\ is\ irreversible}}\end{array}\right. and

{h_{y\to y^{\prime}}}=\left\{\begin{array}[]{ll}{\kappa_{y\to y^{\prime}}}{e^{{T_{.y}}\cdot\mu}}-{\kappa_{y^{\prime}\to y}}{e^{{T_{.y^{\prime}}}\cdot\mu}}&{\rm{if}}\ y\to y^{\prime}\in\mathscr{O}{\rm{\ is\ reversible}}\\ {\kappa_{y\to y^{\prime}}}{e^{{T_{.y}}\cdot\mu}}&{\rm{if}}\ y\to y^{\prime}\in\mathscr{O}{\rm{\ is\ irreversible}}\end{array}\right..

For $i=1,2,...,w$, let ${P_{i}}$ be the equivalence class with ${y_{i}}\to{y_{i}}^{\prime}$ as representative. Moreover, let ${{\rho_{{y_{i}}\to{y_{i}}^{\prime}}}=\dfrac{{{h_{{y_{i}}\to{y_{i}}^{\prime}}}}}{{{g_{{y_{i}}\to{y_{i}}^{\prime}}}}}}$ for nondegenerate fundamental classes ${C_{i}}$.*

• i.

If $y\to y^{\prime}\in{P_{i}}\ (i=1,2,...,w)$ is irreversible then ${g_{y\to y^{\prime}}}>0$, ${h_{y\to y^{\prime}}}>0$, and ${M_{{y_{i}}\to{y_{i}}^{\prime}}}={{{T_{.y}}\cdot\mu}}$.

• ii.

Suppose $y\to y^{\prime}\in{P_{i}}\ (i=1,2,...,w)$ is reversible and ${C_{i}}$ is nondegenerate.

• a.

If ${g_{y_{i}\to y_{i}^{\prime}}}>0$ and $y\to y^{\prime}$ is on the upper shelf, then ${M_{{y_{i}}\to{y_{i}}^{\prime}}}<{{{T_{.y}}\cdot\mu}}<{{{T_{.y^{\prime}}}\cdot\mu}}$.

• b.

If ${g_{y_{i}\to y_{i}^{\prime}}}>0$ and $y\to y^{\prime}$ is on the middle shelf, then ${M_{{y_{i}}\to{y_{i}}^{\prime}}}={{{T_{.y}}\cdot\mu}}={{{T_{.y^{\prime}}}\cdot\mu}}$.

• c.

If ${g_{y_{i}\to y_{i}^{\prime}}}>0$ and $y\to y^{\prime}$ is on the lower shelf, then ${M_{{y_{i}}\to{y_{i}}^{\prime}}}>{{{T_{.y}}\cdot\mu}}>{{{T_{.y^{\prime}}}\cdot\mu}}$.

• d.

If ${g_{y_{i}\to y_{i}^{\prime}}}<0$ and $y\to y^{\prime}$ is on the upper shelf, then ${M_{{y_{i}}\to{y_{i}}^{\prime}}}<{{{T_{.y^{\prime}}}\cdot\mu}}<{{{T_{.y}}\cdot\mu}}$.

• e.

If ${g_{y_{i}\to y_{i}^{\prime}}}<0$ and $y\to y^{\prime}$ is on the middle shelf, then ${M_{{y_{i}}\to{y_{i}}^{\prime}}}={{{T_{.y^{\prime}}}\cdot\mu}}={{{T_{.y}}\cdot\mu}}$.

• f.

If ${g_{y_{i}\to y_{i}^{\prime}}}<0$, and $y\to y^{\prime}$ is on the lower shelf, then ${M_{{y_{i}}\to{y_{i}}^{\prime}}}>{{{T_{.y^{\prime}}}\cdot\mu}}>{{{T_{.y}}\cdot\mu}}$.

• iii.

Suppose $y\to y^{\prime}\in{P_{i}}\ (i=1,2,...,w)$ is reversible and ${C_{i}}$ is degenerate.

• a.

If ${h_{y_{i}\to y_{i}^{\prime}}}>0$ then ${{{T_{.y}}\cdot\mu}}>{{{T_{.y^{\prime}}}\cdot\mu}}$.

• b.

If ${h_{y_{i}\to y_{i}^{\prime}}}=0$ then ${{{T_{.y}}\cdot\mu}}={{{T_{.y^{\prime}}}\cdot\mu}}$.

• c.

If ${h_{y_{i}\to y_{i}^{\prime}}}<0$ then ${{{T_{.y}}\cdot\mu}}<{{{T_{.y^{\prime}}}\cdot\mu}}$.

• iv.

If $y\to y^{\prime}\in{P_{0}}$ is reversible then ${{{T_{.y}}\cdot\mu}}={{{T_{.y^{\prime}}}\cdot\mu}}$.

Lemma 2.10**.**

Let $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ be a PL-RDK system and $\mathscr{O}$ be an orientation. Suppose there exist $\mu\in{\mathbb{R}^{\mathscr{S}}}$, $g,h\in KerL_{\mathscr{O}}$, ${P_{i}}\ (i=1,2,...,w)$ with representative ${{y_{i}}\to{y_{i}}^{\prime}}$ and $\left\{{{\rho_{{y_{i}}\to{y_{i}}^{\prime}}}=\dfrac{{{h_{{y_{i}}\to{y_{i}}^{\prime}}}}}{{{g_{{y_{i}}\to{y_{i}}^{\prime}}}}}|{g_{{y_{i}}\to{y_{i}}^{\prime}}}\neq 0,i=1,...,w}\right\}$ satisfy the conditions given in Lemma 2.9. Then the following holds for a nondegenerate $C_{i}$.

• i.

All irreversible reactions in $C_{i}\ (i\geq 1)$ must belong to the middle shelf.

• ii.

$y\to y^{\prime}\in{C_{i}}$* must belong to the upper shelf if ${{\rho_{{y_{i}}\to{y_{i}}^{\prime}}}}\leq 0$.*

• iii.

If a reaction is reversible, then the reaction and its reversible pair must belong to the same shelf.

• iv.

Any two reactions in $C_{i}$ with the same reactant complex must belong to the same shelf.

• v.

Each reaction whose reactant complex lies in a nonterminal strong linkage class must belong to the middle shelf.

• vi.

Each reaction whose reactant complex lies in a terminal strong linkage class of the fundamental subnetwork must belong to the same shelf.

• vii.

If for a nondegenerate $C_{i}\ (i\geq 1)$, $\mathscr{N}_{i}$ forms a big (undirected) cycle (with at least three vertices), then its reactions are all in a terminal strong linkage class and belong to the middle shelf, where $\mathscr{N}_{i}$ is the subnetwork generated by reactions in $P_{i}$.

Illustration 2.11**.**

*Again, consider ERM0-G in Illustration 2.5. Since irreversible reactions belong to the middle shelf, we obtain the following shelving assignment.

${{\cal U}_{1}}=\left\{{}\right\},{{\cal M}_{1}}=\left\{{R_{1}:X_{2}\to X_{1}+X_{2},R_{2}:X_{1}+X_{5}\to X_{2}+X_{5}}\right\},{{\cal L}_{1}}=\left\{{}\right\}$

${{\cal U}_{2}}=\left\{{}\right\},{{\cal M}_{2}}=\left\{{R_{3}:2X_{5}+X_{1}\to X_{5}+X_{1}}\right\},{{\cal L}_{2}}=\left\{{}\right\}$

${{\cal U}_{3}}=\left\{{}\right\},{{\cal M}_{3}}=\left\{{R_{4}:X_{2}+X_{5}\to X_{3}+X_{5}}\right\},{{\cal L}_{3}}=\left\{{}\right\}$

${{\cal U}_{4}}=\left\{{}\right\},{{\cal M}_{4}}=\left\{{R_{5}:2X_{5}+X_{2}\to X_{5}+X_{2}}\right\},{{\cal L}_{4}}=\left\{{}\right\}$

${{\cal U}_{5}}=\left\{{}\right\},{{\cal M}_{5}}=\left\{{R_{6}:X_{2}+X_{5}\to X_{5}}\right\},{{\cal L}_{5}}=\left\{{}\right\}$

${{\cal U}_{6}}=\left\{{}\right\},{{\cal M}_{6}}=\left\{{R_{7}:X_{2}+X_{5}\to X_{2}}\right\},{{\cal L}_{6}}=\left\{{}\right\}$

${{\cal U}_{7}}=\left\{{}\right\},{{\cal M}_{7}}=\left\{{R_{8}:X_{3}+X_{5}\to X_{4}+X_{5},R_{11}:X_{3}+X_{4}+X_{5}\to X_{3}+X_{5}}\right\},\\ {{\cal L}_{7}}=\left\{{}\right\}$

${{\cal U}_{8}}=\left\{{}\right\},{{\cal M}_{8}}=\left\{{R_{9}:X_{3}+X_{5}\to X_{3}+2X_{5}}\right\},{{\cal L}_{8}}=\left\{{}\right\}$

${{\cal U}_{9}}=\left\{{}\right\},{{\cal M}_{9}}=\left\{{R_{10}:X_{3}+X_{4}+X_{5}\to X_{4}+X_{5}}\right\},{{\cal L}_{9}}=\left\{{}\right\}$

${{\cal U}_{10}}=\left\{{}\right\},{{\cal M}_{10}}=\left\{{R_{12}:X_{3}+X_{4}+X_{5}\to X_{3}+X_{4}+2X_{5}}\right\},{{\cal L}_{10}}=\left\{{}\right\}$

${{\cal U}_{11}}=\left\{{}\right\},{{\cal M}_{11}}=\left\{{R_{13}:2X_{5}\to X_{5}}\right\},{{\cal L}_{11}}=\left\{{}\right\}$*

2.3 Sign Patterns and the Fundamental Theorem of Multistationarity in PL-RDK systems

Suppose $g,h\in Ker{L_{\mathscr{O}}}$ such that $g\neq 0$ but $h=0$. For any nonzero sign pattern which is stoichiometrically compatible with $Ker{L_{\mathscr{O}}}$, we have a solution of nonzero $g\in Ker{L_{\mathscr{O}}}$ with such sign patterns. Define

[TABLE]

with $W=\left\{{{y_{i}}\to{y_{i}}^{\prime}|i=1,...,w}\right\}\in\mathscr{O}$.

We now focus our attention in depicting what is meant by a “valid” pair of sign patterns for ${g_{W}}=g|_{W},{h_{W}}=h|_{W}\in{\mathbb{R}^{\mathscr{O}}}\cap\Gamma_{W}$. A pair of sign patterns for ${g_{W}},{h_{W}}\in{\mathbb{R}^{\mathscr{O}}}\cap\Gamma_{W}$ is said to be compatible with $Ker{L_{\mathscr{O}}}{|_{W}}$ if both of these statements hold.

• i.

If ${C_{i}}\left({i=1,2,...,w}\right)$ is nonreversible then the signs of ${g_{W}}\left({y\to y^{\prime}}\right)$ and ${h_{W}}\left({y\to y^{\prime}}\right)$ are both positive.

• ii.

For every reaction ${y\to y^{\prime}}$ in $W$, $v_{{y_{i}}\to{y_{i}}^{\prime}}^{1}$ and ${g_{W}}\left({y\to y^{\prime}}\right)$ have the same sign, and $v_{{y_{i}}\to{y_{i}}^{\prime}}^{2}$ and ${h_{W}}\left({y\to y^{\prime}}\right)$ have the same sign, for some ${v^{1}},{v^{2}}\in Ker{L_{\mathscr{O}}}$.

In other words, a pair of sign patterns for ${g_{W}}$ and ${h_{W}}$ are said to be valid, if it is nonzero and stoichiometrically compatible with $Ker{L_{\mathscr{O}}}{|_{W}}$.

At this point, we assume that there is a sign pattern for ${g_{W}}\in{\mathbb{R}^{\mathscr{O}}}\cap{\Gamma_{W}}$. Define the sets

[TABLE]

and

[TABLE]

Let $\{b^{j}\}_{j=1}^{q}$ be a basis for $Ker^{\perp}{L_{\mathscr{O}}}\cap{\Gamma_{W}}$, if it exists, where $q$ is the dimension of $Ker^{\perp}{L_{\mathscr{O}}}\cap{\Gamma_{W}}$. For $j=1,2,...,q$, we also define the following sets and consider the following equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Lemma 2.12**.**

[15*]** Suppose a reaction network satisfies the following properties for an orientation $\mathscr{O}$: $P_{i}$ ($i=0,1,2,...,w$) is defined by a representative ${{y_{i}}\to{y_{i}}^{\prime}}$, $W=\{{{y_{i}}\to{y_{i}}^{\prime}}|i=1,2,...,w\}\subseteq\mathscr{O}$, a given basis $\left\{{{b^{j}}}\right\}_{j=1}^{q}$ for $Ker^{\perp}{L_{\mathscr{O}}}\cap{\Gamma_{W}}$, and a set of parameters

$\left\{{{\rho_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)|{g_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)\neq 0,i=1,2,...,w}\right\}$ where the sign of ${{\rho_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}$ is the same as the ratio of the signs of ${{h_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}$ and ${{g_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}$. Further, suppose there exist ${h_{W}},{g_{W}}\in{\mathbb{R}^{\mathscr{O}}}\cap{\Gamma_{W}}$ with a valid pair of sign patterns such that Equations (1), (2), and (3) are satisfied. Then the following holds for $j=1,2,...,q$:*

• i.

If $\sum\limits_{{y_{i}}\to{y_{i}}^{\prime}\in D}{b_{{y_{i}}\to{y_{i}}^{\prime}}^{j}{h_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)>0}$, then one element in $Q_{2}^{j}$ is strictly greater than one element in $Q_{1}^{j}$.

• ii.

If $\sum\limits_{{y_{i}}\to{y_{i}}^{\prime}\in D}{b_{{y_{i}}\to{y_{i}}^{\prime}}^{j}{h_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)<0}$, then one element in $Q_{1}^{j}$ is strictly greater than one element in $Q_{2}^{j}$.

• iii.

If $\sum\limits_{{y_{i}}\to{y_{i}}^{\prime}\in D}{b_{{y_{i}}\to{y_{i}}^{\prime}}^{j}{h_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)=0}$, then $Q_{1}^{j}$ and $Q_{2}^{j}$ are nonsegregated. That is, at least one of these holds:

• a.

There exist an element $a$ from one of multisets $Q_{1}^{j}$ and $Q_{2}^{j}$, and $b<c$ from the other such that a is between $b$ and $c$, i.e., $b<a<c$.

• b.

All the elements in the multisets $Q_{1}^{j}$ and $Q_{2}^{j}$ are equal or there are $a,b\in Q_{1}^{j}$ and $c,d\in Q_{2}^{j}$ where $c=a<b=d$.

Definition 2.13**.**

Let $\left(\mathscr{S},\mathscr{C},\mathscr{R}\right)$ be a chemical reaction network with an orientation $\mathscr{O}$. The chemical reaction network is said to have a forestal property if for given $W\subseteq\mathscr{O}$, $Ker^{\perp}{L_{\mathscr{O}}}\cap{\Gamma_{W}}$ has a forest basis. $Ker^{\perp}{L_{\mathscr{O}}}\cap{\Gamma_{W}}$ has a forest basis if it has a basis such that the graph based on the basis vectors is a forest graph. (The reader may refer to [15] for further details on forest graphs.)

The following theorem is an extension of Theorem 2.11.9 in [15] to PL-RDK system which we shall call the “Fundamental Theorem of Multistationarity in PL-RDK systems”.

Theorem 2.14**.**

Suppose the reaction network $\left(\mathscr{S},\mathscr{C},\mathscr{R}\right)$ has a forestal property and $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ is a PL-RDK system. It has the capacity to admit multiple steady states if and only if the following statements hold:

• i.

$0\neq\mu\in{\mathbb{R}^{\mathscr{S}}}$* exists which is stoichiometrically compatible with $S$,*

• ii.

a valid sign pattern for ${g_{W}},{h_{W}}\in\mathbb{R}^{\mathscr{O}}\cap{\Gamma_{W}}$ exists,

• iii.

a set of parameters $\left\{{{\rho_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)|{g_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)\neq 0}\right\}$ where the sign of ${\rho_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)$ is precisely the ratio of the signs of ${{h_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}$ and

${{g_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}$* that satisfies the conditions in Lemma 2.12, and*

• iv.

a shelving assignment exists for each nondegenerate fundamental class that satisfies the conditions in Lemma 2.10

which together satisfy the conditions in Lemma 2.9 in terms of ${{g_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}$, ${{h_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}$, and ${{\rho_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}$.

3 Applications of the MSA to GMA systems

This section emphasizes how powerful the MSA is, as we present a solution to the problem of monostationarity, for any set of rate constants, of the model of anaerobic fermentation pathway in yeast, and the multistationarity, for particular set of rate constants, of the global carbon cycle model of Heck et al. which are not yet known in literature. Using the algorithm, we can determine whether a PL-RDK, with underlying network of any deficiency, has the capacity to admit multiple steady states.

3.1 Application to the Model of Anaerobic Fermentation Pathway in Yeast

The fermentation pathway in a species of yeast, known as Saccharomyces cerevisiae, has been studied extensively. In particular, Galazzo and Bailey [12, 13] established the experimental basis, and they were able to provide kinetic equations which was used by Curto, Cascante and Sorribas [4] to derive GMA (and S-system) models and performed standard procedures of biochemical systems analysis. It was used as a case study in Chapter 8 of [20] which demonstrates different numerical analyses for biochemical modeling but was not focused on the algebraic aspects of steady-state analyses.

The fermentation pathway is a popular example to be found in numerous BST papers, yet the basic question of its monostationarity or multistationarity is unknown. Different experimental set-ups and conditions under which yeast cells produce ethanol have been investigated which is of high importance for industrial purposes. The pathway in Figure 1 describes how yeast can use glucose to produce ethanol, and also glycerol, glycogen and trehalose [20].

We derive a novel result about the GMA model of anaerobic fermentation pathway of yeast: its monostationarity, i.e., it has at most one steady state for any set of rate constants. The reaction network is given in Illustration 2.5 and its network numbers is provided in Table 2 [2]. The $T$-matrix is given below.

\bordermatrix{&X_{2}&X_{1}+X_{5}&2X_{5}+X_{1}&X_{2}+X_{5}&2X_{5}+X_{2}&X_{3}+X_{5}&X_{3}+X_{4}+X_{5}&2X_{5}\cr X_{1}&0&0.7464&0.7464&0&0&0&0&0\cr X_{2}&-0.2344&0&0&8.6107&0.7318&0&0&0\cr X_{3}&0&0&0&0&0&0.6159&0.05&0\cr X_{4}&0&0&0&0&0&0&0.533&0\cr X_{5}&0&0.0243&0.0243&0&-0.3941&0.1308&-0.0822&1\cr}

STEP 1: CHOOSING AN ORIENTATION

In this step, we refer to Illustration 2.5.

STEP 2: FINDING EQUIVALENCE CLASSES AND FUNDAMENTAL CLASSES

In this step, we refer to Illustration 2.8. If one of these statements (a) and (b) is not satisfied, then the system does not have the capacity to admit multiple equilibria, and we exit the algorithm.

• (a)

All reactions in $P_{0}$ are reversible (with respect to $\mathscr{R}$).

• (b)

For two irreversible reactions (with respect to $\mathscr{R}$), $y\to y^{\prime}$ and ${\overline{y}\to\overline{y}^{\prime}}$ in the same $P_{i}$, there exists $\alpha>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$.

STEP 3: FINDING THE COLINKAGE SETS

In this step, we divide the reaction network into subnetworks in such a way that all reactions belonging from the same fundamental class are in the same subnetwork. We again refer to Table 1.

STEP 4: PICKING $W\subseteq{\mathscr{O}}$

Recall that an equivalence class is reversible if all of its reactions are reversible with respect to the original network $\mathscr{R}$. However, it is nonreversible, if it contains an irreversible reaction. We pick a representative reaction for each of the $P_{i}$’s such that if a class $P_{i}$ is nonreversible, we pick an irreversible reaction. Otherwise, we pick any reversible reaction. The collection of all the representatives from $P_{i}$ where $i=1,2,...,w$ is the set $W$ with $w$ as its number of elements.

In our example, since there are 11 equivalence classes, $w=11$. For $P_{1}$, we choose $R_{1}$. Similarly, we choose $R_{8}$ for $P_{7}$. Since the rest of each $P_{i}$ has only one element, we have no choice but to choose these reactions as representatives. The $i^{\rm th}$ reaction in $W$ is identified as ${y_{i}}\to{y_{i}}^{\prime}$. For instance, $R_{8}$ will be identified as $y_{7}\to y_{7}^{\prime}$.

STEP 5: REALIGNING THE ORIENTATION (if needed)

For each nonzeroth equivalence class $P_{i}$ with $1\leq i\leq w$, for any reaction ${y\to y^{\prime}}$ in $P_{i}$, there exists an $\alpha_{y\to y^{\prime}}>0$ such that $v_{y_{i}\to y^{\prime}_{i}}^{l}=\alpha_{y\to y^{\prime}}v_{y\to y^{\prime}}^{l}$ for all the basis elements $v_{1},v_{2},...,v_{d}$ (for $KerL_{\mathscr{O}}$). If this statement does not hold, we then realign the orientation (or choose another orientation) until it is already satisfied. In the given example, the statement is satisfied so we go to the next step.

STEP 6: FINDING A BASIS FOR $Ker^{\perp}L_{\mathscr{O}}\cap\Gamma_{W}$

To simplify this step, from STEP 2, we just consider the rows of the reactions in $W$. Moreover, a basis for $Ker^{\perp}L_{\mathscr{O}}\cap\Gamma_{W}$ is also given below. \bordermatrix{&v_{1}^{\prime}&v_{2}^{\prime}&v_{3}^{\prime}&v_{4}^{\prime}&v_{5}^{\prime}&v_{6}^{\prime}&v_{7}^{\prime}&v_{8}^{\prime}\cr w=1&0&1&0&0&1&1&0&0\cr w=2&-1&0&-1&1&0&0&1&-1\cr w=3&0&0&0&0&1&1&0&0\cr w=4&1&0&0&0&0&0&0&0\cr w=5&0&1&0&0&0&0&0&0\cr w=6&0&0&1&0&0&0&0&0\cr w=7&0&0&0&0&0&1&0&0\cr w=8&0&0&0&1&0&0&0&0\cr w=9&0&0&0&0&1&0&0&0\cr w=10&0&0&0&0&0&0&1&0\cr w=11&0&0&0&0&0&0&0&1\cr}  \bordermatrix{&a_{1}&a_{2}&a_{3}\cr w=1&-1&0&0\cr w=2&0&0&1\cr w=3&1&-1&0\cr w=4&0&0&1\cr w=5&1&0&0\cr w=6&0&0&1\cr w=7&0&1&0\cr w=8&0&0&-1\cr w=9&0&1&0\cr w={10}&0&0&-1\cr w={11}&0&0&1\cr}

STEP 7: CHECKING THE LINEARITY OF THE SYSTEM OF INEQUALITIES

If there exists a forest basis for $Ker^{\perp}L_{\mathscr{O}}\cap\Gamma_{W}$, then the resulting inequality system is linear. If it does not exist, we may need additional nonlinear equations on the $M_{i}$’s to determine if the kinetic system has the capacity to admit multiple equilibria.

STEP 8: CHOOSING SIGNS FOR ${g_{W}},{h_{W}}\in{\mathbb{R}}^{\mathscr{O}}\cap{\Gamma_{W}}$

Since each $P_{i}$ is nonreversible, the sign patterns for ${g_{W}}$ and ${h_{W}}$ must be positive.

STEP 9: SHELVING REACTIONS IN THE NONDEGENERATE $C_{i}$’s

We refer to Illustration 2.11.

STEP 10: SHELVING EQUALITIES AND INEQUALITIES FROM THE NONDEGENERATE $C_{i}$’s

We refer to Lemma 2.9. For each nondegenerate $C_{i}$, if ${{\rho_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}>0$, let ${M_{i}}=\ln\left({{\rho_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}\right)$. Otherwise, let ${M_{i}}$ to be an arbitrary large and negative number. Suppose $C_{i}$ is a nondegenerate fundamental class. If $y\to y^{\prime}\in{C_{i}}$ is on the middle shelf, then $T_{.y}\cdot\mu={M_{i}}$ is added to the system. If the given reaction is on the upper shelf, then $T_{.y}\cdot\mu>{M_{i}}$ is added to the system. If the given reaction is on the lower shelf, then $T_{.y}\cdot\mu<{M_{i}}$ is added to the system. Now, from a reversible reaction $y\to y^{\prime}$ of $P_{0}$, it is automatic that $T_{.y}\cdot\mu=T_{.y^{\prime}}\cdot\mu$ should be added to the system.

$-0.2344{\mu_{X_{2}}}=M_{1}=M_{2}=0.7464{\mu_{X_{1}}}+0.0243{\mu_{X_{5}}}$

$0.7318{\mu_{X_{2}}}-0.3941{\mu_{X_{5}}}=M_{3}=M_{4}$

$8.6107{\mu_{X_{2}}}=M_{5}=M_{6}$

$0.6159{\mu_{X_{3}}}+0.1308{\mu_{X_{5}}}=M_{7}=M_{8}=M_{9}=M_{10}=0.05{\mu_{X_{3}}}+0.533{\mu_{X_{4}}}-0.0822{\mu_{X_{5}}}$

${\mu_{X_{5}}}=M_{11}$

STEP 11: UPPER AND LOWER SHELVING INEQUALITIES FROM $P_{i}$’s WITH NONDEGENERATE $C_{i}$’s

Suppose $C_{i}$ is nondegenerate and let $y\to y^{\prime}$ in $P_{i}$. If ${g_{W}}\left({y_{i}\to y_{i}^{\prime}}\right)>0$ and ${y\to y^{\prime}}$ is on the upper shelf, or if ${g_{W}}\left({y_{i}\to y_{i}^{\prime}}\right)<0$ and ${y\to y^{\prime}}$ is on the lower shelf, then we add $T_{.y}\cdot\mu<T_{.y^{\prime}}\cdot\mu$ to the system. However, if ${g_{W}}\left({y_{i}\to y_{i}^{\prime}}\right)>0$ and ${y\to y^{\prime}}$ is on the lower shelf, or if ${g_{W}}\left({y_{i}\to y_{i}^{\prime}}\right)<0$ and ${y\to y^{\prime}}$ is on the upper shelf, then we add $T_{.y}\cdot\mu>T_{.y^{\prime}}\cdot\mu$ to the system.

We skip this step since the upper and the lower shelves are both empty.

STEP 12: ADDING EQUALITIES AND INEQUALITIES FROM $P_{i}$’s WITH DEGENERATE $C_{i}$’s

Suppose $C_{i}$ is degenerate and let $y\to y^{\prime}$ in $P_{i}$. If ${h_{W}}\left({y_{i}\to y_{i}^{\prime}}\right)>0$ then we add $T_{.y}\cdot\mu>T_{.y^{\prime}}\cdot\mu$ to the system. However, if ${h_{W}}\left({y_{i}\to y_{i}^{\prime}}\right)=0$ then we add $T_{.y}\cdot\mu=T_{.y^{\prime}}\cdot\mu$ to the system. Otherwise, we add $T_{.y}\cdot\mu<T_{.y^{\prime}}\cdot\mu$ to the system.

Our example has no degenerate $C_{i}$.

STEP 13: ADDING $M$ EQUALITIES AND INEQUALITIES

From STEP 6, we obtained a basis for $Ker^{\perp}L_{\mathscr{O}}\cap\Gamma_{W}$. In particular, $b^{1}=(-1\ 0\ 1\ 0\ 1\ 0\ 0\ 0\ 0\ 0\ 0)$, $b^{2}=(0\ 0\ -1\ 0\ 0\ 0\ 1\ 0\ 1\ 0\ 0)$, and $b^{3}=(0\ 1\ 0\ 1\ 0\ 1\ 0\ -1\ 0\ -1\ 1)$. From $b^{1}$, we have $R^{1}_{+}=\{y_{3}\to y_{3}^{\prime},y_{5}\to y_{5}^{\prime}\}$ and $R^{1}_{-}=\{y_{1}\to y_{1}^{\prime}\}$. From $b^{2}$, we have $R^{2}_{+}=\{y_{7}\to y_{7}^{\prime},y_{9}\to y_{9}^{\prime}\}$ and $R^{2}_{-}=\{y_{3}\to y_{3}^{\prime}\}$. From $b^{3}$, we have $R^{3}_{+}=\{y_{2}\to y_{2}^{\prime},y_{4}\to y_{4}^{\prime},y_{6}\to y_{6}^{\prime},y_{11}\to y_{11}^{\prime}\}$ and $R^{3}_{-}=\{y_{8}\to y_{8}^{\prime},y_{10}\to y_{10}^{\prime}\}$.

Since our example has no degenerate $C_{i}$, we only consider the “nonsegregated case”. We take the required sum as zero if the degenerate set is empty. Recall that two multisets $Q_{1}$ and $Q_{2}$ are nonsegregated if at least one of the following two cases holds:

• i.

There exist $a$ from $Q_{1}$ or $Q_{2}$, and $b<c$ from the other such that $b<a<c$.

• ii.

All elements in the two multisets are equal, or there exist $a,b\in Q_{1}$ and $c,d\in Q_{2}$ such that $c=a<b=d$.

From STEP 10, we already have $M_{1}=M_{2}$, $M_{3}=M_{4}$, $M_{5}=M_{6}$, and $M_{7}=M_{8}=M_{9}=M_{10}$. From $b^{2}$, we obtain $M_{7}=M_{3}=M_{9}$.

CASE 1: From $b^{1}$, we set $M_{3}<M_{1}<M_{5}$. Hence, from $b^{3}$, it can only be $M_{2}>M_{8}>M_{11}$ or $M_{6}>M_{8}>M_{11}$.

CASE 2: From $b^{1}$, we set $M_{3}>M_{1}>M_{5}$. Thus, from $b^{3}$, it can only be $M_{2}<M_{8}<M_{11}$ or $M_{6}<M_{8}<M_{11}$.

CASE 3: From $b^{1}$, we set $M_{3}=M_{1}=M_{5}$. Finally, from $b^{3}$, $M_{i}=M_{j}$ for all $i,j\in\{1,2,...,11\}$.

STEP 14: CHECKING FOR SOLUTION OF THE SYSTEM

If the system of equations and inequalities is linear, then the system obtained from STEP 10 to STEP 13 is complete. If it has a solution, then it is called a signature. However, if the system is nonlinear, then additional nonlinear constraints are needed to make it complete. If a partial linear system we obtained from STEP 10 to STEP 13 has a solution, then it is called a pre-signature. If none of such inequality systems has a solution, then additional nonlinear constraints are not needed, and we conclude that the system does not have the capacity to admit multiple equilibria.

This step checks when one obtains a $\mu$ which is sign-compatible with the stoichiometric subspace $S$. In other words, we check if there is a linear combination of a basis for $S$ that has the same sign with $\mu$.

At this point, we have the following system.

$-0.2344{\mu_{X_{2}}}=M_{1}=0.7464{\mu_{X_{1}}}+0.0243{\mu_{X_{5}}}$

$0.7318{\mu_{X_{2}}}-0.3941{\mu_{X_{5}}}=M_{7}$

$8.6107{\mu_{X_{2}}}=M_{5}$

$0.6159{\mu_{X_{3}}}+0.1308{\mu_{X_{5}}}=M_{7}=0.05{\mu_{X_{3}}}+0.533{\mu_{X_{4}}}-0.0822{\mu_{X_{5}}}$

${\mu_{X_{5}}}=M_{11}$

Ignoring the last equation, we obtain the following values.

$\mu_{X_{1}}=(3442493125M_{1})/2154694458+(101250M_{7})/1225651$

$\mu_{X_{2}}=-(1250M_{1})/293$

$\mu_{X_{3}}=(3988310000M_{1})/2370625789+(52490000M_{7})/24272619$

$\mu_{X_{4}}=(16585421000M_{7})/12937305927-(134085884500M_{1})/97195657349$

$\mu_{X_{5}}=-(9147500M_{1})/1154713-(10000M_{7})/3941$

By considering any of the cases in STEP 13, the system yields a contradiction. Since the system is inconsistent, there is no need to put additional nonlinear equations from STEP 7. Therefore, we cannot find a signature $\mu$ and the system does not have the capacity to admit multiple equilibria. We exit the algorithm.

Remark 3.1**.**

In ERM0-G, not only the network properties are essential for the algorithm but we are also particular with the kinetic order values of the kinetic system. From STEPS 1 to 9, we use the structure of the reaction network alone. Within these steps, the orientation is partitioned into equivalence classes. With respect to this partitioning, consequently, the reaction set is partitioned into fundamental classes. Until the assignment of sign patterns for the fundamental classes and the shelving assignment of the reactions in these classes, there is no use of the kinetic order values in the $T$-matrix. The use of these values takes place from STEPS 10 to 14 when we form the system of equations/inequalities until we solve such system.

3.2 Application to the Heck et al.’s Global Carbon Cycle Model

The global carbon cycle model of Heck et al. [14] is built from the well-cited box model of Anderies et al. [1]. The latter was calibrated according to empirically observed and simulated carbon cycle dynamics. The revised model also incorporates a societal intervention called terrestrial carbon dioxide removal (tCDR). In tCDR, terrestrial carbon is deliberately recovered and dumped into a new sink, referred as carbon engineering sink. A detailed discussion of the model calibration and assumptions can be found in [14].

In the CRN-based analysis of the power-law approximation of this Earth system, Fortun et al. [11] considered pooling the geological carbon pool and the new sink to form a passive carbon pool. Moreover, they separated the atmospheric carbon into two nodes: atmospheric carbon in the pre-industrial state and additional atmospheric carbon due to fossil fuel use. A biochemical map of the system is shown in Figure 2 [11].

The list of network numbers is provided in Table 2. Because of the high deficiency of the network, it becomes a good candidate for the application of MSA to determine the capacity of multiple positive steady states of the system. Multistationarity in global climate implies that there may exist “tipping points” beyond which a return to the original state may be difficult or prolonged. By showing that multistationarity in a global climate model may exist, therefore, makes the search for the tipping point relevant. After the tipping point is determined, appropriate measures may then be set to avoid exceeding it.

The reaction network of the Heck et al.’s global carbon cycle model is given below [11].

[TABLE]

The kinetic system has reactant-determined kinetics and the $T$-matrix is given below.

\bordermatrix{&A_{1}+2A_{2}&A_{1}+A_{2}&A_{2}&A_{3}&A_{4}+A_{5}&A_{1}+2A_{4}&A_{1}+A_{4}&A_{4}&A_{1}+A_{2}+A_{4}\cr A_{1}&199.75&159.84&0&0&0&-43.80&-56.13&0&1\cr A_{2}&-86.03&-63.32&1&0&0&0&0&0&4.44\cr A_{3}&0&0&0&1&0&0&0&0&0\cr A_{4}&0&0&0&0&1&21.42&22.19&1&11.52\cr A_{5}&0&0&0&0&1.54&0&0&0&0\cr}

Choose ${\mathscr{O}}=\{R_{1},R_{2},R_{4},R_{5},R_{6},R_{7},R_{8},R_{10}\}.$ Below is basis for $Ker{L_{\mathscr{O}}}$ obtained by solving $\sum\limits_{y\to y^{\prime}\in\mathscr{O}}{{\alpha_{y\to y^{\prime}}}\left({y^{\prime}-y}\right)}=0$.

[TABLE]

A basis for $Ker^{\perp}L_{\mathscr{O}}\cap\Gamma_{W}$ is also given below.

[TABLE]

We can verify that our example has a forest basis for $Ker^{\perp}L_{\mathscr{O}}\cap\Gamma_{W}$. We will not present the details here. Without this fact, one can obtain a system of equations/inequalities up to STEP 14 and try if $\kappa$ can be solved after getting the value of a solution $\mu$. We pick the sign patterns to be positive and obtain the following shelving assignment. If any of ${{\cal U}_{i}}$, ${{\cal M}_{i}}$, or ${{\cal L}_{i}}$ does not appear in the list, then it is empty.

${{\cal M}_{1}}=\left\{{R_{1}:A_{1}+2A_{2}\to 2A_{1}+A_{2}}\right\}$

${{\cal M}_{2}}=\left\{{R_{2}:A_{1}+A_{2}\to 2A_{2}}\right\}$

${{\cal L}_{3}}\ =\left\{{R_{3}:A_{2}\to A_{3},R_{4}:A_{3}\to A_{2},R_{8}:A_{4}\to A_{3},R_{9}:A_{3}\to A_{4}}\right\}$

${{\cal M}_{4}}=\left\{{R_{5}:A_{4}+A_{5}\to 2A_{4},R_{10}:A_{1}+A_{2}+A_{4}\to A_{5}+A_{2}+A_{4}}\right\}$

${{\cal M}_{5}}=\left\{{R_{6}:A_{1}+2A_{4}\to 2A_{1}+A_{4}}\right\}$

${{\cal M}_{6}}=\left\{{R_{7}:A_{1}+A_{4}\to 2A_{4}}\right\}$

We choose $M_{2}<M_{1}<M_{3}$ from $b^{1}$ and $M_{5}<M_{6}<M_{3}$ from $b^{2}$. Thus, we obtain the following system.

[TABLE]

Since the system of inequalities has a nonzero solution $\mu$, the kinetic system has the capacity for multistationarity. A solution of the system together with the multiple steady states are given in Table 3. The computation of a particular set of rate constants is given in Appendix D.1.

3.3 A Comparison with the Deficiency One Algorithm for PL-RDK systems

We compare the MSA and the DOA on regular networks. Specifically, we apply the MSA on the main example of [9], the Anderies et al.’s pre-industrial carbon cycle model in [1]. Indeed, we can verify from Appendix D.2 that this deficiency-one model has the capacity for multistationarity for some rate constants [9].

Definition 3.2**.**

A pair of complexes $\{i,j\}$ form a cut pair if they are adjacent and the removal of the reaction arrow(s) between these complexes results in a separation of the linkage class containing them.

Let $\mathscr{C}^{1}$, $\mathscr{C}^{2}$, … , $\mathscr{C}^{t}$ denote the complex sets of the terminal strong linkage classes and let $\mathscr{C}^{\prime}=\bigcup_{k=1}^{t}\mathscr{C}^{i}$.

Definition 3.3**.**

A reaction network $\mathscr{N}=\left(\mathscr{S},\mathscr{C},\mathscr{R}\right)$ is said to be regular if the following conditions are satisfied.

• i.

The reaction vectors are positively dependent, i.e., there exists a set of positive numbers $\alpha_{ij}$ for all $(i,j)\in\mathscr{R}$ such that $\sum\limits_{(i,j)\in\mathscr{R}}{{\alpha_{ij}}\left({j-i}\right)}=0$.

• ii.

$\mathscr{N}$* is t-minimal.*

• iii.

The complexes $i$ and $j$ form a cut pair $\forall\ i,j\in\mathscr{C}^{\prime}$ such that $(i,j)\in\mathscr{R}$.

Positive dependency is a necessary condition for the existence of a positive equilibrium. In this section, we present an example where MSA can obtain reasonable results for deficiency one, non-$t$-minimal networks and networks without the “cut pair” property, which are both outside the scope of the DOA for PL-RDK.

Example 3.4**.**

Consider the following reaction network and its kinetic order matrix.

[TABLE]

The network has 2 linkage classes and 3 terminal strong linkage classes, and hence, not $t$-minimal. Its deficiency is 1, its rank is 3 and it is not weakly reversible. Moreover, the system is PL-RDK. Also, the complexes $A_{1}+A_{2}$ and $2A_{3}$ do not form a cut pair. Since neither of last two conditions for regularity of a reaction network is satisfied, the DOA for PL-RDK is not applicable. We will also establish that this system has the capacity to admit multiple steady states using the MSA.

A basis for $KerL_{\mathscr{O}}$ is

[TABLE]

*which induces the following partition:

${P_{1}}=\left\{{{R_{1}}}\right\}={C_{1}}$, ${P_{2}}=\left\{{{R_{2}},{R_{3}}}\right\}={C_{2}}$, and ${P_{3}}=\left\{{{R_{4}},{R_{5}}}\right\}={C_{3}}$.

Each reaction is irreversible so we have the following shelving assignment: ${{\cal M}_{1}}={C_{1}}$, ${{\cal M}_{2}}={C_{2}}$, and ${{\cal M}_{3}}={C_{3}}$. Thus, we have the equations: ${\mu_{{A_{1}}}}+{\mu_{{A_{3}}}}={M_{1}}$ and ${\mu_{{A_{1}}}}+{\mu_{{A_{2}}}}={M_{2}}={M_{3}}={\mu_{{A_{3}}}}$. From the basis $\{(1,-1,1)\}$ for $Ker^{\perp}L_{\mathscr{O}}\cap\Gamma_{W}$, we have $M_{1}=M_{2}=M_{3}$ since we already established ${M_{2}}={M_{3}}$. From these equations, we get $\mu=(1,0,1)$ which is stoichiometrically compatible with $\sigma=(3,0,2)\in S$. (You can choose other $\sigma$ which is stoichiometrically compatible with $\mu$, say $\sigma=\mu=(1,0,1)$.) Thus, the kinetic system has the capacity for multistationarity. We choose $\kappa=(1,2,2,1,1)$, a linear combination of the given basis for $KerL_{\mathscr{O}}$.*

The equilibria $c^{**}$ and $c^{*}$ are given in Table 4. We refer to Table 5 and check that $\sum\limits_{y\to y^{\prime}\in{\mathscr{R}}}{{\kappa_{y\to y^{\prime}}}\left({y^{\prime}-y}\right)=0}$ and $\sum\limits_{y\to y^{\prime}\in{\mathscr{R}}}{{\kappa_{y\to y^{\prime}}}{e^{{T_{.y}\cdot\mu}}}\left({y^{\prime}-y}\right)=0}$ are both satisfied since $\kappa_{y\to y^{\prime}}e^{T_{.y}\cdot\mu}=k_{y\to y^{\prime}}c^{*{T_{.y}}}$ is a multiple of the vector $\kappa$.

4 Multistationarity Algorithm for Power-Law Non-Reactant-Determined Kinetics

Majority of the embedded GMA systems of the 15 identified BST models in [2] are PL-NDK. Hence, there is a need of transforming these PL-NDK systems to dynamically equivalent PL-RDK systems. As a result, we combine the CF-RM method and the extension of the HDA to solve the problem of multistationarity of such systems.

4.1 The CF-RM Method

In this subsection, we present the CF-RM transformation method (transformation of complex factorizable kinetics by reactant multiples). Further discussion can be found in [18]. From a given PL-NDK system, one can construct a PL-RDK system using this method. In this process, at each reactant complex, the branching reactions are partitioned into CF-subsets (a CF-subset contains reactions having the same kinetic order vectors). If a reactant complex has more than one CF-subset, then it is an NF-reactant complex (which makes the system NDK). For each subset, a complex is added to both the reactant and the product complexes of a reaction. This leaves the reaction vectors unchanged. Hence, the stoichiometric subspace remains the same which guarantees the dynamic equivalence of the original and its transform. Under this method, the kinetic order matrix does not change as well.

Let $\mathscr{N}=\left(\mathscr{S},\mathscr{C},\mathscr{R}\right)$ be a reaction network and $\rho:\mathscr{R}\to\mathscr{C}$ be the reactant map. If $y\in\mathscr{C}$, then $\rho^{-1}(y)$ is its reaction set. Let $\iota:\mathscr{R}\to\mathbb{R}^{\mathscr{S}}$ be the interaction map of the system. This map assigns to each reaction its kinetic order row in the kinetic order matrix. If $x\in\iota\left(\rho^{-1}(y)\right)$ then $\{r\in\rho^{-1}(y)|\iota(r)=x\}$ is called a CF-subset of $y$. We denote the number of CF-subsets of $y$ by $N_{R}(y)$. There are $|\iota\left(\rho^{-1}(y)\right)|$ such subsets. Note that a reactant complex is a CF-reactant complex if and only if $N_{R}(y)=1$.

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

Determine the reactant set $\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. If there are several such subsets, then just choose one. Leave this CF-subset unchanged.

• 4.

At this point, there are $N_{R}(y)-1$ remaining CF-subset. For each of these subsets, choose successively a multiple of $y$ which is not among the current set of reactants.

Since we leave each CF-reactant complex unchanged, no new reactant is introduced at a CF-reactant complex. Hence, the total number of new reactants is given by $\sum{\left({{N_{R}}\left(y\right)-1}\right)}$ with the sum taken over all reactants.

We provide an illustration of our computational approach. We show how it can determine multistationarity of a deficiency one PL-RDK.

4.2 A Deficiency One PL-NDK System

In this section, we present a weakly reversible PL-NDK system with deficiency one and we apply the combined CF-RM and HDA approach. Typically, the deficiency will increase. Hence, only the MSA can handle the resulting PL-RDK system. This example further emphasizes the advantages of the MSA over the DOA for deficiency one systems.

Example 4.1**.**

Consider the following reaction network.

[TABLE]

The network has 1 linkage class and 1 strong linkage class which is terminal. Its deficiency is 1 and its rank is 1. Moreover, it is weakly reversible. Consider the following kinetic order values: for $R_{1}:0$, for $R_{2}:0.5$, for $R_{3}:1$, and for $R_{4}:0.5$. Hence, the system is PL-NDK. Using the CF-RM transformation, we modify $R_{3}:3A_{1}\to 4A_{1}$. The deficiency of the CF-transform is 2. We set $\mathscr{O}=\{R_{1},R_{3},R_{4}\}$ and the following are bases for $KerL_{\mathscr{O}}$ and $KerL_{\mathscr{R}}$ are

[TABLE]

The basis for $KerL_{\mathscr{O}}$ given above induces the following partition: ${P_{1}}=\left\{{{R_{1}}}\right\}\subset{C_{1}}=\left\{{{R_{1},R_{2}}}\right\}$, ${P_{2}}=\left\{{{R_{3}}}\right\}={C_{2}}$, and ${P_{3}}=\left\{{{R_{4}}}\right\}={C_{3}}$. In $C_{1}$, consider sign pattern $(+,0)$ for $\left(g_{W},h_{W}\right)$ and the rest are positive since they are nonreversible. Then, we have the following shelving assignment: ${{\cal U}_{1}}={C_{1}}$, ${{\cal M}_{2}}={C_{2}}$, and ${{\cal M}_{3}}={C_{3}}$ which gives $M_{1}<0{\mu_{{A_{1}}}}<0.5{\mu_{{A_{1}}}}$, $1{\mu_{{A_{1}}}}={M_{2}}$, and $0.5{\mu_{{A_{1}}}}={M_{3}}$. Now, a basis for for $Ker^{\perp}L_{\mathscr{O}}\cap\Gamma_{W}$ is $\{(-1,-1,2)\}$ and we choose $M_{1}<M_{3}<M_{2}$. Moreover, take $\kappa=(2,1,1,1)$, a linear combination of the basis elements for $KerL_{\mathscr{R}}$ so $\sum\limits_{y\to y^{\prime}\in{\mathscr{R}}}{{\kappa_{y\to y^{\prime}}}\left({y^{\prime}-y}\right)=0}$. From the sign pattern of $C_{1}$ and we recall from our theory (for reversible reactions) that ${g_{0\to{A_{1}}}}={\kappa_{0\to{A_{1}}}}-{\kappa_{{A_{1}}\to 0}}=2-1>0$ and ${h_{0\to{A_{1}}}}={\kappa_{0\to{A_{1}}}}{e^{0}}-{\kappa_{{A_{1}}\to 0}}{e^{0.5{\mu_{{A_{1}}}}}}=0$. Hence, $2{e^{0}}=1{e^{0.5{\mu_{{A_{1}}}}}}$ and ${\mu_{{A_{1}}}}=\ln 4$ which satisfies the obtained system. Choose $\sigma=3\in S$ so $c^{**}=1$ and $c^{*}=4$. We can verify that

[TABLE]

are both the zero vector.

5 Conclusions and Outlook

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

We modified the higher deficiency algorithm for mass action kinetics of Ji and Feinberg to handle the problem of determining whether a PL-RDK system has the capacity to admit multiple equilibria. This was done by replacing the role of the molecularity matrix by the $T$-matrix.

• 2.

We applied the algorithm to the embedded network of the GMA model of anaerobic fermentation pathway of Saccharomyces cerevisiae and determined that it has no capacity to admit multiple equilibria, no matter what positive rate constants are assumed. On the other hand, we determined that the Heck et al.’s global carbon cycle model has the capacity for multistationarity for particular rate constants.

• 3.

We compared the MSA and the DOA on regular networks by applying it on the Anderies et al.’s pre-industrial carbon cycle model. Indeed, the deficiency-one model has the capacity for multistationarity for some rate constants. This is consistent with the results given in [9]. We also presented an example where MSA can obtain results for deficiency one, non-$t$-minimal networks and networks without the “cut pair” property since these are outside the scope of the DOA.

• 4.

We used the CF-RM transformation to convert a PL-NDK system to a dynamically equivalent PL-RDK system. We also provided a weakly reversible PL-NDK system with deficiency one and we applied the MSA. Under the CF-RM, the deficiency of a network typically increases. Hence, only the MSA can handle the resulting PL-RDK system which further emphasizes its advantage.

• 5.

We can also look into the nonlinear case where there is no forest basis and find an example that would be able to apply the algorithm. In addition, one can explore a possible extension of the algorithm to Rate constant-Interaction map-Decomposable (RID) kinetics with interaction parameter maps.

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

Appendix B Fundamentals of Chemical Reaction Networks and Kinetic Systems

In this section, we present some fundamentals of chemical reaction networks and chemical kinetic systems. These concepts are provided in [6] and [15].

Definition B.1**.**

A chemical reaction network $\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}$.

We can view $\mathscr{C}$ as a subset of $\mathbb{R}^{\mathscr{S}}_{\geq 0}$. The ordered pair $\left({{C_{i}},{C_{j}}}\right)$ corresponds to the reaction ${C_{i}}\to{C_{j}}$.

Definition B.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}$.

Definition B.3**.**

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 B.4**.**

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

[TABLE]

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 B.5**.**

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

Definition B.6**.**

A vector $x\in{\mathbb{R}^{\mathscr{S}}}$ is stoichiometrically compatible with the stoichiometric subspace $S$ if $sign\left({{x_{s}}}\right)=sign\left({{\sigma_{s}}}\right)\forall s\in\mathscr{S}$ for some $\sigma\in S$.

We can view complexes as vertices and reactions as edges. With this, chemical reaction networks 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 chemical reaction network. 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, 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$, and it is said to be $t$-minimal if $t=l$.

Definition B.7**.**

For a chemical reaction network, 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$.

Definition B.8**.**

A kinetics $K$ for a reaction network $\left(\mathscr{S},\mathscr{C},\mathscr{R}\right)$ is an assignment to each reaction $j:y\to y^{\prime}\in\mathscr{R}$ of a rate function ${K_{j}}:{\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_{j}}\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_{j}(c)>0$. The system $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ is called a chemical kinetic system.

Definition B.9**.**

The species formation rate function (SFRF) of a chemical kinetic system is given by

[TABLE]

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$.

Definition B.10**.**

The set of positive equilibria of a chemical kinetic system $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ is given by

[TABLE]

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.

Definition B.11**.**

A kinetics $K$ is a power-law kinetics (PLK) if ${K_{i}}\left(x\right)={k_{i}}{{x^{{F_{i}}}}}$ $\forall 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 B.12**.**

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. That is, ${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.

Definition B.13**.**

[17]*

The $m\times n$ matrix $\widetilde{Y}$ is given by*

[TABLE]

Definition B.14**.**

The $T$-matrix is the $m\times n_{r}$ truncated $\widetilde{Y}$ matrix where the nonreactant columns are removed.

Appendix C Proofs of the Lemmas and the Main Theorem

of Lemma 2.1.

We assume positive equilibria ${c}^{*}$ and ${c}^{**}$ are distinct so ${\mu_{s}}\neq 0$ for some $s\in\mathscr{S}$. Thus, $\mu$ is nonzero. By the monotonicity of $\ln$, $\ln\left({{c_{s}}^{*}}\right)>\ln\left({{c_{s}}^{**}}\right)$ whenever ${c_{s}}^{*}>{c_{s}}^{**}$. Hence, ${c_{s}}^{*}-{c_{s}}^{**}$ and $\ln\left({{c_{s}}^{*}}\right)-\ln\left({{c_{s}}^{*}}\right)$ have the same signs. Now, ${\mu_{s}}=\ln\left({\dfrac{{{c_{s}}^{*}}}{{{c_{s}}^{**}}}}\right)=\ln\left({{c_{s}}^{*}}\right)-\ln\left({{c_{s}}^{**}}\right)$. Since ${c}^{*}$ and ${c}^{**}$ are stoichiometrically compatible, by definition ${c^{*}}-{c^{**}}\in S$. Then, $\mu$ is stoichiometrically compatible with $S$. Finally, take ${\kappa_{y\to y^{\prime}}}={k_{y\to y^{\prime}}}{c}{{}^{**{T_{.y}}}}\ \forall\ y\to y^{\prime}\in\mathscr{R}$ from the set of positive rate constants $\left\{{{k_{y\to y^{\prime}}}|y^{\prime}\to y\in\mathscr{R}}\right\}$. ∎

of Lemma 2.2.

Since $\mu\neq 0$ then $\mu_{s}\neq 0$ for some $s\in\mathscr{S}$. Hence, ${\dfrac{{{c_{s}}^{*}}}{{{c_{s}}^{**}}}={e^{{\mu_{s}}}}\neq 1}$ so ${c}^{*}$ and ${c}^{**}$ are distinct. Note that ${c_{s}}^{*}$ and ${c_{s}}^{**}$ must have the same sign and we set both of them to be positive. Let $\sigma\in S$ be stoichiometrically compatible with $\mu$. If $\mu_{s}\neq 0$, then we set $c_{s}^{*}=\dfrac{{{\sigma_{s}}{e^{{\mu_{s}}}}}}{{{e^{{\mu_{s}}}}-1}}$ and $c_{s}^{**}=\dfrac{{{\sigma_{s}}}}{{{e^{{\mu_{s}}}}-1}}$ for all $s\in\mathscr{S}$. In this case, ${c_{s}}^{*}-{c_{s}}^{**}={\sigma_{s}}\left({\dfrac{{{e^{{\mu_{s}}}}-1}}{{{e^{{\mu_{s}}}}-1}}}\right)={\sigma_{s}}$. If $\mu_{s}=0$, then we set $c_{s}^{*}=c_{s}^{**}=p$ for some positive number $p$. In this case, ${c_{s}}^{*}-{c_{s}}^{**}=p-p=0$. Hence, ${c}^{*}$ and ${c}^{**}$ are stoichiometrically compatible. Finally, take ${k_{y\to y^{\prime}}}=\dfrac{{{\kappa_{y\to y^{\prime}}}}}{{{c}{{{}^{**}}^{T_{.y}}}}}.$ ∎

of Lemma 2.9.

Assume $y\to y^{\prime}$ is irreversible. Then, ${g_{y\to y^{\prime}}}={\kappa_{y\to y^{\prime}}}>0$, ${h_{y\to y^{\prime}}}={\kappa_{y\to y^{\prime}}}{e^{{T_{.y}}\cdot\mu}}>0$, and ${\rho_{y\to y^{\prime}}}=\dfrac{{{h_{y\to y^{\prime}}}}}{{{g_{y\to y^{\prime}}}}}=\dfrac{{{\kappa_{y\to y^{\prime}}}{e^{{T_{.y}}\cdot\mu}}}}{{{\kappa_{y\to y^{\prime}}}}}={e^{{T_{.y}}\cdot\mu}}$. Assume $y\to y^{\prime}$ is reversible. Suppose ${g_{y\to y^{\prime}}}\neq 0$. Now, ${\rho_{y\to y^{\prime}}}{g_{y\to y^{\prime}}}={h_{y\to y^{\prime}}}$, ${g_{y\to y^{\prime}}}={\kappa_{y\to y^{\prime}}}-{\kappa_{y^{\prime}\to y}}$ and ${h_{y\to y^{\prime}}}={\kappa_{y\to y^{\prime}}}{e^{{T_{.y}}\cdot\mu}}-{\kappa_{y^{\prime}\to y}}{e^{{T_{.y^{\prime}}}\cdot\mu}}$. Then, ${g_{y\to y^{\prime}}}{e^{{T_{.y}}\cdot\mu}}={\kappa_{y\to y^{\prime}}}{e^{{T_{.y}}\cdot\mu}}-{\kappa_{y^{\prime}\to y}}{e^{{T_{.y}}\cdot\mu}}$. Also, ${\kappa_{y^{\prime}\to y}}{e^{{T_{.y}}\cdot\mu}}-{\kappa_{y^{\prime}\to y}}{e^{{T_{.y^{\prime}}}\cdot\mu}}={\rho_{y\to y^{\prime}}}{g_{y\to y^{\prime}}}-{g_{y\to y^{\prime}}}{e^{{T_{.y}}\cdot\mu}}$ and ${\kappa_{y^{\prime}\to y}}\left({{e^{{T_{.y}}\cdot\mu}}-{e^{{T_{.y^{\prime}}}\cdot\mu}}}\right)={g_{y\to y^{\prime}}}\left({{\rho_{y\to y^{\prime}}}-{e^{{T_{.y}}\cdot\mu}}}\right)$. Hence,

${\kappa_{y^{\prime}\to y}}=\dfrac{{{\rho_{y\to y^{\prime}}}-{e^{{T_{.y}}\cdot\mu}}}}{{{e^{{T_{.y}}\cdot\mu}}-{e^{{T_{.y^{\prime}}}\cdot\mu}}}}{g_{y\to y^{\prime}}}.$ Similarly, we can solve for the value of ${\kappa_{y\to y^{\prime}}}$ with ${\kappa_{y\to y^{\prime}}}=\dfrac{{{\rho_{y\to y^{\prime}}}-{e^{{T_{.y^{\prime}}}\cdot\mu}}}}{{{e^{{T_{.y}}\cdot\mu}}-{e^{{T_{.y^{\prime}}}\cdot\mu}}}}{g_{y\to y^{\prime}}}.$ Thus, we obtain the following.

• i.

If $y\to y^{\prime}$ is irreversible then ${g_{y\to y^{\prime}}}>0$, ${h_{y\to y^{\prime}}}>0$, and ${\rho_{y\to y^{\prime}}}={e^{{T_{.y}}\cdot\mu}}$.

• ii.

Suppose $y\to y^{\prime}$ is reversible.

• a.

If ${g_{y\to y^{\prime}}}>0$, then either ${\rho_{y\to y^{\prime}}}={e^{{T_{.y}}\cdot\mu}}={e^{{T_{.y^{\prime}}}\cdot\mu}}$, ${\rho_{y\to y^{\prime}}}>{e^{{T_{.y}}\cdot\mu}}>{e^{{T_{.y^{\prime}}}\cdot\mu}}$, or ${\rho_{y\to y^{\prime}}}<{e^{{T_{.y}}\cdot\mu}}<{e^{{T_{.y^{\prime}}}\cdot\mu}}$.

• b.

If ${g_{y\to y^{\prime}}}<0$, then either ${\rho_{y\to y^{\prime}}}={e^{{T_{.y^{\prime}}}\cdot\mu}}={e^{{T_{.y}}\cdot\mu}}$, ${\rho_{y\to y^{\prime}}}>{e^{{T_{.y^{\prime}}}\cdot\mu}}>{e^{{T_{.y}}\cdot\mu}}$, or ${\rho_{y\to y^{\prime}}}<{e^{{T_{.y^{\prime}}}\cdot\mu}}<{e^{{T_{.y}}\cdot\mu}}$.

• c.

If ${g_{y\to y^{\prime}}}=0$ and ${h_{y\to y^{\prime}}}>0$ then ${e^{{T_{.y}}\cdot\mu}}>{e^{{T_{.y^{\prime}}}\cdot\mu}}$.

• d.

If ${g_{y\to y^{\prime}}}=0$ and ${h_{y\to y^{\prime}}}<0$ then ${e^{{T_{.y}}\cdot\mu}}<{e^{{T_{.y^{\prime}}}\cdot\mu}}$.

• e.

If ${g_{y\to y^{\prime}}}={h_{y\to y^{\prime}}}=0$ then ${e^{{T_{.y}}\cdot\mu}}={e^{{T_{.y^{\prime}}}\cdot\mu}}$.

Now, let $i\in\{1,2,...,w\}$ and $y\to y^{\prime}\in P_{i}$. Also, let ${y_{i}\to y_{i}^{\prime}}$ be the representative of $P_{i}$. Then ${\rho_{y_{i}\to y_{i}^{\prime}}}=\dfrac{{{h_{{y_{i}}\to{y_{i}}^{\prime}}}}}{{{g_{{y_{i}}\to{y_{i}}^{\prime}}}}}=\dfrac{{\alpha_{y_{i}\to y_{i}^{\prime}}}{{h_{{y_{i}}\to{y_{i}}^{\prime}}}}}{{\alpha_{y_{i}\to y_{i}^{\prime}}}{{g_{{y_{i}}\to{y_{i}}^{\prime}}}}}=\dfrac{{{h_{{y}\to{y}^{\prime}}}}}{{{g_{{y}\to{y}^{\prime}}}}}={\rho_{y\to y^{\prime}}}$ for some nonzero ${\alpha_{y_{i}\to y_{i}^{\prime}}}$. ∎

of Lemma 2.10.

By considering the columns of the $T$-matrix instead of the columns of the molecularity matrix, we obtain an analogous proof to to one given in [15]. ∎

Lemma C.1**.**

[15*]** Suppose a reaction network satisfies the following properties for an orientation $\mathscr{O}$: $P_{i}$ ($i=0,1,2,...,w$) is defined by a representative ${{y_{i}}\to{y_{i}}^{\prime}}$, $W=\{{{y_{i}}\to{y_{i}}^{\prime}}|i=1,2,...,w\}\subseteq\mathscr{O}$, a given basis $\left\{{{b^{j}}}\right\}_{j=1}^{q}$ for $Ker^{\perp}{L_{\mathscr{O}}}\cap{\Gamma_{W}}$ such that its basis graph is a forest, a valid pair of sign patterns for ${h_{W}},{g_{W}}$, and a set of parameters

$\left\{{{\rho_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)|{g_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)\neq 0,i=1,2,...,w}\right\}$ where the sign of ${{\rho_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}$ is the same as the ratio of the signs of ${{h_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}$ and ${{g_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)}$. Then ${h_{W}},{g_{W}}\in{\mathbb{R}^{\mathscr{O}}}\cap{\Gamma_{W}}$ satisfy Equations (1), (2), and (3) if and only if ${\rho_{W}}\left({{y_{i}}\to{y_{i}}^{\prime}}\right)$’s satisfy the conditions in Lemma 2.12.*

of Theorem 2.14.

By Lemma 2.1, (i) is immediate. Since positive and distinct equilibria exist, $KerL_{\mathscr{O}}$ is nontrivial. Define $g,h\in\mathbb{R}^{\mathscr{O}}$ such that

{g_{y\to y^{\prime}}}=\left\{\begin{array}[]{ll}{\kappa_{y\to y^{\prime}}}-{\kappa_{y^{\prime}\to y}}&{\rm{if}}\ y\to y^{\prime}\in\mathscr{O}{\rm{\ is\ reversible}}\\ {\kappa_{y\to y^{\prime}}}&{\rm{if}}\ y\to y^{\prime}\in\mathscr{O}{\rm{\ is\ irreversible}}\end{array}\right. and

{h_{y\to y^{\prime}}}=\left\{\begin{array}[]{ll}{\kappa_{y\to y^{\prime}}}{e^{{T_{.y}}\cdot\mu}}-{\kappa_{y^{\prime}\to y}}{e^{{T_{.y^{\prime}}}\cdot\mu}}&{\rm{if}}\ y\to y^{\prime}\in\mathscr{O}{\rm{\ is\ reversible}}\\ {\kappa_{y\to y^{\prime}}}{e^{{T_{.y}}\cdot\mu}}&{\rm{if}}\ y\to y^{\prime}\in\mathscr{O}{\rm{\ is\ irreversible}}\end{array}\right..

For $i=1,2,...,w$, let ${P_{i}}$ be the equivalence class with ${y_{i}}\to{y_{i}}^{\prime}$ as representative and $M_{i}$ is associated with the nondegenerate fundamental class. Let $\{b^{j}\}_{j=1}^{q}$ be a basis for $Ker^{\perp}{L_{\mathscr{O}}}\cap{\Gamma_{W}}$, if it exists. Also, let $g,h\in Ker{L_{\mathscr{O}}}$. Hence, $b^{j}\cdot g=\sum\limits_{{y_{i}}\to{y_{i}}^{\prime}\in W}{g_{{y_{i}}\to{y_{i}}^{\prime}}b_{{y_{i}}\to{y_{i}}^{\prime}}^{j}=0}$ and $b^{j}\cdot h=\sum\limits_{{y_{i}}\to{y_{i}}^{\prime}\in W}{h_{{y_{i}}\to{y_{i}}^{\prime}}b_{{y_{i}}\to{y_{i}}^{\prime}}^{j}=0}$ for $j=1,2,...,q$. These two equations can be written as Equations (1), (2), and (3). Hence, (iii) which leads to (iv).

On the other hand, by Lemma 2.2, we only need to solve show the existence of $\kappa$. By Lemma C.1 which requires the existence of the forest basis, Equations (1), (2), and (3) hold. From here, we obtain the values of $g_{y\to y^{\prime}}$ and $h_{y\to y^{\prime}}$. These yield $\kappa$. ∎

Appendix D Computation Details for Section 3

D.1 Computation of a Particular Set of Rate Constants for the Heck et al.’s Carbon Cycle Model

Table 6 is a summary of values for the Heck et al.’s terrestrial carbon recovery model. The rate constant vector $k_{y\to y^{\prime}}=\dfrac{\kappa_{y\to y^{\prime}}}{c^{{**}{T_{.y}}}}$ is indicated. Note that $\kappa_{y\to y^{\prime}}$ was solved by getting the scalars in the linear combination of a basis for $KerL_{\mathscr{R}}$ such that

[TABLE]

D.2 Application of the MSA on the Anderies et al.’s Pre-industrial Carbon Cycle

The reaction network and the $T$-matrix of the Anderies et al.’s pre-industrial carbon cycle in [1] are given below.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We choose the forward reactions so ${\mathscr{O}}=\left\{{{A_{1}}+2{A_{2}}\to 2{A_{1}}+{A_{2}},{A_{1}}+{A_{2}}\to 2{A_{2}},{A_{2}}\to{A_{3}}}\right\}.$ A basis for $Ker{L_{\mathscr{O}}}$, and the corresponding equivalence and fundamental classes are given below.

[TABLE]

${P_{0}}=\left\{{{A_{2}}\to{A_{3}}}\right\}\subset{C_{0}}=\left\{{{A_{2}}\to{A_{3}},{A_{3}}\to{A_{2}}}\right\}$

${P_{1}}=\left\{{{A_{1}}+2{A_{2}}\to 2{A_{1}}+{A_{2}},{A_{1}}+{A_{2}}\to 2{A_{2}}}\right\}={C_{1}}$

We choose $W=\left\{{{A_{1}}+{A_{2}}\to 2{A_{2}}}\right\}$ with $w=1$. The resulting inequality system is automatically linear since $Ker^{\perp}L_{\mathscr{O}}\cap\Gamma_{W}$ is trivial. Since $P_{1}$ is nonreversible, the sign pattern for ${g_{W}}$ and ${h_{W}}$ for ${y_{1}}\to{y_{1}}^{\prime}$ must be positive. We have the following shelving assignment: ${{\cal U}_{1}}=\left\{\right\},{{\cal M}_{1}}=\left\{{{A_{1}}+2{A_{2}}\to 2{A_{1}}+{A_{2}},{A_{1}}+{A_{2}}\to 2{A_{2}}}\right\},{{\cal L}_{1}}=\left\{{}\right\}$. From the middle shelf, we obtain the equation ${p_{1}}{\mu_{{A_{1}}}}+{q_{1}}{\mu_{{A_{2}}}}={M_{1}}={p_{2}}{\mu_{{A_{1}}}}+{q_{2}}{\mu_{{A_{2}}}}.$ From $P_{0}$, we add ${\mu_{{A_{2}}}}={\mu_{{A_{3}}}}$ to the system. Since the $Ker^{\perp}L_{\mathscr{O}}\cap\Gamma_{W}$ is trivial, no additional equation or inequality can be obtained here. We have the following system.

[TABLE]

Hence, ${\mu_{{A_{1}}}}=\dfrac{{{q_{2}}-{q_{1}}}}{{{p_{1}}-{p_{2}}}}{\mu_{{A_{2}}}}$ and $\mu=\left(\dfrac{{{q_{2}}-{q_{1}}}}{{{p_{1}}-{p_{2}}}}{\mu_{{A_{2}}}},{\mu_{{A_{2}}}},{\mu_{{A_{2}}}}\right)$ with ${\mu_{{A_{2}}}}>0$. We take $\left({-2,1,1}\right)$ so that $\mu$ is stoichiometrically compatible with $\sigma\in S$. Thus, $\mu$ is a signature and the reaction network has the capacity to admit multiple equilibria. This is consistent with the results of Fortun et al. when they applied the deficiency-one algorithm for PL-RDK in [9].

D.3 A Remark on the Steps of the MSA

Remark D.1**.**

In the case where the first attempt in adding $M$ equations/inequalities has no solution, we introduce STEP 15 where STEPS 13 and 14 are repeated for every choice of $M$ equations/inequalities. A choice is picking one subcase in STEP 13. If there is still no solution, one proceeds with STEP 16 where STEPS 9 to 15 are repeated for possible shelving assignments in STEP 9. If this happens again, we go to STEP 17 where STEPS 8 to 16 are repeated. All sign pattern choices for $g_{W},h_{W}\in\mathbb{R}^{\mathscr{O}}\cap\Gamma_{W}$ in STEP 8 are repeated in this step. If there is no signature or pre-signature was found, then the system does not have the capacity to admit multiple stady states, no matter what positive values of the rate constants we assume.

Appendix E Additional Useful Propositions

The following propositions are useful in determining if a chemical reaction network does not have the capacity to admit multiple steady states if the reaction has at least one irreversible inflow reaction or irreversible outflow reaction.

Proposition E.1**.**

Let $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ be a PL-RDK system. Suppose it has an irreversible inflow reaction $0\to A$. Let $\mathcal{D}$ be the set of all reactions with $A$ in either reactant or product complex, not including the chosen reaction $0\to A$. If $A$ does not appear for each simplified reaction vector in $\mathcal{D}$, then the system does not have the capacity to admit multiple equilibria.

Proof.

Suppose an irreversible inflow reaction $0\to A$ exists. By definition of orientation, the reaction $0\to A$ must be an element of any orientation. Let $\mathscr{O}$ be an orientation. In STEP 2 of the algorithm, we are getting $KerL_{\mathscr{O}}$. At this point, we solve for the $\alpha_{y\to y^{\prime}}$’s given the equation $\displaystyle\sum\limits_{y\to y^{\prime}\in\mathscr{O}}{{\alpha_{y\to y^{\prime}}}\left({y^{\prime}-y}\right)}=0$. Since after simplifying the reaction vectors, $A$ does not appear in reaction vector $y^{\prime}-y$ for each reaction $y\to y^{\prime}$, the reaction $0\to A$ corresponds to a row with all entries 0. Since $0\to A$ is irreversible, then it must be placed on the zeroth equivalence class $P_{0}$. But for a system to have the capacity to admit multiple equilibria, each element in $P_{0}$ must be reversible (given in STEP 2). ∎

Proposition E.2**.**

Let $\left(\mathscr{S},\mathscr{C},\mathscr{R},K\right)$ be a PL-RDK system. Suppose it has an irreversible outflow reaction $B\to 0$. Let $\mathcal{D}$ be the set of all reactions with $B$ in either reactant or product complex, not including the chosen reaction $B\to 0$. If $B$ does not appear for each simplified reaction vector in $\mathcal{D}$, then the system does not have the capacity to admit multiple equilibria.

Proof.

The proof is similar to the proof of Proposition E.1. ∎

The propositions above are usually applied for total representation of a reaction network where there are several independent variables. In this representation, for each interaction $X_{1}\to X_{2}$ with regularity arrow from each element of $\{X_{j}\}$, the reaction $X_{1}+\sum X_{j}\to X_{1}+\sum X_{j}$ is associated, and any interaction without regulatory arrow are kept as it is [2]. By satisfying the assumptions, one can decide whether a system has the capacity to admit multiple equilibria.