A computational approach to multistationarity of power-law kinetic 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.

Appendices

Nomenclature

1.1 List of abbreviations

Abbreviation	Meaning
ADA	Advanced deficiency algorithm
BST	Biochemical Systems Theory
CKS	Chemical kinetic system
CRN	Chemical reaction network
CRNT	Chemical Reaction Network Theory
DOA	Deficiency one algorithm
DOT	Deficiency one theorem
GMA	Generalized mass action
HDA	Higher deficiency algorithm
MAK	Mass action kinetics
MSA	Multistationarity algorithm
PLK	Power-law kinetics
PL-NDK	Power-law non-reactant-determined kinetics
PL-RDK	Power-law reactant-determined kinetics
SFRF	Species formation rate function

1.2 List of important symbols

Meaning	Symbol
Deficiency	$\delta$
Dimension of the stoichiometric subspace	s
Equivalence class i	$P_i$
Fundamental class i	$C_i$
Lower shelf contained in $C_i$	$\mathcal{L}_i$
Middle shelf contained in $C_i$	$\mathcal{M}_i$
Molecularity matrix	Y
Number of linkage classes	l
Number of strong linkage classes	sl
Number of terminal strong linkage classes	t
Rate constant associated to $y \rightarrow y'$	$k_{y \rightarrow y'}$
Stoichiometric matrix	N
Stoichiometric subspace	S
Upper shelf contained in $C_i$	$\mathcal{U}_i$

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 5

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}}_{\ge 0}$. The ordered pair $\left( {{C_i},{C_j}} \right)$ corresponds to the reaction ${C_i} \rightarrow {C_j}$.

Definition 6

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

$$\begin{aligned} {\left( {{I_a}} \right) _{ij}} = \left\{ \begin{array}{rl} - 1&{}\mathrm{{ if }}{C_i}{\text { is in the reactant complex of reaction }}{R_j},\\ 1&{}\mathrm{{ if }}{C_i}{\text { is in the product complex of reaction }}{R_j},\\ 0&{}\mathrm{{ otherwise}}. \end{array} \right. \end{aligned}$$

The stoichiometric matrix, denoted by N, is the $m\times r$ matrix given by $N=YI_a$.

Definition 7

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 8

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

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

The rank of the network, denoted by s, is given by $s=\dim S$. The set $\left( {x + S} \right) \cap {\mathbb {R}}_{ \ge 0}^{\mathscr {S}}$ is said to be a stoichiometric compatibility class of $x \in {\mathbb {R}}_{ \ge 0}^{\mathscr {S}}$.

Definition 9

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

Definition 10

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 11

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 12

A kineticsK for a reaction network $\left( {\mathscr {S}},{\mathscr {C}},{\mathscr {R}}\right)$ is an assignment to each reaction $j: y \rightarrow y' \in {\mathscr {R}}$ of a rate function ${K_j}:{\varOmega _K} \rightarrow {{\mathbb {R}}_{ \ge 0}}$ such that ${\mathbb {R}}_{ > 0}^{\mathscr {S}} \subseteq {\varOmega _K} \subseteq {\mathbb {R}}_{ \ge 0}^{\mathscr {S}}$, $c \wedge d \in {\varOmega _K}$ if $c,d \in {\varOmega _K}$, and ${K_j}\left( c \right) \ge 0$ for each $c \in {\varOmega _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 13

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} \rightarrow {C_j} \in {\mathscr {R}}} {{K_{{C_i} \rightarrow {C_j}}}\left( x \right) \left( {{C_j} - {C_i}} \right) }.$

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 14

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 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 15

A kinetics K is a power-law kinetics (PLK) if ${K_i}\left( x \right) = {k_i}{{x^{{F_{i}}}}}$$\forall i =1,\ldots ,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 16

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,\ldots ,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 17

[17] The $m \times n$ matrix ${\widetilde{Y}}$ is given by

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

Definition 18

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

Proofs of the lemmas and the main theorem

Proof of Lemma 1

