Introduction

A prominent feature of cancer metabolism is that tumor cells excrete large quantities of fermentation products in the presence of sufficient oxygen (Hanahan and Weinberg, 2011; Liberti and Locasale, 2016; Vander Heiden et al., 2009). This process, discovered by Otto Warburg in the 1920s (Warburg et al., 1924) and known as the Warburg effect, aerobic glycolysis or overflow metabolism, (Basan et al., 2015; Hanahan and Weinberg, 2011; Liberti and Locasale, 2016; Vander Heiden et al., 2009) is ubiquitous for fast-proliferating cells across a broad spectrum of organisms (Vander Heiden et al., 2009), ranging from bacteria (Basan et al., 2015; Holms, 1996; Meyer et al., 1984; Nanchen et al., 2006; Neidhardt et al., 1990) and fungi (De Deken, 1966) to mammalian cells (Hanahan and Weinberg, 2011; Liberti and Locasale, 2016; Vander Heiden et al., 2009). For microbes, cells use standard respiration when nutrients are scarce, while they use the counterintuitive aerobic glycolysis when nutrients are adequate, just analogous to normal tissues and cancer cells, respectively (Vander Heiden et al., 2009). Over the past century, especially with extensive studies in the last two decades (Liberti and Locasale, 2016), various rationales for overflow metabolism have been proposed (Basan et al., 2015; Chen and Nielsen, 2019; Majewski and Domach, 1990; Molenaar et al., 2009; Niebel et al., 2019; Pfeiffer et al., 2001; Shlomi et al., 2011; Vander Heiden et al., 2009; Varma and Palsson, 1994; Vazquez et al., 2010; Vazquez and Oltvai, 2016). In particular, Basan et al. (Basan et al., 2015) provided a systematic characterization of this process, including various types of perturbations in experiments. However, the origin and function of overflow metabolism still remain unclear (DeBerardinis and Chandel, 2020; Hanahan and Weinberg, 2011; Liberti and Locasale, 2016; Vander Heiden et al., 2009).

Why have microbes and cancer cells evolved to possess the seemingly wasteful aerobic glycolysis strategy? For unicellular organisms, there is evolutionary pressure (Vander Heiden et al., 2009) to optimize cellular resources for rapid growth (Dekel and Alon, 2005; Edwards et al., 2001; Hui et al., 2015; Li et al., 2018; Scott et al., 2010; Towbin et al., 2017; Wang et al., 2019; You et al., 2013). In particular, it has been shown that cells allocate protein resources for optimal growth (Hui et al., 2015; Scott et al., 2010; Wang et al., 2019; You et al., 2013), and the most efficient protein allocation corresponds to elementary flux mode (Müller et al., 2014; Wortel et al., 2014). In this study, we extend these approaches in a heterogeneous framework to address the puzzle of aerobic glycolysis. We use Escherichia coli as a typical example and show that overflow metabolism can be understood from optimal protein allocation combined with the heterogeneity in enzyme catalytic rates. The optimal growth strategy varies between respiration and fermentation depending on the concentration and type of the nutrient, and the combination with cell heterogeneity results in the standard picture (Basan et al., 2015; Holms, 1996; Meyer et al., 1984; Nanchen et al., 2006) of overflow metabolism. Our model quantitatively illustrates the growth rate dependence of fermentation/respiration flux and enzyme allocation under various types of perturbations, and can be used to explain the Warburg effect in tumor cells.

Results

Coarse-grained model

Based on the topology of the metabolic network (Neidhardt et al., 1990; Nelson et al., 2008) (Fig. 1A) , we classify the carbon sources that enter from the upper part of glycolysis into Group A (Wang et al., 2019), and the precursors of biomass components (such as amino acids) into five pools (see Appendix 1.2 for details): a1 (entry point: G6P/F6P), a2 (entry point: GA3P/3PG/PEP), b (entry point: pyruvate/acetyl-CoA), c (entry point: a-ketoglutarate), and d (entry point: oxaloacetate). Pools a1 and a2 are also combined as Pool a due to joint synthesis of precursors. Then, the metabolic network for Group A carbon source utilization (Fig. 1A) is coarse-grained into a model shown in Fig. 1B (see Appendix 2.1 for details) , where node A represents an arbitrary carbon source of Group A. Evidently, Fig. 1B is topologically identical to Fig. 1A. Each coarse-grained arrow in Fig. 1B carries a stoichiometric flux Ji, which delivers carbon flux and may consume or produce energy (e.g., J₁, J_a1, see Figs. 1A-B and Appendix-fig. 1A).

Fig 1.

In fact, the stoichiometric flux J_i scales with the cell population. For comparison with experiments, we define the normalized flux , which can be regarded as the flux per unit of biomass (the superscript “(N)” stands for normalized; see Appendix 1.3–1.4 for details). Here, M_carbon represents the carbon mass of the cell population, and m₀ is the weighted average carbon mass of metabolite molecules at the entrance of precursor pools (see Eq. S17). Then, the cell growth rate λ can be represented by (see Appendix 1.4), and the normalized fluxes of respiration and fermentation are and , respectively. In practice, each is characterized by two quantities: the proteomic mass fraction ϕ_i of the enzyme dedicated to carrying the flux and the substrate quality κ_i, such that . We take the Michaelis-Menten form for the enzyme kinetics, and then . (see Eq. S12 and Appendix 1.4 for details) , where [S_i] is the concentration of substrate S_i. For each intermediate node and reaction along the pathway (e.g., node M₁ in J_a1), the substrate quality κ_i can be approximated as a constant (see Appendix 1.5): , where [S_i] ≥ K_i generally holds true in bacteria (Bennett et al., 2009; Park et al., 2016). However, κ_A is a variable that depends on the concentration and nutrient type of a Group A carbon source (see Eq. S27).

Generally, there are three distinct destinies of a Group A carbon source in the metabolic network (Appendix-fig. 1C-E) : fermentation, respiration, and biomass generation. Each draws a proteome fraction of ϕ_f, ϕ_r and ϕ_BM. The net effect of the first two destinies is energy production, while the last one generates precursors of biomass accompanied by energy production. By applying the proteomic constraint(Scott et al., 2010) that there is a maximum fraction ϕ_max for proteome allocation (ϕ_max ≈ 0.48 (Scott et al., 2010)), we have:

In fact, Eq. 1 is equivalent to (see Appendix 2.1) , where ϕ_R stands for the mass fraction of the active ribosome-affiliated proteins. In cell proliferation, ribosomes constitute the majority of enzymes for protein synthesis (Neidhardt et al., 1990; Nelson et al., 2008), while other biomass components such as RNA are optimally produced (Kostinski and Reuveni, 2020) following the growth rate set by protein synthesis. Thus, λ = ϕ_R · κ_t, where κ_t is a parameter determined by the translation rate (Scott et al., 2010) (see Appendix 1.1 for details), which can be approximated as a constant in the growth rate range of concern (Dai et al., 2016).

For balanced cell growth, the energy demand is generally proportional to the biomass production rate. Thus, the normalized energy production rate is proportional to the growth rate λ (see Appendix 2.1 for details):

where η_E is the energy coefficient. By converting all the energy currencies into ATPs, the normalized energy fluxes of respiration and fermentation are and , where and are the stoichiometric coefficients of ATP production per glucose in each pathway (see Appendix-fig. 1C-E and Appendix 2.1 for details). The denominator coefficient (“2”) is derived from the stoichiometry of the coarse-grained reaction M₁→M₂ (see Fig. 1A-B). Using flux balance analysis at each intermediate node (M_i, i = 1, …, 5) and precursor pool (Pool i, i = a1, a2, b, c, d), combined with the constraints of proteome allocation (Eq. 1) and energy demand (Eq. 2), we obtain the relations between energy fluxes and growth rate for a given nutrient condition with fixed κ_A (see Appendix 2.1 for details):

where φ is a constant coefficient mainly determined by η_E (see Eq. S33), and φ·λ represents the normalized energy demand other than the biomass pathway. The coefficients ψ, ε_r and ε_f are all functions of κ_A·. ψ^-1 is the proteome efficiency of the biomass pathway (see Eq. S32), with ψ^-1 = λ/ϕ_BM · ε_r and ε_f are the proteome energy efficiencies of the respiration and fermentation pathways, with and :

where both and are constants, with ≡ 1/κ₁ + 2/κ₂ + 2/κ₃ + 2/κ₄ and ≡ 1/κ₁ + 2/κ₂ + 2/κ₆ (see Appendix 2.1 for details).

Origin of overflow metabolism

The standard picture of overflow metabolism (Basan et al., 2015; Holms, 1996; Meyer et al., 1984; Nanchen et al., 2006) is exemplified by the experimental data (Basan et al., 2015) shown in Fig. 1C: the fermentation flux exhibits a threshold-analog dependence on the growth rate λ. It is well known that respiration is far more efficient than fermentation in terms of energy production per unit carbon (i.e., > ) (Nelson et al., 2008; Vander Heiden et al., 2009). Then, why do cells bother to use the seemingly wasteful fermentation pathway? We proceed to resolve these issues using the optimal protein allocation hypothesis (Scott et al., 2010; You et al., 2013).

For cell proliferation in a given nutrient with fixed κ_A , the values of ε_r, ε_f and ψ are determined (Eqs. 4 and S32). However, the growth rate λ can be influenced by protein allocation between ϕ_r and ϕ_f with the governing equation Eq. 3. If ε_r < ε_r, then . The optimal solution is ε_f = = 0, which means that the cell only uses respiration. Conversely, if ε_f > ε_r, ϕ_r = = 0 is optimal, and the cell only uses fermentation. In either case, the optimal choice between respiration and fermentation is determined by weighing the proteome energy efficiencies.

