Abstract
Efforts to catalogue the structure of metabolic networks have generated highly detailed, genomescale atlases of biochemical reactions in the cell. Unfortunately, these atlases fall short of capturing the kinetic details of metabolic reactions, instead offering only topological information from which to make predictions. As a result, studies frequently consider the extent to which the topological structure of a metabolic network determines its dynamic behaviour, irrespective of kinetic details. Here, we study a class of metabolic networks known as nonautocatalytic metabolic cycles, and analytically prove an open conjecture regarding the stability of their steady states. Importantly, our results are invariant to the choice of kinetic parameters, rate laws, equilibrium fluxes and metabolite concentrations. Unexpectedly, our proof exposes an elementary but apparently open problem of locating the roots of a sum of two polynomials S = P + Q, when the roots of the summand polynomials P and Q are known. We derive two new results named the Stubborn Roots Theorems, which provide sufficient conditions under which the roots of S remain qualitatively identical to the roots of P. Our study illustrates how complementary feedback, from classical fields such as dynamical systems to biology and vice versa, can expose fundamental and potentially overlooked questions.
1. Introduction
Networks of enzymecatalysed metabolic reactions are fundamental to the proliferation of life, using the energy extracted from environmental nutrients to drive the assembly of organic macromolecules and enable the successful reproduction of the cell. The largescale architecture of these networks is rich in structure: they are broadly organized into overlapping pathways (e.g. catabolic glycolysis and anabolic gluconeogenesis), and exhibit powerlawlike degree distributions with highly connected cofactors (e.g. ATP/ADP and NADH/NAD) linking many otherwisedistant metabolites [1]. Perhaps most importantly, these networks are capable of robust operation in spite of heterogeneity in the abundance of crucial enzymes and substrates [2].
To what extent are the robust features of metabolic networks determined by the underlying topological structure of the network itself? This question lies at the centre of many studies precisely because contemporary metabolic models are largely limited to structural information. Using genomic data, it is now possible to reconstruct genomescale models of metabolism, which predict the presence of absence or enzymes (and, by virtue, metabolic reactions) in an organism. However, the kinetic details of these reactions (such as the rate laws they obey as well as rate constants such as K_{m} or V_{max}) are largely unknown owing to the difficulty of measuring them in a highthroughput manner in vivo. As a result, a number of generic results linking the qualitative dynamics of a chemical reaction network (CRN; such as its potential for multistability or sustained oscillations) to its structural organization [3,4] have begun to populate the literature, pointing to a fundamental connection between structure and dynamics. These results make remarkably minimal and physically reasonable assumptions on the generic form of kinetic rate laws that render their conclusions largely independent of the choice of parameters. Perhaps the most wellknown result from studies of this sort is the deficiency zero theorem [3] from CRN theory (CRNT), which gives sufficient conditions for unique equilibria and asymptotic stability of a large class of reaction networks. More recent work extending CRNT, such as that of Shinar & Feinberg [2], has identified structural properties endowing CRNs with absolute concentration robustness (i.e. the steadystate concentration of a molecular species is identical in any steady state the dynamical system admits). Importantly, other existing methodologies, including chemical organization theory [5,6] and the theory of monotone systems [7], take similar ‘topological’ approaches to understanding how dynamics may be inferred from, and in fact directly influenced by, the structure of reaction networks themselves.
In previous work, Reznik & Segré [8] studied a family of nonautocatalytic metabolic cycles, and considered the role that their cyclic topology might play in determining their steadystate properties. Our decision to study a cycle, as opposed to any other structure, was motivated by the prominent role of cycles (such as the tricarboxylic acid and Calvin cycles) in presentday metabolic networks [9]. The approach we took relied on a method known as structural kinetic modelling (SKM), which applied a change of variables to the dynamical system corresponding to the metabolic cycle. This change in variables enabled us to study the dynamics of the metabolic cycle from a structural point of view, with only mild assumptions on the form of the kinetics themselves.
The main outcome of our work in [8] was limited numerical evidence that any steady state of the cycle must be stable to small perturbations, irrespective of equilibrium metabolite concentrations, flux magnitudes or the choice of kinetic parameters. However, we were unable to offer a rigorous analytical proof of this claim. In particular, it was unclear whether small regions of parameter space harbouring unstable equilibria might exist. Perhaps more importantly, computational considerations limited our numerical investigations to relatively short metabolic cycles (including up to eight metabolites), leaving open the possibility that instability appeared as cycles grew longer. As a result, we left the question of stability as an unproven conjecture (herein referred to as the cycle stability conjecture). This conjecture is the object of study in the first part of this article.
The difficulty with proving the cycle stability conjecture reduced to locating the roots of a highorder polynomial in the complex plane. Although a number of classical results from control theory are commonly applied to problems like this, they were rendered largely unusable for the polynomial in [8]. For example, the Routh–Hurwitz criterion, perhaps the best known technique for constraining the locations of a polynomial's roots, requires the precise calculation of coefficients of the polynomial under study. Unfortunately, the calculation of these coefficients for the polynomial in [8] became analytically intractable as the number of reactions in the cycle (which we assumed to be arbitrarily large) grew. Furthermore, because these coefficients were themselves functions of SKM variables, and were only constrained to lie in complicated intervals on the real axis, the problem of proving stability became substantially more difficult. Wellknown methods, such as Kharitonov's theorems [10], exist to study how uncertainty in the coefficients of a polynomial impacts its roots. In addition, the method of resultants, used by Gross & Feudel [11] in the context of studying generalized models, can be used to identify when pairs of imaginary roots cross the imaginary axis (suggesting the onset of instability and sustained oscillations). However, the complexity and scale of the polynomial in [8] again rendered the application of both Kharitonov's theorems and the method of resultants infeasible.
Our difficulty with bringing classical tools to bear on the cycle stability conjecture motivated us to revisit the problem from a completely new perspective. In this study, we resolve the cycle stability conjecture by reformulating it as a question of locating the roots of a sum of two polynomials. By doing so, we are able to apply a classical technique (Rouché's theorem [12]) from complex analysis to resolve the conjecture and prove the stability of nonautocatalytic metabolic cycles. Quite unexpectedly, our proof leads us to a substantially more general question: how do the roots of a sum of polynomials S = P + Q depend on the roots of P and Q themselves? Using a method identical to the one used to prove the cycle stability conjecture, we prove two new theorems, which we call the Stubborn Roots Theorems, which give sufficient conditions for when the roots of S are qualitatively identical to the roots of P. To the best of our knowledge, there are few generic results that provide information regarding the locations of the roots of a sum of two polynomials. Given the fundamental importance of locating the roots of a polynomial in the study of dynamical systems (and in applications of dynamical systems to fields such as in systems biology), we assume that the Stubborn Roots Theorems may find use in other contexts.
2. Results
2.1. Stability of metabolic cycles
We begin by presenting a model of the dynamics of a simple nonautocatalytic metabolic cycle, and then proceed to using SKM to study the stability of its equilibria. First, we describe the generic structure of the metabolic cycle under study, which is identical to the one studied in [8]. The cycle contains n metabolites (M_{1} … M_{n}), and two cofactors (O_{1} and O_{2}) that provide the energetic force to thermodynamically drive the metabolic reactions, and can be illustrated by 2.1
In this cycle, each metabolite M_{i} is converted to metabolite M_{i}_{+}_{1} for i = 1 … n − 1. A constant flux of M_{1} enters the system. A proportion of the last metabolite M_{n} is converted back to M_{1}, whereas the remainder leaves the system. The highenergy cofactor O_{1} is converted to its lowenergy cofactor partner O_{2} in the reaction catalysing the conversion of M_{1} to M_{2}. In a separate reaction, energy is input into the system to drive the reformation of the higher energy molecule O_{1}.
At steady state, the magnitude of the flux through each reaction in the cycle can be calculated by enforcing mass balances on each metabolite. To do so, we assume that a constant flux of generic magnitude αv, 0 < α < 1 of metabolite M_{1} flows into the network. A proportion (1 − α)v of the flux entering M_{n} is channelled back towards M_{1}, whereas the remaining flux αv exits the system. All other reactions carry a steadystate flux of v. It is easily verified that this flux vector is in the nullspace of the stoichiometric matrix S (see appendix A). We assume that the kinetics of each reaction is monotonic, that is, that an increase in the concentration of substrate for any reaction will consequently increase the rate of the reaction. This assumption is quite generic, and is amenable to many wellknown biochemical reaction mechanisms, including the law of mass action as well as Michaelis–Menten and Hill kinetics.
To prove the stability of a steady state of (2.1), we must prove that the Jacobian of (2.1), evaluated at an arbitrary steady state, always has eigenvalues with negative real part. As shown in [8] (and rederived in appendix A, equations (A 8)–(A 10)), the Jacobian J_{n} for the metabolic cycle of size n illustrated above can be calculated using SKM to be 2.2
Crucially, the assumption of monotonic kinetics constrains all the elasticities, θ_{i}, in (2.2) to be greater than zero (see appendix A). In [8], we provided evidence that J_{n} could not have eigenvalues with positive real part. Below, we proceed to analytically prove this conjecture. To simplify some calculations, we elect to work with the negative counterpart of the Jacobian, , and prove that cannot have eigenvalues with negative real part. This is equivalent to proving that J_{n} cannot have eigenvalues with positive real part.
First, we calculate the characteristic polynomial of which we call χ_{n}(λ), explicitly (calculations shown in appendix A) 2.3 2.4 2.5Thus, χ_{n} is the sum of two polynomials, an n + 1th order polynomial P_{n} and a firstorder polynomial Q_{1} (note that the subscript n denotes the size of the cycle, not the degree of the polynomial). Next, we prove three lemmas on the relative location of the roots of P_{n} and Q_{1}, showing that they are strongly constrained. Later on, these constraints will be crucial to proving that χ_{n} can only have roots with positive real part.
Lemma 2.1.
All of the roots of P_{n} are positive and real.
Proof. By inspection, at least n – 1 of P_{n}'s roots, contained in the product term of (2.3), must be positive and real. For the remaining two roots, we must study the quadratic polynomial 2.6
First, we prove that the roots of (2.6) must be real. Calculating the discriminant Δ of this quadratic polynomial, we find
Because Δ > 0, (2.6) cannot have imaginary roots and all of the roots of P_{n} are purely real.
Next, we show that the pair of roots of (2.6) must be positive. If we expand (2.6), we find 2.7
The product of the two roots of (2.7) is θ_{1}θ_{n}_{+}_{3} > 0, and the sum of the roots is θ_{1} + θ_{n}_{+}_{2} + θ_{n}_{+}_{3} > 0. Therefore, both of the roots of (2.7) must be positive.▪
Next, we prove a related lemma regarding the location of the root of Q_{1}.
Lemma 2.2.
The root r_{q} of Q_{1} must be larger than at least one root of P_{n}, and smaller than another root of P_{n}.
Proof. Consider the quadratic factor of P_{n}, (θ_{1} − λ)(θ_{n}_{+}_{2} + θ_{n}_{+}_{3} − λ) − θ_{1}θ_{n}_{+}_{2}. By lemma 2.1, this quadratic polynomial has distinct (because Δ > 0) real, positive roots. Let p_{1}, p_{2} denote these roots, ordered by magnitude so that p_{1} < p_{2}. Because the leading term of the quadratic is positive, the roots of the polynomial divide the real line into three regions: λ ≤ p_{1} and λ ≥ p_{2}, where the polynomial is greater than or equal to zero, and p_{1} < λ < p_{2}, where the polynomial is strictly negative. By inspection, r_{q} = θ_{n}_{+}_{3}. Directly evaluating the value of the quadratic polynomial at λ = r_{q} = θ_{n}_{+}_{3}, we find 2.8
Because the value of the quadratic factor is negative at λ = r_{q}, r_{q} must lie between the roots of (2.7).▪
Finally, we prove a lemma regarding the magnitude of P_{n} and Q_{1} at the origin.
Lemma 2.3.
P_{n}(0) > Q_{1}(0)
Proof. Explicitly calculating P_{n}(0), we find P_{n}(0) = (θ_{n} + θ_{n}_{+}_{1})θ_{n}_{+}_{3}Π_{i} _{=} _{1 … n}_{−}_{1}θ_{i}. This is always greater than Q_{1}(0) = θ_{n}θ_{n}_{+}_{3}Π_{i} _{=} _{1 … n}_{−}_{1}θ_{i}.▪
Now, using lemmas 2.1–2.3, we proceed to prove that the roots of χ_{n} must lie in the positive real half of the complex plane. To do so, we will make use of a wellknown theorem from complex analysis known as Rouché's theorem.
Theorem 2.1 (symmetric Rouché's theorem)
Two holomorphic functions f and g have the same number of roots within a region bounded by some continuous closed contour C (on which neither f nor g have any poles or zeros) if the strict inequality holds on C [12].
In essence, Rouché's theorem offers a way to determine the number of roots of a difference of two functions lying inside a closed contour in the complex plane. We will use f = S = P_{n} + Q_{1} and g = P_{n}. By taking contours bounding larger and larger regions of the lefthalfplane (those complex numbers with negative real part), we will show that Q_{1} < P_{n} on the contour, and thus prove that S has no roots inside this contour, i.e. no roots in the lefthalfplane.
We let our contour C_{R} consist of two parts (figure 1):

— the portion of the circle λ = R centred at the origin in the negative real half of the complex plane and

— the portion of the imaginary axis connecting the two points of intersection of the above circle with the imaginary axis.
First, we will prove that, along an arc of sufficiently large radius, P_{n}/Q_{1} > 1. Because n > 0, given an arc of sufficiently large radius R, Q_{1} < P_{n} simply because P_{n} is of higher order than Q_{1} (i.e. the highestorder term in P_{n} is λ^{n}^{+}^{1} which dominates the highestorder term in Q_{1}, λ^{1}, for very large λ).
Next, we will prove that, along the upper half of the imaginary axis, P_{n}/Q_{1} > 1. To do so, let us consider the behaviour of P_{n}/Q_{1} by substituting λ = iy, y > 0 into (2.5) and taking the modulus. Denoting the roots of P_{n} as r_{i} for i = 1, … , n + 1 and the root of Q_{1} as r_{q}, we have 2.9
and 2.10where c = Q_{1}(0) and r_{q} = c/b. We must show that the following condition holds for all y ≥ 0: 2.11Note that, at y = 0, we know (2.11) is satisfied because P_{n}(0) > Q_{1}(0). If we can show that x(y) is a strictly increasing function, then we know that (2.11) will be satisfied for all y. To do so, let us work with x^{2}(y). Note that x(y) > 0, x^{2}(0) > 1, and if d(x^{2})/dy > 0 for all y, then x^{2}(y) > 1 for all y. This would then imply that x(y) > 1 for all y ≥ 0. We have 2.12
Taking a derivative of x^{2} with respect to y and using the identity where f_{1} = P_{n}^{2}, f_{2} = 1/Q_{1}^{2} (the product rule applied to P_{n}^{2} and 1/Q_{1}^{2}), we find 2.13Without loss of generality, suppose r_{1} is the smallest root of P_{n}. From the first term on the righthand side of (2.13), select the term corresponding to i = 1. Recall that, by lemma 2.2, r_{1} must be smaller than r_{q}. Isolating just this term and the negative term in the parentheses on the righthand side of (2.13) and summing, we find 2.14
Because r_{q} > r_{1}, (2.14) is positive. There are no more negative terms in (2.13), proving that x^{2}(y) is strictly increasing. This proves that P_{n} > Q_{1} on the positive imaginary axis. Furthermore, because real polynomials are symmetric across the real axis, an identical argument shows that P_{n} > Q_{1} on the negative imaginary axis. In particular, setting λ = −iy for y > 0 yields the exact same expressions for P_{n} and Q_{1}.
We have satisfied all of the assumptions of Rouché's theorem, and have proved χ_{n}(λ) contains no roots in the left half of the complex plane. Therefore, the nonautocatalytic metabolic cycle always has stable equilibria.
2.2. The Stubborn Roots Theorem
Can we use the methods illustrated in §2.1 to locate the roots of the sum of two more general polynomials? Our motivation for studying this problem derives from control theory, where it is common to ask whether the roots of a polynomial lie in one half of the complex plane [13]. Such polynomials frequently correspond to the characteristic equation of the Jacobian matrix of a dynamical system. Although we do not provide a generic method for predicting whether a matrix's characteristic equation may be written as the sum of two simpler polynomials, the appearance of such structure in our studies of a metabolic cycle suggests that related ‘wellordered’ systems may exhibit similar properties.
The main question we ask in this section is under what conditions may the roots of a polynomial P be ‘stubborn’ when P is summed with another polynomial Q: in such a case, the roots of the summed polynomial S = P + Q remain qualitatively identical to those of P. By qualitatively identical, we mean specifically that the number of roots of P in the left (right) half of the complex plane is equal to the number of roots of S in the left (right) half of the complex plane. This question follows in the spirit of similar work by Anderson [14]. Our primary result is a theorem, which we call the Stubborn Real Roots Theorem, which gives sufficient conditions under which the location of the roots of a sum of polynomials S = P + Q remains qualitatively unchanged from P.
Theorem 2.2 (Stubborn Real Roots Theorem)
Let P_{n} and Q_{m} be polynomials of order n and m, respectively, and let n > m. Assume that all the roots of P_{n} and Q_{m} are purely real, and that P_{n}(0) > Q_{m}(0). Denote by p_{i} and q_{i} the roots of P_{n} and Q_{m} ordered by magnitude, so that p_{1} < p_{2} < … < p_{n} and If for every j = 1 … m, p_{j} < q_{j}, then the number of roots of S = P_{n} + Q_{m} located in the negative (positive) real half of the complex plane is equal to the number of roots of P_{n} in the negative (positive) real half of the complex plane.
The proof of theorem 2.2 is provided in appendix A, and follows precisely the same line of reasoning as the proof of the cycle stability conjecture in §2.1. The theorem relies on two critical assumptions relating P_{n} and Q_{m}. First, at the origin, P_{n}(0) > Q_{m}(0). Second, it must be possible to assign to each root of q_{i} of Q_{m} a unique root p_{i} of P_{n} such that q_{i} > p_{i}.
The power of theorem 2.2 is that it enables one to qualitatively locate the roots of a polynomial S simply by inspecting the roots of its summands P_{n} and Q_{m}. If the roots of P_{n} and Q_{m} are easily calculated (as in the case of the cycle stability conjecture in §2.1), then the roots of S can be immediately located without resorting to difficult calculations. In many ways, theorem 2.2 is reminiscent of the work of Fisk [15], who describes the behaviour of the roots of sums of polynomials which ‘interlace’. For two polynomials P and Q to interlace, the roots of P and Q alternate when ordered from most negative to most positive, so that p_{1} < q_{1} < p_{2} < q_{2} … . Notably, our result here is more general, and includes interlacing as a special case.
2.3. Stubborn Complex Roots
What happens when matters become complex? In this section, we generalize the Stubborn Roots Theorem to cases when the roots of P are not necessarily all real. This is often the case in dynamical systems, where complex roots indicate oscillatory phenomena such as spiralling or limit cycles [16]. Proceeding along the same lines as before, we find that the roots of P remain stubborn to the addition of Q as long as they remain predominantly real. That is, if the real component of the complex roots of P is larger than their imaginary component, then a more general version of the Stubborn Roots Theorem holds.
Theorem 2.3 (Stubborn Complex Roots Theorem)
Let P_{n} and Q_{m} be polynomials of order n and m, respectively, n > m. Let the m roots of Q_{m} be positive and purely real. Further, let P_{n} have at least m real roots, and let the remainder of the roots be either real or complex. Assume that P_{n}(0) > Q_{m}(0). Furthermore, assume that for each complex root of P, p_{k}, the magnitude of the real component Re(p_{k}) is larger than the magnitude of its imaginary component Im(p_{k}). Denote by p_{i} and q_{i} the real roots of P_{n} and Q_{m} ordered by magnitude, so that p_{1} < p_{2} < … < p_{m} and q_{1} < q_{2} < … < q_{m}. If for every j = 1 … m, p_{j} < q_{j}, then the number of roots of S = P_{n} + Q_{m} located in the negative (positive) real half of the complex plane is equal to the number of roots of P_{n} in the negative (positive) real half of the complex plane (figure 2).
Proof. We prove the theorem for the case when P_{n} has n – 2 real roots, two complex conjugate roots and and Q has one positive real root. The result can be straightforwardly (via wrenching and tedious algebraic calculations) extended to the generic case in theorem 2.3 using an identical argument. As before, we apply Rouché's theorem using a half circle in the negative real half of the complex plane using f = S = P_{n} + Q_{1} and g = P_{n}. First, we consider the behaviour of the two polynomials on the large arc in the negative real half of the complex plane. As before, P_{n} dominates Q_{1} as the radius of the arc grows larger.
Turning our attention to the behaviour of the polynomials on the positive imaginary axis, we substitute z = iy to find 2.15where for simplicity of notation, we have assumed P_{n} has leading coefficient 1. Differentiating and simplifying algebraically, we obtain the expression 2.16where C = ((y − b)^{2} + a^{2})((y + b)^{2} + a^{2}). Using an argument identical to the one used to prove the cycle stability conjecture, it is clear that the top term in (2.16) is positive. We are then left to ensure that the bottom term is positive. Expanding this term, we find 2.17Thus, (2.17) is certain to be positive as long as a > b. In this case, we can once more apply Rouché's theorem to show that S = P_{n} + Q_{1} has the same number of roots in each half of the complex plane as P_{n}.
3. Discussion
A major challenge in systems biology is efficiently studying biological networks through detailed atlases of their topological structures. These atlases, assembled from the cumulative results of many highthroughput, largescale experiments, capture many of the physical links that underlie such networks. However, they fail to describe most of the detailed dynamics taking place on the network itself. Here, we have studied the stability properties of a generic type of metabolic network, and analytically proved that, under quite mild assumptions on the reaction kinetics, any steady state of the network must be stable. To prove this, we reformulated the question of stability as a problem of locating the roots of a sum of two polynomials whose roots were easily calculated. This reformulation exposed a more fundamental problem of locating the roots of the sum of two polynomials, and we proved two new results (the Stubborn Roots Theorems) offering sufficient conditions under which the roots of the sum are not qualitatively different from the roots of one of the summands.
The study of metabolic networks and their stability has played an important role in research into the origin of life. Many of the earliest papers studying simple models of primordial metabolic networks focused on elucidating their stability properties, hypothesizing that molecular selforganization may have arisen through selfsustaining metabolic cycles [17]. Because early metabolism almost certainly lacked the complex regulatory mechanisms and circuits that appear in cells today, these investigations suggested that stability of these cycles to fluctuations in environmental conditions was necessary for their survival. Pursuing this line of thought, Piedrafita et al. [18] recently proposed a simple CRN composed of interlocking cycles that could establish and maintain a stable steady state, even in the face of a sudden loss of some constituent metabolite of the cycle itself. In response to such a catastrophe, the remaining metabolites reproduced the missing metabolite. This notion of closure, in which all of the metabolites necessary for sustaining the system can be produced by the metabolic network itself, also plays an important role in chemical organization theory, a method distinct from SKM for studying CRNs based on structure and mentioned earlier [5,6].
Naturally, one may ask: how important is the stability of equilibria to the robust function of presentday biological networks? A rich and diverse literature, dating back nearly half a century and still expanding today, describes the importance of stabilizing structures in ecological networks [19,20]. However, the extent to which stability plays a role in the fitness of metabolic systems is unclear, and studies investigating this question have failed to produce a definitive conclusion. In particular, if stability endowed metabolic networks with some evolutionary advantage, one should expect an enrichment for stabilizing features and structures in contemporary metabolic networks. While some studies have demonstrated the importance of some key stabilizing edges in metabolic networks, such as the allosteric feedback of ATP onto phosphofructokinase in glycolysis [21,22], others have failed to identify an enrichment of stabilizing structures in metabolism as a whole [23]. In fact, synthetic biologists routinely exploit instability in order to generate circuits exhibiting sustained oscillations, both in metabolic [24] and in transcriptional [25] systems. Importantly, we note that the work presented here only studied dynamics near equilibrium points, and ignored nonlocal dynamics (such as the appearance of periodic orbits arising from global bifurcations). Efforts extending generalized modelling and SKM to understanding nonlocal dynamics are now appearing in the literature [26].
We expect that the results presented here may find useful application in several challenges facing contemporary biology. First, SKM and generalized modelling could be used as coarsegrained techniques for vetting synthetic circuit designs for their potential to exhibit desirable behaviours (such as robust stability). In a previous study [27], we did precisely this, using generalized modelling to identify which topological circuit designs were entirely incapable of oscillations, irrespective of the choice of kinetic parameters or rate laws. Second, a great deal of interest now exists in using highthroughput metabolomics and fluxomics data to identify the role that smallmolecule (e.g. allosteric) regulation of metabolic enzymes plays in shaping metabolic dynamics [28]. Given the difficulty in accurately measuring kinetic rate constants in vivo, we envision that SKM (and our results here highlighting the inherent stability of certain topological motifs) might serve as a useful bridge between detailed mechanistic models and experimental data, highlighting those regulatory interactions that are crucial to the robust function of the network as a whole.
Finally, our results (in particular, the Stubborn Roots Theorems) illustrate the potential for biological questions to reveal interesting and unsolved problems in other fields. What appeared to us initially as a simple problem of locating the roots of the sum of two stable polynomials quickly blossomed into the exploration of widely diverse fields of active research, from control theory to matrix analysis. There is now a growing number of examples of similar feedback from biology to other fields, from classic results in evolutionary optimization [29] to the design of novel algorithms [30]. Interestingly, many of these crossfertilizations of ideas have taken place because of an abundance of biological data, and a need for analytical tools to understand them. Here, it has been quite the contrary: our study of a topological model of a metabolic cycle was motivated by a dearth of data on the kinetics of metabolic reactions. Nevertheless, in both cases, the ultimate outcome is deeper understanding, relevant to both biology and the fields from which it draws new tools and ideas.
Acknowledgements
We are grateful to Daniel Segrè and Gene Wayne for insightful feedback. E.R. was supported in part by grants from the Office of Science (BER), US Department of Energy (DESC0004962), National Science Foundation (NSF DMS0602204 EMSW21RTG Biodynamics at Boston University) and the US Department of Defense Army Research Office W911NF1210390 (UTA12001015). O.C. was supported in part by the National Science Foundation (NSF DMS0908093).
Appendix A
A.1. Structural kinetic modelling
To study the stability of a steady state of a metabolic network, we use a technique known as structural kinetic modelling (SKM) [22]. SKM is a nondimensionalization procedure that replaces conventional kinetic parameters (such as V_{max} and K_{m}) with normalized parameters known as elasticities. In the past, SKM (and its generalization, known as generalized modelling) has been paired with complementary methods studying other dynamic features of a system, such as the effects of noise [21].
As illustrated below, elasticities have several properties that make them powerful tools for studying metabolic dynamics. First, in contrast to kinetic parameters (whose values may be uncertain over many orders of magnitude), elasticities are constrained to lie in welldefined ranges (e.g. between zero and one), and sampling elasticities across this range effectively captures all possible values of kinetic parameters. Second, the value of an elasticity does not depend on the particular choice of kinetic rate law. Instead, an elasticity is simply a normalized measure of the sensitivity of a rate law to infinitesimal changes in a metabolite's concentration.
To study the stability of a steady state, SKM calculates the Jacobian matrix J of the dynamical system corresponding to a metabolic network. If we let C be the mdimensional vector of metabolite concentrations, N be the m × r stoichiometric network and v be the rdimensional vector of metabolic fluxes, then the dynamics of a metabolic network are governed by the system of differential equations A1where k is a vector of parameters and v(C, k) indicates that the vector of fluxes is dependent on both metabolite concentrations and kinetic parameters. Assuming that a nonzero steadystate C^{0} exists, we can make a change of variables and write A2where i = 1 … m and j = 1 … r.
Then, we can write the Jacobian as A3
The stability of the steadystate C^{0} is then dependent on the eigenvalues of J. If the real component of all eigenvalues of J is negative, then the steady state is stable. Thus, the problem of stability reduces to finding the eigenvalues of J.
The element that encodes the effective kinetic dependence of reaction rates on metabolites of SKM is the r × m elasticity matrix Θ. The (i, j)th element of Θ describes the sensitivity of the normalized rate of reaction i to the normalized concentration of metabolite j. This corresponds precisely to the effective kinetic order of the ith reaction with respect to the jth substrate: if the rate of reaction is linear with the amount of substrate, then θ = 1, while if it is zeroth order, θ = 0 [31]. Importantly, we assume that all elasticities in the metabolic cycle are greater than zero; that is, that an increase in the substrate of any reaction will increase the rate of that reaction. The analytical power of SKM comes precisely from the constrained and welldefined ranges of each element of Θ.
To illustrate the utility of elasticities, we derive below the elasticity of a metabolite involved in a Michaelis–Menten reaction. Consider a biochemical reaction governed by the rate law A4
Assuming that this reaction is embedded within a reaction network where metabolite S is at equilibrium concentration S_{0}, we can calculate the normalized reaction rate μ by normalizing (A 4) by its steadystate reaction rate, A5where x = S/S_{0} is the normalized concentration of S. Then, the elasticity is A6Note that, because S_{0} > 0, θ is constrained to the range (0, 1). The outcome of applying SKM to an entire metabolic network is a Jacobian, whose elements are formulated in terms of elasticities with welldefined ranges. Prior studies have used computational surveys [22,32] as well as analytical work [8] to study the role that particular key elasticities play in determining the stability of the network.
Finally, it may be useful to give a bit more intuition regarding the generality of an elasticity. To do so, we consider a reaction, governed by Michaelis–Menten kinetics, which exhibits an elasticity (explicitly calculated in equation (A 6)) of 0.5. First, note that this elasticity may correspond to any combination of S_{0} and K_{m} which satisfy K_{m}/(K_{m} + S_{0}) = 0.5; for example, K_{m} = S_{0} = 1 or K_{m} = S_{0} = 2. Thus, a single value for an elasticity in fact corresponds to a large locus of steadystate concentrations S_{0}. Now, note that, if we consider all elasticities in the range (0, 1), we in fact capture all possible combinations of S_{0} > 0 and K_{m} > 0! Thus, if we prove a theorem using this elasticity, and this theorem holds for all values of the elasticity in the range (0, 1), then it similarly holds for all possible choices of S_{0} and K_{m}. Obviously, the result also holds for any other kinetic rate laws for which this elasticity is valid. Thus, we effectively capture the entire space of possible parameters and steadystate concentrations. This is precisely the approach taken in proving the cycle stability conjecture.
A.2. Characteristic polynomial for the general metabolic cycle
Here, we formulate the structural kinetic model for the metabolic cycle depicted in (2.1) and reproduced below: A7
The analysis presented below is identical to that presented in [8]. The stoichiometric matrix S of the metabolic cycle is A8where the rows correspond to each metabolite in the system. The generalized forms of Λ, the normalized stoichiometric matrix, and Θ, the elasticity matrix, are shown below for the system depicted in (2.1). Note that we make no assumptions on the steadystate concentrations of metabolites, denoted by the vector (M_{1}, M_{2}, … , M_{n}, O_{1}, O_{2}). Furthermore, the steadystate flux through the cycle is equal to a generic magnitude v, v > 0, except for the input and outflow reactions (with flux α v) and the reaction from M_{n} to M_{1} (with flux (1 − α)v), A9and A10
With N + 1 metabolites and N + 3 reactions, Λ is N + 1 × (N + 3) and Θ is (N + 3) × N + 1. Note that the last row (corresponding to cofactor O_{2}) of Λ is omitted because the cofactors come as a conserved pair, and the bottom right element of (corresponding to the dependence of the last reaction on O_{2}) is replaced with a negative element in order to account for this conservation (for more information on modelling of conserved moieties, see the supplementary information text in Steuer et al. [22]). Then, the Jacobian We elect to study the eigenvalues of the negative of J_{n}, which we call Note that the eigenvalues of are precisely the negative of the eigenvalues of J_{n}. Thus, proving that all the eigenvalues of have positive real part is equivalent to proving that all of the eigenvalues of J_{n} have negative real part. The characteristic polynomial of , χ_{n}(λ) is
Above, we have made a change of variables so that , i = 2 … n + 2, , and . Note that this change of variables does not affect our assumption that all values of θ are greater than zero.
Now, we proceed to explicitly calculate χ_{n}(λ) by calculating the determinant. Expanding along the n – 1th column, we have
The important thing to note is that, now, the n – 1th column in the first matrix and the nth row in the second are zero except for a single entry. We can continue expanding along such columns until we arrive at
Finally, we have, for n ≥ 3, A11
If we make a final change of variables so that , i = 1 … n, then we arrive precisely at (2.5).
A.3. Characteristic polynomial for a simple metabolic cycle
To give the reader more intuition regarding the mechanics of SKM calculations, we briefly describe how SKM may be used to model the dynamics of a simple, nonautocatalytic metabolic cycle with two metabolites and one cofactor pair. The reactions of this cycle are A12
Then, the stoichiometric matrix S, the normalized stoichiometric matrix Λ and the elasticity matrix Θ may be written A13 A14 A15where v is an arbitrary unit of flux, , and A^{0}, B^{0} and are the steadystate concentrations of A, B and O_{1}, respectively.
Note that Λ has precisely the same structure as S, but with one fewer row (owing to the mass conservation associated with the cofactor pair, discussed in appendix A.2). Finally, the Jacobian for this system may be straightforwardly written as A16
A.4. Proof of the Real Stubborn Roots Theorem
Consider the sum of two polynomials A17 A18 A19
Above, all roots p_{i} and q_{i} are real and positive. Further, we assume that P_{n}(0) > Q_{1}(0). We claim that P cannot have any roots lying in the left half of the complex plane. To prove this, we make use of the symmetric form of Rouché's theorem (theorem 2.1 in the text). We will use f = S = P_{n} + Q_{m} and g = P_{n}. Our goal will be to show that Q_{m} < P_{n} on the contour, which by Rouché's theorem gives us that P_{n} and S have the same number of zeros inside the contour (which is none, because all of the roots of P_{n} are positive real numbers). By taking contours bounding larger and larger regions contained in and tending towards the entire lefthalfplane of , we will deduce that S has no roots with negative real component.
We let our contour C = C_{R} consist of two parts:
— a large half circle in the left half of the complex plane and
— the portion of the imaginary axis connecting the two intersections of the above circle with the positive imaginary axis.
First, we verify that Q_{m} < P_{n} on the half circle λ = R in the lefthalfplane. Because P_{n} is of higher order than Q_{m}, as R → ∞, P_{n} dominates Q_{m}.
Next, we verify that Q_{m} < P_{n} on the imaginary axis. Again substituting λ = iy for y > 0 into (A 19), we have A20and A21Then, A22and it follows that A23
Rewriting this, we find A24
Ignoring the common denominator in both terms and considering only the first m terms of the first series, we obtain A25
Now, comparing the terms of each series in order, we observe precisely the same pattern as in the main text; because p_{j} < q_{j} for all j = 1 … m, (A 25) is always greater than zero. Then, P_{n} > Q_{m} on the positive imaginary axis, and by symmetry also on the negative imaginary axis. Applying Rouché's theorem once again, we find that the number of roots of S = P_{n} + Q_{m} in the positive (negative) real half of the complex plane is identical to the number of roots of P_{n} in the positive (negative) real half of the complex plane.
 Received January 29, 2013.
 Accepted March 12, 2013.
 © 2013 The Author(s) Published by the Royal Society. All rights reserved.