A New Decomposition of the Graph Laplacian and the Binomial Structure of Mass-Action Systems

We provide a new decomposition of the Laplacian matrix (for labeled directed graphs with strongly connected components), involving an invertible core matrix, the vector of tree constants, and the incidence matrix of an auxiliary graph, representing an order on the vertices. Depending on the particular order, the core matrix has additional properties. Our results are graph-theoretic/algebraic in nature. As a first application, we further clarify the binomial structure of (weakly reversible) mass-action systems, arising from chemical reaction networks. Second, we extend a classical result by Horn and Jackson on the asymptotic stability of special steady states (complex-balanced equilibria). Here, the new decomposition of the graph Laplacian allows us to consider regions in the positive orthant with given monomial evaluation orders (and corresponding polyhedral cones in logarithmic coordinates). As it turns out, all dynamical systems are asymptotically stable that can be embedded in certain binomial differential inclusions. In particular, this holds for complex-balanced mass-action systems, and hence, we also obtain a polyhedral-geometry proof of the classical result.


Introduction
The Laplacian matrix (or graph Laplacian) is a matrix representation of a graph.It can be seen as a discrete version of the Laplace operator defined on graphs.On the one hand, the Laplacian matrix of an undirected graph, its spectrum, and its eigendecomposition have a variety of applications ranging from organic chemistry to signal processing and machine learning [24,22,31,3].On the other hand, labeled, directed graphs underlie dynamical systems ranging from continuous-time Markov processes (linear stochastic models) [23] to mass-action systems (non-linear deterministic models of chemical reaction networks) [20].
In the linear setting, the vertices V of a simple digraph G = (V, E) represent states, and the edges E represent transitions.Moreover, edge labels k represent transition rate constants.The dynamical system for a state variable ψ is given by dψ where A k is the Laplacian matrix of the labeled digraph G k = (V, E, k).That is, (A k ) i,j = k j→i if there is a transition (j → i) ∈ E, (A k ) i,i = − (i→j)∈E k i→j , and (A k ) i,j = 0 otherwise.(As in chemical reaction network theory, we use the letter A for the graph Laplacian and indicate its dependence on the edge labels k by a subscript.)The linear system can be called "Laplacian dynamics", it is equivalent to the stochastic master equation, and it is studied in applications ranging from biochemistry to systems biology [17,23].
In the nonlinear setting, the dynamical system for the species concentrations x is given by dx All notation is defined at the end of this introduction, and mass-action systems are introduced in Section 3. Here, we motivate Eqn. ( 2) in an informal way.As an example, we consider the chemical reaction 1X 1 + 1X 2 → X 3 with "stoichiometric" coefficients equal to 1.Under the assumption of massaction kinetics, its rate is given by v = k (x 1 ) 1 (x 2 ) 1 , where k > 0 is the rate constant, and x 1 , x 2 ≥ 0 are the concentrations of the species X 1 , X 2 .More abstractly, we can write the reaction as y → y ′ with (educt and product) "complexes" y = (1, 1, 0, 0, . ..)T and y ′ = (0, 0, 1, 0, . ..)T , and we can write its rate as v = k x y with the monomial x y := j (x j ) yj = (x 1 ) 1 (x 2 ) 1 (x 3 ) 0 (x 4 ) 0 • • • in the species concentrations x = (x 1 , x 2 , x 3 , x 4 , . ..)T .In a network, an individual reaction y → y ′ contributes the summand k x y (y ′ − y) to the dynamical system for x, where the reaction vector y ′ − y captures the consumption of educts y and the formation of products y ′ .For the example reaction, x y = x 1 x 2 (as stated above) and y ′ − y = (−1, −1, 1, 0, . ..)T .Now, we can introduce a mass-action system as a simple digraph G = (V, E), a map y (assigning complexes to vertices), and edge labels k.In particular, every edge (i → i ′ ) ∈ E defines a reaction y(i) → y(i ′ ) with rate constant k i→i ′ .Hence, the associated dynamical system dx dt = (i→i ′ )∈E k i→i ′ x y(i) (y(i ′ )−y(i)) involves a sum over all edges, and every summand is a product of a reaction rate and a reaction vector.Using the Laplacian matrix A k , the right-hand-side can be decomposed as shown in Eqn.(2).The matrix Y collects the complexes y(i) for i ∈ V , and the vector of monomials x Y is defined via (x Y ) i = x y(i) .Altogether, the dynamical system is polynomial.It is determined by the complex matrix Y (by stoichiometry) as well as by the Laplacian matrix A k (by the graph), and chemical reaction network theory studies the interplay of these two matrices to understand dynamics and steady states of mass-action systems, starting from the foundational 1972 papers [20,18,12] until today.
A steady state x > 0 with A k x Y = 0 is called a positive complex-balanced equilibrium (CBE), also known as vertex-balanced steady state.Indeed, at a CBE, the sum of all "flows" k i→i ′ x y(i) from vertex i/complex y(i) equals the sum of all k i ′ →i x y(i ′ ) to the latter.As shown by Horn [18] and Horn & Jackson [20] in 1972, the existence of a CBE has three important consequences: the components of the graph are strongly connected (the network is "weakly reversible"); all equilibria are complex-balanced and asymptotically stable; and there is a unique equilibrium in every dynamically invariant affine subspace ("stoichiometric compatibility class").More technically, complex-balanced equilibria are given by binomial equations and have a monomial parametrization.
For symmetric digraphs ("reversible" networks), detailed-balanced equilibria are given by binomial equations (by definition).Moreover, the polynomial dynamical system is a sum of binomials.(Just note that every reversible reaction y ⇄ y ′ contributes the summand (k y→y ′ x y −k y ′ →y x y ′ ) (y ′ −y) to the dynamical system for x.)We show that this also holds for weakly reversible networks.To this end, we provide a new decomposition of the graph Laplacian, involving an invertible core matrix, based on an order on the vertices.Further, we extend the classical result by Horn and Jackson on the asymptotic stability of complex-balanced equilibria.In addition to a Lyapunov function (as in classical proofs), we consider regions in the positive orthant with given monomial evaluation orders (and corresponding polyhedral cones in logarithmic coordinates).As it turns out, all dynamical systems are asymptotically stable that can be embedded in certain binomial differential inclusions.In particular, this holds for complex-balanced mass-action systems, and hence we also obtain a polyhedral-geometry proof of the classical result.
Organization of the work.In Section 2, we provide a new decomposition of the graph Laplacian (for labeled directed graphs with strongly connected components), involving an invertible core matrix, based on an order on the vertices.Depending on the particular order, the core matrix has additional properties.
In Section 3, we apply the graph-theoretic/algebraic results to mass-action systems.In Subsection 3.1, we demonstrate their binomial structure, and in 3.2, we introduce monomial evaluation orders and corresponding geometric objects (polyhedra and polyhedral cones).In Subsection 3.3, we embed complex-balanced mass-action systems in binomial differential inclusions and show that all equilibria of the latter are asymptotically stable, and in 3.4, we discuss our results.
In Appendix A, we provide explicit formulas for the vector of tree constants and the Laplacian matrix, using cycle decomposition.In Appendix B, we state auxiliary results used in the new decomposition of the graph Laplacian.In Appendix C, we give another proof of the asymptotic stability of complexbalanced equilibria (and the non-existence of other steady states) without using differential inclusions.
Notation.We denote the positive real numbers by R > and the nonnegative real numbers by R ≥ .Throughout the work, we use index notation: for a finite index set I, we write R I for the real vector space of vectors x = (x i ) i∈I with x i ∈ R, and analogously we write R I ≥ and R I > .(For I = {1, . . ., n}, we have the standard case R I = R n .)We write x > 0 for x ∈ R I > and x ≥ 0 for x ∈ R I ≥ .
For vectors x, y ∈ R I , we denote their scalar product by x • y ∈ R and their componentwise (Hadamard) product by x • y ∈ R I .For x ∈ R I > , y ∈ R I , we define the (generalized) monomial x y = i∈I (x i ) yi ∈ R > , and for x ∈ R I > , Y ∈ R I×J , we define the vector of monomials x Y ∈ R J > via (x Y ) j = x y(j) , where y(j) is the column of Y with index j ∈ J.