In practice, both ε_r and ε_f are functions of κ_A (Eq. 4), and therefore the optimal choice may vary depending on the nutrient conditions. In nutrient-poor conditions where κ_A ≪ and κ_A ≪ , the proteome energy efficiencies are reduced to ε_r ≈ · κ_A and ε_f ≈ · κ_A (see Eq. 4), where ε_r (κ_A)> ε_f (κ_A). In rich media, however, using the parameters of κ_i derived from in vivo/in vitro experimental data of Escherichia coli (Appendix-table 1, see also Appendix 6.1-6.2 and Appendix-table 2), we confirm that ε_r < ε_f with Eq. 4 (see also Eqs. S39–S40), where represents the substrate quality of glucose at saturated concentration (abbreviated as “ST” in the superscript). Indeed, recent experimental studies have validated that ε_r < ε_f for lactose (Basan et al., 2015). In Fig. 1E, we present the growth rate dependence of proteome energy efficiencies ε_r and ε_f in a three-dimensional (3D) form using the data from Appendix-table 1, where ε_r, ε_f and the growth rate λ all vary as functions of κ_A. Therefore, the ratio ε_f/ε_r (defined as Δ(κ_A)) is a monotonically increasing function of κ_A, and there exists a critical value of κ_A (denoted as , see Appendix 2.2 for details) satisfying . Below , the cell grows slowly and the optimal choice is respiration (with ε_f < ε_r, λ = ϕ_max ·(ψ + φ/ε_r)^-1), while above , the cell grows faster and the optimal growth strategy is fermentation (with ε_f > ε_r, λ = ϕ_max ·(ψ + φ/ε_f)^-1). The above analysis qualitatively explains the phenomenon of aerobic glycolysis.

For a quantitative understanding of overflow metabolism, let us first consider the homogeneous case, where all cells share identical biochemical parameters. For optimal protein allocation, the relation between fermentation flux and growth rate is = φ·φ·θ(φ − φ_C), where “θ” represents the Heaviside step function, and λ_C is the critical growth rate for (i.e., λ_C ≡ λ). Similarly, we can determine the growth rate dependence of respiration flux, which is = φ·λ·[1 − θ(λ − λ_C)]. These digital response outcomes are consistent with the numerical simulation findings of Molenaar et al. (Molenaar et al., 2009). However, they are clearly incompatible with the threshold-analog response in the standard picture of overflow metabolism (Basan et al., 2015; Holms, 1996; Meyer et al., 1984; Nanchen et al., 2006).

To address this issue, we take into account cell heterogeneity, which is ubiquitous in both microbes (Ackermann, 2015; Bagamery et al., 2020; Balaban et al., 2004; Nikolic et al., 2013; Solopova et al., 2014; Wallden et al., 2016) and tumor cells (Duraj et al., 2021; Hanahan and Weinberg, 2011; Hensley et al., 2016). For the Warburg effect or overflow metabolism of our concern, experimental studies have reported significant metabolic heterogeneity in the choice between respiration and fermentation within a cell population (Bagamery et al., 2020; Duraj et al., 2021; Hensley et al., 2016; Nikolic et al., 2013). Motivated by the fact that the turnover number (k_cat value) of a catalytic enzyme varies considerably between in vitro and in vivo measurements (Davidi et al., 2016; García - Contreras et al., 2012), we note that the concentrations of potassium and phosphate, which vary from cell to cell, have a significant impact on the k_cat values of the metabolic enzymes (García - Contreras et al., 2012). Therefore, in a cell population, there is a distribution of values for k_cat, which is commonly referred to as extrinsic noise (Elowitz et al., 2002). For simplicity, we assume that each k_cat value follows a Gaussian distribution. This gives the distributions of single-cell growth rate in various types of carbon sources (see Eqs. S155–S157, S163–S165), which are fully verified by recent experiments using isogenic Escherichia coli with single-cell resolution (Wallden et al., 2016) (Appendix-fig. 2B). Then, the critical growth rate λ_C should follow a Gaussian distribution in the cell population (see Appendix 7 for details), where μ_λC is approximated by the deterministic result of λ_C (Eq. S43). We assume that the coefficient of variation (CV) is σ_λC/μ_λC = 12%, or equivalently, the CV for the catalytic rate of each metabolic enzyme is 25%. Thus, we have the growth rate dependence of fermentation and respiration fluxes (see Appendix 2.3 for details):

where “erf” represents the error function. The fermentation flux exhibits a threshold-analog relation with the growth rate (the red curves in Figs. 1C-D, 2B-C and 3B, D, F), while the respiration flux (the blue curve in Fig. 1D) decreases with an increase in fermentation flux. In Fig. 1C-D, we see that the model results (see Eq. 5 and Appendix 8; parameters are set by the experimental data shown in Appendix-table S1) agree quantitatively with the experimental data of Escherichia coli (Basan et al., 2015; Holms, 1996). The fermentation flux was determined by the acetate secretion rate , and the respiration flux was deduced from the carbon dioxide flux (the superscript “(M)” represents the measurable flux in the unit of mM/OD600/h; see Appendix 8.1 for details). Therefore, with the integration of cell heterogeneity, our model of optimal protein allocation quantitatively explains overflow metabolism.

Fig 2.

Fig 3.

Testing the model through perturbations

To further test our model, we systematically investigate model predictions under various types of perturbations and compare them with the experimental data from existing studies (Basan et al., 2015; Holms, 1996) (see Appendices 3 and 4.1 for details).

First, we consider the proteomic perturbation by overexpression of useless proteins encoded by the Lacz gene (i.e., ϕ_Z perturbation) in Escherichia coli. The net effect of the ϕ_Z perturbation is that the maximum fraction of proteome available for resource allocation changes from ϕ_max to ϕ_max − ϕ_Z⁵, where ϕ_Z is the mass fraction of useless proteins. In a cell population, the critical growth rate λ_C(ϕ_Z) still follows a Gaussian distribution N(μ_λC (ϕ_Z), σ_λC (ϕ_Z)²), where the CV of λ_C (ϕ_Z) remains unchanged. Consequently, the growth rate dependence of fermentation flux changes into (see Appendix 3 for model perturbation results of respiration flux), where both the growth rate λ(κ_A, ϕ_Z) and the normalized fermentation flux (κ_A, ϕ_Z) are bivariate functions of κ_A and ϕ_Z (see Eqs. S49 and S56–S57). For each degree of LacZ expression (with fixed ϕ_Z), similar to the wild-type strains, the fermentation flux exhibits a threshold-analog response to growth rate as κ_A varies (Fig. 2C), which agrees quantitatively with experimental results (Basan et al., 2015). The shifts in the critical growth rate λ_C(ϕ_Z) are fully captured by μ_λC (ϕ_Z) = μ_λC(0)(1 − ϕ_Z/ϕ_max). In contrast, for nutrient conditions with each fixed κ_A, since the growth rate changes with ϕ_Z just like λ_C(ϕ_Z): λ(κ_A, ϕ_Z) = λ(κ_A, 0)(1 − ϕ_Z/ϕ_max), the fermentation flux is then proportional to the growth rate for the varying levels of Lacz expression: , where the slope is a monotonically increasing function of the substrate quality κ_A. These scaling relations are well validated by the experimental data (Basan et al., 2015) shown in Fig. 2B. Finally, in the case where both κ_A and ϕ_Z are free to vary, the growth rate dependence of fermentation flux presents a threshold-analog response surface in a 3D plot where ϕ_Z appears explicitly as the y-axis (Fig. 2A). Experimental data points (Basan et al., 2015) lie right on this surface (Fig. 2A), which are highly consistent with the model predictions.

Next, we study the influence of energy dissipation, which introduces an energy dissipation coefficient “w” in Eq. 2: = η_E · λ + w. Similarly, the critical growth rate λ_C (w) follows a Gaussian distribution N(μ_λC(w), σ_λC(w)²) in a cell population. The relation between the growth rate and fermentation flux can be characterized by (see Appendix 3.2 for details): . In Fig. 3A-B, we present a comparison between the model results and experimental data (Basan et al., 2015) in 3D and 2D plots, which demonstrate good agreement. A notable characteristic of energy dissipation, as distinguished from ϕ_Z perturbation, is that the fermentation flux increases despite a decrease in the growth rate when κ_A is fixed.

We proceed to analyze the impact of translation inhibition with different sub-lethal doses of chloramphenicol on Escherichia coli. This type of perturbation introduces an inhibition coefficient “l” in the translation rate, thus turning κ_t into κ_t/(l+1). Still, the critical growth rate λ_C(l) follows a Gaussian distribution N(μ_λC(l), σ_λC(l)²), and then, the growth rate dependence of fermentation flux is: (see Appendix 3.3 for details). In Appendix-fig. 2D-E, we see that the model predictions are generally consistent with the experimental data (Basan et al., 2015). However, there is a noticeable systematic discrepancy when the translation rate is small. Therefore, we turn to maintenance energy, which is tiny and generally negligible for bacteria over the growth rate range of concern (Basan et al., 2015; Neidhardt, 1996). Encouragingly, by assigning a very small value to the maintenance energy coefficient “w₀” (w₀ = 2.5 (h^-1)), the model results for the growth rate-fermentation flux relation agree quantitatively with experiments (Basan et al., 2015) (see Fig. 3C-D and Appendix 3.3 for details).

Finally, we consider the alteration of nutrient categories by switching to a non-Group A carbon source: pyruvate, which enters the metabolic network from the endpoint of glycolysis (Neidhardt et al., 1990; Nelson et al., 2008). The coarse-grained model for pyruvate utilization is shown in Fig. 3E (see also Fig. 1A), which shares identical precursor pools as those for Group A carbon sources (Fig. 1B), yet with several differences in the coarse-grained reactions. The growth rate dependencies of both the proteome energy efficiencies (Appendix-fig. 2H) and energy fluxes are qualitatively similar to those of the Group A carbon source utilization, while there are quantitative differences in the coarse-grained parameters (see Appendices 4.1 and 8 for details). Most notably, the critical growth rate and the ATP production per glucose in the fermentation pathway are noticeably smaller than λ_C and (for Group A sources), respectively. Consequently, the growth rate dependence of fermentation flux in pyruvate should present a quite different curve from that of Group A carbon sources (see Eqs. 5 and S105), which is fully validated by experimental results (Holms, 1996) (Fig. 3F).

Enzyme allocation under perturbations

