The stubborn roots of metabolic cycles

Ed Reznik, Alex Watson, Osman Chaudhary


Efforts to catalogue the structure of metabolic networks have generated highly detailed, genome-scale 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 non-autocatalytic 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 enzyme-catalysed 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 large-scale architecture of these networks is rich in structure: they are broadly organized into overlapping pathways (e.g. catabolic glycolysis and anabolic gluconeogenesis), and exhibit power-law-like degree distributions with highly connected cofactors (e.g. ATP/ADP and NADH/NAD) linking many otherwise-distant 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 genome-scale 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 Km or Vmax) are largely unknown owing to the difficulty of measuring them in a high-throughput 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 multi-stability 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 well-known 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 steady-state 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 non-autocatalytic metabolic cycles, and considered the role that their cyclic topology might play in determining their steady-state 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 present-day 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 high-order 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. Well-known 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 non-autocatalytic 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 non-autocatalytic 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 (M1 … Mn), and two cofactors (O1 and O2) that provide the energetic force to thermodynamically drive the metabolic reactions, and can be illustrated byEmbedded Image 2.1

In this cycle, each metabolite Mi is converted to metabolite Mi+1 for i = 1 … n − 1. A constant flux of M1 enters the system. A proportion of the last metabolite Mn is converted back to M1, whereas the remainder leaves the system. The high-energy cofactor O1 is converted to its low-energy cofactor partner O2 in the reaction catalysing the conversion of M1 to M2. In a separate reaction, energy is input into the system to drive the reformation of the higher energy molecule O1.

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 M1 flows into the network. A proportion (1 − α)v of the flux entering Mn is channelled back towards M1, whereas the remaining flux αv exits the system. All other reactions carry a steady-state 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 well-known 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 re-derived in appendix A, equations (A 8)–(A 10)), the Jacobian Jn for the metabolic cycle of size n illustrated above can be calculated using SKM to beEmbedded Image 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 Jn 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, Embedded Image, and prove that Embedded Image cannot have eigenvalues with negative real part. This is equivalent to proving that Jn cannot have eigenvalues with positive real part.