The graph Laplacian
In the following, we assume that the components of a digraph are strongly connected.For the simplicity of the presentation, we first consider one strongly connected component separately.

One component
We consider a strongly connected, simple, directed graph G = (V, E) with a finite set of vertices V = {1, . . ., m} and a set of edges E ⊆ V × V .Further, we consider positive edge labels k ∈ R E > and the resulting labeled digraph where I E ∈ R V ×E is the incidence matrix and I E,s ∈ R V ×E is the "source matrix".Explicitly, otherwise, and This definition is used in dynamical systems.For example, k is the vector of transition rate constants in the continuous-time, linear process dψ dt = A k ψ (with ψ ∈ R V ≥ and i∈V ψ i = 1).In other fields, the Laplacian matrix is defined as , where 1 ∈ R V is the vector with all entries equal to one.Further, ker with a positive vector K k ∈ R V > (depending on the rate constants).The entries of K k (the tree constants) can be given explicitly in terms of k, where T i is the set of directed spanning trees of G rooted at vertex i ∈ V (and directed towards the root).For a minimal proof of Eqn.(3), see [21, Lemma 1] or Appendix A. We note that the explicit formula is not crucial for our analysis.Finally, the tree constants K k correspond to minors of the matrix −A k which is the content of the matrix-tree theorem (for labeled, directed graphs) [34,Theorem 3.6].
Clearly, the matrix has positive diagonal entries and nonpositive off-diagonal entries.Most importantly, it has zero row and column sums: Indeed, 1T A k diag(K k ) = 0, and also As a consequence, the matrix is diagonally dominant.
The entries of A k diag(K k ) can be given explicitly in terms of k.For a derivation of this formula and a discussion of the Birkhoff/von Neumann Theorem [5,35], see Appendix A. Again, we note that the explicit formula is not crucial for our analysis.
Example.Throughout this section, we consider the labeled directed graph Most importantly, we introduce an auxiliary connected directed graph G E = (V, E) with the same set of vertices V as in G = (V, E), but with an arbitrary set of edges In particular, it has no cycles.Further, G E need not be a subgraph of G nor be directed towards a root.The corresponding incidence matrix Note that the definitions of the incidence matrices I E and I E agree formally.
(Just the sets of edges E and E differ.)Clearly, ker I E = {0} and ker I T E = im 1.
Proposition 1.Let G k = (V, E, k) be a strongly connected, labeled, simple digraph and G E = (V, E) be an auxiliary digraph.Then, there exists a unique invertible matrix A k,E ∈ R E×E , called the core matrix of the graph Laplacian, such that Proof.Since G is strongly connected, for a unique matrix B k,E ∈ R V ×E , where uniqueness follows from ker I E = {0}.
For the same reason, we have Finally, we obtain For an auxiliary digraph In the following two results, we assume G E to be either of the form Proposition 2. Let G k = (V, E, k) be a strongly connected, labeled, simple digraph, and let G E = (V, E) be an auxiliary digraph that is a chain graph.Then A k,E ∈ R E×E , the core matrix of the graph Laplacian, is non-negative with positive diagonal.
Proof.Let G E = (V, E) be the chain graph It induces a natural order on the set of vertices V (and on the set of edges E).
For i, j ∈ V , we write i ≤ j if i = j or i → . . .→ j.An "inverse" of the incidence matrix Explicitly, using the order i 1 , i 2 , . . ., i m on V , , and indeed, J E I E = −I, where I ∈ R E×E is the identity matrix.That is, −J E is a generalized left-inverse of I E .Hence, by Proposition 1,

