Translated chemical reaction networks.

Bull Math Biol (2014) 76:1081–1116 DOI 10.1007/s11538-014-9947-5 O R I G I N A L A RT I C L E

Translated Chemical Reaction Networks Matthew D. Johnston

Received: 23 July 2013 / Accepted: 21 February 2014 / Published online: 8 March 2014 © Society for Mathematical Biology 2014

Abstract Many biochemical and industrial applications involve complicated networks of simultaneously occurring chemical reactions. Under the assumption of mass action kinetics, the dynamics of these chemical reaction networks are governed by systems of polynomial ordinary differential equations. The steady states of these mass action systems have been analyzed via a variety of techniques, including stoichiometric network analysis, deficiency theory, and algebraic techniques (e.g., Gröbner bases). In this paper, we present a novel method for characterizing the steady states of mass action systems. Our method explicitly links a network’s capacity to permit a particular class of steady states, called toric steady states, to topological properties of a generalized network called a translated chemical reaction network. These networks share their reaction vectors with their source network but are permitted to have different complex stoichiometries and different network topologies. We apply the results to examples drawn from the biochemical literature. Keywords Chemical kinetics · Steady state · Mass action system · Complex balancing · Weakly reversible

1 Introduction Chemical reaction networks are given by sets of reactions which transform a set of reactants into a set of products at a given kinetic rate. Under the simplest of kinetic as-

Electronic supplementary material The online version of this article (doi:10.1007/s11538-014-9947-5) contains supplementary material, which is available to authorized users.

B

M.D. Johnston ( ) Department of Mathematics, University of Wisconsin-Madison, 480 Lincoln Dr., Madison, WI 53706, USA e-mail: [email protected]

1082

M.D. Johnston

sumptions, that of mass action kinetics, the dynamics of a continuously-mixed chemical process may be modeled as an autonomous system of polynomial ordinary differential equations called a mass action system. Despite the simplistic formulation of such systems, the resulting dynamical systems may exhibit a wide range of dynamical behaviors, including multistationarity (Craciun and Feinberg 2005, 2006), Hopf bifurcations (Wilhelm and Heinrich 1995, 1996), periodicity and chaos (Érdi and Tóth 1989). Particular attention has been given recently to the nature of the steady states of these mass action systems, and in particular to the positive steady states (that is to say, steady states in Rm >0 ). Such analysis is complicated by two main factors: (i) the nonlinear nature of the steady state equations, and (ii) the partitioning of the positive state space into invariant affine spaces called compatibility classes. The analysis is further complicated by the observation that, for applied chemical processes, many parameter values (i.e., the rate constants associated with each reaction) are typically unknown or only known to a certain precision; consequently, an emphasis has been placed on results which characterize the steady state set regardless of the rate constant values. Nevertheless, many general results about the steady states of mass action systems are well-known. It has been known since the 1970s that two fundamental classes of mass action systems—detailed balanced systems (Vol’pert and Hudjaev 1985) and complex balanced systems (Horn and Jackson 1972)—possess a unique positive steady state within each positive compatibility class. These results were further related to the topological structure of the network’s underlying reaction graph (reversibility and weak reversibility, respectively) in Feinberg (1972), Horn (1972). This network structure approach to characterizing steady states has been continued by Martin Feinberg in a series of papers focusing on network deficiency (Feinberg 1987, 1988, 1995b), network injectivity (Craciun and Feinberg 2005, 2006), and concordance (Shinar and Feinberg 2012). This author, together with Jian Deng, Christopher Jones, and Adrian Nachman, was also instrumental in producing a paper affirming the long-standing conjecture that every weakly reversible network contains a positive steady state (Deng et al. 2011). Beginning with a series of papers published by Karin Gatermann in the early 2000s, interest arose for characterizing the steady state sets of mass action systems by using tools from algebraic geometry (Gatermann 2001; Gatermann and Huber 2002; Gatermann and Wolfrum 2005). Other prominent algebraists, including Alicia Dickenstein and Bernd Sturmfels, have since become involved in adapting chemical reaction network results and terminology to this algebraic setting. These authors, along with Gheorghe Craciun and Anne Shiu, were instrumental in making the connection between toric varieties, Birch’s theorem from algebraic statistics, and complex balanced steady states in Craciun et al. (2009). This paved the way for the introduction of toric steady states, a generalization of complex balanced steady states which no longer shared any direct correspondence on the topological structure of the reaction graph (Pérez Millán et al. 2012). Other related contributions to the study of the steady states of mass action systems have been made in Clarke (1980), Conradi et al. (2008), Dickenstein and Pérez Millán (2011), Feinberg (1989), Flockerzi and Conradi (2008), Markevich et al. (2004).

Translated Chemical Reaction Networks

1083

Research on the steady states of chemical reaction systems has also been conducted for systems which do not possess traditional mass action kinetics. One recent example is that of generalized mass action systems introduced by Müller and Regensburger (2012). Generalized mass action systems maintain the topological structure of standard chemical reaction networks but allow the powers of the monomials appearing in the steady state conditions to differ from those implied by the network stoichiometry. The authors show that a notion of complex balancing is maintained in this generalized setting and that steady state properties can often still be inferred from the topological structure of the generalized reaction graph. This work also led to a generalization of Birch’s theorem. In this paper, we introduce a method for relating the steady states of a mass action system to those of a specially-constructed generalized mass action system. This method, called network translation, allows an explicit connection to be made between systems with toric steady states and generalized mass action systems with complex balanced steady states. It also allows steady state properties to be inferred from generalized network parameters. As such, this paper can be seen as a step toward closing the gap between the network topology approaches to characterizing steady states championed by Martin Feinberg et al., and the approaches of algebraists such as Karin Gatermann and Alicia Dickenstein. We apply the results to several well-studied networks contained in the biochemical literature. While the primary application of this paper is characterizing the steady states of mass action systems, it will be noted that translated chemical reaction networks are interesting objects of study in their own right. We will close with a discussion of some avenues for future research, both within the study of translated chemical reaction networks and generalized chemical reaction networks in general.

2 Background In this section, we present the terminology and notation relevant for the study of chemical reaction networks and mass action systems which will be used throughout this paper. We will present these concepts both in the standard and generalized setting. 2.1 Chemical Reaction Networks A chemical reaction network is given by a triple of sets N = (S, C, R) where: 1. The species set S = {A1 , . . . , Am } consists of the fundamental molecules capable of undergoing chemical change. 2. The complex set C= {C1 , . . . , Cn } consists of linear combinations of the species m of the form Ci = m j =1 yij Aj , i = 1, . . . , n. The terms yij ∈ Z≥0 are called stoichiometric coefficients and allow us to define the stoichiometric vectors yi = (yi1 , . . . , yim ), i = 1, . . . , n. 3. The reaction set R = {R1 , . . . , Rr } consists of interactions of the form Rk = Ci → Cj for some i, j = 1, . . . , n, i = j , for k = 1, . . . , r.

1084

M.D. Johnston

It is also typically assumed in chemical reaction network theory that: (a) every species in S appears in at least one complex in C; (b) every complex in C appears in at least one reaction in R; and (c) there are no self-reactions (i.e., reactions of the form Ci → Ci ). Remark 1 Note on indexing It will be necessary throughout this paper to retain elements of both reaction-centered indexing and complex-centered indexing. This differs from much of chemical reaction network theory literature. To simplify the relationships between these sets, when ambiguity is not a concern we will represent the sets C and R by their corresponding index sets, {1, . . . , n} and {1, . . . , r}, respectively. To further formalize the relationship between the reactions and complexes, we define the set of complexes which appear on the left (right) of at least one reaction to be the reactant (product) complex set CR (CP). We furthermore define the mappings ρ : R → CR and ρ : R → CP so that ρ(i) (ρ (i)) corresponds to the reactant (product) complex of the ith reaction. The mappings ρ and ρ will be called the reactant profile and product profile of N , respectively. This allows the reaction network N to be represented in the form N : Cρ(i) −→ Cρ (i) ,

i = 1, . . . , r.

For example, consider the reaction network 1

3

4

N : C1 C2 −→ C3 ←− C4 . 2

The reactant complex set is CR = {1, 2, 4} and the reaction profile is ρ(1), ρ(2), ρ(3), ρ(4) = (1, 2, 2, 4). Correspondingly, the product complex set is CP = {1, 2, 3} and the product profile is (ρ (1), ρ (2), ρ (3), ρ (4)) = (2, 1, 3, 3). 2.2 Reaction Graph and Deficiency Interpreting chemical reaction networks as interactions between stoichiometrically distinct complexes naturally gives rise to their interpretation as directed graphs G(V , E) where the vertices are the complexes (i.e., V = C) and the edges are the reactions (i.e., E = R). In the literature, this graph has been termed the reaction graph of a network (Horn and Jackson 1972). There are several properties of a network’s reaction graph of which we will need to be aware. We will say that a complex Ci is connected to Cj if there exists a sequence of complexes {Cμ(1) , . . . , Cμ(l) } such that Ci = Cμ(1) , Cj = Cμ(l) , and either Cμ(k) → Cμ(k+1) or Cμ(k+1) → Cμ(k) for all k = 1, . . . , l − 1. We will furthermore say that Ci is strongly connected to Cj if there is such a set where all reactions are in the forward direction, and a (potentially different) set where all reactions are in the backward direction. A linkage class is a maximal set of connected complexes and a


1085

strong linkage class is a maximal set of strongly connected complexes. For instance, consider the chemical reaction network N:

C1 → C2 C3

C4 C5 .

(1)

We have the linkage classes L1 = {1, 2, 3} and L2 = {4, 5}, and the strong linkage classes Λ1 = {1}, Λ2 = {2, 3} and Λ3 = {4, 5}. The number of linkage classes in a network will be denoted by . A chemical reaction network is said to be weakly reversible if the linkage classes and strong linkage classes coincide. For example, we can see that the network (1) is not weakly reversible. By contrast, the following network can easily be seen to be weakly reversible: N:

C1 −→ C2

C3

(2)

To each reaction i ∈ R, we associate the reaction vector yρ (i) − yρ(i) . These vectors keep track of the net stoichiometric change in the individual species as the result of the reaction. They also give rise to the stoichiometric subspace S = span{yρ (i) − yρ(i) | i = 1, . . . , r} ⊂ Rm . The dimension of the stoichiometric subspace is denoted s = dim(S). We may now define the following network parameter, which was introduced by Horn (1972) and Feinberg (1972). Definition 1 The deficiency of a chemical reaction network N is given by δ = n − − s. The deficiency is a nonnegative parameter which can be determined from the structure of the chemical reaction network itself and has been the study of significant research (Feinberg 1987, 1988, 1995a, 1995b). 2.3 Mass Action Systems In order to model how the concentrations of the chemical species evolve over time, we assume that the reaction vessel is spatially homogeneous and that the reacting species are in sufficient quantity to be modeled as chemical concentrations. We will furthermore assume that the system obeys mass action kinetics, so that the rate of each reaction is proportional to the product of concentrations of the reactant species. That is to say, if the ith reaction has the form A1 + A2 → · · · then we have the reaction rate = ki [A1 ][A2 ], where the proportionality constant ki is commonly called the rate constant of the reaction. We define k ∈ Rr>0 to be the vector of rate constants. In order to determine how the concentration vector x = (x1 , x2 , . . . , xm ) ∈ Rm ≥0 evolves over time, it is necessary to introduce the following matrices (notation adapted from Gatermann (2001), Gatermann and Huber (2002), Gatermann and Wolfrum (2005)):

