Complex-balanced equilibria of generalized mass-action systems: necessary conditions for linear stability

1.   Introduction

In their foundational paper from 1972, Horn and Jackson considered chemical reaction networks (CRNs) with mass-action kinetics ^[1]. In particular, they proved that complex-balanced equilibria are asymptotically stable (for all rate constants), by using an entropy-like Lyapunov function. However, mass-action kinetics is an assumption that holds for elementary reactions in homogeneous and dilute solutions. In intracellular environments, which are highly structured and crowded, more general kinetics such as power laws are needed ^[2,3]. Already Horn and Jackson observed that every CRN with power-law kinetics can be written as another CRN with mass-action kinetics (possibly with non-integer stoichiometric coefficients). Typically, the resulting network does not have desired properties such as weak reversibility and zero deficiency. The more recent definition of CRNs with generalized mass-action kinetics (involving both stoichiometric coefficients and kinetic orders) allows to study power-law kinetics without having to rewrite the network ^[4,5]. Moreover, a CRN with classical mass-action kinetics may not have zero deficiency and may not be weakly reversible, but there may be a dynamically equivalent CRN with generalized mass-action kinetics that has the desired properties ^[6].

For the resulting generalized mass-action systems, existence and uniqueness of complex-balanced equilibria (in every stoichiometric class and for all rate constants) is well understood ^[4,7], but not much is known about their stability. Even planar S-systems with a unique complex-balanced equilibrium display rich dynamical behavior, including super/sub-critical or degenerate Hopf bifurcations, centers, and up to three limit cycles; see ^[8,9,10,11].

For chemical reaction networks with classical mass-action kinetics, linear Lyapunov functions have been used in approaches to the global attractor conjecture ^[12,13]. For certain classes of networks, in particular, non-autocatalytic networks with monotonic kinetics, piecewise linear/polyhedral Lyapunov functions have been constructed ^[14,15]. However, no Lyapunov functions are known for complex-balanced equilibria of arbitrary generalized mass-action systems, and hence we first approach the problem by linearization. In other words, instead of asymptotic stability, we investigate linear stability (on the stoichiometric subspace and for all rate constants). To this end, we first discuss several notions of matrix stability (on a linear subspace) such as D-stability and diagonal stability; see [16,Ch. 2], [17,Ch. 26] and the references therein and ^[18] for the first application to chemical reaction networks. In particular, diagonal stability of the Jacobian matrix allows to construct Lyapunov functions with separated variables. Our main results characterize linear stability of complex-balanced equilibria (on the stoichiometric subspace and for all rate constants) for cyclic networks and give necessary conditions for weakly reversible networks.

In the setting of classical mass-action systems, we prove that complex-balanced equilibria are not just asymptotically stable, but even diagonally stable (and hence linearly stable). For an alternative proof, see [19,Thm. 15.2.2.]. As opposed to asymptotic stability, linear stability is robust with respect to small perturbations of the system. In particular, this allows to show the robustness of the classical deficiency zero theorem with respect to small perturbations of the kinetic orders from the stoichiometric coefficients; see [7,Cor. 47].

Organization of the work and main results. In Section 2, we introduce generalized mass-action systems, and in Section 3, we discuss several notions of matrix stability (on a linear subspace). In Section 4, we present our main results: For classical mass-action systems, complex-balanced equilibria are diagonally stable (and hence linearly stable); see Theorem 4.1. For generalized mass-action systems, linear stability of complex balanced equilibria implies uniqueness; see Theorem 4.3. For cyclic networks, linear stability is equivalent to D-stability of the Jacobian matrix; see Theorem 11. For weakly reversible networks, linear stability implies D-semistability of the Jacobian matrices of all cycles in the network; see Theorem 13. In Section 5, we recall and extend characterizations of notions of matrix stability, and finally, in Section 6, we illustrate our results by a series of examples.

2.   Generalized mass-action systems

A generalized chemical reaction network $(G, y, \tilde y)$ is based on a directed graph $G = (V, E)$ without self-loops; every vertex $i \in V = \{1, \ldots, m\}$ is labeled with a (stoichiometric) complex $y(i) \in {\mathbb R}^n_{\ge0}$, and every source vertex $i \in V_s \subseteq V$ is labeled with a kinetic-order complex $\tilde y(i) \in {\mathbb R}^n$. If every component of $G$ is strongly connected, we call $G$ and $(G, y, \tilde y)$ weakly reversible.