For an arbitrary matrix
Explicitly, (σ) is the sum of all entries in the upper left i × j block of A. Now, recall that the matrix A = −A k diag(K k ) has positive diagonal entries and nonpositive off-diagonal as well as zero row and column sums.Hence, the sum (σ) is nonnegative.Finally, recall that the underlying graph G is strongly connected.If i → i ′ equals j → j ′ , then the sum (σ) is positive, since the corresponding subgraph with vertices {i 1 , i 2 , . . ., i} has incoming and outgoing edges.
Example (continued).In the labeled digraph G k = (V, E, k) introduced above, there are 3 vertices and hence 6 possible chain graphs.For example, for ) be a strongly connected, labeled, simple digraph, and let G E = (V, E) be an auxiliary digraph that is a star graph.Then A k,E ∈ R E×E , the core matrix of the graph Laplacian, is (row and column) diagonally dominant with positive diagonal and non-positive off-diagonal entries.
Explicitly, let Proof.Let G E = (V, E) be the star graph An "inverse" of the incidence matrix Explicitly, using the order i 1 , i 2 , . . ., i m on V , , and indeed, That is, (σ ⋆ ) is the sum of all entries of A except the entries in row i and column j.Now, recall that the matrix A = −A k diag(K k ) has zero row and column sums.Hence, (σ ⋆ ) equals the sum of all entries (which is zero) minus the sums of all entries in row i and column j (which are zero) plus the common entry of row i and column j.That is, As claimed, A k,E equals A =−A k diag(K k ) with row i m and column i m removed.Like A, it has positive diagonal entries and nonpositive off-diagonal entries and is (row and column) diagonally dominant.(However, not all row and column sums are zero.) Example (continued).In the labeled digraph G k = (V, E, k) introduced above, there are 3 vertices and hence 3 possible star graphs.For example, for Remark.In applications to mass-action systems in Section 3, we use chain graphs (rather than star graphs).

