Laboratory evolution reveals a two-dimensional rate-yield tradeoff in microbial metabolism

Growth rate and yield are fundamental features of microbial growth. However, we lack a mechanistic and quantitative understanding of the rate-yield relationship. Studies pairing computational predictions with experiments have shown the importance of maintenance energy and proteome allocation in explaining rate-yield tradeoffs and overflow metabolism. Recently, adaptive evolution experiments of Escherichia coli reveal a phenotypic diversity beyond what has been explained using simple models of growth rate versus yield. Here, we identify a two-dimensional rate-yield tradeoff in adapted E. coli strains where the dimensions are (A) a tradeoff between growth rate and yield and (B) a tradeoff between substrate (glucose) uptake rate and growth yield. We employ a multi-scale modeling approach, combining a previously reported coarse-grained small-scale proteome allocation model with a fine-grained genome-scale model of metabolism and gene expression (ME-model), to develop a quantitative description of the full rate-yield relationship for E. coli K-12 MG1655. The multi-scale analysis resolves the complexity of ME-model which hindered its practical use in proteome complexity analysis, and provides a mechanistic explanation of the two-dimensional tradeoff. Further, the analysis identifies modifications to the P/O ratio and the flux allocation between glycolysis and pentose phosphate pathway (PPP) as potential mechanisms that enable the tradeoff between glucose uptake rate and growth yield. Thus, the rate-yield tradeoffs that govern microbial adaptation to new environments are more complex than previously reported, and they can be understood in mechanistic detail using a multi-scale modeling approach.


