Toric Dynamical Systems arXiv:0708.3431v2 [math.DS] 3 Nov 2007 Gheorghe Craciun Dept. of Mathematics, University of Wisconsin, Madison, WI 53706-1388, USA Alicia Dickenstein Dep. de Matemática, FCEN, Universidad de Buenos Aires, (1428), Argentina Anne Shiu Dept. of Mathematics, University of California, Berkeley, CA 94720-3840, USA Bernd Sturmfels Dept. of Mathematics, University of California, Berkeley, CA 94720-3840, USA Abstract Toric dynamical systems are known as complex balancing mass action systems in the mathematical chemistry literature, where many of their remarkable properties have been established. They include as special cases all deficiency zero systems and all detailed balancing systems. One feature is that the steady state locus of a toric dynamical system is a toric variety, which has a unique point within each invariant polyhedron. We develop the basic theory of toric dynamical systems in the context of computational algebraic geometry and show that the associated moduli space is also a toric variety. It is conjectured that the complex balancing state is a global attractor. We prove this for detailed balancing systems whose invariant polyhedron is two-dimensional and bounded. This paper is dedicated to the memory of Karin Gatermann (1961–2005). Key words: chemical reaction network, toric ideal, complex balancing, detailed balancing, deficiency zero, trajectory, Birch’s Theorem, Matrix-Tree Theorem, moduli space, polyhedron Preprint submitted to Elsevier 1 February 2008 1. Introduction Toric dynamical systems describe mass-action kinetics with complex balancing states. These systems have been studied extensively in mathematical chemistry, starting with the work of Horn and Jackson (1972), Feinberg (1972) and Horn (1972, 1973), and continuing with the deficiency theory in (Feinberg, 1979, 1987, 1989, 1995). Mass-action kinetics has a wide range of applications in the physical sciences, and now it is beginning to play a role in systems biology (Craciun, Tang and Feinberg, 2006; Gnacadja et al., 2007; Gunawardena, 2003; Sontag, 2001). Important special cases of these dynamical systems include recombination equations in population genetics (Akin, 1979) and quadratic dynamical systems in computer science (Rabinovich, Sinclair and Wigderson, 1992). Karin Gatermann introduced the connection between mass-action kinetics and computational algebra. Our work drew inspiration both from her publications (Gatermann, 2001; Gatermann and Huber, 2002; Gatermann and Wolfrum, 2005) and from her unpublished research notes on toric dynamical systems. We wholeheartedly agree with her view that “the advantages of toric varieties are well-known” (Gatermann, 2001, page 5). We now review the basic set-up. A chemical reaction network is a finite directed graph whose vertices are labeled by monomials and whose edges are labeled by parameters. The digraph is denoted G = (V, E), with vertex set V = {1, 2, . . . , n} and edge set E ⊆ {(i, j) ∈ V × V : i 6= j}. The node i of G represents the ith chemical complex and is labeled with the monomial cyi = cy1i1 cy2i2 · · · cys is . Here Y = (yij ) is an n × s-matrix of non-negative integers. The unknowns c1 , c2 , . . . , cs represent the concentrations of the s species in the network, and we regard them as functions ci (t) of time t. The monomial labels are the entries in the row vector
Ψ(c) = cy1 , cy2 , . . . , cyn . Each directed edge (i, j) ∈ E is labeled by a positive parameter κij which represents the rate constant in the reaction from the i-th chemical complex to the j-th chemical complex. Note that if there is an edge from i to j and an edge from j to i then we have two unknowns κij and κji . Let Aκ denote the negative of the Laplacian of the digraph G. Hence Aκ is the n × n-matrix whose off-diagonal entries are the κij and whose row sums are zero. Mass-action kinetics specified by the digraph G is the dynamical system dc = Ψ(c) · Aκ · Y. (1) dt A toric dynamical system is a dynamical system (1) for which the algebraic equations Ψ(c) · Aκ = 0 admit a strictly positive solution c∗ ∈ Rs>0 . Such a solution c∗ is a steady state of the system, i.e., the s coordinates of Ψ(c∗ ) · Aκ · Y vanish. The requirement that all n coordinates of Ψ(c∗ ) · Aκ be zero is stronger. The first to study toric dynamical ⋆ GC was supported by was supported by the NSF (DMS-0553687) and DOE (BACTER Institute, DEFG02-04ER25627). AD was supported by UBACYT X042, CONICET PIP 5617 and ANPCyT PICT 20569, Argentina. AS was supported by a Lucent Technologies Bell Labs Graduate Research Fellowship. BS was supported by the NSF (DMS-0456960) and DARPA (Fundamental Laws in Biology, HR001105-1-0057). Email addresses: craciun@math.wisc.edu (Gheorghe Craciun), alidick@dm.uba.ar (Alicia Dickenstein), annejls@math.berkeley.edu (Anne Shiu), bernd@math.berkeley.edu (Bernd Sturmfels). 2 systems, Horn and Jackson (1972), called these systems complex balancing mass action systems and called c∗ a complex balancing steady state. A system (1) being complex balancing (i.e., toric) depends on both the digraph G and the rate constants κij . Example 1. Let s = 2, n = 3 and let G be the complete bidirected graph on three nodes labeled by c21 , c1 c2 and c22 . Here the mass-action kinetics system (1) equals     2 0 −κ12 − κ13 κ12 κ13      
d     = c 1 , c2  ·  1 1  (2) c21 c1 c2 c22 ·  κ21 −κ21 − κ23 κ23 dt     0 2 κ31 κ32 −κ31 − κ32 This is a toric dynamical system if and only if the following algebraic identity holds: (κ21 κ31 + κ32 κ21 + κ23 κ31 )(κ13 κ23 + κ21 κ13 + κ12 κ23 ) = (κ12 κ32 + κ13 κ32 + κ31 κ12 )2 . (3) The equation (3) appears in (Horn, 1973, Equation (3.12)) where it is derived from the necessary and sufficient conditions for complex balancing in mass-action kinetics given by Horn (1972). Our results in Section 2 provide a refinement of these conditions. Let us now replace G by the digraph with four edges (1, 3), (2, 1), (2, 3), (3, 1). This corresponds to setting κ12 = κ32 = 0 in (3). We can check that, for this new G, the system (1) is not toric for any positive rate constants. Note that G is not strongly connected. ✷ Among all chemical reaction networks, toric dynamical systems have very remarkable properties. Some of these properties are explained in (Feinberg, 1979), starting with Proposition 5.3; see also (Gunawardena, 2003, Theorem 6.4). We shall review them in detail in Sections 2 and 3. From our point of view, the foremost among these remarkable properties is that the set Z of all steady states is a toric variety (Gatermann, 2001, §3). Each trajectory of (1) is confined to a certain invariant polyhedron, known to chemists as the stoichiometric compatibility class, which intersects the toric variety Z in precisely one point c∗ . In order to highlight the parallels between toric dynamical systems and toric models in algebraic statistics (Pachter and Sturmfels, 2005, §1.2), we shall refer to the steady state c∗ as the Birch point; see (Sturmfels, 1996, Theorem 8.20). In Example 1, the steady state variety Z is a line through the origin, and the Birch point equals
c∗ = const · κ12 κ32 + κ13 κ32 + κ31 κ12 , κ13 κ23 + κ21 κ13 + κ12 κ23 . Here the constant is determined because c1 + c2 is conserved along trajectories of (2). This article is organized as follows. In Section 2 we develop the basic theory of toric dynamical systems within the context of computational algebraic geometry. For each directed graph G we introduce the moduli space of toric dynamical systems on G. This space parametrizes all rate constants κ for which (1) is toric. In Example 1 this space is the hypersurface (3). Theorem 9 states that this moduli space is itself a toric variety in a suitable system of coordinates. These coordinates are the maximal non-zero minors of the Laplacian of G, and their explicit form as positive polynomials in the κij is given by the Matrix-Tree Theorem (Stanley, 1999, §5.6). Our results in Section 2 furnish a two-fold justification for attaching the adjective “toric” to chemical reaction networks with complex balancing, namely, both the steady state variety and the moduli space are toric. In addition, the subvariety of reaction networks with detailed balancing is toric. 3 In Section 3 we introduce the Global Attractor Conjecture which states that the Birch point is a global attractor for any toric dynamical system. More precisely, we conjecture that all trajectories beginning at strictly positive vectors c0 will converge to the Birch point c∗ in the invariant polyhedron of c0 . The conjecture is currently open, even for deficiency zero systems (cf. Theorem 9). De Leenheer, Angeli and Sontag (2007) found a proof for a class of “monotone” deficiency zero networks where the monomials cyi involve distinct unknowns. We prove the conjecture in Section 5 for toric dynamical systems with detailed balancing that evolve in a bounded polygon in s-dimensional space. The algebraic theory of detailed balancing systems is developed in Section 4. 2. Ideals, Varieties and Chemistry This section concerns the connection between chemical reaction network theory and toric geometry. We use the language of ideals and varieties as in (Cox, Little and O’Shea, 2007). Our reference on toric geometry and its relations with computational algebra is (Sturmfels, 1996). With regard to the dynamical system (1), we use the notation from (Feinberg, 1979, §5) and (Gunawardena, 2003, §3) which has the virtue of separating the roles played by the concentrations ci , the monomials cyi , and the rate constants κij . To study the dynamical system (1) algebraically, we work in the polynomial ring   Q[c, κ] = Q {c1 , c2 , . . . , cs } ∪ {κij : (i, j) ∈ E} , and we introduce various ideals in this polynomial ring. First, there is the steady state ideal hΨ(c) · Aκ · Y i which is generated by the s entries of the row vector on the right hand side of (1). Second, we consider the ideal hΨ(c) · Aκ i which is generated by the n entries of the row vector Ψ(c) · Aκ . The generators of both ideals are linear in the κij but they are usually non-linear in the ci . Next, we define the complex balancing ideal of G to be the following ideal quotient whose generators are usually non-linear in the κij :
CG := hΨ(c) · Aκ i : (c1 c2 · · · cs )∞ . We have thus introduced three ideals in Q[c, κ]. They are related by the inclusions hΨ(c) · Aκ · Y i ⊆ hΨ(c) · Aκ i ⊆ CG . If I is any polynomial ideal then we write V (I) for its complex variety. Likewise, we define the positive variety V>0 (I) and the non-negative variety V≥0 (I). They consist of all points in V (I) whose coordinates are real and positive or, respectively, non-negative. Our algebraic approach to chemical reaction network theory focuses on the study of these varieties. The inclusions of ideals above imply the following inclusions of varieties:

V (CG ) ⊆ V hΨ(c) · Aκ i ⊆ V hΨ(c) · Aκ · Y i . (4) The definition of CG by means of saturation implies that the left hand inclusion becomes equality when we restrict to the points with all coordinates non-zero. In particular,
(5) V>0 (CG ) = V>0 hΨ(c) · Aκ i . Recall from (Sturmfels, 1996) that a toric ideal is a prime ideal which is generated by binomials. We soon will replace CG by a subideal TG which is toric. This is possible by Proposition 5.3 (ii,iv) in (Feinberg, 1979) or Theorem 6.4 (3) in (Gunawardena, 2003), which essentially state that V>0 (CG ) is a positive toric variety. But let us first examine the case when CG is a toric ideal already. 4 Example 2. Suppose that each chemical complex appears in only one reaction, and each reaction is bi-directional. Hence n = 2m is even and, after relabeling, we have E = {(1, 2), (2, 1), (3, 4), (4, 3), . . . , (n−1, n), (n, n−1)}. We start with the binomial ideal hΨ(c) · Aκ i = κ12 cy1 − κ21 cy2 , κ34 cy3 − κ43 cy4 , . . . , κn−1,n cyn−1 − κn,n−1 cyn . The complex balancing ideal CG is a saturation of hΨ(c) · Aκ i, and it coincides with the toric ideal of the extended Cayley matrix in the proof of Theorem 7. There are many programs for computing toric ideals. For instance, the methods in (Sturmfels, 1996, §12.A) are available in maple under the command ToricIdealBasis. Explicitly, the complex balancing ideal CG is generated by all binomials κu+ cv+ − κu− cv− where m X u2i−1,2i (y2i−1 − y2i ) = v and u2i−1,2i + u2i,2i−1 = 0 for i = 1, 2, . . . , m. (6) i=1 Eliminating c1 , . . . , cs from CG , we obtain the ideal of all binomials κu+ − κu− where u ∈ NE satisfies (6) with v = 0. This is the moduli ideal MG to be featured in Theorems 7 and 9 below. It is a prime binomial ideal of Lawrence type (Sturmfels, 1996, §7). ✷ Let us next assume that G = (V, E) is an arbitrary digraph with n nodes which is strongly connected. This means that, for any two nodes i and j, there exists a directed path from i to j. In this case the matrix Aκ has rank n − 1, and all its minors of size (n−1)×(n−1) are non-zero. The next result gives a formula for these comaximal minors. Consider any directed subgraph T of G whose underlying graph is a tree. This means that T has n − 1 edges and contains no cycle. We write κT for the product of all edge labels of the edges in T . This is a squarefree monomial in Q[κ]. Let i be one of the nodes of G. The directed tree T is called an i-tree if the node i is its unique sink, i.e., all edges are directed towards node i. We introduce the following polynomial of degree n − 1: X Ki = κT . (7) T an i-tree The following result is a restatement of the Matrix-Tree Theorem (Stanley, 1999, §5.6). Proposition 3. Consider a submatrix of Aκ obtained by deleting the ith row and any one of the columns. The signed determinant of this (n−1) × (n−1)-matrix equals (−1)n−1 Ki . This minor is independent of the choice of columns because the row sums of Aκ are zero. Combining Proposition 3 with a little linear algebra leads to the following corollary: Corollary 4. The complex balancing ideal CG contains the polynomials Ki cyj − Kj cyi .
We now form the ideal generated by these n2 polynomials and we again saturate with respect to c1 c2 · · · cs . The resulting ideal TG will be called the toric balancing ideal:
TG := hKi cyj − Kj cyi : 1 ≤ i < j ≤ ni : (c1 c2 · · · cs )∞ . It is thus natural to consider TG as an ideal in the polynomial subring Q[c, K] = Q[c1 , . . . , cs , K1 , . . . , Kn ] ⊂ Q[c, κ]. The claim that this is a polynomial ring is the content of the following lemma. 5 Lemma 5. The polynomials K1 , . . . , Kn ∈ Q[κ] are algebraically independent over Q. Proof. Let Ki′ ∈ Q[κ1 , κ2 , . . . , κn ] denote the polynomial obtained from Ki by substituting the new unknown κi for all κij . We need only verify that the Ki′ are algebraically independent, because an algebraic relation among the Ki would be satisfied by the Ki′ as well. Our polynomials are Y Ki′ = (number of i-trees in G) · κt . t6=i Q The n squarefree monomials t6=i κt (for i = 1 . . . n) are algebraically independent because an algebraic dependence among these monomials would specify a dependence among 1/κ1 , 1/κ2 , . . . , 1/κn . Hence, K1′ , K2′ , . . . , Kn′ are algebraically independent. ✷ We new discuss the toric balancing ideal TG . Proposition 6. The toric balancing ideal TG is a toric ideal in Q[c, K]. Moreover, the ideal TG is generated by the binomials K u+ · c(uY )− − K u− · c(uY )+ where u is any row vector in Zn whose coordinate sum u1 + u2 + · · · + un is equal to zero. Proof. Let ∆ denote the
edge-node incidence matrix of the complete directed graph on n nodes. Thus ∆ is the n2 × n-matrix whose rows are ei − ej for 1 ≤ i < j ≤ n. We also
consider the n × (n + s)-matrix −Y In
. The binomials Ki cyj − Kj c
yi which define the ideal TG correspond to the rows of the n2 ×(n+s)-matrix ∆· −Y In , and the binomial
K u+ ·c(uY )− −K u− ·c(uY )+ corresponds to the row vector U ·∆· −Y In , where U is any
row vector of length n2 such that u = U · ∆. The binomial K u+ · c(uY )− − K u− · c(uY )+ yi ±1 yj is a Q[c±1 1 , . . . , cs , K1 , . . . , Kn ]-linear combination of the binomials Ki c − Kj c . This shows that TG is the
lattice ideal in Q[c, K] associated with the lattice spanned by the rows of ∆ · −Y In , i.e., there are no monomial zero-divisors modulo TG . To see that TG is actually a toric ideal, i.e. TG is prime,
it suffices to note that Zn+s modulo the lattice spanned by the rows of ∆ · −Y In is free abelian of rank s + 1. Indeed, the latter matrix has rank n − 1, and its (n − 1) × (n − 1)-minors span the unit ideal in the ring of integers Z, because each (n − 1) × (n − 1)-minor of ∆ is either +1 or −1. ✷ The variety of TG is a toric variety in Spec Q[c, K], but we continue to regard it as a subvariety of Cs × CE (or of Spec Q[c, κ]). In this interpretation we have
V>0 (TG ) = V>0 (CG ) = V>0 hΨ(c) · Aκ i . (8) Thus TG still correctly describes the steady state locus of the toric dynamical system. The equation (8) holds because the matrix Aκ has rank n − 1 over the rational function field Q(κ), and the vector (K1 , K2 , . . . , Kn ) spans its kernel under left multiplication. Finally, the following elimination ideal is called the moduli ideal of the digraph G: MG = TG ∩ Q[κ]. (9) Here Q[κ] is the polynomial ring in only the edge unknowns κij . The generators of MG are obtained from the generators of CG by eliminating the unknown concentrations ci . For instance, if G is the complete bidirected graph on c21 , c1 c2 and c22 as in Example 1 6 then the moduli ideal MG is the principal ideal generated by K1 K3 − K22 . This coincides with condition (3) because K1 = κ21 κ31 + κ32 κ21 + κ23 κ31 , and similarly for K2 , K3 . Suppose now that G is an arbitrary directed graph, and let l be the number of connected components of G. If one of the components Gi fails to be strongly connected, then V>0 (CGi ) is empty and hence V>0 (CG ) is empty, by (Feinberg, 1979, Remark 5.2). In that case we define TG and MG to be the ideal generated by 1. If each connected component Gi of G is strongly connected then we define the toric steady state ideal as
TG := ( TG1 + TG2 + · · · + TGl ) : (c1 c2 · · · cs )∞ . The moduli ideal MG is defined as before in (9). The equality in (8) still holds and this positive variety is in fact non-empty. Here is the first main result of this section: Theorem 7. The equations (1) specify a toric dynamical system if and only if the positive vector of rate constants κij lies in the toric variety V (MG ). In this case, the set of steady states of (1) with all ci > 0 equals the set of positive points on the toric variety V (TG ). Proof. The positive variety V>0 (TG ) consists of all pairs (c, κ) where κ is a strictly positive vector of rate constants and c is a strictly positive solution of the complex balancing equations Ψ(c) · Aκ = 0. The elimination in (9) corresponds to the map of toric varieties V (TG ) → V (MG ) given by (c, κ) 7→ κ. This map is a dominant morphism (by definition of MG ), so its image is Zariski dense in V (MG ). The restriction to real positive points, V>0 (TG ) → V>0 (MG ), is a homomorphism of abelian groups (R>0 )∗ whose image is dense, so it is the monomial map specified by a matrix with maximal row rank. It follows that this restriction is surjective, and this proves our first assertion. The second assertion follows from (Feinberg, 1979, Proposition 5.3). ✷ We now justify calling V (MG ) a toric variety by writing MG explicitly as a toric ideal in Q[K]. As before, G is a directed graph with n nodes labeled by monomials cy1 , . . . , cyn . We assume that each connected component G1 , G2 , . . . , Gl of G is strongly connected, for otherwise MG = h1i. Let Yi denote the matrix with s rows whose columns are the vectors yj where j runs over the nodes of the component Gi . We define the Cayley matrix   Y1 Y2 · · · Yl     0 ··· 0   1     1 ··· 0 . CayG (Y ) =  0   .  ..  .. ..  . . ..  .   0 0 ··· 1 This is an (s + l) × n-matrix. Here 1 and 0 are appropriate row vectors with all entries 1 and 0 respectively. The term “Cayley matrix” comes from geometric combinatorics, and it refers to the Cayley trick in elimination theory (Huber, Rambau and Santos, 2000). Let S denote the linear subspace of Rs which is spanned by the reaction vectors yj − yi where (i, j) ∈ E. This space is known in chemistry as the stoichiometric subspace. We write σ = dim(S) for its dimension. The quantity δ := n−σ−l is known as the deficiency of the chemical reaction network G. For instance, δ = 3 − 1 − 1 = 1 in Example 1. 7 Remark 8. The rank of the Cayley matrix CayG (Y ) equals σ + l. Hence the deficiency δ of the reaction network coincides with the dimension of the kernel of the Cayley matrix. The following theorem is the second main result in this section. Theorem 9. The moduli ideal MG equals the toric ideal of the Cayley matrix CayG (Y ), i.e. MG is the ideal in Q[K] generated by all binomials K u − K v where u, v ∈ Nn satisfy CayG (Y ) · (u − v) = 0. The codimension of this toric ideal equals the deficiency δ. Proof. Let Ids denote the s×s identity matrix and consider the extended Cayley matrix   −Ids Y1 Y2 · · · Yl     1 0 ··· 0   0     0 1 ··· 0 .  0   .. .. .   ..  0 . ..  . .   0 0 0 ··· 1 The toric ideal of this matrix is precisely the toric balancing ideal TG , where the unknowns c1 , c2 , . . . , cs correspond to the first s columns. Deleting these s columns corresponds to forming the elimination ideal MG as in (9). This shows that MG is the toric ideal of the matrix CayG (Y ). The dimension of the affine toric variety V (MG ) in Cn is equal to σ+l = rank(CayG (Y )), and hence its codimension equals the deficiency δ = n−σ−l. ✷ We conclude that V>0 (MG ) is a positive toric variety of codimension δ in Rn>0 . The moment map of toric geometry establishes a natural bijection between V>0 (MG ) and the interior of the Cayley polytope, which is the convex hull of the columns of CayG (Y ). In summary, given any chemical reaction network whose components are strongly connected, we have shown that the positive toric variety of the Cayley polytope equals the moduli space V>0 (MG ) of toric dynamical systems on G. The deficiency δ is precisely the codimension of this moduli space. In particular, if the deficiency is zero then the Cayley polytope is a simplex and (1) is toric for all rate constants κij . Moreover, the positive steady states of a toric dynamical system form a positive toric variety V>0 (TG ). 3. The Global Attractor Conjecture and Some Biological Applications We now consider a fixed toric dynamical system or, equivalently, a chemical reaction network (1) which admits a complex balancing state. The underlying directed graph G = (V, E) has n nodes labeled by monomials cy1 , cy2 , ..., cyn , and we specify positive rate constants by fixing a point κ0 in the moduli space V>0 (MG ). We also fix a strictly positive vector c0 ∈ Rs>0 which represents the initial concentrations of the s species. The equations (1) describe the evolution of these concentrations over time. We seek to understand the long-term behavior of the trajectory which starts at c0 , that is, c(0) = c0 . Let TG (κ0 ) denote the toric ideal in R[c] obtained from TG by substituting the specific rate constants κ0ij ∈ R>0 for the unknowns κij . Then V>0 (TG (κ0 )) coincides with the set of all steady states of the toric dynamical system (1). The following result is well-known: 8 Proposition 10. [Existence and Uniqueness of the Birch Point] The affine subspace c0 + S of Rs intersects the positive toric variety V>0 (TG (κ0 )) in precisely one point c∗ . For a proof and references in the chemistry literature see Horn and Jackson (1972); a different proof can be found in Feinberg (1979, Proposition 5.3) or Gunawardena (2003, Proposition 6.4). We remark that variants of Proposition 10 are ubiquitous across the mathematical sciences, and the result has been rediscovered many times. In statistics, this result is known as Birch’s Theorem; see (Pachter and Sturmfels, 2005, Theorem 1.10). To stress the link with toric models in algebraic statistics we call c∗ the Birch point of the toric dynamical system (1) with starting point c0 . The right hand side of (1) is always a vector in the stoichiometric subspace S = R{yj − yi : (i, j) ∈ E}. Hence the trajectory starting at c0 stays in the affine subspace c0 + S. In fact, concentrations remain non-negative, so the trajectory stays in P := (c0 + S)∩Rs≥0 . We call P the invariant polyhedron. Chemists use the term stoichiometric compatibility class for P . The relative interior of P in c0 + S is denoted by P o := (c0 + S) ∩ Rs>0 . Proposition 11. The Birch point c∗ is the unique point in the invariant polyhedron P for which the transformed entropy function E(c) = s X ci · log(ci ) − ci · log(c∗i ) − ci + c∗i i=1
(10) is a strict Lyapunov function of the toric dynamical system (1). This means the following: (a) For all c ∈ P we have E(c) ≥ 0 and equality holds if and only if c = c∗ , (b) we have dE(c)/dt ≤ 0 along any trajectory c(t) in P , and (c) equality in (b) holds at a point t of any trajectory c(t) in P ◦ if and only if c(t) = c∗ . This proposition was proved by Horn and Jackson (1972). A different proof can be found in (Feinberg, 1979); see especially Proposition 5.3 and its corollaries; see also (Gunawardena, 2003, Theorem 6.4) and the paragraph before it. We suggest comparing this with the proof of (Pachter and Sturmfels, 2005, Theorem 1.10). Any trajectory of the toric dynamical system (1) which starts in the relatively open polyhedron P o = (c0 + S) ∩ Rs>0 will stay in the closed polyhedron P = (c0 + S) ∩ Rs≥0 ; actually, it is not hard to show that P o is an invariant set. The main conjecture below states that any such trajectory converges to the Birch point. This conjecture was first formulated by Horn (1974). A steady state x in P o is called a global attractor if any trajectory that begins in P o converges to x. Global Attractor Conjecture. For any toric dynamical system (1) and any starting point c0 , the Birch point c∗ is a global attractor of the invariant set P o = (c0 + S) ∩ Rs>0 . An important subclass of toric dynamical systems consists of the chemical reaction networks of deficiency zero. If the deficiency δ = n − σ − l is zero then the moduli ideal MG is the zero ideal, by Theorem 9, and (1) is toric for all choices of rate constants. As remarked in the Introduction, the Global Attractor Conjecture is open even for deficiency zero systems. Our last section is devoted to partial results on the conjecture. First, however, we discuss biological examples which illustrate the concepts developed so far. 9 Example 12. [Networks with trivial moduli] We expect that our toric approach will be useful for parametric analyses of chemical reaction networks in systems biology. Analyses of this kind include (Kuepfer, Sauer and Parrilo, 2007), (Gnacadja et al., 2007) and (Sontag, 2001). Many of the explicit examples we found in the literature have trivial toric moduli in the sense that either MG is the unit ideal or MG is the zero ideal. If MG = h1i then (1) is never a toric dynamical system regardless of what the κij are. This happens when components of G are not strongly connected. Examples include Michaelis-Menten kinetics and the covalent modification cycle in (Gunawardena, 2003, §5). If MG = {0} then the network has deficiency zero and (1) is always a toric dynamical system, regardless of what the κij are. Examples include the cycle in (Kuepfer, Sauer and Parrilo, 2007, Equation (9)), the monotone networks in (De Leenheer, Angeli and Sontag, 2007), and the following network which is taken from (Gnacadja et al., 2007). The ligand-receptor-antagonist-trap network has s = 8 species and n = 8 complexes. This network G has four reversible reactions which we write in binomial notation: κ15 · c5 c6 − κ51 · c1 , κ26 · c6 c7 − κ62 · c2 , κ37 · c7 c8 − κ73 · c3 , κ48 · c8 c5 − κ84 · c4 . (11) Here l = 4 and σ = 4, so that δ = 0. In the algebraic notation of Section 2, the toric ideal TG equals the complex balancing ideal hΨ(c)·Aκ i and is generated by the four binomials in (11). Eliminating c1 , c2 , . . . , c8 as prescribed by (9) yields the zero ideal MG = {0}. ✷ Example 13. [DHFR catalysis] Here are some examples from systems biology which show a more complicated dynamical behaviour. We consider the reaction network in (Craciun, Tang and Feinberg, 2006, Figure 5); this reaction network has several positive equilibria for some values of the reaction rate parameters (see Craciun, Tang and Feinberg, 2006, Figure 7). This reaction network allows for inflow and outflow of some chemical species; in the language of deficiency theory, we say that one of the complexes of this reaction network is the zero complex (see Feinberg (1979)), i.e., one of the vectors yi is zero. Note that the group A of reactions in this network has almost the same structure as mechanism 6 in (Craciun, Tang and Feinberg, 2006, Table 1), shown below: E + S1 ⇋ ES1, E + S2 ⇋ ES2, ES1 + S2 ⇋ ES1S2 ⇋ ES2 + S1, (12) ES1S2 → E + P, S1 ⇋ 0, S2 ⇋ 0, P → 0. Like the more complicated DHFR catalysis network, the network (12) also has several positive equilibria for some values of the reaction rate parameters. It is easy to compute the deficiency of this simpler mechanism: the number of complexes is n = 12 (including the zero complex), the number of linkage classes is l = 4 (including the linkage class that contains the inflow and outflow reactions for the substrates S1, S2 and the product P ), and the dimension of its stoichiometric subspace is σ = 6. Therefore the deficiency of the network (12) is δ = 12 − 4 − 6 = 2. This network cannot be toric for any choice of the constant rates because it is not weakly reversible. If we make all reactions reversible in (12), then the complexes, the linkage classes, and the stoichiometric subspace do not change, so the deficiency of the reversible version of (12) is also 2. Example 14. [Recombination on the 3-cube] In population genetics (Akin, 1979, 1982), the evolution of a population is modeled by a dynamical system whose right 10 hand side is the sum of three terms, corresponding to mutation, selection and recombination. The contribution made by recombination alone is a quadratic dynamical system (Rabinovich, Sinclair and Wigderson, 1992) which can be written in the form (1). In our view, toric dynamical systems are particularly well-suited to model recombination. Here we consider a population of three-locus diploids, so the underlying genotope of the haploid gametes is the standard 3-dimensional cube (Beerenwinkel, Pachter and Sturmfels, 2007, Example 3.9). The eight vertices of the cube are the genotypes. They now play the role of the species in chemistry: s=8 genotypes frequencies [000] [001] [010] [011] [100] [101] [110] [111] c1 c2 c3 c4 c5 c6 c7 c8 . The recombination network G has n = 16 nodes corresponding to the pairs of genotypes which are not adjacent on the cube. There are twelve bidirectional edges, representing interactions, and we label them using the notation of (Beerenwinkel, Pachter and Sturmfels, 2007, Example 3.9). Six of the interactions correspond to conditional epistasis: [000] + [110] ↔ [010] + [100] κ1,2 · c1 c7 − κ2,1 · c3 c5 K1 = κ2,1 and K2 = κ1,2 [001] + [111] ↔ [011] + [101] κ3,4 · c2 c8 − κ4,3 · c4 c6 K3 = κ4,3 and K4 = κ3,4 [000] + [101] ↔ [001] + [100] κ5,6 · c1 c6 − κ6,5 · c2 c5 K5 = κ6,5 and K6 = κ5,6 [010] + [111] ↔ [011] + [110] κ7,8 · c3 c8 − κ8,7 · c4 c7 K7 = κ8,7 and K8 = κ7,8 [000] + [011] ↔ [001] + [010] κ9,10 · c1 c4 − κ10,9 · c2 c3 K9 = κ10,9 and K10 = κ9,10 [100] + [111] ↔ [101] + [110] κ11,12 · c5 c8 − κ12,11 · c6 c7 K11 = κ12,11 and K12 = κ11,12 . Secondly, we have marginal epistasis, giving rise to the six pairwise interactions among four complexes [000] + [111] [001] + [110] [010] + [101] [100] + [011] four monomials K13 · c1 c8 K14 · c2 c7 K15 · c3 c6 K16 · c4 c5 . Here K13 , K14 , K15 , K16 are cubic polynomials with 16 terms indexed by trees as in (7). By Proposition 3, they are the 3 × 3 minors of the Laplacian of the complete graph K4 :   κ13,14 +κ13,15 +κ13,16 −κ13,14 −κ13,15 −κ13,16       −κ14,13 κ14,13 +κ14,15 +κ14,16 −κ14,15 −κ14,16  .     −κ15,13 −κ15,14 κ15,13 +κ15,14 +κ15,16 −κ15,16   −κ16,13 −κ16,14 −κ16,15 κ16,13 +κ16,14 +κ16,15 The recombination network G has l = 7 connected components and its deficiency is δ = 5, as there n = 16 complexes, and the stoichiometric subspace S has dimension σ = 4. The 11 moduli ideal MG is minimally generated by 18 binomials. Twelve of them are cubics: K8 K11 K15 − K7 K12 K16 K6 K9 K15 − K5 K10 K16 K4 K11 K14 − K3 K12 K16 K2 K9 K14 − K1 K10 K16 K4 K7 K14 − K3 K8 K15 K2 K5 K14 − K1 K6 K15 K6 K12 K13 − K5 K11 K14 K2 K12 K13 − K1 K11 K15 K8 K10 K13 − K7 K9 K14 K4 K10 K13 − K3 K9 K15 K2 K8 K13 − K1 K7 K16 K4 K6 K13 − K3 K5 K16 . The remaining six generators of MG are quartics: K9 K11 K14 K15 − K10 K12 K13 K16 K6 K8 K13 K15 − K5 K7 K14 K16 K2 K4 K13 K14 − K1 K3 K15 K16 K5 K8 K10 K11 − K6 K7 K9 K12 K1 K4 K10 K11 − K2 K3 K9 K12 K1 K4 K6 K7 − K2 K3 K5 K8 . The moduli space (of toric dynamical systems on G) is the toric variety V (MG ) defined by these 18 binomials. It has codimension 5 and degree 56. For any recombination rates κ0 ∈ V>0 (MG ) and any starting point c0 in the population simplex ∆7 , the trajectory of the toric dynamical system (1) stays in the 4-dimensional polytope (c0 + S) ∩ ∆7 and is conjectured to converge to the Birch point c∗ . Akin (1979) calls c∗ the Wright point. It generalizes the classical Hardy-Weinberg equilibrium in the 2-locus system. ✷ 4. Detailed Balancing Systems In this section we discuss an important subclass of toric dynamical systems called detailed balancing systems. Here, every edge of the digraph G exists in both directions. We can  thus identify G = (V, E) with the underlying undirected graph G̃ = (V, Ẽ), where Ẽ = {i, j} : (i, j) ∈ E . For each undirected edge {i, j} ∈ Ẽ of the graph G̃ we define {i,j} {i,j} an n × n-matrix Aκ as follows. In rows i, j and columns i, j the matrix Aκ equals   −κ κij ,  ij κji −κji {i,j} and all other entries of the matrix Aκ Aκ are 0. The Laplacian of G decomposes as X A{i,j} . (13) = κ {i,j}∈Ẽ A detailed balancing system is a dynamical system (1) for which the algebraic equations {i,j} Ψ(c) · Aκ = 0 for {i, j} ∈ Ẽ admit a strictly positive solution c∗ ∈ Rs>0 . In light of (13), every detailed balancing system is a toric dynamical system, so the positive solution c∗ is unique and coincides with the Birch point. As it is for toric dynamical systems, the condition of being detailed balancing depends on the graph G̃ and the constants κij . We rewrite this condition in terms of binomials in Q[c, κ]. The two non-zero entries of {i,j} are κij cyi − κji cyj and its negative. Moreover, we find the row vector Ψ(c) · Aκ Ψ(c) · A{i,j} ·Y κ = (κij cyi − κji cyj ) · (yj − yi ), 12 and hence the right hand side of the dynamical system (1) can be rewritten as follows: X X (14) (κij cyi − κji cyj ) · (yj − yi ). Ψ(c) · A{i,j} ·Y = Ψ(c) · Aκ · Y = κ {i,j}∈Ẽ {i,j}∈Ẽ For a detailed balancing system, each summand in (14) vanishes at the Birch point c∗ . Example 15. We revisit Example 1. Let s = 2, n = 3 and G̃ the complete graph on three nodes labeled by c21 , c1 c2 and c22 . The dynamical system (2) is now written as d (c1 , c2 ) = (κ12 c21 −κ21 c1 c2 )·(−1, 1) + (κ13 c21 −κ31 c22 )·(−2, 2) + (κ23 c1 c2 −κ32 c22 )·(−1, 1). dt This is a detailed balancing system if and only if the following algebraic identities hold: κ212 κ31 − κ221 κ13 = κ223 κ31 − κ232 κ13 = κ12 κ32 − κ21 κ23 = 0. (15) This defines a toric variety of codimension two which lies in the hypersurface (3). ✷ To fit our discussion into the algebraic framework of Section 2, we now propose the following definitions. The detailed balancing ideal is the following toric ideal in Q[κ, c]:
(16) TeG := h κij cyi − κji cyj | {i, j} ∈ Ẽ i : (c1 c2 · · · cs )∞ . The corresponding elimination ideal in Q[κ] will be called the detailed moduli ideal: fG M := TeG ∩ Q[κ]. The ideal TeG is toric, by the same reasoning as in Proposition 6. The detailed moduli fG is a toric ideal of Lawrence type, as was the ideal in Example 2. Note, however, ideal M fG are toric in the original coordinates κij . Here, we did not that the ideals TeG and M need the transformation to the new coordinates K1 , . . . , Kn in (7). Using the ring inclusion Q[K, c] ⊂ Q[κ, c], we have the following inclusions of ideals: TG ⊆ TeG fG . and MG ⊆ M Here the equality holds precisely in the situation of Example 2, namely, when each chemical complex appears in only one reaction and each reaction is reversible. In general, as seen in Example 15, the corresponding inclusion of moduli spaces will be strict: fG ) ⊂ V>0 (MG ). V>0 (M In words: every detailed balancing system is a toric dynamical system but not vice versa. The following characterization of detailed balancing systems will be used in the next section. If L is any vector in Rs and c the unknown concentration vector then we write L ∗ c := (L1 c1 , L2 c2 , . . . , Ls cs ). Lemma 16. A toric dynamical system is detailed balancing if and only if all the binomials κij cyi −κji cyj in (16) have the form (L∗c)yi −(L∗c)yj , for some positive vector L ∈ Rs>0 . Thus, a detailed balancing system is a toric dynamical system of the special form X
dc = (17) (L ∗ c)yi − (L ∗ c)yj · (yj − yi ). dt {i,j}∈Ẽ 13 Proof. The if-direction is easy: if our binomials have the special form (L ∗ c)yi − (L ∗ c)yj {i,j} = 0. then c∗ = (1/L1 , 1/L2, . . . , 1/Ls ) is a positive solution to the equations Ψ(c)·Aκ Conversely, for the only-if direction, we define L as the reciprocal of the Birch point L = (1/c∗1 , 1/c∗2 , . . . , 1/c∗n ), and the result follows the fact that cyi −yj = (c∗ )yi −yj remains valid for all stationary points c of the system (1) as the starting point c(0) varies. ✷ We now fix a detailed balancing system (17) with a particular starting point c(0). Then the trajectory c(t) evolves inside the invariant polyhedron P = (c(0) + S) ∩ Rs≥0 . Consider any acyclic orientation E ′ ⊂ Ẽ of the graph G̃. This means that E ′ contains one from each pair of directed edges (i, j) and (j, i) in E, in such a way that the resulting directed subgraph of G has no directed cycles. The acyclic orientation E ′ specifies a stratum S inside the relatively open polyhedron P o = (c(0) + S) ∩ Rs>0 as follows:  S := c ∈ P o | (L ∗ c)yi > (L ∗ c)yj for all (i, j) in E ′ . The invariant polyhedron P is partitioned into such strata and their boundaries. We are interested in how the strata meet the boundary of P . Each face of P has the form FI := {c ∈ P | ci = 0 for i ∈ I} where I is subset of {1, 2, . . . , s}. This includes F∅ = P . Lemma 17. Consider a detailed balancing system (17) and fix an acyclic orientation E ′ of the graph G̃. If the closure of the stratum S corresponding to E ′ intersects the relative interior of a face FI of thePinvariant polyhedron P , then there exists a strictly positive ′ vector α ∈ RI>0 such that k∈I (yjk − yik ) · αk ≥ 0 for all directed edges (i, j) in E . P Proof. We proceed by contradiction: assume that the inequalities k∈I (yjk − yik )αk ≥ Programming Duality (Farkas’ 0 have no strictly positive solution α ∈ RI>0 . By Linear P Lemma), there is a non-negative linear combination v = (i,j)∈E ′ λij (yj − yi ) such that the following two conditions on v hold: (a) supp(v + ) ∩ I = ∅, and (b) supp(v − ) contains some j0 ∈ I. We shall prove the following two claims, which give the desired contradiction: Claim One: If c is a point in the relative interior of FI , then (L ∗ c)v+ > (L ∗ c)v− . Since (L ∗ c)i = 0 if and only if i ∈ I, and (L ∗ c)j > 0 for all j ∈ / I, (a) implies that (L ∗ c)v+ is strictly positive, while (b) implies that (L ∗ c)v− = 0, and we are done. Claim Two: If c is a point in the closure of the stratum S, then (L ∗ c)v+ ≤ (L ∗ c)v− . Consider any point s ∈ S. By the construction of v, the following equation holds: P Y
λij λij (yj −yi ) . (18) (L ∗ s)yj −yi = (L ∗ s)v = (L ∗ s) (i,j)∈E′ (i,j)∈E ′ Recall that (L ∗ s)yj −yi ≤ 1 for each oriented edge (i, j) ∈ E ′ . Also, each λij is nonnegative, so ((L ∗ s)yj −yi )λij ≤ 1. Using (18), this implies that (L ∗ s)v ≤ 1, and therefore (L ∗ s)v+ ≤ (L ∗ s)v− . By continuity we can replace s by c in this last inequality. ✷ The vector α ∈ RI>0 in Lemma 17 will play a special role in the next section. In Corollary 18 below we regard α as a vector in Rs≥0 by setting αj = 0 for all j ∈ {1, . . . , s}\I. Corollary 18. Let c(t) be a trajectory of a detailed balancing system (17) on the invariant polyhedron P , and suppose that a point c(t0 ) on this trajectory lies both in the closure of a stratum S and in the relative interior of a face FI of P . Let α ∈ Rs≥0 be the vector obtained as in Lemma 17. Then, the inner product h α, dc dt (t0 ) i is non-negative. 14 Proof. Let E ′ denote the orientation which specifies S. The velocity vector X
(L ∗ c(t0 ))yi − (L ∗ c(t0 ))yj · (yj − yi ). dc dt (t0 ) equals (i,j)∈E ′ Since c(t0 ) is in the closure of the stratum S, we have (L ∗ c(t0 ))yi − (L ∗ c(t0 ))yj ≥ 0. We also have hα, yj − yi i ≥ 0 because α comes from Lemma 17. This implies X
dc (L ∗ c(t0 ))yi − (L ∗ c(t0 ))yj · h α , yj − yi i ≥ 0. (t0 ) i = hα, dt ′ (i,j)∈E This is the claimed inequality. It will be used in the proof of Theorem 23. ✷ 5. Partial Results on the Global Attractor Conjecture This section contains what we presently know about the Global Attractor Conjecture which was stated in Section 3. This conjecture is proved for detailed balancing systems whose invariant polyhedron is bounded and of dimension two. We begin with some general facts on trajectories of toric dynamical systems, which are interesting in their own right. Consider a fixed toric dynamical system (1) with strictly positive starting point c(0) = c0 ∈ Rs>0 . The trajectory c(t) remains in the invariant polyhedron P = (c0 + S) ∩ Rs≥0 . Recall that any face of P has the form FI := {c ∈ P | ci = 0 if i ∈ I}, where I ⊆ {1, . . . s}. The boundary ∂P of P is the union of all faces FI where I is a proper subset of {1, . . . , s}. For positive ε, the ε-neighborhood in P of the boundary of P will be denoted by Vε (∂P ). We note that the transformed entropy function (10) can be extended continuously to the boundary of P , because ci log ci → 0 as ci → 0+ . Equivalent formulations of the following result are well known. For instance, see Siegel and Chen (1994); Sontag (2001). Proposition 19. Suppose that the invariant polyhedron P is bounded and the distance between the boundary of P and the set {c(t) ∈ P | t > 0} is strictly positive. Then the trajectory c(t) converges to the Birch point c∗ of P . Proof. We assume that c(t) does not converge to c∗ . Let ε > 0 be such that c(t) ∈ / Vε (∂P ) for all t > t0 . The strict Lyapunov function (10) ensures that there exists a neighborhood Vε′ (c∗ ) of the Birch point c∗ such that all trajectories that visit Vε′ (c∗ ) converge to c∗ . Then c(t) ∈ / Vε′ (c∗ ) for all t > t0 . Denote the complement of the two open neighborhoods by P0 := P \ (Vε (∂P ) ∪ Vε′ (c∗ )). Then the non-positive and continuous function c 7→ (∇E · dc dt )(c) does not vanish on P0 by Proposition 11, so it is bounded above by some −δ < 0 on P0 . Therefore, the value of E(c(t)) decreases at a rate of at least δ for all t > t0 , which implies that E is unbounded on P0 . This is a contradiction. ✷ Given a trajectory c(t) of (1), a point c̄ ∈ P is called an ω-limit point if there exists a sequence tn → ∞ with limn→∞ c(tn ) = c̄. Proposition 19 says that if the trajectory c(t) does not have any ω-limit points on the boundary of P , then it must converge to the Birch point c∗ . Thus, in order to prove the Global Attractor Conjecture, it would suffice to show that no boundary point of P is an ω-limit point. We first rule out the vertices. 15 Proposition 20. Let r be a vertex of P and consider any ε > 0. Then, there exists a neighborhood W of r such that any trajectory c(t) with starting point c(0) = co satisfying dist(c0 , r) > ε, does not visit W for any t > 0. Proof. The following set is the intersection of a closed cone with a sphere of radius one:   v V := | v ∈ S\{0} and r + v lies in P . kvk  Hence V is compact. We set I = j ∈ {1, . . . , s} : rj = 0 . For each v ∈ V, the ray γv (t) := r + tv extends from the vertex γv (0) = r into the polyhedron P for small t > 0. We consider how the transformed entropy function changes along such a ray: s X X X d log(c∗j vj ) vj log(rj + tvj ) − E(γv (t)) = vj (log(0 + tvj )) + dt i=1 j∈I j ∈I / = (Σj∈I vj ) · log(t) + w(t), where the function w(t) admits a universal upper bound for t close to 0 and v ∈ V. For each j ∈ I we have vj ≥ 0 because rj = 0 and r + tv ∈ P for small t > 0. Also, since v points into P , there exist j ∈ I with vj > 0. Thus, the function Σj∈I vj has a d E(γv (t)) tends to −∞ for t → 0. There exists positive minimum over V. It follows that dt t0 < ε such that for all v ∈ V the function t 7→ E(r + tv) decreases for 0 < t ≤ t0 . So, E(r) > µ := maxv∈V E(r + t0 v). On the other hand, E is continuous, so there is a neighborhood W of the vertex r (contained in {r + tv | t < t0 , v ∈ V}) such that E(c) > (E(r) + µ)/2 for all c ∈ W . Since E decreases along trajectories, we conclude that no trajectory c(t) that starts at distance ≥ ε from the vertex r can enter W . ✷ Remark 21. Chemical reaction networks for which P is bounded are called conservative. For conservative networks, there exists a positive mass assignment for each species that is conserved by all reactions (Feinberg, 1979). On the other hand, if 0 ∈ P , then the reaction network is not conservative. Proposition 20 ensures that, for a toric dynamical system, complete depletion of all the concentrations c1 , c2 , ..., cs is never possible. Lemma 22. Suppose that P is bounded and that the trajectory c(t) has an ω-limit point on the boundary of P . Then for any ε > 0 there exists a positive number tε > 0 such that c(t) belongs to Vε (∂P ) for all t > tε . In other words, the trajectory c(t) approaches the boundary. Proof. Suppose that for some ε > 0 there exists a sequence tn → ∞ such that c(tn ) ∈ / Vε (∂P ) for all n. As P is bounded, the trajectory c(t) has an ω-limit point p ∈ P \Vε (∂P ). On the other hand, c(t) also has an ω-limit point on the boundary of P . Consider a ball B2δ (p) of radius 2δ around p, whose closure lies fully in the relative interior of P . The trajectory c(t) enters and exits the neighborhood Bδ (p) of p infinitely many times, and also enters and exits the neighborhood P \B2δ (p) of the boundary infinitely many times. The trajectory c(t) travels repeatedly between these two sets which are at distance δ from each other. Note that |dc/dt| is bounded above, and ∇E · dc/dt is bounded away from zero on the annulus B2δ (p)\Bδ (p). Then, as in the proof of Proposition 19, each 16 traversal between the neighborhoods decreases the value of E(c(t)) by a positive amount that is bounded away from zero. This contradicts the fact that E is bounded on P . ✷ We shall now prove the main result of this section. Admittedly, Theorem 23 has three rather restrictive hypotheses, namely, “dimension two,” “bounded polyhedron,” and “detailed balancing.” At present we do not know how to remove any of these hypotheses. Theorem 23. Consider a detailed balancing system (17) whose stoichiometric subspace S = R{yj − yi | (i, j) ∈ Ẽ} is two-dimensional and assume that the invariant polygon P = (c0 + S) ∩ Rs≥0 is bounded. Then the Birch point c∗ is a global attractor for P . Proof. By Proposition 19, we need only rule out the possibility that the trajectory c(t) has an ω-limit point on the boundary of P . Proposition 20 gives the existence of open neighborhoods of the vertices such that no trajectory c(t) that starts outside them can visit them. Let V denote the union of these neighborhoods. Suppose now that c(t) has an ω-limit point on ∂P . That limit point lies in the relative interior of some edge F of P . Let Fε denote the set of points in P which have distance at most ε from the edge F . We claim that there exists ε > 0 and tε > 0, such that the trajectory c(t) remains in the subset Fε \V for all t > tε . This is true because c(t) belongs to the neighborhood Vε (∂P ) of the boundary for t ≫ 0, by Lemma 22, and hence c(t) belongs to Vε (∂P )\V for t ≫ 0. But this implies that c(t) belongs to Fε \V for t ≫ 0 because Fε \V is a connected component of Vε (∂P )\V for ε sufficiently small. This uses the dimension two assumption. Consider the closures of all strata S that intersect the relative interior of F . After decreasing ε if necessary, we may assume that the union of these closures contains the set Fε \V , which contains the trajectory c(t) for t > tε . To complete the proof, we will show that the distance from c(t) to the edge F never decreases after c(t) enters Fε \V . Any stratum S whose closure intersects the relative interior of F contributes a vector α = α(S) which satisfies the statement of Lemma 17 for F = FI . The orthogonal projection of α(S) into the two-dimensional stoichiometric subspace is a positive multiple of the unit inner normal α0 ∈ S to F in P . By Corollary 18 we have hα(S), dc dt (t)i ≥ 0 and hence hα0 , dc (t)i ≥ 0 for t > t . Therefore the distance from c(t) to F cannot ε dt decrease. This is a contradiction to the assumption that F contains an ω-limit point. ✷ References Akin, E. The Geometry of Population Genetics, Lecture Notes in Biomathematics, 31, Springer, New York, 1979. Akin, E. Cycling in simple genetic systems, J. Math. Biol., 13:305-324, 1982. Beerenwinkel, N, Pachter, L and Sturmfels, B. Epistasis and the shape of fitness landscapes, Statistica Sinica 17:1317-1342 (2007). Cox, D, Little, J and O’Shea, D. Ideals, Varieties and Algorithms, Undergraduate Texts in Mathematics, Springer Verlag, Third Edition, 2007. Craciun G, Tang Y, and Feinberg, M. Understanding bistability in complex enzymedriven reaction networks, Proc. Natl. Acad. Sci., 103:23, 8697-8702, 2006. De Leenheer, P, Angeli, D and Sontag, E. Monotone chemical reaction networks, J. Math. Chem., 41: 295-314, 2007. 17 Feinberg, M. Complex balancing in general kinetic systems, Arch. Rat. Mech. Anal., 49:3, 187-194, 1972. Feinberg, M. Lectures on chemical reaction networks. Notes of lectures given at the Mathematics Research Center of the University of Wisconsin in 1979, http://www.che.eng.ohio-state.edu/∼FEINBERG/LecturesOnReactionNetworks Feinberg, M. Chemical reaction network structure and the stability of complex isothermal reactors I. The deficiency zero and deficiency one theorems, Chem. Eng. Sci., 42:10, 2229-2268, 1987. Feinberg, M. Necessary and sufficient conditions for detailed balancing in mass action systems of arbitrary complexity, Chem. Eng. Sci., 44:9, 1819-1827, 1989. Feinberg, M. Existence and uniqueness of steady states for a class of chemical reaction networks, Arch. Rational Mech. Anal., 132:311-370, 1995. Gatermann, K. Counting stable solutions of sparse polynomial systems in chemistry, Contemporary Mathematics, Volume 286, Symbolic Computation: Solving Equations in Algebra, Geometry and Engineering, (Editors E. Green et al.), 53–69, 2001. Gatermann, K and Huber, B. A family of sparse polynomial systems arising in chemical reaction systems. J. Symbolic Comput., 33:3, 275-305, 2002. Gatermann, K and Wolfrum, M. Bernstein’s second theorem and Viro’s method for sparse polynomial systems in chemistry, Adv. in Appl. Math., 34:252-294, 2005. Gnacadja, G, Shoshitaishvili, A, Gresser, M, Varnum, B, Balaban, D, Durst, M, Vezina, C, and Li, Y. Monotonicity of interleukin-1 receptor-ligand binding with respect to antagonist in the presence of decoy receptor, J. Theor. Biol., 244:478-488, 2007. Gunawardena, J. Chemical reaction network theory for in-silico biologists. Technical Report, 2003, http://vcp.med.harvard.edu/papers/crnt.pdf Horn, F and Jackson, R. General mass action kinetics. Arch. Rat. Mech. Anal., 47:2, 81-116, 1972. Horn, F. Necessary and sufficient conditions for complex balancing in chemical kinetics, Arch. Rat. Mech. Anal., 49:3, 172-186, 1972. Horn, F. Stability and complex balancing in mass-action systems with three complexes. Proc. Royal Soc. A, 334: 331-342, 1973. Horn, F. The dynamics of open reaction systems. Mathematical aspects of chemical and biochemical problems and quantum chemistry, SIAM-AMS Proceedings, Vol. VIII, 125-137, 1974. Huber, B, Rambau, J and Santos, F. The Cayley trick, lifting subdivisions and the Bohne-Dress theorem on zonotopal tilings, J. Eur. Math. Soc., 2:179-198, 2000. Kuepfer, L, Sauer, U and Parrilo, P. Efficient classification of complete parameter regions based on semidefinite programming, BMC Bioinformatics, 8:12, 2007. Pachter, L and Sturmfels, B. Algebraic Statistics for Computational Biology, Cambridge University Press, Cambridge, 2005. Rabinovich, Y, Sinclair, A and Wigderson, A. Quadratic dynamical systems, Proc. 33rd Annual Symposium on Foundations of Computer Science (FOCS), 1992, 304–313. Siegel, D and Chen, S.F. Global stability of deficiency zero chemical networks, Canadian Appl. Math Quarterly, 2:413–434, 1994. Sontag, E. Structure and stability of certain chemical networks and applications to the kinetic proofreading model of T-cell receptor signal transduction, IEEE Trans. Automat. Control, 46: 1028–1047, 2001. Stanley, R. Enumerative Combinatorics, Volume 2, Cambridge University Press, 1999. Sturmfels, B. Gröbner Bases and Convex Polytopes, American Mathematical Society, University Lectures Series, Vol. 8, Providence, Rhode Island, 1996. 18