As mentioned above, our coarse-grained model is topologically identical to the central metabolic network (Fig. 1A), and thus it can predict enzyme allocation for each gene in glycolysis and the TCA cycle (see Appendix-fig. 1B and Appendix-table 1) under various types of perturbations. In Fig. 1B, the intermediate nodes M₁, M₂, M₃, M₄, and M₅ represent G6P, PEP, acetyl-CoA, α-ketoglutarate, and oxaloacetate, respectively. Then, ϕ₁ and ϕ₂ correspond to enzymes of glycolysis (or at the junction of glycolysis and the TCA cycle), while ϕ₃ and ϕ₄ correspond to enzymes in the TCA cycle.

We first consider enzyme allocation under carbon limitation by varying the nutrient type and concentration of a Group A carbon source (i.e., κ_A perturbation). In fact, this has been extensively studied in more simplified models (Hui et al., 2015; You et al., 2013), where the growth rate dependence of enzyme allocation under κ_A perturbation was generally considered to be a C-line response (Hui et al., 2015; You et al., 2013), i.e., the genes responsible for digesting carbon compounds show a linear increase in gene expression as the growth rate decreases (Hui et al., 2015; You et al., 2013). However, when it comes to enzymes catalyzing reactions between intermediate nodes, we collected experimental data from existing studies (Hui et al., 2015) and found that the enzymes in glycolysis exhibit a completely different response pattern compared to those in the TCA cycle (Appendix-fig. 3A-B). This discrepancy cannot be explained by the C-line response. To address this issue, we apply the coarse-grained model described above (Fig. 1B) to calculate the growth rate dependence of enzyme allocation for each ϕ_i (i = 1, 2, 3, 4) using the model settings for wild-type strains, where no fitting parameters are involved in determining the shape (see Eqs. S118–S119 and Appendix 8). In Figs. 4A-B and Appendix-fig. 3C-D, we see that the model predictions overall match with the experimental data (Hui et al., 2015) for representative genes from either glycolysis or the TCA cycle, and maintenance energy (with w₀ = 2.5 (h^-1)) has a negligible effect on this process. Still, there are minor discrepancies that arise from the basal expression of metabolic genes, which may be attributed to the fact that our model deals with relatively stable growth conditions while microbes need to be prepared for fluctuating environments (Basan et al., 2020; Kussell and Leibler, 2005; Mori et al., 2017).

Fig 4.

We proceed to analyze the influence of ϕ_Z perturbation and energy dissipation. In both cases, our model predicts a linear response to the growth rate reduction for all genes in either glycolysis or the TCA cycle (see Appendix 5.2–5.3 for details). For ϕ_Z perturbation, all predicted slopes are positive, and there are no fitting parameters involved (Eqs. S120–S121). In Figs. 4C-D and Appendix-fig. 3E-J, we show that our model quantitatively illustrates the experimental data (Basan et al., 2015) for representative genes in the central metabolic network, and there is a better agreement with experiments (Basan et al., 2015) by incorporating the maintenance energy (with w₀ = 2.5 (h^-1)). For energy dissipation, however, the predicted slopes of the enzymes corresponding to ϕ₄ are surely negative, and there is a constraint that the slope signs of the enzymes corresponding to the same ϕ_i (i = 1, 2, 3) should be the same. In Appendix-fig. 3K-N, we see that the model results (Eqs. S127 and S123) are consistent with experiments (Basan et al., 2015).

Discussion

The phenomenon of overflow metabolism, or the Warburg effect, has been a long-standing puzzle in cell metabolism. Although many rationales have been proposed (Basan et al., 2015; Chen and Nielsen, 2019; Majewski and Domach, 1990; Molenaar et al., 2009; Niebel et al., 2019; Pfeiffer et al., 2001; Shlomi et al., 2011; Vander Heiden et al., 2009; Varma and Palsson, 1994; Vazquez et al., 2010; Vazquez and Oltvai, 2016) over the past century, the origin and function of this phenomenon remain unclear (DeBerardinis and Chandel, 2020; Hanahan and Weinberg, 2011; Liberti and Locasale, 2016; Vander Heiden et al., 2009). In this study, we use Escherichia coli as a typical example and demonstrate that overflow metabolism can be understood through optimal protein allocation combined with cell heterogeneity. In nutrient-poor conditions, the proteome energy efficiency of respiration is higher than that of fermentation (Fig. 1E), and thus the cell uses respiration to optimize growth. In rich media, however, the proteome energy efficiency of fermentation increases faster and is higher than that of respiration (Fig. 1E), leading the cell to use fermentation to accelerate growth. In further combination with cell heterogeneity in enzyme catalytic rates (Davidi et al., 2016; García - Contreras et al., 2012), our model quantitatively illustrates the thresholdanalog response (Basan et al., 2015; Holms, 1996) in overflow metabolism (Fig. 1C).

Cell heterogeneity is crucial for the threshold-analog response in overflow metabolism. In the homogeneous case, the optimal solution is a digital response (Eq. S44) that corresponds to an elementary flux mode (Müller et al., 2014; Wortel et al., 2014) and agrees with the numerical study of Molenaar et al. (Molenaar et al., 2009). However, this digital response is incompatible with the standard picture of overflow metabolism (Basan et al., 2015; Holms, 1996; Meyer et al., 1984; Nanchen et al., 2006). By incorporating heterogeneity in enzyme catalytic rates (Davidi et al., 2016; García - Contreras et al., 2012), the critical growth rate (i.e., threshold) changes from a single value into a Gaussian distribution (Eq. 45, see Appendix 7 for details; see also Appendix-fig. 4) for a cell population, thus turning a digital response into the threshold-analog response in overflow metabolism (Fig. 1C). Our model results relying on cell heterogeneity are fully validated by the observed distributions of single-cell growth rate (Wallden et al., 2016) (Appendix-fig. 2B) and experiments with various types of perturbations (Basan et al., 2015; Holms, 1996; Hui et al., 2015), both for acetate secretion patterns and gene expression in the central metabolic network (Figs. 2–4 and Appendix-figs. 2D-E and 3).

Finally, our model can be broadly used to address heterogeneity-related challenges in metabolism on a quantitative basis, including the Crabtree effect in yeast (Bagamery et al., 2020; De Deken, 1966), the Warburg effect in cancer (Duraj et al., 2021; Hanahan and Weinberg, 2011; Liberti and Locasale, 2016; Vander Heiden et al., 2009) (see Appendix 6.4 for an explanation of the Warburg effect), and the heterogeneous metabolic strategies of cells in various types of environments (Bagamery et al., 2020; Duraj et al., 2021; Escalante-Chong et al., 2015; Hensley et al., 2016; Liu et al., 2015; Solopova et al., 2014; Wang et al., 2019).

Acknowledgements

The author thanks Chao Tang, Qi Ouyang, Yang-Yu Liu and Kang Xia for helpful discussions. This work was supported by National Natural Science Foundation of China (Grant No.12004443), Guangzhou Municipal Innovation Fund (Grant No.202102020284) and the Hundred Talents Program of Sun Yat-sen University.

Author contributions

X.W. conceived and designed the project, developed the model, carried out the analytical/ numerical calculations, and wrote the paper.

Data, Materials, and Software Availability

All study data are included in the article and/or appendices.

Appendix 1 Model framework

Appendix 1.1 Proteome partition

Here we adopt the proteome partition framework similar to that introduced by Scott et al.(Scott et al., 2010). All proteins in a cell are classified into three classes: the fixed portion Q-class, the active ribosome-affiliated R-class, and the remaining catabolic/anabolic enzymes C-class. Each proteome class has a mass (i = Q, R, C) and mass fraction ϕ_i, where ϕ_Q is a constant, and we define ϕ_max ≡ 1-ϕ_Q. In the exponential growth phase, the ribosome allocation for protein synthesis of each class is f_i, with f_Q + f_R + f_C = 1.

To analyze cell growth optimization, we first consider the homogeneous case where all cells share identical biochemical parameters and simplify the mass accumulation of a cell population into a big cell. Essentially, this approximation would not influence the value of growth rate λ. For bacteria, the protein turnover rate is negligible, and thus the mass accumulation of each class follows:

where m_AA stands for the average molecular weight of amino acids, k_T is the translation rate, is the number of ribosomes, m_R is the protein mass of a single ribosome, and is the protein mass of the ribosomes with (Neidhardt, 1996; Scott et al., 2010). For a specific stable nutrient environment, f_R and k_T are temporal invariants. Then,

where , and the protein mass of the cell population follows:

Over a long period in the exponential growth phase (t → +∞), ϕ_i = f_i (i = Q, R, C) and

where .

Appendix 1.2 Precursor pools

Based on the entry point of the metabolic network, we classify the precursors of biomass components into five pools (Fig. 1A-B): a1 (entry point: G6P/F6P), a2 (entry point: GA3P/3PG/PEP), b (entry point: pyruvate/acetyl-CoA), c (entry point: α-ketoglutarate) and d (entry point: oxaloacetate). For bacteria, these five pools draw roughly r_a1 = 24 %, r_a2 = 24 %, r_b = 28 %, r_c = 12 % and r_d = 12 % of the carbon flux(Nelson et al., 2008; Wang et al., 2019). There are overlapping components between Pools a1 and a2 due to the joint synthesis of some precursors, thus we also use Pool a to represent Pools a1-a2 in the descriptions.

Appendix 1.3 Stoichiometric flux

We consider the following biochemical reaction between substrate S_i and enzyme E_i:

where a_i, d_i and are the reaction parameters, S_i+1 is the product, b_i and c_i are the stoichiometric coefficients. For most of the reactions in the central metabolism, b_i = 1 and c_i = 0. The reaction rate follows Michaelis–Menten kinetics(Nelson et al., 2008):

where , [E_i] and [S_i] are the concentrations of enzymeE_i and substrate E_i, respectively. For this reaction (Eq. S5), and . In the cell populations (the “big cell”), suppose that the cell volume is V_cell, then the stoichiometric flux of the reaction is:

The copy number of enzyme E_i is with a total weight of , where is the molecular weight of E_i. By defining the enzyme cost of an E_i molecule as , where m₀ is a unit mass, then the cost of all E_i molecules is (Wang et al., 2019). By further defining , then,

The mass fraction of E_i is , and thus

Appendix 1.4 Carbon flux and cell growth rate

