Complex and detailed balancing of chemical reaction networks revisited

The characterization of the notions of complex and detailed balancing for mass action kinetics chemical reaction networks is revisited from the perspective of algebraic graph theory, in particular Kirchhoff's Matrix Tree theorem for directed weighted graphs. This yields an elucidation of previously obtained results, in particular with respect to the Wegscheider conditions, and a new necessary and sufficient condition for complex balancing, which can be verified constructively.


Introduction
The notion of complex balancing of mass action kinetics chemical reaction networks, generalizing the classical notion of detailed balancing, dates back at least to the origin of chemical reaction network (CRN) theory; see especially [6,11,12]. The assumption of existence of a complex-balanced equilibrium has powerful consequences for the dynamical behavior, precluding multi-stability and oscillations. In this note we will revisit the notion of complex balancing by a systematic use of notions and results from algebraic graph theory, in particular the Laplacian matrix and Kirchhoff's Matrix Tree theorem (see also [10,14] for other uses of this theorem in chemical reaction dynamics). This will result in a constructive necessary and sufficient condition for complex balancing. Furthermore, motivated by recent work in [4] expanding on [7], we will provide a new perspective and results on the Wegscheider conditions for detailed balancing and formal balancing as introduced in [4].
The structure of this note is as follows. Section 2 gives a brief recap of the basic framework of CRN theory from an algebraic graph theory perspective, based on [18,21,23]. Section 3 introduces Kirchhoff's Matrix Tree theorem and shows how the application of this theorem leads to an improved, and more directly verifiable, condition for complex balancing as compared to [6,11]. Section 4 relates Kirchhoff's Matrix Tree theorem to the notion of formal balancing, cf. [4] and [7], and shows how this leads to an insightful graph-theoretic proof of the result obtained in [4] that complex balancing together with formal balancing implies detailed balancing and conversely. Notation: The space of n-dimensional real vectors consisting of all strictly positive entries is denoted by R n + . The mapping Ln : R n + → R n , x → Ln x, is the elementwise logarithm, and is defined as the mapping whose i-th component is given by ln(x i ). Similarly, Exp : R n → R n + is the mapping whose i-th component is given by exp x i . Furthermore, for two vectors x, y ∈ R n + we let x y denote the vector in R n + with i-th component x i y i . Finally, 1 n denotes the n-dimensional vector with all entries equal to 1. Some graph-theoretic notions [3]: A directed graph 1 G with c vertices and r edges is characterized by its c × r incidence matrix, denoted by D. Each column of D corresponds to an edge of the graph, and contains exactly one element 1 at the position of the head vertex of this edge and exactly one −1 at the position of its tail vertex; all other elements are zero. Clearly, 1 T c D = 0. The graph is connected if any vertex can be reached from any other vertex by following a sequence of edges; direction not taking into account. It holds that rank D = c − , where is the number of connected components of the graph. In particular, G is connected if and only if ker D T = span 1 c . The graph is strongly connected if any vertex can be reached from any other vertex, following a sequence of directed edges. A subgraph of G is a directed graph whose vertex and edge set are subsets of the vertex and edge set of G. A graph is acyclic (or, does not contain cycles) if and only if ker D = 0. A spanning tree of a directed graph G is a connected, acyclic subgraph of G that contains all vertices of G.