1086

M.D. Johnston

– The complex matrix Y ∈ Zm×n ≥0 is the matrix where the j th column is the j th stoichiometric vector yj , i.e., [Y ]·,j = yj , j = 1, . . . , n. – The matrix Ia ∈ Zn×r ≥0 is the matrix with entries [Ia ]j i = −1 if ρ(i) = j , [Ia ]j i = 1 if ρ (i) = j , and [Ia ]j i = 0 otherwise. – The matrix Ik ∈ Rr×n ≥0 is the matrix with entries [Ik ]ij = ki if ρ(i) = j , and [Ik ]ij = 0 otherwise. We will also need the mass action vector Ψ (x) ∈ Rn≥0 , which is the vector with entries [Ψ (x)]j = xyj , j = 1, . . . , n. These definitions allow us to define the mass action system M = (S, C, R, k) governed by dx = Y Ia Ik Ψ (x). dt

(3)

It is worth noting that the stoichiometric matrix Γ := Y Ia ∈ Zm×r contains the reaction vectors yρ (i) − yρ(i) , i = 1, . . . , r, as its columns. It follows that trajectories of (3) are confined to affine translates of the stoichiometric subspace space. We therefore define the stoichiometric compatibility classes to be Cx0 = (S + x0 ) ∩ Rm ≥0 and , then x(t) ∈ C note that, if x(t) is a solution of (3) with x(0) = x0 ∈ Rm x 0 for all ≥0 t ≥ 0 (Horn and Jackson 1972; Vol’pert and Hudjaev 1985). Remark 2 The deficiency can also be defined in terms of Y and Ia . We have δ = dim(ker(Y ) ∩ Im(Ia )). The formula given in Definition 1 is equivalent (see Appendix A) but the formula given here will be more intuitive for the results of Sect. 5.1. 2.4 Generalized Mass Action Systems An alternative but related kinetic form to mass action kinetics is power-law formalism. In this formulation, the kinetic terms still have the product of concentrations form of monomials but are permitted to take (potentially non-integer) powers which do not necessarily correspond to the stoichiometry of the reactant complex (Savageau 1969). This has recently been extended by Müller and Regensburger (2012) to a more network-focused approach called generalized chemical reaction networks. Definition 2 A generalized chemical reaction network N = (S, C, CK , R) is a chemical reaction network (S, C, R) together with a set of kinetic complexes CK which are in one-to-one correspondence with the elements of C. The set (S, C, R) determines the reaction structure and stoichiometry of the generalized chemical reaction network, just as it does for a standard chemical reaction network; however, each complex in C is associated to a kinetic complex in CK . The kinetic complexes CK do not appear directly in the reaction graph but are called upon when determining the kinetics (i.e., the monomials in (3)). Remark 3 The subscript notation used here for the kinetic complex CK differs from ˜ that of Müller and Regensburger (2012), where the kinetic complexes are denoted C.


1087

This modification is made to avoid confusion with translated chemical reaction networks introduced in Sect. 4.1. Since the kinetic complexes are in one-to-one correspondence with the stoichiometric complexes, we may consider properties of a second reaction graph with the kinetic complexes CK in place of the regular complexes C (i.e., the reaction graph of (S, CK , R)). The hypothetical kinetic reaction graph does not determine the stoichiometry of the network but it does play an important role in determining where steady states of the corresponding kinetic model may lie. The kinetic-order subspace SK is defined as SK = span (yK )ρ (i) − (yK )ρ(i) | i = 1, . . . , r , (4) where the vectors (yK )j = ((yK )j 1 , (yK )j 2 , . . . , (yK )j m ), j = 1, . . . , n, correspond to the kinetic complexes (CK )j , j = 1, . . . , n. The kinetic complex matrix YK is the matrix where the j th column is the j th kinetic complex vector (yK )j , i.e., [YK ]·,j = (yK )j , j = 1, . . . , n. The deficiency δ of a generalized chemical reaction network is the same as given by Definition 1. We further define the kinetic deficiency to be the deficiency of the network (S, CK , R) and denote it by δK . When space is not a concern, the correspondence between the stoichiometric and kinetic complexes will be denoted by dotted lines in the reaction graph. For example, we write 7A1 + A3 · · ·

k1

A1 + A2 A3

· · · 5A2

(5)

k2

to imply that the stoichiometric complex C1 = A1 + A2 is associated with the kinetic complex (CK )1 = 7A1 + A3 and that the stoichiometric complex C2 = A3 is associated with the kinetic complex (CK )2 = 5A2 . The kinetic framework for generalized chemical reaction networks is the following. Definition 3 The generalized mass action system M = (S, C, CK , R, k) corresponding to the generalized chemical reaction network N = (S, C, CK , R) is given by dx = Y Ia Ik ΨK (x), dt

(6)

where Y , Ia , and Ik are as in (3), and the generalized mass action vector ΨK (x) has entries [ΨK (x)]j = x(yK )j , j = 1, . . . , n. In other words, a generalized mass action is the mass action system (3) with the monomials xyj replaced by the monomials x(yK )j . The generalized mass action system corresponding to network (5) is dx3 dx1 dx2 = =− = −k1 x17 x3 + k2 x25 , dt dt dt where the stoichiometry of the network comes from the stoichiometric complexes C but the monomials come from the kinetic complexes CK .

1088

M.D. Johnston

Remark 4 Müller and Regensburger also give sufficient conditions on the sign vectors associated with the stoichiometric and kinetic-order subspace S and SK for the existence of a unique (or multiple) generalized complex balanced steady state (states) within compatibility classes Cx0 = (x0 + S) ∪ Rm >0 (Theorem 3.10 and Proposition 3.2 of Müller and Regensburger (2012), respectively). They also give a generalization of the well-known Birch’s Theorem (Proposition 3.9, Müller and Regensburger 2012). These conditions will not be considered in this paper but have received further attention in Müller et al. (2013). Remark 5 It is worth noting that the kinetic complexes CK are not required to have the same support as that of the original complex set C in the generalized chemical reaction network framework. This contrasts with some general kinetic frameworks for chemical reaction systems, e.g., Angeli et al. (2007).

3 Steady States of Mass Action Systems When considering the steady states of a mass action system, we are typically interested in the positive steady state set given by E = x ∈ Rm (7) >0 | Y Ia Ik Ψ (x) = 0 . Characterizing (7) in general is a difficult algebraic task due to the nonlinear nature of the equations. It is somewhat surprising, therefore, that many characterizations exist within the literature which not only guarantee certain properties of the steady state set (7), but guarantee these properties for all compatibility classes and also for all rate constants. We will now consider two classifications of steady states which have appeared prominently in the literature: complex balanced steady states (Feinberg 1972; Horn 1972; Horn and Jackson 1972) and toric steady states (Craciun et al. 2009). 3.1 Complex Balanced Steady States The following class of steady states was introduced by Horn and Jackson (1972) as a generalization of detailed balanced steady states. Definition 4 A positive steady state x ∈ Rm >0 of a mass action system M = (S, C, R, k) is called a complex balanced steady state if Ψ (x) ∈ ker(Ak ), where Ak := Ia Ik ∈ Rn×n . Furthermore, a mass action system will be called a complex balanced system if every steady state is a complex balanced steady state. It is known that if a mass action system has a complex balanced steady state, then all steady states are complex balanced (Lemma 5B, Horn and Jackson 1972). Consequently, all mass action systems with complex balanced steady states are complex balanced systems. It is also known that every positive stoichiometric compatibility


1089

class Cx0 of a complex balanced system contains precisely one steady state (Lemma 5A, Horn and Jackson 1972), and the complex balanced steady state set is given by ⊥ E = x ∈ Rm , >0 | ln(x) − ln(a) ∈ S where a ∈ Rm >0 is an arbitrary complex balanced steady state of the system. The main result of Feinberg (1972) and Horn (1972), popularly called the Deficiency Zero Theorem, relates the capacity of a network to permit complex balanced steady states to properties of the reaction graph. Theorem 1 (Theorem 4A, Horn 1972) Every mass action system M = (S, C, R, k) admitted by a network N = (S, C, R) is complex balanced if and only if N is weakly reversible and has a deficiency of zero (i.e., δ = 0). This result gives computable properties, depending on the network topology alone, which are sufficient to guarantee strong restrictions on the nature, location, and number of steady states of the corresponding mass action systems. Remarkably, these results are guaranteed to hold for all possible choices of positive rate constants and all stoichiometric compatibility classes. A surprising observation of Müller and Regensburger (2012) is that complex balancing may also be meaningfully defined for generalized mass action systems. Definition 5 A positive steady state x ∈ Rm >0 of a generalized mass action system M = (S, C, CK , R, k) is called a generalized complex balanced steady state if ΨK (x) ∈ ker(Ak ), where Ak := Ia Ik ∈ Rn×n . The authors show that a generalized chemical reaction network which permits generalized complex balanced steady states is weakly reversible (Proposition 2.18, Müller and Regensburger 2012) and that the steady state set is given by ⊥ E = x ∈ Rm (8) >0 | ln(x) − ln(a) ∈ SK , where a ∈ Rm >0 is an arbitrary generalized complex balanced steady state of the system and SK is the kinetic-order subspace defined by (4). They also prove the following result. Theorem 2 (Proposition 2.20, Müller and Regensburger 2012) Every generalized mass action system (S, C, CK , R, k) admitted by a generalized chemical reaction network (S, C, CK , R) has at least one generalized complex balanced steady state if the underlying reaction network is weakly reversible and the kinetic deficiency of the network is zero (i.e., δK = 0). Remark 6 It is important to note that not all of the properties of standard complex balanced steady states apply in the generalized setting. In particular, the generalized complex balanced steady state set (8) may intersect a positive stoichiometric compatibility class Cx0 at a unique point, multiple times, or not at all.

1090

M.D. Johnston