Introduction
Growth rate and yield are basic features of microbial life that are widely implicated in cell fitness, adaptation, and evolution [1]. The specific growth rate, μ, represents the number of doublings of bacterial density per unit time [10]. The yield, Y, is the ratio between μ and the rate of substrate consumption [10,11]. The mathematical relation between μ and Y can be written as: where M substrate is the molecular weight of the substrate and q substrate is the substrate uptake rate.
In the context of modeling the phenotypic relation between substrate uptake, metabolism, and biomass growth, there is a great interest in developing quantitative descriptions of the relationship between μ and Y. The wide-ranging measurements of μ and Y ( Fig 1A) across microbial communities and environments raised interest into the exact nature of the μ-Y relationship [1]. At low μ, positive correlations between μ and Y have been observed [8], and these can be explained by non-growth-associated cell maintenance requirements that make slow growth inefficient [11]. At high μ, negative correlations between μ and Y are observed [5], and for E. coli, this can be explained by a tradeoff between metabolic efficiency and enzymatic efficiency that lead to decreasing Y at high μ [12,13]. In particular, E. coli exhibits a tradeoff between respiration, which has higher energy yield per carbon substrate (more metabolicallyefficient), and acetate fermentation, which requires less enzyme per carbon substrate (more proteome-efficient). Therefore, acetate excretion increases linearly with μ above a threshold growth rate (green and blue lines in Fig 1B) [5]. [1] summarized these observations where positive μ-Y correlation at low μ and negative μ-Y correlation at high μ are different parts of a bell-shaped μ-Y curve ( Fig 1A). However, recent experiments suggest that adaptation to new environments can modify the bell-shaped μ-Y tradeoff [3,14,15].
Microorganisms rapidly adapting to environmental niches [17,18], and adaptation mechanisms can be studied directly through adaptive laboratory evolution (ALE) [19]. When strains are adapted through ALE for growth in a liquid minimal medium, they achieve higher μ compared to the wild-type (Fig 1A), ALE-adapted strains have been shown to rapidly acquire regulatory mutations that modify proteome allocation, but they do not acquire new metabolic capabilities within the time frame of reported short-term (4 to 8 weeks) adaptation experiments [3,15,20]. By analyzing ALE-adapted strains, we can reveal the strategies that allow cells to optimize their proteome allocation for growth in an environmental niche, subject to the constraints of their metabolic capabilities (i.e. their repertoire of pathways) and constraints on the kinetic efficiencies of their enzymes [3,20,21].
Contrary to the negative μ-Y relationship at high μ observed in wildtype strains, the endpoint strains of ALE experiments of E. coli selected for high μ in a minimal medium do not have μ-Y data points aligning on the bell-shape curve in Fig 1A, but reveal an uncorrelated  Tables  (Supporting Information). (A-D) Phenotypic data for E. coli strains including Y, μ, q ac , and q glc data. The Y is calculated by m M glc �q glc , where the molecular weight of glucose M glc = 180.156g/mol. Two datasets are presented from literature, for chemostat growth [8] (green triangles) and substrate titration [5] (blue squares). These are compared to strains adapted for maximum growth rate through ALE (this study; red circles; error bars for standard deviation across duplicates). The bell-shaped μ-Y relationship proposed by [1] is included for reference. (E) Diagram of the SSME-model derived from [5]. The model consists of three pathways: respiration (R, res) and fermentation (F, fer) generate different amounts of energy, feeding the biomass (B, bms) pathway to synthesize biomass.(F) Diagram of the genome-scale ME-model that includes a genome-scale reconstruction of metabolic pathways ("M") and protein expression machinery ("E") [9,16]. Only central carbon metabolism is sketched, other biosynthesis pathways (AA amino acids, NT nucleotides and others) are simplified as dashed arrows in the plot, but remain fine-grained in the genome-scale ME-model that we analyze in this study. For the "E" pathway, it starts with the amino acids and nucleotides that are synthesized through "M" pathway, and the end product of "E" pathway is the dilution of the protein (enzyme). The dilution rate of the enzyme determines the corresponding metabolic reaction rate with the factor k eff . More details about k eff constraints can be found in "3 Proteome constraints in the ME-model" in S1 Appendix. relationship between μ and Y [3,15]. In these experiments, strains exhibited little variation in μ but high variation in Y and acetate excretion rate q ac . Previous studies similarly reported that overflow metabolism can be nearly eliminated through genetic engineering without any effect on growth rate in E. coli [2,4]. Thus, the negative μ-Y correlation at high growth rates does not appear to be a fundamental constraint on fast-growing cells. A mechanistic model of the full μ-Y relationship must be able to reconcile the bell-shaped curve observed for individual strains with the uncorrelated μ-Y phenotypes seen in ALE-adapted strains (Fig 1A).
A number of theoretical and computational models have been developed to describe rateyield tradeoffs. For the positive μ-Y correlation, maintenance requirements can be quantitatively described using algebraic growth laws [8,11]. The maintenance requirement has been modeled as non-growth associated maintenance (NGAM) in the genome-scale models (GEMs) of metabolism, which can be simulated as an optimization problem, predicting μ and Y when substrate uptake rates (e.g. q glc ) are known [22]. For the negative μ-Y correlation, quantitative models of overflow metabolism have been developed [5][6][7]. In particular, quantitative measurements of E. coli growth in well-controlled environments revealed a linearthreshold response of acetate excretion (q ac ) with increasing μ [5]. To represent the full range of the μ-Y relationship, a constraint allocation flux balance analysis model (CAFBA) was reported that combines a GEM with proteome allocation constraints [6]. A similar solution can be formulated from a bottom-up reconstruction of metabolism and macromolecular expression (ME-model, [9]) that incorporates the protein synthesis pathways into a GEM and applies coupling constraints related to enzyme kinetics parameters on each individual reaction. However, none of these models have been used to explain experiments where μ and Y are decoupled through laboratory evolution or genetic engineering.
In this study, we show that the wide range of μ-Y observations in E. coli can be explained by a two-dimensional rate-yield tradeoff, where the first dimension is the characteristic μ-Y tradeoff associated with acetate overflow metabolism and the second dimension is a tradeoff between glucose uptake rate (q glc ) and Y that appears during ALE adaptation. We employ a multi-scale modeling approach to provide a mechanistic description of the two-dimensional rate-yield tradeoff. By deriving the relationship between the ME-model and the previously reported small-scale proteome allocation model [5], we are able to develop a workflow for modifying ME-model parameters to fit experimental data, and we achieve quantitative predictions for simulations of μ-Y (the first dimension of the rate-yield tradeoff). This multi-scale modeling approach predicts a two-dimensional rate-yield tradeoff, and it suggests that the second dimension of the tradeoff can be explained by changes in P/O ratio and the flux balance between glycolysis and pentose phosphate pathway. This multi-scale modeling approach predicts the systemic response of the cell to growth selection, representing the relationships between P/O ratio, glycolytic-PPP flux balance, and the two dimensions of the rate-yield tradeoff.