Recall of complex-balanced chemical reaction networks
In this section, in order to set the stage, we will briefly recall the well-established framework of (isothermal) CRN theory, originating in the work of Horn, Jackson and Feinberg in the 1970s [6,9,11,12]. Consider a chemical reaction network with m chemical species (metabolites) with concentrations x ∈ R m + , among which r chemical reactions take place. The left-hand sides of the chemical reactions are called substrate complexes and the right-hand sides the product complexes. To each chemical complex (substrate and/or product) of the reaction network one can associate a vertex of a graph, and to every reaction (from substrate to product complex) a directed edge (with tail vertex the substrate and head vertex the product complex). Let c be the total number of complexes involved in the chemical reaction network, then the resulting directed graph G with c vertices and r edges is called the graph of complexes 2 , and is defined by its c×r incidence matrix D. Furthermore we define the m × c complex composition matrix 3 Z with non-negative integer elements expressing the composition of the complexes in terms of the chemical species: its k-th column expresses the composition of the k-th complex. The dynamics of the chemical reaction network takes the well-known forṁ where v(x) is the vector of reaction rates, and S = Z D the stoichiometric matrix.
The most basic way to define v(x) is mass action kinetics. For example, the mass action kinetics reaction rate of the reaction X 1 + 2X 2 → X 3 is given as v(x) = kx 1 x 2 2 with k > 0 a reaction constant. In general, for a single reaction with substrate complex S specified by its corresponding column Z S = Z S1 . . . Z Sm T of the complex composition matrix Z , the mass action kinetics reaction rate is given by which can be rewritten as k exp(Z T S Ln x), x ∈ R m + . Hence the reaction rates of the total reaction network are given by where S j is the substrate complex of the j-th reaction with reaction constant k j > 0. This yields the following compact description of the rate vector v(x). Define the r × c matrix K as the matrix whose ( j, σ )-th element equals k j if the σ -th complex is the substrate complex for the j-th reaction, and zero otherwise. Then v(x) = K Exp (Z T Ln x), x ∈ R m + , and the dynamics of the mass action reaction network takes the formẋ It can be easily verified that the c × c matrix L := −DK has nonnegative diagonal elements and nonpositive off-diagonal elements. Moreover, since 1 T c D = 0 also 1 T c L = 0, i.e., the column sums of L are all zero. Hence L defines a weighted Laplacian matrix 4 for the graph of complexes G.
The aim of CRN theory, starting with [6,11,12], is to analyze the dynamical properties of (2), and in particular to derive conditions which ensure a dynamical behavior which is independent of the precise values of the reaction constants (which are often poorly known or varying). This has culminated in the deficiency-zero and deficiencyone theorems (see e.g. [8]), while a somewhat complementary approach is based on the assumption of existence of a complex-balanced equilibrium [6,11], generalizing the classical notion of a detailed-balanced equilibrium.
Chemically (3) means that at the complex-balanced equilibrium x * not only all the chemical species, but also the complexes remain constant; i.e., for each complex the total inflow (from the other complexes) equals the total outflow (to the other complexes).
The assumption of complex balancing has been shown to have strong implications for the dynamical properties of (2); see in particular the classical papers [6,11,12]. As detailed in [18], expanding on [21], these properties can be easily proved by defining the diagonal matrix and rewriting the dynamics (2) into the forṁ The key point is that since and thus is a balanced Laplacian matrix (column sums and row sums are zero). Together with convexity of the exponential function this implies the following key fact.
First property which directly follows [18] from Proposition 2.2 is the classical result [6,11,12] that all positive equilibria are in fact complex-balanced, and that given one complex-balanced equilibrium x * the set of all positive equilibria is given by Furthermore, using an elegant result from [8], there exists for every initial condition x 0 ∈ R m + a unique x * * ∈ E such that x * * − x 0 ∈ im S. By using the Lyapunov function Proposition 2.2 then implies that the vector of concentrations x(t) starting from x 0 will converge to x * * ; at least if the reaction network is persistent 6 . The chemical interpretation is that G is (up to a constant) the Gibbs' free energy with gradient vector ∂G ∂ x (x) = Ln x x * * being the chemical potentials. See e.g. [15,21,22] for further information 7 .

