Design principles of autocatalytic cycles constrain enzyme kinetics and force low substrate saturation at flux branch points
 Cited 7
 Views 1,124
 Annotations
Abstract
A set of chemical reactions that require a metabolite to synthesize more of that metabolite is an autocatalytic cycle. Here, we show that most of the reactions in the core of central carbon metabolism are part of compact autocatalytic cycles. Such metabolic designs must meet specific conditions to support stable fluxes, hence avoiding depletion of intermediate metabolites. As such, they are subjected to constraints that may seem counterintuitive: the enzymes of branch reactions out of the cycle must be overexpressed and the affinity of these enzymes to their substrates must be relatively weak. We use recent quantitative proteomics and fluxomics measurements to show that the above conditions hold for functioning cycles in central carbon metabolism of E. coli. This work demonstrates that the topology of a metabolic network can shape kinetic parameters of enzymes and lead to seemingly wasteful enzyme usage.
https://doi.org/10.7554/eLife.20667.001eLife digest
Many bacteria are able to produce all the molecules they need to survive from a limited supply of nutrients. This allows the bacteria to thrive even in harsh environments where other organisms struggle to live. The bacteria act as miniature chemical factories to convert nutrients into the desired molecules via a series of chemical reactions. Some molecules are made in sets of reactions termed autocatalytic cycles. These reaction sets require a molecule to be present in the cell in order to produce more of that molecule; like how a savings account needs to contain some money before it can generate more via interest.
Bacteria have many different enzymes that each drive specific chemical reactions. In order for an autocatalytic cycle to work properly, the cell needs to maintain adequate supplies of the molecule it is trying to make and all of the “intermediate” molecules in the cycle. If less of an intermediate molecule is produced, for example, the cell needs to reduce the demand for that molecule by controlling later chemical reactions in the cycle. Bacteria control chemical reactions by regulating the activities of the enzymes involved, but it is not clear exactly how they regulate the enzymes that drive autocatalytic cycles.
Barenholz et al. combined two approaches called proteomics and fluxomics to study autocatalytic cycles in a bacterium known as E. coli. The experiments suggest several core principles allow autocatalytic cycles to work smoothly in the bacteria. The next step is to apply these principles to different kinds of molecules produced in bacterial cells. A future challenge is to search for other structures that regulate chemical reactions in E. coli and other bacteria. Extending our understanding of autocatalytic cycles and other pathways of chemical reactions is essential for designing and engineering new reactions in bacteria. Such knowledge can be used to modify bacteria to produce valuable chemicals in environmentally friendly ways.
https://doi.org/10.7554/eLife.20667.002Introduction
An essential trait of living systems is their ability to reproduce. This fundamental ability makes all living organisms autocatalytic by definition. Moreover, autocatalytic metabolism is considered to be one of the essential components of life (Ganti et al., 2003).
In this work, we focus on autocatalytic cycles in chemical reaction systems, in the context of metabolic networks. The components we consider are the metabolites of the system, with autocatalytic cycles being formed using the reactions of the metabolic network. An illustrative example for a metabolic autocatalytic cycle is glycolysis. In glycolysis, 2 ATP molecules are consumed in the priming phase, in order to produce 4 ATP molecules in the pay off phase. Therefore, in order to produce ATP in glycolysis, ATP must already be present in the cell. Subsequently, autocatalysis of ATP in glycolysis (also referred to as ‘turbo design’) results in sensitivity to mutations in seemingly irrelevant enzymes (Teusink et al., 1998). Autocatalytic cycles have also been shown to be optimal network topologies that minimize the number of reactions needed for the production of precursor molecules from different nutrient sources (Riehl et al., 2010).
Metabolic networks often require the availability of certain intermediate metabolites, in addition to the nutrients consumed, in order to function. Examples of obligatorily autocatalytic internal metabolites in different organisms, on top of ATP, are NADH, and coenzyme A (Kun et al., 2008). We find that other central metabolites, such as phosphosugars and organic acids, are autocatalytic under common growth conditions. The requirement for availability of certain metabolites in order to consume nutrients implies metabolic processes must be finely controlled to prevent such essential metabolites from running out; in such cases metabolism will come to a halt. Autocatalytic cycles present control challenges because the inherent feedback nature of autocatalytic cycles makes them susceptible to instabilities such as divergence or drainage of their intermediate metabolites (Teusink et al., 1998; Fell et al., 1999; Reznik and Segrè, 2010). The stability criteria typically represent one constraint among the parameters of the cycle enzymes. For large cycles, such as the whole metabolic network, one such constraint adds little information. For compact autocatalytic cycles embedded within metabolism, one such constraint is much more informative. We thus focus our efforts on analyzing small autocatalytic cycles. Finding the unique constraints that metabolic autocatalytic cycles impose is essential for understanding the limitations of existing metabolic networks, as well as for modifying them for synthetic biology and metabolic engineering applications.
A key example of an autocatalytic cycle in carbon metabolism is the CalvinBensonBassham cycle (CBB) (Benson et al., 1950). The carbon fixation CBB cycle, which fixes CO_{2} while transforming fivecarbon compounds into two threecarbon compounds, serves as the main gateway for converting inorganic carbon to organic compounds in nature (Raven et al., 2012). The autocatalytic nature of the CBB cycle stems from the fact that for every 5 fivecarbon compounds the cycle consumes, 6 fivecarbon compounds are produced (by the fixation of 5 CO_{2} molecules). Beyond the CBB cycle, we show that most of the reactions and metabolites in the core of central carbon metabolism are part of compact (i.e. consisting of around 10 reactions or fewer) metabolic autocatalytic cycles. Some of the autocatalytic cycles we find are not usually considered as such. The span of autocatalytic cycles in central carbon metabolism suggests that the constraints underlying their stable operation have networkwide biological consequences.
In this study, we present the specific requirements metabolic autocatalytic cycles must meet in order to achieve at least one, nonzero, steady state which is stable in respect to fluctuations of either metabolites or enzyme levels close to the steady state point. The mathematical tools we use are part of dynamical systems theory (Strogatz, 2014). We identify the kinetic parameters of enzymes at metabolic branch points out of an autocatalytic cycle as critical values that determine whether the cycle can operate stably. We show that the affinity of enzymes consuming intermediate metabolites of autocatalytic cycles must be limited to prevent depletion of these metabolites. Moreover, we show that the stable operation of such cycles requires low saturation, and thus excess expression, of these enzymes. Low saturation of enzymes has previously been suggested to stem from a number of reasons in different contexts: (A) to achieve a desired flux in reactions close to equilibrium, for example in glycolysis (Staples and Suarez, 1997; Eanes et al., 2006; Flamholz et al., 2013); (B) to provide safety factors in the face of varying nutrient availability, for example in the brushborder of the mouse intestine (Weiss et al., 1998); (C) to accommodate rapid shifts in demand from the metabolic networks in muscles with low glycolytic flux (Suarez et al., 1997); (D) to allow fast response times, for example to pulses of oxidative load in erythrocytes, resulting from their adherence to phagocytes (Salvador and Savageau, 2003). Our findings add to these reasons the essential stabilizing effect of low saturation of branch reactions on the stability of fluxes through autocatalytic cycles.
We use recent fluxomics and proteomics data to test the predictions we make. We find them to hold in all cases tested where autocatalytic cycles support flux. Our analysis demonstrates how the requirement for stable operation of autocatalytic cycles results in design principles that are followed by autocatalytic cycles invivo. The results and design principles presented here can be further used in synthetic metabolic engineering applications that require proper functioning of autocatalytic cycles.
Results
Compact autocatalytic cycles are common and play important roles in the core of central carbon metabolism
Different definitions exist for autocatalytic sets in the context of chemical reaction networks (Hordijk and Steel, 2004; Eigen and Schuster, 2012; Kun et al., 2008). Here we define an autocatalytic cycle as a set of reactions and metabolites that form a cycle, and that, when the reactions are applied to the substrates at the given stoichiometric ratios, increase the amount of the intermediate metabolites. A minimal example of a metabolic autocatalytic cycle is shown in Figure 1, where an internal metabolite joins with an external assimilated metabolite to give rise to $1+\delta $ copies of the internal metabolite, representing an increase by $\delta $ copies. For stable operation, $\delta $ copies have to branch out of the cycle, and this consumption must be robust to small fluctuations in enzyme levels and metabolite concentrations. For a formal, mathematical definition, see Materials and methods section "Formal definition of an autocatalytic metabolic cycle".
While rarely discussed as such, a systematic search in the central carbon metabolism core model of E. coli (see Materials and methods section "Systematic identification of autocatalytic cycles in metabolic networks") shows the ubiquity of compact autocatalytic cycles. On top of the previously discussed CBB cycle (Figure 2, example I), we show two other prominent examples:
The glyoxylate cycle within the TCA cycle, which turns an internal malate and two external acetylCoAs into two malate molecules. This is achieved by transforming malate to isocitrate, while assimilating acetylCoA, and then splitting the isocitrate to produce two malate molecules, assimilating another acetylCoA (Kornberg, 1966) (Figure 2, example II).
A cycle formed by the glucose phosphotransferase system (PTS) in bacteria. This transport system imports a glucose molecule using phosphoenolpyruvate (pep) as a cofactor. The imported glucose is further catabolized, producing two pep molecules via glycolysis (Figure 2, example III).
Two additional examples are presented in Figure 2—figure supplements 1 and 2 and discussed below.
The ubiquity of compact autocatalytic cycles in the core of central carbon metabolism motivates the study of unique features of autocatalytic cycles, as derived below, which may constrain and shape the kinetic parameters of a broad set of enzymes at the heart of metabolism.
Steady state existence and stability analysis of a simple autocatalytic cycle
To explore general principles governing the dynamic behavior of autocatalytic cycles, we consider the simple autocatalytic cycle depicted in Figure 3A. This cycle has a single intermediate metabolite, $X$. We denote the flux through the autocatalytic reaction of the cycle by ${f}_{a}$, such that for any unit of $X$ consumed, it produces two units of $X$. The autocatalytic reaction assimilates an external metabolite (denoted $A$), which we assume to be at a constant concentration. We denote the flux through the reaction branching out of the cycle by ${f}_{b}$. Biologically, ${f}_{b}$ represents the consumption of the intermediate metabolite $X$. In the cycles we find in central carbon metabolism, the branch reactions provide precursors that support growth through subsequent reactions. We thus also sometimes consider ${f}_{b}$ to represent biomass generation.
For simplicity in the derivation, we assume irreversible MichaelisMenten kinetics for the two reactions. Even though ${f}_{a}$ should follow bisubstrate velocity equation, assuming constant concentration of $A$ reduces the bisubstrate equation to a simple MichaelisMenten equation. The apparent kinetic constants of the equation depend on the constant value of $A$ (see Materials and methods section "Connecting bisubstrate reaction kinetic constants with simple MichaelisMenten constants"). We extend our analysis to bisubstrate reaction equations in the next section. We therefore assume that:
where ${V}_{\mathrm{max}}$ is the maximal flux each reaction can carry and ${K}_{M}$ is the substrate concentration at which half the maximal flux is attained. Physiologically, these kinetic parameters must be positive. Using these simple forms allows us to obtain an analytic solution. We discuss more general cases below.
We characterize the metabolic state of this system by the concentration of the metabolite $X$. We note that knowing the concentration of $X$ suffices in order to calculate the fluxes originating from it, ${f}_{a}$ and ${f}_{b}$, thus fully defining the state of the system. A steady state of the system is defined as a concentration, ${X}^{*}$, which induces fluxes that keep the concentration constant, such that the total influx to $X$ is exactly equal to the total outflux from it. In our example, the outgoing flux from $X$ is ${f}_{a}+{f}_{b}$ and the incoming flux to $X$ is $2{f}_{a}$, so at steady state it holds that:
Intuitively, at steady state, the branch reaction must consume all the excess intermediate metabolite that is produced by the autocatalytic reaction. Indeed, expanding the condition above gives:
which is satisfied either if ${X}^{*}=0$ or if:
implying that:
The concentration of $X$ cannot be negative, and thus we get a constraint on the kinetic parameters for which a positive steady state exists. Either both the numerator and the denominator of Equation 3 are positive, such that:
or both are negative, such that:
These constraints are graphically illustrated in Figure 3B, cases (III) and (I).
In order to gain intuition for this relationship we note that $\frac{{V}_{\mathrm{max}}}{{K}_{m}}$ is the slope of the Michaelis Menten function at $X=0$. The existence of a positive steady state can be used to get that:
The last inequality above implies that in order for a positive steady state to exist, the reaction with higher maximal flux must have a shallower slope at $X=0$. Mathematically, the constraint states that if $V}_{max,a}\text{}\text{}{V}_{max,b$ then $\frac{{V}_{max,a}}{{K}_{M,a}}\text{}\text{}\frac{{V}_{max,b}}{{K}_{M,b}}$. Alternatively, if $V}_{max,a}\text{}\text{}{V}_{max,b$ then $\frac{{V}_{max,a}}{{K}_{M,a}}\text{}\text{}\frac{{V}_{max,b}}{{K}_{M,b}}$. This condition can be intuitively understood, as the reaction with shallower slope at $X=0$ has smaller fluxes for small values of $X$ compared with the other reaction, so unless it has higher fluxes than the other reaction for large values of $X$ (meaning that its maximal flux is higher), the two will never intersect (see Figure 3B).
While having a steady state at positive concentration is an essential condition to sustain flux, it is not sufficient in terms of biological function. The steady state at positive concentration must also be stable to small perturbations. Stability with respect to small perturbations is determined by the response of the system to small deviations from the steady state, ${X}^{*}$ (at which, by definition $\dot{X}=0$). Assuming $X={X}^{*}+\mathrm{\Delta}X$, stability implies that if $\mathrm{\Delta}X$ is positive then $\dot{X}$ needs to be negative at ${X}^{*}+\mathrm{\Delta}X$, reducing $X$ back to ${X}^{*}$, and if $\mathrm{\Delta}X$ is negative, $\dot{X}$ will need to be positive, increasing $X$ back to ${X}^{*}$. It then follows that in order for ${X}^{*}$ to be stable, $\frac{d\dot{X}}{dX}\text{}\text{}0$ at $X={X}^{*}$, implying that upon a small deviation from the steady state ${X}^{*}$ (where $\dot{X}=0$), the net flux $\dot{X}$ will oppose the direction of the deviation.
For the simple kinetics we chose, the stability condition dictates that:
The analysis is straightforward for the case of ${X}^{*}=0$, yielding that 0 is a stable steady state concentration if $\frac{{V}_{max,b}}{{K}_{M,b}}\text{}\text{}\frac{{V}_{max,a}}{{K}_{M,a}}$, corresponding to the area above the diagonal in Figure 3B, where $\frac{{V}_{max,b}}{{V}_{max,a}}\text{}\text{}\frac{{K}_{M,b}}{{K}_{M,a}}$. These cases are denoted as cases (II) and (III). If the relation is reversed (i.e. $\frac{{V}_{max,b}}{{K}_{M,b}}\text{}\text{}\frac{{V}_{max,a}}{{K}_{M,a}}$), then 0 is an unstable steady state. The criterion that is of interest, however, is the criterion for stability of the nonzero steady state, ${X}^{*}=\frac{{V}_{\mathrm{max},b}{K}_{M,a}{V}_{\mathrm{max},a}{K}_{M,b}}{{V}_{\mathrm{max},a}{V}_{\mathrm{max},b}}$. In this case, substituting ${X}^{*}$ in Equation 4 gives the opposite condition to that of ${X}^{*}=0$. This steady state is thus stable if $\frac{{V}_{\mathrm{max},b}}{{K}_{M,b}}<\frac{{V}_{\mathrm{max},a}}{{K}_{M,a}}$, corresponding to the magenta domain in Figure 3B, containing cases (I) and (IV), and unstable otherwise.
The stability criterion can be generally stated in metabolic control terms (Fell, 1997) using the notion of elasticity coefficients of reactions, defined as:
In these terms, stability is obtained if and only if the elasticity of the branch reaction at the positive steady state concentration is greater than the elasticity of the autocatalytic reaction:
The complete analysis is summarized in Figure 3B. Domain (I) is the only domain where a positive, stable steady state exists. Domains (I) and (III) are the domains at which a positive steady state concentration exists, but in domain (III) that steady state is not stable. The domains below the diagonal (cases (I) and (IV)) are the domains where ${X}^{*}=0$ is an unstable steady state concentration, so that if another steady state exists, it is stable, but in domain (IV) no positive steady state exists. The domains above the diagonal (cases (II) and (III)) are the domains where ${X}^{*}=0$ is a stable steady state concentration, so that the other steady state, if it exists, is unstable.
Aside from existence and stability, a quantitative relationship between the affinity of the biomass generating, branching reaction and the flux it carries can be obtained. This relationship is opposite to the standard one, meaning that unlike the common case where the flux $f$ increases when the affinity becomes stronger, in this case, because the steady state concentration increases when ${K}_{M,b}$ becomes weaker (Equation 7 in Materials and methods section "Steady state concentration dependence on kinetic parameters of autocatalytic and branch reactions"), ${f}_{b}$ also increases when ${K}_{M,b}$ becomes weaker.
To conclude, for this simple cycle, we get that in order for a positiveconcentration stable steady state to exist (case (I)), two conditions must be satisfied:
The first requirement states that the maximal flux of the biomass generating, branching reaction should be higher than the maximal flux of the autocatalytic reaction. This requirement ensures a stable solution exists, as large concentrations of $X$ will result in its reduction by the branch reaction. The second requirement states that for concentrations of $X$ that are close enough to 0, the autocatalytic reaction is higher than the branch reaction (as can be inferred from the slopes). This requirement implies that the two fluxes will be equal for some positive concentration of $X$, ensuring a positive steady state exists. As this requirement further implies that below the positive steady state the branch reaction will carry less flux than the autocatalytic reaction, it follows that small deviations of the concentration of $X$ below the steady state will result in an increase in its concentration by the autocatalytic reaction, driving it back to the steady state. Meeting the second constraint has another consequence.
Interestingly, these conditions imply that if $K}_{M,b}\text{}\text{}{K}_{M,a$ then no positive stable steady state can be achieved. Specifically, changes to the expression levels of the enzymes catalyzing ${f}_{a}$ and ${f}_{b}$ only affect ${V}_{\mathrm{max},a}$ and ${V}_{\mathrm{max},b}$, and therefore do not suffice to attain a stable positive steady state. This indicates that stability of autocatalytic cycles, that are represented by the model analyzed above, depends on inherent kinetic properties of the enzymes involved and cannot always be achieved by modulating expression levels. We suggest this property to be a design principle that can be critical in metabolic engineering.
Integrating the bisubstrate nature of the autocatalytic reaction into the simple model
In the model above, to keep the analysis concise, we neglected the bisubstrate nature of the autocatalytic reaction. We extend the analysis to the most common classes of bisubstrate reaction mechanisms in the Materials and methods, sections "Connecting bisubstrate reaction kinetic constants with simple MichaelisMenten constants", "Constraints on concentration of assimilated metabolite and kinetic constants of bisubstrate reactions", and "Dependence of steady state concentration on assimilated metabolite". All bisubstrate reaction schemes analyzed take a MichaelisMenten like form once the concentration of the assimilated metabolite is kept constant (Equations 8, 10, 12, and 14 in Materials and methods section "Connecting bisubstrate reaction kinetic constants with simple MichaelisMenten constants").
For any set of kinetic parameters, under all ternary enzyme complex schemes, a lower bound on the concentration of $A$ exists, under which the conditions for the existence and stability of a positive steady state cannot be satisfied (Equations 18 and 23 in Materials and methods section "Constraints on concentration of assimilated metabolite and kinetic constants of bisubstrate reactions"). The exact value of the minimal concentration of $A$ depends on the specific bisubstrate reaction scheme and the kinetic parameters of it.
In the simplified model analyzed above, stability implied the affinity of the branch reaction towards its substrate was limited. A similar limit exists in most cases of bisubstrate reaction schemes (Equations 16, 19, and 21 in Materials and methods section "Constraints on concentration of assimilated metabolite and kinetic constants of bisubstrate reactions"). Interestingly, if the bisubstrate reaction is ordered with the internal metabolite binding first, then no strict constraints exist on ${K}_{M,b}$ and a stable steady state solution can always be achieved by setting appropriate values to ${V}_{\mathrm{max},b}$ and ${V}_{\mathrm{max}}$, the maximal flux of the bisubstrate autocatalytic reaction.
Finally, regarding the dynamic behavior of the system when the concentration of $A$ varies, we note that in all three ternary enzyme complex cases, as the concentration of $A$ approaches its lower bound, the steady state concentration of $X$ approaches 0, reducing both the autocatalytic and the branch fluxes (Equations 24, 25, and 26 in Materials and methods section "Dependence of steady state concentration on assimilated metabolite"). In the substituted enzyme mechanism, the lower bound on the concentration of $A$ is 0, at which the steady state concentration of $X$ is trivially 0 as well. In all cases, if the maximal flux of the autocatalytic reaction is higher than the maximal flux of the branch reaction, an upper bound on the concentration of $A$ may also exist, to satisfy the condition that $V}_{max,a}\text{}\text{}{V}_{max,b$. However, this bound can be removed by increasing ${V}_{\mathrm{max},b}$ or reducing ${V}_{\mathrm{max}}$.
Extensions of the simple autocatalytic cycle model
Generalizing for different autocatalytic stoichiometries
Our didactic analysis considered an autocatalytic reaction with 1:2 stoichiometry, such that for every substrate molecule consumed, two are produced. Realworld autocatalytic cycles may have different stoichiometries. For example, the CBB cycle has a stoichiometry of 5:6 so that for every 5 molecules of fivecarbon sugar that the autocatalytic reaction consumes, 6 fivecarbon molecules are produced. We can generalize our analysis by defining a positive $\delta $ such that for every molecule of $X$ that ${f}_{a}$ consumes, it produces $1+\delta $ molecules of $X$, where $\delta $ may be a fraction. This extension implies that Equation (1) becomes:
Therefore, all of the results above can be extended to different stoichiometries by replacing ${V}_{\mathrm{max},a}$ with $\delta \cdot {V}_{\mathrm{max},a}$. As a result, the qualitative conditions and observations from the 1:2 stoichiometry case remain valid but with a constant factor that changes the quantitative relations according to the relevant stoichiometry.
Input flux increases the range of parameters for which a stable steady state solution exists
Autocatalytic cycles are embedded within a larger metabolic network. Specifically, such cycles may have independent input fluxes to some of their intermediate metabolites, not requiring the use of other intermediate metabolites of the cycle. For example, in the glucose based, PTSdependent autocatalytic cycle, the existence of alternative transporters can generate flux of glucose 6phosphate into the cycle without the use of pep (Ferenci, 1996).
When adding a constant input flux, ${f}_{i}$ to our simple system (Figure 4A) the steady state condition changes to include this flux, giving:
In this situation, at $X=0$, $\dot{X}={f}_{i}\text{}\text{}0$ so the concentration of $X$ increases and there is no steady state at zero. If $V}_{max,b}\text{}\text{}{f}_{i}+{V}_{max,a$ then at a large enough value of $X$, $\dot{X}$ will be negative, implying that at some value of $X$ between these two extremes, $\dot{X}$ attains the value of zero, such that under this condition a positive stable steady state concentration exists (Figure 4I). This case therefore differs from the case with no input flux analyzed above, as now a positive stable steady state can always be achieved by modifying only ${V}_{\mathrm{max},a}$ and/or ${V}_{\mathrm{max},b}$. In this setup, cells can therefore tune the expression levels of enzymes to meet the needs of a stable steady state flux.
In cases where $V}_{max,b}\text{}\text{}{f}_{i}+{V}_{max,a$ either no steady states exist (Figure 4II), or two positive steady states exist (Figure 4III). The latter case implies that there exists a positive concentration $X$ that satisfies:
In this case, the lower concentration steady state will be stable.
To conclude, input fluxes change the steady state(s) of autocatalytic cycles. When an input flux is present, an autocatalytic cycle can always achieve a non zero, stable steady state by tuning the expression levels of the enzymes forming the cycle.
Interestingly, we find that in the two autocatalytic cycles shown in Figure 2—figure supplements 1 and 2, reactions that generate direct input flux into the cycle exist. In the ribose5P assimilating autocatalytic cycle (Figure 2—figure supplement 1), the rpi reaction serves as a shortcut, allowing input flux directly from ribose5P into the cycle. In the glyceronephosphate assimilating cycle (Figure 2—figure supplement 2), the tpi reaction similarly serves as such a shortcut. Interestingly, in these two cases, these shortcuts relax the constraints imposed by the strict use of the corresponding autocatalytic cycles as they prevent zero from being a stable steady state concentration. Another example of the effects of the addition of an input flux to an autocatalytic cycle is the input flux of fructose6phosphate from the catabolism of starch into the CBB cycle. This input flux can be used to ’kick start’ the cycle even without using the intermediate metabolites of the cycle.
Reversible branch reaction can either be far from equilibrium, resulting in the simple case, or near equilibrium, pushing the stability conditions down the branch pathway
The simple model assumed both the autocatalytic and the branch reactions are irreversible under physiological conditions. Assuming the branch reaction, ${f}_{b}$, can be reversible, with a product $Y$, the system can be analyzed in two extreme cases.
If $Y$ is consumed very rapidly by subsequent reactions, keeping its concentration low, then ${f}_{b}$ operates far from equilibrium. In this case, the reversible reaction equation reduces to an irreversible MichaelisMenten equation, resulting in the same constraints as in the simple, irreversible case analyzed above.
If $Y$ is consumed very slowly, and if the maximal consumption of $Y$ is larger than ${V}_{\mathrm{max},a}$, then, as long as $V}_{max,b}\text{}\text{}{V}_{max,a$, a stable steady state exists both when ${V}_{\mathrm{max},b}\to \mathrm{\infty}$, making ${f}_{b}$ operate near equilibrium, and when ${V}_{\mathrm{max},b}\to {V}_{\mathrm{max},a}$. A mathematical analysis is provided in the Materials and methods section "Reversible branch reaction analysis". The assumptions on the consumption of $Y$ in this case are analogous to the constraints in Equation 5, namely that the reaction downstream of $Y$ is less saturated than the autocatalytic reaction, and that it consumes $Y$ at a lower rate than the rate at which the autocatalytic reaction produces $X$ near $X=0$.
Analysis of a reversible autocatalytic cycle reaction
The simple model assumed both the autocatalytic and the branch reactions are irreversible under physiological conditions. The autocatalytic reaction in the simple model represents an effective overall reaction for all of the steps in autocatalytic cycles found in real metabolic networks. In order for the combined autocatalytic reaction to be physiologically reversible, all of the reactions in the real metabolic network must be reversible under physiological conditions. We note that this is not the case in any of the cycles we identify in central carbon metabolism. Nevertheless, this case can be mathematically analyzed.
If the autocatalytic reaction is reversible, then it must be driven by the displacement from thermodynamic equilibrium of the concentration of $A$ versus the concentration of $X$. Therefore, for any fixed concentration of $A$, $\widehat{A}$, a concentration of $X$ exists such that ${f}_{a}(X,\hat{A})=0\text{}\text{}{f}_{b}(X)$. It then follows that a sufficient condition for a positive steady state to exist is that at $X=0$, $\dot{X}(\hat{A})\text{}\text{}0$, which implies that
This condition can always be satisfied by high expression of the autocatalytic enzyme, increasing ${V}_{\mathrm{max},a}$. For this case, it therefore follows that for any concentration of $A$, a minimal value for ${V}_{\mathrm{max},a}$ exists, above which a positive steady state is achieved.
Stability analysis for multiplereaction cycles
Even the most compact realworld autocatalytic cycles are composed of several reactions. It is thus useful to extend the simple criteria we derived to more complex autocatalytic cycles. In such cycles the criteria for the existence of a steady state become nuanced and detail specific. We therefore focus on evaluating stability of such cycles, under the assumption that a nonzero steady state exists, which is usually known based on experimental measurements.
We analyze the stability criteria for the autocatalytic cycles depicted in Figure 5A and B in the Materials and methods, section "Extending the stability analysis from single to multiple reaction cycles". The analysis is performed for autocatalytic ratios up to 1:2, which is the case for all the autocatalytic cycles we identify. We find that in the multiple reaction case a steady state is stable if there exists $i$ such that $\beta}_{i}\text{}\text{}{\alpha}_{i$ (where $i$ can be any number in the range $1\mathrm{\dots}n$, ${\alpha}_{i}=\frac{\partial {f}_{{a}_{i}}}{\partial {X}_{i}}$, and ${\beta}_{i}=\frac{\partial {f}_{{b}_{i}}}{\partial {X}_{i}}$), Materials and methods section "Limits on derivatives of branch reactions for complex autocatalytic cycles".
Using the connection between derivatives of reactions and saturation levels (Materials and methods section "Inverse relationship between derivatives, affinities, and saturation levels"), we conclude that if ${\beta}_{i}>{\alpha}_{i}$ for some $i$ at the steady state point, then the saturation of the branch reaction, denoted $S(f)$, must be lower than the saturation of the corresponding cycle reaction at ${X}_{i}$:
This condition also dictates that the affinity of the branch reaction to the intermediate metabolite of the cycle it consumes must be weaker than the affinity of the corresponding recycling reaction of the cycle.
While having a single branch point at which $\beta}_{i}\text{}\text{}{\alpha}_{i$ is a sufficient condition for stability, we note that the larger the number of branch points satisfying this condition, the more robust the steady state point will be to perturbations, as such branch points reduce the propagation of deviations along the cycle (see Materials and methods section "Multiple unsaturated branch reactions increase convergence speed and dampen oscillations"). As we show below, these predictions hold for functioning autocatalytic cycles.
Using different kinetic equations
Although we utilized the widelyused irreversible MichaelisMenten kinetics equation to model enzyme kinetics, our results can be extended to different kinetic equations. Generally, two conditions must be met for a stable flux through an autocatalytic cycle to exist: (A) there should be a positive concentration of the intermediate metabolites for which the outgoing fluxes balance the autocatalytic fluxes, resulting in a steady state, and, (B) at the steady state point at least one derivative of an outgoing reaction out of the cycle should be higher than the derivative of the corresponding cycle reaction, as is implied by Equation 31, to enforce stability in the presence of small perturbations. Therefore, these two conditions should be explicitly evaluated for every case with different kinetic equations and autocatalytic cycles topologies to assert whether it can carry stable fluxes or not.
Testing the predictions of the analysis with experimental data on functioning autocatalytic cycles
To evaluate the validity of our analysis of autocatalytic cycles we searched for growth conditions under which the autocatalytic cycles we identified in central carbon metabolism carry substantial flux invivo. We used recent invivo flux measurements in E. coli from Gerosa et al. (2015). According to the data, two autocatalytic cycles carry substantial flux under at least one of the growth conditions measured: a cycle using the PTS carries significant fluxes in growth on glucose and on fructose; the glyoxylate cycle carries significant flux in growth on acetate and on galactose.
As noted above, we predict a design principle for functioning autocatalytic cycles: that at least one branch reaction should have a steeper response than the corresponding autocatalytic reaction at steady state. This requirement is sufficient, but not necessary, for the autocatalytic cycle to be at a stable steady state point. Moreover, having more than one branch point at which the branch reaction has a steeper response than the autocatalytic reaction increases the robustness of the steady state flux in the cycle as shown in the Materials and methods section "Multiple unsaturated branch reactions increase convergence speed and dampen oscillations". An outcome of the relationship between the steepness’s of responses is a reverse relationship between the saturation levels of the corresponding reactions (Equation 35). In order to evaluate the saturation level of a reaction under a given condition, two values must be obtained:
The maximal flux capacity of the reaction under the given condition, ${V}_{\mathrm{max}}$.
The actual flux through the reaction at the steady state, $f$.
To estimate the maximal capacity of a reaction we followed the procedure described in (Davidi et al., 201622) (see Materials and methods section "Evaluating maximal flux capacity of reactions under a given condition"). We used the data from (Gerosa et al., 2015) to identify the major branch points in each functioning cycle and the actual flux in them under each of the relevant conditions. The results are presented in Figure 6 and are provided, with the relevant calculations, in Supplementary file 1.
Our results show that for any of the four functioning autocatalytic cycle cases, in at least one branch point the biomass generating branch reaction has a larger maximal flux capacity, and is considerably less saturated than the respective autocatalytic reaction, in accordance with our predictions. Moreover, out of nine branch points analyzed, in six branch points the branching reactions were significantly less saturated than the autocatalytic reactions, in two branch points the saturation levels were similar, and only in one branch point the autocatalytic reaction was less saturated than the branching reaction.
The branch point at which the autocatalytic reaction is less saturated than the branch reaction is the branch point from fructose1,6bisphosphate in growth on fructose as the carbon source. The high saturation of the branch reaction arises as a large flux is reported for the fbp reaction, whereas the corresponding enzyme is not highly expressed under this condition. The large reported flux through fbp arises due to assuming a single transport pathway for fructose, entering the cycle as fructose1,6bisphosphate. However, an alternative fructose transport pathway is known to occur for the concentration at which the measurements were made (Kornberg, 1990). The alternative transport pathway produces fructose6phosphate from external fructose. Therefore, any flux through the alternative transport pathway should be directly deducted from the flux through fbp. Assuming $20\%$ of the consumed fructose uses this pathway suffices in order to balance the saturation levels at the fructose1,6bisphosphate branch point.
We made two negative control analyses to examine whether other reasons do not underlie the trend we find. First, we compared the saturation levels at the same branch points in growth conditions at which the autocatalytic cycles do not function, but the reactions carry flux. We find that for these cases, only 4 out of 9 cases satisfy the low branch saturation condition (Supplementary file 1). Second, we searched for branch points out of nonautocatalytic cycles and tested whether in such points branch reactions are also consistently less saturated than their corresponding cycle reactions. We found two fluxcarrying cycles: the TCA cycle, carrying flux in glucose, fructose, and glycerol growth, and a cycle consisting of the pentosephosphate pathway combined with gluconeogenesis, carrying flux in acetate, glycerol, and succinate growth. Out of the total six conditionsbranch points cases, in three the branch reaction was less saturated than the cycle reaction, and in three the cycle reaction was less saturated than the branch reaction (Supplementary file 1). We therefore conclude that, for cases that do not involve autocatalysis, the saturation of branch versus cycle reactions seems evenly distributed.
The consistently lower saturation values of biomass generating branch reactions demonstrate that the expressed enzymes have enough capacity to prevent the autocatalytic cycle from increasing the concentration of intermediate metabolites infinitely. Moreover, the lower saturation values of the biomass generating reactions suggest that at the steady state point their derivatives are higher, ensuring stable operation of the cycle.
Another demonstration of the autocatalytic mechanism being at play is in the CBB cycle, which is not a part of the metabolic network of wild type E. coli, and for which no flux measurements are available. This cycle has been recently introduced synthetically into E. coli and was shown to carry flux in it, given further metabolic engineering of central carbon metabolism (Antonovsky et al., 2016). The experimentally observed key evolutionary event enabling the functioning of the CBB cycle, was a mutation affecting the kinetic properties of the main branching reaction out of the CBB pathway, prs, weakening its affinity to its substrate, ribose5p. The observed weakening of affinity of prs is directly in line with our predictions on the relationship between the affinity of branch reactions and the affinity of the corresponding cycle reactions (see Materials and methods section of Antonovsky et al., 2016).
The other examples of autocatalytic cycles we found did not carry flux in any of the conditions for which data were available. The pentosephosphate cycle variants do not carry flux in any of the measured conditions, which is expected given that growth on ribose was not measured. The gluconeogenic FBA with ED pathway cycle also did not carry flux in any of the measured conditions. Although glycerol could have been a potential carbon source to use this pathway, the metabolic network allows for a more energy efficient growth by using the tpi reaction, as was indeed observed.
To conclude, existing data supports predictions made by our model, given the requirement for stable steady state operation of autocatalytic cycles. This agreement between predictions and measurements is especially encouraging given the highly limited information on kinetic properties, concentrations, and fluxes under various growth conditions.
Analysis of allosteric regulation potential for cycle improvement
Allosteric regulation can modulate the kinetic properties of enzymes at branch points, and of the cycle in general. As such, the relevant condition for the existence of a stable positive steady state should hold for the updated kinetic properties as defined following the effect of allosteric regulation.
We further analyze the ability of specific allosteric interactions to support fast convergence and stability of autocatalytic cycles in the Materials and methods section "Allosteric regulation can improve network performance". We compare the expected beneficial allosteric interactions against the allosteric regulation network of the two functioning autocatalytic cycles we identified, the PTSusing autocatalytic cycle and the glyoxylate cycle (Supplementary file 1, regulation data were taken from Keseler et al. (2013) and Schomburg et al. (2004)).
For the PTS cycle, we find that there are a total of 12 allosteric interactions, 7 inhibitions and five activations. Out of these 12 interactions, 11 interactions follow our expectations in terms of the type of the regulating metabolite (assimilated metabolite, cycle intermediate, or branch product), the regulated reaction (cycle reaction, branch reaction, or the reverse of a branch reaction), and the direction of the regulation (activation or inhibition). One interaction, the activation of fba by pep, does not follow our expectation.
For the glyoxylate cycle, we find that there are a total of 13 interactions, 12 inhibitions and one activation. Out of these 13 interactions, 8 interactions follow our expectations and 5 do not. The lack of significant agreement between the expected regulation direction and the actual regulation found for this cycle is consistent with the observation in Gerosa et al. (2015) that TCA cycle fluxes are regulated mainly by transcription and not by reactants levels.
It is important to note that allosteric regulation serves many roles, and that the metabolic network faces many more challenges than just the support of stable autocatalysis. Therefore, the agreement we find between existing allosteric interactions and the expected regulation scheme supporting autocatalysis does not suggest that the autocatalytic nature of the PTS is the only, or even main underlying reason for these allosteric interactions.
Discussion
Our study into the dynamics and stability of autocatalytic cycles suggests design principles applicable to both systems biology, that aims to understand the function of natural networks, and in the context of synthetic biology, in the effort to express novel heterologous cycles.
While autocatalytic cycles are often overlooked in the study of metabolism, we find that such cycles are at the heart of central carbon metabolism. Our autocatalytic modeling framework gives concrete predictions on saturation levels of branch reactions for operating autocatalytic cycles. We find these predictions agree well with empirically measured fluxomics and proteomics data sets. Given that there are other suggestions (Staples and Suarez, 1997; Weiss et al., 1998; Suarez et al., 1997) that may underlie the low saturation of branch reactions, we compare the saturation levels of branch reactions versus their corresponding cycle reactions both under conditions when the autocatalytic cycle does not function, and for branch points out of nonautocatalytic cycles. Both tests show no bias towards low saturation of branch reactions out of nonautocatalytic cycles, contrary to the clear trend we find for reactions branching out of autocatalytic cycles. Our findings thus support the addition of stability of intermediate metabolites of autocatalytic cycles as an explanation for the seemingly wasteful expression of enzymes (Salvador and Savageau, 2003, 2006). The model we present can also highlight metabolic branch points at which the kinetic efficiency of enzymes is constrained due to stability requirements of a corresponding autocatalytic cycle.
A common concept in synthetic biology is that the successful implementation of novel pathways requires the expression of functional enzymes in the correct amounts in the target organism. Here we show that in the context of autocatalytic cycles, such expression modulation may not suffice. Specifically, changes to the substrate affinity of enzymes at branch points of the cycle may be required in order for the novel pathway to function.
Another aspect of our findings is that while it is common to assume that strong affinity and high catalytic rate are desirable traits for enzymes, such seeming improvements may actually lead to instability and thus to non functional metabolic cycles. Furthermore, for reactions branching out of autocatalytic cycles, weaker affinities increase the steady state concentration of intermediate metabolites, resulting in higher fluxes both through the cycle, and through the branch reaction, suggesting an unconventional strategy for optimizing fluxes through such reactions. We note that because allosteric regulators modify the affinity of the enzymes they target, such regulators can potentially be used to restrict the affinity of branch reactions only when the autocatalytic cycle functions.
An experimental demonstration of these principles invivo is the recent implementation of a functional CBB cycle in E. coli by introducing the two genes missing for its function (Antonovsky et al., 2016). The successful introduction of the genes did not suffice to make the cycle function, and further directed evolution was needed in order to achieve successful operation of the cycle. Strikingly, most evolutionary changes occurred in branch points from the cycle (Antonovsky et al., 2016). The change which was biochemically characterized in the evolutionary process was the decrease of the value of $\frac{{k}_{\text{cat}}}{{K}_{M}}$ of phosphoribosylpyrophosphate synthetase (prs), one of the enzymes responsible for flux out of the CBB cycle, corresponding to the branch reaction in our simple model. This is beautifully in line with the predictions of our analysis that suggest that decreasing $\frac{{V}_{max}}{{K}_{M}}$ of branch reactions can lead to the existence of a stable flux solution.
Our observation regarding the stabilizing effect of input fluxes into an autocatalytic cycle can provide means to mitigate the stability issue in synthetic biology metabolic engineering setups. In such setups, conducting directed evolution under gradually decreasing input fluxes, such as those achieved in a chemostat, allows for a pathway to gradually evolve towards sustainable, substantial flux.
Finally, while our work focuses on cycles increasing the amount of carbon in the system, we note that autocatalysis can be defined with respect to other quantities such as energy (e.g. ATP investment and production in glycolysis [Teusink et al., 1998]), noncarbon atoms, reducing power, or other moieties (Reich and Selkov, 1981). As autocatalysis is often studied with relation to the origin of life, our analysis may be useful in studying synthetic autocatalytic systems such as the one recently described in Semenov et al. (2016). The analysis we present here can thus be of relevance for the analysis of metabolic networks in existing organisms and for the design of novel synthetic systems.
Materials and methods
Formal definition of an autocatalytic metabolic cycle
Given a metabolic network composed of a set of reactions and metabolites, the following criteria can be used to define a subset of the network that is an autocatalytic cycle: First we define a metabolic cycle. A set of irreversible reactions (for reversible reactions only one direction can be included in the set) and metabolites forms a cycle if every metabolite of the set can be converted, by sequential application of reactions in the set (where two reactions can be chained if a metabolite in the set is a product of the first reaction and a substrate of the second reaction), to every other metabolite in the set. A cycle is autocatalytic if the reactions of the cycle can be applied, each reaction at an appropriate, positive number of times, such that the resulting change in the amount of each of the metabolites forming the cycle is nonnegative, with at least one metabolite of the cycle having a strictly positive change.
The same definition can be stated in terms of reaction vectors and a stoichiometric matrix. If a metabolic network has $n$ metabolites, indicated by the numbers $1$ to $n$, then every reaction, $r$, in the network can be described as a vector ${V}_{r}$ in ℤ^{n}, such that the $i$’th coordinate of ${V}_{r}$ specifies how much of metabolite $i$ the reaction $r$ produces (if $r$ consumes a metabolite, then the value at the coordinate representing the metabolite is negative).
With this notation, a set of metabolites: $M={m}_{1}\mathrm{\cdots}{m}_{j}$ and a set of reactions, $R={r}_{1}\mathrm{\cdots}{r}_{k}$ form an autocatalytic cycle if:
Every row and every column of the stoichiometric matrix have at least one positive and one negative number.
There is a set of positive integers, ${i}_{1}\mathrm{\cdots}{i}_{k}$ such that the total reaction vector ${r}^{*}={\sum}_{l=1}^{k}{i}_{l}{r}_{l}$ is non negative at all the coordinates ${m}_{1}\mathrm{\cdots}{m}_{j}$ and is strictly positive for at least one coordinate in this range.
The condition that the total reaction vector ${r}^{*}={\sum}_{l=1}^{k}{i}_{l}{r}_{l}$ is non negative at all the coordinates ${m}_{1}\mathrm{\cdots}{m}_{j}$ and is strictly positive for at least one coordinate in this range cannot be satisfied by a set of nonnegative integers, ${i}_{1}\mathrm{\cdots}{i}_{k}$, if this set includes values that are 0 (this condition eliminates the addition of disjoint cycles to an autocatalytic cycle).
Systematic identification of autocatalytic cycles in metabolic networks
We implemented an algorithm to systematically search for autocatalytic cycles in metabolic networks. The algorithm is not comprehensive, in the sense that there may be autocatalytic cycles that will not be identified by it. Further work will enable a more advanced algorithm to identify additional autocatalytic cycles in full metabolic networks. We used the algorithm on the core carbon metabolism network of E. coli (Orth et al., 2010).
In our framework, a metabolic network is defined by a set of reactions, $(\overline{R})$. Each reaction is defined by a set of substrates and a set of products, with corresponding stoichiometries ${R}_{i}=(S,P,{N}^{S},{N}^{P})$, such that ${R}_{i}$ describes the reaction ${\sum}_{j}{N}_{j}^{S}{S}_{j}\to {\sum}_{k}{N}_{k}^{P}{P}_{k}$. The algorithm works as follows:
All cofactors are removed from the description of the metabolic network.
The metabolic network is converted to a directed graph, G: The nodes of G are all the metabolites and all the reactions of the network. For each metabolite, $M$, and each reaction, $R$, if $M$ is a substrate of $R$ then the edge $(M,R)$ is added to the graph, and if $M$ is a product of $R$, then the edge $(R,M)$ is added to the graph.
The Tarjan cycle identification algorithm is used to enumerate all the cycles in the graph (Tarjan, 1973).
For every cycle identified by the Tarjan algorithm, $C$, the algorithm checks if the cycle can be the backbone of an autocatalytic cycle as follows:
For every reaction in the cycle, $R$, the algorithm checks if it consumes more than one intermediate metabolite of the cycle. If so, $C$ is assumed not to be autocatalytic and the algorithm continues to evaluate the next cycle.
Otherwise, for every reaction in the cycle, $R$, the algorithm checks if it has more than one product that is an intermediate metabolite of the cycle. If so, then the algorithm lists $C$ as an autocatalytic cycle.
Finally, the algorithm checks, for every reaction in the cycle, if it has a product that is not an intermediate metabolite of the cycle. If so, denote by ${M}_{E}$ such a metabolite. The algorithm proceeds to check if, for every intermediate metabolite of the cycle, ${M}_{i}$ a reaction exists from ${M}_{E}$ to ${M}_{i}$ that does not use any of the reactions of the cycle, and does not consume any of the intermediate metabolites of the cycle. If so then the algorithm lists $C$ as an autocatalytic cycle.
The algorithm assumes reactions consume exactly one molecule of any of their substrates and produce exactly one molecule of any of their products, an assumption that holds for the core model of E.coli, but not in metabolic networks in general.
Steady state concentration dependence on kinetic parameters of autocatalytic and branch reactions
The simple cycle steady state concentration, ${X}^{*}$, is given in Equation 2. Taking the derivative of this expression with respect to ${K}_{M,a}$, ${K}_{M,b}$, ${V}_{\mathrm{max},a}$, and ${V}_{\mathrm{max},b}$, under the assumption that the kinetic parameters satisfy the stability conditions in Equation 5 gives:
So that ${X}^{*}$ increases when ${K}_{M,a}$ decreases or ${V}_{\mathrm{max},a}$ increases (or both) corresponding to activation of ${f}_{a}$. On the other hand, ${X}^{*}$ decreases when ${K}_{M,b}$ decreases or ${V}_{\mathrm{max},b}$ increases (or both) corresponding to activation of ${f}_{b}$.
Connecting bisubstrate reaction kinetic constants with simple MichaelisMenten constants
Three standard equations are commonly used to describe the flux through irreversible bisubstrate reactions (Leskovac, 2003). We show that, under the assumption that the assimilated metabolite maintains constant concentration, these equations reduce to simple MichaelisMenten equations. We derive the expressions for the apparent MichaelisMenten constants, ${K}_{M}$ and ${V}_{\mathrm{max}}$, as functions of the kinetic constants of the bisubstrate reaction and the concentration of the assimilated metabolite. While the substrates in these equations are generally denoted as $A$ and $B$, here, to keep the notation consistent, we will denote by $A$ the assimilated metabolite and by $X$ the internal metabolite of the cycle.
The simplest equation describing a bisubstrate reaction assumes substituted enzyme (Ping Pong) mechanism (Imperial and Centelles, 2014). As this equation is symmetric with respect to the two substrates, we can arbitrarily decide which of the two substrates is the assimilated metabolite, and which is the internal metabolite. We get that the flux through the reaction is:
Rearranging to get the dependence of the flux on $X$ in a MichaelisMenten like form we get that:
which gives apparent MichaelisMenten kinetic constants of:
The second bisubstrate reaction mechanism we consider is the ternary enzyme complex with random binding order of the two substrates. As this equation is also symmetric with respect to the two substrates, we can again arbitrarily decide which of the two substrates is the assimilated metabolite, and which is the internal metabolite. We get that the flux through the reaction is:
Rearranging to get the dependence of the flux on $X$ in a MichaelisMenten like form we get that:
which gives apparent MichaelisMenten kinetic constants of:
The other equation describing a ternary enzyme complex bisubstrate reaction assumes ordered binding of the substrates. Because in ordered binding the equation is asymmetric with respect to the two substrates, analyzing this reaction is further split according to which of the two substrates is assumed to be the assimilated metabolite with constant concentration.
If the first binding metabolite is assumed to be the assimilated metabolite we get that:
which gives apparent MichaelisMenten kinetic constants of:
If the first binding metabolite is assumed to be the internal metabolite we get that:
which gives apparent MichaelisMenten kinetic constants of:
To summarize, the most common equations describing bisubstrate reactions reduce to equations of the same form as MichaelisMenten equations, under the assumption that one of the metabolites maintains a constant concentration. The apparent kinetic constants of the MichaelisMenten equivalent equations depend on the kinetic constants of the bisubstrate reactions, as well as on the concentration of the assimilated metabolite.
Constraints on concentration of assimilated metabolite and kinetic constants of bisubstrate reactions
In Equation 5 we obtain constraints on the kinetic parameters of MichaelisMenten reactions that ensure the existence and stability of a positive steady state. We observe that these constraints imply that even if the maximal rates of the two reactions can be easily modified, if ${K}_{M,b}<{K}_{M,a}$ then such changes cannot suffice in order to satisfy the existence and stability constraints.
Here, we map the same constraints from Equation 5 onto bisubstrate autocatalytic reactions. This mapping results in constraints on the assimilated metabolite concentration, as well as on the kinetic parameters of the bisubstrate autocatalytic reactions. We show that in all ternary enzyme complex bisubstrate reaction schemes, there is a lower bound on the concentration of the assimilated metabolite, below which the system cannot attain a stable positive steady state. We further show that the nature of the bisubstrate reaction qualitatively affects the ability to satisfy the stability constraints by changing expression levels alone. In the cases of substituted enzyme mechanism, random binding order ternary complex, and ordered binding ternary complex, with the assimilated metabolite binding first, unless the kinetic parameters of the participating enzymes satisfy specific inequalities, changes to the maximal reaction rates alone cannot suffice in order to satisfy the existence and stability constraints. However, in the case of ordered binding ternary complex with the internal metabolite binding first, changes to the maximal reaction rates alone suffice in order to allow for stable steady state to occur, given high enough concentration of the assimilated metabolite. We analyze each of the four possible bisubstrate reaction schemes separately below.
Substituted enzyme (Ping Pong) mechanism
The case of substituted enzyme mechanism is the simplest case to analyze. We can substitute Equation 9 into the conditions from Equation 5 to get:
As $\frac{A}{{K}_{A}+A}\text{}\text{}1$, the first inequality can always be satisfied if $V}_{max,b}\text{}\text{}{V}_{max$, which is equivalent to the first condition in Equation 5. The second condition is identical to the second condition from Equation 5. Therefore, this case imposes equivalent conditions to those derived for the simple, single substrate case.
Random binding order
In the case of a random binding order, we can substitute Equation 11 into the conditions from Equation 5 to get:
We first note that from the second inequality we get that:
Giving a lower bound on the concentration of the assimilated metabolite for which a stable steady state is attainable.
We now wish to obtain a lower bound on ${K}_{M,b}$. In order to obtain such a lower bound, we need an upper bound on ${V}_{\mathrm{max}}\frac{A}{{K}_{i,A}+A}$. However, we only have an upper bound on ${V}_{\mathrm{max}}\frac{A}{{K}_{A}+A}$. We use the first inequality in Equation 17 to get that:
We note that for positive $A$, $\frac{{K}_{A}+A}{{K}_{i,A}+A}\text{}\text{}max(1,{K}_{A}/{K}_{i,A})$ and therefore:
Substituting this inequality in the second inequality of Equation 17 therefore gives us that:
We have therefore obtained a lower bound on the affinity of the branch reaction, ${K}_{M,b}$, in this case.
For the random binding order we can thus conclude that, like in the singlesubstrate case, a lower bound exists on the affinity of the branch reaction, below which a positive steady state is not attainable, even if the expression levels of the enzymes, and the concentration of the assimilated metabolite are modified. Furthermore, for any set of kinetic parameters, there is a lower bound on the concentration of $A$ for which a positive steady state is attainable.
Ordered binding with the assimilated metabolite binding first
In the case of ordered binding, with the assimilated metabolite binding first, we can substitute Equation 13 into the conditions from Equation 5 to get:
As the second inequality is identical to the one in the random binding order case, we can immediately conclude that the same lower bound on the concentration of $A$ from Equation 18 holds in this case as well.
Regarding a lower bound on ${K}_{M,b}$, following a similar reasoning as in the previous case, we first note that for any value of $A$:
Therefore, we can deduce, by using the first inequality from Equation 20 in the second inequality from that equation, that:
which immediately yields:
setting an absolute lower bound on ${K}_{M,b}$.
We thus arrive at the same conclusions in this case, as we have arrived to in the previous case, namely that a lower bound exists on the affinity of the branch reaction, and that, for any set of kinetic parameters, there is a lower bound on the concentration of the assimilated metabolite, below which a positive stable steady state cannot be obtained.
Ordered binding with the internal metabolite binding first
In the case of ordered binding, with the internal metabolite binding first, we can substitute Equation 15 into the conditions from Equation 5 to get:
As in the previous two cases, the second inequality can be used to obtain a lower bound on the concentration of $A$:
However, unlike in the previous two cases, in this case if $V}_{max}\text{}\text{}{V}_{max,b$, then the first inequality in Equation 22 holds for any concentration of $A$, and, for any concentration of $A$ that is larger than its lower bound, the second inequality is also satisfied, resulting in a stable steady state. This case is therefore more robust than the other cases as it allows for the conditions to be satisfied, at least for high concentrations of $A$, given any set of kinetic parameters.
Dependence of steady state concentration on assimilated metabolite
Equation 2 shows the dependency between the steady state concentration of the internal metabolite $X$, ${X}^{*}$, and the kinetic parameters of the reactions in the system. Substituting the dependencies of the apparent kinetic parameters from Equations 11, 13, and 15 into Equation 2 gives the dependency of ${X}^{*}$ on the kinetic parameters of the bisubstrate reactions and the concentration of the assimilated metabolite, $A$. We get for these three cases respectively that:
Assuming the kinetic parameters satisfy the stable steady state conditions derived in Equations 17, 20, and 22, we note that when $A$ is equal to its lower bound, the numerator in all three cases is 0, resulting in ${X}^{*}=0$. Furthermore, as $A$ decreases towards its lower bound, ${X}^{*}$ decreases resulting in a decrease in both ${f}_{b}$ and ${f}_{a}$ (for the two latter cases this is trivial to show, as the terms involving $A$ increase and decrease monotonically in accordance with their effect on ${X}^{*}$. In the first case, taking the derivative of the numerator w.r.t. $A$ shows the derivative is always positive, resulting in the same conclusion). Interestingly, in the first and last cases, if $V}_{max}\text{}\text{}{V}_{max,b$, then an upper bound on the concentration of $A$ also exists. As the concentration of $A$ approaches this upper bound, the denominator approaches 0 resulting in an increase in the concentration of ${X}^{*}$ towards infinity.
Reversible branch reaction analysis
The simple model assumed both the autocatalytic and the branch reactions are irreversible. Here we assume the branch reaction is reversible, and let $Y$ denote its product. For simplicity, we further assume that ${K}_{\text{eq}}=1$, noting that this assumption can always be satisfied by measuring the concentration of $Y$ in units of ${K}_{\text{eq}}X$. We recall that the reversible MichaelisMenten equation states that:
We assume that a third reaction, ${f}_{c}$, irreversibly consumes $Y$. While assuming ${f}_{c}$ follows irreversible MichaelisMenten kinetics is analytically tractable, the analysis is simpler, and as informative, under the assumption that ${f}_{c}=DY$ for some constant $D$. This simplification is equivalent to assuming ${f}_{c}$ follows MichaelisMenten kinetics with $\frac{{V}_{\mathrm{max},c}}{{K}_{M,c}}\approx D$, and ${V}_{max,c}\text{}\text{}max({V}_{max,a},{V}_{max,b})$.
We start by deriving the necessary conditions for steady state existence. Because at steady state ${f}_{a}={f}_{c}$, it follows that:
Furthermore, as at the steady state ${f}_{a}={f}_{b}$, we get that:
Substituting ${Y}^{*}$ from Equation 27 gives:
Which is satisfied when ${X}^{*}=0$ (implying that ${X}^{*}={Y}^{*}=0$ is a steady state), or when ${X}^{*}$ satisfied the quadratic equation:
Albeit intimidating, this quadratic equation can be used to derive the conditions for the existence of a positive steady state. Only if both of the roots of this equation are negative, no positive steady state exists. We recall that the two roots of a quadratic equation of the form $0=a{X}^{2}+bX+c$ are negative iff:
As in the irreversible case, the sign of ${V}_{\mathrm{max},b}{V}_{\mathrm{max},a}$ determines the required condition on the numerators. We assume that $V}_{max,b}\text{}\text{}{V}_{max,a$, noting that if $V}_{max,b}\text{}\text{}{V}_{max,a$, a positive steady state cannot be globally stable because for $X$ such that $f}_{a}(X)\text{}\text{}{V}_{max,b$, the system will diverge regardless of the value of $Y$.
Under the assumption that $V}_{max,b}\text{}\text{}{V}_{max,a$, the denominator of both $b$ and $c$ is positive, meaning a positive steady state exists only if the nominators of $b$ or $c$ (or both) are negative. Thus, two options may arise.
If ${K}_{M,a}\text{}\text{}{V}_{max,a}/D$ (implying $D\text{}\text{}{V}_{max,a}/{K}_{M,a}$, qualitatively suggesting rapid removal of $Y$) then an upper bound on ${V}_{\mathrm{max},b}$ exists, above which the two solutions are negative, implying no positive steady state exists. A sufficient condition for existence in this case is that $\frac{{V}_{max,b}}{{K}_{X}}\text{}\text{}\frac{{V}_{max,a}}{{K}_{M,a}}$, ensuring that $c\text{}\text{}0$. This condition is equivalent to the condition in the irreversible case. We further show below that for large enough $D$, the resulting steady state is stable.
If $D\text{}\text{}{V}_{max,a}/{K}_{M,a}$, then for any $V}_{max,b}\text{}\text{}{V}_{max,a$, $c\text{}\text{}0$ implying a positive steady state exists. As we show below, in this case both when ${V}_{\mathrm{max},b}\to {V}_{\mathrm{max},a}$, and when ${V}_{\mathrm{max},b}\to \mathrm{\infty}$, the steady state is stable.
We now turn to analyze the stability of the steady state. For a steady state to be stable, the eigenvalues of the Jacobian matrix must have negative real values. In our system it holds that
We use the following notation:
We can use this notation to write the Jacobian matrix as:
which gives a characteristic polynomial of:
In order for the real values of the roots of the characteristic polynomial to be negative it must hold that $b\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ and $c\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$, where $b$ and $c$ are now the coefficients of the quadratic equation $a{\lambda}^{2}+b\lambda +c=0$. We therefore get that:
We denote by ${f}^{*}$ the steady state flux in the system, such that ${f}^{*}={f}_{a}={f}_{b}={f}_{c}$ We note that for MM kinetics and positive concentrations it holds that:
First, we note that if $\alpha \ge {\beta}_{x}$ then the steady state cannot be stable as, looking at the value of $c$, we see that in such a case $({\beta}_{x}\alpha )D\text{}\text{}0$ and since $\alpha {\beta}_{y}\text{}\text{}0$, $c\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}0$ violating the stability conditions. However, because we assume that $V}_{max,b}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{V}_{max,a$, then for $Y={Y}^{*}$, at $X=0$, $f}_{b}\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{f}_{a$, but for $X\to \mathrm{\infty}$, ${f}_{b}\to {V}_{\mathrm{max},b}$ and ${f}_{a}\to {V}_{\mathrm{max},a}$, so that $f}_{b}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{f}_{a$. It then follows that, because the two fluxes can only intersect once for positive $X$ and fixed $Y$, at the steady state point, where ${f}_{a}={f}_{b}$, $\alpha \phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{\beta}_{x}$, so this condition is satisfied. We note that this condition is sufficient to ensure that $b\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$. We also note that as $\alpha \phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{\beta}_{x}$, a large enough value of $D$ exists at which the steady state is stable, concluding that if $D$ is large enough, then a stable steady state exists if:
If $D$ is small, such that $D\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{V}_{max,a}/{K}_{M,a}$, and $V}_{max,b}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{V}_{max,a$ (implying that $\alpha \phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{\beta}_{x}$), we need to check what other conditions are necessary in order to ensure that ${\beta}_{x}D+\alpha {\beta}_{y}\alpha D\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$. We look at the limit ${V}_{\mathrm{max},b}\to \mathrm{\infty}$. At this limit, the quadratic equation for ${X}^{*}$ converges to:
For this equation, $c\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}0$, implying that one of the roots is negative and one is positive. The positive root is:
As this ${X}^{*}$ is finite, we get that when ${V}_{\mathrm{max},b}\to \mathrm{\infty}$, ${Y}^{*}$ also converges to $\frac{{V}_{\mathrm{max},a}}{D}{K}_{M,a}$. At this limit, ${\beta}_{x}$ increases infinitely and ${\beta}_{y}$ decreases infinitely, but ${\beta}_{x}+{\beta}_{y}$ converges to:
that is constant. Therefore, rearranging c such that:
we note that as ${V}_{\mathrm{max},b}$ increases, the dominant term becomes ${\beta}_{y}(D\alpha )\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ ensuring that $c\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ and therefore stability.
On the other hand, when ${V}_{\mathrm{max},b}\to {V}_{\mathrm{max},a}$, we note that because $f}_{c}\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}{V}_{max,a$, ${Y}^{*}$ is bounded by $Y}^{\ast}\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}\frac{{V}_{max,a}}{D$, but ${X}^{*}\to \mathrm{\infty}$. Thus, both $\alpha $ and ${\beta}_{x}$ diminish like ${\frac{1}{{X}^{*}}}^{2}$, and ${\beta}_{y}$ diminishes like $\frac{1}{X}$. The dominant term in $c={\beta}_{x}D+{\beta}_{y}\alpha \alpha D$ therefore becomes $({\beta}_{x}\alpha )D\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ so again stability is maintained.
Therefore, for small values of $D$, as long as $V}_{max,b}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{V}_{max,a$, a positive stable steady state exists both in the lower limit of ${V}_{\mathrm{max},b}\to {V}_{\mathrm{max},a}$, and in the upper limit of ${V}_{\mathrm{max},b}\to \mathrm{\infty}$.
Our conclusions are therefore as follows: As in the irreversible case, $V}_{max,b}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{V}_{max,a$ is a necessary condition for the existence of a globally stable steady state. For large values of $D$, the reversible reaction is far from equilibrium, resulting in an additional condition, equivalent to the condition we obtained for the irreversible case, namely that ${V}_{\mathrm{max},b}/{K}_{X}$ is upper bounded by a term that is larger than ${V}_{\mathrm{max},a}/{K}_{M,a}$, but approaches it as $D$ increases. This condition is sufficient for the existence and stability of the steady state. For small values of $D$, a steady state always exists (given that $V}_{max,b}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{V}_{max,a$). We can show that this steady state is stable both when ${V}_{\mathrm{max},b}\to \mathrm{\infty}$, and when ${V}_{\mathrm{max},b}\to {V}_{\mathrm{max},a}$. We therefore conclude that in this case, no further restrictions apply on ${K}_{X}$, ${K}_{Y}$, or ${K}_{M,a}$ but rather that a steady state can always be achieved at most by changing ${V}_{\mathrm{max},b}$.
Qualitatively, the cases we analyze show that, on top of the required $V}_{max,b}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{V}_{max,a$ condition, the second condition is that either the slope of ${f}_{c}=D$ is smaller than ${V}_{\mathrm{max},a}/{K}_{M,a}$, or that the maximal slope of ${f}_{b}$, ${V}_{\mathrm{max},b}/{K}_{X}$, is smaller than ${V}_{\mathrm{max},a}/{K}_{M,a}$.
Extending the stability analysis from single to multiple reaction cycles
We analyze the stability criteria for the autocatalytic cycles depicted in Figure 5A and B. We start by writing the relevant equations for the autocatalytic cycle depicted in Figure 5A. In this system, there are two intermediate metabolites, ${X}_{1}$ and ${X}_{2}$, two reactions that form the cycle, ${f}_{{a}_{1}}$ and ${f}_{{a}_{2}}$, and two branch reactions, ${f}_{{b}_{1}}$ and ${f}_{{b}_{2}}$. We assume, without loss of generality, that the autocatalytic reaction (the reaction that increases the amount of carbon in the cycle) is ${f}_{{a}_{2}}$ and that the autocatalysis is in a 1:2 ratio. The equations describing the dynamics of the system are thus:
We note that in steady state, where $\stackrel{.}{X}1=\stackrel{.}{X}2=0$, because the autocatalysis is in a 1:2 ratio, it must hold that ${f}_{{b}_{1}}+{f}_{{b}_{2}}={f}_{{a}_{2}}$, meaning the total outgoing flux balances the total increase of intermediate metabolites due to autocatalysis. Given that a steady state of the system exists for some value $({X}_{1}^{*},{X}_{2}^{*})$, we can evaluate the condition for stability. In multivariable systems, stability dictates that the real part of the eigenvalues of the Jacobian matrix must all be negative. We define ${\alpha}_{i}=\frac{\partial {f}_{{a}_{i}}}{\partial {X}_{i}}$ and ${\beta}_{i}=\frac{\partial {f}_{{b}_{i}}}{\partial {X}_{i}}$ for $i=1,2$. We note that as we assume Michaelis Menten kinetics, ${\alpha}_{i}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$ and ${\beta}_{i}\phantom{\rule{thinmathspace}{0ex}}\ge \phantom{\rule{thinmathspace}{0ex}}0$, where ${\beta}_{i}=0$ is the case where there is no flux branching out at $i$. We then get that the Jacobian matrix is:
Solving for the characteristic polynomial gives:
that has two negative roots when:
which is satisfied if $\beta}_{1}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{\alpha}_{1$ or $\beta}_{2}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{\alpha}_{2$. Therefore, if either $\beta}_{1}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{\alpha}_{1$ or $\beta}_{2}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{\alpha}_{2$ at the steady state, then the steady state is stable.
The twometabolites cycle case can be easily extended to a larger number of intermediate metabolites and reactions, as is depicted in Figure 5B. For this extension, we again assume, without loss of generality, that the autocatalytic reaction is the last reaction, ${f}_{{a}_{n}}$, and that the autocatalysis is in a 1:2 ratio.
In this case, steady state implies that the concentration of each intermediate metabolite is conserved, meaning that for all $i\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1$:
(for $i=1$, as ${f}_{{a}_{n}}$ is the autocatalytic reaction, we get that $2\cdot {f}_{{a}_{n}}\phantom{\rule{thinmathspace}{0ex}}\ge \phantom{\rule{thinmathspace}{0ex}}{f}_{{a}_{1}}$). Also, because at steady state the total outgoing flux from the cycle must balance the total incoming flux into the system, which is the amount of autocatalysis carried out by ${f}_{{a}_{n}}$, we get that:
(due to our assumption of a 1:2 autocatalytic ratio) implying that for all $i$:
We stress that Equation 29 is only valid if the autocatalysis is in up to a 1:2 ratio. Deriving a stability criterion for the multiplereaction case, we get that in this case a steady state is stable if there exists $i$ such that $\beta}_{i}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{\alpha}_{i$ (see section 9 below).
To conclude, for the straightforward extension of the simple model to multiple reactions with a single autocatalytic reaction, steady state implies that for all $i$:
Where the left inequality is due to Equation 29 and the right inequality is due to Equation 28.
A sufficient condition for such a steady state point to be stable is that at the steady state point there exists at least one branching point $i$ at which the derivative of the branch reaction is larger than the derivative of the equivalent autocatalytic reaction:
Limits on derivatives of branch reactions for complex autocatalytic cycles
Stability analysis of a model complex autocatalytic cycle with $n$ reactions in the cycle results in the following Jacobian matrix:
The characteristic polynomial of this matrix is given by:
To extract the conditions under which all the roots of the characteristic polynomial have negative real parts we use Rouche’s theorem. Our strategy will be as follows: We will define a contour that contains only numbers with negative real parts. We will show that all the roots of the polynomial $0={\prod}_{i=1}^{n}(\lambda +{\alpha}_{i}+{\beta}_{i})$ lie within the area this contour encloses. We will find the conditions for which ${\prod}_{i=1}^{n}(\lambda +{\alpha}_{i}+{\beta}_{i})>2{\prod}_{i=1}^{n}{\alpha}_{i}$ on the contour, satisfying the premise of Rouche’s theorem. We will then claim that under these conditions all the roots of the polynomial in Equation 33 must also lie inside the contour, and therefore must have negative real parts. Given that all the roots of this polynomial have negative real parts, we will conclude that the eigenvalues of the Jacobian matrix at the steady state point all have negative real parts, making the steady state point stable.
Proof
We pick a large parameter $R$, such that $R\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}3\text{}{\mathrm{m}\mathrm{a}\mathrm{x}}_{j}({\alpha}_{j}+{\beta}_{j})$. We look at the closed half circle contour, $K$, composed of the segment $[(0,iR),(0,iR)]$ and the half circle arc $(x,iy)$ such that $x\phantom{\rule{thinmathspace}{0ex}}\le \phantom{\rule{thinmathspace}{0ex}}0$ and ${x}^{2}+{y}^{2}={R}^{2}$. We define
noting that it is constant over all of ℂ and specifically over $K$. We define
noting that all of $f$’s roots lie inside $K$ as the roots are $0\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}({\alpha}_{j}+{\beta}_{j})\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}R$ for all $j$. We check the conditions under which $f(\lambda )>g(\lambda )$ over the contour $K$.
For the arc segment we note that, as for complex numbers it holds that $xy=xy$, then
From the triangle inequality we know that $x+y\ge xy$ and therefore for all $j$ it holds that
As we picked $R$ such that $R>3\text{}{\mathrm{m}\mathrm{a}\mathrm{x}}_{j}({\alpha}_{j}+{\beta}_{j})$ we get that
and therefore
concluding that over the arc, $f(\lambda )\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}g(\lambda )$.
For the part of $K$ on the imaginary axis, we note that $\lambda =iy$ where $y\in [R,R]$. For this segment we therefore get that
and, as before, that
To meet the condition that $f(\lambda )\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}g(\lambda )$, which is equivalent to: $\frac{f(\lambda )}{g(\lambda )}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}1$, it is sufficient to find the conditions under which:
Simplifying this inequality gives:
A sufficient condition to satisfy this inequality, given that all the ${\alpha}_{j}$ are positive and all the ${\beta}_{j}$’s are non negative, is that there exists $j$ such that $\beta}_{j}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{\alpha}_{j$.
We therefore get that if there exists $j$ such that $\beta}_{j}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}{\alpha}_{j$, then $f(\lambda )\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}g(\lambda )$ over the contour $K$. In this case, by Rouche’s theorem, we deduce that, as all of $f$’s roots lie inside $K$, then it follows that all of $fg$’s roots lie inside $K$, concluding that the real part of all of the eigenvalues of the characteristic polynomial of the Jacobian matrix of the complex autocatalytic cycle are negative, making any steady state that meets this criterion stable.
Multiple unsaturated branch reactions increase convergence speed and dampen oscillations
Using the Jacobian matrix from Equation 32 we can analyze the effect of multiple low saturation branch points on convergence to steady state. The analysis shows that the more $i$’s exist for which ${\beta}_{i}\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}0$, and the larger ${\beta}_{i}$ is (resulting in lower saturation of ${f}_{{b}_{i}}$), the faster the convergence of the cycle to steady state will be.
We denote by $\overrightarrow{{X}^{*}}$ the steady state vector of the concentrations of the intermediate metabolites. We denote by $\overrightarrow{X}=\overrightarrow{{X}^{*}}+\mathrm{\Delta}{X}_{j}$ a state where for all the intermediate metabolites that are not ${X}_{j}$, their concentration is the same as the steady state concentration, and ${X}_{j}$ differs by a small amount, $\mathrm{\Delta}{X}_{j}$, from its steady state concentration. We let $F$ denote the fluxes function of the system such that $F(X)={\dot{X}}_{\overrightarrow{X}}$. Evaluating the dynamics of the system at $\overrightarrow{X}$ by noting that $F(\overrightarrow{X})\approx F(\overrightarrow{{X}^{*}})+J\cdot \mathrm{\Delta}{X}_{j}=J\cdot \mathrm{\Delta}{X}_{j}$ (where $F(\overrightarrow{{X}^{*}})=0$ as $\overrightarrow{{X}^{*}}$ is a steady state) results in $F{(\overrightarrow{X})}_{k}=0$ for all $k\ne j,j+1$. For ${X}_{j}$ such that $j\ne n$ we get $F{(\overrightarrow{X})}_{j}\approx ({\alpha}_{j}+{\beta}_{j})\mathrm{\Delta}{X}_{j}$ and for ${X}_{j+1}$ we get $F{(\overrightarrow{X})}_{j+1}\approx {\alpha}_{j}\mathrm{\Delta}{X}_{j}$. Therefore, the difference from the steady state decreases proportionally to ${\beta}_{j}$ (and cycles to the next intermediate metabolite, ${X}_{j+1}$). For $j=n$, we get that $F{(\overrightarrow{X})}_{j}\approx ({\alpha}_{j}+{\beta}_{j})\mathrm{\Delta}{X}_{j}$, as for $j\ne n$, but $F{(\overrightarrow{X})}_{1}\approx 2{\alpha}_{j}\mathrm{\Delta}{X}_{j}$ where the factor of $2$ is due to the effect of the assimilating reaction, that causes an amplification of the deviation from steady state (an amplification that is dampened by subsequent reactions along the cycle if the conditions for stable steady state are satisfied).
It therefore follows that any increase in ${\beta}_{j}$, for any $j$, increases the speed of convergence to steady state and reduces the propagation of deviations from steady state for ${X}_{j}$. Because of the linearity of matrix multiplication, an arbitrary deviation from $\overrightarrow{{X}^{*}}$ can always be decomposed to individual deviations with respect to every intermediate metabolite, making the analysis above valid for such deviations as well. Thus, to keep deviations from steady state at check, it is beneficial to increase ${\beta}_{j}$, for all $j$, which implies decreasing the saturation of ${f}_{{b}_{j}}$.
Inverse relationship between derivatives, affinities, and saturation levels
It turns out that for the MichaelisMenten kinetics equations, the following useful lemma can be used to connect theoretical observations on the relationships of derivatives to physiological observations on affinities and saturation levels.
We define the saturation level of a reaction as the ratio between the flux it carries, and the maximal flux it can carry, given the expression level of the relevant enzyme, that is:
Given this definition we can show that if two MichaelisMenten reactions consume the same metabolite, $X$, and at a given concentration, ${X}^{*}$, it holds that ${f}_{a}({X}^{*})\ge {f}_{b}({X}^{*})$, then if:
then it follows that:
Proof: expanding the condition that ${f}_{a}({X}^{*})\ge {f}_{b}({X}^{*})$, we get that:
Expanding the premise of the lemma in Equation 34 gives us that:
Because Equation 36 holds, it follows that:
setting the affinity of the autocatalytic enzyme as a lower bound for the affinity of the branch enzyme. Finally, given this relation of affinities it follows that:
concluding the proof.
We note that a multiple reaction autocatalytic cycle at a stable steady state point satisfies Equations 30 and 31, so the lemma applies.
Evaluating maximal flux capacity of reactions under a given condition
To evaluate the maximal flux capacity of a reaction under a prescribed growth condition, given expression level and flux data for a set of conditions, we follow the procedure described in Davidi et al. (2016). For each reaction, under every condition, we divide the flux the reaction carries (obtained from Gerosa et al., 2015) by the amount of the corresponding enzyme expressed under that condition (obtained from Schmidt et al., 2016). We thus get a flux per enzyme estimate for the given reaction under each of the conditions. We define the enzyme maximal invivo catalytic rate as the maximum flux per unit enzyme it carries across all conditions analyzed (noting that this is actually only a lower bound on this rate). Multiplying the enzyme maximal catalytic rate by the enzyme amount at each condition results in an estimate of the maximal possible flux through the given reaction under the relevant condition.
Allosteric regulation can improve network performance
In this section we touch upon the potential (in somewhat simplified and naively nonrigorous terms) of allosteric regulation to improve the properties of autocatalytic cycles. The constraint on the affinity of the branch reaction imposed by the stability requirement (Equation 35) may be suboptimal under other flux modes. Furthermore, allosteric regulation can be used to accelerate the rate at which an autocatalytic cycle converges to its stable steady state mode. While many allosteric regulation schemes exist (Leskovac, 2003), all of these schemes affect the affinity of the regulated enzyme, and some of these schemes also affect the maximal rate. We qualitatively analyze the expected regulation benefits for autocatalytic cycles.
From the perspective of the simple model, we recall that $\dot{X}={f}_{a}{f}_{b}$. If the cycle is such that some steady state concentration, ${X}^{*}$, is the desired value for biological function, then for levels of $X$ below ${X}^{*}$ convergence will be faster if ${f}_{a}$ is increased and ${f}_{b}$ is decreased, compared with their values at ${X}^{*}$. Conversely, for levels of $X$ above ${X}^{*}$, convergence will be faster if ${f}_{a}$ is decreased and ${f}_{b}$ is increased, compared with their values at ${X}^{*}$. Convergence to ${X}^{*}$ can therefore be accelerated if, for example, $X$ activates the branch reactions and inhibits the cycle reactions.
The assimilated metabolite can also allosterically regulate the reactions of the cycle. We assume that the desired steady state, denoted $\widehat{X}$, does not depend on the concentration of the assimilated metabolite, $A$. Under this assumption, we further assume that $\widehat{X}$ is attained for some constant concentration of the assimilated metabolite, $\widehat{A}$. It then follows that because the autocatalytic activity is higher when $A\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}\hat{A}$, then in order to maintain ${X}^{*}$ close to its desired level, when $A\phantom{\rule{thinmathspace}{0ex}}>\phantom{\rule{thinmathspace}{0ex}}\hat{A}$, ${f}_{a}$ should be inhibited, and ${f}_{b}$ should be activated, but when $A\phantom{\rule{thinmathspace}{0ex}}<\phantom{\rule{thinmathspace}{0ex}}\hat{A}$, ${f}_{a}$ should be activated, and ${f}_{b}$ should be inhibited. Therefore, to increase the robustness of the steady state concentration to changes in the concentration of the assimilated metabolite, the assimilated metabolite should inhibit the cycle reactions and activate the branch reactions.
Another possible class of regulators are the products of the branch reactions. Taking a somewhat simplistic view, if the level of $Y$, the product of a branch reaction is low, this can indicate that the cycle does not carry sufficient flux to supply the demand for $Y$. Regulation can then be used to increase ${X}^{*}$. From Equation 7, we get that the steady state concentration, ${X}^{*}$, increases as ${K}_{M,b}$ increases and ${V}_{\mathrm{max},b}$ decreases, corresponding to inhibition of ${f}_{b}$, and that ${X}^{*}$ decreases as ${K}_{M,a}$ increases and ${V}_{\mathrm{max},a}$ decreases, corresponding to inhibition of ${f}_{a}$. So, to tune autocatalytic fluxes to match the demands of $Y$, regulation should increase ${X}^{*}$ when $Y$ is low, by activating the recycling and autocatalytic reactions and inhibiting the branch reactions. On the other hand, regulation should decrease ${X}^{*}$ when $Y$ is high, by inhibiting the autocatalytic reactions and activating the branch reactions. Therefore, to synchronize the demand of the cycle product with the cycle flux, the cycle branch products should inhibit the cycle reactions and activate the branch reactions.
Finally, we note that in the autocatalytic cycles we identify in central carbon metabolism, there are also reactions that operate in the reverse direction to the branch reactions, such that they consume products of the cycle and produce intermediate metabolites of the cycle. As such reactions are mirror images of branch reactions, we expect them to be oppositely regulated to branch reactions.
We find that these predictions hold for the cycle using the PTS, that is known to be allosterically controlled, but not for the glyoxylate cycle, which is known to be transcriptionally controlled (Gerosa et al., 2015).
References
 1