Adaptive laboratory evolution reveals a two-dimensional rate-yield tradeoff
To explore the metabolic constraints on E. coli growth, adaptive laboratory evolution (ALE) was used to adapt E. coli K-12 MG1655 to maximize growth at 37˚C in a liquid culture with a minimal medium containing glucose [3]. Eight independent experiments were performed on an automated ALE platform to achieve 8.3 × 10 12 to 18.3 × 10 12 cumulative cell divisions [23]. Phenotype characterization was performed on eight ALE endpoint strains, including quantitative measurements of μ, q glc , q ac , and other common metabolic byproducts of E. coli (Materials and methods).
A diversity of metabolic phenotypes was observed in the ALE endpoint strains. Through ALE, μ increased from 0.7 h -1 for wild-type (red triangles in Fig 1A-1D) to 0.95-1.10 h -1 (red circles with error bars in Fig 1A-1D). Based on previous reports, we expected a linear relationship between μ and q ac . However, ALE endpoint strains achieved a wide ranging q ac from 3.9-11.4 mmol gDW -1 h -1 (where wild-type q ac was 3.9 mmol gDW -1 h -1 . While we did not observe a correlation between μ and q ac in these strains (Fig 1B), there was a clear correlation between q glc and q ac (Fig 1D).
Two of the ALE endpoint strains with similar μ (3% difference) but distinct Y (30% difference) have been processed for 13 C metabolic flux measurements (S11 Table). The measured metabolic fluxes using 13 C metabolic flux analysis ( 13 C MFA, see "Materials and methods") shows the positive correlation between TCA fluxes, q TCA , and Y for the ALE strains. For the ALE strain with larger Y, q glc and q ac are lower and q TCA is higher. Therefore, for a fixed μ, Y increases as q TCA increases and q ac decreases, indicating a pathway switch between the TCA cycle and acetate overflow depending on q glc . Therefore, combining with the referenced study [5], for a wild-type strain, there is a μ-Y tradeoff. And for the isogenic ALE strains, a q glc -Y tradeoff appears. For both tradeoffs, Y varies with the pathway switch between TCA cycle and acetate overflow. In this paper, we call the μ-Y and q glc -Y tradeoffs a two-dimension rate-yield tradeoff, since they are tradeoffs between different "rates" (growth rate μ and glucose uptake rate q glc ) and the same yield (glucose yield Y), and they share the same phenotypic behavior of TCA-aceate overflow pathway switch.
Correlations between q glc , Y, and q ac have been observed previously for E. coli strains [2][3][4], and moreover, a bacterial engineering approach has been reported to vary q ac by manipulating the substrate uptake system [24]. In one of these studies, [4] showed that switching electron transport chain (ETC) enzyme selection (and thereby modifying the P/O ratio) can cause a q glc -Y tradeoff at a low μ of 0.15 h -1 . ALE gained q ac and Y decoupled from μ, which seemingly differs from the reported correlation between μ-Y and μ-q ac [5,8]. The ME-model used in this study simulates the relationships between these q glc -q ac and q glc -Y tradeoffs, connecting to the mechanisms of μ-Y tradeoffs (the bell-shaped curve in Fig 1A) by established models [5,6].
To enable our analysis, it is important to note that ALE endpoint strains rapidly acquire regulatory mutations, but they do not acquire new metabolic capabilities within the time frame of these experiments [3,15,20]. The linear correlation between q ac and μ reported previously was identified for an isogenic strain [5,8]. In contrast, our observations of a decoupling between q ac and μ appear when comparing adapted strains. However, because these adapted strains have only regulatory mutations, their phenotypes represent the limits of what E. coli cells can achieve while bounded by metabolic and proteomic constraints (but not by regulation). This type of adaptation and the associated phenotypic tradeoffs are useful for understanding cellular adaptation to ecological niches where regulatory adaptation can occur rapidly [18].

ME-model data fitting with a multi-scale modeling approach
To explain these experimental observations, we sought a modeling approach that could quantitatively predict the μ-Y and μ-q ac relationships. Our modeling approach starts with fitting the linear-threshold (blue line in Fig 1B) mu-q ac relation [5] using the framework of ME-model [9,16]. We first considered a previously reported coarse-grained model of proteome allocation [5] that describes E. coli overflow metabolism ( Fig 1E). [5] solves Y and q ac as functions of μ, and assuming that the cells pick the maximum Y under each particular μ. This indicates that high-Y growth strategies have a fitness benefit in spatially structured environments, for instance, the wild-type cultures collected from colonies, that has been demonstrated through a Y-selection system [14], and more efficient strategies also leave more resources for cells that are hedging against future stresses [20]. The evolutionary history of E. coli includes growth in structured environments and a wide range of stresses that could have placed a selection pressure on increasing Y. Therefore, we focused on fitting the observed wild-type chemostat [8] and uptake titration [5] data for the Y-maximized growth solution (green and blue data points in Fig 2). Simulations were fit to experimental data for each of the three datasets, K-12 MG1655 chemostat [8] (green triangles), NCM3722 substrate titration [5] (blue squares), and strains adapted from wild-type K-12 MG1655 (red triangle, [3] for maximum growth rate through ALE (this study, red circles, error bars for standard deviation across duplicates). The Y-maximized solutions are displayed as solid lines in all plots. Solution spaces are simulated by taking the feasible range between maximum (Y max -Y min in A and C, q ac,max -q ac,min in B and D)For both models, fitting was performed by manipulating three global parameters: unmodeled protein fraction (UPF), growth-associated maintenance (GAM), and non-growth associated maintenance (NGAM). Details of the fitting approach are provided in Materials and Methods. https://doi.org/10.1371/journal.pcbi.1007066.g002 Laboratory evolution reveals a two-dimensional rate-yield tradeoff in microbial metabolism The coarse-grained proteome allocation model [5] was intended to make predictions at high μ and thus only captures the negative μ-Y relation (Fig 2A). The parameters in the coarse-grained model have a strong experimental basis in fine-grained protein abundances measurements in high growths, and the model simulates accurate predictions of μ-q ac [5].
We also considered the genome-scale ME-model iJL678-ME [16]. With the default parameter settings in the ME-model, simulations had a poor quantitative prediction [9] of μ-q ac to the uptake titration data (S3F Fig). As [5] has shown experimentally that the overflow metabolism is fundamentally caused by the tradeoff between metabolic efficiency (reaction stoichiometry) and protein efficiency (enzyme turnover rate). Since the reaction stoichiometry in the MEmodel has been mass-balanced and well established, we suspect that this poor fit from the MEmodel can be explained by inaccurate genome-wide enzyme turnover rates (k eff s) that MEmodel researchers have been seeking to improve [16,25,26]. We sought to modify the k eff s to fit the μ-q ac data. However, since each of the 5266 reactions in the genome-scale ME-model has a k eff parameter, it is difficult to directly fit the parameters to measured data.
Therefore, we pursued a multi-scale modeling approach where the coarse-grained model was used to analyze the effects of proteome-efficiency at the level of coarse-grained pathways instead of each individual reaction, which helps to tune the fine-grained parameters in the ME-model. To connect the coarse-grained and fine-grained models, we first found that the proteome efficiency (ε) parameters in the coarse-grained model share a conceptual basis with the enzyme efficiency parameter k eff s in ME-models ("3 Proteome constraints in the MEmodel" in S1 Appendix). Thus, we were able to reformulate the coarse-grained model within the framework as the ME-model (S1 Fig). The resulting small-scale ME-model (SSME-model) has parameters directly analogous to those in the genome-scale ME-model (See "5 SSMEmodel parameters derivation" and "6 Matlab and COBRAme implementation" in S1 Appendix). The resulting SSME-model generates identical μ-Y and μ-q ac predictions to the proteome allocation model. The SSME-model is a good tool for k eff parameter sensitivity analysis [27], which provides insights on how to modify the k eff of the ME-model to achieve quantitative fit. As a result, we gained predictions for μ-Y and μ-q ac from both the SSME-and ME-models (blue curves in Fig 2). Details of the ME-model modifications are in the S1 Appendix ("8 Experimental data fitting"). In summary, with the multi-scale modeling approach, we identified the reactions whose enzymes turnover rates are too high to match the observed phenotypes. Those reactions are involved in different pathways, including the TCA cycle, Entner-Doudoroff pathway, glyoxylate shunt, nucleotide salvage, and fatty acids metabolism (S3 Table). Three global parameters, unmodeled protein fraction (UPF), growth-associated maintenance (GAM), and non-growth-associated maintenance(NGAM) (S2 Table) were then used to predict the phenotype from different strains (green, blue and red curves in Fig 2).
The reason for only modifying global parameters to simulate the ALE adaptation is that the mutations in the ALE strains do not directly related to the enzyme turnover rate (k eff value) of a particular metabolic reaction. According to previous ALE studies [3], most mutations occur in genes associated with regulations or translations. Even in the cases where mutations might directly change a k eff , this is hard to model. Therefore, rather than exploring the mechanistic effects of ALE mutations, we focused on the phenotypic changes in the endpoint strains. Some recent studies have shown how individual mutations can have wide-reaching effects on gene expression, metabolic pathway activity, and cell phenotype [3,20].
The most obvious difference between the SSME-model derived from [5] and ME-model for these phenotypic predictions is the expanded solution space of the ME-model (Fig 2). However, much of the ME-model solution space corresponds to very low yield metabolic solutions. If Y is maximized during simulations of the SSME-model and ME-model (achieved by minimizing q glc at a given μ), the resulting predictions are more similar between the models and lie closer to experimental data (solid blue curves in Fig 2). For the growth-rate-dependent Y-maximized solutions (solid curves in all 4 panels), though the q ac lines looks similar, the Y lines looks very different in the low growth regime. Where the SSME-model predicts a constant high Y, the ME-model predicts an initially low Y that increases rapidly with μ. This is because the additional non-growth maintenance energy (NGAM) added in the ME-model. We can also see that in Fig 2C, among the three different solution curves from the ME-model, the curve with lower NGAM has a higher Y at low μ. The NGAM parameters for the different curves are shown in S1 Table in Supporting information. In fact, by adding complexity to the SSME-model or simplifying the ME-model, many intermediate models can be built.

The ME-model predicts phenotypic diversity in ALE strains
As a result of data fitting, we achieved a quantitative fit of chemostat [8] and batch [5] uptake titration data with the Y-maximized ME-model solutions (blue and green curves in Fig 2C and  2D). The ALE-adapted strains (red circles in Fig 2) do not align well with the Y-maximizing solutions (red curves in Fig 2), but they are encompassed by the ME-model solution space. Further analysis of these ALE data points and the corresponding ME-model solutions were used to understand the phenotypic diversity of these adapted strains.
Feasible solutions other than the Y-maximized solution are achieved through the activation of alternative metabolic pathways which are sub-optimal. The SSME-model does not capture the ALE data points with high q ac (red region in Fig 2B), while the genome-scale ME-model does (red region in Fig 2D). Moreover, the ME-model predicts feasible growth at lower Y in the μ-Y solution space than the SSME-model. We sought to determine which pathways are responsible for the lower Y and higher q ac in ME-model that was not captured by the SSMEmodel.
Removing reactions from the ME-model can decrease the size of the solution space (S5, S6 and S7 Figs, "8 Solution space variation" in S1 Appendix), making the solution space more similar to the SSME-model solution space. We employed a workflow to identify 24 reactions (S6 Table) that are not activated in the Y-maximized solutions but are used to enable higher q ac at lower Y. We observed that these 24 reactions are part of metabolically inefficient pathways that are alternatives to the Y-optimal pathways. By extension, metabolically inefficient pathways can be added to the SSME-model to increase the size of the solution space (S7 Fig), making it more similar to the ME-model solution space. Thus, the modified SSME-model can achieve low Y (S7A Fig) at high q ac (S7C Fig). Therefore, the difference in predictions of MEmodel from the SSME-model is a result of the greater range of metabolic capabilities of the genome-scale model.

The two-dimensional rate-yield tradeoff
We can now put forward a theory to connect the correlations in μ-Y (Fig 1A) (and the associated acetate curve in q ac --Y, Fig 1B) with the negative correlation in q glc --Y ( Fig 1C) and positive q glc -q ac correlations (Fig 1D).
To see the relationship between the three variables μ, q glc , and Y we generated ME-model solution spaces in q glc and Y at increasing lower bounds of μ (Fig 3A). These solution spaces represent the flexibility in the model to achieve a particular growth rate. At the Y-maximized limit of these solution space, we see the established negative μ-Y tradeoff where increasing growth rate requires increasing q glc and decreasing Y (dashed arrow marked as "d1" in Fig 3A) coupling with increasing q ac (top edges of solution spaces in Fig 3B). This is the first dimension of the rate-yield tradeoff, "d1". , the growth rate of each ALE endpoint strain is labeled as black number beside the corresponding data point. (A) Two dimensions of the rate-yield tradeoff. The first dimension "d1" is the negative μ-Y correlation at maximum Y, and the second dimension "d2" is the negative q glc -Y correlation at a fixed μ. These correlations are observed in ME-model simulations and experimental data from ALE strains. (B) A correlation between q glc and q ac is also observed at fixed μ in both the ME-model and ALE endpoint data. Linear fits for the experimental data at similar growth rates are shown as dash-dotted (μ = 0.95-0.97 h -1 ), dashed (μ = 1.00-1.04 h -1 ), and solid (μ = 1.08-1.10 h -1 ) orange lines. These fits are described by the upper edges of the q glc -q ac solution space at fixed μ. For growth between 1.00 and 1.04 h -1 , r 2 = 0.931 and p = 0.035. For growth between 1.08 and 1. Laboratory evolution reveals a two-dimensional rate-yield tradeoff in microbial metabolism Considering only "d1", one would expect the acetate production rate in all strains to be fully defined by the growth rate. In the case of this 1-dimensional tradeoff, all points in Fig 3A  would appear on the line at the top of the blue solution spaces (parallel to the dotted "d1" line). However, we observed another degree of freedom in the phenotypic space. At a given μ, evolved strains can acquire higher q glc , higher q ac , and lower Y. The "d2" tradeoff is defined by a linear correlation in q glc -q ac (Fig 3B) and a corresponding inverse proportional tradeoff in q glc -Y (Fig 3A).
The "d2" tradeoff is also predicted by ME-model simulations. At a given μ, the ME-model solution spaces extend toward lower Y and higher q glc , revealing this inverse proportional relationship in q glc -Y. The second dimension "d2" can also be seen in q ac -q glc where the MEmodel predicts the q ac -q glc correlation observed in ALE endpoint as the q ac -maximized edges of the solution spaces (Fig 3B). The solution spaces predicted by ME-model show broad feasible ranges of acetate production q ac at a given q glc and μ ("bold" solution spaces in Fig 3B), so the q glc -Y tradeoff is not required by the model. On the other hand, the relationship between q glc and Y is a strict tradeoff in the model ("thin" solution spaces in Fig 3A). Therefore, the ME-model suggests that q glc -Y is the more fundamental second dimension of the rate-yield tradeoff. To verify that hypothesis, one would look for mutant strains where q glc increased while the other three phenotypic variables remained fixed (a shift to the right in Fig 3B).

Mechanisms for the additional rate-yield tradeoff
We sought to identify the particular alternate metabolic strategies in the ME-model that could enable a q glc -Y tradeoff by identifying the differential pathway usage at a fixed high μ (1.05 h -1 in the ME-model (Fig 3C). The model predicts that when q glc increases from the Y-maximized state (minimum q glc ), flux through the proton-coupled NAD(P) transhydrogenase increases (reaction THD2pp, catalyzed by pntAB. In addition, a pathway switch between two different NADH dehydrogenase reactions, NADH5 (ndh and NADH16pp (nuo, appears at high q glc . In fact, each of or any combination of the 24 reactions in S6 Table can be activated in the MEmodel to achieve high q glc , high q ac , and low Y. There are two common threads among these pathway activations. First, they all reduce the P/O ratio in the simulations (Fig 3C). NADH5 transports fewer protons to the periplasm per electron than NADH16pp. And increasing THD2pp flux drains the proton gradient without contributing to ATP production, thereby reducing P/O ratio (Fig 3D). Second, with the activation of those 24 reactions, glycolytic flux increases ( Fig 3D) and pentose phosphate pathway flux decreases (Fig 3C). By comparing to the 13 C metabolic flux analysis (Fig 3E and 3F), the ME-model shows quantitative predictive power for the second-order rate-yield tradeoff.
Experiments that introduce proton leakage have shown a shift towards high q ac and low Y [5] in the same μ. It has also been shown that the variation of P/O ratio can uncouple the regulation of cytochrome oxidase from the cellular ATP demand [4]. More broadly, energy dissipation through proton leakage is known to be a method of metabolic control in bacteria [28,29]. To clarify the effect of decreasing of P/O ratio in the ME-model, we added a reaction in the model representing proton leakage (Methods). As a result, we see the Y-maximized solution with decreased P/O ratios have higher q glc , higher q ac , and lower Y at a given μ (Fig 4). Finally, experiments have shown that knocking out gnd leads to increased q glc and q ac and decreased Y with little change in μ [30]. The ME-model also predicts that gnd knockout mutants ("gnd knockout simulation" Methods) will have increased q glc , q ac and decreased Y (Fig 4). Since the ALE experiments do not introduce leaky proton or knock out any genes, it is also possible that multiple mechanisms working together, where the ME-model points to the systemic mechanisms for this fundamental second-order tradeoff. The exact pathways involved can be determined in future experiments.
Alternative explanations of the rate-yield tradeoff have been proposed, including membrane [31,32] and cytosolic crowding [33,34]. It is difficult to rule out these alternative constraints on cell growth, and it may be that multiple constraints operate at the same time. However, it is encouraging to see that the ME-model can explain the complex relationship between μ, Y, q ac , and q glc with only metabolic and proteome allocation constraints. In the future, it will be possible to extend ME-models with additional constraints. For example, it has been proposed that the unmodeled protein fraction(UPF) is growth-rate dependent, and thus existing proteome allocation models with fixed UPF are inaccurate [34]. If this is indeed the case, then SSME-and ME-models with cytosolic crowding constraints can be developed to fully represent the interplay between crowding, proteome allocation, and pathway selection.

Discussion
The E. coli ME-model provides a mechanistic and predictive model of rate-yield tradeoffs. It successfully reconciles several experimental data sets: i) uptake titration at low growth [8], ii) batch culture at higher growth rates [5], and iii) ALE endpoint strains (this study). These data sets, when analyzed with the ME-model, show the existence of a two-dimensional rate-yield tradeoff. The first dimension ("d1") rate-yield tradeoff is μ-Y tradeoff and the second dimension ("d2") is q glc -Y tradeoff.
From a mathematical perspective, one can describe the observed tradeoffs as correlations between any pair of the four variables μ, Y, q glc , and q ac . The two particular dimensions of the tradeoff that we describe, μ-Y ("d1") and q glc -Y ("d2"), are motivated by two different trends in our physiological observations. First, the previously-reported strong linear correlation between μ and Y [5] occurs for isogenic cultures under carbon limitation. The second dimension "d2" appears when comparing laboratory evolution endpoint strains, where q glc , q ac , and Y are observed to vary at a fixed μ, with a linear relationship in q glc -q ac and a corresponding inverse proportional relationship in q glc -Y. This two-dimensional tradeoff cannot be deciphered from simpler intuitive models, but it can be derived from the comprehensive set of metabolic and gene expression pathways represented by the ME-model. Furthermore, this study employed a multi-scale modeling approach where a smallscale model was used to guide parameter estimation in the genome-scale ME-model. This Laboratory evolution reveals a two-dimensional rate-yield tradeoff in microbial metabolism approach-which has been termed Tunable Resolution (TR) modeling [35]-was essential to the success of the study, and we expect that both small-scale and genome-scale models will continue to play an important role in understanding the genotype-phenotype relationship.
The two-dimensional rate-yield tradeoff appears as a result of ALE selection for μ when alternative pathway selection strategies achieve the same growth rate. Proton leakage and alternative ETC pathway selection are plausible mechanisms for modifying the P/O ratio and creating the q glc -Y tradeoff. In addition, the flux ratio between glycolysis (GAPD, gapA) and the pentose phosphate pathway (GND, gnd) might play a significant role in the q glc -Y tradeoff. Those mechanisms can be tested experimentally. Finally, revealing the underlying regulation would be of great interest for establishing a deeper understanding of rate-yield tradeoffs. Combining ME-models with known regulatory mechanisms to explain cellular choices would achieve a long-standing goal in systems biology [36].

Materials and methods
Phenotypic data including μ, q glc , q ac , and excretion rates of other metabolic byproducts were collected for ALE endpoint strains ("1 Phenotypic characterization of E. coli strains" in S1 Appendix). In addition, 13 C fluxes were measured from two of the strains with different growth rate and different glucose yield ("2 13 C metabolic flux analysis" in S1 Appendix). Reference data points of rate-yield, growth-acetate relations of wild-type MG1655 and NCM3722 E. coli strains were collected from published studies [5,8]. The coarse-grained proteome allocation model from [5] was reformulated as a small-scale ME-model (SSME-model, detail in "5 SSME-model parameters derivation" in S1 Appendix) and implemented by the COBRAme framework [16]. The genome-scale model iJL1678-ME was modified to fit experimental data by modifying the k eff s (enzyme turnover rate) of TCA cycle reactions, blocking target reactions, and modifying UPF (unmodeled protein fraction), GAM (growth associate maintenance energy), and NGAM (non-growth associate maintenance energy) ("8 Experimental data fitting" in S1 Appendix). Solution spaces were generated using flux balance analysis (incorporated in COBRAme) in the ME-model ("7 Solution space of the ME-model" in S1 Appendix).
To determine the effect of modifying P/O ratio on ME-model solution spaces, a reaction representing proton leakage was added to the ME-model ("10 P/O ratio manipulation" in S1 Appendix). The effect of the gnd knockout was demonstrated by blocking the reaction GND in ME-model simulations ("11 gnd knockout simulation" in S1 Appendix).
Supporting information S1 Appendix. Detailed introduction and discussions of materials and methods. (PDF) S1 Table. Parameters comparison between two coarse-grained models. Comparison between the coarse-grained proteome allocation model [5] and SSME-model. The derivation in detail is shown in "5 SSME-model parameter derivation" in S1 Appendix. (PDF) S2 Table. ME-model parameters. Global parameter selection in iJL1678-ME model to fit the μ-Y, μ-q ac data as in Fig 2C and 2D. (PDF) S3 Table. iJL1678-ME model modification (blocked reactions). Reactions that need to be turned off in the model to get quantitative fit of μ-Y, μ-q ac data as in Fig 2C and 2D.  [8]. (XLSX) S11 Table. 13 C metabolic flux analysis data. Metabolic fluxes distribution of the highest μ strain and highest Y strain among the ALE endpoint strains.
• Tab "Reactions" 13 C MFA model and carbon mapping network.
• Tab "MS_data" Measured mass distribution vectors (MDVs) by LC-MS/MS and their associated carbon mappings used for MFA calculations.
• Tab "Flux_data" Measured uptake and secretion rates by HPLC. In the genome-scale MEmodel, some TCA cycle reactions appeared as bp1 reactions, but, because they belong to the major respiration pathway of the cell, we will decreased their k eff s rather than blocking them entirely.
(TIF)  Table) that are blocked where maximum q ac s in high μ get lower, where the upper edge of the yellow region is below the upper edge of the pink region. The activation of one of these 24 reactions thus corresponding to higher q ac with lower Y. (B) 11 target reactions (S5 Table) corresponding to lower q ac with lower Y, blocking those reactions will get the minimum q ac (lower edge of the yellow region) closed to the Y-maximized q ac solution. (C) The method of picking reactions to block: Looking for the reactions that are not activated in the yield-maximized solution but activated at the maximal and minimal of the μq ac solution space, where the principal is to keep the Y-maximized solutions unchanged.  (1) and (3) expand the solution space to low-q ac at high μ. (C) Reactions (2) and (4) expand the solution space to high-q ac across all μ. (D) Model reactions that are added in the SSME-model for expanding the original solution space, all those reactions are guaranteed not be activated in the Y-maximized solutions so that the Y-optimal solution remains the same to fit data from [5]. Reaction (1) corresponds to the reactions that would generate products other than acetate such as pyruvate excretion, lactate excretion, etc. Reaction (2) is representative to the reactions that would generate other products, but at the same time generating acetate, such as pyruvate formate lyase (PFL), which produce formate and acetyl-CoA (precursor of acetate) from pyruvate. Reaction (3) and (4) could both be referred from the futile cycle in energy production and consumption, where (3) are the reactions that are less efficient than the optimal pathway, such as the alternative reactions in ETC which are less efficient in transporting electrons, while (4) are the reactions that would waste more energy in the same growth comparing to the optimal state, such as the reactions that would cause proton leakage. (TIF)