A generalized mass-action system $(G_k, y, \tilde y)$ is a generalized chemical reaction network $(G, y, \tilde y)$ together with edge labels $k \in {\mathbb R}^E_+$, resulting in the labeled digraph $G_k$. Every edge $(i \to i') \in E$, representing the chemical reaction $y(i) \to y(i')$, is labeled with a rate constant $k_{i \to i'} > 0$.

The ODE system for the species concentrations $x \in {\mathbb R}^n_+$, associated with the generalized mass-action system $(G_k, y, \tilde y)$, is given by

The sum ranges over all reactions, and every summand is a product of the reaction rate and the difference of product and educt complexes. Thereby, for $x \in {\mathbb R}^n_+$ and $y \in {\mathbb R}^n$, we define $x^y = x_1^{y_1} \cdots x_n^{y_n} \in {\mathbb R}_+$. Moreover, for $Y = (y^1, \ldots, y^m) \in {\mathbb R}^{n \times m}$, we define $x^Y \in {\mathbb R}^m_+$ via $(x^Y)_j = x^{y^j}$ for $j = 1, \ldots, m$.

Example 1. Consider the Lotka reactions $\mathsf{0} \to \mathsf{X}$, $\mathsf{X} \to \mathsf{Y}$, and $\mathsf{Y} \to \mathsf{0}$ with generalized mass-action kinetics

Thereby, the first vertex is labeled with the stoichiometric complex $\mathsf{0}$ and the kinetic-order complex $\alpha_1 \mathsf{X} + \beta_1 \mathsf{Y}$ in brackets. (Analogously, for the second and third vertex.) Hence, the reaction rate of $\mathsf{0} \to \mathsf{X}$ is given by $k_{1\to2} \, x^{\alpha_1} y^{\beta_1}$.

The associated ODE system is given by

The right-hand-side of the ODE system (1) can be decomposed into stoichiometric, graphical, and kinetic-order contributions,

where $Y, \, \widetilde Y \in {\mathbb R}^{n \times V}$ are the matrices of stoichiometric and kinetic-order complexes, $I_E, I_E^s \in {\mathbb R}^{V \times E}$ are the incidence and source matrices^* of the digraph $G$, and $A_k = I_E {\rm{diag}}(k) (I_E^s)^ \mathsf{T} \in {\mathbb R}^{V \times V}$ is the Laplacian matrix of the labeled digraph $G_k$. (We note that columns of $\widetilde Y$ corresponding to non-source vertices can be chosen arbitrarily.)

^*Explicitly, $(I_E)_{i, j \to j'} = -1$ if $i = j$, $(I_E)_{i, j \to j'} = 1$ if $i = j'$, and $(I_E)_{i, j \to j'} = 0$ otherwise. Further, $(I_E^s)_{i, j \to j'} = 1$ if $i = j$ and $(I_E^s)_{i, j \to j'} = 0$ otherwise. Finally, $(A_k)_{i', i} = k_{i \to i'}$ if $(i \to i') \in E$, $(A_k)_{i, i} = - \sum_{i \to i'} k_{i \to i'}$, and $(A_k)_{i', i} = 0$ otherwise.

Clearly, the change over time lies in the stoichiometric subspace

that is, $\frac{\mathrm{d} {x}}{\mathrm{d} {t}} \in S$. Equivalently, trajectories are confined to cosets of $S$, that is, $x(t) \in x(0)+S$. For $x' \in {\mathbb R}^n_+$, the set $(x'+S) \cap {\mathbb R}^n_+$ is called a stoichiometric class.

Analogously, we introduce the kinetic-order subspace

A positive equilibrium $x \in {\mathbb R}^n_+$ of the ODE system (2) that fulfills $A_k \, x^{\widetilde Y} = 0$ is called a complex-balanced equilibrium. The set of all complex-balanced equilibria is given by

The Jacobian matrix of the right-hand-side of (2) is given by

For given generalized chemical reaction network $(G, y, \tilde y)$, we study whether, for all rate constants $k\in {\mathbb R}^E_+$ (such that $Z_k\neq\emptyset$), all complex-balanced equilibria $x^* \in Z_k$ are linearly stable on $S$, that is, the corresponding Jacobian matrices $J(x^*)$ are stable on $S$.

Before we turn to stability, we recall the well-known fact that every positive vector is a complex-balanced equilibrium for some rate constant (see e.g. the proof of Lemma 1 in ^[5]). Further, we state a characterization of the uniqueness of complex-balanced equilibria in terms of sign vectors of the stoichiometric and kinetic-order subspaces (see Proposition 3.1 in ^[4]).

Proposition 2. Consider a weakly reversible generalized chemical reaction network, and let $x^* \in {\mathbb R}_+^n$. Then, there exist rate constants $k \in {\mathbb R}^E_+$ such that $x^* \in Z_k$.

Proposition 3. Consider a weakly reversible generalized chemical reaction network. There exists at most one complex-balanced equilibrium in every stoichiometric class and for all rate constants if and only if ${\rm{sign}}(S) \cap {\rm{sign}}(\widetilde S^\perp) \neq \{0\}$.

For surveys on the uniqueness of equilibria and related injectivity results, see ^[20,21,22].

3.   Notions of matrix stability