2
The path of carbon in photosynthesis. v. paper chromatography and radioautography of the products^{1}Journal of the American Chemical Society 72:1710–1718.https://doi.org/10.1021/ja01160a080
 3
 4

5
Springer Science & Business MediaThe hypercycle: a principle of natural selforganization, Springer Science & Business Media.

6
Portland Press LondonUnderstanding the control of metabolism, Portland Press London.

7
Plant Carbohydrate Biochemistry17–27, Modelling metabolic pathways and analysing control, Plant Carbohydrate Biochemistry.
 8
 9
 10
 11

12
Detecting autocatalytic, selfsustaining sets in chemical reaction systemsJournal of Theoretical Biology 227:451–461.https://doi.org/10.1016/j.jtbi.2003.11.020

13
Enzyme kinetic equations of irreversible and reversible reactions in metabolismJournal of Biosciences and Medicines 02:24–29.https://doi.org/10.4236/jbm.2014.24005

14
EcoCyc: fusing model organism databases with systems biologyNucleic Acids Research 41:D605–D612.https://doi.org/10.1093/nar/gks1027

15
The role and control of the glyoxylate cycle in Escherichia coliThe Biochemical Journal 99:1–11.https://doi.org/10.1042/bj0990001

16
Fructose transport by Escherichia coliPhilosophical Transactions of the Royal Society B: Biological Sciences 326:505–513.https://doi.org/10.1098/rstb.1990.0028
 17
 18