To clarify the relation between the stoichiometric flux J_i and growth rate λ, we consider the carbon flux in the biomass production. The carbon mass of the cell population (the “big cell”) is given by M_carbon = M_protein · r_carbon/r_protein, where r_carbon and r_protein represent the mass fraction of carbon and protein within a cell. In the exponential growth phase, the carbon flux of the biomass production is given by:

where m_carbon is the mass of a carbon atom. In fact, the carbon mass flux per stoichiometry varies depending on the entry point of the precursor pool. Taking Pool b as an example, there are three carbon atoms in a molecule of the entry point metabolite (i.e., pyruvate). Assuming that carbon atoms are conserved from pyruvate to Pool b, then the carbon flux of Pool b is given by , where J_b is the stoichiometric flux from pyruvate to Pool b (Fig. 1A-B)and stands for the carbon number of a pyruvate molecule. Combining with Eq. S10 and noting that , we get . Similarly, for each precursor pool, we have:

where the subscript “EP_i” represents the entry point of Pool i, and stands for the number of carbon atoms in a molecule of the entry point metabolite.

For each substrate in intermediate steps of the metabolic network, we define κ_i as the substrate quality:

and for each precursor pool, we define:

Combining Eqs. S8, S9 and S11, we have

Then, we define the normalized flux, which can be regarded as the flux per unit of biomass:

where the superscript “(N)” stands for normalized. Combined with Eqs. S8, S9 and S12, we have:

Since , by setting

then,

and we have , and

Appendix 1.5 Intermediate nodes

In a metabolic network, the metabolites between the carbon source and precursor pools are the intermediate nodes. As specified in Wang et al.(Wang et al., 2019), to optimize cell growth rate, the substrate of each intermediate node is nearly saturated, and thus

The real cases could be more complicated because of other metabolic regulations. Recent quantitative studies(Bennett et al., 2009; Park et al., 2016)have shown that at least in E. coli, for most of the substrate-enzyme pairs [S_i] > K_i, which implies .

Appendix 2 Model and analysis

Appendix 2.1 Coarse-grained model

In the coarse-grained model shown in Fig. 1B, node A represents an arbitrary carbon source of Group A(Wang et al., 2019), which joins at the upper part of glycolysis. Nodes M1, M2, M3, M4, and M5 stand for G6P, PEP, acetyl-CoA, α-ketoglutarate, and oxaloacetate, respectively. In the analysis of carbon supply into precursor pools, we lump sum G6P/F6P as M1, GA3P/3PG/PEP as M2, and pyruvate/acetyl-CoA as M3 for approximation. For the biochemical reactions, each follows Eq. S5 with b_i = 1 except that M1→2M2 and M3+M5→M4. Basically, there are 3 possible destinies of a Group A carbon source (e.g., glucose, see Appendix-fig. 1C-E): energy contributions in the fermentation and respiration pathways (Appendix-fig.1C-D), or biomass components accompanied by energy production in the biomass pathway (Appendix-fig.1E).

By applying flux balance to the stoichiometric fluxes and combining with Eq. S8, we have:

Obviously, the stoichiometric fluxes of respiration J_r and fermentation J_f (Appendix-fig. 1C-D) are:

We further assume that the carbon atoms are conserved from each entry point metabolite to the precursor pool, and then,

In terms of energy production for the relevant reactions, for convenience, we convert all the energy currencies into ATPs, namely, NADH → 2ATP(Neidhardt et al., 1990), NADPH → 2ATP(Neidhardt et al., 1990; Sauer et al., 2004), FADH₂ → 1ATP(Neidhardt et al., 1990). Then, we have

where J_E represents the stoichiometric flux of ATPs, and β_i is the stoichiometric coefficient with β₁ = 4, β₂ = 3, β₃ = 2, β₄ = 6, and β_a1 = 4 (Neidhardt et al., 1990; Sauer et al., 2004). Generally, the energy demand is proportional to the carbon flux infused into biomass production, thus,

where r_E is the ratio and also a constant.

By applying the substitutions specified in Eqs. S9, S12, S14–S18, combined with Eqs. S4, S10, S21–S25, and the constraint of proteome resource allocation ϕ_R + ϕ_C = ϕ_max, we have:

where . Here, for each intermediate node, κ_i follows Eq. S20, which can be approximated as a constant. The substrate quality of the Group A carbon source κ_A varies with the identity and concentration of the Group A carbon source:

which is determined externally by the culture condition. From Eq. S26, all ϕ_i can be expressed by , , and λ:

By substituting Eq. S28 into Eq. S26, we have:

Here, and stand for the normalized energy fluxes of respiration and fermentation, with

The coefficients ε_r and ε_f represent the proteome energy efficiencies of the respiration and fermentation pathways (Appendix-fig.1C-D), respectively, with

ψ^-1 is the proteome efficiency of biomass pathway (Appendix-fig. 1E), with

φ is the energy demand coefficient (a constant), with

and φ · λ stands for the normalized energy demand other than the accompanying energy production from the biomass pathway.

Appendix 2.2 The reason for overflow metabolism

Microbes optimizetheir growth rate to survive in the evolutionary process(Vander Heiden et al., 2009).Basically, this also applies to tumor cells, which proliferate rapidly ignoring signals of growth restriction(Vander Heiden et al., 2009). To optimize cell growth, we first consider the best strategy for a single cell. The coarse-grained model is summarized in Eq. S26 and further simplified into Eq. S29. Here, ε_r, ε_f and ψ are functions of κ_A (see Eqs. S31, S32), so we also denote them as ε_r(κ_A), ε_f(κ_A), ψ(κ_A). Apparently, the fluxes of both respiration and fermentation take non-negative values, i.e., , and all the coefficients are positive: ε_r(κ_A), ε_f(κ_A), ψ(κ_A), φ > 0.

Thus, if ε_r > ε_f, then . Obviously, the optimal solution is:

Similarly, if ε_f > ε_f, then the optimal solution is:

In both cases, the growth rate λ takes the maximum value for a given nutrient condition (i.e., given κ_A):

So, why do microbes use the wasteful fermentation pathway when the growth rate is large under aerobic conditions? An intuitive speculation is that the fermentation pathway is more efficient in terms of the proteome energy efficiency, i.e., ε_f > ε_r. If so, then why do microbes still use the normal respiration pathway when the growth rate is small? The answer lies in that both ε_r(κ_A) and ε_f(κ_A) are not constants, but are dependent on nutrient conditions. In Eq. S31, when κ_A is small, just consider the extreme case of κ_A → 0, and then

Since β₃ + β₄ ≫ β₆, clearly,

Combined with Eq. S36, thus cells would certainly use the respiration pathway when the growth rate is very small. Meanwhile, suppose that is the maximum value of κ_A available across different Group A carbon sources, and if there exists a κ_A (with ) satisfying ε_r(κ_A) < ε_f(κ_A), specifically,

then Δ(κ_A) ≡ ε_f(κ_A)/ε_r(κ_A) is a monotonously increasing function of κ_A. Thus,

and cells would use the fermentation pathway when the growth rate is large.

In practice, experimental studies(Basan et al., 2015)in E. coli have reported that the proteome energy efficiency in fermentation is higher than that in respiration when the Group A carbon source is lactose at saturated concentration(Molenaar et al., 2009), i.e., . Here, represents the substrate quality of lactose and the superscript “(ST)” signifies saturated concentration. In fact, E. coli grows much faster in G6p than lactose(Basan et al., 2015), thus, . The above example verifies that Eq. S40 holds for E. coli. In fact, more recent studies(Chen and Nielsen, 2019)supported that Eq. S40 holds generically for many microbial species. From the theoretical side, we can verify Eq. S39 and thus Eq. S40 using Eq. S20 combined with the in vivo/in vitro biochemical parameters obtained from experimental data (see Appendix-tables 1–2). For example, it is easy to check that with this method (see Appendix 8.2), which also confirms the validity of Eqs. S39–S40.

Now that Eqs. S38–S40 are all valid, then there exists a critical value of κ_A (denoted as ) satisfying . Thus,

Combined with Eq. S31, we have:

By substituting Eq. S42 into Eqs. S31, S32 and S36, we obtain the expressions for , and the critical growth rate at the transition point (i.e., ):

where ε_r/f represents either ε_r or ε_f. In Fig. 1E, we show the dependencies of ε_r(κ_A), ε_f(κ_A) and λ(κ_A) on κ_A in a 3-dimensional form, as κ_A changes.

Appendix 2.3 The relation betweenrespiration/fermentationfluxes and growth rate

We proceed to study the relation between the respiration/fermentation flux and the cell growth rate. From Eqs. S16 and S30, we see that the stoichiometric fluxes J_r, J_f, the normalized fluxes , and the normalized energy fluxes , are all interconvertible. For convenience, we first analyze the relations between , , and λ under growth rate optimization. In fact, all these terms are merely functions of κ_A (see Eqs. S34–S36), which is determined by the nutrient condition (Eq. S27).

In the homogeneous case, i.e., all microbes share identical biochemical parameters, as λ(κ_A) increases with κ_A, show up abruptly and vanish simultaneously as κ_A right exceeds (Fig. 1E, see also Eqs. S34–S35, S41).Combining Eqs. S34–S36 and S43, we have:

where “θ” stands for the Heaviside step function. Defining , and then, [0, λ_max] is the relevant range of the x axis. In fact, the digital responses in Eq. S44 are consistent with the numerical simulation results in Molenaar et al.(Molenaar et al., 2009). However, these results are incompatible with the threshold-analog response in the standard picture of overflow metabolism(Basan et al., 2015; Holms, 1996).

In practice, the values of can be greatly influenced by the concentrations of potassium and phosphate(García‐Contreras et al., 2012), which vary from cell to cell. Then, there is a distribution of the values for among the cell populations, commonly we call it extrinsic noise(Elowitz et al., 2002). For convenience, we assume that each (and thus κ_i) follows a Gaussian distribution with the coefficient of variation (CV) taken as 25%. Then, the distribution of λ_C can be approximated by a Gaussian distribution (see Appendix 7.1):

where and stand for the mean and standard deviation of λ_C, with the CV calculated to be 12% (see Appendix 8.2 for details). Note that λ is κ_A dependent, while λ_C is independent of κ_A. Thus, given the growth rate of microbes in a culturing medium (e.g., in a chemostat), the normalized energy fluxes are:

