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

It is well known that, for mass-action systems, complex-balanced equilibria are asymptotically stable. For generalized mass-action systems, even if there exists a unique complex-balanced equilibrium (in every stoichiometric class and for all rate constants), it need not be stable. We first discuss several notions of matrix stability (on a linear subspace) such as D-stability and diagonal stability, and then we apply our abstract results to complex-balanced equilibria of generalized mass-action systems. In particular, we show that linear stability (on the stoichiometric subspace and for all rate constants) implies uniqueness. For cyclic networks, we characterize linear stability (in terms of D-stability of the Jacobian matrix); and for weakly reversible networks, we give necessary conditions for linear stability (in terms of D-semistability of the Jacobian matrices of all cycles in the network). Moreover, we show that, for classical mass-action systems, complex-balanced equilibria are not just asymptotically stable, but even diagonally stable (and hence linearly stable). Finally, we recall and extend characterizations of D-stability and diagonal stability for matrices of dimension up to three, and we illustrate our results by examples of irreversible cycles (of dimension up to three) and of reversible chains and S-systems (of arbitrary dimension).


Introduction
In their foundational paper from 1972, Horn and Jackson considered chemical reaction networks (CRNs) with mass-action kinetics [13]. In particular, they proved that complexbalanced equilibria are asymptotically stable (for all rate constants), by using an entropylike Lyapunov function. Moreover, they 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 [17,18].
For the resulting generalized mass-action systems, existence and uniqueness of complexbalanced equilibria (in every stoichiometric class and for all rate constants) is well understood [17,16], 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 [3,4,5,2].
Since it is hard to find Lyapunov functions for generalized mass-action systems, we 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 [14,Ch. 2], [12,Ch. 26], and the references therein. 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 [9,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 coeffcients; see [16,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: 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.

Generalized mass-action systems
A generalized chemical reaction network (G, y,ỹ) is based on a directed graph G = (V, E) without self-loops; every vertex i ∈ V = {1, . . . , m} is labeled with a (stoichiometric) complex y(i) ∈ R n ≥0 , and every source vertex i ∈ V s ⊆ V is labeled with a kinetic-order complexỹ(i) ∈ R n . If every component of G is strongly connected, we call G and (G, y,ỹ) weakly reversible.
A generalized mass-action system (G k , y,ỹ) is a generalized chemical reaction network (G, y,ỹ) together with edge labels k ∈ R E + , resulting in the labeled digraph G k . Every edge (i → i ′ ) ∈ E, representing the chemical reaction y(i) → y(i ′ ), is labeled with a rate constant k i→i ′ > 0.
The ODE system for the species concentrations x ∈ R n + , associated with the generalized mass-action system (G k , 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 ∈ R n + and y ∈ R n , we define . . , m. The right-hand-side of the ODE system (1) can be decomposed into stoichiometric, graphical, and kinetic-order contributions, where Y, Y ∈ R n×V are the matrices of stoichiometric and kinetic-order complexes, I E , I s E ∈ R V ×E are the incidence and source matrices 1 of the digraph G, and A k = I E diag(k)(I s E ) T ∈ R V ×V is the Laplacian matrix of the labeled digraph G k . (We note that columns of Y corresponding to non-source vertices can be chosen arbitrarily.) Clearly, the change over time lies in the stoichiometric subspace Equivalently, trajectories are confined to cosets of S, that is, x(t) ∈ x(0) + S. For x ′ ∈ R n + , the set (x ′ + S) ∩ R n + is called a stoichiometric class. Analogously, we introduce the kinetic-order subspace A positive equilibrium x ∈ R n + of the ODE system (2) that fulfills A k x 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,ỹ), we study whether, for all rate constants k ∈ R E + (such that Z k = ∅), all complex-balanced equilibria x * ∈ 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 [18]). 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 [17]). Proposition 1. Consider a weakly reversible generalized chemical reaction network, and let x * ∈ R n + . Then, there exist rate constants k ∈ R E + such that x * ∈ Z k . For surveys on the uniqueness of equilibria and related injectivity results, see [15,1,10].