19
Reconstruction and use of microbial metabolic networks: the core Escherichia coli metabolic model as an educational guideEcoSal Plus, 4, 10.1128/ecosalplus.10.2.1, 26443778.

20
Algal evolution in relation to atmospheric CO2: carboxylases, carbonconcentrating mechanisms and carbon oxidation cyclesPhilosophical Transactions of the Royal Society B: Biological Sciences 367:493–507.https://doi.org/10.1098/rstb.2011.0212
 21

22
On the stability of metabolic cyclesJournal of Theoretical Biology 266:536–549.https://doi.org/10.1016/j.jtbi.2010.07.023

23
Signatures of arithmetic simplicity in metabolic network architecturePLoS Computational Biology 6:e1000725.https://doi.org/10.1371/journal.pcbi.1000725
 24
 25

26
The quantitative and conditiondependent Escherichia coli proteomeNature Biotechnology 34:104–110.https://doi.org/10.1038/nbt.3418

27
BRENDA, the enzyme database: updates and major new developmentsNucleic Acids Research 32:431D–433.https://doi.org/10.1093/nar/gkh081
 28

29
Honeybee flight muscle phosphoglucose isomerase: matching enzyme capacities to flux requirements at a nearequilibrium reactionThe Journal of Experimental Biology 200:1247–1254.