where “erf” represents the error function. In practice, given a culturing medium, there is also a probability distribution for the growth rate (Appendix-fig.2B, see also Eq. S157). For approximation, in plotting the flux-growth rate relations, we use the deterministic (noise-free) value of the growth rate as the proxy. To compare with experiments, basically, we are comparing the normalized fluxes , (see Appendix 8.1 for details). Combining Eqs. S30 and S46, we get:

In Fig. 1C-D, we see that Eq. S47 quantitatively illustrates the experimental data(Basan et al., 2015), where the model parameters were obtained using the biochemical data for the catalytic enzymes(see Appendix-table1 for details).

Appendix 3 Model perturbations

Appendix 3.1 Overexpression of useless proteins

Here we consider the case of overexpression of the protein encoded by Lacz gene in E. coli. Effectively, this limits the proteome by altering ϕ_max:

where ϕ_Z stands for the fraction of useless proteins, which is controllable in experiments. Then, the growth rate changes into a bivariate function of κ_A and ϕ_Z:

and thus,

Obviously, is still a constant (following Eq. S42), while and become functions of ϕ_Z:

In the homogeneous case, and follow:

Combined with Eqs. S50–S51, we have:

To compare with experiments, we assume that each and κ_i follow the extrinsic noise with a CV of 25% specified in Appendix 2.3, and neglect the noise on ϕ_Z and ϕ_max. Combining Eqs. S45 and S51, we see that λ_C(ϕ_Z) approximately follows a Gaussian distribution:

where and represent the mean and standard deviation of λ_C(ϕ_Z), with

Here, , , λ_C(0), λ_max(0) and λ(κ_A, 0) represent the parameters or variables free from Lacz perturbation, just as those in Appendix 2.3. Since the noise on the multiplier term (i.e., 1 − ϕ_Z/ϕ_max) is negligible, the CV of λ_C(ϕ_Z) (i.e., ) is unaffected by ϕ_Z. By combining Eqs. S46 and S48, we obtain the relations between the normalized energy fluxes and growth rate:

where λ(κ_A, ϕ_Z), and follow Eqs. S50 and S55 accordingly. For a given value of ϕ_Z, i.e., ϕ_Z is fixed, then, λ(κ_A, ϕ_Z) changes monotonously with κ_A. Combining Eqs. S55–S56 and S30, we have the relation between the normalized fluxes , and growth rate (here ϕ_Z is a parameter):

In Fig. 2C. we show that the model predictions (Eq. S57) quantitatively agree with the experiments(Basan et al., 2015).

Meanwhile, we can also perturb the growth rate by tuning ϕ_Z in a stable culturing environment with fixed concentration of a Group A carbon source (i.e., given [A]). In fact, for this case there is a distribution of κ_A values due to the extrinsic noise in , yet this distribution is fixed. For convenience of description, we still referred to it as fixed κ_A. Then, combining Eqs. S30, S50, S55 and S56, we get:

Here, λ(κ_A, 0) remains unaltered as κ_A is fixed. Therefore, in this case, and are proportional to λ, where the slopes are both functions of κ_A. More specifically, the slope of is a monotonically increasing function of κ_A, while that of is a monotonically decreasing function of κ_A. In Fig. 2B, we see that the model predictions (Eq. S58) agree quantitatively with the experiments(Basan et al., 2015).

In fact, the growth rate can be altered by tuning ϕ_Z and κ_A simultaneously. Then, the relations among the energy fluxes, growth rate and ϕ_Z still follow Eq. S57 (here ϕ_Z is a variable). In a 3-D representation, these relations correspond to a surface. In Fig. 2A, we show that the model predictions (Eq. S57) match well with the experimental data(Basan et al., 2015).

Appendix 3.2 Energy dissipation

In practice, energy dissipation breaks the proportional relationship between energy demand and biomass production. Thus, Eq. S25 changes to:

where W represents the dissipation coefficient. In fact, maintenance energy contributes to energy dissipation, and we define the maintenance energy coefficient as w₀. In bacteria, the impact of w₀ is often negligible, particularly for all the analysis in the sections above. While in tumor cell, w₀ plays a much more significant role.

The introduction of energy dissipation leads to a modification to Eq. S26: combining Eq. S59 and Eq. S16, we have:

Then, Eq. S29 changes to:

Consequently, if ε_r > ε_f, the best strategy for the cell is:

and if ε_f > ε_r, the best strategy is:

Then, the growth rate turns into a bivariate function of both κ_A and w:

Clearly, is still a constant, while and changes into functions of w:

For a cell population, in the homogeneous case, and follow:

To compare with experiments, we assume the same extent of extrinsic noise in (and thus κ_i) as that specified in Appendix 2.3. Combining Eqs. S45 and S65, then, λ_C(w) approximately follows a Gaussian distribution:

where and represent the mean and standard deviation of λ_C(w), and

Here, , , λ_C(0), λ_max(0), and λ(κ_A, 0) represent the parameters or variables free from energy dissipation. In fact, there is a distribution of values for . For approximation, we use the deterministic value of in Eq. S68, and then the CV of λ_C(w) is roughly unperturbed by w. Combining Eqs. S46, S66 and S67, we have:

Since the dissipation coefficient w is tunable in experiments, for a given value of w, λ(κ_A, w) changes monotonously with κ_A. Combining Eqs. S68–S69 and S30, we have (here w is a parameter):

The comparison between model predictions (Eq. S70) and experimental results(Basan et al., 2015) is shown in Fig.3B, which agrees quantitatively. Meanwhile, the growth rate can also be perturbed by changing κ_A and w simultaneously. Then, the relations among the energy fluxes, growth rate and follow Eq. S70 (here w is a variable). In a 3D representation, these relations correspond to a surface. As shown in Fig. 3A, the model predictions (Eq. S70) agree quantitatively with the experimental results(Basan et al., 2015).

Appendix 3.3 Translation inhibition

In E. coli, the translation rate can be modified by adding different concentrations of translation inhibitors, e.g., chloramphenicol (Cm). The net effect of this perturbation is represented as:

where ι stands for the inhibition coefficient with ι > 0, and (1+ ι)^-1 represents the translation efficiency. Thus, Eq. S32 changes to:

First, we consider the case of neglecting the maintenance energy, i.e., w₀ = 0. Then, the growth rate takes the following form:

where λ(κ_A, 0) and ψ(κ_A, 0) represent the terms free from translation inhibition. Thus, and change into functions of ι:

In the homogeneous case, and follow:

To compare with experiments, we assume that there exists extrinsic noise in and κ_i as specified in Appendix 2.3. Combining Eqs. S45 and S74, λ_C(ι) can be approximated by a Gaussian distribution:

where and represent the mean and standard deviation of λ_C(ι), with

Here, , , , λ_C(0) and λ_max(0) stand for the terms free from translation inhibition. Basically, there are distributions of values for , and . For approximation, we use the deterministic values of these terms in Eq. S77, and then the CV of λ_C(ι) can be approximated by λ_C(0). Combining Eqs. S46, S75 and S76, we have:

In the experiments, the inhibition coefficient ι is controllable by tuning the concentration of translation inhibitor. For a given value of ι, λ(κ_A, ι) changes monotonously with κ_A. Combining Eqs. S30 and S78, we have (ι is a parameter here):

where and follow Eq. S77. The growth rate can also be perturbed by altering κ_A and ι simultaneously. Then, the relations among the energy fluxes, growth rate and ι still follow Eq. S79 (here ι is a variable).The comparison between Eq. S79 and experimental data(Basan et al., 2015) is shown in Appendix-fig. 2D (3-D) and 2E (2-D). Overall, it shows good consistency, however, there is still a noticeable discrepancy when ι is large (i.e., with high concentration of translation inhibitor). Then, it led us to consider the maintenance energy w₀, which is small but may account for this discrepancy. Then, λ(κ_A, ι) changes into:

while and still follow Eq. S74, although the forms of λ_C(0) and λ_max(0) change into:

In the homogeneous case, and follow:

To compare with experiments, we assume that the extrinsic noise follows that specified in Appendix 2.3. Combining Eqs. S45, S74 and S81, then, λ_C(ι) approximately follows a Gaussian distribution:

Here and still follow Eq. S77, while and change accordingly with Eq. S81. For approximation, we use the deterministic values of the relevant terms in Eq. S77, and then the CV of λ_C(ι) is roughly the same as λ_C(0). Combining Eqs. S46, S82 and S83, we have:

Thus, for a given ι, λ(κ_A, ι) changes monotonously with κ_A. Combining Eqs. S30 and S84, we have (here ι is a parameter):

The growth rate and fluxes can also be perturbed by altering κ_A and ι simultaneously. The relations among the energy fluxes, growth rate and ι would still follow Eq. S85 other than that now ι is regarded as a variable. Assuming that there is a tiny amount of maintenance energy. Basically, we assign w₀ = 2.5(h^-1). Then, we see that the experimental results(Basan et al., 2015) agree quantitatively well with the model predictions (Fig.3C-D).

Appendix 4 Overflow metabolism in substrates other than Group A carbon sources

Due to the topology of metabolic network, Group A carbon sources follow the equation (Eq. S47) of overflow metabolism upon κ_A perturbation (i.e., varying the type or concentration of a Group A carbon source). This has been demonstrated clearly in the above analysis, which agrees quantitatively with experiments. However, further analysis is required for substrates other than Group A sources due to the topological differences in carbon utilization(Wang et al., 2019). Basically, substrates entering from glycolysis or the points before acetyl-CoA are potentially involved in overflow metabolism, while those join from the TCA cycle are not relevant to this behavior. Still, mixed carbon sources are likely to induce a different profile of overflow metabolism, so long as there is a carbon source coming from glycolysis.

Appendix 4.1 Pyruvate

The coarse-grained model for pyruvate utilization is shown in Fig. 3E. Here, nodes M₁, M₂, M₃, M₄, M₅ follow everything depicted in Appendix 2.1. Each biochemical reaction follows Eq. S5 with b_i = 1 except that 2M₂→M₁ and M₃+M₅→M₄. By applying flux balance to the stoichiometric fluxes, combining with Eq. S8, we have:

For energy production, we convert all the energy currencies into ATPs, and then,