Notions of matrix stability
Let S be a linear subspace of R n . For an ODE systemẋ = f (x) with f : R n → R n and im f ⊆ S, cosets of S are forward invariant. Hence, given an equilibrium x ∈ 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) = ∂f ∂x ∈ R n×n , more precisely, the linear map J(x)| S : S → 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 systeṁ y = J(x)y on S; see e.g. [20,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 [ Proposition 3. Let A ∈ R n×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 T > 0 with P A + A T P ≤ 0, then the origin is Lyapunov stable for the linear ODEẋ = Ax (by the Lyapunov function x → x 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 im A ⊆ S is stable on S (respectively, semistable on S) if all eigenvalues of the linear map A| S : S → S have negative (respectively, non-positive) real part. We extend Lyapunov's Theorem to stability on a linear subspace.
Proposition 4. Let A ∈ R n×n and S ⊆ R n be a linear subspace with im A ⊆ 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 = im B with B ∈ R n×s and B T B = I ∈ R s×s . Then, for every x ∈ S, there exists a unique y ∈ R s such that x = By. Further, y = B T x and, usingẋ = Ax (on S), we haveẏ = B T ABy (on R s ). Thus, A is stable on S if and only if B T AB ∈ R s×s is stable on R s . By Theorem 3, the latter is equivalent to the existence of Q ∈ R s×s with Q = Q T > 0 such that y T (QB T AB + B T A T BQ)y < 0 for all 0 = y ∈ R s . This is further equivalent to the existence of P = BQB T ∈ R n×n with P = P T > 0 on S such that Finally, let P = P T > 0 on S and P A + A T P ≤ 0 on S. With B as above, let Q = B T P B ∈ R s×s and note that A = BB T A (because x = BB T x for all x ∈ S). Then Q = Q T > 0 and QB T AB + B T A T BQ ≤ 0. By Proposition 3, B T AB 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 D + ⊆ R n×n denote the set of diagonal matrices with positive diagonal. A matrix A ∈ R n×n is diagonally stable (respectively, diagonally semistable) if there exists P ∈ D + such that P A + A T P < 0 (respectively, P A + A T P ≤ 0). We note that diagonal stability is also known as Lyapunov diagonal stability or Volterra-Lyapunov stability. A matrix A ∈ R n×n is D-stable (respectively, D-semistable) if AD is stable (respectively, semistable) for all D ∈ D + . The following result summarizes the relations between these notions.
Proposition 5. Let A ∈ R n×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 ∈ D + be such that P A + A T P < 0 (respectively, P A + A T P ≤ 0). Then, for any D ∈ D + , we have (DP )(AD) + (AD) T (DP ) = D(P A + A T P )D < 0 (respectively, ≤ 0), and by Proposition 3, 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 ∈ D + , then this holds for D = I.) Let S be a linear subspace. We say that for all D ∈ D + . Finally, we introduce an even stronger notion. We say that A ∈ R n×n with im A ⊆ S is diagonally D-stable on S (respectively, diagonally D-semistable on S) if, for all D ∈ D + , there exists P ∈ D + such that P AD + DA T P < 0 on S (respectively, P AD+DA T P ≤ 0 on S). We note that diagonal D-stability on S trivially implies diagonal stability on S, and the two notions agree for S = R n , see also [8, p. 257]. Moreover, we have the following relations.
Proposition 6. Let A ∈ R n×n and S ⊆ R n be a linear subspace with im A ⊆ 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 4.
Finally, the implications from the second row to the third row are trivial.

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

Linear stability for mass-action systems
We start by showing that, for classical mass-action systems (whereỹ = y), complexbalanced 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).
In fact, the author proved diagonal stability, without noting it and using an ad-hoc inequality. Here, we even show diagonal D-stability, thereby using the negative semidefiniteness of the Laplacian matrix of an undirected graph.
Theorem 7. 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 ∈ R E + and x * ∈ Z k . We show that the corresponding Jacobian matrix J = J(x * ) is diagonally D-stable on S. By Proposition 6, all other conclusions follow.
Let D ∈ D + . We show that there exists P = diag((x * ) −1 )D ∈ D + such that H = P JD+DJ T P < 0 on S. Let k * ∈ R E + be defined by k Now, let Ak = A k * + A T k * and hence H = P Y Ak Y T P . The symmetric matrix Ak is the Laplacian matrix of a labeled undirected graph Gk with G = (V, E) andk ∈ R E + . In particular, Ak = −I E ′ diag(k) I T E ′ for some directed version E ′ of the undirected edges E and the corresponding incidence matrix I E ′ (with im I E ′ = im I E ); see e.g. [6]. Hence, , and hence P v ∈ S ⊥ . Clearly, v 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 P J + J T P < 0 on S. This property implies the existence of a function with separated variables L(x) = n i=1 L i (x i ) that serves as a Lyapunov function for showing the asymptotic stability of the complex-balanced equilibrium x * ∈ R n + (in its stoichiometric class).
As an example, consider a scaled version of the classical entropy-like Lyapunov function where P = diag(p 1 , . . . , p n ), and the function . Then, the gradient of g at x * vanishes, and the Hessian matrix of g at x * is H = P J + J 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).

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 8. Consider a weakly reversible generalized mass-action system, and let x * ∈ Z k . Then, ker(A k diag((x * ) Y )) = ker I T E .
Proof. The dimension of ker I T E 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.
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. [18,Section 4] and the references therein.
By definition, x * ∈ Z k if A k (x * ) 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 diag((x * ) Y ). (As an example, consider a strongly connected graph G, and let 1 V be the characteristic vector of V . Then, Theorem 9. 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 ′ ∈ R E + , there exists a stoichiometric class with at least two complex-balanced equilibria. Then, by Proposition 2, there exist vectors 0 = u ∈ S ⊥ and 0 = v ∈ S with sign(u) = sign(v). Now, S = im( Y I E ), S ⊥ = ker(I T E Y T ), and hence I T E Y T u = 0. Let x * ∈ R n + be such that u = diag((x * ) −1 )v, and let k ∈ R E + be such that x * ∈ Z k , see Proposition 1. Then, , the latter by Lemma 8. Using (3) for the Jacobian matrix J, we have Hence, for all P = P T > 0, we have v T (P J + J T P )v = 0, and by Proposition 4, J is not stable on S.
As can be demonstrated by the generalized reaction network X(X) ⇋ Y(0), linear stability implies uniqueness, but not existence in every stoichiometric class.

The network is a cycle
In the following result, we characterize diagonal stability and linear stability of complexbalanced equilibria, provided that the network is a cycle.  A is diagonally D-stable on S.
⇒ A is D-stable on S.
Proof. The implication in the first (respectively, second) row follows immediately from Proposition 4 (respectively, Proposition 6).
To prove the two equivalences, let k ∈ R E + and x * ∈ Z k . Since G = (V, E) is a cycle, the quantity c = k i→i ′ (x * )ỹ (i) is the same for all (i → i ′ ) ∈ E. In particular, A k diag((x * ) Y ) = c A k=1 . Using (3) for the Jacobian matrix J, we have . Finally, recall that every x * ∈ R n + is a complex-balanced equilibrium for some k ∈ R E + , see Proposition 1. Hence, every D ∈ D + is of the form D = diag((x * ) −1 ) for some k, and the equivalences follow.

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 11. Consider a weakly reversible generalized chemical reaction network, let C be a cycle of the network, and let x * ∈ R n + . Then, there exists a family of rate constants is the Laplacian matrix of the cycle C with all rate constants set to 1.
Proof. For a cycle C ′ , let k C ′ ∈ R E ≥0 be defined by For ε > 0, let k ε = k C + ε C ′ =C k C ′ , where the summation is over all cycles C ′ of G, except for the given cycle C. Then, x * ∈ Z k ε and Using (3) for the Jacobian matrix J ε , we obtain the desired limit as ε → 0.
Theorem 12. Consider a weakly reversible generalized chemical reaction network, and, for a cycle C of the network, let A C = Y A C k=1 Y T and S C be the corresponding stoichiometric subspace. The following implications hold: For all rate constants, complex-balanced equilibria are diagonally stable (in their stoichiometric classes).

⇒
For all rate constants, complex-balanced equilibria are linearly stable (in their stoichiometric classes).

⇓ ⇓
For all cycles C, A C is diagonally D-semistable on S C . ⇒ For all cycles C, A C is D-semistable on S C .
Proof. The implications in the first (respectively, second) row follows immediately from Proposition 4 (respectively, Proposition 6).
Next, we prove the implication in the right column. Assume there exists a cycle C in G and a matrix D ∈ D + such that A C D has an eigenvalue with positive real part. Let x * ∈ R n + be such that D = diag((x * ) −1 ). Further, let k ε ∈ R E + (with ε > 0) be a family of rate constants as in Lemma 11, and let J ε denote the corresponding Jacobian matrix. Since J ε → A C D as ε → 0 and, in general, the spectrum of a matrix depends continuously on its entries, the matrix J ε has an eigenvalue with positive real part for ε > 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 ∈ D + such that for all P ∈ D + there exists a w ∈ S C with Clearly, the term v T (P B + B T P )v is continuous (in all arguments v ∈ S, P ∈ D + , B ∈ R n×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 ∈ S), the inequality (4) implies g(A C D) > 0.
As above, let x * ∈ R n + be such that D = diag((x * ) −1 ). Further, let k ε ∈ R E + (with ε > 0) be a family of rate constants as in Lemma 11, and let J ε denote the corresponding Jacobian matrix. Since, by assumption, all complex-balanced equilibria are diagonally stable, there exists P ε ∈ D + such that P ε J ε + (J ε ) T P ε < 0 on S. As a consequence, g(J ε ) < 0 for all ε > 0, and

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. [8, Section 2]).
Proposition 13. Let A ∈ R n×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.
(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 [7] and [8, Theorem 4(c)]). Thereby, for A ∈ R n×n , let M ij denote the principal minor (of order two) a ii a jj − a ij a ji . Proof. For D ∈ D + , the characteristic polynomial of AD is Thus, A is D-stable on S if and only if b > 0 and c > 0 for all d 1 , . . . , d n > 0 which is equivalent to (a) and (b).