30
Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and EngineeringWestview Press.
 31

32
Enumeration of the elementary circuits of a directed graphSIAM Journal on Computing 2:211–216.https://doi.org/10.1137/0202017

33
The danger of metabolic pathways with turbo designTrends in Biochemical Sciences 23:162–169.https://doi.org/10.1016/S09680004(98)012055
 34
Decision letter

Daniel SegreReviewing Editor; Boston University, United States
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Design principles of autocatalytic cycles constrain enzyme kinetics and force overexpression at flux branch points" for consideration by eLife. Your article has been reviewed by four peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Aviv Regev as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Frank Bruggeman (Reviewer #2).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
This article addresses a very interesting and fundamental question about the nature of autocatalytic cycles within cellular metabolic networks, namely whether the existence and stability of steady states in such cycles imposes constraints on the kinetic parameters of the reactions involved. Starting from simple reactions, and gradually looking at slightly more complicated cases, the authors identify general rules that enzymes should follow in order for the autocatalytic cycles to support a stable flux.
All four reviewers agreed that the article is well written and original, and that it provides insightful and valuable results that would be in principle worth presenting to the audience of eLife. However, they also requested a number of important clarifications and revisions, as listed cumulatively below:
Essential revisions:
1) In the equation for "fa", there is no mention of the bisubstrate nature of the reaction. It is understood that the authors consider A constant, which seems a reasonable assumption to start with (though it may be interesting to entertain simple alternatives, e.g. as a constant influx of A – see also next point). However, it is not completely obvious that, even if the kinetic parameters relative to A in the M.M. equations are incorporated into the constants of the equation as shown, the bisubstrate nature of the reaction would not affect the generality of the conclusions drawn throughout the paper.
2) Both in the simple model, and in the extended ones, the authors consider A to be constant. This can be a convenient assumption for the simple model, but makes conclusions for realistic models less relevant, since consumption of A (lowering its level to a critical amount) can be an important part of the system's dynamics. The authors should either explore this alternative possibility, or provide a justification for why this would not be expected to affect significantly the main conclusions.
3) Even though the simple example that is treated in the paper nicely sets the stage and helps the reader in understanding the analytical approach, it is based on irreversible enzyme kinetics only. How do the requirements for proper functioning change when the reactions are catalysed by productsensitive and reversible kinetics? Does product sensitivity prevent loss of stability? How does the displacement from thermodynamic equilibrium of the enzymes at the branch point influence the stability?
4) Allosteric regulation, within or onto, the autocatalytic subnetwork can also improve the robust functioning of the autocatalytic subnetwork. Have the authors analysed this? What did they find? This is particularly relevant in relation to the remark, which suggests that enzyme expression adjustment  so, in fact rate adjustment  can always lead to a nonzero, stable steady state. If so, then is it true that allosteric regulation could always lead to a steady state? Could one then make predictions about the kinetic design of autocatalytic cycles based on their allosteric control?
5) In subsection “Stability analysis for multiplereaction cycles”, the authors mention that the affinity of the branch reaction to the intermediate metabolite of the cycle it consumes must be lower than the affinity of the corresponding recycling reaction of the cycle. Is this statement really correct, given that it is based on either a unscaled sensitivity (dv/ds and not dlnv/dlns) or a saturation degree, which generally also depend on the expression levels of the enzymes, i.e. their Vmaxs?
6) The theoretical predictions could and should be tested more strictly. Although the experimental results that are discussed are consistent with the predictions from the mathematical analysis, they could have alternative explanations. Low saturation of enzymes with their substrates is not uncommon in metabolic networks, and a number of hypotheses have been put forward to explain the apparently wasteful investment in high enzyme concentrations. See for instance refs. [16], not to mention some from this manuscript's authors. In fact, while according to the theory a single branching enzyme being less saturated than the cycle enzyme consuming the same substrate would warrant stability, the experimental data point to most of the branching enzymes being less saturated than the cycle enzymes, which makes the concern about alternative explanations pertinent. The authors may thus wish to discuss these potential alternative explanations, and explain why they eventually don't apply. As a negative control, the authors may wish to examine if enzymes branching from nonautocatalytic cycles tend to be more saturated with substrates from the cycle than those branching from autocatalytic cycles.
7) If lower expression of the enzymes catalyzing reactions that branch out from autocatalytic cycles leads to higher production rates of metabolic precursors, one would expect a negative correlation between the cellular abundance of some of these enzymes and growth rate. This prediction should be relatively straightforward to test based on published quantitative proteomics datasets. Again, a comparison with enzymes branching from nonautocatalytic cycles could be made as control.
8) An examination of the regulation of the branching enzymes may also provide additional support for the theory. E.g., these enzymes may be subject to competitive inhibition (to increase the apparent KM) under conditions where there is a higher demand for the cycle's intermediates.
9) It would also be very useful to have experimental data showing that forcing a lower activity or higher saturation of key branching enzymes would cause instability as predicted. However, this may require dedicated experiments, and should not be viewed as a necessary condition for acceptance.
10) While most of the paper seemed clearly written and unambiguous, several statements in the Introduction are imprecise or unclear. (a) One confusing aspect is the ambiguity between collectively autocatalytic systems and autocatalysis of their components. The component of a system that is collectively autocatalytic is in general not autocatalytic. Thus the statement in the second sentence of the Introduction is inaccurate. In fact, the authors themselves contradict that statement in paragraph two of the Introduction. (b) The authors mention that "autocatalytic systems have an inherent potential to be unstable as their operation changes the amount of their components". This last part of the sentence could be true for many other nonautocatalytic dynamical systems, so it is not clear that this statement is justified. (c)The concept of a system being catalytic "with respect to something" is not clearly defined (and, perhaps, unnecessary for clarifying the concept of a purely metabolic autocatalytic cycle); (d) It is not clear how the statement at the end of paragraph two ("Therefore….") follows as a logical consequence of the prior statement.
11) In the interest of reproducibility, the authors should provide supplementary material tables containing the data used for deriving Figure 6, as well as a description of how exactly the data was processed to generate the final flux estimates.
References
[1]. Weiss, S.L., Lee, E.A. & Diamond, J. (1998) Evolutionary matches of enzyme and transporter capacities to dietary substrate loads in the intestinal brush border. Proceedings Of The National Academy Of Sciences Of The United States Of America, 95, 21172121.
[2]. Suarez, R.K., Staples, J.F., Lighton, J.R.B. & West, T.G. (1997) Relationships between enzymatic flux capacities and metabolic flux rates: Nonequilibrium reactions in muscle glycolysis. Proceedings Of The National Academy Of Sciences Of The United States Of America, 94, 70657069.
[3]. Staples, J.F. & Suarez, R.K. (1997) Honeybee flight muscle phosphoglucose isomerase: Matching enzyme capacities to flux requirements at a nearequilibrium reaction. Journal of Experimental Biology, 200, 12471254
[4]. Salvador, A. & Savageau, M.A. (2003) Quantitative evolutionary design of glucose 6phosphate dehydrogenase expression in human erythrocytes. Proceedings Of The National Academy Of Sciences Of The United States Of America, 100, 1446314468
[5]. Salvador, A. & Savageau, M.A. (2006) Evolution of enzymes in a series is driven by dissimilar functional demands. Proceedings Of The National Academy Of Sciences Of The United States Of America, 103, 22262231
[6]. Eanes, W.F., Merritt, T.J.S., Flowers, J.M., Kumagai, S., Sezgin, E. & Zhu, C.T. (2006) Flux control and excess capacity in the enzymes of glycolysis and their relationship to flight metabolism in Drosophila melanogaster. Proceedings Of The National Academy Of Sciences Of The United States Of America, 103, 1941319418
[7]. Tibor Gánti (2003) "The principles of life" Oxford, Oxford University Press
[8]. Semenov, S.N., Kraft, L.J., Ainla, A., Zhao, M., Baghbanzadeh, M., Campbell, V.E., Kang, K., Fox, J.M. & Whitesides, G.M. (2016) Autocatalytic, bistable, oscillatory networks of biologically relevant organic reactions. Nature, 537, 656660.
https://doi.org/10.7554/eLife.20667.012Author response
Essential revisions:
1) In the equation for "fa", there is no mention of the bisubstrate nature of the reaction. It is understood that the authors consider A constant, which seems a reasonable assumption to start with (though it may be interesting to entertain simple alternatives, e.g. as a constant influx of A – see also next point). However, it is not completely obvious that, even if the kinetic parameters relative to A in the M.M. equations are incorporated into the constants of the equation as shown, the bisubstrate nature of the reaction would not affect the generality of the conclusions drawn throughout the paper.
We thank the reviewers for this (and the following) important points. The reviewers correctly point at an important gap in the paper regarding the bisubstrate nature of the autocatalytic reaction.
We updated the manuscript (subsection “Integrating the bisubstrate nature of the autocatalytic reaction into the simple model”) to include, beyond the unisubstrate case, also four cases of bisubstrate reactions: substituted enzyme (PingPong) scheme, random binding order, ordered binding, with the assimilated metabolite binding first, and ordered binding with the internal metabolite binding first (Materials and methods section "Connecting bisubstrate reaction kinetic constants with simple MichaelisMenten constants"). We show that in all four cases, keeping A constant results in a M.M. like equation, with kinetic parameters relative to A, maintaining the main result of the paper. In three cases a critical amount for A exists, below which a positive stable steady state cannot be attained (Materials and methods section "Constraints on concentration of assimilated metabolite and kinetic constants of bisubstrate reactions").
In the first three cases of bisubstrate binding, bounds on the affinity of the bisubstrate enzyme to the internal metabolite X, exist (compare, e.g. to equation 5), supporting the conclusions drawn throughout the paper. In the last bisubstrate case (ordered binding with the internal metabolite binding first) the expression levels of the autocatalytic and branch enzymes can be modified to achieve positive stable steady state for any set of affinity parameters.
We therefore conclude that the generality of the conclusions drawn throughout the paper hold also with bisubstrate reactions, with a modification in the case of ordered binding with the internal metabolite binding first. Even for this case the conclusions apply, but do not impose constraints on the affinity of the bisubstrate reaction.
2) Both in the simple model, and in the extended ones, the authors consider A to be constant. This can be a convenient assumption for the simple model, but makes conclusions for realistic models less relevant, since consumption of A (lowering its level to a critical amount) can be an important part of the system's dynamics. The authors should either explore this alternative possibility, or provide a justification for why this would not be expected to affect significantly the main conclusions.
As noted in the previous comment, we find that (in all cases but the PingPong reaction scheme) a critical amount for A exists, below which the system loses the positive stable steady state. Interestingly, under some circumstances, an upper bound on the level of A also exists, but this upper bound can be removed by lowering the expression level of the bisubstrate enzyme. The exact critical level of A depends on the bisubstrate reaction scheme and on the specific kinetic parameters and expression levels of the bisubstrate enzyme and the branch enzyme. Nevertheless, qualitatively, the system’s dynamics is identical in all of the reactions schemes we analyze (Materials and methods section "Dependence of steady state concentration on assimilated metabolite"). The concentration of the intermediate metabolite, X, at the st.st. decreases as the level of A decreases (decreasing both the autocatalytic flux, and the branching flux), and reaches 0 when A hits its critical amount (or 0, in the PingPong scheme). Thus, the effect of draining A does not affect the main conclusions, but does add another observation regarding the system’s dynamics when the assimilated metabolite is depleted. We have updated the text to reflect these observations.
3) Even though the simple example that is treated in the paper nicely sets the stage and helps the reader in understanding the analytical approach, it is based on irreversible enzyme kinetics only. How do the requirements for proper functioning change when the reactions are catalysed by productsensitive and reversible kinetics? Does product sensitivity prevent loss of stability? How does the displacement from thermodynamic equilibrium of the enzymes at the branch point influence the stability?
The reviewers address a biologically important extension to the simple example treated in the paper. We added analysis of the requirements for proper functioning under two cases of productsensitive, reversible kinetics:
1) When the branch reaction is reversible (subsection “Reversible branch reaction can either be far from equilibrium, resulting in the simple case, or near equilibrium, pushing the stability conditions down the branch pathway”), the system behavior depends on the assumptions made regarding the consumption of the product of the branch reaction, Y. If Y is removed rapidly, the branch reaction operates far from thermodynamic equilibrium and the same requirements as in the simple example hold. If Y is removed slowly, then we show that the requirements from the irreversible case are ‘forwarded’ to the reaction that consumes Y, requiring its maximal flux to be higher than that of the autocatalytic reaction, and its flux for small values of Y to be smaller than the autocatalytic flux for small X (Materials and methods section "Reversible branch reaction analysis").
2) When the autocatalytic reaction is productsensitive (subsection “Analysis of a Reversible autocatalytic cycle reaction”) we find that the level of X for which the reaction proceeds in the autocatalytic direction is bounded (regardless of the exact bisubstrate reaction mechanism). Therefore, in this case the only requirement for proper functioning is that, for constant A, the initial rate of the autocatalytic reaction for small X is larger than the rate of the branch reaction for the same X, a requirement that can always be satisfied by increasing the expression level of the autocatalytic enzyme. Therefore, in this case, product sensitivity prevents loss of stability.
4) Allosteric regulation, within or onto, the autocatalytic subnetwork can also improve the robust functioning of the autocatalytic subnetwork. Have the authors analysed this? What did they find? This is particularly relevant in relation to the remark, which suggests that enzyme expression adjustment  so, in fact rate adjustment  can always lead to a nonzero, stable steady state. If so, then is it true that allosteric regulation could always lead to a steady state? Could one then make predictions about the kinetic design of autocatalytic cycles based on their allosteric control?
We thank the reviewers for highlighting this important aspect of metabolic control, overlooked in our work. We find allosteric regulation to be valuable in optimizing native autocatalytic cycles enabling both faster convergence to st.st. and modulation of affinities such that branch enzymes can be more efficient, by having higher affinity, when the corresponding autocatalytic subnetwork is not functioning (Materials methods section 13). We note that, as our work also targets synthetic design and metabolic engineering leading to the introduction of novel autocatalytic cycles, we find allosteric regulation implementation to be more complicated in such contexts and thus less appealing a strategy compared with, for example, the introduction of external input fluxes.
We added two sections to the paper analyzing allosteric regulation. One section (Materials and methods section "Allosteric regulation can improve network performance") describes the allosteric regulation interactions (inhibitions or activations) that are expected to improve the function of autocatalytic cycles. Another section (” Analysis of Allosteric regulation potential for cycle improvement”) compares the predicted interactions with the actual direction (activation or inhibition) found for the two autocatalytic cycles we identify and find to be functional in the central carbon metabolism of native E. coli (the PTS using cycle and the glyoxylate cycle). Briefly, for the PTS using cycle we find good agreement between our predicted direction of regulation and the regulation found (11/12 interactions agree with our prediction). For the glyoxylate cycle, the agreement is much weaker (8/13 interactions agree with our prediction). We note that this finding is in agreement with the work of Gerosa et al. (Cell Systems 2015) that suggests that TCA cycle fluxes are mostly regulated by transcription and not reactant levels.
Clearly, allosteric regulation serves many more roles and the metabolic network faces many more challenges than just the support of stable autocatalysis. We do not therefore suggest that proper autocatalytic function is the underlying reason for the interactions that agree with our predictions.
5) In subsection “Stability analysis for multiplereaction cycles”, the authors mention that the affinity of the branch reaction to the intermediate metabolite of the cycle it consumes must be lower than the affinity of the corresponding recycling reaction of the cycle. Is this statement really correct, given that it is based on either a unscaled sensitivity (dv/ds and not dlnv/dlns) or a saturation degree, which generally also depend on the expression levels of the enzymes, i.e. their Vmaxs?
We thank the reviewers for showing us we did not make this delicate point clear. We now address the issue of the connection between the unscaled sensitivities of reactions, their saturation degrees, and their affinities in the Materials and methods, section "Inverse relationship between derivatives, affinities, and saturation levels". According to the lemma we prove, we conclude that if f_a(X*)>f_b(X*) and β>α, then it holds that S_a(X*)>S_b(X*) and K_(M,a)<K_(M,b), which correspond to the statements we make in subsection “Analysis of a Reversible autocatalytic cycle reaction”.
6) The theoretical predictions could and should be tested more strictly. Although the experimental results that are discussed are consistent with the predictions from the mathematical analysis, they could have alternative explanations. Low saturation of enzymes with their substrates is not uncommon in metabolic networks, and a number of hypotheses have been put forward to explain the apparently wasteful investment in high enzyme concentrations. See for instance refs. [16], not to mention some from this manuscript's authors. In fact, while according to the theory a single branching enzyme being less saturated than the cycle enzyme consuming the same substrate would warrant stability, the experimental data point to most of the branching enzymes being less saturated than the cycle enzymes, which makes the concern about alternative explanations pertinent. The authors may thus wish to discuss these potential alternative explanations, and explain why they eventually don't apply. As a negative control, the authors may wish to examine if enzymes branching from nonautocatalytic cycles tend to be more saturated with substrates from the cycle than those branching from autocatalytic cycles.
We thank the reviewers for this point. We agree that alternative explanations may underlie the lower saturation of branching enzymes from autocatalytic cycles.
We added references to the excellent works the reviewers pointed at.
We address this issue in three independent ways:
1) We show that having multiple unsaturated reactions increases the robustness of the cycle to perturbations, as is suggested in minor point 7 below (Materials and methods section "Multiple unsaturated branch reactions increase convergence speed and dampen oscillations").
2) We compare the saturation levels of the same branch points as those presented in Figure 6 for conditions at which the autocatalytic cycles do not function, but the relevant reactions carry flux. We find that in 5 out of 9 cases the saturation levels are reversed, meaning that in these cases the cycle reactions are less saturated than the branch reactions (subsection “Testing the predictions of the analysis with experimental data on functioning autocatalytic cycles” and supplementary file 1).
3) We consider branching points from two nonautocatalytic cycles – the TCA cycle and the PP + gluconeogenesis cycle (composed of the PP pathway, coupled with gluconeogenic flux in upper glycolysis). For each such cycle we identify the major branch points (akg consumption from the TCA cycle and r5p consumption from the PP cycle). We find that out of 6 branch cases, in 3 cases the branch reaction is less saturated than the cycle reaction, and in 3 cases, the situation is reversed (subsection “Testing the predictions of the analysis with experimental data on functioning autocatalytic cycles” and supplementary file 1).
To conclude, we find that the bias towards lower saturation of the branch reaction is specific to cases of functioning autocatalytic cycles and is not present neither for the same pairs of reactions in the absence of autocatalysis, nor in branches out of non autocatalytic cycles. We further provide a mathematical analysis highlighting the advantage of multiple unsaturated branch reactions that enable more robust functioning of the cycle.
7) If lower expression of the enzymes catalyzing reactions that branch out from autocatalytic cycles leads to higher production rates of metabolic precursors, one would expect a negative correlation between the cellular abundance of some of these enzymes and growth rate. This prediction should be relatively straightforward to test based on published quantitative proteomics datasets. Again, a comparison with enzymes branching from nonautocatalytic cycles could be made as control.
The reviewers suggest an interesting, additional validation of the findings of the paper. However, we do not find this suggestion to be straightforward to test. In order to be tested, it requires multiple conditions with different growth rates, under all of which the same autocatalytic cycle functions and the expression levels of the branch or cycle reactions differ. Glucoselimited chemostat data may seem to be an adequate candidate, but there are a few caveats for using it:
1) The level of the assimilated metabolite, glucose, changes dramatically between the different chemostat conditions.
2) Different transporters exist and function under different concentrations of glucose, not necessarily following the autocatalytic properties of the PTS.
3) As noted in our response to the points regarding allosteric regulation, it may play a major role in controlling fluxes in the PTSusing autocatalytic cycle, which is the cycle used in growth on glucose.
Indeed, evaluating the correlations for two data sets (from Schmidt et al., used also for the expression levels in our article, and from Peebo et al. Mol. Biosyst. 2015) shows that the branch and corresponding recycling enzymes display various dependencies on the growth rate, mostly even nonuniform (that is not monotonically increasing or monotonically decreasing). We believe a specific experimental system, titrating the expression levels of various enzymes along the cycle while maintaining the level of external metabolite fixed, is probably the best way to investigate this phenomenon, but find it to lie outside the scope of this work.
8) An examination of the regulation of the branching enzymes may also provide additional support for the theory. E.g., these enzymes may be subject to competitive inhibition (to increase the apparent KM) under conditions where there is a higher demand for the cycle's intermediates.
We thank the reviewers for highlighting this potential additional support. We address this issue together with addressing point 4 above and find that, in the PTSusing cycle, such evidence indeed exists. In the glyoxylate cycle such regulation is not statistically significant, suggesting other control mechanisms may be at play (subsection “Analysis of Allosteric regulation potential for cycle improvement” and supplementary file 1).
9) It would also be very useful to have experimental data showing that forcing a lower activity or higher saturation of key branching enzymes would cause instability as predicted. However, this may require dedicated experiments, and should not be viewed as a necessary condition for acceptance.
We agree with the suggestion the reviewers make here, and find such a system to also be useful for investigating point 7 above further.
10) While most of the paper seemed clearly written and unambiguous, several statements in the Introduction are imprecise or unclear. (a) One confusing aspect is the ambiguity between collectively autocatalytic systems and autocatalysis of their components. The component of a system that is collectively autocatalytic is in general not autocatalytic. Thus the statement in the second sentence of the Introduction is inaccurate. In fact, the authors themselves contradict that statement in paragraph two of the Introduction. (b) The authors mention that "autocatalytic systems have an inherent potential to be unstable as their operation changes the amount of their components". This last part of the sentence could be true for many other nonautocatalytic dynamical systems, so it is not clear that this statement is justified. (c)The concept of a system being catalytic "with respect to something" is not clearly defined (and, perhaps, unnecessary for clarifying the concept of a purely metabolic autocatalytic cycle); (d) It is not clear how the statement at the end of paragraph two ("Therefore….") follows as a logical consequence of the prior statement.
We thank the reviewers for pointing out these confusing points regarding the general introduction to autocatalysis. As we found this part of the Introduction does not contribute to the main theme of the paper, we removed it from the revised version.
11) In the interest of reproducibility, the authors should provide supplementary material tables containing the data used for deriving Figure 6, as well as a description of how exactly the data was processed to generate the final flux estimates.
The reviewers correctly point at the missing tables, which we added in the revised version.
https://doi.org/10.7554/eLife.20667.013Article and author information
Author details
Funding
Israel Science Foundation (740/16)
 Uri Barenholz
 Dan Davidi
 Yinon BarOn
 Niv Antonovsky
 Ron Milo