Let $S$ be a linear subspace of ${\mathbb R}^n$. For an ODE system $\dot x = f(x)$ with $f \colon {\mathbb R}^n \to {\mathbb R}^n$ and ${\rm{im}} f \subseteq S$, cosets of $S$ are forward invariant. Hence, given an equilibrium $x \in {\mathbb R}^n$, one is interested in its asymptotic stability on the coset $x+S$. To approach the problem via linearization, one considers the Jacobian matrix $J(x) = \frac{\partial {f}}{\partial {x}} \in {\mathbb R}^{n \times n}$, more precisely, the linear map $J(x)|_S \colon S \to S$. By the Hartman-Grobman Theorem, if $x$ is hyperbolic (that is, if all eigenvalues of $J(x)|_S$ have non-zero real part), then the original ODE system is dynamically equivalent (technically: topologically conjugate) to the linear ODE system $\dot y = J(x) y$ on $S$; see e.g. [23,Thm. 9.9]. In the following, we recall notions of stability of a square matrix and extend them to stability on a linear subspace.

A square matrix is stable (respectively, semistable) if all its eigenvalues have negative (respectively, non-positive) real part. We start with a matrix formulation of Lyapunov's Theorem; see [24,Ch. XV,Thm. 3'] or [16,Thm. 2.2.1].

Proposition 4. Let $A\in {\mathbb R}^{n \times n}$. The following implications hold:

Proof. Obviously, the implications from left to right hold.

The equivalence on the left is Lyapunov's Theorem.

Finally, if there exists $P = P^ \mathsf{T} > 0$ with $PA+A^ \mathsf{T} P\leq 0$, then the origin is Lyapunov stable for the linear ODE $\dot x = Ax$ (by the Lyapunov function $x \mapsto x^ \mathsf{T} P x$), and thus $A$ cannot have an eigenvalue with positive real part.

Let $S$ be a linear subspace. We say that a square matrix $A$ with ${\rm{im}} A \subseteq S$ is stable on $S$ (respectively, semistable on $S$) if all eigenvalues of the linear map $A|_S\colon S\to S$ have negative (respectively, non-positive) real part. We extend Lyapunov's Theorem to stability on a linear subspace.

Proposition 5. Let $A \in {\mathbb R}^{n \times n}$ and $S \subseteq {\mathbb R}^n$ be a linear subspace with ${\rm{im}} A \subseteq S$. The following implications hold:

Proof. Obviously, the implications from left to right hold.

To prove the equivalence on the left, let $s = \dim S$ and $S = {\rm{im}} B$ with $B \in {\mathbb R}^{n\times s}$ and $B^ \mathsf{T} B = I \in {\mathbb R}^{s \times s}$. Then, for every $x\in S$, there exists a unique $y\in {\mathbb R}^s$ such that $x = By$. Further, $y = B^ \mathsf{T} x$ and, using $\dot x = Ax$ (on $S$), we have $\dot y = B^ \mathsf{T} A B y$ (on ${\mathbb R}^s$). Thus, $A$ is stable on $S$ if and only if $B^ \mathsf{T} A B \in {\mathbb R}^{s \times s}$ is stable on ${\mathbb R}^s$. By Proposition 3.1, the latter is equivalent to the existence of $Q \in {\mathbb R}^{s\times s}$ with $Q = Q^ \mathsf{T} > 0$ such that $y^ \mathsf{T}(Q B^ \mathsf{T} A B + B^ \mathsf{T} A^ \mathsf{T} BQ) y < 0$ for all $0\neq y \in {\mathbb R}^s$. This is further equivalent to the existence of $P = B Q B^ \mathsf{T} \in {\mathbb R}^{n \times n}$ with $P = P^ \mathsf{T} > 0$ on $S$ such that $x^ \mathsf{T} (PA+A^ \mathsf{T} P) x < 0$ for all $0\neq x \in S$.

Finally, let $P = P^ \mathsf{T} > 0$ on $S$ and $PA+A^ \mathsf{T} P\leq0$ on $S$. With $B$ as above, let $Q = B^ \mathsf{T} P B \in {\mathbb R}^{s \times s}$ and note that $A = B B^ \mathsf{T} A$ (because $x = B B^ \mathsf{T} x$ for all $x \in S$). Then $Q = Q^ \mathsf{T} > 0$ and $QB^ \mathsf{T} A B + B^ \mathsf{T} A^ \mathsf{T} B Q \le 0$. By Proposition 3.1, $B^ \mathsf{T} A B$ is semistable, and, as in the previous paragraph, the latter is equivalent to $A$ being semistable on $S$.

We recall more notions of stability of square matrices. For convenience, let $\mathcal{D}_+ \subseteq {\mathbb R}^{n \times n}$ denote the set of diagonal matrices with positive diagonal. A matrix $A \in {\mathbb R}^{n\times n}$ is diagonally stable (respectively, diagonally semistable) if there exists $P \in \mathcal{D}_+$ such that $PA+A^ \mathsf{T} P < 0$ (respectively, $PA+A^ \mathsf{T} P\leq0$). We note that diagonal stability is also known as Lyapunov diagonal stability or Volterra-Lyapunov stability. A matrix $A \in {\mathbb R}^{n\times n}$ is D-stable (respectively, D-semistable) if $AD$ is stable (respectively, semistable) for all $D \in \mathcal{D}_+$. The following statement summarizes the relations between these notions.