Examples
To illustrate our main results, Theorems 10 and 12, we present a series of examples. In particular, we consider the following networks, leading to special ODE systems: • irreversible three-cycle (dimension two) → trinomial ODE system • irreversible three-cycle (dimension three, stoichiometric subspace of dimension two) → binomial ODE system • irreversible four-cycle (dimension three) → binomial ODE system • reversible chain (arbitrary dimension) Example 18. Consider the generalized mass-action system k 23 k 31 and the resulting ODE system dx dt = k 12 (a 2 − a 1 )x α1 y β1 + k 23 (a 3 − a 2 )x α2 y β2 + k 31 (a 1 − a 3 )x α3 y β3 , The matrix Y A k=1 Y T is given by Now, assume d ab = 0, that is, dim S = 2. By Theorem 10 and Proposition 14(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 19. Consider the generalized mass-action system and the resulting ODE system By Theorem 10 and Proposition 17, all the 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 20. Consider the generalized mass-action system For simplicity, we consider particular kinetic orders 0 (γZ) In the following table, we choose particular values for the kinetic orders α, β, and γ, leading to different situations regarding stability, D-stability, diagonal stability, and being a P + 0 -matrix. α β γ Y A k=1 Y T 0 0 0 diagonally stable 5 0 −3 D-stable, but not diagonally stable 3 4 −4 stable P + 0 -matrix, but not D-stable 2 −2 1 stable, but not P + 0 -matrix 0 −2 −3 unstable 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 10 and Proposition 15(i) and unique by Theorem 9. In the third and fourth case, the unique complex-balanced equilibrium is linearly stable for some rate constants, but not for all.
x n = α n x gn1 1 · · · x gnn n − β n x hn1 1 · · · x hnn n on R n + is called an S-system [19,21]. This ODE is associated to the generalized massaction system 0 (g 11 X 1 + · · · + g 1n X n ) X 1 (h 11 X 1 + · · · + h 1n X n ) α 1 β 1 . . . 0 (g n1 X 1 + · · · + g nn X n ) X n (h n1 X 1 + · · · + h nn X n ) α n β n For every cycle C (corresponding to a reversible reaction 0 ⇄ X i ), the matrix Y A C k=1 Y T has (g ·i − h ·i ) T as its i-th row, and all other rows are zero. By Proposition 16, this (rankone) matrix is D-semistable on S C = im(e i ) if and only if the diagonal entry g ii − h ii is non-positive. By Theorem 12, if there is a reversible reaction 0 ⇄ X i such that g ii > h ii , then there is a complex-balanced equilibrium (for some rate constants) that is not linearly stable.