3.2 Toric Steady States Algebraic techniques have recently become prominent in the study of steady state properties of mass action systems (Conradi et al. 2008; Craciun et al. 2009; Pérez Millán et al. 2012; Shiu 2010). We omit here background on the algebraic objects of interest, such as varieties, ideals, and Gröbner bases. The interested reader is directed to the accessible textbook of Cox, Little, and O’Shea (Cox et al. 2007). The connection between toric geometry, Birch’s theorem in algebraic statistics, and complex balanced steady states is made in Craciun et al. (2009). In particular, they show that the steady state ideal for any complex balanced system is a toric ideal. That is to say, it is a prime ideal which is generated by binomials. This justified the authors’ choice to refashion complex balanced systems as toric dynamical systems. The authors furthermore showed the following inclusion. (For a detailed introduction to the tree constants Ki and Kj (31), and the relationship between Theorem 3 and Definition 4, see Appendix B.) Theorem 3 (Corollary 4, Craciun et al. 2009) The steady state ideal of a complex balanced system contains the binomials Ki xyj − Kj xyi where i, j ∈ Lk for some k = 1, . . . , and Ki and Kj are the corresponding tree constants. There are a number of desirable features which follow from a system having a toric ideal. It allows, for instance, an easy parametrization of the associated variety. It was noted in Pérez Millán et al. (2012) that many mass action systems which do not admit complex balanced steady states nevertheless have steady state ideals which are generated by binomials. The authors say that such systems have toric steady states. In order to derive sufficient conditions for a system to have toric steady states, they introduce the complex-to-species matrix Σ := Y Ia Ik ∈ Rm×n . It is a surprising result that, for many non-complex balanced systems, ker(Σ) can, in fact, be decomposed in the same way as ker(Ak ) can be for complex balanced systems (see Appendix B). The authors show the following inclusion. Theorem 4 (Theorem 3.3, Pérez Millán et al. 2012) Consider a chemical reaction network N = (S, C, R). Suppose that ker(Σ) has dimension d and that there exists a partition Λ1 , Λ2 , . . . , Λd of {1, . . . , n} and a basis bk , k = 1, . . . , d, of ker(Σ) with supp(bk ) = Λk . Then the steady state ideal is generated by the binomials [bk ]i xyj − [bk ]j xyi for i, j ∈ Λk , k = 1, . . . , d. It is striking that the binomials in Theorem 3 and Theorem 4 are both constructed by first partitioning the complexes of the chemical reaction network into disjoint components and then computing a basis for a specific kernel restricted to the support of these components. Furthermore, the components in Theorem 3 have a clear interpretation in terms of the reaction network: they are the linkage classes of the network’s reaction graph. The components in Theorem 4 are less well-understood and are left as computational constructs in Pérez Millán et al. (2012). It is the clarification of the connection between Theorem 3 and Theorem 4, and of complex balanced steady states and toric steady states in general, which will be


1091

the main concern of this paper. We will show that the supports of the components derived in Theorem 4 can often be corresponded to linkage classes just as they are in Theorem 3. These linkage classes, however, will not be those of the original reaction network; rather, they will be the linkage classes of a related generalized reaction network which we will call a translated chemical reaction network. 4 Main Results In this section, we introduce the notion of network translation and show how this concept can be used to characterize mass action systems with toric steady states. To motivate the results, throughout this section we will consider the following pairs of example networks. In both cases, the first is a standard chemical reaction network, and the second is a generalized one. Example 1 Consider the chemical reaction network N1 and generalized chemical reaction network N˜ 1 given, respectively, by k1

A1 −→ 2A1 , k2

A1 + A2 −→ 2A2 ,

(Regular Network N1 )

k3

A2 −→ 0, and A1 · · ·

k1

0 −→ A1 k3 k2 A2 .. .

· · · A1 + A2 (Generalized Network N˜ 1 )