Proposition 6. Let $A\in {\mathbb R}^{n \times n}$. The following implications hold:

Proof. Obviously, the implications from left to right hold.

To prove the implications from the first to the second row, let $P \in \mathcal{D}_+$ be such that $PA+A^ \mathsf{T} P < 0$ (respectively, $PA+A^ \mathsf{T} P\leq0$). Then, for any $D \in \mathcal{D}_+$, we have $(DP)(AD)+(AD)^ \mathsf{T}(DP) = D(PA+A^ \mathsf{T} P)D < 0$ (respectively, $\le 0$), and by Proposition 3.1, $AD$ is stable (respectively, semistable).

Finally, the implications from the second row to the third row are trivial. (If $AD$ is (semi)stable for any $D\in\mathcal{D}_+$, then this holds for $D = I$.)

Let $S$ be a linear subspace. We say that $A \in {\mathbb R}^{n \times n}$ with ${\rm{im}} A \subseteq S$ is diagonally stable on $S$ (respectively, diagonally semistable on $S$) if there exists $P \in \mathcal{D}_+$ such that $PA+A^ \mathsf{T} P < 0$ on $S$ (respectively, $PA+A^ \mathsf{T} P\le0$ on $S$). We say that $A \in {\mathbb R}^{n \times n}$ with ${\rm{im}} A \subseteq S$ is D-stable on $S$ (respectively, D-semistable on $S$) if $AD$ is stable on $S$ (respectively, semistable on $S$) for all $D \in \mathcal{D}_+$. Finally, we introduce an even stronger notion. We say that $A \in {\mathbb R}^{n \times n}$ with ${\rm{im}} A \subseteq S$ is diagonally D-stable on $S$ (respectively, diagonally D-semistable on $S$) if, for all $D \in \mathcal{D}_+$, there exists $P \in \mathcal{D}_+$ such that $PAD+DA^ \mathsf{T} P < 0$ on $S$ (respectively, $PAD+DA^ \mathsf{T} P\le0$ on $S$). We note that diagonal D-stability on $S$ trivially implies diagonal stability on $S$, and the two notions agree for $S = {\mathbb R}^n$, see also [25,p. 257]. Moreover, we have the following relations.

Proposition 7. Let $A \in {\mathbb R}^{n \times n}$ and $S \subseteq {\mathbb R}^n$ be a linear subspace with ${\rm{im}} A \subseteq S$. The following implications hold:

Proof. Obviously, the implications from left to right hold.

The implications from the first to the second row follow immediately from Proposition 3.2.

Finally, the implications from the second row to the third row are trivial.

4.   Main results

In the following, we say that an equilibrium $x^*$ is linearly stable (in its stoichiometric class $x^*+S$) if the Jacobian matrix $J(x^*)$ is stable on $S$. In short:

Analogously, we say that $x^*$ is diagonally D-stable/diagonally stable/D-stable (in $x^*+S$) if $J(x^*)$ is diagonally D-stable/diagonally stable/D-stable on $S$.

4.1.   Linear stability for mass-action systems

We start by showing that, for classical mass-action systems (where $\tilde{y} = y$), complex-balanced equilibria are diagonally D-stable (and hence diagonally stable/D-stable/linearly stable) in their stoichiometric classes, and not just locally asymptotically stable (as shown via the classical entropy-like Lyapunov function).

Linear stability of complex-balanced equilibria was shown in [26,Theorem 4.3.2] and [19,Theorem 15.2.2.]. In the latter work, the author also proved diagonal stability, without noting it and using an ad-hoc inequality. Here, we even show diagonal D-stability, thereby using the negative semi-definiteness of the Laplacian matrix of an undirected graph.

Theorem 8. Consider a weakly reversible chemical reaction network. Then, for all rate constants, complex-balanced equilibria are linearly stable (in their stoichiometric classes). In fact, they are diagonally D-stable (and hence also diagonally stable, D-stable, and linearly stable).

Proof. Let $k \in {\mathbb R}^E_+$ and $x^* \in Z_k$. We show that the corresponding Jacobian matrix $J = J(x^*)$ is diagonally D-stable on $S$. By Proposition 3.4, all other conclusions follow.