We assume positive equilibria ${c}^{*}$ and ${c}^{**}$ are distinct so ${\mu _s} \ne 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 \rightarrow y'}} = {k_{y \rightarrow y'}}{c}{^{**{T_{.y}}}} \ \forall \ y \rightarrow y' \in {\mathscr {R}}$ from the set of positive rate constants $\left\{ {{k_{y \rightarrow y'}}|y' \rightarrow y \in {\mathscr {R}}} \right\}$. $\square$

Proof of Lemma 2

Since $\mu \ne 0$ then $\mu _s \ne 0$ for some $s\in {\mathscr {S}}$. Hence, ${\dfrac{{{c_s}^*}}{{{c_s}^{**}}} = {e^{{\mu _s}}} \ne 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 \ne 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 \rightarrow y'}} = \dfrac{{{\kappa _{y \rightarrow y'}}}}{{{c}{{^{**}}^{T_{.y}}}}}.$$\square$

Proof of Lemma 4

Assume $y \rightarrow y'$ is irreversible. Then, ${g_{y \rightarrow y'}} = {\kappa _{y \rightarrow y'}} > 0$, ${h_{y \rightarrow y'}} = {\kappa _{y \rightarrow y'}}{e^{{T_{.y}} \cdot \mu }} > 0$, and ${\rho _{y \rightarrow y'}} = \dfrac{{{h_{y \rightarrow y'}}}}{{{g_{y \rightarrow y'}}}} = \dfrac{{{\kappa _{y \rightarrow y'}}{e^{{T_{.y}} \cdot \mu }}}}{{{\kappa _{y \rightarrow y'}}}} = {e^{{T_{.y}} \cdot \mu }}$. Assume $y \rightarrow y'$ is reversible. Suppose ${g_{y \rightarrow y'}} \ne 0$. Now, ${\rho _{y \rightarrow y'}}{g_{y \rightarrow y'}} = {h_{y \rightarrow y'}}$, ${g_{y \rightarrow y'}} = {\kappa _{y \rightarrow y'}} - {\kappa _{y' \rightarrow y}}$ and ${h_{y \rightarrow y'}} = {\kappa _{y \rightarrow y'}}{e^{{T_{.y}} \cdot \mu }} - {\kappa _{y' \rightarrow y}}{e^{{T_{.y'}} \cdot \mu }}$. Then, ${g_{y \rightarrow y'}}{e^{{T_{.y}} \cdot \mu }} = {\kappa _{y \rightarrow y'}}{e^{{T_{.y}} \cdot \mu }} - {\kappa _{y' \rightarrow y}}{e^{{T_{.y}} \cdot \mu }}$. Also, ${\kappa _{y' \rightarrow y}}{e^{{T_{.y}} \cdot \mu }} - {\kappa _{y' \rightarrow y}}{e^{{T_{.y'}} \cdot \mu }} = {\rho _{y \rightarrow y'}}{g_{y \rightarrow y'}} - {g_{y \rightarrow y'}}{e^{{T_{.y}} \cdot \mu }}$ and ${\kappa _{y' \rightarrow y}}\left( {{e^{{T_{.y}} \cdot \mu }} - {e^{{T_{.y'}} \cdot \mu }}} \right) = {g_{y \rightarrow y'}}\left( {{\rho _{y \rightarrow y'}} - {e^{{T_{.y}} \cdot \mu }}} \right)$. Hence, ${\kappa _{y' \rightarrow y}} = \dfrac{{{\rho _{y \rightarrow y'}} - {e^{{T_{.y}} \cdot \mu }}}}{{{e^{{T_{.y}} \cdot \mu }} - {e^{{T_{.y'}} \cdot \mu }}}}{g_{y \rightarrow y'}}.$ Similarly, we can solve for the value of ${\kappa _{y \rightarrow y'}}$ with ${\kappa _{y \rightarrow y'}} = \dfrac{{{\rho _{y \rightarrow y'}} - {e^{{T_{.y'}} \cdot \mu }}}}{{{e^{{T_{.y}} \cdot \mu }} - {e^{{T_{.y'}} \cdot \mu }}}}{g_{y \rightarrow y'}}.$ Thus, we obtain the following.

• i. If $y \rightarrow y'$ is irreversible then ${g_{y \rightarrow y'}} > 0$, ${h_{y \rightarrow y'}} > 0$, and ${\rho _{y \rightarrow y'}} = {e^{{T_{.y}} \cdot \mu }}$.

• ii. Suppose $y \rightarrow y'$ is reversible.

  • a. If ${g_{y \rightarrow y'}} > 0$, then either ${\rho _{y \rightarrow y'}} = {e^{{T_{.y}} \cdot \mu }} = {e^{{T_{.y'}} \cdot \mu }}$, ${\rho _{y \rightarrow y'}}> {e^{{T_{.y}} \cdot \mu }} > {e^{{T_{.y'}} \cdot \mu }}$, or ${\rho _{y \rightarrow y'}}< {e^{{T_{.y}} \cdot \mu }} < {e^{{T_{.y'}} \cdot \mu }}$.

  • b. If ${g_{y \rightarrow y'}} < 0$, then either ${\rho _{y \rightarrow y'}} = {e^{{T_{.y'}} \cdot \mu }} = {e^{{T_{.y}} \cdot \mu }}$, ${\rho _{y \rightarrow y'}}> {e^{{T_{.y'}} \cdot \mu }} > {e^{{T_{.y}} \cdot \mu }}$, or ${\rho _{y \rightarrow y'}}< {e^{{T_{.y'}} \cdot \mu }} < {e^{{T_{.y}} \cdot \mu }}$.

  • c. If ${g_{y \rightarrow y'}} = 0$ and ${h_{y \rightarrow y'}} > 0$ then ${e^{{T_{.y}} \cdot \mu }} > {e^{{T_{.y'}} \cdot \mu }}$.

  • d. If ${g_{y \rightarrow y'}} = 0$ and ${h_{y \rightarrow y'}} < 0$ then ${e^{{T_{.y}} \cdot \mu }} < {e^{{T_{.y'}} \cdot \mu }}$.

  • e. If ${g_{y \rightarrow y'}} = {h_{y \rightarrow y'}} = 0$ then ${e^{{T_{.y}} \cdot \mu }} = {e^{{T_{.y'}} \cdot \mu }}$.

Now, let $i \in \{1,2,\ldots ,w\}$ and $y \rightarrow y' \in P_i$. Also, let ${y_i\rightarrow y_{i}'}$ be the representative of $P_i$. Then ${\rho _{y_i \rightarrow y_i'}}=\dfrac{{{h_{{y_i} \rightarrow {y_i}'}}}}{{{g_{{y_i} \rightarrow {y_i}'}}}}=\dfrac{{\alpha _{y_i \rightarrow y_i'}}{{h_{{y_i} \rightarrow {y_i}'}}}}{{\alpha _{y_i \rightarrow y_i'}}{{g_{{y_i} \rightarrow {y_i}'}}}}=\dfrac{{{h_{{y} \rightarrow {y}'}}}}{{{g_{{y} \rightarrow {y}'}}}}={\rho _{y \rightarrow y'}}$ for some nonzero ${\alpha _{y_i \rightarrow y_i'}}$. $\square$

Proof of Lemma 5

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

Lemma 7

[15] Suppose a reaction network satisfies the following properties for an orientation ${\mathscr {O}}$: $P_i$ ($i=0,1,2,\ldots ,w$) is defined by a representative ${{y_i} \rightarrow {y_i}'}$, $W=\{{{y_i} \rightarrow {y_i}'}|i=1,2,\ldots ,w\}\subseteq {\mathscr {O}}$, a given basis $\left\{ {{b^j}} \right\} _{j = 1}^q$ for $Ker^\perp {L_{\mathscr {O}}}\cap {\varGamma _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} \rightarrow {y_i}'} \right) |{g_W}\left( {{y_i} \rightarrow {y_i}'} \right) \ne 0,i = 1,2,\ldots ,w} \right\}$ where the sign of ${{\rho _W}\left( {{y_i} \rightarrow {y_i}'} \right) }$ is the same as the ratio of the signs of ${{h_W}\left( {{y_i} \rightarrow {y_i}'} \right) }$ and ${{g_W}\left( {{y_i} \rightarrow {y_i}'} \right) }$. Then ${h_W},{g_W} \in {{\mathbb {R}}^{\mathscr {O}}} \cap {\varGamma _W}$ satisfy Eqs. (1), (2), and (3) if and only if ${\rho _W}\left( {{y_i} \rightarrow {y_i}'} \right)$’s satisfy the conditions in Lemma 6.

Proof of Theorem 4

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

$$\begin{aligned} {g_{y \rightarrow y'}}= & {} \left\{ \begin{array}{ll} {\kappa _{y \rightarrow y'}} - {\kappa _{y' \rightarrow y}}&{}\mathrm{{ if }} \ y \rightarrow y' \in {\mathscr {O}}{\text { is reversible}}\\ {\kappa _{y \rightarrow y'}}&{}\mathrm{{ if }} \ y \rightarrow y' \in {\mathscr {O}}{\text { is irreversible}} \end{array} \right. and\\ {h_{y \rightarrow y'}}= & {} \left\{ \begin{array}{ll} {\kappa _{y \rightarrow y'}}{e^{{T_{.y}} \cdot \mu }} - {\kappa _{y' \rightarrow y}}{e^{{T_{.y'}} \cdot \mu }} &{} \mathrm{{ if }} \ y \rightarrow y' \in {\mathscr {O}}\mathrm{{\ is \ reversible}}\\ {\kappa _{y \rightarrow y'}}{e^{{T_{.y}} \cdot \mu }} &{} \mathrm{{ if }} \ y \rightarrow y' \in {\mathscr {O}}\mathrm{{\ is \ irreversible}} \end{array} \right. . \end{aligned}$$

For $i=1,2,\ldots ,w$, let ${P_i}$ be the equivalence class with ${y_i} \rightarrow {y_i}'$ 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 {\varGamma _W}$, if it exists. Also, let $g,h \in Ker {L_{\mathscr {O}}}$. Hence, $b^j \cdot g = \sum \limits _{{y_i} \rightarrow {y_i}' \in W} {g_{{y_i} \rightarrow {y_i}'} b_{{y_i} \rightarrow {y_i}'}^j= 0}$ and $b^j \cdot h = \sum \limits _{{y_i} \rightarrow {y_i}' \in W} {h_{{y_i} \rightarrow {y_i}'} b_{{y_i} \rightarrow {y_i}'}^j= 0}$ for $j=1,2,\ldots ,q$. These two equations can be written as Eqs. (1), (2), and (3). Hence, (iii) which leads to (iv).

On the other hand, by Lemma 2, we only need to show the existence of $\kappa$. By Lemma 7 which requires the existence of the forest basis, Eqs. (1), (2), and (3) hold. From here, we obtain the values of $g_{y \rightarrow y'}$ and $h_{y \rightarrow y'}$. These yield $\kappa$. $\square$

Computation details for Sect. 3

1.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 \rightarrow y'}=\dfrac{\kappa _{y \rightarrow y'}}{c^{{**}{T_{.y}}}}$ is indicated. Note that $\kappa _{y \rightarrow y'}$ was solved by getting the scalars in the linear combination of a basis for $Ker L_{{\mathscr {R}}}$ such that

$$\begin{aligned} \sum \limits _{y \rightarrow y' \in {{\mathscr {R}}}} {{\kappa _{y \rightarrow y'}}\left( {y' - y} \right) = 0} \mathrm{{ \ and \ }} \sum \limits _{y \rightarrow y' \in {{\mathscr {R}}}} {{\kappa _{y \rightarrow y'}}{e^{{T_{.y }\cdot \mu }}}\left( {y' - y} \right) = 0}. \end{aligned}$$

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

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

We choose the forward reactions so ${{\mathscr {O}}}=\{ {A_1} + 2{A_2} \rightarrow 2{A_1} + {A_2},{A_1} + {A_2} \rightarrow 2{A_2},{A_2} \rightarrow {A_3} \}.$ A basis for $Ker{L_{\mathscr {O}}}$, and the corresponding equivalence and fundamental classes are given below.

$$\begin{aligned} \begin{array}{*{20}{c}} {{A_1} + 2{A_2} \rightarrow 2{A_1} + {A_2}}\\ {{A_1} + {A_2} \rightarrow 2{A_2}}\\ {{A_2} \rightarrow {A_3}} \end{array}\left( {\begin{array}{*{20}{c}} 1\\ 1\\ 0 \end{array}} \right) \end{aligned}$$

${P_0} = \left\{ {{A_2}\rightarrow {A_3}} \right\} \subset {C_0} = \left\{ {{A_2} \rightarrow {A_3},{A_3} \rightarrow {A_2}} \right\}$

${P_1} = \left\{ {{A_1} + 2{A_2} \rightarrow 2{A_1} + {A_2},{A_1} + {A_2} \rightarrow 2{A_2}} \right\} ={C_1}$

We choose $W = \left\{ {{A_1} + {A_2} \rightarrow 2{A_2}} \right\}$ with $w=1$. The resulting inequality system is automatically linear since $Ker^{\perp }L_{{\mathscr {O}}}\cap \varGamma _W$ is trivial. Since $P_1$ is nonreversible, the sign pattern for ${g_W}$ and ${h_W}$ for ${y_1} \rightarrow {y_1}'$ must be positive. We have the following shelving assignment: ${\mathcal{U}_1} = \left\{ \right\} ,{\mathcal{M}_1} = \left\{ {{A_1} + 2{A_2} \rightarrow 2{A_1} + {A_2},{A_1} + {A_2} \rightarrow 2{A_2}} \right\} ,{\mathcal{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 \varGamma _W$ is trivial, no additional equation or inequality can be obtained here. We have the following system.

$$\begin{aligned} \begin{aligned} {p_1}{\mu _{{A_1}}} + {q_1}{\mu _{{A_2}}}&= {M_1} = {p_2}{\mu _{{A_1}}} + {q_2}{\mu _{{A_2}}}\\ {\mu _{{A_2}}}&= {\mu _{{A_3}}} \end{aligned} \end{aligned}$$

(4)

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

1.3 A remark on the steps of the MSA

Remark 5

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 \varGamma _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.

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 1

Let $\left( {\mathscr {S}},{\mathscr {C}},{\mathscr {R}},K\right)$ be a PL-RDK system. Suppose it has an irreversible inflow reaction $0 \rightarrow A$. Let ${\mathcal {D}}$ be the set of all reactions with A in either reactant or product complex, not including the chosen reaction $0 \rightarrow 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 \rightarrow A$ exists. By definition of orientation, the reaction $0 \rightarrow A$ must be an element of any orientation. Let ${\mathscr {O}}$ be an orientation. In STEP 2 of the algorithm, we are getting $Ker L_{\mathscr {O}}$. At this point, we solve for the $\alpha _{y \rightarrow y'}$’s given the equation $\displaystyle \sum \limits _{y \rightarrow y' \in {\mathscr {O}}}^{} {{\alpha _{y \rightarrow y'}}\left( {y' - y} \right) }=0$. Since after simplifying the reaction vectors, A does not appear in reaction vector $y'-y$ for each reaction $y \rightarrow y'$, the reaction $0 \rightarrow A$ corresponds to a row with all entries 0. Since $0 \rightarrow 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). $\square$

Proposition 2

Let $\left( {\mathscr {S}},{\mathscr {C}},{\mathscr {R}},K\right)$ be a PL-RDK system. Suppose it has an irreversible outflow reaction $B \rightarrow 0$. Let ${\mathcal {D}}$ be the set of all reactions with B in either reactant or product complex, not including the chosen reaction $B \rightarrow 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 1. $\square$

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 \rightarrow X_2$ with regularity arrow from each element of $\{X_j\}$, the reaction $X_1+\sum X_j \rightarrow 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.