A graph-theoretic characterization of complex-balancing
In this section we will expand on earlier investigations to characterize the existence of a complex-balanced equilibrium (see in particular [4,6,11], and the references quoted therein), and derive a new necessary and sufficient condition for complex balancing which can be constructively verified.
First note that by the definition of Ln : or equivalently, Exp (Z T μ * ) ∈ ker L. Furthermore, note that Exp (Z T μ * ) ∈ R c + . In case the graph G is connected the kernel of L is 1-dimensional, and a vector ρ ∈R c + (the closure of the positive orthant) with ρ ∈ ker L can be computed by what is sometimes called Kirchhoff's Matrix Tree theorem 8 , which for our purposes can be summarized as follows. Denote the (i, j)-th cofactor of L by C i j = (−1) i+ j M i, j , where M i, j is the determinant of the (i, j)-th minor of L, which is the matrix obtained from L by deleting its i-th row and j-th column. Define the adjoint matrix adj(L) as the matrix with (i, j)-th element given by C ji . It is well-known that 6 The reaction network is called persistent if for every x 0 ∈ R m + the ω-limit set of the dynamics (2) does not intersect the boundary ofR m + . It is generally believed that most reaction networks are persistent. However, up to now this persistence conjecture has been only proved in special cases (cf. [1,2,19] and the references quoted in there). 7 The form (5) also provides a useful starting point for structure-preserving model reduction [17,18,21].
L · adj(L) = (det L)I c = 0 ( 1 0 ) Furthermore, since 1 T c L = 0 the sum of the rows of L is zero, and hence by the properties of the determinant function it directly follows that C i j does not depend on i; implying that C i j = ρ j , j = 1, . . . , c. Therefore by defining ρ := (ρ 1 , . . . , ρ c ) T , it follows from (10) that Lρ = 0. Furthermore, cf. [3, Theorem14 on p.58], ρ i is equal to the sum of the products of weights of all the spanning trees of G directed towards vertex i. In particular, it follows that ρ j ≥ 0, j = 1, . . . , c. In fact, ρ = 0 if and only if G has a spanning tree. Furthermore, since for every vertex i there exists at least one spanning tree directed towards i if and only if the graph is strongly connected, we may conclude that ρ ∈ R c + if and only if the graph is strongly connected. Example 3.1 Consider the cyclic reaction network The Laplacian matrix is given as By Kirchhoff's Matrix Tree theorem the corresponding vector ρ satisfying Lρ = 0 is given as where each term corresponds to one of the three weighted spanning trees pointed towards the three vertices.
In case the graph G is not connected the same analysis can be performed on any of its connected components.
Returning to the existence of μ * ∈ R m satisfying LExp (Z T μ * ) = 0 this implies the following. Let G j , j = 1, . . . , , be the connected components of the graph of complexes G. For each connected component, define the vectors ρ 1 , . . . , ρ as above by Kirchhoff's Matrix Tree theorem (i.e., as cofactors of the corresponding diagonal sub-blocks of L or as sums of products of weights along spanning trees). Define the total vector ρ as the stacked column vector ρ := col(ρ 1 , . . . , ρ ). Partition correspondingly the composition matrix Z as Z = [Z 1 . . . Z ]. Then there exists μ * ∈ R m satisfying LExp (Z T μ * ) = 0 if and only if each connected component G j , j = 1, . . . , , is strongly connected and This in its turn is equivalent to strong connectedness of each connected component G j and the existence of constants β j such that Z T j μ * = Ln ρ j + β j 1, j = 1, . . . , . Furthermore, this is equivalent to strong connectedness of each connected component of G, and Ln ρ ∈ im Z T + ker D T Finally, (12) is equivalent to Summarizing we have obtained Remark 3.4 The above theorem is a restatement of Theorem 3C in [11]; the main difference being that in [11] the positive vector ρ ∈ ker L remains unspecified, while in our case it is explicitly given by Kirchhoff's Matrix Tree theorem.
Example 3.6 Consider Example 3.1 for the special case k − 1 = k − 2 = k − 3 = 0 (irreversible reactions). Then the vector ρ reduces to The reaction network with complex composition matrix Z is complex-balanced if and only k + i > 0, i = 1, 2, 3, and This last condition can be further written out as ⎡ As a mathematical example, take the complex composition matrix Z = 1 0 2 1 1 1 (corresponding to the complexes X 1 + X 2 , X 2 , 2X 1 + X 2 ). In this case the network is complex-balanced if and only if k + then there also is a directed edge from j to i (and in the case of multiple edges from i to j there are as many edges from i to j as edges from j to i). This means that the connected components of G are always strongly connected, or equivalently, that ρ as obtained from Kirchhoff's Matrix Tree theorem is in R c + . Define the undirected graphḠ as having the same vertices as G but half its number of edges, by replacing every pair of oppositely directed edges of G by one undirected edge ofḠ. Denote the number of edges ofḠ byr = 1 2 r . Endow subsequentlyḠ with an arbitrary orientation (all results in the sequel will be independent of this orientation), and denote the resulting incidence matrix byD. Clearly, after possible reordering of the edges,D is related to the incidence matrix D of G as To the j-th edge ofḠ there now correspond two reaction constants k + j , k − j (the forward and reverse reaction constants with respect to the chosen orientation ofḠ). Then define the equilibrium constants K eq j := k + j k − j , j = 1, . . . ,r , and the vector K eq := (K eq 1 , . . . , K eq r ) T . Recall [7,24] that the reaction network is called detailed-balanced 11

if and only if it satisfies
Ln whereS := ZD is the stoichiometric matrix of the reversible network with graphḠ. This is equivalent to σ 1 lnK eq 1 + · · · + σr lnK eq r = 0 for all σ = (σ 1 , . . . , σr ) such that σ TST = 0. Writing out K eq j = k + j k − j this is seen to be equivalent to for all σ such thatSσ = 0, known as the (generalized) Wegscheider conditions. Recently in [4] the notion of formally balanced was introduced, based on Feinberg's circuit conditions in [7], and weakening the above Wegscheider conditions. In our setup this notion is defined as follows.