Let $D \in \mathcal{D}_+$. We show that there exists $P = {\rm{diag}}((x^*)^{-1})D \in \mathcal{D}_+$ such that $H = P J D + DJ^ \mathsf{T} P < 0$ on $S$. Let $k^* \in {\mathbb R}^E_+$ be defined by $k^*_{i \to i'} = k_{i \to i'} \, (x^*)^{y(i)}$ for $(i \to i') \in E$. Then, $A_k {\rm{diag}}((x^*)^Y) = A_{k^*}$. Using (3) for the Jacobian matrix $J$, we have

Now, let $A_{\bar k} = A_{k^*} + A_{k^*}^ \mathsf{T}$ and hence $H = P \, Y A_{\bar k} \, Y^ \mathsf{T} P$. The symmetric matrix $A_{\bar k}$ is the Laplacian matrix of a labeled undirected graph ${\bar G}_{\bar k}$ with ${\bar G} = (V, {\bar E})$ and $\bar k \in {\mathbb R}^{{\bar E}}_+$. In particular, $A_{\bar{k}} = -I_{E'} {\rm{diag}}(\bar k) \, I_{E'}^ \mathsf{T}$ for some directed version $E'$ of the undirected edges ${\bar E}$ and the corresponding incidence matrix $I_{E'}$ (with ${\rm{im}} I_{E'} = {\rm{im}} I_E$); see e.g. ^[27]. Hence, $A_{\bar k} \le 0$ and $H = - P \, Y I_{E'} {\rm{diag}}(\bar k) \, I_{E'}^ \mathsf{T} \, Y^ \mathsf{T} P \le 0$ on $S$.

Suppose $v^ \mathsf{T} Hv = 0$ for some $v \in S$. Then, $I_{E'}^ \mathsf{T} \, Y^ \mathsf{T} Pv = 0$. Now, $S = {\rm{im}} (Y I_E) = {\rm{im}} (Y I_{E'})$, $S^\perp = \ker (I_{E'}^ \mathsf{T} \, Y^ \mathsf{T})$, and hence $Pv \in S^\perp$. Clearly, $v^ \mathsf{T} P v = 0$ which finally implies $v = 0$, since $P > 0$.

Hence, $H < 0$ on $S$, and $J$ is diagonally D-stable on $S$.

In the result above, we prove diagonal stability of $J$ on $S$, by providing a diagonal, positive definite matrix $P$ such that $PJ+J^ \mathsf{T} P < 0$ on $S$. This property implies the existence of a function with separated variables $L(x) = \sum_{i = 1}^nL_i(x_i)$ that serves as a Lyapunov function for showing the asymptotic stability of the complex-balanced equilibrium $x^*\in {\mathbb R}^n_+$ (in its stoichiometric class).

As an example, consider a scaled version of the classical entropy-like Lyapunov function $L(x) = \sum_{i = 1}^n p_i \, x_i^* \, [x_i(\log(x_i/x_i^*)-1)+x_i^*]$, where $P = {\rm{diag}}(p_1, \ldots, p_n)$, and the function $g(x) = \nabla L(x) f(x)$, where $\frac{\mathrm{d} {x}}{\mathrm{d} {t}} = f(x)$. Then, the gradient of $g$ at $x^*$ vanishes, and the Hessian matrix of $g$ at $x^*$ is $H = PJ+J^ \mathsf{T} P$ with $H < 0$ on $S$. Hence, in $x^* + S$, the function $g$ has a local maximum at $x^*$ with $g(x^*) = 0$. That is, $L$ is a Lyapunov function at $x^*$, and $x^*$ is asymptotically stable (in $x^* + S$).

4.2.   Linear stability implies uniqueness

Now, we turn to generalized mass-action systems. We show that linear stability of complex-balanced equilibria implies uniqueness (in their stoichiometric classes). We start with a useful technical result.

Lemma 9. Consider a weakly reversible generalized mass-action system, and let $x^* \in Z_k$. Then, $\ker (A_k {\rm{diag}}((x^*)^{\widetilde Y})) = \ker I_E^ \mathsf{T}$.

Proof. The dimension of $\ker I_E^ \mathsf{T}$ equals the number of connected components of the graph $G = (V, E)$, and the characteristic vectors of the vertex sets of the connected components form a basis.

Let $k^* \in {\mathbb R}^E_+$ be defined by $k^*_{i \to i'} = k_{i \to i'} \, (x^*)^{\tilde{y}(i)}$ for $(i \to i') \in E$. Then, $A_k {\rm{diag}}((x^*)^{\widetilde Y}) = A_{k^*}$. Clearly, $A_{k^*}$ is the Laplacian matrix of the labeled directed graph $G_{k^*}$, and the dimension of its kernel also equals the number of connected components of $G$; see e.g. [5,Section 4] and the references therein.

By definition, $x^* \in Z_k$ if $A_k (x^*)^{\widetilde Y} = 0$. Due to the block structure of the Laplacian matrix, the characteristic vectors of the vertex sets of the connected components are in the kernel of $A_k {\rm{diag}}((x^*)^{\widetilde Y})$. (As an example, consider a strongly connected graph $G$, and let $1^V$ be the characteristic vector of $V$. Then, $A_k {\rm{diag}}((x^*)^{\widetilde Y}) \, 1^V = A_k (x^*)^{\widetilde Y} = 0$.) Altogether, $\ker (A_k {\rm{diag}}((x^*)^{\widetilde Y})) = \ker I_E^ \mathsf{T}$.

Theorem 10. Consider a weakly reversible generalized chemical reaction network. If, for all rate constants, complex-balanced equilibria are linearly stable (in their stoichiometric classes), then they are unique (in their stoichiometric classes).

Proof. Assume that, for some $k' \in {\mathbb R}^E_+$, there exists a stoichiometric class with at least two complex-balanced equilibria. Then, by Proposition 2.3, there exist vectors $0\neq u \in \widetilde S^\perp$ and $0\neq v \in S$ with ${\rm{sign}}(u) = {\rm{sign}}(v)$. Now, $\widetilde S = {\rm{im}} (\widetilde Y I_E)$, $\widetilde S^\perp = \ker (I_E^ \mathsf{T} \, \widetilde Y^ \mathsf{T})$, and hence $I_E^ \mathsf{T} \, \widetilde Y^ \mathsf{T} u = 0$. Let $x^*\in {\mathbb R}^n_+$ be such that $u = {\rm{diag}}((x^*)^{-1}) v$, and let $k\in {\mathbb R}^E_+$ be such that $x^*\in Z_k$, see Proposition 2.2. Then,

that is, $\widetilde Y^ \mathsf{T} {\rm{diag}}((x^*)^{-1}) v \in \ker I_E^ \mathsf{T} = \ker (A_k {\rm{diag}}((x^*)^{\widetilde Y}))$, the latter by Lemma 4.2. Using (3) for the Jacobian matrix $J$, we have

Hence, for all $P = P^ \mathsf{T} > 0$, we have $v^ \mathsf{T}(PJ+J^ \mathsf{T} P)v = 0$, and by Proposition 3.2, $J$ is not stable on $S$.

As can be demonstrated by the generalized reaction network $\mathsf{X} (\mathsf{X}) \rightleftharpoons \mathsf{Y} (\mathsf{0})$, linear stability implies uniqueness, but not existence in every stoichiometric class.

4.3.   The network is a cycle

In the following result, we characterize diagonal stability and linear stability of complex-balanced equilibria, provided that the network is a cycle.

Theorem 11. Consider a cyclic generalized chemical reaction network, and let $A = Y A_{k = 1} \widetilde{Y}^ \mathsf{T}$. The following implications hold:

Proof. The implication in the first (respectively, second) row follows immediately from Proposition 3.2 (respectively, Proposition 3.4).

To prove the two equivalences, let $k \in {\mathbb R}^E_+$ and $x^* \in Z_k$. Since $G = (V, E)$ is a cycle, the quantity $c = k_{i \to i'} \, (x^*)^{\tilde{y}(i)}$ is the same for all $(i \to i')\in E$. In particular, $A_k {\rm{diag}}((x^*)^{\widetilde{Y}}) = c \, A_{k = 1}$. Using (3) for the Jacobian matrix $J$, we have

with $D = {\rm{diag}}((x^*)^{-1})$. Finally, recall that every $x^*\in {\mathbb R}^n_+$ is a complex-balanced equilibrium for some $k\in {\mathbb R}^E_+$, see Proposition 2.2. Hence, every $D \in \mathcal{D}_+$ is of the form $D = {\rm{diag}}((x^*)^{-1})$ for some $k$, and the equivalences follow.

4.4.   The network is weakly reversible

In the previous subsection, we characterized diagonal stability and linear stability of complex-balanced equilibria, provided that the network is a cycle. In this subsection, we extend those results to weakly reversible networks, but instead of equivalent we only give necessary conditions for diagonal and linear stability. We start with a useful technical result.

Lemma 12. Consider a weakly reversible generalized chemical reaction network, let $C$ be a cycle of the network, and let $x^* \in {\mathbb R}_+^n$. Then, there exists a family of rate constants $k^\varepsilon \in {\mathbb R}^E_+$ (with $\varepsilon > 0$) such that $x^* \in Z_{k^\varepsilon}$ for every $\varepsilon > 0$ and $J^\varepsilon \to Y A_{k = 1}^C \widetilde{Y}^ \mathsf{T} {\rm{diag}}((x^*)^{-1})$ as $\varepsilon \to 0$, where $A_{k = 1}^C$ is the Laplacian matrix of the cycle $C$ with all rate constants set to $1$.

Proof. For a cycle $C'$, let $k^{C'} \in {\mathbb R}^E_{\geq0}$ be defined by

For $\varepsilon > 0$, let $k^\varepsilon = k^C + \varepsilon \sum_{C'\neq C} k^{C'}$, where the summation is over all cycles $C'$ of $G$, except for the given cycle $C$. Then, $x^* \in Z_{k^\varepsilon}$ and

Using (3) for the Jacobian matrix $J^\varepsilon$, we obtain the desired limit as $\varepsilon \to 0$.

Theorem 13. Consider a weakly reversible generalized chemical reaction network, and, for a cycle $C$ of the network, let $A^C = Y A^C_{k = 1} \widetilde{Y}^ \mathsf{T}$ and $S^C$ be the corresponding stoichiometric subspace. The following implications hold:

Proof. The implications in the first (respectively, second) row follows immediately from Proposition 3.2 (respectively, Proposition 7).

Next, we prove the implication in the right column. Assume there exists a cycle $C$ in $G$ and a matrix $D \in \mathcal{D}_+$ such that $A^CD$ has an eigenvalue with positive real part. Let $x^*\in {\mathbb R}^n_+$ be such that $D = {\rm{diag}}((x^*)^{-1})$. Further, let $k^\varepsilon \in {\mathbb R}^E_+$ (with $\varepsilon > 0$) be a family of rate constants as in Lemma 4.5, and let $J^\varepsilon$ denote the corresponding Jacobian matrix. Since $J^\varepsilon \to A^CD$ as $\varepsilon \to 0$ and, in general, the spectrum of a matrix depends continuously on its entries, the matrix $J^\varepsilon$ has an eigenvalue with positive real part for $\varepsilon > 0$ small enough. That is, the complex-balanced equilibrium $x^*$ is not linearly stable.

Finally, we prove the implication in the left column (by contradiction). Assume there exists a cycle $C$ in $G$ and a matrix $D \in \mathcal{D}_+$ such that for all $P \in \mathcal{D}_+$ there exists a $w \in S^C$ with

Clearly, the term $v^ \mathsf{T}(PB+B^ \mathsf{T} P)v$ is continuous (in all arguments $v \in S$, $P \in \mathcal{D}_+$, $B \in {\mathbb R}^{n \times n}$), and hence the map $g$, defined by

is also continuous, since maximum and minimum are taken over compact sets. Since $S^C$ is a subspace of $S$ (and hence $w \in S$), the inequality (4) implies $g(A^CD) > 0$.

As above, let $x^*\in {\mathbb R}^n_+$ be such that $D = {\rm{diag}}((x^*)^{-1})$. Further, let $k^\varepsilon \in {\mathbb R}^E_+$ (with $\varepsilon > 0$) be a family of rate constants as in Lemma 4.5, and let $J^\varepsilon$ denote the corresponding Jacobian matrix. Since, by assumption, all complex-balanced equilibria are diagonally stable, there exists $P^\varepsilon \in \mathcal{D}_+$ such that $P^\varepsilon J^\varepsilon+(J^\varepsilon)^ \mathsf{T} P^\varepsilon < 0$ on $S$. As a consequence, $g(J^\varepsilon) < 0$ for all $\varepsilon > 0$, and

a contradiction.

5.   Characterization of D-stability and diagonal stability

For both D-stability and diagonal stability, an explicit characterization is available only up to dimension three. For arbitrary dimension, we recall the following necessary conditions (see e.g. [25,Section 2]).

Proposition 14. Let $A \in {\mathbb R}^{n\times n}$. The following statements hold.

$(i)$ If $A$ is D-stable, then $A$ is a $P_0^+$-matrix (that is, all its signed principal minors are non-negative and at least one of each order is positive).

$(ii)$ If $A$ is diagonally stable, then $A$ is a $P$-matrix (that is, all its signed principal minors are positive).

It is relatively easy to check that, for dimension two, equivalence holds in the previous result.

Proposition 15. Let $A \in {\mathbb R}^{2 \times 2}$. The following statements hold.

$(i)$ The matrix $A$ is D-stable if and only if it is a $P_0^+$-matrix.

$(ii)$ The matrix $A$ is diagonally stable if and only if it is a $P$-matrix.

For dimension three, we have the following characterizations (see ^[28] and [25,Theorem 4(c)]). Thereby, for $A \in {\mathbb R}^{n \times n}$, let $M_{ij}$ denote the principal minor (of order two) $a_{ii}a_{jj}-a_{ij}a_{ji}$.

Proposition 16. Let $A \in {\mathbb R}^{3 \times 3}$. The following statements hold.

$(i)$ The matrix $A$ is D-stable if and only if

$(a)$ $A$ is a $P_0^+$-matrix and

$(b)$ the three pairs $(-a_{11}, M_{23})$, $(-a_{22}, M_{13})$, $(-a_{33}, M_{12})$ dominate $-\det A$, that is,

with equality implying that at least one of the three pairs has exactly one member equal to zero.

(ii) The matrix $A$ is diagonally stable if and only if

$(a)$ $A$ is a $P$-matrix and

$(b)$ there exists $y\in {\mathbb R}$ such that

where $b_1 = a_{12}a_{23}-a_{22}a_{13}$ and $b_2 = a_{21}a_{32}-a_{22}a_{31}$.

Finally, as new results, we characterize D-(semi)stability on a linear subspace with dimension one or two.

Proposition 17. Let $A \in {\mathbb R}^{n \times n}$ and $S \subseteq {\mathbb R}^n$ be a linear subspace with ${\rm{im}} A \subseteq S$ and $\dim S = 1$. Then, $A$ is D-semistable on $S$ if and only if $a_{ii} \leq 0$ for all $i$.

Proof. For $D \in \mathcal{D}_+$, the characteristic polynomial of $AD$ is

Thus, $A$ is D-semistable on $S$ if and only if $b\geq 0$ for all $d_1, \ldots, d_n > 0$ which is equivalent to $a_{ii} \leq 0$ for all $i$.

Proposition 18. Let $A \in {\mathbb R}^{n \times n}$ and $S \subseteq {\mathbb R}^n$ be a linear subspace with ${\rm{im}} A \subseteq S$ and $\dim S = 2$. Then, $A$ is D-stable on $S$ if and only if

$(a)$ $a_{ii} \leq 0$ for all $i$ and $a_{ii} < 0$ for some $i$ and

$(b)$ $M_{ij} \geq 0$ for all $i\neq j$ and $M_{ij} > 0$ for some $i\neq j$.

Proof. For $D \in \mathcal{D}_+$, the characteristic polynomial of $AD$ is

with

Thus, $A$ is D-stable on $S$ if and only if $b > 0$ and $c > 0$ for all $d_1, \ldots, d_n > 0$ which is equivalent to $(a)$ and $(b)$.

6.   Examples

To illustrate our main results, Theorems 4.4 and 4.6, we present a series of examples. In particular, we consider the following networks, leading to special ODE systems:

● irreversible three-cycle (dimension two) $\to$ trinomial ODE system

● irreversible three-cycle (dimension three, stoichiometric subspace of dimension two) $\to$ binomial ODE system

● irreversible four-cycle (dimension three) $\to$ binomial ODE system

● reversible chain (arbitrary dimension)

● S-system (arbitrary dimension) $\to$ binomial ODE system

For the two- and three-dimensional examples, we use the characterizations of D-stability and diagonal stability given in Propositions 5.2, 5.3, and 5.5. For the examples with arbitrary dimension, we use Proposition 17.

Example 19. Consider the generalized mass-action system

and the resulting ODE system

The matrix $YA_{k = 1}\widetilde Y^ \mathsf{T}$ is given by

and

Now, assume $d_{ab} \neq 0$, that is, $\dim S = 2$. By Theorem 11 and Proposition 5.2(i), the unique complex-balanced equilibrium is linearly stable for all $k$ if and only if

and the latter two are not both zero.

Example 20. Consider the generalized mass-action system

and the resulting ODE system

Clearly, $\dim S = 2$. The matrix $YA_{k = 1}\widetilde Y^ \mathsf{T}$ is given by

By Theorem 11 and Proposition 18, all complex-balanced equilibria are linearly stable if and only if

with at least one of the three inequalities satisfied strictly and

with at least one of the three inequalities satisfied strictly.

Example 21. Consider the generalized mass-action system

and the resulting ODE system

Clearly, $\dim S = 3$. The matrix $YA_{k = 1}\widetilde Y^ \mathsf{T}$ is given by

For simplicity, we consider particular kinetic orders

and the resulting matrix

In the following table, we choose particular values for the kinetic orders $\alpha$, $\beta$, and $\gamma$, leading to different situations regarding stability, D-stability, diagonal stability, and being a $P_0^+$-matrix.

In the first and second case, for all rate constants, there exists a complex-balanced equilibrium which is linearly stable by Theorem 11 and Proposition 16(i) and unique by Theorem 10. In the third and fourth case, the unique complex-balanced equilibrium is linearly stable for some rate constants, but not for all.

Example 22. Consider the generalized chemical reaction network

where the graph $G$ is a reversible chain. In particular, all cycles of $G$ correspond to reversible reactions. For every cycle $i \rightleftarrows i+1$ (denoted by $C$), we have

a dyadic product. By Proposition 17, this (rank-one) matrix is D-semistable on $S^C = {\rm{im}} (y(i+1)-y(i))$ if and only if all diagonal entries are non-positive. By Theorem 13, if there is a cycle $i \rightleftarrows i+1$ and an index (species) $s \in \{1, \ldots, n\}$ with

then there is a complex-balanced equilibrium (for some rate constants) that is not linearly stable.

Example 23. For the real matrices $G, H \in {\mathbb R}^{n \times n}$ and the positive vectors $\alpha, \beta \in {\mathbb R}^n_+$, the ODE system

on ${\mathbb R}^n_+$ is called an S-system ^[2,3].

This ODE is associated to the generalized mass-action system

For every cycle $C$ (corresponding to a reversible reaction $\mathsf{0} \rightleftarrows \mathsf{X}_i$), the matrix $YA_{k = 1}^C\widetilde{Y}^ \mathsf{T}$ has $(g_{\cdot i}-h_{\cdot i})^ \mathsf{T}$ as its $i$-th row, and all other rows are zero. By Proposition 17, this (rank-one) matrix is D-semistable on $S^C = {\rm{im}} (e_i)$ if and only if the diagonal entry $g_{ii}-h_{ii}$ is non-positive. By Theorem 13, if there is a reversible reaction $\mathsf{0} \rightleftarrows \mathsf{X}_i$ such that $g_{ii} > h_{ii}$, then there is a complex-balanced equilibrium (for some rate constants) that is not linearly stable.

Acknowledgments

We thank three anonymous referees for helpful comments. BB and SM were supported by the Austrian Science Fund (FWF), project 28406. GR was supported by the FWF, project 27229.

Conflict of interest

The authors declare there is no conflicts of interest.