A2 The network N1 is a common network representation of the Lotka–Volterra predator–prey model, where A1 represents the prey and A2 represents the predator. It can be easily checked that the governing equations (3) of the mass action system M1 = (S, C, R, k) and the governing equations (6) of the generalized mass ˜ C, ˜ C˜K , R, ˜ k) ˜ 1 = (S, ˜ coincide. That is to say, they are dynamically action system M equivalent. We will give general conditions for such a property in Lemma 2. Example 2 Consider the chemical reaction network N2 and generalized chemical reaction network N˜ 2 given, respectively, by k1

k3

A1 A2 −→ A3 , k2 k4

A1 + A3 −→ A1 + A2 , k5

A2 + A3 −→ 2A2 ,

(Regular Network N2 )

1092

M.D. Johnston

and A1 · · ·

k1

A1

A2

k2

k3

A3

· · · A1 + A3

k 2

k4 + k1 k5

.. . A2

(Generalized Network N˜ 2 )

While the governing equations (3) of the mass action system M2 and the governing ˜ 2 do not coincide, it can be equations (6) of the generalized mass action system M easily checked that the systems have the same steady state set. In particular, they both yield the binomial Gröbner basis

k1 x1 − k2 x2 , (k1 k5 + k2 k4 )x1 x3 − k1 k3 x2

(9)

in the lexicographical ordering x3 > x2 > x1 . We will give general conditions for this correspondence of steady states in Lemma 4. In particular, we will give conditions for when such a correspondence can be made, and how the nontrivial rate constant k4 + (k1 /k2 )k5 was chosen. 4.1 Translated Chemical Reaction Networks The following is the fundamental new concept of this paper. Definition 6 Consider a chemical reaction network N = (S, C, R) with source com˜ with ˜ C, ˜ C˜K , R) plex set CR and a generalized chemical reaction network N˜ = (S, ˜ We will say N˜ is a translation of N if: source complex set CR. ˜ so that y˜ρ˜ (h (i)) − y˜ρ(h 1. There is a surjection h1 : R → R ˜ 1 (i)) = yρ (i) − yρ(i) for 1 all i = 1, . . . , r; ˜ so that h2 (ρ(i)) = ρ(h 2. There is a surjection h2 : CR → CR ˜ 1 (i)) for all i = 1, . . . , r; and ˜ → CR so that h2 (h3 (j )) = j and (C˜K )j = h3 (j ) for 3. There is an injection h3 : CR ˜ all j ∈ CR. The process of finding a generalized network N˜ which is a translation of N will be called network translation. Remark 7 We will use the tilde notation exclusively to denote properties of translations N˜ . For instance, we will let C˜ denote the stoichiometric complexes and C˜K denote the kinetic complexes of the translation N˜ . We will correspondingly let δ˜ and δ˜K denote the regular and kinetic deficiency of N˜ . This differs from the definitions given in Müller and Regensburger (2012), where C˜ was used to denote the set of kinetic complexes and δ˜ was use to represent the kinetic deficiency of a generalized chemical reaction network.


1093

˜ CR, and CR ˜ through the mappings Remark 8 The relationship between R, R, h1 , h2 , ρ, and ρ˜ can be visualized as

N:

h1 ˜ R −→ R ρ↓ ↓ ρ˜ h2 ˜ CR −→ CR,

: N˜

This representation is useful when interpreting property 2 of Definition 6. These formal conditions will be required for the technical results contained in Sects. 4.2 and 4.3. It will be helpful, however, to use the following intuition about the three properties given in Definition 6 when considering examples: 1. There is a correspondence between the reactions in N and N˜ which preserves the reaction vectors. 2. Every reaction from a common source complex in N must be mapped to a reaction from a common source complex in N˜ . In other words, reactions from a single source complex may not be broken apart by the translation. 3. The kinetic complexes are chosen from the set of source complexes of N which are mapped to the corresponding source complex of N˜ . Another useful way to think about a network translations N˜ is that it is a network produced by translating the reactions N by adding or subtracting species to the leftand right-hand sides of each reaction, while preserving the original complexes as the kinetic complexes of the new generalized network. This operation satisfies the three conditions of Definition 6. To investigate some subtleties which may arise, consider the examples introduced at the beginning of the section. Example 1 Consider translating the individual reactions of N1 in the following way: A1 −→ 2A1

(−A1 )

A1 + A2 −→ 2A2

(−A2 )

A2 −→ 0

0 −→ A1 , =⇒

A1 −→ A2 , A2 −→ 0.

(+0)

This satisfies properties 1 and 2 of Definition 6. In order to satisfy property 3, we associate to each source complex of the post-translation network the source complex of the original reaction. That is to say, we associate A1 to 0, A1 + A2 to A1 , and A2 to A2 . This produces the generalized network N˜ 1 so that we may say that N˜ 1 is a translation of N . Example 2 Consider translating the individual reactions of N2 in the following way: A1 A2 −→ A3

(+0)

A1 + A3 −→ A1 + A2

(−A1 )

A2 + A3 −→ 2A2

(−A2 )

A1 A2 −→ A3 , =⇒

A3 −→ A2 , A3 −→ A2 .

1094

M.D. Johnston

While this produces the reaction structure of N˜ 2 , there are two important differences. The first is that we have two reactions in N2 , A1 + A3 → A1 + A2 and A2 + A3 → 2A2 , which are associated with the reaction A3 → A2 in N˜ 2 . We have therefore reduced the number of reactions. This is allowed by property 1 but will require some additional consideration when transferring rate constants between the networks. A second important difference is that we may no longer unambiguously assign kinetic complexes to all of the source complex in N˜ 2 . This can be seen by noting that A1 + A3 and A2 + A3 are both mapped to A3 . This is allowed by property 2 of Definition 6, and property 3 allows us to choose either complex as the kinetic complex. Choosing A1 + A3 completes the correspondence, so that N˜ 2 is a translation of N2 by Definition 6. These examples demonstrate a fundamental difference between translations, namely, whether or not there is a one-to-one correspondence between the source complexes of the two networks. We therefore introduce the following important classification of translations. Definition 7 Consider a chemical reaction network N = (S, C, R) and a translation ˜ C, ˜ C˜K , R). ˜ We will say N˜ is a proper translation of N if h2 is injective as N˜ = (S, well as surjective. A translation N˜ will be called improper if it is not proper. ˜ Remark 9 If N˜ is proper then h3 = h−1 2 so that the kinetic complexes CK are uniquely defined in property 3 of Definition 6. A translation is improper if the kinetic complexes are not uniquely defined by the mapping h2 . Remark 10 Note that, for proper translations, h1 is necessarily injective as well as surjective. This follows from the observation that if two reactions in N are translated to the same reaction in N˜ there will necessarily be two source complexes in N mapped to the same source complex in N˜ . Examples 1& 2 It can be readily seen that N˜ 1 is a proper translation of N1 while N˜ 2 is an improper translation of N2 . It will also be important to understand how the stoichiometric and kinetic-order subspaces associated with N and its translation N˜ are related, and in particular in the case when N˜ is weakly reversible. Lemma 1 Consider a chemical reaction network N = (S, C, R) and a translation ˜ C, ˜ C˜K , R). ˜ Then the stoichiometric subspaces S of N and S˜ of N˜ coincide N˜ = (S, ˜ and, if N is weakly reversible, the kinetic-order subspace S˜K of N˜ is given by ˜ S˜K = span{yh3 (i) − yh3 (j ) | i, j ∈ L˜ k , k = 1, . . . , },

(10)

˜ are the linkage classes of N˜ . where L˜ k , k = 1, . . . , , ˜ be ˜ C, ˜ C˜K , R) Proof Let N = (S, C, R) be a chemical reaction network and N˜ = (S, ˜ a strong translation of N . By property 1 of Definition 6, N and N have the same


1095

reaction vectors and therefore have the same stoichiometric subspace S. This proves the first claim. To the second claim, it is well-known that the span of the reaction vectors of a chemical reaction network is the same as the span of the stoichiometric differences of complexes on the same connected component (for example, see Feinberg 1972, p. 189). Since the kinetic complexes are drawn from CR by h3 , this completes the proof. 4.2 Properly Translated Mass Action Systems We assign a kinetics to a proper translation in the following way. ˜ is a proper translation of a chemical re˜ C, ˜ C˜K , R) Definition 8 Suppose N˜ = (S, action network N = (S, C, R) and M = (S, C, R, k) is a mass action system corresponding to N . Then we define the properly translated mass action system of M ˜ k) ˜ = (S, ˜ C, ˜ C˜K , R, ˜ where k˜j = ki if to be the generalized mass action system M j = h1 (i). This relationship is the natural correspondence since we make the same correspondence for rate constants as we make for reactions. In other words, for proper translations, we will simply transfer the rate constant along with the reaction in the translation process. The following relates the dynamics of a properly translated mass action system ˜ to the original mass action system M. M ˜ k) ˜ = (S, ˜ C, ˜ C˜K , R, ˜ is a properly translated mass action sysLemma 2 Suppose M ˜ tem of M = (S, C, R, k). Then the generalized mass action system (6) governing M is identical to the mass action system (3) governing M. In particular, the two systems have the same steady states. Proof Consider a chemical reaction network N = (S, C, R) with corresponding ˜ ˜ C, ˜ C˜K , R) mass action system M = (S, C, R, k) and a proper translation N˜ = (S, ˜ ˜ ˜ ˜ ˜ ˜ of N . Let M = (S, C, CK , R, k) be a properly translated mass action system of M defined by Definition 8. Without loss of generality, we will index the reactions of N˜ so that h1 may be taken to be the identity. Since N˜ is a translation of N , it follows from property 1 of Definition 6 that the system (3) governing M is given by dx = Y Ia Ik Ψ (x) = Y˜ Iã Ik Ψ (x), dt where Y˜ and Iã correspond to the translation N˜ . It remains to relate the rate vector ˜ Notice that we have: ˜ := I˜k˜ Ψ˜ K (x) corresponding to M. R(x) := Ik Ψ (x) to R(x) – – –

˜ by property 3 of Definition 6; [ΨK (x)]j = xyh3 (j ) for all j ∈ CR ˜ki = ki , i = 1, . . . , r, by Definition 8; ˜ for all i = 1, . . . , r, by property 2 of Definition 6. h2 (ρ(i)) = ρ(i)

1096

M.D. Johnston

– h3 (h2 (j )) = j for all j ∈ CR by Condition 1 of Definition 7. ˜ ˜ = ki xyh3 (ρ(i)) = ki xyh3 (h2 (ρ(i))) = It follows that, for all i = 1, . . . , r, R˜ i (x) = k˜i xyh3 (ρ(i)) y ρ(i) ki x = Ri (x). Consequently, we have

dx = Y˜ Iã Ik Ψ (x) = Y˜ Iã I˜k˜ Ψ˜ K (x) dt ˜ have the same dynamics, and we are done. so that M and M

(11)

Example 1 We can now see that the rate constant choices for the proper translation N˜ 1 yield a translated mass action system by Definition 8. We have already observed that the mass action system (3) corresponding to M1 and the generalized mass action ˜ 1 coincide. This is exactly the conclusion guaranteed system (6) corresponding to M by Lemma 2. 4.3 Improperly Translated Mass Action Systems ˜ = It is more difficult to sensibly define a generalized mass action system M ˜ C, ˜ C˜K , R, ˜ k) ˜ C, ˜ C˜K , R) ˜ than a proper trans˜ for an improper translation N˜ = (S, (S, lation. This is because there is necessarily at least one source complex in N which does not appear as the kinetic complex of any source complex in N˜ . As a result, an ˜ will necessarily have fewer improperly translated generalized mass action system M monomials than the original mass action system M. A comprehensive dynamical result like Lemma 2 is consequently not possible. In this section, we introduce conditions, which we call resolvability conditions, ˜ correspondunder which the steady state set of a generalized mass action system M ˜ ing to an improper translation N may still be shown to coincide with the steady state set of the original mass action system M. We will use the networks N2 and N˜ 2 of Example 2 to illustrate the main concepts. We will need first of all to correspond the reactant complexes in N which are not used as kinetic complexes in N˜ to those which are. To that end, we introduce the following definitions. ˜ C, ˜ C˜K , R) ˜ is an improper translation of a chemical Definition 9 Suppose N˜ = (S, reaction network N = (S, C, R). Then: 1. The kinetically relevant complex of a reaction i ∈ R will be denoted ρ(i)K = h3 h2 ρ(i) .

(12)

2. A reaction i ∈ R will be said to be improperly translated if ρ(i) = ρ(i)K and the set of improperly translated reactions will be denoted RI = i ∈ R | ρ(i) = ρ(i)K . (13) 3. The improper subspace S˜I of N˜ will be defined to be S˜I = span{yρ(i) − yρ(i)K | i ∈ RI }.

(14)


1097

These definitions clarify several relationships unique to improper translations. A source complex is said to be kinetically relevant to a reaction if it is the kinetic complex in the translation. A reaction will then be called improperly translated if this kinetically relevant complex differs from its own source complex in N . Finally, the stoichiometric differences between the original and kinetically relevant complexes forms a basis of the improper subspace S˜I . The following result gives conditions under which the monomials corresponding to source complexes of improper reaction may be explicitly related to the corresponding kinetically relevant complexes. ˜ C, ˜ C˜K , R) ˜ is an improper translation of a chemical reLemma 3 Suppose N˜ = (S, action network N = (S, C, R). Suppose furthermore that N˜ is weakly reversible and that S˜I ⊆ S˜K where S˜K is the kinetic-order subspace of N˜ . Then, for every i ∈ RI , ˜ j = 1, . . . , s˜ , such that there exist constants cij and pairs pj , qj ∈ CR, x

yρ(i)

=

s˜ yh3 (pj ) cij

x j =1

x

yh3 (qj )

xyρ(i)K .

(15)

Proof Consider an arbitrary i ∈ RI and the difference yρ(i) − yρ(i)K . Since N˜ is weakly reversible, by Lemma 1 we can define a basis for S˜K by {yh3 (pj ) − yh3 (qj ) }sj˜ =1 ˜ by removing linearly dependent vectors from the generating set where pj , qj ∈ CR (10). Since S˜I ⊆ S˜K , it follows from (14) that we can write yρ(i) − yρ(i)K =

s˜

cij (yh3 (pj ) − yh3 (qj ) ),

(16)

j =1

where cij , j = 1, . . . , s˜ , are constants. The form (15) follows directly by rearranging (16) and raising the terms into the exponent of x ∈ Rm >0 . Lemma 3 gives conditions under which the monomials corresponding to source complexes in the original network may be explicitly related to the corresponding kinetic complexes. The result is deficient, however, in that the factors which relate the monomials are state dependent. The following stronger condition will allow us to simplify this dependence. ˜ C, ˜ C˜K , R) ˜ is an improper translation of a chemical Definition 10 Suppose N˜ = (S, reaction network N = (S, C, R). Suppose that N˜ is weakly reversible and that S˜I ⊆ S˜K where S˜K is the kinetic-order subspace of N˜ . Define K˜ i , i = 1, . . . , n, ˜ to be the ˜ C, ˜ R, ˜ k) ˜ with the reaction weights tree constants (31) associated with the network (S, k˜j = ki . (17) {i|h1 (i)=j }

Then:

1098

M.D. Johnston

1. For every i ∈ RI , define the kinetic adjustment factor of ρ(i) and ρ(i)K K˜ ρ(i),ρ(i)K =

s˜ ˜ cij

Kpj j =1

K˜ qj

,

(18)

˜ are determined according to (15) of where the constants cij and pairs pj , qj ∈ CR Lemma 3. 2. The translation N˜ will be called resolvable if, for every i ∈ RI , the kinetic adjustment factors K˜ ρ(i),ρ(i)K given in (18) do not depend explicitly on any kj , j ∈ RI . The condition of resolvability will be the key to removing the state dependence in (15). It will allow us to define a translated mass action system for an improper translation in the following way. Definition 11 Consider a chemical reaction network N = (S, C, R) and an associ˜ is a resolv˜ C, ˜ C˜K , R) ated mass action system M = (S, C, R, k). Suppose N˜ = (S, able improper translation of N . Consider the adjusted rate constants ki , for i ∈ / RI , (19) ki = (K˜ ρ(i),ρ(i)K )ki , for i ∈ RI . We define the improperly translated mass action system to be the generalized mass ˜ k) ˜ C, ˜ C˜K , R, ˜ with rate constants action system (S, ki . (20) k˜j = {i|h1 (i)=j }

This defines a generalized mass action system for improper translations which are also resolvable. In this case, we may sensibly define such a system by adjusting the rate constants by a carefully constructed kinetic adjustment factor which depends on the remainder of the rate constants. In order to see how Definition 11 applies in practice, reconsider Example 2. Example 2 We know the N˜ 2 is an improper translation of N2 . We now want to determine if it is resolvable and, if it is, what the corresponding improperly translated ˜ 2 is by Definition 11. mass action system M We start by defining the concepts in Definition 9. We have that the kinetic relevant complex ρ(i)K for the first four reactions of N2 coincides with their source complex ρ(i) in N2 . The fifth reaction, however, is assigned the kinetic complex A1 + A3 in the translation N˜ 2 rather than the original source complex A2 + A3 so that ρ(5) = ρ(5)K . It follows that RI = {5}. The improper subspace is given by S˜I = span{(0, 1, 1) − (1, 0, 1)} = span{(−1, 1, 0)}. We now want to apply Lemma 3. It is clear that N˜ 2 is weakly reversible. We have already determined S˜I , so it only remains to find the kinetic-order subspace S˜K . This ˜ ˜ C˜K , R): is given by the span of the reaction vectors of the network (S, A1 A2 A1 + A3 .


1099

It can be easily seen that S˜K = span{(−1, 1, 0), (1, −1, 1)} so that S˜I ⊂ S˜K . Allowing the basis elements to be represented in terms of kinetic reaction vectors corresponding to C˜K , we have that (16) gives yρ(5) − yρ(5)K = (0, 1, 1) − (1, 0, 1) = (1) (0, 1, 0) − (1, 0, 0) . After raising this into a power of x = (x1 , x2 , x3 ) ∈ R3>0 and rearranging, we arrive at the identity x2 x1 x3 x2 x3 = (21) x1 corresponding to (15). The important feature of (21) is that it explicitly relates the monomial x2 x3 , which corresponds to the unused complex A2 + A3 , to the monomial x1 x3 , which corresponds to the kinetically relevant complex A1 + A3 of the fifth reaction. The adjustment factor x2 /x1 is a ratio of monomials in the remainder of the set C˜K (the complexes A1 and A2 ). We now want to determine if the translation N˜ 2 is resolvable. We have already seen that the network N˜ 2 is weakly reversible and S˜I ⊆ S˜K . To compute the kinetic ˜ C, ˜ R, ˜ k) ˜ with reaction weights k˜i , adjustment factor (18) we consider the network (S, i = 1, . . . , 4, given by (17): k1

k3

k2

k4 +k5

A1 A 2 A 3 .

(22)

Setting C˜1 = A1 , C˜2 = A2 , and C˜3 = A3 , the required tree constants (31) are K˜ 1 = k2 (k4 + k5 )

and K˜ 2 = k1 (k4 + k5 )

so that K˜ ρ(5),ρ(5)K = K˜ 2 /K˜ 1 = k1 /k2 . Since this does not depend explicitly on the rate constant for the improper reaction, k5 , we have that N˜ 2 is a resolvable improper translation of N2 . We may now apply Definition 11 to construct the improperly translated mass ac˜ k). ˜ = (S, ˜ C, ˜ C˜K , R, ˜ The only difference between the rate constants tion system M of this network and (23) is that the rate constant k5 must be adjusted by the kinetic adjustment factor K˜ ρ(5),ρ(5)K = k1 /k2 . By (19) and (20), we define k˜1 = k1 , k˜2 = k2 , k˜3 = k3 , and k˜4 = k4 + (k1 /k2 )k5 to get: A1 · · ·

k1

A1

A2

k2

k3

A3

· · · A1 + A3

k 2

k4 + k1 k5

.. . A2 This coincides with the generalized mass action system given at the beginning of the section. We noted in the introduction of this section that the mass action system associated with N2 and the generalized mass action system associated with N˜ 2 shared the same

1100

M.D. Johnston

steady states. This is a general property of improperly translated mass action systems defined by Definition 11 with the additional condition δ˜ = 0 on the stoichiometric deficiency of the translation. This condition allows the steady state set to be completely characterized by tree constants, a requirement of the result. The following result is proved in Appendix C. ˜ C, ˜ C˜K , R) ˜ of a chemical reLemma 4 Consider an improper translation N˜ = (S, ˜ action network N = (S, C, R). Suppose that N is resolvable and δ˜ = 0. Let M = ˜ = (S, ˜ C, ˜ C˜K , R, ˜ k) ˜ (S, C, R, k) be a mass action system corresponding to N and M ˜ be an improperly translated mass action system corresponding to N and defined by ˜ coincide with Definition 11. Then the steady states of the system (6) governing M those of the system (3) governing M. Example 2 We can easily determine that the structural deficiency of N˜ 2 is zero ˜ C, ˜ C˜K , R, ˜ k) ˜ 2 = (S, ˜ de(i.e., δ˜ = 0) so that the generalized mass action system M fined by Definition 11 has the same steady state set as the mass action system M2 = (S, C, R, k) associated with N2 . 4.4 Connection with Toric Steady States We now make the following connection between translated mass action systems, complex balanced generalized mass action systems, and systems with toric steady states. The following is the main result of the paper and is proved in Appendix D. ˜ C, ˜ C˜K , R) ˜ is a weakly reversible translation of a chemTheorem 5 Suppose N˜ = (S, ical reaction network N = (S, C, R) which is either proper or improper and re˜ = (S, ˜ C, ˜ C˜K , R, ˜ k) ˜ denote a properly or improperly translated mass solvable. Let M action system corresponding to M = (S, C, R, k) (Definition 8 or Definition 11). Suppose furthermore that N˜ satisfies δ˜ = δ˜K = 0. Then: 1. The mass action system M has toric steady states for all rate constant vectors k ∈ Rr>0 . 2. The toric steady state ideal of M is generated by the binomials K˜ i xyh3 (j ) − K˜ j xyh3 (i) , ˜ and K˜ i , i = 1, . . . , n, where i, j ∈ L˜ k for some k = 1, . . . , , ˜ are the tree constants ˜ C, ˜ R, ˜ k). ˜ (31) corresponding to the reaction graph of (S, 3. The toric steady states of M correspond to the complex balanced steady states of ˜ and can be represented by M ˜⊥ E = x ∈ Rm >0 | ln(x) − ln(a) ∈ SK , ˜ ˜ where a ∈ Rm >0 is any generalized complex balanced steady state of M and SK is given by (10). Example 1 We have that N˜ 1 is weakly reversible and δ˜ = δ˜K = 0. Since N˜ 1 is a proper translation of N1 , it follows that M1 has toric steady states. To characterize


1101

˜ C, ˜ R, ˜ k) ˜ where the steady state ideal, we construct the weighted reaction graph (S, the reaction weights come from Definition 8: k1

0 −→ A1 k3 k2 A2 Setting C˜1 = 0, C˜2 = A1 , and C˜3 = A2 gives K˜ 1 = k2 k3 , K˜ 2 = k1 k3 , and K˜ 3 = k1 k2 . The monomials come from the mapping h3 which corresponds source complexes ˜ to their kinetic complexes C˜K . We have Ch (1) = A1 , Ch (2) = A1 + A2 , and CR 3 3 Ch3 (3) = A2 . It follows that the steady state ideal is K˜ 1 xyh3 (2) − K˜ 2 xyh3 (1) , K˜ 1 xyh3 (3) − K˜ 3 xyh3 (1) = k2 k3 x1 x2 − k1 k3 x1 , k2 k3 x2 − k1 k2 x1 .

Direct evaluation yields x1 = k3 /k1 and x2 = k3 /k2 , consistent with the known equilibrium of M1 . Example 2 We have that N˜ 2 is weakly reversible and δ˜ = δ˜K = 0. Since N˜ 2 is a resolvable improper translation of N2 , it follows by Theorem 5 that the mass action system M2 corresponding to N2 has toric steady states. In order to characterize the ˜ C, ˜ R, ˜ k) ˜ where the steady state ideal, we consider the weighted reaction graph (S, reaction weights come from Definition 11: k1

k3

A1 A2

k2

k4 + k1 k5

A3 .

(23)

k

2

Letting C˜1 = A1 , C˜2 = A2 , and C˜3 = A3 , the tree constants can be computed to be K˜ 1 = k2 k4 + k1 k5 , K˜ 2 = (k1 /k2 )(k2 k4 + k1 k5 ), and K˜ 3 = k1 k3 . We also have Ch3 (1) = A1 , Ch3 (2) = A2 , and Ch3 (3) = A1 + A3 . It follows that the ideal can be generated by K˜ 1 xyh3 (2) − K˜ 2 xyh3 (1) , K˜ 1 xyh3 (3) − K˜ 3 xyh3 (1) , which gives (k2 k4 + k1 k5 )x2 −

k1 (k2 k4 + k1 k5 )x1 , (k2 k4 + k1 k5 )x1 x3 − k1 k3 x1 . k2

This simplifies to (9). Also, since S˜K = span{(−1, 1, 0), (1, −1, 1)}, it follows that ⊥ = span{(1, 1, 0)} so that the steady set can be given by S˜K E = x ∈ R3>0 | ln(x) − ln(a) ∈ span (1, 1, 0) ˜ 2 . This can be directly where a ∈ R3>0 is any complex balanced steady state of M rearranged to give the parametrization E = (x1 , x2 , x3 ) ∈ R3>0 | x1 = et a1 , x2 = et a2 , x3 = a3 , t ∈ R .

1102

M.D. Johnston

5 Techniques and Applications We would like to apply Theorem 5 to characterize mass action systems which have toric steady states. This requires having a translation—either proper, or improper and resolvable—for which specific structural conditions are satisfied. In general, however, we do not have a translation N˜ of N given to us; rather, we must find the translation. Although not the primary purpose of this paper, in this section we present a heuristic method for generating translations which will then be successfully employed on several models from the literature. The method introduced here requires the introduction of a few further concepts (Clarke 1980; Gatermann and Wolfrum 2005). Definition 12 Consider a chemical reaction network N = (S, C, R) with associated complex matrix Y and matrix Ia as defined in Sect. 2.3. Then: 1. The flux cone of the network is given by Cv = {v ∈ Rm ≥0 | v ∈ ker(Y Ia )}. 2. A vector E ∈ Cv is called an elementary flux mode if E = λv + (1 − λ)w for any λ ∈ (0, 1) and v, w ∈ Cv , v = E, w = E. 3. An elementary flux mode E ∈ Cv is called a cyclic generator if E ∈ ker(Ia ). / 4. An elementary flux mode E ∈ Cv is called a stoichiometric generator if E ∈ ker(Ia ). The flux modes represent modes of positive stoichiometric flux balance in the reaction network while the cyclic generators have the additional interpretation of corresponding to cycles in the reaction graph of N . For example, the network 1

A1 2A1 2

3

5

4

6

A1 + A2 A3 A2

(24)

considered in Feinberg (1988) and Gatermann and Wolfrum (2005) has the elementary flux modes E1 = (1, 1, 0, 0, 0, 0), E2 = (0, 0, 1, 1, 0, 0), E3 = (0, 0, 0, 0, 1, 1), and E4 = (1, 0, 1, 0, 1, 0), of which E1 , E2 , and E3 are cyclic generators (reversible cycles in the reaction graph), and E4 is a stoichiometric generator. Every vector in the flux cone can be represented as a convex combination of the elementary flux modes (i.e., the cyclic and stoichiometric generators). Remark 11 Notice that the condition that a network has a deficiency of zero (i.e., δ = dim(ker(Y ) ∩ Im(Ia )) = 0) is equivalent to the condition that the network permits no stoichiometric generators. That is to say, if δ = 0, then every elementary flux mode of the network is a cyclic generator. Notice that (24) has δ = 1 as a result of the stoichiometric generator E4 . 5.1 Translation Algorithm In order to apply the Theorem 5, we require that the stoichiometric deficiency of the translation be zero (i.e., δ˜ = dim(ker(Y˜ ) ∩ Im(Iã )) = 0). This means that the translation N˜ does not have any stoichiometric generators whereas, if δ > 0, then


1103

original network N necessarily has stoichiometric generators. When attempting to apply Theorem 5, therefore, we have the following critical intuition: – The stoichiometric generators of N must be mapped into cyclic generators of N˜ through the translation process. We therefore propose the following algorithm for constructing translated chemical reaction networks. Translation Algorithm: 1. Identify the stoichiometric generators of the flux cone. 2. Assign an ordering to the reactions on the support of each stoichiometric generator found in part. If possible, add and/or subtract species to the left- and right-hand side of each reaction so that the product complex of each reaction coincides with the reactant complex of the next to create a cycle. 3. If possible, translate all reactions from the same source complex, or those already on the support of a cyclic generator, in unison (i.e., add/subtract the same term). 4. Repeat steps 2 and 3 until either successful, or all orderings of all stoichiometric generators has been exhausted. 5. If successful, retain the source complex of each reaction as the kinetic complex in the translation. In the case of multiple source complexes being assigned to a single translated complex, any of those source complexes may be chosen. 6. If the translation is proper, assign rate constants by Definition 8. If improper, determine whether it is resolvable (Definition 10) and, if so, assign rate constants by Definition 11. If successful, this algorithm produces a weakly reversible translated chemical reaction network by Definition 6 with no stoichiometric generators (i.e., δ˜ = 0) and an associated properly or improperly translated mass action system by either Definition 8 or 11. The algorithm is deficient in several ways. Notably, even if the algorithm works for some ordering of the reactions of a stoichiometric generator in step 2, it may not work for all choices. Certain orderings may allow multiple stoichiometric generators to fit together while certain others may not. Furthermore, even when the stoichiometric generators can be translated into cycles, they may fail to construct a consistent network by requirements of step 3. We will see in this section that the algorithm may nevertheless provide a valuable means by which to construct translations, and for characterizing the steady states of mass action systems by Theorem 5. The key observation is that this procedure allows an explicit characterization of the steady set (7) of M based on knowledge of the topological network structure of a related generalized network, N˜ (Müller and Regensburger 2012). This clarifies the connection between established deficiency-based approaches and the newer algebraic work contained in Pérez Millán et al. (2012). Further consideration of algorithmic approaches to network translation will be the subject of future work.

1104

M.D. Johnston

5.2 Application I: Futile Cycle Consider the enzymatic network N given by k1+

k2

S + E SE → P + E, k1−

(25)

k3+

k4

P + F P F → S + F. k3−

This network describes an enzymatic mechanism where one enzyme E catalyzes the transition of a substrate S into a product P , and a separate enzyme F catalyzes the reverse transition. Due to the appearance that the two enzymes are competing to undo the work of the other, this network is often called the futile cycle (Angeli and Sontag 2008; Pérez Millán et al. 2012; Wang and Sontag 2008). The steady states of this network under mass action (and more general) kinetics has been well-studied in the mathematical literature. The most thorough discussion is given in Angeli and Sontag (2008), where the authors show through a monotonicity argument that, for every choice of rate constants, every stoichiometric compatibility class Cx0 of (25) contains precisely one positive steady state and that this steady state is globally asymptotically stable relative to Cx0 . It has also be shown by the deficiency one algorithm (Feinberg 1988), the main argument of Wang and Sontag (2008), and Theorem 5.5 of Pérez Millán et al. (2012) that the network may not permit multistationarity. We show in the Supplemental Material that the Translation Algorithm can be successfully applied by using the translations k1+

k2

S + E SE → P + E

(+F ),

k1−

(26)

k3+

k4

P +F PF → S +F

(+E),

k3−

˜ C, ˜ C˜K , R) ˜ to give the proper translation N˜ = (S, S +E

k1+

· · · S + E + F SE + F k4

PF

↑

k1−

k3−

↓k2

··· PF + E P + E + F k3+

· · · SE,

··· P + F

˜ with the given rate constants. Since N˜ satiswith associated mass action system M ˜ ˜ fies δ = δK = 0, by claim 1 of Theorem 5 we have that M has toric steady states for all values of the rate constants. Further details characterizing the steady state set by Theorem 5 are contained in the Supplemental Material (ESM).


1105

5.3 Application II: Multiple Futile Cycle Now consider the multiple futile cycle N given by kon0

kcat0

S0 + E ES0 −→ S1 + E

lon0

koff0

loff0

.. . Sn−1 + E

konn−1

lcat0

S1 + F F S1 −→ S0 + F, .. .

kcatn−1

ESn−1 −→ Sn + E

Sn + F

koffn−1

lonn−1

lcatn−1

F Sn −→ Sn−1 + F.

loffn−1

(27) This network is a generalization of the futile cycle analyzed in Sect. 5.2. In this network, one enzyme E facilitates a forward cascade of transitions from substrate S0 to substrate Sn while another enzyme F facilitates the reverse transitions. This network has been frequently used in the literature to model multisite sequentially distributed phosphorylation networks of unspecified length (Gunawardena 2005; Holstein et al. 2013; Manrai and Gunawardena 2009). Despite the structural similarities with (25), there are notable dynamical differences in the corresponding mass action systems M. It was first shown in Markevich et al. (2004) that, even for the simple case n = 2, the system (27) admits rate constant vectors k ∈ Rr>0 for which M exhibits multistationarity. A subsequent focused study of the case n = 2 in Conradi et al. (2008) provided a detailed characterization of the rate parameters which permit multistationarity. The network (27) has also been studied for arbitrary values of n ≥ 1 in Gunawardena (2005, 2007), Pérez Millán et al. (2012), Wang and Sontag (2008). It is known that, for all n ≥ 2, the associated mass action systems M admit rate constant vectors k ∈ Rr>0 for which multistationarity is permitted and that an upper bound on the number of steady states in each compatibility class is 2n − 1 (Wang and Sontag 2008). It is furthermore conjectured that this upper bound can be tightened to n + 1 for odd n, and n for even n. In Pérez Millán et al. (2012), the authors prove that, for all rate constant vectors k ∈ Rr>0 , the mass action system M associated with this network has toric steady states. The authors then explicitly calculate the steady state ideal in terms of those rate constants. In the Supplemental Material we show that the Translation Algorithm can be applied by breaking the network into subcomponents of the form Si−1 + E

koni−1

kcati−1

ESi−1 −→ Si + E,

koffi−1

Si + F

loni−1

lcati−1

F Si −→ Si−1 + F,

loffi−1

1106

M.D. Johnston

and using the translations Si−1 + E

koni−1

kcati−1

ESi−1 −→ Si + E

(+iE + F ),

koffi−1

Si + F

loni−1

+(i + 1)E .

lcati−1

F Si −→ Si−1 + F

loffi−1

˜ = ˜ C, ˜ C˜K , R) This gives the proper translation N˜ = (S, by Si−1 + (i + 1)E + F

koni−1

n

˜ where N˜ i is given

i=1 Ni

ESi−1 + iE + F,

koffi−1 lcati−1

↑

↓kcati−1

(28)

loni−1

F Si + (i + 1)E Si + (i + 1)E + F, loffi−1

for i = 1, . . . , n. For each N˜ i , we associate the kinetic complexes Si−1 + E, ESi−1 , Si + F , and F Si , respectively, to the translated complexes in (28), starting in the top left and rotating clockwise. It can be computed that δ˜ = δ˜K = 0 for the translation N˜ (see Supplemental Material). Since the network is clearly weakly reversible, it follows by claim 1 of Theorem 5 that M has toric steady states for all rate constant vectors k ∈ Rr>0 . Further details characterizing the steady state set are contained in the Supplemental Material. 5.4 Application III: EnvZ/OmpR Signaling System Consider the signaling network N given by k1

k3

k2

k4

k5

XD X XT → Xp , k6

k8

Xp + Y Xp Y → X + Yp , k7 k9

(29) k11

XT + Yp XT Yp → XT + Y, k10 k12

k14

XD + Yp XDYp → XD + Y. k13

This network was first considered in Example (S60) of the Supporting Online Material of Shinar and Feinberg (2010) as a model for the EnvZ/OmpR signaling system is Escherichia coli. The mass action systems M associated with N were shown in that paper to exhibit concentration robustness in the species Yp . That is to say, the steady state value of [Yp ] were shown to be independent of the overall molar concentrations


1107

of X, Y , and their derivatives. The network was reproduced in Pérez Millán et al. (2012) where the authors showed that the systems M have toric steady states for all rate constant values. We now show that the steady states of M can be characterized by appealing to network translation and Theorem 5. Details of application of the Translation Algorithm are provided in the Supplemental Material. We start by relabeling the species as in Pérez Millán et al. (2012): A1 = XD,

A2 = X,

A3 = XT ,

A6 = Xp Y,

A7 = Yp ,

A4 = Xp ,

A8 = XT Yp ,

A5 = Y,

A9 = XDYp .

We then translate each linkage class in the following way: 1

3

2

4

5

(+A1 + A3 + A5 ),

8

(+A1 + A3 ),

11

(+A1 + A2 ),

14

(+A2 + A3 )

A1 A2 A3 → A4 6

A4 + A5 A6 → A2 + A 7 7

9

A3 + A7 A8 → A3 + A 5 10

12

A1 + A7 A9 → A1 + A 5 13

˜ where we have labeled the reac˜ C, ˜ C˜K , R), to get the following translation N˜ = (S, tions as they correspond to (29): 1

3

2

4

2A1 + A3 + A5 A1 + A2 + A3 + A5 A1 + 2A3 + A5 14 A2 + A3 + A9 13 12

↑11

↓5

A1 + A2 + A8 9

A1 + A3 + A4 + A5

↑↓ 10

7

(30)

↑↓ 6

8

A1 + A2 + A3 + A7 ← A1 + A3 + A6 For all complexes except A1 + A2 + A3 + A7 , we assign the kinetic complexes associated to each complex in N˜ to be the pre-translation source complex in N . We notice, however, that the source complexes of R9 and R12 (A3 + A7 and A1 + A7 , respectively) are both translated to A1 + A2 + A3 + A7 . We may choose either original reactant complex to be the associated kinetic complex so we will arbitrarily choose A3 + A7 . Since A1 + A7 does not appear as the kinetic complex for any source complex in N˜ , this is an improper translation. We show in the Supplementary Material that N˜ is a resolvable improper translation and that the rate constants from (19) and (20) of

1108

M.D. Johnston

Definition 11 are k˜i =

ki , 5) ( k2 (kk14k+k )k12 , 3

for i ∈ R \ {12}, for i = 12.

Since N˜ is a resolvable improper translation of N , is weakly reversible, and that δ˜ = δ˜K = 0 (easily checked), by claim 1 of Theorem 5, M has toric steady states for all rate constant values. Further details characterizing the steady state set are contained in the Supplemental Material.

6 Conclusions and Future Work In this paper, we introduced the notion of a translated chemical reaction network as a method for characterizing the steady states of mass action systems. The method of network translation relates a chemical reaction network N = ˜ C, ˜ C˜K , R), ˜ called a (S, C, R) to a generalized chemical reaction network N˜ = (S, translation of N , which has the same reaction vectors as N but different complexes and consequently different connectivity properties in the translated reaction graph. We defined two classes of translations, proper translations (Definition 7) and resolvable improper translations (Definition 10), which allowed a translated mass action ˜ = (S, ˜ C, ˜ C˜K , R, ˜ k) ˜ to be defined (Definition 8 and Definition 11, respecsystem M tively). We then presented conditions on the network topology of N˜ which allowed ˜ and an explicit connection to be made between complex balanced steady states of M toric steady states of M (Theorem 5). Finally, in Sect. 5 and the corresponding Supplemental Material, we applied the results to a series of examples drawn from the literature. The study of translated chemical reaction networks specifically, and generalized chemical reaction networks in general, is very new and there are consequently many aspects of the theory which have not be fully investigated. A few of the key points of future work include: 1. The Translation Algorithm presented in Sect. 5.1 depends heavily on intuition which may be lacking for large-scale biochemical networks. A stronger algorithm, and computational implementation, is required for broad-based application. 2. There is notable room for improvement in the conditions for resolvability of improper translations (Definition 10). In particular, computing the ratios K˜ pj /K˜ qj in (18) may be tedious. The author suspects that there are simpler sufficient conditions for resolvability of improper translations. 3. Translated chemical reaction networks are generalized chemical reaction networks, and consequently conclusions may only be drawn as far as they are justified by this underlying theory. The author suspects that, as this nascent theory becomes more fully developed, there will be increased application for the process of network translations in characterizing the steady states of mass action systems.


1109

Acknowledgements The author is supported by NSF grant DMS-1009275 and NIH grant R01GM086881. The author is also grateful for the numerous constructive conversations with Anne Shiu, Carsten Conradi, Casian Pantea, Stefan Müller, and others, over email and at the AIM workshop “Mathematical problems arising from biochemical reaction networks,” which pointed him toward the strong connection between toric steady states and complex balancing in generalized mass action systems. The author also thanks the two anonymous referees whose suggestions have significantly improved the paper.

Appendix A: Deficiency Result Lemma 5 The deficiency δ = n − − s of a chemical reaction network N , where n is the number of stoichiometrically distinct complexes, is the number of linkage classes, and s = dim(S), also satisfies δ = dim(ker(Y ) ∩ Im(Ia )). Proof It follows from basic dimensional considerations that dim ker(Y Ia ) = dim ker(Ia ) + dim ker(Y ) ∩ Im(Ia ) . From the rank–nullity theorem we have dim ker(Y Ia ) = r − dim Im(Y Ia ) = r − s. The rank of Ia corresponds to the number of complexes minus the number of linkage classes, so that dim(Im(Ia )) = n − . It follows that dim ker(Ia ) = r − (n − ) = r + − n. It follows that δ = dim ker(Y ) ∩ Im(Ia ) = dim ker(Y Ia ) − dim ker(Ia ) = (r − s) − (r + − n) = n − − s,

and we are done. Appendix B: Kernel of Ak

In this appendix, we present a more detailed characterization of ker(Ak ) for a mass action system M = (S, C, R, k). Consider a weakly reversible chemical reaction network N = (S, C, R) and let Lk , k = 1, . . . , , denote the network’s linkage classes. Define a subgraph T ⊆ R to be a spanning i-tree if T spans all of the complexes in some linkage class Lk , contains no cycles, and has the unique sink i ∈ C. Let Ti denote the set of all spanning i-trees for i = 1, . . . , n. We define the following network constants. Definition 13 Consider a weakly reversible chemical reaction network N = (S, C, R) with reaction weights kj , j = 1, . . . , r. Then the tree constant for i = 1, . . . , n is given by

kj . (31) Ki = T ∈T i j ∈T i

1110

M.D. Johnston

Remark 12 To compute the tree constants Ki , we restrict ourselves to the linkage class containing the complex i ∈ C. We then determine all of the spanning trees which contain this complex as the unique sink, multiply across all the weighted edges in each tree, and then sum over all such trees. The terms Ki can also be computed by computing specific minors of the kinetic matrix Ak restricted to the support of the linkage classes (Proposition 3, Craciun et al. 2009). Note that the term “tree constant” is our own. The following result characterizes ker(Ak ) in terms of the tree constants (31). This result appears in various forms within the chemical reaction network literature. A basic form, just concerned with the signs of the individual components, can be found in Feinberg (1979) (Proposition 4.1) and Gatermann and Huber (2002) (Theorem 3.1). A more specific result can be obtained by the Matrix-Tree Theorem (Stanley 1999). This form is explicitly connected with the reaction graph of a chemical reaction network in Craciun et al. (2009) (Corollary 4). A direct argument is also contained in Sect. 3.4 of Johnston (2011). We defer to these references for the proof. Theorem 6 Let N = (S, C, R) denote a weakly reversible chemical reaction network. Let Ki denote the tree constants (31) corresponding to i = 1, . . . , n. Then ker(Ak ) = span{K1 , K2 , . . . , K }, where Kj = ([Kj ]1 , [Kj ]2 , . . . , [Kj ]n ) has entries Ki , if i ∈ Lj , [Kj ]i = 0 otherwise. Remark 13 This theorem may be extended to networks which are not weakly reversible by considering the terminal strongly linked components of a chemical reaction network. As all the relevant networks considered in this paper are weakly reversible, however, Theorem 6 will suffice for our purposes here.

Appendix C: Proof of Lemma 4 ˜ C, ˜ C˜K , R) ˜ of a chemical reaction Proof Consider an improper translation N˜ = (S, network N = (S, C, R) which is resolvable. Let M = (S, C, R, k) be a mass action ˜ = (S, ˜ C, ˜ C˜K , R, ˜ k) ˜ be an improperly translated system corresponding to N and M mass action system corresponding to N˜ and defined by Definition 11. We will write the steady state condition for M as Y Ia Ik Ψ (x) = 0

(32)

˜ as and the steady state condition for M Y˜ Iã I˜k˜ Ψ˜ K (x) = 0.

(33)


1111

Since N˜ is improper, the vector Ψ˜ K (x) contains a subset of the monomials in Ψ (x) by property 3 of Definition 6. Consequently, to relate (32) and (33), we need to remove explicit dependence on the monomials in Ψ (x) corresponding to complexes not in C˜K . We will accomplish this by rewriting the unused monomials in Ψ (x) in terms of the monomial in Ψ˜ (x) corresponding to the kinetically relevant complex, and absorbing the required adjustment factor into the matrix Ik . (This will produce a state-dependent matrix I˜k(x) ˜ . We will resolve the state-dependency at a later stage.) ˜ Since N is resolvable, it follows that it is weakly reversible and S˜I ⊆ S˜K . Consequently, by Lemma 3, it follows that, for every i ∈ RI , there are constants cij and ˜ j = 1, . . . , s˜ , such that pairs pj , qj ∈ CR, x

yρ(i)

=

s˜ yh3 (pj ) cij

x j =1

x

yh3 (qj )

xyρ(i)K .

(34)

We now introduce ki (x) =

ki , s˜

xyh3 (pj ) cij ki , j =1 yh3 (qj ) x

and define k˜j (x) =

ki (x)

for i ∈ / RI , for i ∈ RI

(35)

(36)

{i|h1 (i)=j } ˜ ×n˜ ∈ Rr≥0 for j = 1, . . . , n. ˜ This gives rise to the state dependent kinetic matrix I˜k(x) ˜ with entries [I˜˜ ]ij = k˜i (x) if ρ(i) ˜ = j and [I˜˜ ]ij = 0 otherwise. k(x)

k(x)

Ψ˜ K (x) by showing that, for all We now will prove that Y Ia Ik Ψ (x) = Y˜ Iã I˜k(x) ˜ j = 1, . . . , r˜ , [Γ˜ ]·,j · R˜ j (x) = [Γ ]·,i · Ri (x), {i|h1 (i)=j }

˜ Ψ˜ K (x). We first where Γ := Y Ia , Γ˜ := Y˜ Iã , R(x) := Ik Ψ (x), and R(x) := I˜k(x) ˜ make several observations, listed in the order they will be used: – By definition, [Γ ]·,i = yρ (i) − yρ(i) and [Γ˜ ]·,j = y˜ρ˜ (j ) − y˜ρ(j ˜ ). – For h1 (i) = j , we have ρ(j ˜ ) = h3 (ρ(h ˜ 1 (i))) = h3 (h2 (ρ(i))) = ρ(i)K (since yρ(j ˜ ) = ρ(h ˜ 1 (i)) = h2 (ρ(i)) by property 2 of Definition 6) so that [Ψ˜ K (x)]ρ(j ˜ )=x yρ(i)K x . – By the constructions (34) and (35), we have ki (x) xyρ(i)K = ki xyρ(i) = Ri (x) for all i = 1, . . . , r. – For every i = 1, . . . , r, we have y˜ρ˜ (h1 (i)) − y˜ρ(h ˜ 1 (i)) = yρ (i) − yρ(i) by property 1 of Definition 6. It follows that, for every j = 1, . . . , r˜ , we have

1112

M.D. Johnston

˜ ) = (y ˜ yρ(j [Γ˜ ]·,j · R˜ j (x) = (y˜ρ˜ (j ) − y˜ρ(j ˜ρ˜ (j ) − y˜ρ(j ˜ ) ) kj x ˜ ))

=

{i|h1 (i)=j }

{i|h1 (i)=j }

(yρ (i) − yρ(i) ) ki xyρ(i) =

ki (x) xyρ(i)K

[Γ ]·,i · Ri (x).

{i|h1 (i)=j }

It follows that we have Y Ia Ik Ψ (x) = Y˜ Iã I˜k(x) Ψ˜ K (x) = Y˜ A˜ k(x) Ψ˜ K (x), ˜ ˜

(37)

˜ n˜ := Iã I˜k(x) ∈ Rn× where A˜ k(x) ˜ ˜ >0 is a state dependent kinetic matrix with positive offdiagonal entries corresponding to the structure of the translation N˜ and rates given by (35) and (36). ˜ C, ˜ R, ˜ k(x)) ˜ Consider the reaction graph of the network (S, with state dependent r ˜ ˜ edge weights k(x) ∈ R≥0 given by (35). In order to remove the state dependence in ˜ ˜ ˜ A˜ k(x) ˜ , we consider the system at steady state. Since δ = dim(ker(Y ) ∩ Im(Ia )) = 0, it follows that

Y˜ A˜ k(x) Ψ˜ K (x) = 0 ˜

A˜ k(x) Ψ˜ K (x) = 0. ˜

⇐⇒

(38)

˜ denote the state dependent tree constants (31) of the Now let K˜ j (x), j = 1, . . . , n, ˜ ˜ ˜ ˜ reaction graph of (S, C, R, k(x)). Since N˜ is weakly reversible, by Theorem 6, we have that ˜ ˜ ˜ ker(A˜ k(x) ˜ ) = span K1 (x), K2 (x), . . . , K˜ (x) , ˜ j (x) = ([K˜ j (x)]1 , [K˜ j (x)]2 , . . . , [K˜ j (x)]n˜ ) has entries where K K˜ i (x) if i ∈ L˜ j , K˜ j (x) i = 0 otherwise. ˜ we have It follows that, if i, j ∈ L˜ k for some k = 1, . . . , , xyh3 (j ) xyh3 (i) = K˜ i (x) K˜ j (x)

xyh3 (i) K˜ i (x) . yh3 (j ) = ˜ x Kj (x)

=⇒

(39)

Now consider i ∈ RI and let cij and pj , qj , j = 1, . . . , s˜ , denote the values guaranteed by Lemma 3. By (39), for every i ∈ RI , we have s˜ yh3 (pj ) cij

x j =1

x

yh3 (qj )

=

s˜ ˜

Kpj (x) cij j =1

K˜ qj (x)

.

(40)

It follows by (35) and the assumption that N˜ is resolvable that (40) only depends / RI . It follows that (40) may be written on the state-independent rates ki (x) = ki , i ∈ s˜ yh3 (pj ) cij

x j =1

x

yh3 (qj )

=

s˜ ˜ cij

Kpj j =1

K˜ qj

,

(41)


1113

where the tree constants K˜ i are determined with respect to the reaction graph of ˜ C, ˜ R, ˜ k) ˜ with the rate constants given by (35) and (36) for k = ki for i ∈ (S, / RI and i ki arbitrary for i ∈ RI (since the product (40) does not depend on these rates). We may now substitute (41) into (35) to get ⎧ for i ∈ / RI , ⎨ ki , ki (x) = ki = s˜ K˜ pj cij (42) ki , for i ∈ RI ⎩ j =1 K˜ qj

and

k˜j (x) = k˜j =

ki .

(43)

{i|h1 (i)=j }

Notice that, while the tree constant pairs K˜ pi and K˜ qi depend upon a choice for the rate constants ki for i ∈ RI , their ratios do not so that (42) has been defined consistently. It follows from (18) that (42) and (43) correspond to the choice of rate constants ˜ k) ˜ = (S, ˜ C, ˜ C˜K , R, ˜ defined by for the improperly translated mass action system M Definition 11. Consequently, from (37) we have that Y Ia Ik Ψ (x) = 0

⇐⇒

Y˜ Iã I˜k˜ Ψ˜ K (x) = 0

(44)

so that the steady states of the system (3) governing M and the steady states of the ˜ defined by Definition 11 coincide, and we are done. system (6) governing M

Appendix D: Proof of Theorem 5

˜ C, ˜ C˜K , R) ˜ be Proof Let N = (S, C, R) be a chemical reaction network and N˜ = (S, a weakly reversible translation of N which is either proper or improper and resolvable. Suppose M = (S, C, R, k) is a mass action system corresponding to N . We ˜ = (S, ˜ C, ˜ C˜K , R, ˜ k) ˜ according to Definidefine the translated mass action system M ˜ ˜ tion 8 if N is proper and by Definition 11 if N is resolvable and improper. ˜ correFrom either Lemma 2 and Lemma 4 we have that the steady state set of M sponds to the steady state set of M. Correspondingly, by either (11) or (44) we have that Y Ia Ik Ψ (x) = 0

⇐⇒

Y˜ A˜ k˜ Ψ˜ K (x) = 0,

˜ where A˜ k˜ := Iã I˜k˜ and Ψ˜ K (x) has entries [Ψ˜ K (x)]j = xyh3 (j ) for j ∈ CR. Since δ˜K = 0, we may conclude by Proposition 2.20 of (Müller and Regensburger ˜ has a complex balanced steady state. 2012) that the translated mass action system M m That is to say, there is a point a ∈ R>0 which satisfies Ψ˜ K (a) ∈ ker(A˜ k ).

(45)

1114

M.D. Johnston

Furthermore, since δ˜ = dim(ker(Y˜ ) ∩ Im(Iã )) = 0, we have from (38) that all steady states are complex balanced steady states. It follows from Proposition 2.21 of Müller and Regensburger (2012) that the set of such steady states may be parametrized by ˜⊥ E = x ∈ Rm >0 | ln(x) − ln(a) ∈ SK , where S˜K is given by (10) of Lemma 1. This is sufficient to prove claim 3. Now consider claims 1 and 2. Since N˜ is weakly reversible it follows by Theorem 6 that ˜ 1, K ˜ 2, . . . , K ˜ ˜}, ker(A˜ k ) = span{K ˜ j = ([K˜ j ]1 , [K˜ j ]2 , . . . , [K˜ j ]n˜ ) has entries where K K˜ i if i ∈ Lj , ˜ [Kj ]i = 0 otherwise

(46)

(47)

˜ are the tree constants corresponding to the reaction graph where K˜ i , i = 1, . . . , n, ˜ C, ˜ R, ˜ k). ˜ (S, ˜ It follows from (45), (46), and (47) that, for every i, j ∈ L˜ k for some k = 1, . . . , , satisfy the steady states x ∈ Rm >0 xyh3 (j ) xyh3 (i) = K˜ i K˜ j

⇐⇒ K˜ j xyh3 (i) − K˜ i xyh3 (j ) = 0.

Since this set corresponds to the steady states of M by either Lemma 2 or Lemma 4, we have shown that M has toric steady states generated by binomials of the form required by claim 2. Since the choice of rate constants in the definition of M was arbitrary, claim 1 follows and we are done. References Angeli, D., & Sontag, E. (2008). Translation-invariant monotone systems, and a global convergence result for enzymatic futile cycles. Nonlinear Anal., Real World Appl., 9, 128–140. Angeli, D., Leenheer, P., & Sontag, E. (2007). A Petri net approach to the study of persistence in chemical reaction networks. Math. Biosci., 210(2), 598–618. Clarke, B. L. (1980). Stability of complex reaction networks. Adv. Chem. Phys., 43, 1–215. Conradi, C., Flockerzi, D., & Raisch, J. (2008). Multistationarity in the activation of a MAPK: parametrizing the relevant region in parameter space. Math. Biosci., 211, 105–131. Cox, D., Little, J., & O’Shea, D. (2007). Undergraduate texts in mathematics. Ideals, varieties and algorithms (3rd ed.). Berlin: Springer. Craciun, G., & Feinberg, M. (2005). Multiple equilibria in complex chemical reaction networks: I. The injectivity property. SIAM J. Appl. Math., 65(5), 1526–1546. Craciun, G., & Feinberg, M. (2006). Multiple equilibria in complex chemical reaction networks: II. The species-reaction graph. SIAM J. Appl. Math., 66(4), 1321–1338. Craciun, G., Dickenstein, A., Shiu, A., & Sturmfels, B. (2009). Toric dynamical systems. J. Symb. Comput., 44(11), 1551–1565. Deng, J., Feinberg, M., Jones, C., & Nachman, A. (2011). On the steady states of weakly reversible chemical reaction networks. Preprint available on the arXiv:1111.2386. Dickenstein, A., & Pérez Millán, M. (2011). How far is complex balancing from detailed balancing? Bull. Math. Biol., 73, 811–828.


1115

Érdi, P., & Tóth, J. (1989). Mathematical models of chemical reactions. Princeton: Princeton University Press. Feinberg, M. (1972). Complex balancing in general kinetic systems. Arch. Ration. Mech. Anal., 49, 187– 194. Feinberg, M. (1979). Lectures on chemical reaction networks. Unpublished written versions of lectures given at the Mathematics Research Center, University of Wisconsin. Available at http://www. chbmeng.ohio-state.edu/~feinberg/LecturesOnReactionNetworks/. Feinberg, M. (1987). 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. Feinberg, M. (1988). Chemical reaction network structure and the stability of complex isothermal reactors: II. Multiple steady states for networks of deficiency one. Chem. Eng. Sci., 43(1), 1–25. Feinberg, M. (1989). Necessary and sufficient conditions for detailed balancing in mass action systems of arbitrary complexity. Chem. Eng. Sci., 44(9), 1819–1827. Feinberg, M. (1995a). The existence and uniqueness of steady states for a class of chemical reaction networks. Arch. Ration. Mech. Anal., 132, 311–370. Feinberg, M. (1995b). Multiple steady states for chemical reaction networks of deficiency one. Arch. Ration. Mech. Anal., 132, 371–406. Flockerzi, D., & Conradi, C. (2008). Subnetwork analysis for multistationarity in mass-action kinetics. J. Phys. Conf. Ser., 138(1). Gatermann, K. (2001). Counting stable solutions of sparse polynomial systems in chemistry. In E. L. Green, S. Hosten, R. C. Laubenbacher, & V. A. Powers (Eds.), Contemporary math: Vol. 286. Symbolic computation: solving equations in algebra, geometry and engineering (pp. 53–69). Gatermann, K., & Huber, B. (2002). A family of sparse polynomial systems arising in chemical reaction systems. J. Symb. Comput., 33(3), 275–305. Gatermann, K., & Wolfrum, M. (2005). Bernstein’s second theorem and Viro’s method for sparse polynomial systems in chemistry. Adv. Appl. Math., 34(2), 252–294. Gunawardena, J. (2005). Multisite protein phosphorylation makes a good threshold but can be a poor switch. Proc. Natl. Acad. Sci. USA, 102, 14617–14622. Gunawardena, J. (2007). Distributivity and processivity in multisite phosphorylation can be distinguished through steady-state invariants. Biophys. J., 93, 3828–3834. Holstein, K., Flockerzi, D., & Conradi, C. (2013). Multistationarity in sequentially distributed multisite phosphorylation networks. Bull. Math. Biol., 75, 2028–2058. Horn, F. (1972). Necessary and sufficient conditions for complex balancing in chemical kinetics. Arch. Ration. Mech. Anal., 49, 172–186. Horn, F., & Jackson, R. (1972). General mass action kinetics. Arch. Ration. Mech. Anal., 47, 187–194. Johnston, M. D. (2011). Topics in chemical reaction network theory. PhD thesis, University of Waterloo. Manrai, A., & Gunawardena, J. (2009). The geometry of multisite phosphorylation. Biophys. J., 95, 5533– 5543. Markevich, N. I., Hoek, J. B., & Kholodenko, B. N. (2004). Signaling switches and bistability arising from multisite phosphorylation in protein kinase cascades. J. Cell Biol., 164(3), 353–359. Müller, S., & Regensburger, G. (2012). Generalized mass action systems: complex balancing equilibria and sign vectors of the stoichiometric and kinetic-order subspaces. SIAM J. Appl. Math., 72(6), 1926– 1947. Müller, S., Feliu, E., Regensburger, G., Conradi, C., Shiu, A., & Dickenstein, A. (2013). Sign conditions for injectivity of generalized polynomial maps with applications to chemical reaction networks and real algebraic geometry. Preprint available on the arXiv:1311.5492. Pérez Millán, M., Dickenstein, A., Shiu, A., & Conradi, C. (2012). Chemical reaction systems with toric steady states. Bull. Math. Biol., 74(5), 1027–1065. Savageau, M. A. (1969). Biochemical systems analysis II. The steady-state solutions for an n-pool system using a power-law approximation. J. Theor. Biol., 25, 370–379. Shinar, G., & Feinberg, M. (2010). Structural sources of robustness in biochemical reaction networks. Science, 327(5971), 1389–1391. Shinar, G., & Feinberg, M. (2012). Concordant chemical reaction networks. Math. Biosci., 240(2), 92– 113. Shiu, A. J. (2010). Algebraic methods for biochemical reaction network theory. PhD thesis, University of California, Berkeley. Stanley, R. (1999). Enumerative combinatorics (Vol. 2). Cambridge: Cambridge University Press. Vol’pert, A. I., & Hudjaev, S. I. (1985). Analysis in classes of discontinuous functions and equations of mathematical physics. Dordrecht: Martinus Nijhoff.

1116

M.D. Johnston

Wang, L., & Sontag, E. (2008). On the number of steady states in a multiple futile cycle. J. Math. Biol., 57(1), 25–52. Wilhelm, T., & Heinrich, R. (1995). Smallest chemical reaction system with Hopf bifurcations. J. Math. Chem., 17(1), 1–14. Wilhelm, T., & Heinrich, R. (1996). Mathematical analysis of the smallest chemical reaction system with Hopf bifurcation. J. Math. Chem., 19(2), 111–130.

Stochastic flux analysis of chemical reaction networks.

Grip on complexity in chemical reaction networks.

Chemical reaction networks: Colour by number.

Complex Chemical Reaction Networks from Heuristics-Aided Quantum Chemistry.

Solution of Chemical Master Equations for Nonlinear Stochastic Reaction Networks.

Simple methods for calculating reaction equilibria in single reaction or single moiety chemical networks.

Plasmid-derived DNA Strand Displacement Gates for Implementing Chemical Reaction Networks.

Sensitivity of chemical reaction networks: a structural approach. 1. Examples and the carbon metabolic network.

Exact probability distributions of selected species in stochastic chemical reaction networks.

A technique for determining the signs of sensitivities of steady states in chemical reaction networks.

Thermodynamics of random reaction networks.

Efficiently gap-filling reaction networks.

Single-molecule chemical reaction reveals molecular reaction kinetics and dynamics.

Correction: Thermodynamics of random reaction networks.

Absolute asymmetric synthesis in enantioselective autocatalytic reaction networks: theoretical games, speculations on chemical evolution and perhaps a synthetic option.

Perspective: Bimolecular chemical reaction dynamics in liquids.

Chemical reaction-induced multi-molecular polarization (CRIMP).

Reaction efficiency effects on binary chemical reactions.

Nature-Inspired Chemical Reaction Optimisation Algorithms.

On correlated reaction sets and coupled reaction sets in metabolic networks.

Quantum-chemical insights from deep tensor neural networks.

Cluster synchronization in networks of neurons with chemical synapses.

Molecular Engineering of Robustness and Resilience in Enzymatic Reaction Networks.

Inference for reaction networks using the linear noise approximation.