Definition 4.1
The reversible reaction networkḠ with incidence matrixD, and vector of equilibrium constants K eq is called formally balanced if Ln K eq ∈ imD T SinceS = ZD 'formally balanced' is trivially implied by 'detailed-balanced', while if imS T = imD T (zero-deficiency) the reverse holds. Furthermore, any reaction network with acyclicḠ (and thus kerD = 0) is automatically formally balanced.
As above, the notion of 'formally balanced' is seen to be equivalent to σ 1 lnK eq 1 + · · · + σ r lnK eq r = 0 for all σ = (σ 1 , . . . , σr ) T such that σ TDT = 0, which in turn is equivalent to (19) for all σ such thatDσ = 0 (that is, for all cycles σ ). We will refer to (19) as the weak Wegscheider conditions. Note that the weak Wegscheider conditions only depend on the structure of the graphḠ (i.e., its cycles) and the equilibrium constants, and not on the complex composition matrix Z as in the case of the 'strong' Wegscheider conditions (18).

Theorem 4.2
Consider a reversible chemical reaction network given by the graphḠ with incidence matrixD, and with ρ determined by L. The following statements are equivalent Let us first prove the equivalence between (1) and (2). Consider the i-th and j-th vertex ofḠ, and suppose that the orientation has been taken such that the α-th edge between i and j is such that i is the tail vertex and j is the head vertex. Then the (i, j)-th element of Ldiag (ρ 1 , . . . , ρ c ) is given by k − α ρ j , while the ( j, i)-th element equals k + α ρ i . Symmetry of Ldiag (ρ 1 , . . . , ρ c ) thus amounts to for all pairs of vertices i, j. On the other hand, the α-th element of the vectorD T Ln ρ is given by ln ρ j − ln ρ i while the α-th element of Ln (K eq ) is given by Equality of ln ρ j − ln ρ i and ln k + α − ln k − α is thus equivalent to ln ρ j + ln k − α = ln ρ i + ln k + α which in its turn is equivalent to k + α ρ i = k − α ρ j as above.
Obviously (2) implies (3). For the reverse implication, consider any pair of vertices linked by an edge of the graphḠ. Depending on the orientation ofḠ refer to one vertex as the tail vertex t and the other vertex as the head vertex h. Refer to the positive reaction constant from t to h by k + and to the negative reaction constant by k − . Now consider a spanning tree directed towards t with product of weights denoted by τ t . In case the edge between t and h is part of this spanning tree then it follows that by reversing the orientation of this edge it defines a spanning tree directed towards to h with product of weights denoted by τ h . It follows directly that In case the edge between t and h is not part of this spanning tree then divide the edges of the spanning tree into two sets; the set E 1 containing the edges of the part of the spanning tree from t to h (containing say edges) and the set E 2 containing the remaining edges of the spanning tree. Observe that E 1 together with the edge from t to h forms a cycle of the graphḠ. Since Ln (K eq ) ∈ imD T it follows that for the reaction constants along this cycle (choosing an appropriate orientation), Now within the spanning tree directed towards t, if the orientation of each of the edges of E 1 is reversed, we obtain another spanning tree directed towards h with product of weights denoted again by τ h . By using (21) it is readily verified that also in this case we obtain the same relation (20). Summing up over all spanning trees we thus obtain the equality which can be equivalently written as ln k + k − = ln ρ h − ln ρ t . Doing this for all adjacent vertices t and h this exactly amounts to the required equality Ln (K eq ) =D T Ln ρ. Example 4.3 Consider again the reaction network described in Example 3.1 (without specifying the complexes C 1 , C 2 , C 3 ). The transformed Laplacian matrix is computed as On the other hand, Ln K eq ∈D T amounts to = 0, and hence to the same condition (23).
Thus the reaction network is formally balanced if and only if (23) holds.
Now let us relate all this to the necessary and sufficient conditions for complex balancing obtained before, cf. (13). Note that by (16) D T Ln ρ ∈ im D T Z T ⇔D T Ln ρ ∈ imD T Z T Hence a reversible reaction network is complex-balanced if and only ifD T Ln ρ ∈ imD T Z T = imS T .
We directly obtain the following corollary proved by other methods in [4] 12 .

Conclusions
By a systematic use of notions from algebraic graph theory, in particular the Laplacian matrix and Kirchhoff's Matrix Tree theorem, previously derived results on complex, detailed and formal balancing have been proved in a simple manner. Furthermore, it has resulted in a new necessary and sufficient condition for complex balancing, which can be verified constructively.
The results obtained in this note can be immediately extended to mass action kinetic reaction networks with constant inflows and mass action outflow exploiting the classical idea of adding a 'zero complex'; see [9] and [23] for further details.
Current research is concerned with the application of the developed framework to questions of occurrence of multi-stability and structure-preserving model reduction; see for the latter also [17,18,21].