Several components
In general, we consider a labeled, simple digraph Accordingly, an auxiliary digraph a chain graph, and analogously for a star graph.Propositions 1, 2, and 3 imply the main result of this section.Theorem 4. Let G k = (V, E, k) be a labeled, simple digraph with strongly connected components, and let G E = (V, E) be an auxiliary digraph.Then, there exists an invertible, block-diagonal matrix A k,E ∈ R E×E , called the core matrix of the graph Laplacian, such that is diagonally dominant with positive diagonal and non-positive off-diagonal entries.

Mass-action systems
We apply the graph-theoretic/algebraic results from the previous section to mass-action systems.We start with a brief summary of fundamental concepts and results.
A chemical reaction network (G, y) is given by a simple directed graph G = (V, E) with a finite set of vertices V = {1, . . ., m} and a set of edges (reactions) E ⊆ V × V together with an injective map y : ) If the components of G (the linkage classes) are strongly connected, then the network is called weakly reversible.
A mass-action system (G k , y) is a chemical reaction network (G, y) where every edge (i (If the network is weakly reversible, then also the mass-action system is called weakly reversible.) The resulting dynamical system for x ∈ R n ≥ (the concentrations of n molecular species) is given by The right-hand side of the ODE can be decomposed as where I E ∈ R V ×E is the incidence matrix, I E,s ∈ R V ×E is the "source matrix", and is the resulting Laplacian matrix of the labeled, simple digraph G k .In the following, we consider the dynamical system in the form The stoichiometric subspace is given by S = im(Y I E ).Clearly, dx dt = f k (x) ∈ S, and hence x(t) ∈ x(0) + S.
then it is a positive complex-balanced equilibrium (CBE), also known as vertexbalanced steady state.
Remark.In the linear setting, the Laplacian matrix captures state transitions on a graph.Let ψ = x Y be the state variable, given by the vector of monomials.
If A k ψ = 0, then transitions are balanced (at every vertex of the graph), and As shown by Horn [18] and Horn & Jackson [20], if there exists a positive CBE (in some stoichiometric class), then 1. the mass-action system is weakly reversible [18, Theorem 3C], 2. the equilibrium is asymptotically stable, and all equilibria are complexbalanced [20, Theorem 6A], and, 3. there exists a unique positive (necessarily complex-balanced) equilibrium in every stoichiometric class [20,Lemma 4B].
In the following remarks, we elaborate on results 1, 2, and 3.
Remark (result 1).Let G be weakly reversible and G E = (V, E) be some auxiliary digraph.By Theorem 4, Given a particular positive CBE x * ∈ R n > , Eqn. ( 6) is equivalent to and further to (y(i ′ ) − y(i)) T ln(x/x * ) = 0 for (i , the set of all positive CBEs is given by the monomial parametrization x = x * • e S ⊥ . Remark (result 2).In Section 3.3, we extend the classical stability result.As it turns out, it holds not only for complex-balanced equilibria of mass-action systems, but for all equilibria of binomial differential inclusions.
In Appendix C, we give another proof for the asymptotic stability of complexbalanced equilibria (and the non-existence of other steady states) without using differential inclusions.

Binomial structure
Given that the network is weakly reversible (the components of the graph G are strongly connected), our main graph-theoretic/algebraic result, Theorem 4, implies that the dynamical system (4) for the mass-action system (G k , y) can be decomposed as where G E = (V, E) is some auxiliary digraph.
Again, we have a closer look at the term That is, the right-hand side of the dynamical system is a sum of binomials.This is obvious for symmetric digraphs (reversible networks); cf.[9, Eqn. ( 14)].
By Theorem 4, it also holds for digraphs with strongly connected components (weakly reversible networks).
In particular, for a complex-balanced equilibrium, not just the right-hand side of ( 7) is zero, but every individual binomial is zero.In this sense, the ODE (7) does not only have binomial steady states (positive complex-balanced equilibria, given by binomial equations), but truly is a binomial dynamical system.

Monomial evaluation orders and corresponding polyhedra/polyhedral cones
Let (G k , y) be a mass-action system based on the labeled, simple digraph G k = (V, E, k) and the map y (the matrix Y ).
For fixed x ∈ R n > , the values of the monomials x y(i) with i ∈ V are ordered (using the order on R).For simplicity, we first consider a connected graph G = (V, E).Obviously, the total order x y(i1) ≤ x y(i2) ≤ . . .≤ x y(im) can be represented by a chain graph, If the order is non-strict (if some monomials have the same value), then the representation is not unique.Analogously, the partial order x y(i1) ≤ x y(im) , x y(i2) ≤ x y(im) , . . ., x y(im−1) ≤ x y(im) can be represented by a star graph, In general, every auxiliary graph G E = (V, E) represents a partial order on the vertices of G and hence on the values of the monomials.
In the following, we will consider monomials with coefficients: , for weakly reversible networks with tree constants K k ∈ R V > , and

Weak reversibility
Let (G k , y) be a weakly reversible mass-action system, and fix x ∈ R n > .We call an order on the entries of x Y K k ∈ R V > that is total within connected components, but does not relate entries in different components, a monomial evaluation order (since the notion monomial order(-ing) has a different meaning in algebra).We represent the order by a chain graph G E = (V, E) and often just by the set of edges E. Explicitly, (i Thereby, the vertices i, i ′ ∈ V are necessarily in the same component.If the order is non-strict, then E is not unique. Analogously, the maximal entries of x Y K k ∈ R V > within connected components are greater or equal than all other entries in the respective components.We represent this order by a star graph G E = (V, E).If there is more than one maximal entry within a component, then E is not unique.
Conversely, fix an auxiliary graph G E = (V, E), for example, a chain graph or a star graph.The subset of R n > with monomial evaluation order represented by E is given by By the monotonicity of the logarithm, with the polyhedron

Complex balancing
If there exists a positive CBE x * ∈ R n > , then the polyhedra become polyhedral cones.
Fix an auxiliary graph G E = (V, E).Using complex balancing (6) for x * , the subset (8) can be written as By the monotonicity of the logarithm, which does not depend on k. (Of course, x * depends on k.)The lineality space of C E does not even depend on E, Obviously, S ⊥ = lineal C E ⊆ C E .For fixed E, there are two possibilities: • C E = S ⊥ .Then, all defining (non-strict) inequalities of C E (and S k,E ) are fulfilled with equality, and S k,E = x * • e S ⊥ equals the set of complexbalanced equilibria.
• C E ⊃ S ⊥ .Then C E and S k,E are full-dimensional, and the monomial evaluation order is strict in the interior of S k,E and non-strict on the boundary (where some monomials have the same value).
In the following study of complex-balanced mass-action systems (and their extension to binomial differential inclusions), we use chain graphs G E , representing monomial evaluation orders.In this setting, a full-dimensional subset S k,E is called a stratum, cf.[32].This term has also been used for partial orders related to the original graph, rather than to an auxiliary graph, cf.[9].
Remark.As stated above, for every x ∈ R n > , there is a (non-unique) E such that x ∈ S k,E .In particular, R n > is a union of strata which intersect only on their boundaries.Correspondingly, R n is a union of polyhedral cones C E .Indeed, by the monotonicity of the logarithm, an order on the entries of ( (within components) is equivalent to an order on the entries of Y T z ∈ R V with z = ln x x * , and the set of pairs of vertices within components, induces an arrangement of central hyperplanes, The central hyperplane arrangement decomposes R n into open polyhedral cones called faces; full dimensional faces are called cells.In our terminology, a cell is the interior of a polyhedral cone C E and hence corresponds to the interior of a stratum S k,E .
x * S k,E The positive orthant is a union of strata corresponding to monomial evaluation orders.In particular, consider the stratum given by the order x y(1) ≤ x y(2) ≤ x y (3) , that is, S k,E with E = {1 → 2, 2 → 3}, bounded by the green and blue lines.The green line specifies x y(1) = x y (2) ; above it, x y(2) > x y (1) , as indicated by the corresponding vertices 2 and 1.The blue line specifies x y(2) = x y(3) ; below it, x y(3) > x y (2) .(The dashed black line specifies x y(1) = x y(3) , which does not bound the particular stratum.)In the interior of S k,E , the order is strict.In logarithmic coordinates z = ln(x/x * ), the stratum corresponds to the polyhedral cone C E .
In general, S k,E = x * • e S ⊥ (equals the set of complex-balanced equilibria) if and only if C E = S ⊥ .In the example, S = R 2 and S ⊥ = {0}.

Binomial differential inclusions
Finally, we extend a classical result by Horn and Jackson from 1972.
Theorem 5 (cf.[20], Theorem 6A).Let (G k , y) be a mass-action system and x * ∈ R n > be a positive CBE of the dynamical system (4).Then, for all x ∈ R n > that are not complex-balanced equilibria.Hence, (i) all positive equilibria are complex-balanced, and (ii) x * is asymptotically stable.
All proofs are based on the entropy-like Lyapunov function L : R n > → R, For T , and hence T f k (x) ≤ 0 with "=" if and only if x = x * , then L(x) is a strict Lyapunov function, and x * is asymptotically stable.
Previous proofs further use inequalities for the exponential function or the logarithm and cycle decomposition of the graph, cf.[20,33,1,15].For a new proof using monomial evaluation orders and corresponding geometric objects (strata and polyhedral cones), see Appendix C.
In the following, we extend the stability result and provide a maximally transparent, polyhedral-geometry proof.First, we relate the dynamics in a given stratum to the corresponding polyhedral cone.

Polar cone
In Proposition 6 below, we use the concept of the polar cone Proposition 6.Let (G k , y) be a complex-balanced mass-action system, G E be a chain graph, and S k,E ⊂ R n > be a stratum.Then, for all x ∈ S k,E that are not positive complex-balanced equilibria, Proof.Let x ∈ S k,E and u ∈ C E .Using the dynamical system (4) and Theorem 4, we have Using S k,E and C E as in Eqns.( 8) and ( 10), we have b ≥ 0 and a ≥ 0.
By Theorem 4, the core matrix of the graph Laplacian, A k,E ∈ R E×E , is nonnegative with positive diagonal.Hence, Now, let (G k , y) be a mass-action system and x * ∈ R n > be a positive CBE of the dynamical system (4).Proposition 6 suggests to introduce a corresponding piece-wise constant binomial differential inclusion as thereby explicitly specifying the set of positive equilibria x * • e S ⊥ .Equivalently, Proposition 6 immediately implies the following result.
Theorem 7. Let (G k , y) be a mass-action system and x * ∈ R n > be a positive CBE of the dynamical system (4).Then, the mass-action system can be embedded in the binomial differential inclusion (12).
Finally, we extend Theorem 5 (from complex-balanced mass-action systems to binomial differential inclusions).Theorem 8. Let x * ∈ R n > be a positive equilibrium of the binomial differential inclusion (12).Then, for all x ∈ R n > that are not positive equilibria and all f ∈ F (ln x x * ).Hence, x * is asymptotically stable.
Proof.Let S k,E ⊂ R n > be a stratum and x ∈ S k,E not be a positive equilibrium.On the one hand, that is, ln x x * lies in C E , but not in the lineality space lineal C E = S ⊥ .On the other hand, and T f < 0, and L(x) is a strict Lyapunov function.
Remark 9.Even if a weakly reversible mass-action system (G k , y) does not admit a complex-balanced equilibrium x * , it can be embedded in a piece-wise constant differential inclusion.Technically, the absence of a CBE x * does not allow to pass from the polyhedron P k,E (with given monomial evaluation order) to the cone C E , cf.Eqns.( 9) and (10).That is, instead of a central hyperplane arrangement (that defines the cones C E ), one considers a non-central hyperplane arrangement (that defines the polyhedra P k,E ).In analogy to Proposition 6, one can show that, for a chain graph G E and a stratum S k,E , it holds that f k (x) ∈ int(rec(P k,E ) pol ), for all x ∈ S k,E .Here, rec(C) denotes the recession cone of a set C.

Discussion
As Horn and Jackson in 1972 [20, Theorem 6A], we have shown that, in massaction systems with a positive complex-balanced equilibrium, every positive equilibrium is complex-balanced and asymptotically stable.For a proof using the new decomposition of the graph Laplacian, monomial evaluation orders, and corresponding geometric objects (strata and polyhedral cones), see Appendix C.
In fact, we have extended the result to binomial differential inclusions (BDIs), introduced in this work.Every positive equilibrium of a BDI is asymptotically stable, see Theorem 8.

Binomial and toric differential inclusions
Given a reaction network (G, y) with graph G = (V, E) and "complex" map y : V → R n ≥ , a BDI depends on the components of the graph (but not on the exact edge set E) and on some positive equilibrium x * (but not explicitly on the rate constants).In fact, it is mainly determined by stoichiometry, namely by pairwise differences of complexes, defining a hyperplane arrangement.In particular, monomial evaluation orders correspond to polyhedral cones (in logarithmic coordinates) and strata (in the original positive variables).More formally, a BDI is given by a hyperplane arrangement (with lineality space S ⊥ ) and a positive equilibrium x * , see Equation (12).Most importantly, complex-balanced massaction systems can be embedded in BDIs.
Recently, toric differential inclusions (TDIs) have been used in a proposed proof [7,8] of the global attractor conjecture [19], stating that complex-balanced equilibria are not just asymptotically, but also globally stable.In fact, TDIs also allow to tackle the persistence and permanence conjectures for (weakly reversible) mass-action systems with (time-)variable rate constants.In the classical setting, rate constants k > 0 are fixed, whereas, in the study of the conjectures mentioned above, rate constants ǫ ≤ k(t) ≤ 1/ǫ may vary over time, but are bounded [1,11].To address this complication, "uncertainty regions" with thickness δ(ǫ) around the boundaries of "regions with definite monomial order" are introduced.On the one hand, BDIs are special cases of TDIs with δ → 0 (modulo a translation of the hyperplane arrangement by log x * ), and also the piece-wise constant differential inclusions mentioned in Remark 9 can be embedded in TDIs (with δ > 0).On the other hand, BDIs allow to consider (the asymptotic stability of) positive equilibria, whereas TDIs capture the dynamics close to the boundary of the positive orthant without being explicit about equilibria.

Generalized mass-action systems
In previous work, we have studied generalized mass-action systems [27,28,25,26,10,6].In order to motivate the setting, we consider the reaction 1X 1 + 1X 2 → X 3 with "stoichiometric" coefficients equal to 1.Under the assumption of generalized mass-action kinetics, its rate is given by v = k (x 1 ) a (x 2 ) b with arbitrary "kinetic orders" a, b > 0 (in particular, different from 1).Using the complexes y = (1, 1, 0, 0, . ..)T , y ′ = (0, 0, 1, 0, . ..)T , and the kinetic-order complex ỹ = (a, b, 0, 0, . ..)T , we can write the reaction as y → y ′ with rate v = k x ỹ .For a network, the resulting dynamical system, is determined by the matrices Y (by stoichiometry), Ỹ (by kinetics), and A k (by a graph).For generalized mass-action systems, asymptotic stability of complexbalanced equilibria and non-existence of other steady states are not guaranteed (as for classical mass-action systems, cf.Theorem 5).We have already provided necessary conditions for linear stability of complex-balanced equilibria [6].In parallel work [29], we use the new decomposition of the graph Laplacian and monomial evaluation orders to study sufficient conditions for linear stability of complex-balanced equilibria and non-existence of other steady states.
On the other hand, every subgraph S ∈ G i gives rise to a spanning tree T ∈ T i ′ and vice versa (by removing/adding the edge i ′ → i that is in the cycle).For i ∈ V , Hence, where the sum is over all cycles C contained in G, and A C is the Laplacian matrix of the cycle C with k = 1 ∈ R E > (all edge labels set to 1).
Proof.Both matrices, A k diag(K k ) and C λ k,C A C , have zero row and column sums.Hence, it is sufficient to compare the off-diagonal entries.
On the one hand, every spanning tree in T i gives rise to a subgraph in G C that contains the edge i → i ′ in the cycle and vice versa (by adding/removing the edge i On the other hand, and the two matrices, A k diag(K k ) and C λ k,C A C , agree.
Remark.In a time-discrete, linear process otherwise, the edge labels k ∈ R n > do not represent transition rates, but transition probabilities.Then, (i→i ′ )∈E k i→i ′ ≤ 1, and B k is simply the matrix of transition probabilities with "k i→i "= 1 − (i→i ′ )∈E k i→i ′ column sums equal to one.That is, B k = A k + I, the identity matrix.Obviously, ψ = B k ψ if and only if A k ψ = 0. Whereas A k diag(K k ) always has zero row and column sums, B k may (or may not) be doubly stochastic (have column and row sums equal to one).
The Birkhoff/von Neumann Theorem [5,35] states that every doubly stochastic (d.s.) matrix B ∈ R n×n ≥ is the convex sum of permutation matrices; however, this decomposition is not unique.In fact, there are n! permutation matrices.Still, the polytope of d.s.matrices lies in an (n− 1) 2 -dimensional affine subspace of R n×n ≥ , and hence every d.s.matrix can be written as the sum of at most (n − 1) 2 + 1 permutation matrices.
On the contrary, the matrix A k diag(K k ) is the unique sum of all Laplacian matrices of cycles.However, there are more than (n − 1) 2 + 1 cycles, in general.

B Auxiliary graph-theoretic results
Lemma 12 (cf.[13], Lemma 2).Let G k = (V, E, k) be a connected, labeled, simple digraph with one absorbing strong component, and A k and I E be the corresponding Laplacian and incidence matrices.Then, im(A k ) = im(I E ).Lemma 13 (cf.[28], Proposition 5).Let G = (V, E) be a connected, simple digraph, G E = (V, E) be an auxiliary digraph, and I E and I E be the corresponding incidence matrices.Then, im(I E ) = im(I E ).
Proof.From graph theory and the definition of an auxiliary graph, we know that dim im(I E ) = dim im(I E ) = |V | − 1.In the rest of the proof, we show that im(I E ) ⊆ im(I E ).We consider the edge (i → j) ∈ E and the corresponding column e j − e i of I E , where e i denotes the ith standard basis vector in R V .
Since G E is a directed tree, there is a path from i to j in the undirected version of G E , that is, i = i 1 − − i 2 − − . . .− − i l = j with either (i k → i k+1 ) ∈ E or (i k ← i k+1 ) ∈ E for k = 1, . . ., l − 1.Hence, where α k ∈ {−1, 1} and e i k+1 − e i k is the column of I E corresponding to either the edge (i k → i k+1 ) ∈ E or (i k ← i k+1 ) ∈ E.
C A proof of Theorem 5 We provide a proof of Theorem 5 in the main text, based on the entropy-like Lyapunov function.Previous proofs further use inequalities for the exponential function or the logarithm and cycle decomposition of the graph, cf.[20,33,2,15].We use monomial evaluation orders and corresponding geometric objects (strata and polyhedral cones).
Theorem.Let (G k , y) be a mass-action system and x * ∈ R n > be a positive CBE of the dynamical system (4).Then, ln x x * T f k (x) < 0 for all x ∈ R n > that are not complex-balanced equilibria.Hence, (i) all positive equilibria are complex-balanced, and (ii) x * is asymptotically stable.
Proof.Let x ∈ R n > not be a CBE.Then there is a full-dimensional subset (a stratum) S k,E ⊂ R n > for some chain graph G E = (V, E) such that x ∈ S k,E , that is, ln x x * ∈ C E .Using the dynamical system (4) and Theorem 4, we have Using S k,E and C E as in Eqns.( 8) and ( 10), we have b ≥ 0 and a ≥ 0.
Since x is not be a CBE, b = 0, that is, there is i → i ′ ∈ E such that By complex balancing (6), (i) If there is a positive equilibrium x ∈ R n > that is not complex-balanced, then f k (x) = 0, contradicting ln x x * T f k (x) < 0.
Lemma 12 in Appendix B, and further im I E = im I E , cf.Lemma 13 in Appendix B. Altogether, we have im B k,E = im I E and hence B k,E = −I E A k,E for a unique matrix A k,E ∈ R E×E .(The minus sign ensures positive diagonal entries of A k,E for particular auxiliary graphs; see below.

Proof.
From graph theory, we know that dim im(I E ) = |V | − 1 and ker(A k ) = im ξ, where ξ ∈ R V ≥ has support on the absorbing strong component of G. Hence, also dim im(A k ) = |V | − 1.By definition, im(A k ) ⊆ im(I E ) and hence im(A k ) = im(I E ).
a i→i ′ = (y(i ′ ) − y(i)) T ln x x * > 0.By Theorem 4, the core matrix of the graph Laplacian, A k,E ∈ R E×E is nonnegative with positive diagonal.Hence,ln x x * T f k (x) = −a T A k,E b < 0.

(
ii) Recall that a positive CBE x * is the unique steady state in its stoichiometric compatibility class (forward invariant set).Hence,d dt L(x(t)) = ln x x * T f k (x) ≤ 0 with "=" if and only if x = x * ,and L(x) is a strict Lyapunov function.
and y ∈ int C pol if and only if y • x < 0 for all x ∈ C \ lineal C.