where β₇ = 1, β₈ = 2, β₃ = 2, β₄ = 6, β₆ = 1, β₉ = 6, β_a1 = 4 (Neidhardt et al., 1990; Sauer et al., 2004), and J_E follows Eq. S25. By applying the substitutions specified in Eqs. S9, S12, S14–S18, combined with Eqs. S4, S10, S22, S23, S25, S86–S87, and the constraint of proteome resource allocation, we have:

where . κ_i is approximately a constant which follows Eq. S20 for each of the intermediate node. The substrate quality of κ_py varies with the external concentration of pyruvate ([py]),

From Eq. S88, all ϕ_i can be expressed by , , and λ:

By substituting Eq. S90 into Eq. S88, we have:

Here, and stand for the normalized energy fluxes of respiration and fermentation, respectively, with

The coefficients and represent the proteome energy efficiencies of respiration and fermentation, respectively (Appendix-fig. 1C-D), with

is the proteome efficiency of biomass pathway (Appendix-fig. 1E), with

φ_py is the energy demand coefficient (a constant), with

Evidently, Eq. S91 is identical in form with Eq. S29. The growth rate changes into κ_py dependent:

When κ_py is very small, combined with Eq. S93, then,

Obviously,

As long as

where the superscript “(ST)” stands for the saturated concentration, then,

and there exists a critical value of κ_py, denoted as , with

Here, is the growth rate at the transition point, and stands for either or . In Appendix-fig. 2H, we show the dependencies of , and λ(κ_py) on κ_py in a 3-dimensional form. In the homogeneous case, and follow:

Defining , and then, is the relevant range of the x axis. To compare with experiments, we assume the same extent of extrinsic noise in as specified in Appendix 2.3. Then, approximately follows a Gaussian distribution:

where and stand for the mean and standard deviation of . Then, the relations between the normalized energy fluxes and growth rate are:

Combined with Eq. S92, we have:

In Fig. 3F, we show that the model predictions (Eq. S105) agree quantitatively with the experimental results(Holms, 1996).

Appendix 4.2 Mixture of a Group A carbon source with extracellular amino acids

In the case of a Group A carbon source mixed with amino acids, the coarse-grained model is shown in Appendix-fig. 2A.In fact, this model can be used to analyze mixtures with one or multiple types of extracellular amino acids. Here, Eqs. S21, S22, S24 and S25 still apply, but Eq. S23 changes to (the case of i remains the same as Eq. S23):

Here, represents the number of carbon atoms in a molecule of Pool i. For simplicity, we assume:

In the case where all 21 types of amino acids are present and each in saturated concentration (denoted as “21AA”), we have:

where ϕ_i and κ_i are defined following Eqs. S9 and S12. Since the cell growth rate elevates significantly with the mixture of amino acids, we deduce that Pools a2-d are supplied by amino acids in growth optimization, with

Basically, amino acids should be more efficient in the supply of biomass production than the Group A carbon source for Pools a2-d, i.e.,

In practice, the requirement is even higher for proteome efficiency using amino acids, since the biomass production pathway is accompanied by energy production in the case of Group A carbon sources, yet not for amino acids. Combining Eqs. S108 and S109, we have:

where , follow Eq. S30, while ε_r and ε_f follow Eq. S31. is the proteome efficiency of biomass pathway under this nutrient condition, with

φ_21AA is the energy demand coefficient, with

Combining Eqs. S111 and S31, it is easy to obtain the formula for the growth rate:

In fact, Eqs. S37–S42 still apply. satisfies Eq. S43, while and are:

When extrinsic noise is taken into account, approximately follows a Gaussian distribution:

and the normalized fluxes , change into:

In fact, the above analysis can be extended to cases where a Group A carbon source is mixed with arbitrary combinations of amino acids. Eqs. S111, S114–S117 would remain in a similar form, while Eqs. S112–S113 would change depending on the amino acid combinations.In Appendix-fig.2B-C, we show the comparisons between model predictions (see also Appendix7.2 and Eq. S157) and experimental data(Basan et al., 2015; Wallden et al., 2016) in mixtures of 21 or 7 types of amino acids together with a Group A carbon source, which agree quantitatively.

Appendix 5 Enzyme allocation upon perturbations

Appendix 5.1 Carbon limitationwithin Group A carbon sources

In Eq. S28, we present the model predictions of the dependencies of enzyme protein fractions on growth rate and energy fluxes. To compare with experiments, we assume the same extent of extrinsic noise in as specified in Appendix 2.3. In fact, relative protein expression data for enzymes within glycolysis and the TCA cycle are available from existing studies, which are comparable to the ϕ₁-ϕ₄ enzymes of our model (Fig. 1B). Upon κ_A perturbation, κ_A is a variable while w₀ = 0 is fixed (see Appendix 1.5). Combining Eqs. S28 and S47 (with w₀ = 0), we obtain:

In Appendix-fig.3C-D, we show the comparisons between model predictions (Eq. S118, w₀ = 0) and experimental data(Hui et al., 2015), which are consistent overall. We proceed to consider the influence of maintenance energy as specified in Appendix3.2. Here, we still choose w₀ = 2.5 (h^-1) as previously adopted in Appendix 3.3. Then, Eq. S28 still holds, combined with Eq. S85 in the condition that ι = 0, we have:

In Fig. 4A-B, we show that the model predictions (Eq. S119, w₀ = 2.5 (h^-1)) generally agree with the experiments(Hui et al., 2015). However, there are different basal expressions of these enzymes, which are probably due to living demands other than cell proliferation, such as preparation for starvation(Mori et al., 2017) or alteration in the type of the nutrient(Basan et al., 2020; Kussell and Leibler, 2005).

Appendix 5.2 Overexpression of useless proteins

In the case of ϕ_Z perturbation under each nutrient condition with fixed κ_A (see Appendix 3.1), we consider the same extent of extrinsic noise in as specified in Appendix 2.3. The relation between enzyme allocation and growth rate can be obtained by combining Eqs. S28 and S58 (with w₀ = 0):

Here λ(κ_A, 0) is the growth rate for ϕ_Z = 0 and thus is a parameter rather than a variable. The growth rate is defined as λ(κ_A, ϕ_Z), which follows Eq. S50. Thus, ϕ_i is proportional to the growth rate λ.In Appendix-fig. 3E-F, we see that the model predictions (Eq. S120) agree with the experiments(Basan et al., 2015) overall. Next, we consider the influence of maintenance energy with w₀ = 2.5 (h^-1). Combining Eqs. S28, S58 and S85 (with ι = 0), we get:

Here, the growth rate is defined as λ(κ_A, ϕ_Z), and λ(κ_A, 0) is a parameter rather than a variable. Thus, ϕ_i is a linear function of the growth rate λ, with a positive slope and a positive y-intercept.In Fig. 4C-D and Appendix-fig.3I-J, we show that the model predictions (Eq. S121) agree quantitively with the experimental data(Basan et al., 2015).

Appendix 5.3 Energy dissipation

In the case of energy dissipation under each nutrient condition, w is perturbed while κ_A is fixed. The relation between protein allocation and growth rate can be obtained by combining Eqs. S28 and S70. However, since is explicitly present in Eq. S70, we need to reduce this variable to obtain the growth rate dependence of enzyme allocation. In fact, from Eq. S64, we have:

Here, (satisfying Eq. S64), which is a parameter rather than a variable. “θ” stands for the Heaviside step function. Thus, we have:

where the energy dissipation coefficient w is regarded as a function of the growth rate.

Combining Eqs. S28, S70 and S123, we get:

where w(λ) follows Eq. S123. When κ_A lies in the vicinity of or w is small so that

then we have:

and thus,

Note that in Eq. S123, w is a linear function of λ with a negative slope. Thus ϕ_i exhibits a linear relation with λ when Eq. S125 is satisfied (see Eq. S127). In fact, the slope of ϕ₄ is surely negative (combining Eqs. S64, S123 and S127), while the slope sign of other ϕ_i depends on parameters. For a given nutrient, the enzymes corresponding to the same ϕ_i should exhibit the same slope sign. Another restriction is that if the slope sign of ϕ₁ is negative, then the slope sign of ϕ₂ is surely negative. In Appendix-fig. 3K-N, we show that our model results agree well with the experimental data(Basan et al., 2015) (Eq. S127).

Appendix 6 Other aspects of the model

Appendix 6.1 A coarse-grained model with more details

To compare with experiments, we consider a coarse-grained model with more details as shown in Appendix-fig. 2F.Here, nodes M₆, M₇ represent GA3P and DHAP, respectively. Other nodes follow the descriptions specified in Appendix 2.1. Each biochemical reaction follows Eq. S5 with b_i = 1 except that M₁→M₆+M₇ and M₃+M₅→M₄. By applying flux balance to the stoichiometric fluxes, combined with Eq. S8, we obtain:

While Eqs. S22–S25 still hold. By applying the substitutions specified in Eqs. S9, S12, S14–S18, combined with Eqs. S4, S10, S22–S25, S128, and the constraint of proteome resource allocation, we get:

Then, Eq. S28 still hold, while ϕ₁₀ and ϕ₁₁ are:

By substituting Eqs. S28 and S130 into Eq. S129, we get:

where “dt” stands for details. Eqs. S30 and S33 still hold. and represent the proteome energy efficiencies of respiration and fermentation, respectively, with

is the proteome efficiency of biomass pathway, with

Appendix 6.2 Estimation of the in vivo enzyme catalytic rates

We use the method introduced by Davidi et al.(Davidi et al., 2016),combined with proteome experimental data(Basan et al., 2015) (Appendix-table2), to estimate the in vivo enzyme catalytic rates. Combining Eqs. S28 and S130, we have:

Here, , , λ and ϕ_i (i = 1–6, 10–11) are measurable from experiments (see Appendix 8.1 and Appendix-table2). Thus, we can obtain the in vivo values of κ_i from Eq. S134. Combined with Eqs. S17 and S20, we have

Eq. S135 is the in vivo result for the enzyme catalytic rate. In Appendix-fig. 2G, we show a comparison between in vivo and in intro results for k_cat values of enzymes within glycolysis and the TCA cycle, which are roughly consistent. In the applications, we prioritized the use of in vivo results for enzyme catalytic rates, and use in intro data as a substitute when there were vacancies.