BeckCanadian Center for Alternative Energy Research
 Uri Barenholz
 Dan Davidi
 Yinon BarOn
 Niv Antonovsky
 Ron Milo
Dana and Yossie Hollander
 Uri Barenholz
 Dan Davidi
 Yinon BarOn
 Niv Antonovsky
 Ron Milo
Leona M. and Harry B. Helmsley Charitable Trust
 Uri Barenholz
 Dan Davidi
 Yinon BarOn
 Niv Antonovsky
 Ron Milo
The Larson Charitable Foundation
 Uri Barenholz
 Dan Davidi
 Yinon BarOn
 Niv Antonovsky
 Ron Milo
Wolfson Family Charitable Trust
 Uri Barenholz
 Dan Davidi
 Yinon BarOn
 Niv Antonovsky
 Ron Milo
Charles Rothchild
 Uri Barenholz
 Dan Davidi
 Yinon BarOn
 Niv Antonovsky
 Ron Milo
Selmo Nussenbaum
 Uri Barenholz
 Dan Davidi
 Yinon BarOn
 Niv Antonovsky
 Ron Milo
Alternative Sustainable Energy Research Initiative (Graduate Student Fellowship)
 Uri Barenholz
European Research Council (NOVCARBFIX 646827)
 Uri Barenholz
 Dan Davidi
 Yinon BarOn
 Niv Antonovsky
 Ron Milo
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We would like to thank Arren BarEven, Katja Tummler, Matthias Heinemann, David Fell, Patrick Shih, Noam Prywes, Leeat Keren, David Wernick, Wolfram Liebermeister, Rami Pugatch and Ron Sender for fruitful discussions and valuable insights contributing to this work. We would also like to thank Frank Bruggeman and the other anonymous reviewers for their significant contributions to the completeness and coherence of this work.
Reviewing Editor
 Daniel Segre, Boston University, United States
Publication history
 Received: August 15, 2016
 Accepted: February 1, 2017
 Accepted Manuscript published: February 7, 2017 (version 1)
 Version of Record published: March 2, 2017 (version 2)
Copyright
© 2017, Barenholz et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,124
 Page views

 312
 Downloads

 7
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Download citations (links to download the citations from this article in formats compatible with various reference manager tools)
Open citations (links to open the citations from this article in various online reference manager services)
Further reading

 Cell Biology
 Computational and Systems Biology

 Computational and Systems Biology
 Ecology