First, we calculate the characteristic polynomial of Embedded Image which we call χn(λ), explicitly (calculations shown in appendix A)Embedded Image 2.3Embedded Image 2.4Embedded Image 2.5Thus, χn is the sum of two polynomials, an n + 1th order polynomial Pn and a first-order polynomial Q1 (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 Pn and Q1, 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 Pn are positive and real.

Proof. By inspection, at least n – 1 of Pn'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 polynomialEmbedded Image 2.6

First, we prove that the roots of (2.6) must be real. Calculating the discriminant Δ of this quadratic polynomial, we findEmbedded Image

Because Δ > 0, (2.6) cannot have imaginary roots and all of the roots of Pn are purely real.

Next, we show that the pair of roots of (2.6) must be positive. If we expand (2.6), we findEmbedded Image 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 Q1.

Lemma 2.2.

The root rq of Q1 must be larger than at least one root of Pn, and smaller than another root of Pn.

Proof. Consider the quadratic factor of Pn, (θ1λ)(θn+2 + θn+3λ) − θ1θn+2. By lemma 2.1, this quadratic polynomial has distinct (because Δ > 0) real, positive roots. Let p1, p2 denote these roots, ordered by magnitude so that p1 < p2. Because the leading term of the quadratic is positive, the roots of the polynomial divide the real line into three regions: λp1 and λp2, where the polynomial is greater than or equal to zero, and p1 < λ < p2, where the polynomial is strictly negative. By inspection, rq = θn+3. Directly evaluating the value of the quadratic polynomial at λ = rq = θn+3, we findEmbedded Image 2.8

Because the value of the quadratic factor is negative at λ = rq, rq must lie between the roots of (2.7).▪

Finally, we prove a lemma regarding the magnitude of Pn and Q1 at the origin.

Lemma 2.3.

|Pn(0)| > |Q1(0)|

Proof. Explicitly calculating Pn(0), we find |Pn(0)| = |(θn + θn+1)θn+3Πi = 1 … n1θi|. This is always greater than |Q1(0)| = |θnθn+3Πi = 1 … n1θ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 well-known 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 inequalityEmbedded Image 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 = Pn + Q1 and g = Pn. By taking contours bounding larger and larger regions of the left-half-plane (those complex numbers with negative real part), we will show that |Q1| < |Pn| on the contour, and thus prove that S has no roots inside this contour, i.e. no roots in the left-half-plane.

We let our contour CR 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.

Figure 1.

The semicirular contour CR (dashed line) and the region it encloses (grey region) used in the proof of the cycle stability conjecture. By allowing the radius R of the semicircle to go to infinity, we are able to encompass the entire left half of the complex plane.

First, we will prove that, along an arc of sufficiently large radius, |Pn|/|Q1| > 1. Because n > 0, given an arc of sufficiently large radius R, |Q1| < |Pn| simply because Pn is of higher order than Q1 (i.e. the highest-order term in Pn is λn+1 which dominates the highest-order term in Q1, λ1, for very large λ).

Next, we will prove that, along the upper half of the imaginary axis, |Pn|/|Q1| > 1. To do so, let us consider the behaviour of |Pn|/|Q1| by substituting λ = iy, y > 0 into (2.5) and taking the modulus. Denoting the roots of Pn as ri for i = 1, , n + 1 and the root of Q1 as rq, we haveEmbedded Image 2.9

andEmbedded Image 2.10where c = |Q1(0)| and rq = c/b. We must show that the following condition holds for all y ≥ 0:Embedded Image 2.11Note that, at y = 0, we know (2.11) is satisfied because |Pn(0)| > |Q1(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 x2(y). Note that x(y) > 0, x2(0) > 1, and if d(x2)/dy > 0 for all y, then x2(y) > 1 for all y. This would then imply that x(y) > 1 for all y ≥ 0. We haveEmbedded Image 2.12

Taking a derivative of x2 with respect to y and using the identityEmbedded Image where f1 = |Pn|2, f2 = 1/|Q1|2 (the product rule applied to |Pn|2 and 1/|Q1|2), we findEmbedded Image 2.13Without loss of generality, suppose r1 is the smallest root of Pn. From the first term on the right-hand side of (2.13), select the term corresponding to i = 1. Recall that, by lemma 2.2, r1 must be smaller than rq. Isolating just this term and the negative term in the parentheses on the right-hand side of (2.13) and summing, we findEmbedded Image 2.14

Because rq > r1, (2.14) is positive. There are no more negative terms in (2.13), proving that x2(y) is strictly increasing. This proves that |Pn| > |Q1| on the positive imaginary axis. Furthermore, because real polynomials are symmetric across the real axis, an identical argument shows that |Pn| > |Q1| on the negative imaginary axis. In particular, setting λ = −iy for y > 0 yields the exact same expressions for |Pn| and |Q1|.

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 non-autocatalytic 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 ‘well-ordered’ 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 Pn and Qm be polynomials of order n and m, respectively, and let n > m. Assume that all the roots of Pn and Qm are purely real, and that |Pn(0)| > |Qm(0)|. Denote by pi and qi the roots of Pn and Qm ordered by magnitude, so that |p1| < |p2| < … < |pn| and Embedded Image If for every j = 1 … m, |pj| < |qj|, then the number of roots of S = Pn + Qm located in the negative (positive) real half of the complex plane is equal to the number of roots of Pn 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 Pn and Qm. First, at the origin, |Pn(0)| > |Qm(0)|. Second, it must be possible to assign to each root of qi of Qm a unique root pi of Pn such that qi > pi.

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 Pn and Qm. If the roots of Pn and Qm 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 p1 < q1 < p2 < q2 … . 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 Pn and Qm be polynomials of order n and m, respectively, n > m. Let the m roots of Qm be positive and purely real. Further, let Pn have at least m real roots, and let the remainder of the roots be either real or complex. Assume that |Pn(0)| > |Qm(0)|. Furthermore, assume that for each complex root of P, pk, the magnitude of the real component |Re(pk)| is larger than the magnitude of its imaginary component |Im(pk)|. Denote by pi and qi the real roots of Pn and Qm ordered by magnitude, so that |p1| < |p2| < … < |pm| and |q1| < |q2| < … < |qm|. If for every j = 1 … m, |pj| < |qj|, then the number of roots of S = Pn + Qm located in the negative (positive) real half of the complex plane is equal to the number of roots of Pn in the negative (positive) real half of the complex plane (figure 2).

Figure 2.

Region of validity (in grey) for the Stubborn Complex Roots Theorem. To apply the Stubborn Complex Roots Theorem, the complex roots of Pn must have a smaller imaginary component than real component.

Proof. We prove the theorem for the case when Pn has n – 2 real roots, two complex conjugate roots Embedded Image and Embedded Image 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 = Pn + Q1 and g = Pn. First, we consider the behaviour of the two polynomials on the large arc in the negative real half of the complex plane. As before, |Pn| dominates |Q1| 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 findEmbedded Image 2.15where for simplicity of notation, we have assumed Pn has leading coefficient 1. Differentiating and simplifying algebraically, we obtain the expressionEmbedded Image 2.16where C = ((yb)2 + a2)((y + b)2 + a2). 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 findEmbedded Image 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 = Pn + Q1 has the same number of roots in each half of the complex plane as Pn.

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 high-throughput, large-scale 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 re-formulated 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 self-organization may have arisen through self-sustaining 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 re-produced 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 present-day 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 non-local dynamics (such as the appearance of periodic orbits arising from global bifurcations). Efforts extending generalized modelling and SKM to understanding non-local 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 coarse-grained 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 high-throughput metabolomics and fluxomics data to identify the role that small-molecule (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 cross-fertilizations 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.


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 (DE-SC0004962), National Science Foundation (NSF DMS-0602204 EMSW21-RTG Biodynamics at Boston University) and the US Department of Defense Army Research Office W911NF-12-1-0390 (UTA12-001015). O.C. was supported in part by the National Science Foundation (NSF DMS-0908093).

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 non-dimensionalization procedure that replaces conventional kinetic parameters (such as Vmax and Km) 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 well-defined 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 m-dimensional vector of metabolite concentrations, N be the m × r stoichiometric network and v be the r-dimensional vector of metabolic fluxes, then the dynamics of a metabolic network are governed by the system of differential equationsEmbedded Image 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 non-zero steady-state C0 exists, we can make a change of variables and writeEmbedded Image A2where i = 1 … m and j = 1 … r.

Then, we can write the Jacobian asEmbedded Image A3

The stability of the steady-state C0 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 well-defined 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 lawEmbedded Image A4

Assuming that this reaction is embedded within a reaction network where metabolite S is at equilibrium concentration S0, we can calculate the normalized reaction rate μ by normalizing (A 4) by its steady-state reaction rate,Embedded Image A5where x = S/S0 is the normalized concentration of S. Then, the elasticity isEmbedded Image A6Note that, because S0 > 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 well-defined 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 S0 and Km which satisfy Km/(Km + S0) = 0.5; for example, Km = S0 = 1 or Km = S0 = 2. Thus, a single value for an elasticity in fact corresponds to a large locus of steady-state concentrations S0. Now, note that, if we consider all elasticities in the range (0, 1), we in fact capture all possible combinations of S0 > 0 and Km > 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 S0 and Km. 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 steady-state 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:Embedded Image A7

The analysis presented below is identical to that presented in [8]. The stoichiometric matrix S of the metabolic cycle isEmbedded Image 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 steady-state concentrations of metabolites, denoted by the vector (M1, M2, , Mn, O1, O2). Furthermore, the steady-state 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 Mn to M1 (with flux (1 − α)v),Embedded Image A9andEmbedded Image 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 O2) of Λ is omitted because the cofactors come as a conserved pair, and the bottom right element of Embedded Image (corresponding to the dependence of the last reaction on O2) 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 Embedded Image We elect to study the eigenvalues of the negative of Jn, which we call Embedded Image Note that the eigenvalues of Embedded Image are precisely the negative of the eigenvalues of Jn. Thus, proving that all the eigenvalues of Embedded Image have positive real part is equivalent to proving that all of the eigenvalues of Jn have negative real part. The characteristic polynomial of Embedded Image, χn(λ) isEmbedded Image

Above, we have made a change of variables so that Embedded Image, i = 2 … n + 2, Embedded Image, and Embedded Image. 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

Embedded Image

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

Embedded Image

Finally, we have, for n ≥ 3,Embedded Image A11

If we make a final change of variables so that Embedded Image, i = 1 … n, Embedded Image Embedded Image Embedded Image 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, non-autocatalytic metabolic cycle with two metabolites and one cofactor pair. The reactions of this cycle areEmbedded Image A12

Then, the stoichiometric matrix S, the normalized stoichiometric matrix Λ and the elasticity matrix Θ may be writtenEmbedded Image A13Embedded Image A14Embedded Image A15where v is an arbitrary unit of flux, Embedded Image, and A0, B0 and Embedded Image are the steady-state concentrations of A, B and O1, 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 asEmbedded Image A16

A.4. Proof of the Real Stubborn Roots Theorem

Consider the sum of two polynomialsEmbedded Image A17Embedded Image A18Embedded Image A19

Above, all roots pi and qi are real and positive. Further, we assume that |Pn(0)| > |Q1(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 = Pn + Qm and g = Pn. Our goal will be to show that |Qm| < |Pn| on the contour, which by Rouché's theorem gives us that Pn and S have the same number of zeros inside the contour (which is none, because all of the roots of Pn are positive real numbers). By taking contours bounding larger and larger regions contained in and tending towards the entire left-half-plane of Embedded Image, we will deduce that S has no roots with negative real component.

We let our contour C = CR 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 |Qm| < |Pn| on the half circle |λ| = R in the left-half-plane. Because Pn is of higher order than Qm, as R, |Pn| dominates |Qm|.

Next, we verify that |Qm| < |Pn| on the imaginary axis. Again substituting λ = iy for y > 0 into (A 19), we haveEmbedded Image A20andEmbedded Image A21Then,Embedded Image A22and it follows thatEmbedded Image A23

Rewriting this, we findEmbedded Image A24

Ignoring the common denominator in both terms and considering only the first m terms of the first series, we obtainEmbedded Image A25

Now, comparing the terms of each series in order, we observe precisely the same pattern as in the main text; because pj < qj for all j = 1 … m, (A 25) is always greater than zero. Then, |Pn| > |Qm| 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 = Pn + Qm in the positive (negative) real half of the complex plane is identical to the number of roots of Pn in the positive (negative) real half of the complex plane.

  • Received January 29, 2013.
  • Accepted March 12, 2013.


View Abstract