Appendix 6.3 Comparison with existing models thatillustrate experimental results

For the coarse-grained model described in Appendix 2, the normalized stoichiometric influx of a Group A carbon source is given by:

Combined with the first equation in Eq. S28, we obtain:

where e_r = β₁ + 2(β₂ + β₃ + β₄), e_f = β₁ + 2(β₂ + β₆), and . Evidently, e_r, e_f and ϑ are constant parameters. In this subsection, we highlight the major differences between our model presented in Appendix 2 and existing models that may illustrate the growth rate dependence of fermentation flux in the standard picture(Basan et al., 2015; Holms, 1996; Meyer et al., 1984; Nanchen et al., 2006) of overflow metabolism.

Based on the modeling principles rather than the detailed mechanisms, there are two major classes of existing models that can illustrate experimental results. In fact, both classes of models regard the proteome energy efficiencies ε_r and ε_f as constants, with ε_f > ε_r if used, or follow functionally equivalent propositions. In our model, however, ε_r and ε_f are both functions of κ_A, which vary significantly upon nutrient perturbation, where and (see Eqs. S38, S40–S41). Furthermore, there are significant differences in the modeling/optimization principles, as listed below.

The first class of models(Chen and Nielsen, 2019; Majewski and Domach, 1990; Niebel et al., 2019; Shlomi et al., 2011; Varma and Palsson, 1994; Vazquez et al., 2010; Vazquez and Oltvai, 2016)optimize the ratio of biomass outflow to carbon influx: . Either to optimize the growth rate for a given carbon influx, or to minimize the carbon influx for a given growth rate. Since respiration is far more efficient than fermentation in terms of the energy production per unit carbon, to optimize the ratio , cells would preferentially use respiration when the carbon influx is small. As carbon influx increases above a certain threshold, factors such as proteome allocation direct cells towards fermentation in an analog way since they consider ε_f > ε_r. Our model is significantly different from this class of models in the optimization principle, where we purely optimize the cell growth rate for a given nutrient condition, without any further constraints on the carbon influx.

The second class of models, represented by Basan et al.(Basan et al., 2015), also adopt the optimization of in the interpretation of their model results. However, the growth rate dependence of fermentation flux was derived prior to the application of growth rate optimization. In fact, Eqs. S29 and S137 in our model are very similar form to Basan et al.(Basan et al., 2015), yet there are critical differences as we list below. In Eq. S29, by regarding and as the two variables in a system of linear equations, we obtain the following expressions:

In Basan et al.(Basan et al., 2015), Eq. S138 is considered to be the relation between and λ upon nutrient (and thus ) perturbation, while ε_r and ε_f are regarded as constants through the perturbation. By contrast, in our model, Eq. S138 is the constraint under a given nutrient condition with fixed κ_A, yet not relevant to nutrient perturbation. For wild-type strains, if ε_r (κ_A) > ε_f (κ_A) (or vice versa), then the optimal solution is and , with . This solution (which satisfies Eq. S138) corresponds to a point rather than a line in the relation between and λ upon κ_A perturbation.

Appendix 6.4 Explanation of the Warburg effect in tumor cells

Our model and analysis shown in Appendix 2 can be naturally extended to explain the Warburg effect in tumor cell metabolism with the following modifications in the model settings. In the applications for tumor cell metabolism, the fermentation flux changes from acetate secretion rate into the lactate secretion rate, and thus the stoichiometric coefficients (β_i) for ATP production change accordingly. Consequently, in the coarse-grained model shown in Fig. 1B, M₃ stands for pyruvate, and β₂, β₆ change into β₂ = 1, β₆ = -2 (Nelson et al., 2008).

Evidently, Eqs. S37–S38 still hold. As long as (Eq. S40), then there exists a critical switching point for κ_A (denoted as , see Eq. S41), below which the respiration pathway is more efficient, while above , the fermentation pathway is more efficient in the proteome energy production. To maximize proliferation rate, the tumor cells would preferentially use the respiration pathway upon starvation while use the aerobic glycolysis pathway with abundant nutrient (Warburg effect). Based on the facts that the rate of glucose metabolism is 10–100 times faster in aerobic glycolysis than that of respiration(Liberti and Locasale, 2016), and oncogenes induce Warburg effect intrinsically to promote cell proliferation(Hanahan and Weinberg, 2011; Liberti and Locasale, 2016), it is suggested that . Thus, the aerobic glycolysis with abundant nutrients provides a fitness advantage for tumor cell proliferation.

Appendix 7 Probability density functions of variables and parameters

Appendix 7.1 Probability density function of κ_i

Enzyme catalysis is crucial for the survival of living organisms, as it can significantly accelerate the rate of a biochemical reaction by moderating the energy barrier between the substrate and product(Nelson et al., 2008).However, the maximal turnover rate of enzymes, k_cat values, vary notably between the in vivo and in vitro measurements(Davidi et al., 2016). Recent studies suggest thatdifferences in the aquatic medium should be the causes(Davidi et al., 2016; García‐Contreras et al., 2012). In particular, the concentrations of potassium and phosphate have a biginfluence on k_cat (García‐Contreras et al., 2012), which possess a certain degree of variation among the cellpopulations under intracellular conditions(García‐Contreras et al., 2012).For simplicity, we assume that the turnover rate of each E_i enzyme follows a Gaussian distribution with among cells (extrinsic noise(Elowitz et al., 2002), denoted as η_ext),and then the probability density function of is given by:

When the CV of the distribution (i.e., ) is less than 1/3, is almost identical to . In this regard, follows the positive inverse of Gaussian distribution (IOG), and the probability density function is:

where and .

Meanwhile, due to the stochastic nature of biochemical reactions, we apply Gillespie’s chemical Langevin equation(Gillespie, 2000) to account for the intrinsic noise(Elowitz et al., 2002)(denoted as η_int). For cell size regulation of E. coli within a cell cycle, the cell mass at initiation of DNA per chromosome origin remains constant(Donachie, 1968). Thus, the duration of enzyme E_i to finish a catalytic job (with a timescale of ) can be approximated as the first passage time of a stochastic process, with

Here , Θ is proportional to the cell volume, and Γ_i(t) is the independent temporally uncorrelated Gaussian white noise. Then, for a given value of , the first passage time T_Γ follows an Inverse Gaussian distribution (IG)(Folks and Chhikara, 1978):

where , , and the variance is . Thus, we can obtain the CV:

which is inversely scaled with the square root of cell volume. Evidently, the intrinsic and extrinsic noise make orthogonal contributions to the total noise(Elowitz et al., 2002)(denoted as θ_tot):

In fact, when the CV is small (i.e., CV<<1), both IOG and IG distributions converge into Gaussian distributions (Appendix-fig.4). In the back-of-the-envelope calculations, we approximate x in all the denominator terms of IOG(x, μ, ζ) and IG(x, μ, ζ) as μ (since CV<<1). Then, both IOG and IG distributions can be approximated as follows:

with , and

with . Rigorously, we show below that IG(x, μ, ζ) shrinks to be N(μ, μ³/ζ) when the CV is small. Here, the characteristic function of the variable x is(Folks and Chhikara, 1978; Van Kampen, 1992):

and thus,

When the variance σ² ≡ μ³/ζ is very small, basically, we require 2μ²k/ζ = 2σ²k/μ ≪ 1, and then . Thus,

Then, we have:

In fact, although intrinsic noise affects the short-term measurement of enzyme catalytic rate and growth rate at the single-cell level, its contribution in the long-term is negligible. Thus, we approximate η_tot ≈ η_ext. Combined with Eqs. S145–S146, it is easy to check that shares roughly identical CV as :

For convenience, in the model analysis, we approximate both IOG and IG distributions as Gaussian distributions. Then, all are normally distributed independent random variables following Gaussian distributions:

Using the properties of Gaussian distributions, for a series of constant real numbers γ_i, the summation of , where we defined it as , follows a Gaussian distribution(Van Kampen, 1992):

with and . The relation between κ_i and is shown in Eq. S12. To optimize growth rate, each κ_i of the intermediate nodes satisfies Eq. S20, while κ_A satisfies Eq. S27. Thus, for a given nutrient condition ([A] is fixed), all the ratios are constants. Combined with Eqs. S139, S145–S146, S152, then, the distributions of all κ_i and 1/κ_i can be approximated as Gaussian distributions:

where and are the means of κ_i and 1/κ_i, while and are their standard deviations. Using the properties of Gaussian distributions combined with Eq. S31, S32, S36, S42–S43, S145–S146 and S153, then ε_r, ε_f, ψ, λ_r, λ_f, and λ_C also roughly follow Gaussian distributions.

Appendix 7.2 Probability density function of the growth rate λ

From Appendix 7.1, we note that λ_r and λ_f (see Eq. S36) roughly follow Gaussian distributions, with

where and represent the mean and standard deviation, respectively. We further assume that the correlation between λ_r and λ_f is ρ_rf. From Eq. S36, we see that the growth rate λ takes the maximum of λ_r and λ_f, i.e.,

Then, the cumulative distribution function of λ is , where . Thus, the probability density function of the growth rate λ is given by:

In Appendix-fig. 2B, we show that Eq. S157 quantitatively illustrates the experimentaldata of E.coli under the relevant conditions.

Appendix 8 Model comparison withexperiments

Appendix 8.1 Flux comparison with experiments

In Appendix 6.2, we see that the values of and are required to calculate the in vivo enzyme catalytic rates of the intermediate nodes. Here, we use J_acetate and to represent the stoichiometric fluxes of acetate from the fermentation pathway and CO₂ from the respiration pathway, respectively. Combined with the stoichiometric coefficients of both pathways, we have:

By further combining with Eqs. S16–S17, we get:

In fact, the values of J_acetate and scale with the mass of the “big cell”, which increases over time. In experiments, the measurable fluxes are typically expressed in the unit of mM/OD₆₀₀/h(Basan et al., 2015). Thus, we define and as the fluxes of J_acetate and in the unit of mM/OD₆₀₀/h, respectively. The superscript “(M)” represents the measurable flux in this unit. For E. coli, we use the following biochemical data collected from published literature: 1 OD₆₀₀ roughly corresponds to 6×10⁸ cells/mL(Stevenson et al., 2016), the average mass of a cell is 1pg(Milo and Phillips, 2015), the biomass percentage of the cell weight is 30%(Neidhardt et al., 1990), the molar mass of carbon is 12g(Nelson et al., 2008), r_carbon = 0.48 (Neidhardt et al., 1990) and r_protein = 0.55 (Neidhardt et al., 1990).Combined with the values of r_i (see Appendix 1.2) and , where , , , , and (Nelson et al., 2008), we have:

From Eq. S18, we obtain the values of η_i for each precursor pool : η_a1 = 0.15, η_a2 = 0.30, η_b = 0.35, η_c = 0.09, η_d = 0.11. Still, the value of η_E is required to compare the growth rate dependence of fermentation/respiration fluxes between model results and experiments, which we will specify in Appendix 8.2.

Appendix 8.2 Model parameter settings using experimental data

We have collected biochemical data of E. coli shown in Appendix-tables 1–2 to set the model parameters. This includes the molecular weight (MW) and in vitro k_cat values of the catalytic enzymes, as well as the proteome and flux data used to calculate the in vivo turnover numbers. To reduce measurement noise, we take the average rather than the maximum value of in vivo k_cat from calculations with data from four cultures (see Appendix-table 2). Here, we prioritize the use of in vivo k_cat wherever applicable unless there is a vacancy (see Appendix-table1).

Note that our models are coarse grained. For example, the flux J₃ shown in Fig. 1B actually corresponds to three different reactions in the metabolic network (see Fig.1A and Appendix-table1),which we label as (i=1, 2, 3). For each , there are corresponding variables/parameters of , , , satisfying Eqs. S8, S9 and S12, Evidently, (i=1, 2, 3), and it is straightforward to derive the following relation between and κ₃:

In fact, Eq. S161 can be generalized to determine the values of other κ_i in the coarse-grained models combined with the biochemical data. For the coarse-grained model of Group A carbon source utilization shown in Fig. 1B, we have the values for parameters κ_i (i=1, …, 6), and then . Evidently, , , and thus . For pyruvate, we have (see Eqs. S43 and S101), and it is easy to check that .

For the remaining model parameters, note that we have classified the inactive ribosomal-affiliated proteins into the Q-class, and then ϕ_max = 48 % (Scott et al., 2010). The values of κ_t is obtainable from experiments: the translational speed is 20.1aa/s(Scott et al., 2010), with 7336 amino acids per ribosome(Neidhardt, 1996) and ϛ ≈ 1.67 (Neidhardt, 1996; Scott et al., 2010)(see Appendix 1.1), hence κ_t = 1/610 (s^-1).However, there are insufficient data to determine the values of κ_i (i=a1, a2, b, c, d) from the entry point metabolites to the precursor pools. Basically, it involves many steps, and thus these values should be considerably large. Here, we lump sum the contributions of κ_t and κ_i (i=a1, a2, b, c, d) by defining a composite parameter:

We proceed to estimate the values Ω of φ and using experimental data(Basan et al., 2015)for wild-type strainson the relation(Fig. 1C),and then all the rest of model parameters are set accordingly.

For the case of w₀ = 0, where all k_cat values follow a Gaussian distribution with an extrinsic noise of 25% CV, which is the general setting we use unless otherwise specified, then φ = 10.8 and Ω = 1345 (s). Accordingly, we have η_E = 14.78, , and , where the CV of the extrinsic noise for Ω is estimated using the averaged CV of other κ_i. For the translation inhibition effect of Cm, we estimate the values for ι as , , , where the superscript stands for the concentration of Cm, and the subscript represents the choice of w₀.

For pyruvate, with the value of η_E, we get φ_py = 14.82. However, there is still a lack of proteome data to determine the value of κ₉, which actually involves many steps in the metabolic network and thus can be considerably large. Here we define another composite parameter and estimate its value as from the growth rate data for E. coli measured under the relevant nutrient conditions(Basan et al., 2015),where the subscript “Gg” stands for glucogenesis. Then, , and , where the same CV of extrinsic noise for Ω applies to Ω’_Gg.

For the case of a Group A carbon source mixed with 21 amino acids (21AA, with saturated concentrations), then φ_21AA = 14.2. Comparing Eq. S32 with Eq. S112, the parameter Ω should change to .Obviously, 1/κ_t < Ω_21AA < Ω, and we estimate Ω_21AA = 1000 (s) from the growth rate data for E. coli measured under the relevant nutrient conditions(Wallden et al., 2016). Then, we have , and .

For the case of a Group A carbon source mixed with 7 amino acids (7AA: His, Iso, Leu, Lys, Met, Phe, and Val), similar to the roles of φ_21AA and Ω_21AA, we define φ_7AA and Ω_7AA. Using the mass fraction of the 7AA combined with Eq. S18, we have φ_7AA = 11.6. For the value of Ω_7AA, evidently, Ω_21AA < Ω_7AA < Ω, and we estimate φ_7AA = 1215 (s) from the growth rate data for E. coli measured under the relevant culture media(Basan et al., 2015). Then, , and .

For the case of w₀ = 2.5 (h^-1), we have φ = 8.3, and thus η_E = 12.28,, while other parameters such as Ω, and remain the same values as for w₀ = 0. Nevertheless, the values for ι upon translation inhibition with Cm are influenced by the choice of w₀ = 0, where the values of ι change into , , .

From Appendix 7.1–7.2, combined with Eq. S114, the distributions of and can be approximated by Gaussian distributions:

where and stand for the mean values, while and represent the standard deviations. For the case of glucose mixed with 21AA (labeled as “Glucose+21AA”), the distribution of the growth rate follows Eq. S157. With Ω_21AA = 1000 (s), we have , (both definitions follow Eq. S163), and ρ_rf ≈ 1.0 (obtained with numerical results).

For the case of succinate mixed with 21AA (labeled as “Succinate+21AA”), the respiration pathway is always more efficient since succinate lies within the TCA cycle, then the cell growth rate (defined as ) would just take the value of the respiration one and thus follows a Gaussian distribution:

For the case that acetate is the sole carbon source, evidently, the cells only use the respiration pathway, and thus the growth rate (defined as λ_acetate) follows a Gaussian distribution:

With the measured growth rate data(Wallden et al., 2016), we estimate and . To illustrate the distribution of growth rates , and λ_acetate shown in Appendix-fig. 2B, if there was no other source of noise, extrinsic noise with a CV of 40% would be required for each k_cat value. Then, , , , and . Allowing for that intrinsic noise may also play a non-negligible role for the observed single-cell growth rate (which is not a long-term average), we still use extrinsic noise with a CV of 25% other than those shown in Appendix-fig. 2B.

Appendix 8.3 Notes on the application of reference data

Data calibration

Throughout our manuscript, we use experimental data from the original references except for two calibrations. The first calibration is shown in the footnote to Appendix-table2. With this calibration, the data(Basan et al., 2015)in Appendix-table2 align with the curve shown in Fig. 1C, which includes experimental data from other sources. The second calibration is for the data shown in Figs. 3F and 1C (chemostat data). The unit in the original reference(Holms, 1996) is mmol/(dry mass)g/h.To convert this into the unit mM/OD₆₀₀/h applied in our text, the conversion factor should be 0.18. Here, we deduce that only 60% of the measured dry biomass in centrifuged material is effective to calibrate with other experimental results. Therefore, there is a calibration factor of 0.6, and the conversion factor changes into 0.3.

Data of the inducible strains

We note that part of the experiment data in the original references(Basan et al., 2015; Hui et al., 2015)were obtained using strains with titratable systems (e.g. titratable ptsG, LacY). Basically, the relation of these inducible strains align with the same curve as that of the wild-type E.coli (Fig. 1C). Given that evolution treatment is not involved for the inducible strain, we regard the titration perturbation as a technique to fool the strains. That is, these inducible strains behave as if they were cultured in a less efficient Group A carbon source.

Experimental data sources

The batch culture data shown in Fig. 1C(labeled with minimum/rich media or inducible strains) and Appendix-fig.2C were taken from the source data of the reference’s figure 1(Basan et al., 2015). The chemostat data shown in Fig. 1C were taken from the reference’s table 7(Holms, 1996). The data shown in Fig. 1D were taken from the reference’s extended data figure 3a(Basan et al., 2015) with the calibration specified in the footnote to Appendix-table2.

The data shown in Fig. 2A were adopted from the reference’s extended data figure 4a-b(Basan et al., 2015). The data shown in Fig. 2B were taken from the source data of the reference’s figure 2a(Basan et al., 2015). The data shown in Fig. 2C were taken from the source data of the reference’s figure 3a(Basan et al., 2015). The data shown in Fig. 3A-B were taken from the source data of the reference’s figure 3d(Basan et al., 2015).

The data shown in Fig. 3C-D and Appendix-fig.2D-E were taken from the source data of the reference’s figure 3c(Basan et al., 2015). The data shown in Fig. 3F were taken from the reference’s table 7(Holms, 1996), with a calibration factor specified in the above paragraph (“Data calibration”).

The data shown in Fig.4A-B and Appendix-fig.3A-D were taken from the reference’s table S2 with the label “C-lim”(Hui et al., 2015).We excludedthe reference’s data withλ=0.45205 (h^-1)as there are other unconsidered factors involved during slow growth(Dai et al., 2016) (for λ<0.5 h^-1), and we suspect that there may be unknown calibration factors. The data shown in Fig. 4C-D and Appendix-fig.3E-N were adopted from the reference’s extended data figure 6–7(Basan et al., 2015).

The gene names depicted in Appendix-fig. 1B were identified using KEGG database. The data shown in Appendix-fig. 2G were drawn from Appendix-table1, which includes the original references themselves. The flux data presented in Appendix-table2 were obtained from the reference’s extended data figure 3a(Basan et al., 2015), with the calibration specified in the footnote. The proteome data shown in Appendix-table2 were taken from the reference’s supplementary Table N5(Basan et al., 2015).

Appendix-table 1.

Appendix-table 2.

Appendix-figure 1.

Appendix-figure 2.

Appendix-figure 3.

Appendix-figure 4.