{'fraction_non_alphanumeric': 0.025728701888995453, 'fraction_numerical': 0.005983903923939785, 'mean_word_length': 4.714643155010127, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Stability of Metabolic Cycles

We investigate the stability properties of two different classes of metabolic cycles using a combination of analytical and computational methods. Using principles from structural kinetic modeling (SKM), we show that the stability of metabolic networks with certain structural regularities can be studied using a combination of analytical and computational techniques. We then apply these techniques to a class of single input, single output metabolic cycles, and find that the cycles are stable under all conditions tested. Next, we extend our analysis to a small autocatalytic cycle, and determine parameter regimes within which the cycle is very likely to be stable. We demonstrate that analytical methods can be used to understand the relationship between kinetic parameters and stability, and that results from these analytical methods can be confirmed with computational experiments. In addition, our results suggest that elevated metabolite concentrations and certain crucial saturation parameters can strongly affect the stability of the entire metabolic cycle. We discuss our results in light of the possibility that evolutionary forces may select for metabolic network topologies with a high intrinsic probability of being stable. Furthermore, our conclusions support the hypothesis that certain types of metabolic cycles may have played a role in the development of primitive metabolism despite the absence of regulatory mechanisms.


Introduction
Cycles are at the heart of the metabolic networks of organisms spanning the entire tree of life [1,2,3]. For example, the tricarboxylic acid (TCA) cycle sits at the core of energy production for many species and additionally plays the role of regenerating essential cellular nutrients and components. A great deal of research has been devoted to understanding the stability properties of such cycles, which have often been represented as chemical reaction networks (CRNs). Many of these studies have drawn conclusions about the steady state of a CRN without regard to details about the rate constants or chemical concentrations of the particular CRN (for instance, see [4,5,6]). Instead, these approaches rely on the topology of the network and commonly assume mass-action kinetics laws to limit (and sometimes determine) the possible behaviors the system may exhibit at steady state. Can one hope to extend these results to general, potentially nonlinear, rate laws?
In this work, we use a generalized modeling framework known as structural kinetic modeling (SKM) to extend previous results on two classes of metabolic cycles. The first class has a special but relatively flexible structure: it contains a single input, a single output, and unlimited length. This relatively organized structure makes it amenable to analytical methods of study. For cycles belonging to this first class, we prove that the cycle can only lose dynamical stability in an oscillatory manner, regardless of absolute metabolite concentrations, flux values, kinetic rate constants, and form of monotonic kinetic rate law. We further investigate the cycle computationally, and find that it remains completely stable under all conditions tested, leading us to make a conjecture about its stability under all conditions. We then apply the same SKM approach to a second class of cycles, autocatalytic cycles, and demonstrate that these are not stable under all combinations of kinetic parameters and metabolic conditions. Instead, we use analytical calculations to derive relationships between key parameters which can enhance or attenuate the stability properties of the cycle. These relationships are then computationally tested and verified.
Nearly all prior published work using SKM [7,8] employs computational experiments and statistical methods in order to draw conclusions about the system of interest. In this work, we show that metabolic systems with relatively well-organized structures are not only analytically tractable but furthermore that their analysis can actually lead to interesting and novel conclusions. Finally, we discuss the implications of our results on the evolution of primitive metabolic networks. The SKM technique used throughout the paper highlights essential parameters which directly impact the stability of the cycle under study. Particular combinations of these parameters may have conferred an evolutionary advantage to their respective metabolic networks by making them robust to small environmental perturbations.

Structural Kinetic Modeling (SKM)
The methods developed in this paper are based on the SKM framework [9]. SKM is a specific application of generalized modeling [10] in which normalized parameters replace conventional parameters (such as V max or K M in the modeling of metabolic networks). The normalized parameters have a direct connection to the original kinetic parameters, but are much easier to work with. As will be shown below, these parameters usually have well-defined and limited ranges (e.g. [0,1]), and sampling them across this range effectively samples all possible values of the original kinetic parameters.
The goal of SKM is to capture the local stability properties of a biochemical system. In this sense, it bridges the gap between genome-scale steady state modeling [11] and explicit kinetic modeling of a metabolic process. To study stability, one usually determines the Jacobian J of the system of interest and evaluates it at the system's steady states. Assuming knowledge of the form of kinetic rate laws for a CRN, it is quite possible to write down the corresponding J. However, in many cases J will be unnecessarily complicated and quite difficult to work with. By performing a change of variables, SKM actually simplifies the mathematical form of each entry in J. The new entries are in almost all cases easier to work with. The specific method by which the change of variables occurs is briefly described below and more completely in Appendix A. Much of this material is based on [9] and we refer readers to that reference and its supplementary materials for more information.
If we let S be the m-dimensional vector of metabolite concentrations, N be the m ×r stoichiometric matrix, and v be the r -dimensional vector of reaction rates, then we can describe the dynamics of the system with the equation where v(S, k) denotes that the reaction rates are dependent on both metabolite concentrations S and kinetic parameters (such as Michaelis-Menten constants) k. If we assume that a non-negative steady state S 0 exists, then we can redefine our system using the definitions where i = 1...m and j = 1...r, x is a vector of metabolite concentrations normalized with respect to their steady state concentrations and µ represents flux normalized with respect to steady state flux values. The matrix Λ represents the stoichiometric matrix normalized with respect to steady state fluxes and steady state metabolite concentrations.
The Jacobian, evaluated at x 0 = 1 (which, because of the way x is defined, is actually the equilibrium of the system), can be written as Note that the equations above were derived without regard to the actual form of the kinetic equations that determined the ODE system. The matrix θ µ x contains elements which represent the degree of saturation of normalized flux µ j with respect to normalized substrate concentration x i . In terms of derivatives, each element of θ represents the degree of change in a flux as a particular metabolite is incrementally increased. This is analogous to the concept of elasticity in metabolic control analysis [12].
What does the θ matrix look like? Its columns correspond to each metabolite, and its rows to each flux. A non-zero element θ i j in the matrix represents the effect a small change in a metabolite j has on flux i. In the case of Michaelis-Menten kinetics, this element in the matrix may take values ranging from [0,1]. In the case of standard competitive inhibition (e.g. allosteric inhibition by a product), the element takes values in [-1,0].
To usefully illustrate the meaning of a single θ element, consider an equation following Michaelis-Menten kinetics shown in (2.4). First, we write µ(x), which we recall is the flux normalized by the flux at the steady state. Here, S 0 is the concentration of the substrate S at steady state. Then, we manipulate the equation by substituting xS 0 for S, where x is the normalized steady state concentration of substrate S. The result is shown in (2.5).
Finally, we take a derivative with respect to x and evaluate it at x = 1 to obtain θ in (2.6). Notice that θ can only take values between 0 and 1 for any positive value of S 0 .
The power of SKM comes from the parametrization illustrated above. Each element in θ has a precise correspondence to some combination of kinetic parameters in the original model. However, it is far more tractable to study a system using θ parameters rather than the original kinetic parameters. To see this more clearly, recall that in most cases, biochemical kinetic constants are poorly estimated. In order to build precise ODE models for biochemical systems, it is usually necessary to actually choose values for these constants. While the chosen values may be estimated from experimental measurements, hidden dependencies in the model may actually result in non-obvious correlations between parameters that can strongly affect the output of the model.
On the other hand, if one is only interested in the stability of a characterized steady-state, it may not be necessary to actually have knowledge of kinetic parameters. Experimental measurements can provide data on absolute metabolite concentrations and flux values, and the stoichiometry is often known a priori. Then, one can parameterize the system as shown above and sample many possible combinations of θ parameters. For each unique set of θ parameters, the stability of the Jacobian is determined and recorded. Analysis of the results can lead to several important conclusions such as which θ parameters have the strongest correlation with stability of the entire metabolic system. Thus, the benefit of employing SKM over other techniques is that it provides the means to analyze and make sense of a large number of possible cases of a metabolic network, rather than a single instance.
At first glance, it is unclear whether this will make things easier or harder; values for fluxes and concentrations are almost never known. However, our work to follow shows that the absence of this knowledge does not make the problem of stability intractable. Instead, we operate with the following idea: we assume that the steady state flux vector v and the metabolite concentration vector X are variables, and we do not specify them. Instead, we simply determine the stability of the system in terms of these and other (θ) variables. Then, we can find trends in stability as the flux and metabolite concentrations change: it may be that as a concentration goes up, the stability of the system tends to fall. Isolating these trends is the motivation of our analysis.

Previous results on stability of metabolic cycles
Our work is preceded by a great deal of work aimed at understanding the stability of metabolic networks [3,13,14,15,16]. One flavor of such work is referred to as chemical reaction network theory (CRNT). CRNT is an elegant and powerful framework for understanding how the structure of chemical networks dictates their equilibrium and stability properties. By applying CRNT to the cycles studied in this paper, one can (in some cases) prove the stability of their steady states. This result can be easily obtained, for example, by feeding the stoichiometric matrices of the corresponding networks into the ERNEST MATLAB Toolbox [17]. Notably, these results cannot be easily extended beyond the case of mass-action kinetics (where the rate of reaction is proportional to the concentration of its reactants). The SKM method presented in this article extends the stability results for these cycles to general monotonic kinetics.
Dynamical systems approaches have also been used in understanding the stability properties of autocatalytic cycles (the second type of cycle presented in the paper). Previous research by King [14] (focusing on dynamics) and by Eigen [3] (focusing on steady-state stability properties) have studied autocatalytic cycles assuming first-order (non-saturating) mass-action kinetic rate laws. Extending such work, our analysis takes into account all forms of saturating rate laws up to order one, and derives simple expressions for the stability of metabolic cycles given adequate information about rate constants and steady-state concentrations. To put this in context, the prior work on autocatalytic cycles represents a single point (at which all θ parameters equal 1) in the large space of θ parameters within which the actual metabolic network resides. Our approach searches over this entire space, and identifies trends in stability across it.

A Simple Example
To illustrate the main concepts behind analytically characterizing the stability of a metabolic cycle using SKM, we will start with a simple example. Consider the cycle shown in Figure 1. To analyze the conditions under which it is stable, we will write the Λ and θ matrices, multiply them together to obtain the Jacobian J, and then find the eigenvalues of J by writing out its characteristic equation.
The cycle contains 3 metabolites and 6 reactions. We impose that the reactions can only proceed in the forward direction. Furthermore, the 6 reactions are constrained by mass balance, and we note that two linearly independent rate vectors will characterize the fluxes of all 6 reactions in the network. Assume that the steady states metabolite concentrations (A 0 , B 0 , C 0 , O 0 1 , O 0 2 ) of the system are (1, 1, 1, 1, 1). Furthermore, assume that a flux of magnitude αF enters the cycle through reaction v 1 , and (1 − α)F flux returns through reaction v 5 (where 0 < α < 1).
Finally, we impose that there is no complicated activation or inhibition present in the system. This means that v 2 uses only metabolite A and cofactor O 1 as its substrates, v 3 uses only metabolite B, v 4 and v 5 use only metabolite C, v 6 uses only O 2 , and v 1 is constant and uses none of the metabolites as substrates. Notice that this constrains all θ parameters to be greater than or equal to zero. Constraining θ to be positive ensures that any kinetic rate laws we consider in our analysis are monotonic, i.e. that an increase in the concentration of the substrate of a reaction will never decrease the rate of a reaction. This will be essential in the analysis to follow. Now, we can write our system in terms of SKM variables: Because the cofactors O 1 and O 2 are conserved, we can remove them from the two matrices above (for details and justification, see the Supplementary Material of [9]). Doing this, we omit the last row of Λ, and θ becomes Then, to obtain the Jacobian we compute the product Λθ. In order to determine the stability of the system, we find the eigenvalues of the Jacobian and determine whether any of them are positive. We find the eigenvalues by imposing |J − λI| = 0, i.e. by setting to zero the following determinant: This corresponds to the equation The stability properties will be determined by the signs of the real parts of the eigenvalues obtained from solving this equation. For all eigenvalues with only real components, it is possible to determine the signs without explicitly solving the equation, by using Descartes' Rule of Signs. If we arrange the characteristic polynomial in order of decreasing variable exponents, then the number of positive real roots is equal to the number of signs changes in the coefficients of the ordered polynomial, or less than the number of size changes by a multiple of 2. In our case, there are no sign changes (all coefficients are greater than zero, notice that any negative terms cancel), so there are no positive roots. Note that it is not possible to have a zero eigenvalue either, because there is a nonzero coefficient for the λ 0 term. Then, we must consider two alternative cases. In one, the Jacobian has only real eigenvalues, in which case they are all negative and the system is stable. In the other, the Jacobian may contain pairs of complex conjugate eigenvalues whose real part will determine the stability of the system. To investigate systems which may display such complex eigenvalues, we sampled many instances of the Jacobian computationally and computed the corresponding eigenvalues. Any eigenvalues with real part greater than zero would indicate instability. For this particular cycle, we randomly sampled each θ in the range (0,1] for a value of α in the open set (0,1). Then, the eigenvalues of the Jacobian were calculated. This process was repeated 4 × 10 6 times and in no cases were any unstable Jacobians identified.

A More General Treatment
In this section we generalize our previous results to a metabolic cycle of arbitrary length, relying on the following assumptions: 1. There exists at least a single steady state for the system. It is trivial to show that such a steady state exists because rate constants (such as V max ) are not related to any θ parameters which determine the stability, and can be chosen independently and arbitrarily in order to satisfy the steady-state condition. Furthermore, results from [4] and computational analysis using the ERNEST MATLAB package [17] demonstrate that the cycle under study is strongly sign determined (SSD) and injective. It therefore admits only a single unique equilibrium for a given set of kinetic parameters (e.g. rate constants, saturation constants, etc.). It is worthwhile to note that if multiple equilibria did exist, the results presented in this section would still be valid. 2. The input and output reactions from the cycle are separated by one metabolite. In other words, if the output reaction of the cycle leaves through metabolite x i in the cycle, then the input reaction to the cycle enters through metabolite x i+1 . In analogy to the simple example above, the output reaction leaves through metabolite C, and the input reaction enters in the next metabolite in the chain, A. This assumption is made to simplify the analysis, but computational experiments indicate it does not appear to be necessary. However, this has not been explicitly proven. 3. The j th reaction (j = 2, 3, ..., N + 1) uses only the (j − 1) th metabolite as a substrate. There is no feedback inhibition or activation of any kind. This means that all θ parameters are constrained to be greater than zero. A negative θ would correspond to a metabolite having an inhibitory effect on the rate of reaction. 4. Cofactors such as ATP and NADH are often used as the energetic driving forces in metabolic reactions in order to make these reactions thermodynamically feasible. We assume that a single cofactor pair (O 1 and O 2 ) is involved in the metabolic cycle, and that it is not involved in the reaction creating the first or last metabolite in the cycle. The reformation of cofactors is assumed to be driven by some energy input into the system. Our goal is to study the stability of the cycle shown in figure 2, which generalizes the cycle of Figure 1.
The generalized form of Λ and θ, analogous to the ones from (3.1) and (3.3), is shown in Appendix Appendix C. Note that we make no further assumptions about the steady states of the system: X 0 = (X 0 1 , X 0 2 , . . . , X 0 N ) and the flux through the system remains equal to F . With N metabolites and N + 3 reactions, Λ is N × (N + 3) and θ is (N + 3) × N . The complete analysis is shown in Appendix C, and the results are summarized below.
The characteristic equation for the eigenvalues of the system reduces to The terms on the bottom two lines will cancel with terms from the top two lines, resulting in a polynomial equation of order N + 1 containing all nonzero coefficients with exactly the same sign. Furthermore, a constant λ 0 term will remain, indicating that there is no zero eigenvalue. Therefore, all eigenvalues with only real part are negative, and the system can only lose stability in an oscillatory manner.
We proceeded with computational experiments in which SKM parameters were randomly sampled to search for cases of systems whose eigenvalues lie in the positive real part of the complex plane. In these experiments, we generated random samples of θ parameters (in the range (0,10] to account for high Hill coefficients), α in (0,1), metabolite concentrations in (0,100), and calculated the eigenvalues of the Jacobian. Approximately 4 × 10 6 different instances were simulated for each particular cycle length from 3 to 8, and for each possible cofactor location. In none of the sampled instances was the real part of any eigenvalue of the Jacobian greater than zero, indicating that all samples were stable.
In a separate computational experiment, we tested the dependency of our stability results on the structure of the cycle. An example is shown in the lower part Figure 3 in which the first metabolite activates the reaction creating the third metabolite (similar to a feedforward loop). This introduces an instability into the system which does not appear in our prototypical system.
Thus, both computational experiments and analytical results support the conjecture that under general kinetic rate laws, cycles of this form are stable under all conditions. We believe that a proof of this conjecture may follow from methods similar to those used to establish stability results for linear chains with positive feedback, as presented in [16].

Autocatalytic Cycles
The next question we want to address is whether the methods and results we derived for the cycles of metabolic reactions seen in the previous sections can be extended to cycles of fundamentally different topology, namely autocatalytic cycles. Autocatalytic cycles consist of a set of metabolites which reproduce themselves with each turn of the cycle. Their simplified stoichiometry often takes the form A + X → 2A, where a substrate X transforms into a product A, through catalysis of A itself [1,18]. The formose reaction and the reductive citric acid cycle are two cases of autocatalytic cycles; both have been hypothesized to play roles in early metabolism, and significant work has been devoted to understanding the conditions under which they and other autocatalytic cycles may have naturally appeared and proliferated [2,19].
Shown in Figure 4, our simple example contains 3 metabolites and 5 reactions. It is distinct from the cycles shown in prior sections because both the first and last metabolite (A and C) are required for reaction v 2 to occur. Furthermore, we assume that reactions exiting from the cycle transport metabolites to a separate compartment, where they are not accessible to the reactions within the cycle. Notably, we will show that this system does not contain a universally stable steady state. Instead, the steady state remains stable under certain parameter regimes, to be defined in the analysis below. Although this example is relatively simple, the results can be extended to more complicated examples using a method similar to the one described in Appendix C.
The same assumptions are made as before, and all reactions are assumed to proceed in only the forward direction. There are 2 independent flux vectors which span the null space of the stoichiometric matrix for this system. Therefore, we assign the flux vector where γF is the outflow from the cycle through the reaction v 4 and F (x − γ − 1) is the outflow from the cycle through the reaction v 5 . Furthermore, it is noted that for flux balance considerations, x > γ. Since all reactions are assumed to go in the forward direction, x > γ + 1.
Following the same analysis as before, we can write the matrices Λ and θ: Figure 3: (a) Histogram of values of the maximal eigenvalue for a case of the system shown in Figure 2. The specific system contains 6 metabolites, and the reaction v 2 is coupled to a cofactor. All eigenvalues are stable, indicating the system is stable. (b) Histogram of values of the maximal eigenvalue for a system in which one of the assumptions is broken. In this case, metabolite X 1 activates reaction v 3 . Notice that not all maximal eigenvalues are negative, indicating that the system may be unstable for certain combinations of parameters. This suggests that global stability strongly relies on the topology of the network, and that additional activating or inhibitory reactions may prevent a steady state from being stable.
Then, the Jacobian J is the product of Λ and θ: By imposing |J − λI| = 0, we get the characteristic equation Expanding this equation, the coefficients of the third-order polynomial in λ are: To analyze the signs of the roots of this polynomial, we use the Routh-Hurwitz criterion. According to this criterion, a third-order polynomial in z with coefficients, a, b, c of the form z 3 + az 2 + bz + c = 0 has roots lying in the negative real half of the complex plane if the conditions a > 0, c > 0, ab − c > 0 are satisfied. For our polynomial, it is obvious that both a and c are always greater than zero as long as x > γ + 1, so we are left to find conditions where the last inequality is satisfied. Note that if x < γ + 1, then c < 0, there is a positive root to the characteristic equation, and the system is unstable.
Rather than deriving a precise condition for stability in terms of any single variable, we can use simple observations regarding the particular forms of the coefficients to infer trends in stability. In particular, notice that all negative terms in −c are precisely canceled out by corresponding terms in ab. Then, the only term which may cause ab − c to be less than zero is (x−γ)θ2θ4(1−x)
We can computationally test these predictions by following a procedure similar to [9] and randomly sampling the parameters that define the autocatalytic cycle. Specifically, we can calculate the eigenvalues of the Jacobian matrix as a single parameter of interest (such as a θ parameter) is varied, while all other θ parameters and metabolite concentrations are sampled randomly. By doing so, we sample over many possible instances of the metabolic network, and can identify correlations between individual parameters and stability. For each combination of dependent variable (e.g.: x, θ 1 , θ 2 , A 0 ) and γ shown in the figures below, 5 × 10 5 instances of the Jacobian were sampled and characterized. The value of the maximal eigenvalue for each instance indicated whether the network is stable or not.
From the terms of the λ 0 coefficient derived above, one can see that as long as x > γ + 1, it is ideal for γ to be as large as possible and x to be as small as possible. This minimizes the value of . In Figure 5(a), the stability of the system is plotted against values of the stoichiometric coefficient x of metabolite B in reaction v 2 . Each line indicates a different value for the outflow γ. As x increases, the stability of the system falls. Additionally, as the value of γ increases, the stability of the system increases, while the general shape of the plot remains identical. Note that if x < γ + 1, then the system is universally unstable, in agreement with the condition derived before. It is particularly interesting that the system's optimally stable operating point is located at the boundary of a region where it is entirely unstable due to dilution effects. Similarly, Figure 5(b) confirms that as θ 1 increases, the stability of the system rises. In Figure 5(c), it is shown that as θ 2 increases, the stability of the system falls. Furthermore, the relationship between stability and γ is retained in all these plots.
The combination of analytical and computational results presented so far in this section indicates that the tuning of a single parameter can greatly increase the likelihood of a metabolic system being stable. One may wonder which parameters, among those that can affect stability, are the ones that a cellular system would be most likely to adjust. An intriguing possibility supported by the current findings is that a cell could approach a stable steady state by appropriately tuning metabolite concentrations. A sufficiently high steady state concentration of B or C would nearly guarantee that the simple autocatalytic cycle would be stable. In the future it may be interesting to verify experimentally whether certain metabolic pathways -and autocatalytic cycles in particular -may display steady state metabolite concentrations that fall into the range of high stability.

Discussion
We have used structural kinetic modeling (SKM) to assess the dynamical stability of certain classes of metabolic cycles. In prior work, it was demonstrated that SKM can be used as a powerful computational tool for identifying stabilizing sites in metabolic networks by sampling large ranges of parameters [7]. The main thrust of our work was to show that SKM can also be used analytically to characterize stability. Drawing from basic concepts in nonlinear dynamics and linear algebra, we demonstrated that it is possible to make a strong assertion about the stability of a class of non-autocatalytic metabolic cycles, and to identify trends in stability as a function of model parameters for a type of autocatalytic cycles.
For the non-autocatalytic cycles of the type illustrated in Figure 2, we have proven that the topology of the networks ensures that they can only lose stability in an oscillatory manner, regardless of metabolite concentrations, flux magnitude, or form of kinetic equations. The generality of the result with respect to kinetic laws is particularly noteworthy. Computational experiments support the conjecture that these cycles are completely stable and do not produce Jacobians with positive complex eigenvalues. The stability of non-autocatalytic cycles similar to the one studied in this section could be easily proven using the Deficiency One Theorem of Chemical Reaction Network Theory (CRNT), [6]. However the CRNT result, in this case, would be restricted to cases of simple mass action kinetics. Hence, our work shows that an analytical approach completely distinct from CRNT converges to the same strong result, and extends it to a large category of kinetic rate laws. By approaching the question of stability with an analytical flavor of SKM, we demonstrate that the dynamics of a well-organized metabolic network, even in the presence of nonlinear kinetic rate laws, is not necessarily beyond "pencil and paper" methods, and that useful and significant conclusions can be drawn from such analyses. Additionally, recent work has shown that a broad category of regulatory and biochemical networks can be analyzed using the theory of monotone systems [20]. In the future it will be interesting to determine whether this theory applies to the cycles discussed here, and conversely how our analytical use of SKM may tie into the theory of monotone systems.
In applying the same SKM techniques to a simple autocatalytic cycle (Figure 4), we found that a combination of analytically derived conditions and computational experiments could highlight the dependence of stability on metabolite concentrations, fluxes, and degrees of saturation. In this case, the stability of the cycle was not guaranteed, and the analysis suggested possible directions that evolution may take in order to steer the network towards an even more robust steady state. While it was possible to derive expressions equivalent to Eq.(5.4) without employing the SKM non-dimensionalization, it should be obvious that the algebra required to do so would be complicated and non-trivial assuming the rate laws were nonlinear. SKM provided an appealing (but nevertheless rigorous and useful) workaround, and still produced the correct characteristic equation for the Jacobian. To date, we are not aware of a similar effort to characterize the dynamic properties of autocatalytic cycles for general kinetic rate laws.
Our results suggest that a system does not necessarily require complex forms of regulation in order to maintain a stable steady state. In fact, in many cases it seems plausible that the elevated concentration of a single metabolite may be sufficient to endow stability. Why is it then that we observe many complex forms of allosteric (and genetic) enzyme regulation among biochemical networks? Certainly, such regulation may play a crucial role in stabilizing more complex enzymatic networks, but perhaps there is more to the story. One possibility may be that these forms of regulation actually control the transition of the system from one steady state to another. When conditions in the environment change, it may be favorable for a cell to alter its metabolism to react to these changes. Allosteric and transcriptional regulation may determine the manifolds along which the transitions between two steady states actually occur.
It is also interesting to consider the possibility that the inherent stability of certain types of pathways may influence the evolution of biochemical networks. One possible link between dynamic stability and evolution might be the adaptive fine tuning of kinetic parameters to guarantee stability of a given pathway topology. Conversely, we are suggesting that natural selection might gradually amplify the preference for network topologies with an intrinsically high probability of being stable. Such topologies would be robust and not require substantial fine tuning of kinetic parameters [21,22]. Defining the precise parameter ranges and topological structures which promote stability among these processes may enable one to actually search for these common motifs across a broad range of well-studied metabolic networks. As more data is collected characterizing the rates of reactions and magnitudes of concentrations in metabolism, the possibility of such a search becomes more likely. Searching for dynamical behaviors by varying topologies and parameter ranges is quickly becoming a prototypical problem in systems biology. Work such as [23,24] demonstrate that solutions to these problems can have quite fruitful repercussions in our ability to understand naturally occurring biological systems and synthetically design novel genetic and signalling networks.
Finally, the observation that metabolic cycles have robust steady state properties may have ramifications in the study of early metabolism. Several studies have suggested that the early stages of molecular self-organization and self-reproduction may have involved the emergence of self-sustaining proto-metabolic cycles [1,2]. Arguments such as those presented in [25], however, remind us that the earliest biological reactions may have taken place in environments which siphoned off the products of these reactions. In the case of cycles, these dilutions would have potentially precluded a replenishment of substrates in the cycle, and thus prevented a stable steady state. Irrespective of this open question, our work suggests that when autocatalytic cycles did first appear, they would have benefitted from inherent robustness to environmental perturbations. In general, these cycles may be only a single instance of a large array of other metabolic topologies which exhibit equally robust stability properties. Future studies may address whether these stability properties played a role in establishing the core of current cellular metabolism.

Appendix A. Structural Kinetic Modeling (SKM)
The methods used in this paper are based on the SKM framework [9]. SKM is a specific application of generalized modeling [10] in which normalized parameters replace conventional kinetic parameters (such as V max or K M ) in the modeling of metabolic networks. The normalized parameters have a direct connection to the original kinetic parameters, but are much easier to work with. As will be shown below, these parameters usually have well-defined and extremely limited ranges (e.g. [0,1]), and sampling them across this range effectively samples all possible values of the original kinetic parametsrs.
The goal of SKM is to capture the local stability properties of a biochemical system. In this sense, it bridges the gap between genome-scale steady state modeling [11] and explicit kinetic modeling of a metabolic process. To accomplish this, one usually determines the Jacobian J of the system of interest and evalutes it at its equilibrium. Assuming knowledge of the form of kinetic rate laws for a CRN, it is quite possible to write down the corresponding J. However, in many cases J will be unnecessarily complicated and quite difficult to work with. By performing a change of variables, SKM actually simplifies the mathematical form of each entry in J. The new entries are in almost all cases easier to work with. Much of this material is taken directly from [9] and we refer readers to that reference and its supplementary materials for more information.
If we let S be the m-dimensional vector of metabolite concentrations, N be the m ×r stoichiometric matrix, and v be the r -dimensional vector of reaction rates, then we can describe the dynamics of the system with the equation where v(S, k) denotes that the reaction rates are dependent on both metabolite concentrations S and unknown parameters (such as Michaelis-Menten constants) k. If we assume that a non-negative steady state S 0 exists, then we can redefine our system using the definitions where i = 1...m and j = 1...r. Now, x is a vector of metabolite concentrations normalized with respect to their steady state concentrations and µ represents flux normalized with respect to steady state flux values. The matrix Λ represents the stoichiometric matrix normalized with respect to steady state fluxes and steady state metabolite concentrations. We can now rewrite the system of differential equations as Now, we can write the Jacobian and evaluate it at x 0 = 1 (which, because of the way x is defined, is actually the equilibrium of the system). Calculating the Jacobian, we find that If we evaluate (A.4) at x 0 = 1 we find Note that the equations above were derived without regard to the actual form of the kinetic equations that determined the ODE system. The matrix θ µ x contains elements which represent the degree of saturation of normalized flux µ j with respect to normalized substrate concentration x i . In terms of derivatives, each element of θ represents the degree of change in a flux as a particular metabolite is incrementally increased. This is analogous to the concept of elasticity in metabolic control analysis [12].
What does the θ matrix look like? Its columns correspond to each metabolite, and its rows to each flux. A non-zero element θ i j in the matrix represents the effect a small change in a metabolite j has on flux i. In the case of Michaelis-Menten kinetics, this element in the matrix may take values ranging from [0,1]. In the case of standard competitive inhibition (e.g. allosteric inhibition by a product), the element takes values in [-1,0].
To usefully illustrate the meaning of a single θ element, consider an equation following Michaelis-Menten kinetics shown in (2.4). First, we write µ(x), which we recall is the flux normalized by the flux at the steady state. Here, S 0 is the concentration of the substrate S at steady state. Then, we manipulate the equation by substituting xS 0 for S, where x is the normalized steady state concentration of substrate S. The result is shown in (2.5).
Finally, we take a derivative with respect to x and evaluate it at x = 1 to obtain θ in (2.6). Notice that θ can only take values between 0 and 1 for any positive value of S 0 .
It is also possible to approximate the value of θ for different extremes of saturation. In the regime of extremely low substrate concentrations when the enzyme catalyzing the reaction is not saturated at all, formation of product can be approximated by which is approximately linear in S. In this case, θ i j would be approximately 1. This is identically the case for mass-action kinetics, where the rate of product formation is linear with the amount of substrate. In contrast, at extremely high substrate concentrations, In this case, S has nearly no effect on the rate of product formation, and so θ i j is approximately 0. For all other substrate concentrations, θ falls between 0 and 1. Note that the above analysis can be performed for other forms of kinetic equations including mass-action and Hill-type equations. In these cases, the values for θ i j are similarly constrained. For example, in the case of Hill kinetics, θ i j falls in the range [0,n], where n is the Hill constant. For more information on analyzing other types of kinetics, see the supplementary materials in [9]. The power of SKM comes from the parameterization illustrated above. Each element in θ has a precise correspondence to some combination of kinetic parameters in the original model. However, it is far more tractable to study a system using θ parameters rather than the original kinetic parameters. To see this more clearly, recall that in most cases, biochemical kinetic constants are poorly estimated. In order to build precise ODE models for biochemical systems, it is usually necessary to actually choose values for these constants. While the chosen values may be estimated from experimental measurements, hidden dependencies in the model may actually results in non-obvious correlations between parameters that can strongly affect the output of the model.
On the other hand, if one is only interested in the stability of a characterized steady-state, it may not be necessary to actually have knowledge of kinetic parameters. Experimental measurements can provide data on absolute metabolite concentrations and flux values, and the stoichiometry is often known a priori. Then, one can parameterize the system as shown above and sample many possible combinations of θ parameters. For each unique set of θ parameters, the stability of the Jacobian is determined and recorded. Analysis of the results can lead to several important conclusions such as which θ parameters have the strongest correlation with stability of the entire metabolic system. Thus, the benefit of employing SKM over other techniques is that it provides the means to analyze and make sense of a large number of possible cases of a metabolic network, rather than a single instance.
Traditionally, SKM has been used as a computational tool to reveal relationships between the stability of a metabolic module and the parameters which govern its kinetics. However, many aspects of the SKM parameterization lend themselves to analytical methods. By using θ parameters, it is no longer necessary to decide on the kinetic form (e.g. mass action, Michaelis-Menten, Hill) of the chemical reactions taking place. Instead, the problem of studying stability is reduced to finding the signs of the roots of a polynomial equation (the characterstic equation of the Jacobian J), whose coefficients are determined by the fluxes, absolute metabolite concentrations, and θ's of the system at steady state. Moreover, it is in many cases possible to limit the values of θ (e.g.-if simple Michaelis-Menten kinetics are assumed without inhibition, θ falls between 0 and 1). If the signs of the roots are all negative, then the eigenvalues of the Jacobian are negative and the current steady state is stable.
At first glance, it is unclear whether this will make things easier or harder; values for fluxes and concentrations are almost never known. However, our work shows that the absence of this knowledge does not make the problem of stability intractable. Instead, we operate with the following idea: we assume that the steady state flux vector v and the metabolite concentration vector X are variables, and we do not specify them. Instead, we simply determine the stability of the system in terms of these variables. Then, we can find trends in stability as the flux and metabolite concentrations change: it may be that as a concentration goes up, the stability of the system tends to fall. Isolating these trends is the motivation of our analysis.
It is worthwhile to clarify what it means to change a steady-state parameter like a flux or a θ. At first glance, it may not be obvious that it is even possible to alter such a parameter. In other words, suppose that at steady state, the concentration of metabolite A in some network is 1. What does it mean to increase the concentration of A? If we do so without altering any other parameters, the metabolic network will in all likelihood not be in steady state any longer. However, there may exist (and in fact for the systems described in this paper there does exist) another set of concentrations, fluxes, and θ's which does lead to a steady state at the increased concentration of A. The problem, again, is that we do not know their actual values. However, we may be able to infer, based on changes in the coefficients of the characteristic equation, that if we were to sample a large number of unique instances of the Jacobian for both concentrations of A, the higher concentration may actually be more likely to be stable. This is in fact the method we will use to study autocatalytic cycles later in the paper. However, in the next two sections, we will demonstrate that for certain types of cycles, it is not necessary to go to such great lengths in order to characterize stability; it will come much more naturally.

Appendix B. Linear Chains
Up to this point, linear metabolic networks have been ignored. The reason for this is that, using the same analysis as above, it is trivial to show such a network is stable. In this appendix, we provide a short example illustrating this point. Consider the metabolic network shown in Figure 3 with flux F flowing through it and observed steady states (A 0 B 0 ).
If we follow the same analysis as before (writing Λ and θ and calculating the Jacobian J), we will find The eigenvalues of the Jacobian are then −F A 0 and −F B 0 , both of which are negative. Therefore, the system is stable.

Appendix C. Cofactors
The reaction system shown below is the complete, more general version of the simple system analyzed above. Now, we can write a system very similar to the one in Section 2. Consider the example shown in the figure below. The generalized form of Λ and θ is shown below. Note that we make no further assumptions about the steady states of the system: X 0 = (X 0 1 , X 0 2 , . . . , X 0 N ) and the flux through the system remains equal to F . With N metabolites and N + 3 reactions, Λ is N × (N + 3) and θ is (N + 3) × N .
Note that because the cofactors come in a conserved pair, one of them can be removed from the system. This modifies the entry in the bottom right element of the Jacobian.
Note that the last row of Λ is omitted because the cofactors come as a conserved pair. Furthermore, note that the bottom right element of θ is replaced with a negative element in order to account for this conservation. Aftering factor out F , the Jacobian becomes Next, we evaluate take the determinant |J − λI| and take the Laplace expansion over the last column. This yields three smaller determinants. The first is −θN+2 This determinant is easily evaluated and results in The second matrix resulting from the Laplace expansion of the original Jacobian is (−1) N +R+1 −θN+2 If another Laplace expansion is now taken along the N th row, then the resulting matrix is The resulting term for this determinant is The third matrix resulting from the Laplace expansion is (−1) N +R+2 θN+2 X 0 R+1 multiplied by the determinant of all nonzero coefficients with exactly the same sign. Furthermore, a constant λ 0 term will remain, indicating that there is no zero eigenvalue.