Overflow metabolism originates from growth optimization and cell heterogeneity

  1. Xin Wang  Is a corresponding author
  1. School of Physics, Sun Yat-sen University, China

eLife Assessment

This valuable study tackles the well-established overflow metabolism issue by applying a coarse-grained metabolic flux model to predict how individual cells execute various energy strategies, such as respiration versus fermentation. The model's population average is convincing enough to align with experimental observations on overflow metabolism. The potential source of metabolic or proteomic heterogeneity of individual cells remains an open question to be studied. How individual cells adjust their metabolic strategies also requires future study of the underlying mechanisms. Overall, this work provides a key aspect on cell-to-cell variability on general metabolic response.

https://doi.org/10.7554/eLife.94586.4.sa0

Abstract

A classic problem in metabolism is that fast-proliferating cells use seemingly wasteful fermentation for energy biogenesis in the presence of sufficient oxygen. This counterintuitive phenomenon, known as overflow metabolism or the Warburg effect, is universal across various organisms. Despite extensive research, its origin and function remain unclear. Here, we show that overflow metabolism can be understood through growth optimization combined with cell heterogeneity. A model of optimal protein allocation, coupled with heterogeneity in enzyme catalytic rates among cells, quantitatively explains why and how cells choose between respiration and fermentation under different nutrient conditions. Our model quantitatively illustrates the growth rate dependence of fermentation flux and enzyme allocation under various perturbations and is fully validated by experimental results in Escherichia coli. Our work provides a quantitative explanation for the Crabtree effect in yeast and the Warburg effect in cancer cells and can be broadly used to address heterogeneity-related challenges in metabolism.

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, 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 among 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, and particularly through 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; Peebo et al., 2015; 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; Zhuang et al., 2011). Notably, Basan et al., 2015 provided a systematic characterization of this process, including various types of experimental perturbations. Currently, prevalent explanations (Basan et al., 2015; Chen and Nielsen, 2019) hold that overflow metabolism arises from the proteome efficiency in fermentation being consistently higher than that in respiration. However, recent studies have shown that the measured proteome efficiency in respiration is actually higher than in fermentation for many yeast and cancer cells (Shen et al., 2024), even though these cells generate fermentation products through aerobic glycolysis. This finding (Shen et al., 2024) apparently contradicts the prevalent explanations (Basan et al., 2015; Chen and Nielsen, 2019). Furthermore, most explanations (Basan et al., 2015; Chen and Nielsen, 2019; Majewski and Domach, 1990; Shlomi et al., 2011; Varma and Palsson, 1994; Vazquez et al., 2010; Vazquez and Oltvai, 2016; Zhuang et al., 2011) rely on the assumption that cells optimize their growth rate for a given rate of carbon influx (i.e. nutrient uptake rate) under each nutrient condition (or its equivalents). However, this assumption remains open to further scrutiny, as the given factors in a nutrient condition are the identities and concentrations of the carbon sources (Molenaar et al., 2009; Scott et al., 2010; Wang et al., 2019), rather than the carbon influx. Therefore, 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 strategy of aerobic glycolysis? 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). For cancer cells, disrupting the growth control system and evading immune destruction from the host are prominent hallmarks of their survival (Hanahan and Weinberg, 2011), which in certain ways mimic the evolutionary pressure on microbes to optimize cell growth rate. In this study, we apply the optimal growth principle of microbes, which also roughly holds for cancer cells, to a heterogeneous framework to address the puzzle of aerobic glycolysis. We use Escherichia coli as a typical example to show that overflow metabolism can be understood from optimal protein allocation combined with 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; van Hoek et al., 1998) of overflow metabolism. Our model quantitatively illustrates the growth rate dependence of fermentation/respiration flux and enzyme allocation under various types of perturbations in E. coli. Furthermore, it provides a quantitative explanation for the data on the Crabtree effect in yeast and the Warburg effect in cancer cells (Bartman et al., 2023; Shen et al., 2024).

Results

Coarse-grained model

Based on the topology of the metabolic network (Neidhardt et al., 1990; Nelson and Cox, 2008) (see Figure 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. Specifically, each pool is designated according to its entry point (see Figure 1A and Appendix 2.2 for details): 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). Pools a1 and a2 are also combined as Pool a due to the joint synthesis of precursors. Then, the metabolic network for Group A carbon source utilization (see Figure 1A) can be coarse-grained into a model shown in Figure 1B (see Appendix 3.1 for details), where node A represents an arbitrary carbon source of Group A. Evidently, Figure 1B is topologically identical to Figure 1A. Each coarse-grained arrow in Figure 1B represents a stoichiometric flux Ji, which delivers carbon flux and may be accompanied by energy consumption or biogenesis (e.g. J1, Ja1; see Figure 1A–B and Appendix 1—figure 1A).

Model and results of overflow metabolism in E. coli.

(A) The central metabolic network of carbon source utilization. The Group A carbon sources (Wang et al., 2019) are labeled with green squares. (B) Coarse-grained model for Group A carbon source utilization. (C) Model predictions (see Equations S47 and S160) and experimental results (Basan et al., 2015; Holms, 1996) of overflow metabolism, covering the data for all the Group A carbon sources shown in (A). (D) Growth rate dependence of respiration and fermentation fluxes (see Equations S47 and S160). (E) The proteome efficiencies for energy biogenesis in the respiration and fermentation pathways vary with growth rate as functions of the nutrient quality of a Group A carbon source (see Equations S31 and S36). See Appendices 9 and 11 for model parameter settings and experimental data sources (Basan et al., 2015; Holms, 1996; Hui et al., 2015) for Figures 14 of E. coli.

In fact, the stoichiometric flux Ji scales with the cell population. For comparison with experiments, we define the normalized flux Ji(N)Jim0/m0McarbonMcarbon, which can be regarded as the flux per unit of biomass (the superscript ‘(N)’ stands for normalized; see Appendix 2.3–2.4 for details). Here, Mcarbon represents the carbon mass of the cell population, and m0 is the weighted average carbon mass of metabolite molecules at the entry of precursor pools (see Equation S17). Then, the cell growth rate λ can be represented by the total outflow of the normalized fluxes: λ=ia1,a2,b,c,dJi(N) (see Appendix 2.4). The normalized fluxes of respiration and fermentation are Jr(N)J4(N) and Jf(N)J6(N), respectively (see Figure 1A and B). In practice, each Ji(N) 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 Ji(N)=ϕiκi. We take the Michaelis-Menten form for the enzyme kinetics (Nelson and Cox, 2008), and then κiki[Si][Si]+Ki (see Equation S12 and Appendix 2.4 for details), where [Si] is the concentration of substrate Si, and Ki is the Michaelis constant. For each intermediate node and reaction along the pathway (e.g. node M1 in Ja1), the substrate quality κi can be approximated as a constant (see Appendix 2.5): κiki[Si][Si]+Kiki, where [Si]Ki generally holds true in bacteria (Bennett et al., 2009; Park et al., 2016). However, the nutrient quality κA is a variable that depends on the nutrient type and concentration of a Group A carbon source (see Equation S27).

Generally, there are three independent fates for a Group A carbon source in the metabolic network (Chen and Nielsen, 2019): fermentation, respiration, and biomass generation (see Appendix 1—figure 1C-E). Each draws a distinct proteome fraction of ϕf, ϕr and ϕBM, with no overlap between them (see Appendix 3.1). The net effect of the first two fates is energy biogenesis, while the last one generates precursors for biomass, accompanied by energy biogenesis. By applying the proteomic constraint that there is a maximum fraction, ϕmax, for proteome allocation: ϕmax0.48 (Scott et al., 2010), we have:

(1) ϕf+ϕr+ϕBM=ϕmax.

In fact, Equation 1 is equivalent to ϕR+ϕA+j=16ϕj+ia1,a2,b,c,dϕi=ϕmax (see Appendix 3.1 for derivation details), where ϕR and ϕA represent the proteomic mass fractions of the active ribosome-affiliated proteins and the cargo proteins responsible for the uptake of the Group A carbon source, respectively. During cell proliferation, ribosomes serve as the factories for protein synthesis and are primarily composed of proteins (Neidhardt et al., 1990; Nelson and Cox, 2008), while other biomass components, such as RNA, are optimally produced (Kostinski and Reuveni, 2020) in accordance with the growth rate determined by protein synthesis. Thus, the cell growth rate is proportional to ϕRλ=ϕRκt, where κt is a parameter set by the translation rate (Scott et al., 2010) (see Appendix 2.1 for details), which can be approximated as a constant within the growth rate range of interest (Dai et al., 2017).

For balanced cell growth in bacteria, the energy demand JE, expressed as the stoichiometric energy flux in ATP, is generally proportional to the biomass production rate (Ebenhöh et al., 2024), since the proportion of maintenance energy is roughly negligible (Locasale and Cantley, 2010) (see Appendix 10 for the cases of yeast and tumor cells). Thus, the normalized flux of energy demand in ATP, denoted as JE(N), representing the energy demand per unit of biomass, is proportional to the growth rate λ (see Appendix 3.1 for details):

(2) JE(N)=ηEλ,

where ηE is an energy coefficient (see Equations S25 and S26 for details). By converting all energy currencies (such as NADH, FADH2, etc.) into ATP, the normalized energy fluxes for respiration and fermentation are given by Jr(E)=βr(A)Jr(N)/2 and Jf(E)=βf(A)Jf(N)/2, where βr(A) and βf(A) are the stoichiometric coefficients of ATP production per glucose in each pathway (see Appendix 1—figure 1C-E and Appendix 3.1 for details). The denominator coefficient of ‘2’ is derived from the stoichiometry of the coarse-grained reaction M12M2 (see Figure 1A and B). Applying the criteria of flux balance (i.e. mass conservation; see Appendix 2.3) at each intermediate node (Mi, i = 1, …, 5) and precursor pool (Pool i, i= a1, a2, b, c, d), along with the constraints of proteome allocation (see Equation 1) and energy demand (see Equation 2), we obtain the relations between normalized energy fluxes and growth rate for a given nutrient condition with a fixed κA (see Appendix 3.1 for details):

(3) {Jr(E)+Jf(E)=φλ,Jr(E)εr+Jf(E)εf=ϕmaxψλ,

where φ is a constant coefficient primarily determined by the coefficient ηE (see Equation S33), and φλ represents the normalized flux of energy demand, excluding energy biogenesis from the biomass synthesis pathway. The coefficients ψ, εr, and εf are functions of κA, such that their values are highly dependent on nutrient conditions. ψ1 denotes the proteome efficiency for biomass generation in the biomass synthesis pathway (see Equation S32), defined as ψ1λ/ϕBM (see Appendix 3.1). εr and εf represent the proteome efficiencies for energy biogenesis in the respiration and fermentation pathways, respectively, defined as the normalized energy fluxes expressed in ATP generated per proteomic mass fraction, with εrJr(E)/ϕr and εfJf(E)/ϕf. Hence,

(4) {εr=βr(A)1/κA+1/κr(A),εf=βf(A)1/κA+1/κf(A),

where both κr(A) and κf(A) are composite parameters that can be approximated as constants, with 1/κr(A)1/κ1+2/κ2+2/κ3+2/κ4 and 1/κf(A)1/κ1+2/κ2+2/κ6 (see Appendices 2.5 and 3.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; van Hoek et al., 1998) is exemplified by the experimental data (Basan et al., 2015) presented in Figure 1C, where the fermentation flux exhibits a threshold-analog dependence on the growth rate λ. It is well established that respiration is significantly more efficient than fermentation in terms of energy biogenesis per unit of carbon (i.e. βr(A)>βf(A)) (Nelson and Cox, 2008; Vander Heiden et al., 2009). Then, why do cells bother to use the seemingly wasteful fermentation pathway? We proceed to address this issue by applying optimal protein allocation (Scott et al., 2010; Wang et al., 2019) within the framework of optimal growth.

For cell proliferation in a given nutrient condition (i.e. with a fixed κA), the values of εr, εf, and ψ are determined (see Equations 4 and S32). However, the growth rate λ can be influenced by protein allocation between respiration and fermentation, specifically ϕr and ϕf, according to the governing equation (Equation 3). If εr>εf, that is, if the proteome efficiency in respiration is higher than that in fermentation, then λ=ϕmaxJf(E)(1/εf1/εr)ψ+φ/εrϕmaxψ+φ/εr. The optimal growth strategy is ϕf=Jf(E)=0, meaning that the cell exclusively uses respiration. Conversely, if εf>εr, then ϕr=Jr(E)=0 is optimal, and the cell solely uses fermentation. In either case, the choice between respiration and fermentation for growth optimization is determined by comparing their proteome efficiencies.

In practice, both proteome efficiencies εr and εf are functions of nutrient quality κA, which can be significantly influenced by the nutrient type and concentration of the carbon source (see Equations 4 and S27). Therefore, the optimal growth strategy may vary depending on the nutrient conditions. In nutrient-poor conditions where κAκr(A) and κAκf(A), the proteome efficiencies can be approximated by εrβr(A)κA and εfβf(A)κA (see Equation 4), and hence εr(κA)>εf(κA) (since βr(A)>βf(A)), meaning that the proteome efficiency of respiration is higher than that of fermentation under these conditions. In contrast, in rich media, using parameters for κi derived from in vivo/in vitro experimental data for E. coli (see Appendix 1—table 1, Appendix 1—table 2 and Appendix 7.1–7.2), we obtain εr(κglucose(ST))<εf(κglucose(ST)) with Equation 4 (see also Equations S39-S40), where κglucose(ST) represents the substrate quality of glucose at saturated concentration (abbreviated as ‘ST’ in the superscript). This indicates that the proteome efficiency in fermentation is higher than that in respiration for bacteria in rich media. Indeed, recent studies have validated that the measured proteome efficiency in fermentation is higher than in respiration for E. coli in lactose at saturated concentration (Basan et al., 2015), i.e., εr(κlactose(ST))<εf(κlactose(ST)). In Figure 1E, we present the growth rate dependence of proteome efficiencies εr and εf in a three-dimensional (3D) format using the collected data shown in Appendix 1—table 1, where εr, εf and the growth rate λ all vary as functions of nutrient quality κA. Furthermore, the ratio Δ (defined as Δ(κA)εf(κA)/εr(κA)) is a monotonically increasing function of κA, and there exists a critical value of κA (denoted as κA(C); see Appendix 3.2 for details) satisfying Δ(κA(C))=1. Below κA(C), where the nutrient is poorer and the cell grows slowly, the proteome efficiency of fermentation is lower than that of respiration (i.e. εf<εr), hence respiration is the optimal choice (with λ=ϕmax(ψ+φ/εr)1). Above κA(C), where the nutrient is richer and the cell grows faster, fermentation is more efficient than respiration in terms of proteome efficiency (i.e. εf>εr) and becomes the optimal growth strategy (with λ=ϕmax(ψ+φ/εf)1). This 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 under nutrient variation (with significantly varying κA) is given by Jf(E)=φλθ(λλC), where ‘θ’ represents the Heaviside step function, and λC denotes the critical growth rate corresponding to the nutrient condition with nutrient quality κA(C) (i.e. λCλ(κA(C))). Similarly, the growth rate dependence of respiration flux is Jr(E)=φλ[1θ(λλC)]. These digital response outcomes are consistent with the numerical simulation findings of Molenaar et al., 2009. However, they are clearly incompatible with the threshold-analog response observed in the standard picture of overflow metabolism (Basan et al., 2015; Holms, 1996; Meyer et al., 1984; Nanchen et al., 2006; van Hoek et al., 1998).

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; Yaginuma et al., 2014; Zhang et al., 2018) and tumor cells (Duraj et al., 2021; Shibao et al., 2018; Hanahan and Weinberg, 2011; Hensley et al., 2016). In the context of the Warburg effect or overflow metabolism, 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; Shibao et al., 2018; Hensley et al., 2016; Nikolic et al., 2013). Motivated by the observation that the turnover number (kcat 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 kcat values of metabolic enzymes (García-Contreras et al., 2012). Therefore, within a cell population, there is a distribution of kcat values for a catalytic enzyme, commonly referred to as extrinsic noise (Elowitz et al., 2002). For simplicity, we assume that the kcat values for each enzyme follow a Gaussian distribution. Consequently, the proteome efficiencies εr and εf, which are crucial for determining the choice between respiration and fermentation, also follow Gaussian distributions (see Appendix 8 for details). This variability leads to diverse distributions of single-cell growth rates across different carbon sources (see Equations S155-S157 and S163-S165), which has been fully verified by recent experiments using isogenic E. coli at single-cell resolution (Wallden et al., 2016; see Appendix 1—figure 2B). Accordingly, the critical growth rate λC is expected to follow a Gaussian distribution N(μλC,σλC2) within a cell population (see Appendix 8 for details), where μλC is approximated by the deterministic result of λC (Equation S43). Assuming the coefficient of variation (CV) of λC is σλC/μC=12%, or equivalently that the CV for the catalytic rate of each metabolic enzyme is 25%, we derive the growth rate dependence of fermentation and respiration fluxes (see Appendix 3.3 for details):

(5) {Jf(N)(λ)=φλβf(A)[erf(λμλC2σλC)+1],Jr(N)(λ)=φλβr(A)[1erf(λμλC2σλC)],

where ‘erf’ represents the error function. The fermentation flux exhibits a threshold-analog relation with the growth rate (the red curves in Figures 1C–D3B, D and F), while the respiration flux (the blue curve in Figure 1D) decreases as the fermentation flux increases. In Figure 1C–D, we observe that the model results (see Equation 5 and Appendix 9 for details; parameters are set based on the experimental data shown in Appendix 1—table 1) quantitatively agree with the experimental data from E. coli (Basan et al., 2015; Holms, 1996). The fermentation flux is represented by the acetate secretion rate Jactate(M)=2Jf(N), and the respiration flux is exemplified by the carbon dioxide flux JCO2,r(M)=6Jr(N) (the superscript ‘(M)’ represents the measurable flux in the unit of mM/OD600/h; see Appendix 9.1 for details). By incorporating cell heterogeneity, our model of optimal protein allocation quantitatively explains overflow metabolism.

Testing the model through perturbations

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

First, we consider the proteomic perturbation caused by overexpression of useless proteins encoded by the lacZ gene (i.e. ϕZ perturbation) in E. coli. The net effect of the ϕZ perturbation is that the maximum fraction of the proteome available for resource allocation changes from ϕmax to ϕmaxϕZ (Basan et al., 2015), where ϕZ is the proteomic 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)2), where the CV of λC(ϕZ) remains unchanged. Consequently, the growth rate dependence of fermentation flux changes to Jf(N)=φλβf(A)[erf(λμλC(ϕZ)2σλC(ϕZ))+1] (see Appendix 4 for model perturbation results regarding respiration flux), where both the growth rate λ(κA,ϕZ) and the normalized fermentation flux Jf(N)(κA,ϕZ) are bivariate functions of κA and ϕZ (see Equations S49, S56 and S57). For each degree of LacZ expression (with fixed ϕZ), similar to wild-type strains, the fermentation flux exhibits a threshold-analog response to growth rate as κA varies (see Figure 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/ϕZϕmaxϕ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: Jf(N)=φβf(A)[erf(λ(κA,0)μλC(0)2σλC(0))+1]λ, 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 Figure 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 (see Figure 2A). Experimental data points (Basan et al., 2015) lie right on this surface, which is highly consistent with the model predictions.

Influence of protein overexpression on overflow metabolism in E. coli.

(A) A 3D plot of the relations among fermentation flux, growth rate, and the expression level of useless proteins. In this plot, both the acetate excretion rate and growth rate vary as bivariate functions of the nutrient quality of a Group A carbon source (denoted as κA) and the useless protein expression encoded by lacZ gene (denoted as ϕZ perturbation; see Equations S57 and S160). (B) Growth rate dependence of the acetate excretion rate upon ϕZ perturbation for each fixed nutrient condition (see Equations S58 and S160). (C) Growth rate dependence of the acetate excretion rate as κA varies (see Equations S58 and S160), with each fixed expression level of LacZ.

Next, we study the influence of energy dissipation, which introduces an energy dissipation coefficient w to Equation 2: JE(N)=ηEλ+w. Similarly, the critical growth rate in this case, λC(w), follows a Gaussian distribution N(μλC(w),σλC(w)2) in a cell population. The relation between the growth rate and fermentation flux can be characterized by: Jf(N)=φλ+wβf(A)[erf(λμλC(w)2σλC(w)2)+1] (see Appendix 4.2 for details). In Figure 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.

Influence of energy dissipation, translation inhibition, and carbon source category alteration on overflow metabolism in E. coli.

(A) A 3D plot of the relations among fermentation flux, growth rate, and the energy dissipation coefficient (see Equations S70 and S160). (B) Growth rate dependence of the acetate excretion rate as the nutrient quality κA varies, with each fixed energy dissipation coefficient determined by or fitted from experimental data. (C) A 3D plot of the relations among fermentation flux, growth rate, and the translation efficiency (see Equations 85 and S160). Here, the translation efficiency is adjusted by the dose of chloramphenicol (Cm). (D) Growth rate dependence of the acetate excretion rate as κA varies, with each fixed dose of Cm. (E) Coarse-grained model for pyruvate utilization. (F) The growth rate dependence of fermentation flux in pyruvate (see Equations 105 and S160) significantly differs from that of the Group A carbon sources (see Equations 47 and S160).

We proceed to analyze the impact of translation inhibition with different sub-lethal doses of chloramphenicol on E. coli. This type of perturbation introduces an inhibition coefficient ι to the translation rate, thus turning κt into κt/ι+1. Still, the critical growth rate λC(ι) follows a Gaussian distribution N(μλC(ι),σλC(ι)2), and then, the growth rate dependence of fermentation flux is given by: Jf(N)=φλβf(A)[erf(λμλC(ι)2σλC(ι))+1] (see Appendix 4.3 for details). In Appendix 1—figure 2D and E, we observe that the model predictions are generally consistent with the experimental data (Basan et al., 2015). However, a noticeable systematic discrepancy arises when the translation rate is low. Therefore, we consider maintenance energy, which is typically tiny and generally negligible for bacteria over the growth rate range of interest (Basan et al., 2015; Locasale and Cantley, 2010; Neidhardt, 1996). Encouragingly, by assigning a very small value to the maintenance energy coefficient w0 (where w0=2.5(h1)), the model results for the growth rate-fermentation flux relation Jf(N)=φλ+w0βf(A)[erf(λμλC(ι)2σλC(ι))+1] quantitatively agree with experiments (Basan et al., 2015) (see Figure 3C–D and Appendix 4.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 and Cox, 2008). The coarse-grained model for pyruvate utilization is shown in Figure 3E (see also Figure 1A), which shares identical precursor pools with those for Group A carbon sources, yet has several differences in the coarse-grained reactions. The growth rate dependencies of both the proteome efficiencies (see Appendix 1—figure 2H) and energy fluxes (see Figure 3F) are qualitatively similar to those of Group A carbon source utilization, while there are quantitative differences in the coarse-grained parameters (see Appendices 5.1 and 9 for derivation details). Most notably, the critical growth rate λC(py) and the ATP production per glucose in the fermentation pathway βf(py) for pyruvate utilization are noticeably smaller than those for Group A sources (i.e. λC and βf(A), respectively). Consequently, the growth rate dependence of fermentation flux in pyruvate should present a distinctly different curve from that of Group A carbon sources (see Equations 5 and S105), which is fully validated by experimental results (Holms, 1996; see Figure 3F).

Enzyme allocation under perturbations

As mentioned above, our coarse-grained model is topologically identical to the central metabolic network (see Figure 1A) and can thus predict enzyme allocation for each gene in glycolysis and the TCA cycle (see Appendix 1—figure 1B and Appendix 1—table 1) under various types of perturbations. In Figure 1B, the intermediate nodes M1, M2, M3, M4, and M5 represent G6P, PEP, acetyl-CoA, α-ketoglutarate, and oxaloacetate, respectively. Therefore, ϕ1 and ϕ2 correspond to enzymes involved in glycolysis (or at the junction of glycolysis and the TCA cycle), while ϕ3 and ϕ4 correspond to enzymes in the TCA cycle (see Figure 1A–B and Appendix 3.1).

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). 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 is generally described by a C-line response (Hui et al., 2015; You et al., 2013). Specifically, the genes responsible for digesting carbon compounds exhibit 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 gathered 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 (see Appendix 1—figure 3A and B). This discrepancy cannot be explained by the C-line response. To address this issue, we apply the coarse-grained model described above (see Figure 1B) to calculate the growth rate dependence of enzyme allocation for each ϕi (i=1,2,3,4) using model settings for wild-type strains, with no fitting parameters influencing the shape (see Equations S118-S119 and Appendix 9). In Figure 4A–B and Appendix 1—figure 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 w0=2.5(h1)) 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).

We proceed to analyze the influence of ϕZ perturbation and energy dissipation. In both cases, our model predicts a linear response to growth rate reduction for all genes in either glycolysis or the TCA cycle (see Appendix 6.2–6.3 for details). For ϕZ perturbation, all predicted slopes are positive, and there are no fitting parameters involved (Equations S120-S121). In Figure 4C–D and Appendix 1—figure 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 w0=2.5(h1) as aforementioned). For energy dissipation, however, the predicted slopes of the enzymes corresponding to ϕ4 are 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 1—figure 3K-N, we see that the model results (Equations S127 and S123) are consistent with experiments (Basan et al., 2015).

Relative protein expression of central metabolic enzymes in E. coli under carbon limitation and proteomic perturbation.

(A, C) Relative protein expression of representative genes from glycolysis. (B, D) Relative protein expression of representative genes from the TCA cycle. (A, B) Results of the perturbation through changes in nutrient quality κA (see Equation S119). (C, D) Results of proteomic perturbation via varied levels of expression of the useless protein LacZ (i.e. ϕZ perturbation; see Equation S121).

Explanation of the Crabtree effect in yeast and the Warburg effect in cancer cells

We proceed to apply our model to explain the Crabtree effect in yeast (Bagamery et al., 2020; De Deken, 1966; Shen et al., 2024) and the Warburg effect in tumors (Bartman et al., 2023; Duraj et al., 2021; Hanahan and Weinberg, 2011; Shen et al., 2024; Vander Heiden et al., 2009) with slight modifications using the optimal growth principle combined with cell heterogeneity (see Appendix 10 and Appendix 1—figure 5). For yeast and tumors, similar to the case of E. coli, the proteome efficiencies εr and εf are both increasing functions of nutrient quality κA (see Equation S170). Under poor nutrient conditions (i.e. κA is small), the proteome efficiency in respiration is higher than that in fermentation: εr>εf (see Equations S174-S175), making respiration the optimal choice for growth optimization (see Equation S171). Conversely, when nutrients are abundant and εf>εr, aerobic glycolysis (i.e. fermentation) becomes the optimal growth strategy (see Equation S172). Further combination with cell heterogeneity results in the standard picture of overflow metabolism, which has indeed been observed in yeast (van Hoek et al., 1998). However, it remains challenging to tune the growth rate of cancer cells in vivo.

Recently, Shen et al., 2024 discovered that the proteome efficiency measured at the cell population level in respiration (i.e. εr ; where ‘’ denotes the population average) is higher than that in fermentation (i.e. εf) for many yeast and cancer cells, despite the presence of fermentation fluxes through aerobic glycolysis. Evidently, this finding (Shen et al., 2024) contradicts prevalent explanations (Basan et al., 2015; Chen and Nielsen, 2019), which hold that overflow metabolism arises because the proteome efficiency in fermentation is consistently higher than in respiration. Nevertheless, our model may resolve this puzzle due to the incorporation of two important features. First, our model predicts that the proteome efficiency in respiration is larger than that in fermentation when nutrient quality is low (see Equations S174-S175). Second, and crucially, by accounting for cell heterogeneity, our model allows a proportion of cells to have a higher proteome efficiency in fermentation than in respiration, even when the overall proteome efficiency in respiration at the cell population level is greater than that in fermentation (i.e. εr>εf).

To compare our model results quantitatively with experimental data on yeast and tumors (Shen et al., 2024), we define PrfJf(E)Jf(E)+Jr(E) as the fraction of ATP produced through fermentation. To account for cell heterogeneity, we apply Gaussian distributions to enzyme turnover numbers, as described above. This yields the relationship between Prf (i.e.Jf(E)Jf(E)+Jr(E)) and εr and εf through derivations (see Equations S180-S190 and Appendix 10 for details):

(6) Jf(E)Jf(E)+Jr(E)=12[erf(1εr/εf2χεr2+χεf2(εr/εf)2)+1],

where χεr and χεf represent the CVs of proteome efficiencies εr and εf, respectively. Due to the higher levels of cell heterogeneity in yeast (Bagamery et al., 2020) and cancer cells (Duraj et al., 2021; Shibao et al., 2018; Hanahan and Weinberg, 2011; Hensley et al., 2016), the CVs of εr and εf (i.e. χεr and χεf) in these cells are expected to be significantly higher than those in E. coli, although their precise values are unknown. The values for the variables shown in Equation 6 can be obtained from experiments. Therefore, we plot the theoretical results from Equation 6 using χεr and χεf values of 0.25, 0.40, and 0.58 to compare with experimental data from yeast and in vivo mouse tumors (Bartman et al., 2023; Shen et al., 2024). As shown in Figure 5A–B, the theoretical results with χεr=χεf=0.58 align quantitatively with the experimental data (Bartman et al., 2023; Shen et al., 2024) on both logarithmic and linear scales, demonstrating that our model has the potential to quantitatively explain the Crabtree effect in yeast and the Warburg effect in cancer cells.

Model comparison with data on the Crabtree effect in yeast and the Warburg effect in tumors.

(A) A linear scale representation on the y-axis. (B) A log scale representation on the y-axis. In (A–B), εr and εf represent the population averages of εr and εf, while χεr and χεf are the coefficients of variation (CVs) of εr and εfεr/εf represents the ratio of proteome efficiency between respiration and fermentation at the population-averaged level, while Jf(E)/(Jf(E)+Jr(E)) stands for the fraction of energy flux generated by the fermentation pathway (see Equation 6). The data for yeast in batch culture and chemostat were calculated from experimental data of S. cerevisiae and I. orientalis (Shen et al., 2024). The data for mouse tumors were calculated from in vivo experimental data of pancreatic ductal adenocarcinoma (PDAC) and leukemic spleen of mice (Bartman et al., 2023; Shen et al., 2024). See Appendix 11 for detailed information on the experimental data sources (Bartman et al., 2023; Shen et al., 2024).

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 over the past century (Basan et al., 2015; Chen and Nielsen, 2019; Majewski and Domach, 1990; Molenaar et al., 2009; Niebel et al., 2019; Peebo et al., 2015; Pfeiffer et al., 2001; Shlomi et al., 2011; Vander Heiden et al., 2009; Varma and Palsson, 1994; Vazquez et al., 2010; Zhuang et al., 2011), contradictions persist (Shen et al., 2024), leaving the origin and function of this phenomenon unclear (DeBerardinis and Chandel, 2020; Hanahan and Weinberg, 2011; Vander Heiden et al., 2009). In this study, we use E. coli as a typical example and demonstrate that overflow metabolism can be understood through optimal protein allocation combined with cell heterogeneity. Under nutrient-poor conditions, the proteome efficiency of respiration is higher than that of fermentation (see Figure 1E), and thus the cell uses respiration to optimize growth. In rich media, however, the proteome efficiency of fermentation increases more rapidly and surpasses that of respiration (see Figure 1E), leading the cell to adopt fermentation as the optimal growth strategy. 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 threshold-analog response (Basan et al., 2015; Holms, 1996) in overflow metabolism (see Figure 1C). Furthermore, it quantitatively explains the data on the Crabtree effect in yeast and the Warburg effect in cancer cells (Bartman et al., 2023; Shen et al., 2024).

Mechanistically, the optimal growth strategy for the binary choice between respiration and fermentation can be facilitated by the direct sensing and comparison of proteome efficiencies between the two pathways (see Appendix 3.4). A growing body of evidence suggests that the cyclic AMP (cAMP)-cAMP receptor protein (CRP) system plays a crucial role in sensing proteome efficiency and executing the optimal strategy (Basan et al., 2015; Towbin et al., 2017; Valgepea et al., 2010; Wehrens et al., 2023). However, it has also been suggested that the cAMP-CRP system alone is insufficient, and that additional regulators remain to be identified to fully elucidate this mechanism (Basan et al., 2015; Valgepea et al., 2010). Furthermore, since the binary choice between respiration and fermentation is driven by the comparison of proteome efficiencies, the optimal growth principle in our model can be relaxed to the case where efficient protein allocation is required only for enzymes, rather than ribosomes. This allows our model to remain applicable under suboptimal growth conditions (see Appendix 3.4 for details), where recent experimental studies have shown that the inactive portion of ribosomes (i.e. ribosomes not bound to mRNAs) may vary with culturing conditions (Dai et al., 2017; Li et al., 2018) and between individual cells within the same culture (Pavlou et al., 2025), despite an overall trend toward growth optimization.

In existing rationales (Basan et al., 2015; Chen and Nielsen, 2019; Majewski and Domach, 1990; Shlomi et al., 2011; Varma and Palsson, 1994; Vazquez et al., 2010; Vazquez and Oltvai, 2016), the standard picture of overflow metabolism (Basan et al., 2015; Holms, 1996; Meyer et al., 1984; Nanchen et al., 2006; van Hoek et al., 1998) has primarily been illustrated by a threshold-linear response, which largely relies on the assumption that cells optimize their growth rate for a given rate of carbon influx under each nutrient condition (or similar equivalents; see Appendix 7.3). However, in practice, for microbes or tumor cells grown in vitro or in vivo, the given factors are the identity and concentration of the nutrient (Molenaar et al., 2009; Scott et al., 2010; Wang et al., 2019), rather than the rate of carbon influx. Additionally, prevalent explanations (Basan et al., 2015; Chen and Nielsen, 2019) suggest that overflow metabolism originates from the proteome efficiency in fermentation always being higher than that in respiration (see Appendix 7.3 for details). While it has been observed in E. coli that proteome efficiency in fermentation is higher than that in respiration for cells cultured in lactose at saturated concentration (Basan et al., 2015), Shen et al., 2024 reported that for many yeast and cancer cells, the proteome efficiency in fermentation is noticeably lower than that in respiration, despite the presence of aerobic glycolytic fermentation flux. This observation (Shen et al., 2024) evidently contradicts the prevalent explanations (Basan et al., 2015; Chen and Nielsen, 2019). Our model resolves this puzzle by significantly differing from existing rationales in its optimization principle, where we optimize cell growth rate purely through protein allocation without imposing a special constraint on carbon influx (see Appendix 7.3 for details). More importantly, our model incorporates cell heterogeneity, which is crucial for both explaining the threshold-analog response in overflow metabolism and for resolving this puzzle raised by Shen et al., 2024.

In the homogeneous case, the optimal growth strategy for growth rate dependent fermentation flux results in a digital response (see Equation S44), corresponding to an elementary flux mode (Müller et al., 2014; Wortel et al., 2014), which aligns with the numerical study by Molenaar et al., 2009 but is incompatible with the standard picture of overflow metabolism (Basan et al., 2015; Holms, 1996; Meyer et al., 1984; Nanchen et al., 2006; van Hoek et al., 1998). Furthermore, in this case, cells would not generate fermentation flux if the proteome efficiency in fermentation were lower than that in respiration, under the optimal growth framework. By incorporating heterogeneity in enzyme catalytic rates (Davidi et al., 2016; García-Contreras et al., 2012), the critical growth rate (i.e. threshold) shifts from a single value to a Gaussian distribution (see Equation S45 and Appendix 8 for details; see also Appendix 1—figure 4) across a cell population, thereby turning a digital response into the threshold-analog response observed in overflow metabolism (see Figure 1C). Moreover, cell heterogeneity allows a fraction of cells to possess a larger proteome efficiency in fermentation than in respiration despite the overall proteome efficiency in respiration at the cell population level is higher than in fermentation. This mechanism facilitates the fermentation flux in yeast and cancer cells observed by Shen et al., 2024 (see Figure 5A–B).

Our model results, based on cell heterogeneity, are further supported by observed distributions of single-cell growth rates in E. coli (Wallden et al., 2016) (see Appendix 1—figure 2B), as well as by experiments involving various types of perturbations (Basan et al., 2015; Holms, 1996; Hui et al., 2015), both in terms of acetate secretion patterns and gene expression in the central metabolic network (see Figures 24, Appendix 1—figures 2D and E and 3). Furthermore, the heterogeneity patterns predicted by our model for fermentation and respiration modes in an isogenic cell population under the same culturing conditions are highly consistent with the non-genetic heterogeneity observed in single-cell experiments with E. coli (Nikolic et al., 2013) and S. cerevisiae (Bagamery et al., 2020), and align with experiments on intra-tumor heterogeneity in glioblastoma (Duraj et al., 2021; Shibao et al., 2018). Finally, our model can be broadly applied to address heterogeneity-related challenges in metabolism on a quantitative basis, including diverse metabolic strategies of cells in various 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).

Appendix 1

Appendix 1—table 1
Molecular weight (MW) and in vivo/in vitro kcat data for E. coli.
No.*ReactionEnzymeGene nameECMW (kDa)In vitro kcat (s-1)ReferencesIn vivokcat (s-1)Selected kcat (s-1)
J1Glucose-6P ↔ Fructose-6PGlucose-6-phosphate isomerasepgiEC:5.3.1.91.2×1022.6×102PMID:7004378;
DOI:https://doi.org/10.1016/j.ijms.2004.09.017
8.7×1028.7×102
Fructose-6P → Fructose-1,6PPhosphofructokin-asepfkAEC:2.7.1.111.4×1024.4×102PMID:6218375; 702261.7×1031.7×103
Fructose-1,6P ↔ Glyceraldehyde 3-phosphate+Dihydroxyacetone phosphateFructose-bisphosphate aldolasefbaAEC:4.1.2.137.8×101.4×10PMID:8939754; 155316271.6×1021.6×102
Dihydroxyacetone phosphate ↔ Glyceraldehyde 3-phosphateTriosephosphate IsomerasetpiAEC:5.3.1.15.4×104.3×102PMID:3887397; 60928572.7×1022.7×102
Glyceraldehyde 3-phosphate ↔ 1,3-BisphosphoglycerateGlyceraldehyde-3-phosphate dehydrogenasegapAEC:1.2.1.121.4×1029.5×10PMID:4932978; 22009291.5×1021.5×102
1,3-Bisphosphoglycerate ↔ 3-PhosphoglyceratePhosphoglycerate kinasepgkEC:2.7.2.34.4×103.5×102PMID:367367; 1662741.9×1021.9×102
3-Phosphoglycerate ↔ 2-PhosphoglyceratePhosphoglycerate mutasegpmAEC:5.4.2.114.9×103.3×102PMID:104378014.5×1024.5×102
2-Phosphoglycerate ↔ PhosphoenolpyruvateEnolaseenoEC:4.2.1.119.0×102.2×102PMID:1094232; 49423261.7×1021.7×102
J2Phosphoenolpyruvate → PyruvatePyruvate kinasepykFEC:2.7.1.402.4×1025.0×102PMID:67598521.6×1031.6×103
Pyruvate → Acetyl-CoAPyruvate dehydrogenaseaceEEC:1.2.4.11.0×1021.2×102PMID:230884223.4×1023.4×102
J3Oxaloacetate +Acetyl CoA → CitrateCitrate synthasegltAEC:2.3.3.19.7×102.4×102PMID:4900996; 239543057.1×107.1×10
Citrate ↔ IsocitrateAconitate hydrataseacnBEC:4.2.1.39.4×107.0×10PMID:15963579; 159635796.3×106.3×10
Isocitrate→ α-KetoglutarateIsocitrate dehydrogenaseicdEC:1.1.1.429.5×102.0×102PMID:8141; 36923; 22009293.3×103.3×10
J4α-Ketoglutarate → Succinyl-CoAα-Ketoglutarate dehydrogenase complex E1 componentsuc A suc BEC:1.2.4.2, EC:2.3.1.611.9×1021.5×102PMID:6380583; 45886791.3×1021.3×102
Succinyl-CoA ↔ SuccinateSuccinyl-CoA synthetasesuc C suc DEC:6.2.1.51.6×1029.1×10PMID:53381301.0×1021.0×102
Succinate → FumarateSuccinate dehydrogenasesdh A sdh BEC:1.3.5.11.0×1021.1×102PMID:4334990; 164842321.1×1021.1×102
Fumarate ↔ MalateFumarasefumAEC:4.2.1.22.0×1021.2×103PMID:3282546; 120214534.9×1024.9×102
Malate ↔ OxaloacetateMalate dehydrogenasemdhEC:1.1.1.376.1×105.5×102doi:https://doi.org/10.1016/0076-6879(69)13029-36.6×106.6×10
J5Phosphoenolpyruvate →OxaloacetatePhosphoenolpyru-vate carboxylaseppcEC:4.1.1.314.0×1021.5×102PMID:9927652; 4932977/1.5×102
J6Acetyl-CoA ↔ Acetyl phosphatePhosphate acetyltransferaseptaEC:2.3.1.87.7×103.0×10PMID:202363193.7×1023.7×102
Acetyl phosphate↔ AcetateAcetate kinaseackAEC:2.7.2.14.3×103.6×103EcoCyc:
EG10027;
PMID:24801996
3.3×1023.3×102
Acetate (intracellular) ↔ Acetate (extracellular)Acetate transporteractP/2×104.7×102PMID:31405984 (Estimated)/4.7×102
J7Pyruvate → PhosphoenolpyruvatePyruvate, water dikinaseppsAEC:2.7.9.22.5×1023.5×10PMID:4319237/3.5×10
JAGlucose-6P (extracellular) → Glucose-6P (intracellular)Glucose-6-phosphate transporterUhpT/5×102×102PMID:3283129; 2197272;
20018695
(Estimated)
/2×102
Glucose (extracellular) → Glucose-6PGlucose-specific PTS enzymeptsGEC:
2.7.1.199
5×101×102PMID:9575173; 20018695; 12146972/1×102
Lactose (extracellular) → Lactose (intracellular)Lactose transporterlacY/4.6×106×10PMID:6444453; 20018695/6×10
Lactose →Glucose +Galactoseβ-galactosidaselacZEC:3.2.1.234.6×1026.4×102PMID:8008071;
23011886
(Estimated)
/6.4×102
JpyPyruvate (extracellular) → Pyruvate (intracellular)Pyruvate transporterbtsT CstA/8×106×10PMID:20018695; 33260635;
EcoCyc: G7942; EG10167
(Estimated)
/6×10
  1. *

    The classification of Ji follows the coarse-grained models shown in Figures 1B and 3E.

  2. In vivo kcat values were obtained using the experimental data shown in Appendix 1—table 2, combined with Equations S134-S135.

  3. See Appendix 1—figure 1B for additional genes that may play a secondary role.

Appendix 1—table 2
Proteome and flux data (Basan et al., 2015) used to calculate the in vivo kcat of E. coli.
Culture 1Culture 2Culture 3Culture 4
Growth rate λ (h–1)*0.820.870.971.03
Jacetate (mM OD600–1 h–1)0.391.182.682.84
JCO2, r (mM OD600–1 h–1) 7.446.054.303.04
Gene nameProteomic mass fractions obtained using absolute abundance (ϕi)
pgi0.09%0.09%0.10%0.11%
pfkA0.06%0.06%0.06%0.06%
fbaA0.32%0.35%0.35%0.39%
tpiA0.12%0.15%0.13%0.18%
gapA1.19%1.29%1.33%1.47%
pgk0.30%0.31%0.32%0.36%
gpmA0.15%0.15%0.15%0.16%
eno0.63%0.70%0.75%0.83%
pykF0.15%0.15%0.18%0.21%
aceE0.30%0.32%0.34%0.41%
gltA0.88%0.80%0.61%0.48%
acnB0.92%0.84%0.66%0.57%
icd1.55%1.55%1.31%1.39%
suc A suc B0.71%0.75%0.64%0.55%
suc C suc D0.88%0.84%0.66%0.52%
sdh A sdh B0.49%0.45%0.42%0.35%
fumA0.24%0.21%0.17%0.13%
mdh0.45%0.45%0.41%0.39%
pta0.10%0.10%0.10%0.10%
ackA0.06%0.07%0.06%0.06%
  1. *

    For calibration purposes, a factor of 1.03/0.97 was multiplied by the reference data (Basan et al., 2015).

  2. For calibration purposes, a factor of 2.84/3.24 was multiplied by the reference data (Basan et al., 2015).

  3. Here, (1.03, 2.84) and (0.97, 3.24) are both the data points for (λ h-1, Jacetate mM OD600-1 h-1) for E. coli strain NCM3722 cultured with lactose in the same reference (Basan et al., 2015). The former is specified in the source data of the reference’s figure 1 (Basan et al., 2015), while the latter is recorded in the reference’s extended data figure 3a (Basan et al., 2015). With the calibrations above, the data for the Jacetate(M)λ relation shown here align with the curve depicted in Figure 1C.

Appendix 1—table 3
Illustrations of symbols in this manuscript.
SymbolsIllustrations/DefinitionsModel variable/parameter settings for E. coli*
A (in the figures)A Group A carbon source joining the metabolic network from the upper part of glycolysis.NA
Mi (in the figures)A metabolite in the metabolic network that serve as intermediate node.NA
Ji (in the figures)The stoichiometric flux delivering carbon flux, an extensive variable; see Equation S7.see Equations S7-S8.
ri (in the figures)The mass fraction of carbon flux drawn from a precursor pool.ra1=24%; ra2=24%; rb = 28%; rc = 12%; rd = 12%
(Nelson and Cox, 2008).
λGrowth rate of the cell population; see Equation S36 for the optimal model solution.see Equations S4 and S36.
Jr, JfJr and Jf are stoichiometric fluxes of respiration and fermentation, extensive variables.Jr = J4; Jf = J6 (see Equation S22)
m0 The weighted average carbon mass of metabolite molecules at the entrance of precursor pools.See Equation S17.
McarbonThe carbon mass of the cell population, an extensive variable.NA
MproteinThe protein mass of the cell population; an extensive variable.NA
MQ(P), MR(P), MC(P)The mass of Q-class, R-class, or C-class proteome.See Equation S2.
fQ, fR, fCThe ribosome allocation fraction for protein synthesis of Q-class, R-class, or C-class.fQ=ϕQ.
mAAThe average molecular weight of amino acids.A reducible parameter for the results.
kTTranslation speed of ribosomes.kT=20.1 aa/s (Scott et al., 2010).
ϕQ, ϕR, ϕCThe mass fraction of Q-class, R-class, or C-class proteome; see Appendix 2.1.ϕQ=52% (Scott et al., 2010).
ϕmaxThe maximum proteomic mass fraction of proteome allocation for fermentation, respiration, and biomass generation, with ϕmax1ϕQ.ϕmax=48% (Scott et al., 2010).
mRThe protein mass of a single ribosome.mR=7336mAA
(Neidhardt et al., 1990).
VcellThe cell volume of the cell population (the ‘big cell’); an extensive variable.NA
NR, Mrp(P)The number or the total protein mass of ribosomes in the big cell; extensive variables.NA
ςThe ratio of the mass of R-class proteome to the protein mass of ribosomes: ςMR(P)/Mrp(P).ς=1.67 (Scott et al., 2010).
[Ei], [Si]The concentration of enzyme Ei or substrate Si; intensive variables.NA
ai, di, bi, ciai and di are reaction parameters; bi and ci are stoichiometric coefficients. See Appendix 2.3.NA
KiThe Michaelis constant, defined as Ki≡(di+kicat)/ ai.Obtainable from Bennett et al., 2009, yet unused in practice since [Si]>Ki
(see Appendix 2.5).
viThe reaction rate per volume of a biochemical reaction catalyzed by Ei; an intensive variable.See Equation S6.
NEi, MEiThe copy number or the total weight enzyme Ei in the cell population; extensive variables.NEi=Vcell[Ei];
MEi=NEimEi.
mcarbonThe mass of a carbon atom.mcarbon=12NAvogadrog, where g represents gram and NAvogadro is the Avogadro constant.
ΦiThe enzyme cost of all Ei molecules in the cell population; an extensive variable.ΦiNEinEi.
ξiξi is defined such that ξi=Ji/Φi.ξikicatnEi[Si][Si]+Ki.
Ji(N)The normalized flux, i.e., flux per unit of biomass; an intensive variable§Ji(N)Jim0/Mcarbon
see Equations S15-S16.
Jr(N), Jf(N)Jr(N) and Jf(N) are the normalized fluxes of respiration and fermentation, intensive variables.Jr(N)=J4(N); Jf(N)=J6(N).
NEPicarbonThe number of carbon atoms in the entry point metabolite molecule of Precursor Pool i.NEPa1carbon=6; NEPa2carbon=3; NEPbcarbon=3; NEPccarbon=5; NEPdcarbon=4 (Nelson and Cox, 2008).
kcat, kicatThe turnover number of a catalytic enzyme.See Appendix 1—table 1.
mEi, nEimEi and nEi are the molecular weight and the enzyme cost of an Ei molecule, respectively.See Appendix 1—table 1.
rcarbon, rproteinrcarbon and rprotein are the mass fractions of all carbon and protein within a cell, respectively.rprotein=0.55; rcarbon=0.48
(Neidhardt et al., 1990).
κiSubstrate quality of a metabolite in a biochemical reaction; see Equation S12 and S20.Calculated from the values of kicat, mEi, m0, rprotein, rcarbon.
κASubstrate quality of a Group A carbon source; see Equation S27.Calculated from the values of kAcat, mEA, m0, rprotein, rcarbon, KA and the concentration of the Group A carbon source [A].
ϕiThe proteomic mass fraction of enzyme Ei: ϕiMEi/Mprotein; an intensive variable.See Equation S9.
ηiThe fraction of stoichiometric flux drawn from a precursor pool; see Equations S13, S14 and S18.ηa1=15%; ηa2=30%; ηb=35%; ηc=9%; ηd=11% (calculated from the values of ri and NEPicarbon).
ϕr, ϕf, ϕBMϕf, ϕf, ϕBM are the proteomic mass fraction of enzymes dedicated to fermentation, respiration, and biomass generation, respectively.NA
κtA parameter determined by the translation rate, defined as κtkTmAA/(ςmR).κt=1/610 (s–1) (calculated from the values of kT, ς and mR).
JBMThe carbon flux of biomass production; an extensive variable.See Equation S10.
JEThe energy demand for cell growth, expressed as the stoichiometric energy flux in ATP; an extensive variable.See Equation S25.
JE(N)The normalized flux of energy demand in ATP; an intensive variable.JE(N)JEm0/Mcarbon.
rE, ηErE and ηE are energy coefficients. rE is the slope of JE versus JBM; ηE=rE[iri/NEPicarbon].See Appendix 9.2.
βiThe stoichiometric coefficient of ATPs in biochemical reactions shown in Figures 1B and 3E (for E. coli) or Appendix 1—figure 5E and F (for yeast and mammalian cells).β1=4, β2=3, β3=2, β4=6, β6=1, βa1=4, β7=1, β8=2, β9=6 (E. coli); β1=5, β2=1, β3=5, β4=7.5, β6=2.5, βa1=5 (eukaryotic cells) (Neidhardt et al., 1990; Sauer et al., 2004).
βr(A), βf(A)βr(A) and βf(A) are the stoichiometric coefficients of ATP production per glucose in respiration and fermentation, respectively.βr(A)=26, βf(A)=12 (E. coli); βr(A)=32, βf(A)=2 (eukaryotic cells) (Neidhardt et al., 1990).
Jr(E), Jf(E)Jr(E) and Jf(E) are normalized energy fluxes of respiration and fermentation, intensive variables.Jr(E)βr(A)2Jr(N);Jf(E)βf(A)2Jf(N).
εr, εf
εr(dt), εf(dt)
εr (or εr(dt)) and εf (or εf(dt)) are the proteome efficiencies for energy biogenesis in the respiration and fermentation pathways: εrJr(E)/ϕr and εfJf(E)/ϕf.Calculated from the values of κA, κi, βr(A) and βf(A) with Equations S132 and S161.
φφ is an energy demand coefficient, defined in Equation S33 and mainly determined by ηE.Calculated from the values of ηE, βi, ηi with Equation S33. See Appendix 9.2.
ψ, ψdtψ-1 (or ψdt1) is the proteome efficiency for biomass generation in the biomass pathway, with ψ1/λ/ϕBM.Calculated from the values of ηi, κA, κi, Ω, κt with Equations S133 and S162.
κrA, κf(A)κrA and κfA are parameters defined as κr(A)[1κ1+2κ2+2κ3+2κ4]1 and
κf(A)[1κ1+2κ2+2κ6]1.
Calculated from the values of κi.
ΩΩ is a composite parameter defined as Ω1/κt+ia1,a2,b,c,dηi/κi.See Appendix 9.2.
κglucose(ST),κlactose(ST)The substrate quality of glucose or lactose at saturated concentration.Calculated using Equation S27 and the approximation used in Equation S20.
ΔΔ is a function of κA defined as Δ(κA)εf(κA)/εr(κA) .Δεf/εr.
κA(C)The critical value of κA which satisfy Δ(κA)=1 and thus εf(κA)=εr(κA);
See Equation S42 (for E. coli) and S176 (for yeast and mammalian cells).
Calculated from the values of βi and κi with Equation S42.
λCThe critical growth rate at the transition point: λCλ(κA(C)); See Equations S43 and S177.Calculated from the values of ϕmax, φ, βi, κi, κA(C), Ω, ηi with Equations S43, S32 and S162.
θThe Heaviside step function.NA
Jacetate,JCO2,rJacetate and JCO2,r are the stoichiometric fluxes of acetate from the fermentation pathway and CO2 from the respiration pathway; extensive variables.Jacetate=Jf; JCO2,r=3Jr. See Appendix 9.1 and Equations S158.
Jacetate(M),JCO2,r(M)Jacetate(M) and JCO2,r(M) are the fluxes of Jacetate and JCO2,r (per biomass) in the unit of mM/OD600/h, which are measurable in experiment. Intensive variables.Jacetate(M)2Jf(N); JCO2,r(M)6Jr(N). See Appendix 9.1 and Equation S160.
κAmaxThe maximum value of κA available across different Group A carbon sources.Approximated by the max κA across Group A carbon sources, calculated with Equation S27 and the approximation used in Equation S20.
λmaxThe population cell growth rate for the maximum value of κA: λmax=λκAmax.Calculated from the maximum of Equation S36 with the values of βi, κi, κAmax, φ, Ω, κt, and Equations S32, S132, Equation S161 and S162.
N(μ,σ2)A Gaussian distribution with a mean of μ and a standard deviation of σ.The probability density function is f(x)=1σ2πe12(xμσ)2.
μλC, σλCμλC and σλC are the mean and standard deviation of λC, respectively.μλC is approximated by the deterministic value of λC; see Appendix 3.3 for σλC settings. See Appendix 9.2 for the values.
erfThe error function in mathematics.erf(x)=2π0xexp(t2)dt
ϕZThe proteomic mass fraction of useless proteins encoded by the LacZ gene.See Appendix 4.1.
wAn energy dissipation coefficient.See Appendix 4.2.
w0The maintenance energy coefficient.w0=0 or 2.5 (h–1) as specified in Figures 34, Appendix 1—figures 2 and 3. See Appendices 4.3 and 9.2.
ιι is the inhibition coefficient such that (1+ι)1 represents the translation efficiency.See Appendices 4.3 and 9.2
ιw0=0(2μmCm), ιw0=0(4μmCm), ιw0=0(8μmCm), ιw0=2.5(2μmCm), ιw0=2.5(4μmCm), ιw0=2.5(8μmCm)The values for ι in the cases with 2 μm , 4 μm, or 8 μm of chloramphenicol and the maintenance energy coefficient w0 chosen as 0 or 2.5 (h–1).ιw0=0(2μmCm)=1.15;ιw0=0(4μmCm)=2.33; ιw0=0(8μmCm)=6.25; ιw0=2.5(2μmCm)=1.05; ιw0=2.5(4μmCm)=2.00; ιw0=2.5(8μmCm)=5.40. See Appendix 9.2.
κpyThe substrate quality of pyruvate; see Equation S89.Calculated from the values of kpycat, mEpy, m0, rprotein, rcarbon, Kpy and the external concentration of pyruvate [py].
βr(py), βf(py)βr(py) and βf(py) are the stoichiometric coefficients of ATP production per pyruvate in respiration and fermentation, respectively.βr(py)=10; βf(py)=3.
(Neidhardt et al., 1990).
Jr(E,py), Jf(E,py)Jr(E,py) and Jf(E,py) are the normalized energy fluxes of respiration and fermentation for pyruvate utilization; intensive variables.The corresponding variables of Jr(E) and Jf(E) in the case of pyruvate utilization.
εr(py),εf(py)εr(py) and εf(py) are the proteome efficiencies for energy biogenesis using pyruvate in the respiration and fermentation pathways.The corresponding variables of εr and εf in the case of pyruvate utilization.
ΩGgΩ`Gg is a composite parameter defined as ΩGg(ηb+ηc)/κ8+ηa1/κ9.See Appendix 9.2.
ψpy, φpy, κpy(ST)κpy(C), λmax(py)ψpy, φpy, κpy(ST), κpy(C) and λmax(py) are the corresponding variables/parameters of ψ, φ, κAmax, κAC and λmax in the case of pyruvate utilization.See Appendices 5.1 and 9.2.
λC(py), μλC(py), σλC(py)λC(py), μλC(py) and σλC(py) are the corresponding variables/parameters of λC, μλC and σλC in the case of pyruvate utilization.See Appendices 5.1 and 9.2.
NPicarbonThe number of carbon atoms in a molecule of Pool i.The value of NPicarbon is approximated by NEPicarbon (Equation S107).
κi(21AA)The substrate quality of the external supplied amino acids identical to those in Pool i.See Appendices 5.2 and 9.2.
Ω21AAΩ21AA is a composite parameter defined as Ω21AA1/κt+ηa1/κa1+ia2,b,c,dηi/κi(21AA).See Appendices 5.2 and 9.2.
ψ21AA, φ21AA, λmax(21AA), λC(21AA), μλC(21AA), σλC(21AA)ψ21AA, φ21AA, λmax(21AA), λC(21AA), μλC(21AA) and σλC(21AA) are the corresponding variables/parameters of ψ, φ, λmax, λC, μλC and σλC in the case of a Group A carbon source is mixed with 21 types of amino acids at saturated concentrations.See Appendices 5.2 and 9.2.
Ω7AA, φ7AA, μλC(7AA),σλC(7AA)Ω7AA, φ7AA, μλC(7AA) and σλC(7AA) are the corresponding parameters of Ω, φ, μλC and σλC in the case of a Group A carbon source is mixed with 7 types of amino acids.See Appendices 5.2 and 9.2.
Jin(N), ϑJin(N) is the normalized stoichiometric influx of a Group A carbon source (Equation S136). ϑ is a parameter defined as ϑ=ηa1+ηc+(ηa2+ηb+ηd)/2 for the model shown in Figure 1B.See Appendix 7.3
χext, χint, χtotχext, χint and χtot are the level of extrinsic noise, intrinsic noise and total noise in a system.See Appendix 8.1
μkicat, σkicat, μ1/kicat, σ1/kicat, μ1/kicat, σ1/kicatμkicat and σkicat are the mean and standard deviation of kicat. μ1/kicat (or μ1/kicat) and σ1/kicat (or σ1/kicat) are the mean and standard deviation of 1/kicat. See Appendix 8.1.μkicat is approximated by the deterministic value of kicat. The CV of kicat is set to 25%. μ1/kicat≈1/μkicat; σ1/kicat/μ1/kicatσkicat/μkicat.
IG(x;μ,ζ)The inverse Gaussian (IG) distribution: variable x>0 with parameters μ and ζ. See Equation S142.The probability density function is ζ2πx3exp-ζx-μ22μ2x.
IOGx;μ,ζThe positive inverse of Gaussian (IOG) distribution: variable x>0 with parameters μ and ζ. See Equation S140 and Appendix 8.1.The probability density function is ζ2πx4exp-ζx-μ22μ2x2.
ζ1/kicat, ζ1/kicatDistributional parameters of 1/kicat corresponding to ζ in an IG or IOG distribution.See Appendix 8.1
GkThe characteristic function of IG distribution. See Equation S147.Gk=-eikxIGx;μ,ζdx
Xi, αi, Θ, TΘ,Γi(t)Xi, αi, Θ and Γi(t) are variables and parameters used to calculate the first passage time TΘ of a stochastic process that mimics the duration of an enzyme to finishing a catalytic job.See Appendix 8.1.
γi, Ξ, μΞ, σΞγi is a real number; Ξ is a variable defined as Ξi=1nγi/kicat; μΞ and σΞ are the mean and standard deviation of Ξ.See Equation S153 and Appendix 8.1.
μκi, σκi, μ1/κi, σ1/κiμκi and σκi are the mean and standard deviation of κi; μ1/κi and σ1/κi are the mean and standard deviation of 1/κi.See Equation S154 and Appendices 8.1 and 9.2.
λr, λf, μλr, σλr, μλf, σλf, ρrfλr and λf are the growth rates when cells choose respiration or fermentation; μλr, μλf and σλr, σλf are the means and standard deviations of λr and λf; ρrf is the correlation of λr and λf.See Equation S36 and Appendices 8.1 and 9.2.
λsuccinate(21AA), λacetate,μλsuccinate(21AA), μλacetate,σλsuccinate(21AA), σλacetateλsuccinate(21AA) and λacetate are the growth rates for succinate mixed with 21AA or acetate as the sole carbon source; μλsuccinate(21AA), μλacetate and σλsuccinate(21AA), σλacetate are the means and standard deviations of λsuccinate(21AA) and λacetate.See Appendix 9.2.
ϕMT, κMTϕMT and κMT are the proteomic mass fraction of the enzymes and the effective substrate quality of related metabolites in the mitochondria for yeast and mammalian cells, respectively.NA
PrfThe proportion of ATP generated from fermentation: PrfJf(E)Jf(E)+Jr(E).See Equations S180, S189 and Appendix 10.
Δ¯The proteome efficiency difference between respiration and fermentation: Δ¯1/εr1/εf.See Equations S181, S187 and Appendix 10.
μεr, μεf, μ1/εr, μ1/εfμεr, μεf, μ1/εr and μ1/εf are the mean values of εr , εf, 1/εr and 1/εf, respectively.See Equations S182-S184 and Appendix 10.
σεr, σεf, σ1/εr, σ1/εfσεr, σεf, σ1/εr, and σ1/εf are the standard deviations of εr , εf, 1/εr and 1/εf, respectively.See Equations S182, S185 and Appendix 10.
χεr, χεf, χ1/εr, χ1/εfχεr, χεf, χ1/εr, and χ1/εf are the coefficients of variation of εr , εf, 1/εr and 1/εf, respectively.See Equations S185-S186 and Appendix 10.
μΔ¯, σΔ¯μΔ¯ and σΔ¯ are the mean and standard deviation of Δ¯, respectively.See Equations S187-S188 and Appendix 10.
εr,εfεr and εf are the population-averaged values of εr and εf, respectively.Measurable from experiments.
See Equations S183-S184 and Appendix 10.
  1. *

    Parameter settings for yeast and mammalian cells are specifically labeled as ‘eukaryotic cells.’

  2. ‘NA’ represents ‘Not applicable.’

  3. Extensive variables scale with the size of the cell population.

  4. §

    Intensive variables are scale-invariant with respect to the cell population.

Appendix 1—figure 1
Central metabolic network and carbon utilization pathways of E. coli.

(A) Energy biogenesis details in the central metabolic network. In E. coli, NADPH and NADH are interconvertible (Sauer et al., 2004), and all energy carriers can be converted to ATP through ADP phosphorylation. The conversion factors are: NADH = 2 ATP, NADPH = 2 ATP, FADH2=1 ATP (Neidhardt et al., 1990). (B) Relevant genes encoding enzymes in the central metabolic network of E. coli. (C–E) Three independent fates of glucose metabolism in E. coli. (C) For energy biogenesis through fermentation, a molecule of glucose generates 12 ATPs. (D) For energy biogenesis via respiration, a molecule of glucose generates 26 ATPs. (E) For biomass synthesis, glucose is converted into precursors of biomass. Note that biomass synthesis is accompanied by ATP production (see Appendix 3.1).

Appendix 1—figure 2
Model and results for experimental comparison of E. coli.

(A–C) Model analysis for carbon utilization in mixtures with amino acids. (A) Coarse-grained model for the case of a Group A carbon source mixed with extracellular amino acids. (B) Model predictions (Equations S157, S164-S165) and single-cell reference experimental results (Wallden et al., 2016) showing growth rate distributions for E. coli in three culturing conditions. (C) Comparison of the growth rate-fermentation flux relation for E. coli in Group A carbon sources between minimal media and enriched media (those with 7AA). (D–E) Influence of translation inhibition on overflow metabolism in E. coli. (D) A 3D plot illustrating the relations among fermentation flux, growth rate, and translation efficiency (Equations S79 and S160). (E) Growth rate dependence of acetate excretion rate as κA varies, with each fixed dose of Cm. Translation efficiency is tuned by the dose of Cm, and the maintenance energy coefficient is set to 0 (i.e. w0=0). (F) Coarse-grained model for Group A carbon source utilization, which includes more details to compare with experiments. (G) Comparison of in vivo and in vitro catalytic rates for enzymes of E. coli within glycolysis and the TCA cycle (see Appendix 1—table 1 for details). (H) The proteome efficiencies for energy biogenesis in the respiration and fermentation pathways vary with growth rate as functions of the substrate quality of pyruvate (Equations S93 and S96)

Appendix 1—figure 3
Relative protein expression of central metabolic enzymes in E. coli under various types of perturbations.

(A–D) Relative protein expression under κA perturbation. (A) Experimental data (Hui et al., 2015) for the catalytic enzymes at each step of glycolysis. (B) Experimental data (Hui et al., 2015) for the catalytic enzymes at each step of the TCA cycle. (C) Model predictions (Equation S118, with w0=0) and experimental data (Hui et al., 2015) for representative glycolytic genes. (D) Model predictions (Equation S118, with w0=0) and experimental data (Hui et al., 2015) for representative genes from the TCA cycle. (E–J) Relative protein expression under ϕZ perturbation. (E, F, I) Model predictions and experimental data (Basan et al., 2015) for representative glycolytic genes. (G, H, J) Model predictions and experimental data (Basan et al., 2015) for representative genes from the TCA cycle. (E–H) Results of ϕZ perturbation with w0=0 (Equation S120). (I–J) Results of ϕZ perturbation with w0=2.5(h1) (Equation S121). (K–N) Relative protein expression upon energy dissipation. (K–L) Model fits (Equations S127 and S123) and experimental data (Basan et al., 2015) for representative glycolytic genes. (M–N) Model fits (Equations S127 and S123) and experimental data (Basan et al., 2015) for representative genes from the TCA cycle.

Appendix 1—figure 4
Asymptotic distributions of inverse Gaussian distribution and the inverse of Gaussian distribution.

(A) Comparison between the inverse of Gaussian distribution and the corresponding Gaussian distribution for various values of the coefficient of variation (CV) (Equations S140 and S145). (B) Comparison between the inverse Gaussian distribution and the corresponding Gaussian distribution for various values of CV (Equations S142 and S146). Both the inverse Gaussian distribution and the inverse of Gaussian distribution converge to Gaussian distributions when CV is small.

Appendix 1—figure 5
Carbon utilization in yeast and mammalian cells.

(A–D) Three independent fates of glucose metabolism in yeast and mammalian cells. (A–B) For energy biogenesis through fermentation, one molecule of glucose generates 2 ATPs. (C) For energy biogenesis through respiration, one molecule of glucose generates 32 ATPs. (D) For biomass synthesis, glucose is converted into biomass precursors, with ATP produced as a byproduct. In yeast and mammalian cells, the energy stored in NADH and FADH2 converts ADP into ATP in the mitochondria, with higher conversion factors than in E. coli: NADH = 2.5 ATP, FADH2=1.5 ATP (Nelson and Cox, 2008). (E) Coarse-grained model for Group A carbon source utilization in yeast. (F) Coarse-grained model for Group A carbon source utilization in mammalian cells.

Appendix 2

Model framework

2.1 Proteome partition

Here, we adopt the proteome partition framework similar to that introduced by 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 Mi(P) (i=Q,R,C) and mass fraction ϕi, where ϕQ is a constant, and we define ϕmax1ϕQ. In the exponential growth phase, the ribosome allocation for protein synthesis of each class is fi, with fQ+fR+fC=1.

To analyze cell growth optimization, we first consider the homogeneous case where all cells share identical biochemical parameters, simplifying the mass accumulation of the cell population into a ‘big cell.’ This simplification does not affect the value of growth rate λ. For bacteria, the protein turnover is negligible, so the mass accumulation of each class follows:

(S1) dMi(P)/dt=fikTNRmAA(i=Q,R,C),

where mAA stands for the average molecular weight of amino acids, kT is the translation rate, NR=Mrp(P)/mR is the number of ribosomes, mR is the protein mass of a single ribosome, and Mrp(P) is the total protein mass of ribosomes, with MR(P)/Mrp(P)=ς1.67 (Neidhardt, 1996; Scott et al., 2010). For a specific stable nutrient environment, fR and kT are temporal invariants. Thus,

(S2) Mi(P)(t)=Mi(P)(0)+fi/fRMR(P)(0)[exp(λt)1](i=Q,R,C),

where λ=fRkTmAA/(ςmR), and the total protein mass of the cell population MproteiniQ,R,CMi(P) follows:

(S3) Mprotein(t)=Mprotein(0)+MR(0)[exp(λt)1]/fR

Over a long period in the exponential growth phase (i.e. t+), we have ϕi=fi (i=Q,R,C), and

(S4) λ=ϕRκt,

where κt=kTmAA/(ςmR).

2.2 Precursor pools

Based on the entry points of the metabolic network, we classify the precursors of biomass components into five pools (Figure 1A and 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). These five pools draw approximately ra1=24%, ra2=24%, rb=28%, rc=12%, and rd=12% of the carbon flux (Nelson and Cox, 2008; Wang et al., 2019). There are overlapping components between Pools a1 and a2 due to the joint synthesis of some precursors. Therefore, we use Pool a to represent both Pools a1 and a2 in the descriptions.

2.3 Stoichiometric flux

We consider the following biochemical reaction between substrate Si and enzyme Ei:

(S5) Ei+SidiaiEiSikicatEi+biSi+1+ciCO2,

where ai, di and kicat are the reaction parameters, Si+1 is the product, bi and ci are the stoichiometric coefficients. For most of the reactions in the central metabolism, bi=1 and ci=0. The reaction rate follows Michaelis–Menten kinetics (Nelson and Cox, 2008):

(S6) vi=kicat[Ei][Si][Si]+Ki,

where Ki(di+kicat)/ai, [Ei] and [Si] are the Michaelis constant, and the concentrations of enzyme Ei and substrate Si, respectively. For this reaction (Equation S5), d[Si+1]/dt=bivi and d[Si]/dt=vi. In the cell population (the ‘big cell’), suppose that the cell volume is Vcell, then the stoichiometric flux of the reaction is:

(S7) JiVcellvi.

The copy number of enzyme Ei is NEi=Vcell[Ei] with a total weight of MEi=NEimEi, where mEi is the molecular weight of Ei. By defining the enzyme cost of an Ei molecule as nEimEi/m0, where m0 is a unit mass, then the cost of all Ei molecules is ΦiNEinEi (Wang et al., 2019). By further defining ξikicatnEi[Si][Si]+Ki, then:

(S8) Ji=Φiξi.

The mass fraction of enzyme Ei in the proteome is ϕi=MEi/Mprotein, and thus:

(S9) ϕi=Φim0Mprotein.

2.4 Carbon flux and cell growth rate

To clarify the relation between the stoichiometric flux Ji 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 Mcarbon=Mproteinrcarbon/rprotein , where rcarbon and rprotein 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:

(S10) JBM=1mcarbondMcarbondt=λMcarbonmcarbon,

where mcarbon 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 Jbcarbon=JbNpycarbon, where Jb is the stoichiometric flux from pyruvate to Pool b (Figure 1A and B) and Npycarbon stands for the carbon number of a pyruvate molecule. Combining with Equation S10 and noting that Jbcarbon=rbJBM, we get JbNpycarbonmcarbon=rbλMcarbon. Similarly, for each precursor pool, we have:

(S11) JiNEPicarbonmcarbon=riλMcarbon(i=a1,a2,b,c,d),

where the subscript ‘EPi’ represents the entry point of Pool i, and NEPi is 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:

(S12) κiξirproteinrcarbon=rproteinrcarbonkicatnEi[Si][Si]+Ki,

and for each precursor pool, we define:

(S13) ηirim0/(NEPicarbonmcarbon)(i=a1,a2,b,c,d).

Combining Equations S8, S9 and S11, we have

(S14) ϕκi=ηiλ(i=a1,a2,b,c,d).

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

(S15) Ji(N)ϕiκi,

where the superscript ‘(N)’ stands for normalized. Combined with Equations S8, S9 and S12, we have:

(S16) Ji(N)Jim0Mcarbon.

Since ia1,a2,b,c,dri=1, by setting

(S17) m0=[iri/NEPicarbon]1mcarbon,

we then obtain:

(S18) ηi=riNEPicarbon[ja1,a2,b,c,drjNEPjcarbon]1(i=a1,a2,b,c,d),

and we have ia1,a2,b,c,dηi=1, and

(S19) ia1,a2,b,c,dϕiκi=λ.

2.5 Intermediate nodes

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

(S20) κirproteinrcarbonkicatnEi

Real cases could be more complicated due to other forms of 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, Ki is lower than the substrate concentration (i.e. [Si]>Ki), which implies κirproteinrcarbonkicatnEi.

Appendix 3

Model and analysis

3.1 Coarse-grained model

In the coarse-grained model shown in Figure 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 Equation S5 with bi=1 except for M1→2M2 and M3 +M5→M4. Basically, there are three independent possible fates for a Group A carbon source (e.g. glucose; see Appendix 1—figure 1C-E; Chen and Nielsen, 2019): energy biogenesis through fermentation; energy biogenesis via respiration (Appendix 1—figure 1C and D), or conversion into biomass components accompanied by energy biogenesis in the biomass pathway. Each fate involves a distinct fraction of the proteome, with no overlap between them (Appendix 1—figure 1).

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

(S21) {ΦAξA=Φ1ξ1+Φa1ξa1,2Φ1ξ1=Φ2ξ2+Φ5ξ5+Φa2ξa2,Φ2ξ2=Φ3ξ3+Φ6ξ6+Φbξb,Φ5ξ5+Φ4ξ4=Φ3ξ3+Φdξd,Φ3ξ3=Φ4ξ4+Φcξc.

Obviously, the stoichiometric fluxes of respiration Jr and fermentation Jf (Appendix 1—figure 1C and D) are:

(S22) {JrJ4=Φ4ξ4,JfJ6=Φ6ξ6.

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

(S23) ΦiξiNEPicarbon=riJBM(i=a1,a2,b,c,d).

In terms of energy biogenesis 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), FADH2→1ATP (Neidhardt et al., 1990). Then, we have

(S24) β1Φ1ξ1+β2Φ2ξ2+β3Φ3ξ3+β4Φ4ξ4+β6Φ6ξ6+βa1Φa1ξa1=JE,

where JE represents the energy demand for cell proliferation, expressed as the stoichiometric energy flux in ATP. βi is the stoichiometric coefficient with β1=4, β2=3, β3=2, β4=6, β6=1, and βa1=4 for E. coli (Neidhardt et al., 1990; Sauer et al., 2004). For bacteria, the energy demand is generally proportional to the carbon flux infused into biomass production, as the proportion of maintenance energy is roughly negligible (Locasale and Cantley, 2010). Thus,

(S25) JE=rEJBM,

where rE is the ratio and also a constant.

By applying the substitutions specified in Equations S9, S12, S14-S18, combined with Equations S4, S10, S21-S25, and the constraint of proteome resource allocation ϕR+ϕC=ϕmax, we have:

(S26) {ϕAκA=ϕ1κ1+ϕa1κa1,2ϕ1κ1=ϕ2κ2+ϕ5κ5+ϕa2κa2,ϕ2κ2=ϕ3κ3+ϕ6κ6+ϕbκb,ϕ5κ5+ϕ4κ4=ϕ3κ3+ϕdκd,ϕ3κ3=ϕ4κ4+ϕcκc,ϕa1κa1=ηa1λ,ϕa2κa2=ηa2λ,ϕbκb=ηbλ,ϕcκc=ηcλ,ϕdκd=ηdλ,β1ϕ1κ1+β2ϕ2κ2+β3ϕ3κ3+β4ϕ4κ4+β6ϕ6κ6+βa1ϕa1κa1=JE(N),JE(N)=ηEλ,λ=ϕRκt,Jr(N)=ϕ4κ4,Jf(N)=ϕ6κ6,ϕR+ϕA+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ϕ6+ϕa1+ϕa2+ϕb+ϕc+ϕd=ϕmax,

where JE(N) and ηE are defined as JE(N)JEm0Mcarbon and ηErE[iri/NEPicarbon]1, respectively. Here, for each intermediate node, κi follows Equation 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:

(S27) κArproteinrcarbonkAcatmEA[A][A]+KAm0,

which is determined externally by the culture condition. From Equation S26, all ϕi and ϕR can be expressed in terms of Jr(N), Jf(N), and λ:

(S28) {ϕA=[Jr(N)+Jf(N)+(2ηa1+ηa2+ηb+2ηc+ηd)λ]/(2κA),ϕ1=[Jr(N)+Jf(N)+(ηa2+ηb+2ηc+ηd)λ]/(2κ1),ϕ2=[Jr(N)+Jf(N)+(ηb+ηc)λ]/κ2,ϕ3=(Jr(N)+ηcλ)/κ3,ϕ4=Jr(N)/κ4,ϕ5=(ηc+ηd)λ/κ5,ϕ6=Jf(N)/κ6,ϕR=λ/κt,ϕi=ηiλ/κi(i=a1,a2,b,c,d).

In Equation S28, for each ϕi or ϕR, the Jr(N)- and Jf(N)-related proteome fraction terms belong to the fractions of the proteome dedicated to respiration (denoted as ϕr) and fermentation (denoted as ϕf), respectively. The λ-related proteome fraction terms belong to those involved in the biomass synthesis pathway (denoted as ϕBM). Thus, ϕr=Jr(N)[1/(2κA) +1/(2κ1) +1/κ2 +1/κ3 +1/κ4], ϕf=Jf(N) [1/(2κA) +1/(2κ1) +1/κ2 +1/κ6], and ϕBM=λ(1κt +1+ηa1+ηc2κA +1ηa1+ηc2κ1 +ηb+ηcκ2 +ηcκ3 +ηc+ηdκ5 +ia1,a2,b,c,dηiκi). By substituting Equation S28 into Equation S26, we have:

(S29) {Jr(E)+Jf(E)=φλ,Jr(E)εr+Jf(E)εf=ϕmaxψλ.

Here, Jr(E) and Jf(E) stand for the normalized energy fluxes of respiration and fermentation, with

(S30) {Jr(E)=βr(A)Jr(N)/2,Jf(E)=βf(A)Jf(N)/2,

where βr(A)=β1+2(β2+β3+β4) and βf(A)=β1+2(β2+β6), with βr(A)=26 and βf(A)=12 for E. coli. εr and εf represent the proteome efficiencies for energy biogenesis in the respiration and fermentation pathways (Appendix 1—figure 1C-D), defined as εrJr(E)/ϕr and εfJf(E)/ϕf; that is, the normalized energy fluxes expressed in ATP generated per proteomic mass fraction dedicated to respiration and fermentation, respectively. Hence,

(S31) {lεr=βr(A)1/κA+1/κ1+2/κ2+2/κ3+2/κ4,εf=βf(A)1/κA+1/κ1+2/κ2+2/κ6.

ψ1 is the proteome efficiency for biomass generation in the biomass synthesis pathway (Appendix 1—figure 1E), defined as ψ1λ/ϕBM=ia1,a2,b,c,dJi(N)/ϕBM (see Equations S15 and S19); that is, the normalized flux (which differs from the normalized energy flux used to define εr and εf) generated per proteomic mass fraction dedicated to biomass synthesis. Hence

(S32) ψ=1κt+1+ηa1+ηc2κA+ηa2+ηb+2ηc+ηd2κ1+ηb+ηcκ2+ηcκ3+ηc+ηdκ5+ia1,a2,b,c,dηiκi.

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

(S33) φηEβ1(ηa2+ηb+2ηc+ηd)/2β2(ηb+ηc)β3ηcβa1ηa1,

and φλ stands for the normalized flux of energy demand other than the accompanying energy biogenesis from the biomass synthesis pathway.

3.2 The reason for overflow metabolism

Microbes optimize their growth rate to survive through the evolutionary process (Vander Heiden et al., 2009). The optimal growth principle also roughly holds for tumor cells, which proliferate while ignoring growth restriction signals and evading immune destruction by the host (Vander Heiden et al., 2009). First, we consider the optimal growth strategy for a single cell. The coarse-grained model for bacteria is summarized in Equation S26 and further simplified in Equation S29. Here, εr, εf and ψ are functions of κA (see Equations S31, S32), so we also denote them as εr(κA), εf(κA), ψ(κA). Evidently, the fluxes of both respiration and fermentation take non-negative values, i.e., Jr(E),Jf(E)0, and all the coefficients are positive: εr(κA),εf(κA),ψ(κA),φ>0.

Thus, if εr>εf, then (ψ+φ/εr)λ=ϕmaxJf(E)(1/εf1/εr)ϕmax. Obviously, the solution for optimal growth is:

(S34) {Jf(E)=0,Jr(E)=φλ.εr>εf.

Similarly, if εf>εr, then the optimal growth solution is:

(S35) {Jf(E)=φλ,Jr(E)=0.εr<εf.

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

(S36) λ={λr=ϕmaxφ/εr(κA)+ψ(κA)εr(κA)>εf(κA),λf=ϕmaxφ/εf(κA)+ψ(κA)εr(κA)<εf(κA).

So, why do microbes use the seemingly wasteful fermentation pathway when the growth rate is large under aerobic conditions? Prevalent explanations (Basan et al., 2015; Chen and Nielsen, 2019) suggest that it originates from that the proteome efficiency in fermentation is consistently higher than in respiration (i.e. εf>εr). If this is the case, why do microbes still use the normal respiration pathway when the growth rate is small? The answer lies in the fact that both εr(κA) and εf(κA) are not constants, but are dependent on nutrient conditions. In Equation S31, when κA is small, consider the extreme case of κA0, and then

(S37) {εr(κA0)βr(A)κA,εf(κA0)βf(A)κA.

Since βr(A)βf(A), clearly,

(S38) εr(κA0)>εf(κA0).

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

(S39) βr(A)βf(A)κA<βf(A)(1κ1+2κ2+2κ3+2κ4)βr(A)(1κ1+2κ2+2κ6),

then Δ(κA)εf(κA)/εr(κA) is a monotonically increasing function of κA. Thus,

(S40) εr(κAmax)<εf(κAmax),

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

In practice, experimental studies using E. coli (Basan et al., 2015) have demonstrated that proteome efficiency in fermentation is higher than in respiration when the Group A carbon source is lactose at a saturated concentration, i.e., εr(κlactose(ST))<εf(κlactose(ST)). Here, κlactose(ST) 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, κAmax>κlactose(ST), and hence, Equation S40 holds for E. coli. From a theoretical perspective, we can verify Equation S39 and consequently Equation S40 using Equation S20, combined with the in vivo/in vitro biochemical parameters obtained from experimental data (see Appendix 1—table 1; Appendix 1—table 2). For example, it is straightforward to confirm that εr(κglucose(ST))<εf(κglucose(ST)) using this method (see Appendix 9.2), further supporting the validity of Equations S39-S40 (see also Appendix 10).

Now that Equations S38-S40 are all valid, a critical value of κA, denoted as κA(C), exists, satisfying Δ(κA(C))=1. Thus,

(S41) {εf(κA)>εr(κA),κA>κA(C);εf(κA)=εr(κA),κA=κA(C);εf(κA)<εr(κA),κA<κA(C).

Combined with Equation S31, we have:

(S42) κA(C)=βr(A)βf(A)βf(A)(1/κ1+2/κ2+2/κ3+2/κ4)βr(A)(1/κ1+2/κ2+2/κ6).

By substituting Equation S42 into Equations S31, S32 and S36, we obtain the expressions for εr(κA(C)), εf(κA(C)) and the critical growth rate at the transition point (i.e. λCλ(κA(C))):

(S43) {εr(κA(C))=εf(κA(C))=βr(A)βf(A)2(1/κ3+1/κ41/κ6)=β3+β4β61/κ3+1/κ41/κ6,λC=ϕmaxφ/εr/f(κA(C))+ψ(κA(C)),

where εr/f represents either εr or εf. In Figure 1E, we show the dependencies of εr(κA), εf(κA), and λ(κA) on κA in a three-dimensional form, as κA changes.

3.3 The relation between respiration/fermentation flux and growth rate

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

In the homogeneous case, where all microbes share identical biochemical parameters, as λ(κA) increases with κA, Jf(E) appear abruptly and Jr(E) vanish simultaneously as κA exceeds κA(C) (Figure 1E; see also Equations S34-S35, S41). Combining Equations S34-S36 and S43, we obtain:

(S44) {Jf(E)=φλθ(λλC),Jr(E)=φλ[1θ(λλC)].

where ‘θ’ stands for the Heaviside step function. Defining λmax=λ(κAmax), and then, [0,λmax] is the relevant range of the x-axis. In fact, the digital responses in Equation S44 are consistent with the numerical simulation results of 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 kicat can be greatly influenced by the concentrations of potassium and phosphate (García-Contreras et al., 2012), which vary from cell to cell. Consequently, there is a distribution of values for kicat among cell populations, commonly referred to as extrinsic noise (Elowitz et al., 2002). For convenience, we assume that each kicat (and thus κi) follows a Gaussian distribution with a coefficient of variation (CV) of 25%. Therefore, the distributions of proteome efficiencies that determine the choice between respiration and fermentation, εr and εf, and the critical growth rate for the transition, λC, can be approximated by Gaussian distributions for a cell population (see Appendix 8.1 for details). Specifically, λC follows:

(S45) λCN(μλC,σλC2),

where μλC and σλC represent the mean and standard deviation of λC, with the CV σλC/μλC calculated to be 12% (see Appendix 9.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:

(S46) {Jf(E)(λ)=12φλ[erf(λμλC2σλC)+1],Jr(E)(λ)=12φλ[1erf(λμλC2σλC)],

where ‘erf’ represents the error function. In practice, given a culturing medium, there is also a probability distribution for the growth rate (Appendix 1—figure 2B; see also Equation S157). For approximation, in plotting the growth rate-respiration/fermentation flux relations, we use the deterministic (noise-free) value of the growth rate as a proxy. To compare with experiments, we essentially compare the normalized fluxes, Jr(N) and Jf(N) (see Appendix 9.1 for details). Combining Equations S30 and S46, we obtain:

(S47) {Jf(N)(λ)=φβf(A)λ[erf(λμλC2σλC)+1],Jr(N)(λ)=φβr(A)λ[1erf(λμλC2σλC)].

In Figure 1C–D, we see that Equation S47 quantitatively illustrates the experimental data (Basan et al., 2015), where the model parameters were obtained using biochemical data for the catalytic enzymes (see Appendix 1—table 1 for details).

3.4 Dependence of the model on optimization principles

In the derivation of the growth rate dependence of respiration/fermentation flux described above (Equation S44 for the single-cell level and Equation S47 for the population-averaged level), we applied the principles of optimal growth, incorporating both efficient protein allocation to enzymes and ribosomes (through ribosomal proteins). However, recent experimental studies show that the inactive portion of ribosomes (i.e. ribosomes not bound to mRNAs) may vary with culturing conditions (Dai et al., 2017; Li et al., 2018) and between individual cells within the same culture (Pavlou et al., 2025), despite an overall trend toward growth optimization. These findings (Dai et al., 2017; Li et al., 2018; Pavlou et al., 2025) suggest that ribosome allocation may be suboptimal under many culturing conditions, likely as a result of cells preparing for potential environmental changes (Li et al., 2018). Nevertheless, since our model’s predictions regarding the binary choice between respiration and fermentation rely solely on comparing proteome efficiency between these two pathways, which involves only efficient protein allocation to enzymes, and because the active portion of ribosomes and the translation elongation rate can be approximated as constants within the growth rate range of interest for cells exhibiting overflow metabolism (Dai et al., 2017), our model remains applicable to suboptimal growth conditions. This can be achieved by incorporating suboptimal ribosome allocation factors, lowering the parameter κt (which results in a larger ψ through Equation S32), to account for these influences. For convenience, we present results for optimal growth below, while all model results can be extended to cases of suboptimal ribosome allocation.

Regarding the mechanism by which cells sense and choose between respiration and fermentation, although the standard picture of overflow metabolism (Basan et al., 2015; Holms, 1996) presents a growth rate dependence of fermentation flux, it is the proteome efficiency of respiration and fermentation, rather than the growth rate, that a cell should sense directly. Due to stochasticity in gene expression and metabolic reactions, the cell growth rate may fluctuate within a cell cycle (Kiviet et al., 2014; Pavlou et al., 2025), and suboptimal factors related to ribosome allocation (Dai et al., 2017; Li et al., 2018) would further complicate the scheme if cells were sensing via growth rate. Essentially, to expedite cell growth and survive under evolutionary pressure, cells should adopt the optimal strategy by directly sensing and comparing proteome efficiencies between respiration and fermentation, choosing the pathway with higher efficiency. This is analogous to how microbes choose between two types of carbon sources in a mixture for nutrient uptake (Wang et al., 2019). Mechanistically, the cyclic AMP (cAMP)-cAMP receptor protein (CRP) system plays an important role in sensing proteome efficiency and executing the optimal strategy between respiration and fermentation (Basan et al., 2015; Towbin et al., 2017; Valgepea et al., 2010; Wehrens et al., 2023). However, the roles of additional unidentified regulators are required to fully elucidate this mechanism (Basan et al., 2015; Valgepea et al., 2010).

Appendix 4

Model perturbations

4.1 Overexpression of useless proteins

Here, we consider the case of overexpression of the protein encoded by the lacZ gene (i.e. ϕZ perturbation) in E. coli. Effectively, this limits the proteome by altering ϕmax:

(S48) ϕmaxLacZ overexpressionϕmaxϕZ,

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

(S49) λ(κA,ϕZ)={ϕmaxϕZφ/εr(κA)+ψ(κA)εr(κA)>εf(κA),ϕmaxϕZφ/εf(κA)+ψ(κA)εr(κA)<εf(κA),

and thus,

(S50) λ(κA,ϕZ)=λ(κA,0)(1ϕZ/ϕmax).

Obviously, κA(C) remains a constant (following Equation S42), while λC(ϕZ)λ(κA(C),ϕZ) and λmax(ϕZ)λ(κAmax,ϕZ) become functions of ϕZ:

(S51) {λC(ϕZ)=λC(0)(1ϕZ/ϕmax),λmax(ϕZ)=λmax(0)(1ϕZ/ϕmax).

In the homogeneous case, Jf(E) and Jr(E) follow:

(S52) {Jf(E)(κA,ϕZ)=φλ(κA,ϕZ)θ(λ(κA,ϕZ)λC(ϕZ)),Jr(E)(κA,ϕZ)=φλ(κA,ϕZ)[1θ(λ(κA,ϕZ)λC(ϕZ))].

Combined with Equations S50-S51, we have:

(S53) {Jf(E)(κA,ϕZ)=φλ(κA,ϕZ)θ(λ(κA,0)λC(0)),Jr(E)(κA,ϕZ)=φλ(κA,ϕZ)[1θ(λ(κA,0)λC(0))].

To compare with experiments, we assume that each kicat and κi follow the extrinsic noise with a CV of 25% specified in Appendix 3.3, and we neglect the noise on ϕZ and ϕmax. Combining Equations S45 and S51, λC(ϕZ) approximately follows a Gaussian distribution:

(S54) λC(ϕZ)N(μλC(ϕZ),σλC(ϕZ)2),

where μλC(ϕZ) and σλC(ϕZ) represent the mean and standard deviation of λC(ϕZ), with

(S55) {μλC(ϕZ)=μλC(0)(1ϕZ/ϕmax),σλC(ϕZ)=σλC(0)(1ϕZ/ϕmax).

Here, μλC(0) , σλC(0), λC(0) , λmax(0), and λ(κA,0) represent the parameters or variables free from ϕZ perturbation, just as those in Appendix 3.3. Since the noise on the multiplier term (i.e. 1ϕZ/ϕmax) is negligible, the CV of λC(ϕZ) (i.e. σλC(ϕZ)/μλC(ϕZ)) is unaffected by ϕZ. By combining Equations S46 and S48, we obtain the relations between the normalized energy fluxes and growth rate:

(S56) {Jf(E)(λ(κA,ϕZ),ϕZ)=12φλ(κA,ϕZ)[erf(λ(κA,ϕZ)μλC(ϕZ)2σλC(ϕZ))+1],Jr(E)(λ(κA,ϕZ),ϕZ)=12φλ(κA,ϕZ)[1erf(λ(κA,ϕZ)μλC(ϕZ)2σλC(ϕZ))],

where λ(κA,ϕZ), μλC(ϕZ), and σλC(ϕZ) follow Equations S50 and S55 accordingly. For a given value of ϕZ, i.e., ϕZ is fixed, then, λ(κA,ϕZ) changes monotonically with κA. Combining Equations S55-S56 and S30, we obtain the relation between the normalized fluxes Jr(N), Jf(N), and the growth rate (where ϕZ is a parameter):

(S57) {Jf(N)(λ,ϕZ)=φβf(A)λ[erf(λμλC(0)(1ϕZ/ϕmax)2σλC(0)(1ϕZ/ϕmax))+1],Jr(N)(λ,ϕZ)=φβr(A)λ[1erf(λμλC(0)(1ϕZ/ϕmax)2σλC(0)(1ϕZ/ϕmax))].

In Figure 2C. we show that the model predictions (Equation 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 kAcat, yet this distribution is fixed. For convenience of description, we still referred to it as fixed κA. Then, combining Equations S30, S50, S55 and S56, we get:

(S58) {Jf(N)(λ,ϕZ)=φβf(A)[erf(λ(κA,0)μλC(0)2σλC(0))+1]λ,Jr(N)(λ,ϕZ)=φβr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))]λ.

Here, λ(κA,0) remains unaltered as κA is fixed. Therefore, in this case, Jf(N) and Jr(N) are proportional to λ, where the slopes are both functions of κA. More specifically, the slope of Jf(N) is a monotonically increasing function of κA, while that of Jr(N) is a monotonically decreasing function of κA. In Figure 2B, we see that the model predictions (Equation 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 Equation S57 (where ϕZ is a variable). In a 3-D representation, these relations correspond to a surface. In Figure 2A, we show that the model predictions (Equation S57) match well with the experimental data (Basan et al., 2015).

4.2 Energy dissipation

In practice, energy dissipation disrupts the proportional relationship between energy demand and biomass production. Thus, Equation S25 becomes:

(S59) JE=rEJBM+wMcarbonm0,

where w represents the dissipation coefficient. In fact, maintenance energy contributes to energy dissipation, and we define the maintenance energy coefficient as w0. In bacteria, the impact of maintenance energy is roughly negligible, yet in tumor cells, it plays a much more significant role (Locasale and Cantley, 2010).

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

(S60) JE(N)=ηEλ+w.

Then, Equation S29 changes to:

(S61) {Jr(E)+Jf(E)=φλ+w,Jr(E)εr+Jf(E)εf=ϕmaxψλ.

Consequently, if εr>εf, the optimal growth strategy for the cell is:

(S62) {Jf(E)=0,Jr(E)=φλ+w,εr>εf,

and if εf>εr, the optimal growth strategy is:

(S63) {Jf(E)=φλ+w,Jr(E)=0.εr<εf.

Then, the growth rate becomes a bivariate function of both κA and w:

(S64) λ(κA,w)={ϕmaxw/εr(κA)φ/εr(κA)+ψ(κA)εr(κA)>εf(κA),ϕmaxw/εf(κA)φ/εf(κA)+ψ(κA)εr(κA)<εf(κA).

Clearly, κA(C) is still a constant, while λC(w)λ(κA(C),w) and λmax(w)λ(κAmax,w) become functions of w:

(S65) {λC(w)=λC(0){1w/[εr/f(κA(C))ϕmax]},λmax(w)=λmax(0){1w/[εf(κAmax)ϕmax]}.

For a cell population, in the homogeneous case, Jf(E) and Jr(E) follow:

(S66) {Jf(E)(κA,w)=[φλ(κA,w)+w]θ(λ(κA,w)λC(w)),Jr(E)(κA,w)=[φλ(κA,w)+w][1θ(λ(κA,w)λC(w))].

To compare with experiments, we assume the same extent of extrinsic noise in kicat (and thus κi) as that specified in Appendix 3.3. Combining Equations S45 and S65, λC(w) approximately follows a Gaussian distribution:

(S67) λC(w)N(μλC(w),σλC(w)2),

where μλC(w) and σλC(w) represent the mean and standard deviation of λC(w), and

(S68) {μλC(w)=μλC(0){1w/[εr/f(κA(C))ϕmax]},σλC(w)σλC(0){1w/[εr/f(κA(C))ϕmax]}.

Here, μλC(0) , σλC(0), λC(0) , λmax(0), and λ(κA,0) represent parameters or variables unaffected by energy dissipation. In fact, there is a distribution of values for εr/f(κA(C)). For approximation, we use the deterministic value of εr/f(κA(C)) in Equation S68, and then the CV of λC(w) remains largely unperturbed by w. Combining Equations S46, S66 and S67, we have:

(S69) {Jf(E)(λ(κA,w),w)=12(φλ(κA,w)+w)[erf(λ(κA,w)μλC(w)2σλC(w))+1],Jr(E)(λ(κA,w),w)=12(φλ(κA,w)+w)[1erf(λ(κA,w)μλC(w)2σλC(w))].

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

(S70) {Jf(N)(λ,w)=φλ+wβf(A)[erf(λμλC(0){1w/[εr/f(κA(C))ϕmax]}2σλC(0){1w/[εr/f(κA(C))ϕmax]})+1],Jr(N)(λ,w)=φλ+wβr(A)[1erf(λμλC(0){1w/[εr/f(κA(C))ϕmax]}2σλC(0){1w/[εr/f(κA(C))ϕmax]})].

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

4.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:

(S71) κtTranslationinhibitionκt/(ι+1),

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

(S72) ψ(κA,ι)=ι+1κt+1+ηa1+ηc2κA+ηa2+ηb+2ηc+ηd2κ1+ηb+ηcκ2+ηcκ3+ηc+ηdκ5+ia1,a2,b,c,dηiκi

First, we consider the case where maintenance energy is neglected, i.e., w0=0. In this case, the growth rate takes the following form:

(S73) λ(κA,ι)={ϕmaxφ/εr(κA)+ψ(κA,ι)εr(κA)>εf(κA),ϕmaxφ/εf(κA)+ψ(κA,ι)εr(κA)<εf(κA),

where λ(κA,0) and ψ(κA,0) represent the terms unaffected by translation inhibition. Thus, λC(ι)λ(κA(C),ι) and λmax(ι)λ(κAmax,ι) become functions of ι:

(S74) {λC(ι)=λC(0)φ/εr/f(κA(C))+ψ(κA(C),0)φ/εr/f(κA(C))+ψ(κA(C),ι),λmax(ι)=λmax(0)φ/εf(κAmax)+ψ(κAmax,0)φ/εf(κAmax)+ψ(κAmax,ι).

In the homogeneous case, Jf(E) and Jr(E) follow:

(S75) {Jf(E)(κA,ι)=φλ(κA,ι)θ(λ(κA,ι)λC(ι)),Jr(E)(κA,ι)=φλ(κA,ι)[1θ(λ(κA,ι)λC(ι))].

To compare with experiments, we assume that extrinsic noise exists in kicat and κi as specified in Appendix 3.3. Combining Equations S45 and S74, λC(ι) can be approximated by a Gaussian distribution:

(S76) λC(ι)N(μλC(ι),σλC(ι)2),

where μλC(ι) and σλC(ι) represent the mean and standard deviation of λC(ι), with

(S77) {μλC(ι)=μλC(0)φ/εr/f(κA(C))+ψ(κA(C),0)φ/εr/f(κA(C))+ψ(κA(C),ι),σλC(ι)σλC(0)φ/εr/f(κA(C))+ψ(κA(C),0)φ/εr/f(κA(C))+ψ(κA(C),ι).

Here, μλC(0) , σλC(0), ψ(κA(C),0), λC(0) and λmax(0) stand for the terms unaffected by translation inhibition. Essentially, there are distributions of values for εr/f(κA(C)), ψ(κA(C),0) and ψ(κA(C),ι). For approximation, we use the deterministic values of these terms in Equation S77, and then the CV of λC(ι) can be approximated by λC(0). Combining Equations S46, S75 and S76, we have:

(S78) {Jf(E)(λ(κA,ι),ι)=12φλ(κA,ι)[erf(λ(κA,ι)μλC(ι)2σλC(ι))+1],Jr(E)(λ(κA,ι),ι)=12φλ(κA,ι)[1erf(λ(κA,ι)μλC(ι)2σλC(ι))].

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

(S79) {Jf(N)(λ,ι)=φλβf(A)[erf(λμλC(ι)2σλC(ι))+1],Jr(N)(λ,ι)=φλβr(A)[1erf(λμλC(ι)2σλC(ι))],

where μλC(ι) and σλC(ι) follow Equation S77. The growth rate can also be perturbed by altering both κA and ι simultaneously. In this case, the relations among the energy fluxes, growth rate and ι still follow Equation S79 (here ι is a variable). The comparison between Equation S79 and experimental data (Basan et al., 2015) is shown in Appendix 1—figure 2D (3-D) and E(2-D). Overall, there is good consistency; however, there remains a noticeable discrepancy when ι is large (i.e. at high concentration of the translation inhibitor). This led us to consider the maintenance energy through the coefficient w0, which is small but may account for this discrepancy. Then, λ(κA,ι) changes into:

(S80) λ(κA,ι)={ϕmaxw0/εr(κA)φ/εr(κA)+ψ(κA,ι)εr(κA)>εf(κA),ϕmaxw0/εf(κA)φ/εf(κA)+ψ(κA,ι)εr(κA)<εf(κA),

while λC(ι)λ(κA(C),ι) and λmax(ι)λ(κAmax,ι) still follow Equation S74, though the forms of λC(0) and λmax(0) change to:

(S81) {λC(0)=ϕmaxw0/εr/f(κA(C))φ/εr/f(κA(C))+ψ(κA(C),0),λmax(0)=ϕmaxw0/εf(κAmax)φ/εf(κAmax)+ψ(κAmax,0).

In the homogeneous case, Jf(E) and Jr(E) follow:

(S82) {Jf(E)(κA,ι)=[φλ(κA,ι)+w0]θ(λ(κA,ι)λC(ι)),Jr(E)(κA,ι)=[φλ(κA,ι)+w0][1θ(λ(κA,ι)λC(ι))].

To compare with experiments, we assume that the extrinsic noise follows the specification in Appendix 3.3. Combining Equations S45, S74 and S81, λC(ι) approximately follows a Gaussian distribution:

(S83) λC(ι)N(μλC(ι),σλC(ι)2)

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

(S84) {Jf(E)(λ(κA,ι),ι)=12(φλ(κA,ι)+w0)[erf(λ(κA,ι)μλC(ι)2σλC(ι))+1],Jr(E)(λ(κA,ι),ι)=12(φλ(κA,ι)+w0)[1erf(λ(κA,ι)μλC(ι)2σλC(ι))].

Thus, for a given ι, λ(κA,ι) changes monotonically with κA. Combining Equations S30 and S84, we have (here ι is a parameter):

(S85) {Jf(N)(λ,ι)=φλ+w0βf(A)[erf(λμλC(ι)2σλC(ι))+1].Jr(N)(λ,ι)=φλ+w0βr(A)[1erf(λμλC(ι)2σλC(ι))].

The growth rate and fluxes can also be perturbed by altering both κA and ι simultaneously. The relations among the energy fluxes, growth rate, and ι would still follow Equation S85, except that ι is now regarded as a variable. Assuming a small amount of maintenance energy by assigning w0=2.5 (h1), we find that the experimental results (Basan et al., 2015) agree quantitatively well with the model predictions (Figure 3C and D).

Appendix 5

Overflow metabolism in substrates other than Group A carbon sources

Due to the topology of the metabolic network, for cells using Group A carbon sources, the behavior of overflow metabolism follows Equation 5 (or Equation S47) 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 and agrees quantitatively with experiments. However, further analysis is required for cells using substrates other than Group A sources due to the topological differences in carbon utilization (Wang et al., 2019). In principle, substrates entering from glycolysis or the points before acetyl-CoA are potentially involved in overflow metabolism, while those joining from the TCA cycle are not relevant to this behavior. Still, mixed carbon sources are likely to induce a different profile of overflow metabolism, as long as there is a carbon source derived from glycolysis.

5.1 Pyruvate

The coarse-grained model for pyruvate utilization is shown in Figure 3E. Here, nodes M1, M2, M3, M4, M5 follow the descriptions in Appendix 3.1. Each biochemical reaction follows Equation S5 with bi=1 except that 2M2→M1 and M3+M5→M4. By applying flux balance to the stoichiometric fluxes, combining with Equation S8, we have:

(S86) {Φpyξpy=Φ7ξ7+Φ8ξ8,Φ7ξ7=2Φ9ξ9+Φ5ξ5+Φa2ξa2,Φ9ξ9=Φa1ξa1,Φ8ξ8=Φ3ξ3+Φ6ξ6+Φbξb,Φ5ξ5+Φ4ξ4=Φ3ξ3+Φdξd,Φ3ξ3=Φ4ξ4+Φcξc.

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

(S87) β8Φ8ξ8+β3Φ3ξ3+β4Φ4ξ4+β6Φ6ξ6+βa1Φa1ξa1β7Φ7ξ7β9Φ9ξ9=JE

where β7=1, β8=2, β3=2, β4=6, β6=1, β9=6,βa1=4 for E. coli (Neidhardt et al., 1990; Sauer et al., 2004), and JE follows Equation S25. By applying the substitutions specified in Equations S9, S12, S14-S18, combined with Equations S4, S10, S22, S23, S25, S86-S87, and the constraint of proteome resource allocation, we have:

(S88) {ϕpyκpy=ϕ7κ7+ϕ8κ8,ϕ7κ7=2ϕ9κ9+ϕ5κ5+ϕa2κa2,ϕ9κ9=ϕa1κa1ϕ8κ8=ϕ3κ3+ϕ6κ6+ϕbκbϕ3κ3=ϕ4κ4+ϕcκcϕ5κ5+ϕ4κ4=ϕ3κ3+ϕdκdϕa1κa1=ηa1λ,ϕa2κa2=ηa2λ,ϕbκb=ηbλ,ϕcκc=ηcλ,ϕdκd=ηdλ,β8ϕ8κ8+β3ϕ3κ3+β4ϕ4κ4+β6ϕ6κ6+βa1ϕa1κa1β7ϕ7κ7β9ϕ9κ9=JE(N),JE(N)=ηEλ,λ=ϕRκt,Jr(N)=ϕ4κ4,Jf(N)=ϕ6κ6,ϕR+ϕpy+ϕ3+ϕ4+ϕ5+ϕ6+ϕ7+ϕ8+ϕ9+ϕa1+ϕa2+ϕb+ϕc+ϕd=ϕmax,

where ηE=rE[iri/NEPicarbon]1. κi is approximately a constant which follows Equation S20 for each of the intermediate node. The substrate quality of κpy varies with the external concentration of pyruvate ([py]),

(S89) κpyrproteinrcarbonkpycatmEpy[py][py]+Kpym0.

From Equation S88, all ϕi can be expressed by Jr(N), Jf(N), and λ:

(S90) {ϕpy=[(2ηa1+ηa2+ηb+2ηc+ηd)λ+Jr(N)+Jf(N)]/κpy,ϕ7=(2ηa1+ηa2+ηc+ηd)λ/κ7,ϕ9=ηa1λ/κ9ϕ8=[Jr(N)+Jf(N)+(ηb+ηc)λ]/κ8ϕ3=(Jr(N)+ηcλ)/κ3,ϕ4=Jr(N)/κ4,ϕ5=(ηc+ηd)λ/κ5,ϕ6=Jf(N)/κ6,ϕi=ηiλ/κi(i=a1,a2,b,c,d).

By substituting Equation S90 into Equation S88, we have:

(S91) {Jr(E,py)+Jf(E,py)=φpyλ,Jr(E,py)εr(py)+Jf(E,py)εf(py)=ϕmaxψpyλ.

Here, Jr(E,py) and Jf(E,py) stand for the normalized energy fluxes of respiration and fermentation, respectively, with

(S92) {Jr(E,py)=βr(py)Jr(N),Jf(E,py)=βf(py)Jf(N).

where βr(py)=β3+β4+β8 and βf(py)=β6+β8, with βr(py)=10 and βf(py)=3 for E. coli. The coefficients εr(py) and εf(py) represent the proteome efficiencies for energy biogenesis using pyruvate in respiration and fermentation pathways, respectively, with

(S93) {εr(py)=βr(py)1/κpy+1/κ8+1/κ3+1/κ4,εf(py)=βf(py)1/κpy+1/κ8+1/κ6.

ψpy1 is the proteome efficiency for biomass generation using pyruvate in the biomass synthesis pathway, with

(S94) ψpy=1κt+1+ηa1+ηcκpy+1ηb+ηa1κ7+ηb+ηcκ8+ηa1κ9+ηcκ3+ηc+ηdκ5+ia1,a2,b,c,dηiκi

φpy is an energy demand coefficient (a constant), with

(S95) φpyηE+β7(1ηb+ηa1)+β9ηa1β8(ηc+ηb)β3ηcβa1ηa1,

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

(S96) λ(κpy)={ϕmaxφpy/εr(py)(κpy)+ψpy(κpy)εr(py)(κpy)>εf(py)(κpy),ϕmaxφpy/εf(py)(κpy)+ψpy(κpy)εr(py)(κpy)<εf(py)(κpy).

When κpy is very small, combined with Equation S93, then,

(S97) {εr(py)(κpy0)βr(py)κpy,εf(py)(κpy0)βf(py)κpy.

Obviously, βr(py)βf(py), and hence

(S98) εr(py)(κpy0)>εf(py)(κpy0).

As long as

(S99) βr(py)βf(py)κpy(ST)<βf(py)(1κ8+1κ3+1κ4)βr(py)(1κ8+1κ6),

where the superscript ‘(ST)’ stands for the saturated concentration, then,

(S100) εr(py)(κpy(ST))<εf(py)(κpy(ST)),

and there exists a critical value of κpy, denoted as κpy(C), with

(S101) {εr(py)(κpy(C))=εf(py)(κpy(C))=βr(py)βf(py)1/κ3+1/κ41/κ6=β3+β4β61/κ3+1/κ41/κ6,λC(py)λ(κpy(C))=ϕmaxφpy/εr/f(py)(κpy(C))+ψpy(κpy(C)).

Here, λC(py) is the growth rate at the transition point, and εr/f(py) stands for either εr(py) or εf(py). In Appendix 1—figure 2H, we show the dependencies of εr(py)(κpy), εf(py)(κpy) and λ(κpy) on κpy in a 3-dimensional form. In the homogeneous case, Jf(E,py) and Jr(E,py) follow:

(S102) {Jf(E,py)=φpyλθ(λλC(py)),Jr(E,py)=φpyλ[1θ(λλC(py))].

Defining λmax(py)=λ(κpy(ST)), and then, [0,λmax(py)] is the relevant range of the x axis. To compare with experiments, we assume the same extent of extrinsic noise in kicat as specified in Appendix 3.3. Then, λC(py) approximately follows a Gaussian distribution:

(S103) λC(py)N(μλC(py),σλC(py)2),

where μλC(py) and σλC(py) stand for the mean and standard deviation of λC(py). Then, the relations between the normalized energy fluxes and growth rate are:

(S104) {Jf(E,py)(λ)=12φpyλ[erf(λμλC(py)2σλC(py))+1],Jr(E,py)(λ)=12φpyλ[1erf(λμλC(py)2σλC(py))].

Combined with Equation S92, we have:

(S105) {Jf(N)(λ)=φpy2βf(py)λ[erf(λμλC(py)2σλC(py))+1],Jr(N)(λ)=φpy2βr(py)λ[1erf(λμλC(py)2σλC(py))].

In Figure 3F, we show that the model predictions (Equation S105) align quantitatively with the experimental results (Holms, 1996).

5.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 1—figure 2A. This model can be used to analyze mixtures with one or multiple types of extracellular amino acids. Here, Equations S21, S22, S24 and S25 still apply, but Equation S23 changes to (the case of i=a1 remains the same as Equation S23):

(S106) ΦiξiNEPicarbon+ΦiξiNPicarbon=riJBM(i=a2,b,c,d).

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

(S107) NPicarbonNEPicarbon.

In the case where all 21 types of amino acids are present and each is at saturated concentration (denoted as ‘21AA’), we have:

(S108) {ϕAκA=ϕ1κ1+ϕa1κa1,2ϕ1κ1=ϕ2κ2+ϕ5κ5+ϕa2κa2,ϕ2κ2=ϕ3κ3+ϕ6κ6+ϕbκb,ϕ5κ5+ϕ4κ4=ϕ3κ3+ϕdκd,ϕ3κ3=ϕ4κ4+ϕcκc,ϕa1κa1=ηa1λ,ϕa2κa2+ϕa2(21AA)κa2(21AA)=ηa2λ,ϕbκb+ϕb(21AA)κb(21AA)=ηbλ,ϕcκc+ϕc(21AA)κc(21AA)=ηcλ,ϕdκd+ϕd(21AA)κd(21AA)=ηdλ,β1ϕ1κ1+β2ϕ2κ2+β3ϕ3κ3+β4ϕ4κ4+β6ϕ6κ6+βa1ϕa1κa1=JE(N),JE(N)=ηEλ,λ=ϕRκt,Jr(N)=ϕ4κ4,Jf(N)=ϕ6κ6,ϕR+ϕA+i6ϕi+ja1,a2,b,c,dϕj+ϕa2(21AA)+ϕb(21AA)+ϕc(21AA)+ϕd(21AA)=ϕmax,

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

(S109) ϕi=0(i=a2,b,c,d)

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

(S110) {1/κa2(21AA)<1/κa2+1/(2κ1)+1/(2κA),1/κb(21AA)<1/κb+1/κ2+1/(2κ1)+1/(2κA),1/κc(21AA)<1/κc+1/κ5+1/κ3+1/κ2+1/κ1+1/κA,1/κd(21AA)<1/κd+1/κ5+1/(2κ1)+1/(2κA).

In practice, the requirement for proteome efficiency when using amino acids is even higher, since the biomass synthesis pathway is accompanied by energy biogenesis for Group A carbon sources, but not for amino acids. Combining Equations S108 and S109, we have:

(S111) {Jr(E)+Jf(E)=φ21AAλ,Jr(E)εr+Jf(E)εf=ϕmaxψ21AAλ,

where Jr(E), Jf(E) follow Equation S30, while εr and εf follow Equation S31. ψ21AA1 is the proteome efficiency for biomass generation in the biomass synthesis pathway under this nutrient condition, with

(S112) ψ21AA=1κt+ηa1κA+ηa1κa1+ηa2κa2(21AA)+ηbκb(21AA)+ηcκc(21AA)+ηdκd(21AA)

φ21AA is an energy demand coefficient, with

(S113) φ21AAηEβa1ηa1

Combining Equations S111 and S31, the formula for the growth rate is:

(S114) λ(κA)={λr(21AA)=ϕmaxφ21AA/εr(κA)+ψ21AA(κA)εr(κA)>εf(κA),λf(21AA)=ϕmaxφ21AA/εf(κA)+ψ21AA(κA)εr(κA)<εf(κA).

In fact, Equations S37-S42 still apply. εr/f(κA(C)) satisfies Equation S43, while λC(21AA)λ(κA(C)) and λmax(21AA)λ(κAmax) are:

(S115) {λC(21AA)=ϕmaxφ21AAεr/f(κA(C))+ψ21AA(κA(C)),λmax(21AA)=ϕmaxφ21AAεf(κAmax)+ψ21AA(κAmax).

When extrinsic noise is taken into account, λC(21AA) approximately follows a Gaussian distribution:

(S116) λC(21AA)N(μλC(21AA),σλC(21AA)2),

and the normalized fluxes Jr(N), Jf(N) change to:

(S117) {Jf(N)(λ)=φ21AAβf(A)λ[erf(λμλC(21AA)2σλC(21AA))+1],Jr(N)(λ)=φ21AAβr(A)λ[1erf(λμλC(21AA)2σλC(21AA))].

The above analysis can be extended to cases where a Group A carbon source is mixed with arbitrary combinations of amino acids. Equations S111, S114-S117 would remain in a similar form, while Equations S112-S113 would change depending on the combinations of amino acid. In Appendix 1—figure 2B and C, we compare model predictions (see also Appendix 8.2 and Equation S157) with experimental data (Basan et al., 2015; Wallden et al., 2016) from mixtures of 21 or 7 types of amino acids along with a Group A carbon source, demonstrating quantitative agreement. Additionally, the increase in the critical threshold of growth rate for the growth rate-dependent fermentation flux in mixtures with extracellular amino acids (i.e. λC(21AA),λC(7AA)>λC; see Appendix 1—figure 2C) has also been observed in other experimental findings (Peebo et al., 2015).

Appendix 6

Enzyme allocation upon perturbations

6.1 Carbon limitation within Group A carbon sources

In Equation S28, we present the model predictions for the dependencies of enzyme proteomic mass fractions on growth rate and energy fluxes. To compare with experiments, we assume the same extent of extrinsic noise in kicat as specified in Appendix 3.3. Relative protein expression data for enzymes within glycolysis and the TCA cycle are available from existing studies and are comparable to the ϕ1-ϕ4 enzymes of our model (Figure 1B). Upon κA perturbation, κA is a variable while w0 is fixed (see Appendix 2.5). Combining Equations S28 and S47 (with w0=0), we obtain:

(S118) {ϕ1=λκ1{φ(βr(A)βf(A))2βr(A)βf(A)[erf(λμλC2σλC)+1]+φβr(A)+ηa2+ηb+2ηc+ηd2},ϕ2=λκ2{φ(βr(A)βf(A))βr(A)βf(A)[erf(λμλC2σλC)+1]+2φβr(A)+ηb+ηc},ϕ3=λκ3{φβr(A)[1erf(λμλC2σλC)]+ηc},ϕ4=λκ4φβr(A)[1erf(λμλC2σλC)].

In Appendix 1—figure 3C and D, we show the comparisons between model predictions (Equation S118, w0=0) and experimental data (Hui et al., 2015), which are consistent overall. We then consider the influence of maintenance energy as specified in Appendix 4.2. Here, we continue to choose w0=2.5 (h1) as previously adopted in Appendix 4.3. Thus, Equation S28 still holds. Combined with Equation S85 under the condition that ι=0, we have:

(S119) {ϕ1=12κ1{φλ+w0βf(A)[erf(λμλC2σλC)+1]+φλ+w0βr(A)[1erf(λμλC2σλC)]+(ηa2+ηb+2ηc+ηd)λ},ϕ2=1κ2{φλ+w0βf(A)[erf(λμλC2σλC)+1]+φλ+w0βr(A)[1erf(λμλC2σλC)]+(ηb+ηc)λ},ϕ3=1κ3{φλ+w0βr(A)[1erf(λμλC2σλC)]+ηcλ},ϕ4=1κ4φλ+w0βr(A)[1erf(λμλC2σλC)].

In Figure 4A–B, we show that the model predictions (Equation S119, w0=2.5(h1)) generally agree with the experiments (Hui et al., 2015). However, there are different basal expressions of these enzymes, likely due to living demands other than cell proliferation, such as preparation for starvation (Mori et al., 2017) or changes in the type of the nutrient (Basan et al., 2020; Kussell and Leibler, 2005).

6.2 Overexpression of useless proteins

In the case of ϕZ perturbation under each nutrient condition with fixed κA (see Appendix 4.1), we consider the same extent of extrinsic noise in kicat as specified in Appendix 3.3. The relation between enzyme allocation and growth rate can be obtained by combining Equations S28 and S58 (with w0=0):

(S120) {ϕ1=λ2κ1{(ηa2+ηb+2ηc+ηd)+φβr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))]+φβf(A)[erf(λ(κA,0)μλC(0)2σλC(0))+1]},ϕ2=λκ2{(ηb+ηc)+φβr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))]+φβf(A)[erf(λ(κA,0)μλC(0)2σλC(0))+1]},ϕ3=λκ3{φβr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))]+ηc},ϕ4=λκ4{φβr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))]}.

Here λ(κA,0) is the growth rate for ϕZ=0, and thus it is a parameter rather than a variable. The growth rate is defined as λ(κA,ϕZ), which follows Equation S50. Thus, ϕi is proportional to the growth rate λ. In Appendix 1—figure 3E and F, we observe that the model predictions (Equation S120) generally agree with the experiments (Basan et al., 2015). Next, we consider the influence of maintenance energy with w0=2.5 (h1). Combining Equations S28, S58 and S85 (with ι=0), we get:

(S121) {ϕ1=w02κ1{1βf(A)[erf(λ(κA,0)μλC(0)2σλC(0))+1]+1βr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))]}+λ2κ1{φβf(A)[erf(λ(κA,0)μλC(0)2σλC(0))+1]+φβr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))]+(ηa2+ηb+2ηc+ηd)},ϕ2=w0κ2{1βf(A)[erf(λ(κA,0)μλC(0)2σλC(0))+1]+1βr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))]}+λκ2{φβf(A)[erf(λ(κA,0)μλC(0)2σλC(0))+1]+φβr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))]+(ηb+ηc)},ϕ3=λκ3{φβr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))]+ηc}+w0κ31βr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))],ϕ4=λκ4φβr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))]+w0κ41βr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))].

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 Figure 4C–D and Appendix 1—figure 3I-J, we show that the model predictions (Equation S121) agree quantitively with the experimental data (Basan et al., 2015).

6.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 Equations S28 and S70. However, since w is explicitly present in Equation S70, it is necessary to reduce this variable to obtain the growth rate dependence of enzyme allocation. From Equation S64, we have:

(S122) λ(κA,w)=λ(κA,0){1wϕmax[1εr(κA)θ(εf(κA)εr(κA))(1εr(κA)1εf(κA))]}.

Here, λ(κA,0)λ(κA,w=0) (satisfying Equation S64) is a parameter rather than a variable. ‘θ’ stands for the Heaviside step function. Thus, we have:

(S123) w(λ)=ϕmax[1λ/λ(κA,0)][1/εr(κA)θ(εf(κA)εr(κA))(1/εr(κA)1/εf(κA))],

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

Combining Equations S28, S70 and S123, we get:

(S124) {ϕ1=12κ1{erf(λμλC(0)[1w(λ)εr/f(κA(C))ϕmax]2σλC(0)[1w(λ)εr/f(κA(C))ϕmax])[φλ+w(λ)βf(A)φλ+w(λ)βr(A)]+[φλ+w(λ)](1βf(A)1βr(A)+2βr(A))+(ηa2+ηb+2ηc+ηd)λ},ϕ2=1κ2{erf(λμλC(0)[1w(λ)εr/f(κA(C))ϕmax]2σλC(0)[1w(λ)εr/f(κA(C))ϕmax])[φλ+w(λ)βf(A)φλ+w(λ)βr(A)]+[φλ+w(λ)](1βf(A)1βr(A)+2βr(A))+(ηb+ηc)λ},ϕ3=1κ3{[φλ+w(λ)]βr(A)[1erf(λμλC(0)[1w(λ)εr/f(κA(C))ϕmax]2σλC(0)[1w(λ)εr/f(κA(C))ϕmax])]+ηcλ},ϕ4=1κ4[φλ+w(λ)]βr(A)[1erf(λμλC(0)[1w(λ)εr/f(κA(C))ϕmax]2σλC(0)[1w(λ)εr/f(κA(C))ϕmax])],

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

(S125) (1wεr/f(κA)ϕmax)/(1wεr/f(κA)ϕmax)(1wεr/f(κA(C))ϕmax)(1wεr/f(κA(C))ϕmax)1

then we have:

(S126) {Jf(N)(λ,w)=φλ+wβf(A)[erf(λ(κA,0)μλC(0)2σλC(0))+1],Jr(N)(λ,w)=φλ+wβr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))],

and thus:

(s127) {ϕ1=12κ1{[φλ+w(λ)](1βf(A)1βr(A))[erf(λ(κA,0)μλC(0)2σλC(0))+1]+2βr(A)[φλ+w(λ)]+(ηa2+ηb+2ηc+ηd)λ},ϕ2=1κ2{[φλ+w(λ)](1βf(A)1βr(A))[erf(λ(κA,0)μλC(0)2σλC(0))+1]+2βr(A)[φλ+w(λ)]+(ηb+ηc)λ},ϕ3=1κ3{[φλ+w(λ)]βr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))]+ηcλ},ϕ4=1κ4[φλ+w(λ)]βr(A)[1erf(λ(κA,0)μλC(0)2σλC(0))],

Note that in Equation S123, w is a linear function of λ with a negative slope. Thus ϕi exhibits a linear relation with λ when Equation S125 is satisfied (see Equation S127). In fact, the slope of ϕ4 is certainly negative (combining Equations S64, S123 and S127), while the sign of the slope for 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 ϕ1 is negative, then the slope sign of ϕ2 is surely negative. In Appendix 1—figure 3K-N, we show that our model results agree well with the experimental data (Basan et al., 2015; Equation S127).

Appendix 7

Other aspects of the model

7.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 1—figure 2F. Here, nodes M6, M7 represent GA3P and DHAP, respectively. Other nodes follow the descriptions specified in Appendix 3.1. Each biochemical reaction follows Equation S5 with bi=1 except that M1→M6+M7 and M3+M5→M4. By applying flux balance to the stoichiometric fluxes, combined with Equation S8, we obtain:

(S128) {ΦAξA=Φ1ξ1+Φa1ξa1,Φ11ξ11=Φ10ξ10+Φ1ξ1,Φ10ξ10=Φ1ξ1,Φ11ξ11=Φ2ξ2+Φ5ξ5+Φa2ξa2,Φ2ξ2=Φ3ξ3+Φ6ξ6+Φbξb,Φ5ξ5+Φ4ξ4=Φ3ξ3+Φdξd,Φ3ξ3=Φ4ξ4+Φcξc.

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

(S129) {ϕAκA=ϕ1κ1+ϕa1κa1,ϕ11κ11=ϕ10κ10+ϕ1κ1,ϕ10κ10=ϕ1κ1,ϕ11κ11=ϕ2κ2+ϕ5κ5+ϕa2κa2,ϕ2κ2=ϕ3κ3+ϕ6κ6+ϕbκb,ϕ5κ5+ϕ4κ4=ϕ3κ3+ϕdκd,ϕ3κ3=ϕ4κ4+ϕcκc,ϕa1κa1=ηa1λ,ϕa2κa2=ηa2λ,ϕbκb=ηbλ,ϕcκc=ηcλ,ϕdκd=ηdλ,β1ϕ1κ1+β2ϕ2κ2+β3ϕ3κ3+β4ϕ4κ4+β6ϕ6κ6+βa1ϕa1κa1=JE(N),JE(N)=ηEλ,λ=ϕRκt,Jr(N)=ϕ4κ4,Jf(N)=ϕ6κ6,ϕR+ϕA+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ϕ6+ϕ7+ϕ8+ϕa1+ϕa2+ϕb+ϕc+ϕd=ϕmax.

Then, Equation S28 still holds, while ϕ10 and ϕ11 are:

(S130) {ϕ10=[Jr(N)+Jf(N)+(ηa2+ηb+2ηc+ηd)λ]/2κ10(2κ10),ϕ11=[Jr(N)+Jf(N)+(ηa2+ηb+2ηc+ηd)λ]/κ11.

By substituting Equations S28 and S130 into Equation S129, we get:

(S131) {Jr(E)+Jf(E)=φλ,Jr(E)εr(dt)+Jf(E)εf(dt)=ϕmaxψdtλ,

where ‘dt’ stands for details. Equations S30 and S33 still hold. εr(dt) and εf(dt) represent the proteome efficiencies for energy biogenesis in the respiration and fermentation pathways, respectively, with

(S132) {εr(dt)=βr(A)1/κA+1/κ1+1/κ10+2/κ11+2/κ2+2/κ3+2/κ4,εf(dt)=βf(A)1/κA+1/κ1+1/κ10+2/κ11+2/κ2+2/κ6.

ψdt1 is the proteome efficiency for biomass generation in the biomass synthesis pathway, with

(S133) ψdt=1κt+1+ηa1+ηc2κA+(ηa2+ηb+2ηc+ηd)(12κ1+12κ10+1κ11)+ηb+ηcκ2+ηcκ3+ηc+ηdκ5+ia1,a2,b,c,dηiκi.

7.2 Estimation of the in vivo enzyme catalytic rates

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

(S134) {κ1=[Jr(N)+Jf(N)+(ηa2+ηb+2ηc+ηd)λ]/(2ϕ1),κ2=[Jr(N)+Jf(N)+(ηb+ηc)λ]/ϕ2,κ3=(Jr(N)+ηcλ)/ϕ3,κ4=Jr(N)/ϕ4,κ5=(ηc+ηd)λ/ϕ5,κ6=Jf(N)/ϕ6κ10=[Jr(N)+Jf(N)+(ηa2+ηb+2ηc+ηd)λ]/(2ϕ10),κ11=[Jr(N)+Jf(N)+(ηa2+ηb+2ηc+ηd)λ]/ϕ11.

Here, Jr(N), Jf(N), λ and ϕi (i=1-6,10-11) are measurable from experiments (see Appendix 9.1 and Appendix 1—table 2). Thus, we can obtain the in vivo values of κi from Equation S134. Combined with Equations S17 and S20, we have

(S135) kicat=rcarbonrproteinmEimcarbonκi[iri/NEPicarbon]

Equation S135 is the in vivo result for the enzyme catalytic rate. In Appendix 1—figure 2G, we show a comparison between in vivo and in vitro results for kcat 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 vitro data as a substitute when there were gaps.

7.3 Comparison with existing models that illustrate experimental results

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

(S136) Jin(N)JA(N)=ϕAκA.

Combined with the first equation in Equation S28 and Equation S30, we obtain:

(S137) Jin(N)ϑλ=Jr(E)βr(A)+Jf(E)βf(A),

where ϑ=ηa1+ηc+(ηa2+ηb+ηd)/2. Evidently, βr(A), βf(A) and ϑ are constant parameters. In this subsection, we highlight the major differences between our model presented in Appendix 3 and existing models that illustrate the growth rate dependence of fermentation flux in the standard picture of overflow metabolism (Basan et al., 2015; Holms, 1996; Meyer et al., 1984; Nanchen et al., 2006).

Based on the modeling principles rather than the detailed mechanisms, there are two major classes of existing models that can illustrate experimental results. Both classes of models regard the proteome efficiencies εr and εf as constants, with εf>εr if used, or follow functionally equivalent propositions. However, in our model, εr and εf are both functions of κA, which vary significantly upon nutrient perturbation, with εr(κA0)>εf(κA0) and εr(κAmax)<εf(κAmax) (see Equations S38, S40-S41). Furthermore, there are significant differences in the modeling and optimization principles, as listed below.

The first class of models (Chen and Nielsen, 2019; Majewski and Domach, 1990; Shlomi et al., 2011; Varma and Palsson, 1994; Vazquez et al., 2010; Vazquez and Oltvai, 2016; Zhuang et al., 2011) optimize the ratio of biomass outflow to carbon influx λ/Jin(N), 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 energy biogenesis per unit carbon, to optimize the ratio λ/Jin(N), 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 toward fermentation in a threshold-linear response, since they consider εf>εr. Our model is significantly different from this class of models in the optimization principle, as we purely optimize the cell growth rate for a given nutrient condition, without imposing a special constraint on the carbon influx.

The second class of models, represented by Basan et al., 2015, also adopt the optimization of λ/Jin(N) 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 (although it can be derived by optimizing λ/Jin(N)). In fact, Equations S29 and S137 in our model are very similar in form to those in Basan et al., 2015, yet there are critical differences, which we list below. In Equation S29, by regarding Jr(E) and Jf(E) as the two variables in a system of linear equations, we obtain the following expressions:

(S138) {Jr(E)=ϕmax(ψ+φεf)λ1εr1εf,Jf(E)=(ψ+φεr)λϕmax1εr1εf.

In Basan et al., 2015, Equation S138 is considered to be the relation between Jr/f(E) and λ upon nutrient (and thus Jin(N)) perturbation, while εr and εf are regarded as constants throughout the perturbation. By contrast, in our model, Equation S138 serves as a constraint under a given nutrient condition with fixed κA, and is not relevant to nutrient perturbation. For wild-type strains, if εr(κA)>εf(κA) (or vice versa), then the solution for optimal growth is Jr(E)(κA)=φλ(κA) and Jf(E)(κA)=0, with λ(κA)=εr(κA)ϕmaxφ+εr(κA)ψ(κA). This solution, which satisfies Equation S138, corresponds to a point rather than a line in the relation between growth rate λ and normalized energy flux Jr/f(E) upon κA perturbation.

Appendix 8

Probability density functions of variables and parameters

8.1 Probability density function of κi

Enzyme catalysis is crucial for the survival of living organisms, as it significantly accelerates biochemical reactions by reducing the energy barrier between the substrate and product (Nelson and Cox, 2008). However, the maximal turnover rate of enzymes, kcat, varies notably between in vivo and in vitro measurements (Davidi et al., 2016). Recent studies suggest that differences in the aquatic medium are the primary cause of this variation (Davidi et al., 2016; García-Contreras et al., 2012). In particular, potassium and phosphate concentrations have a significant influence on kcat (García-Contreras et al., 2012), and these concentrations exhibit some degree of variation among cell populations under intracellular conditions (García-Contreras et al., 2012). For simplicity, we assume that the turnover rate of each enzyme Ei, kicat, follows a Gaussian distribution N(μkicat,σkicat2) with kicat>0 among cells (representing extrinsic noise [Elowitz et al., 2002], denoted as χext). The probability density function of kicat is then given by:

(S139) kicatN(x;μkicat,σkicat2)={l1σkicat2πe12(xμkicatσkicat)2,x0.0,x<0.

When the CV of the kicat distribution (i.e. σkicat/μkicat) is less than 1/3 , N(x;μkicat,σkicat2) is almost identical to N(μkicat,σkicat2). In this case, 1/kicat follows the positive inverse of Gaussian (IOG) distribution, and the probability density function is:

(S140) IOG(x;μ1/kicat,ζ1/kicat)={ζ1/kicat2πx4exp(12ζ1/kicat(xμ1/kicat)2x2μ1/kicat2),x0,0,x<0,

where ζ1/kicat=1/σkicat2 is the shape parameter, and μ1/kicat=1/μkicat is the mean.

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

(S141) {Xi(t=0)=0,dXi/dt=αi+αiΓi(t),TΘ=inf{t>0|Xi(t)=Θ}.

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

(S142) IG(x;μ1/kicat,ζ1/kicat)={ζ1/kicat2πx3exp(12ζ1/kicat(xμ1/kicat)2xμ1/kicat2),x0,0,x<0,

where ζ1/kicat=Θ/kicat is the shape parameter, and μ1/kicat=1/kicat represents the mean. The variance of this distribution is σ1/kicat2μ1/kicat3/ζ1/kicat=1/[Θ(kicat)2]. Thus, we can obtain the CV:

(S143) σ1/kicat/μ1/kicat=Θ12,

which is inversely proportional to 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):

(S144) χtot2=χint2+χext2.

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

(S145) IOG(x;μ1/kicat,ζ1/kicat)CV1N(μ1/kicat,σ1/kicat2),

with a variance of σ1/kicat2=μ1/kicat4/ζ1/kicat, and

(S146) IG(x;μ1/kicat,ζ1/kicat)CV1N(μ1/kicat,σ1/kicat2),

with a variance of σ1/kicat2=μ1/kicat3/ζ1/kicat. Rigorously, we show below that IG(x;μ,ζ) shrinks to be N(μ,μ3/ζ) when the CV is small. For the IG distribution, the characteristic function of the variable x is given by Folks and Chhikara, 1978; Kampen, 1992:

(S147) G(k)=eikxIG(x;μ,ζ)dx=exp{ζμ[112iμ2kζ]},

and therefore,

(S148) IG(x;μ,ζ)=12πeikxG(k)dk.

When the variance σ2μ3/ζ is very small, we essentially require 2μ2k/ζ=2σ2k/μ1, and then 12iμ2kζ1μ2ζki+μ42ζ2k2. Thus,

(S149) {G(k)exp(μkiμ32ζk2),IG(x;μ,ζ)ζ2πμ3exp(ζ(xμ)22μ3)=N(μ,μ3ζ).

This leads to:

(S150) limσ0IG(x;μ,ζ)=N(μ,μ3/ζ).

In fact, intrinsic noise does affect the short-term measurement of enzyme catalytic rate and growth rate at the single-cell level. However, its contribution in the long term is averaged out and thus becomes negligible. For simplicity, we approximate χtotχext. Combined with Equations S145-S146, it is straightforward to verify that 1/kicat shares roughly the same CV as kicat:

(S151) σ1/kicat/μ1/kicat=σkicat/μkicat.

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

(S152) 1/kicatN(μ1/kicat,σ1/kicat2).

Using the properties of Gaussian distributions, for a series of constant real numbers γi, the summation of γi/kicat, which we define as Ξi=1nγi/kicat , follows a Gaussian distribution (Kampen, 1992):

(S153) ΞN(μΞ,σΞ2),

with μΞ=i=1nγiμ1/kicat and σΞ2=i=1n(γiσ1/kicat)2. The relation between κi and kicat is shown in Equation S12. To optimize cell growth rate, each κi of the intermediate nodes satisfies Equation S20, while κA satisfies Equation S27. Thus, for a given nutrient condition ([A] is fixed), all the ratios kicat/κi are constants. Combined with Equations S139, S145-S146, and S152, the distributions of all κi and 1/κi can be approximated as Gaussian distributions:

(S154) {κiN(μκi,σκi2),1κiN(μ1/κi,σ1/κi2),

where μκi and μ1/κi are the means of κi and 1/κi, and σκi and σ1/κi are their standard deviations. Using the properties of Gaussian distributions, combined with Equation S31, S32, S36, S42-S43, S145-S146 and S153, εr, εf, ψ, λr, λf, κA(C) and λC also roughly follow Gaussian distributions.

8.2 Probability density function of the growth rate λ

From Appendix 8.1, we note that λr and λf (see Equation S36) roughly follow Gaussian distributions, with

(S155) {λrN(μλr,σλr2),λfN(μλf,σλf2),

where μλr/f and σλr/f represent the mean and standard deviation, respectively. We further assume that the correlation between λr and λf is ρrf. From Equation S36, we see that the growth rate λ takes the maximum of λr and λf, i.e.,

(S156) λ=max(λr,λf).

Then, the cumulative distribution function of λ is P(λx)=xxf(x1,x2)dx1dx2, where

f(x1,x2)=(1ρrf2)122πσλrσλfexp(12(1ρrf2)[(x1μλrσλr)22ρrf(x1μλrσλr)(x2μλfσλf)+(x2μλfσλf)2]).

Thus, the probability density function of the growth rate λ is given by:

(S157) fλ(x)=122πσλre12(xμλrσλr)2[erf((xμλf)σλrρrfσλf(xμλr)σλrσλf2(1ρrf2))+1]+122πσλfe12(xμλfσλf)2[erf((xμλr)σλfρrfσλr(xμλf)σλrσλf2(1ρrf2))+1].

In Appendix 1—figure 2B, we show that Equation S157 quantitatively matches the experimental data for E. coli under the relevant conditions.

Appendix 9

Model comparison with experiments on E. coli

9.1 Flux comparison with experiments on E. coli

In Appendix 7.2, we see that the values of Jf(N) and Jr(N) are required to calculate the in vivo enzyme catalytic rates of the intermediate nodes. Here, we use Jacetate and JCO2,r to represent the stoichiometric fluxes of acetate from the fermentation pathway and CO2 from the respiration pathway, respectively. Combined with the stoichiometric coefficients of both pathways, we have:

(S158) {Jacetate=Jf,JCO2,r=3Jr.

By further combining with Equations S16-S17, we get:

(S159) {Jf(N)=JacetatemcarbonMcarbon[iri/NEPicarbon]1,Jr(N)=13JCO2,rmcarbonMcarbon[iri/NEPicarbon]1.

In fact, the values of Jacetate and JCO2,r 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/OD600/h (Basan et al., 2015). Thus, we define Jacetate(M) and JCO2,r(M) as the fluxes of Jacetate and JCO2,r (per biomass) in the unit of mM/OD600/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 OD600 roughly corresponds to 6×108 cells/mL (Stevenson et al., 2016), the average mass of a cell is 1 pg (Milo and Phillips, 2015), the biomass percentage of the cell weight is 30% (Neidhardt et al., 1990), the molar mass of carbon is 12 g (Nelson and Cox, 2008), rcarbon=0.48 (Neidhardt et al., 1990) and rprotein=0.55 (Neidhardt et al., 1990). Combined with the values of ri (see Appendix 2.2) and NEPicarbon, where NEPa1carbon=6, NEPa2carbon=3, NEPbcarbon=3, NEPccarbon=5, and NEPdcarbon=4 (Nelson and Cox, 2008), we have:

(S160) {Jf(N)Jacetate(M)/2,Jr(N)JCO2,r(M)/6.

From Equation S18, we obtain the values of ηi for each precursor pool: ηa1=0.15, ηa2=0.30,ηb=0.35, ηc=0.09, and η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 9.2.

9.2 Model parameter settings using experimental data of E. coli

We have collected biochemical data for E. coli, as shown in Appendix 1—table 1 and Appendix 1—table 2, to set the model parameters. This includes the molecular weight (MW) and in vitro kcat 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 kcat from calculations using data from four cultures (see Appendix 1—table 2). Here, we prioritize the use of in vivo kcat wherever applicable unless there is a gap in the in vivo data (see Appendix 1—table 1).

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

(S161) 1/κ3=i=131/κ3(i).

In fact, Equation 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 Figure 1B, we have the values for parameters κi (i=1, …, 6), and then εr/f(κA(C))=122(h1). Evidently, εr(κglucose(ST))<εf(κglucose(ST)), εr(κlactose(ST))<εf(κlactose(ST)), and thus εr(κAmax)<εf(κAmax). For pyruvate, we have εr/f(py)(κpy(C))=εr/f(κA(C))=122(h1) (see Equations S43 and S101), and it is easy to check that εr(κpy(ST))<εf(κpy(ST)).

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 value of κt is obtainable from experiments: the translation 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 2.1), hence κt=1/610(s1). However, there are insufficient data to determine the values of κi (i=a1, a2, b, c, d) for the metabolites between the entry point metabolites shown in Figure 1A to the precursor pools. These processes involves many steps, so these values are expected to be quite large. Here, we combine the contributions of κt and κi (i=a1, a2, b, c, d) by defining a composite parameter:

(S162) Ω1/κt+ia1,a2,b,c,dηi/κi.

We proceed to estimate the values of Ω and φ using experimental data (Basan et al., 2015) for wild-type strains on the Jacetate(M)-λ relation (Figure 1C), and then all the remaining model parameters are set accordingly.

For the case of w0=0, where all kcat values follow a Gaussian distribution with an extrinsic noise of 25% CV (which is the general setting we use unless otherwise specified), we have φ=10.8 and Ω=1345 (s). Accordingly, we obtain ηE=14.78, μλC=0.92 (h1), and σλC=0.12μλC, 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 ιw0=0(2μm Cm)=1.15, ιw0=0(4μm Cm)=2.33, and ιw0=0(8μm Cm)=6.25, where the superscript stands for the concentration of Cm, and the subscript represents the choice of w0.

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 κ9, which involves many steps in the metabolic network and thus can be considerably large. Here we define another composite parameter, ΩGg(ηb+ηc)/κ8+ηa1/κ9, and estimate its value as ΩGg=690 (s) from growth rate data for E. coli measured under the relevant nutrient conditions (Basan et al., 2015), where the subscript ‘Gg’ stands for glucogenesis. Then, μλC(py)=0.67 (h1), and σλC(py)=0.10μλC(py), 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), we have φ21AA=14.2. Comparing Equation S32 with Equation S112, the parameter Ω should change to Ω21AA1/κt+ηa1/κa1+ia2,b,c,dηi/κi(21AA). 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 μλC(21AA)=1.13 (h1), and σλC(21AA)=0.12μλC(21AA).

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 Equation S18, we have φ7AA=11.6. For the value of Ω7AA, evidently, Ω21AA<Ω7AA<Ω, and we estimate Ω7AA=1215 (s) from growth rate data for E. coli measured under the relevant culture media (Basan et al., 2015). Then, μλC(7AA)=0.98 (h1), and σλC(7AA)=0.12μλC(7AA).

For the case of w0=2.5 (h1), we have φ=8.3, and thus ηE=12.28, while other parameters such as Ω, μλC and σλC remain the same as for w0=0. Nevertheless, the values for ι under translation inhibition by Cm are influenced by the choice of w0, where the values of ι change to ιw0=2.5(2μm Cm)=1.05, ιw0=2.5(4μm Cm)=2.00, and ιw0=2.5(8μm Cm)=5.40.

From Appendix 8.1–8.2, combined with Equation S114, the distributions of λr(21AA) and λf(21AA) can be approximated by Gaussian distributions:

(S163) {λr(21AA)N(μλr(21AA),σλr(21AA)2),λf(21AA)N(μλf(21AA),σλf(21AA)2),

where μλr(21AA) and μλf(21AA) stand for the mean values, while σλr(21AA) and σλf(21AA) represent the standard deviations. For the case of glucose mixed with 21AA (labeled as ‘Glucose + 21AA’), the distribution of the growth rate λglucose(21AA) follows Equation S157. With Ω21AA=1000 (s), we have μλglucose, r(21AA)=1.34 (h1) , μλglucose,f(21AA)=1.46(h1) (both definitions follow Equation S163), and ρrf1.0 (obtained from 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. Thus, the cell growth rate (defined as λsuccinate(21AA)) would take the value of the respiration one and follows a Gaussian distribution:

(S164) λsuccinate(21AA)N(μλsuccinate(21AA),σλsuccinate(21AA)2).

For the case where acetate is the sole carbon source, the cells exclusively use the respiration pathway, and the growth rate (defined as λacetate) follows a Gaussian distribution:

(S165) λacetateN(μλacetate,σλacetate2).

Using the measured growth rate data (Wallden et al., 2016), we estimate μλsuccinate(21AA)=0.67(h1) and μλacetate=0.253(h1). To illustrate the distribution of growth rates λglucose(21AA), λsuccinate(21AA) and λacetate shown in Appendix 1—figure 2B, if no other source of noise existed, extrinsic noise with a CV of 40% would be required for each kcat value. Then, σλglucose, r(21AA)0.21μλglucose, r(21AA), σλglucose, f(21AA)0.23μλglucose, f(21AA), σλsuccinate(21AA)=0.22μλsuccinate(21AA), and σλacetate=0.22μλacetate. Allowing for the possibility that intrinsic noise may also play a non-negligible role in the observed single-cell growth rate (which is not a long-term average), we still use extrinsic noise with a CV of 25% for the model results of E. coli, except for those shown in Appendix 1—figure 2B.

Appendix 10

Explanation of the Crabtree effect in yeast and the Warburg effect in tumors

Our model, along with the analysis presented in Appendix 3, can be extended with modifications to explain the Crabtree effect in yeast and the Warburg effect in tumors. In both cases, the optimization objective remains maximizing the cell growth rate. Consequently, yeast and tumor cells use the most efficient pathway for ATP production at the single-cell level.

For model applications in yeast or tumor cell metabolism, the fermentation flux shifts from acetate secretion to ethanol and lactate secretion, respectively (see Appendix 1—figure 5A and B). The respiration and biomass generation pathways remain largely similar to those of E. coli, except that the biochemical reactions within the TCA cycle and respiratory chain occur in the mitochondria (see Appendix 1—figure 5C and D). This leads to an increased enzyme cost for the respiration pathway due to energy currency exchanges between NADH or FADH2 and ATP in the mitochondria. The coarse-grained models for Group A carbon source utilization in yeast and mammalian cells are shown in Appendix 1—figure 5E and F, where M3 represents pyruvate. In yeast and mammalian cells, the stoichiometric coefficients for ATP production (i.e. βi) are identical to each other but differ from those of E. coli (see Figure 1B and Appendix 1—figure 5C and D), with β1=5, β2=1, β3=5, β4=7.5, β6=2.5, and βa1=5 (Nelson and Cox, 2008). Hence, the stoichiometric coefficients of ATP production per glucose in each pathway are βr(A)=32 and βf(A)=2, respectively, where βr(A)=β1+2(β2+β3+β4) and βf(A)=β1+2(β2+β6).

The impact of maintenance energy in yeast and tumor cells is significantly higher than that in E. coli (Locasale and Cantley, 2010). Therefore, Equation S25 changes to (see Equation S59):

(S166) JE=rEJBM+w0Mcarbonm0,

where w0 is the aforementioned maintenance energy coefficient. Thus, we have (see Equation S60):

(S167) JE(N)=ηEλ+w0.

To account for the protein cost of energy currency exchanges in the mitochondria, we introduce ϕMT and κMT to represent the proteomic mass fraction of the enzymes and the effective substrate quality of related metabolites in the mitochondria, respectively. Note that the energy currency exchanges between NADH or FADH2 and ATP only occur during respiration, as there is no net NADH or FADH2 generation during fermentation (see Appendix 1—figure 5C and D). Combined with Equation S167, Equation S25 changes to:

(S168) {ϕAκA=ϕ1κ1+ϕa1κa1,2ϕ1κ1=ϕ2κ2+ϕ5κ5+ϕa2κa2,ϕ2κ2=ϕ3κ3+ϕ6κ6+ϕbκb,ϕ5κ5+ϕ4κ4=ϕ3κ3+ϕdκd,ϕ3κ3=ϕ4κ4+ϕcκc,ϕa1κa1=ηa1λ,ϕa2κa2=ηa2λ,ϕbκb=ηbλ,ϕcκc=ηcλ,ϕdκd=ηdλ,β1ϕ1κ1+β2ϕ2κ2+β3ϕ3κ3+β4ϕ4κ4+β6ϕ6κ6+βa1ϕa1κa1=JE(N),JE(N)=ηEλ+w0,λ=ϕRκt,Jr(N)=ϕ4κ4=ϕMTκMT,Jf(N)=ϕ6κ6,ϕR+ϕA+ϕ1+ϕ2+ϕ3+ϕ4+ϕ5+ϕ6+ϕMT+ϕa1+ϕa2+ϕb+ϕc+ϕd=ϕmax.

Here, Equation S28 still holds, and we have:

(S169) {Jr(E)+Jf(E)=φλ+w0,Jr(E)εr+Jf(E)εf=ϕmaxψλ,

where Jr(E) and Jf(E) follow Equation S30, and ψ and φ satisfy Equation S32 and S33, respectively. The expression for εf follows Equation S31. However, the expression for εr differs from that in Equation S31. For yeast and mammalian cells, we have:

(S170) {εr=βr(A)1/κA+1/κ1+2/κ2+2/κ3+2/κ4+2/κMT,εf=βf(A)1/κA+1/κ1+2/κ2+2/κ6.

At the single-cell level, from Equation S169, and similar to Equation S61-S63, if εr>εf, the optimal growth strategy is:

(S171) {Jf(E)=0,Jr(E)=φλ+w0,εr>εf,

while if εf>εr, the optimal growth strategy is:

(S172) {Jf(E)=φλ+w0,Jr(E)=0.εr<εf.

In both cases, the growth rate λ reaches its maximum value for a given nutrient condition with fixed κA:

(S173) λ(κA)={ϕmaxw0/εr(κA)φ/εr(κA)+ψ(κA)εr(κA)>εf(κA),ϕmaxw0/εf(κA)φ/εf(κA)+ψ(κA)εr(κA)<εf(κA).

From Equation S170, when κA is very small such that κA0, it is evident that for yeast and mammalian cells, we still have:

(S174) {εr(κA0)βr(A)κA,εf(κA0)βf(A)κA.

Thus,

(S175) εr(κA0)>εf(κA0),

since βr(A)βf(A) still holds. Then, as long as εr(κAmax)<εf(κAmax), there exists a critical switching point for κA (denoted as κA(C); see Equation S41), below which respiration is more efficient, while above κA(C), fermentation becomes more efficient in ATP production per proteome. Combined with Equation S170, we have:

(S176) κA(C)=βr(A)βf(A)βf(A)(1/κ1+2/κ2+2/κ3+2/κ4+2/κMT)βr(A)(1/κ1+2/κ2+2/κ6).

Accordingly, we obtain the expressions for εr(κA(C)), εf(κA(C)) and λC (i.e. λ(κA(C))):

(S177) {εr(κA(C))=εf(κA(C))=βr(A)βf(A)2(1/κ3+1/κ4+1/κMT1/κ6),λC=ϕmaxw0/εr/f(κA(C))φ/εr/f(κA(C))+ψ(κA(C)),

Consequently, yeast and tumor cells would preferentially use respiration under starvation conditions (where εr>εf), yet switch to aerobic glycolysis when nutrients are abundant (where εr<εf) for optimal cell growth. This qualitatively illustrates the Crabtree effect in yeast and the Warburg effect in tumors.

At the cell population level, cell heterogeneity resulting from intrinsic and extrinsic noise causes the turnover numbers (i.e. kcat) of enzymes and the critical growth rates at the transition point (λC) to follow distributions, which we assume to be Gaussian (see Equation S45, Appendices 3.3 and 8.1). Due to the higher level of heterogeneity observed in tumor cells (Duraj et al., 2021; Shibao et al., 2018; Hanahan and Weinberg, 2011; Hensley et al., 2016) and yeast (Bagamery et al., 2020) compared to E. coli, the extent of noise—and thus the CVs of kcat and λC—in yeast and tumor cells are expected to be larger than those in E. coli. The growth rate dependence of the normalized energy fluxes is as follows:

(S178) {Jf(E)(λ)=12(φλ+w0)[erf(λμλC2σλC)+1],Jr(E)(λ)=12(φλ+w0)[1erf(λμλC2σλC)],

where μλC and σλC are the mean and standard deviation of λC, respectively, similar to the case of E. coli. Therefore, the growth rate dependence of the normalized fluxes is:

(S179) {Jf(N)(λ)=φλ+w0βf(A)[erf(λμλC2σλC)+1],Jr(N)(λ)=φλ+w0βr(A)[1erf(λμλC2σλC)].

Combined with Equation S160, Equation S179 can be compared to experimental results, although in practice, it is difficult to tune the growth rate of tumor cells in vivo in experiments.

Recently, Shen et al., 2024 reported that in many yeast and tumor cells, the measured proteome efficiencies in respiration at the cell population level are higher than the corresponding proteome efficiencies in fermentation, even though aerobic glycolysis fermentation fluxes still occur. This finding apparently contradicts prevalent explanations (Basan et al., 2015; Chen and Nielsen, 2019), which assert that overflow metabolism originates from the proteome efficiency in fermentation always being higher than in respiration.

Our model can resolve the puzzle above based on two important features: First, our model predicts that as long as ATP generation per glucose in respiration is higher than in fermentation (i.e. βr(A)>βf(A)), which definitely holds true for all organisms, the proteome efficiency in respiration is higher than that in fermentation when the nutrient quality κA is low (see Equations S37-S38 and S174-S175). Second, and importantly, due to cell heterogeneity at the population level, a subset of cells exhibiting greater proteome efficiency in fermentation compared to respiration could exist, even if the proteome efficiency at the cell population level in respiration is higher than in fermentation.

To facilitate comparison between our model and the experiments of Shen et al., 2024, we define Prf as the proportion of ATP generated from fermentation, and Δ¯ as the proteome efficiency difference between respiration and fermentation, with

(S180) PrfJf(E)Jf(E)+Jr(E),

and

(S181) Δ¯1/εr1/εf.

At the cell population level, εr, εf,1/εr, and 1/εf roughly follow Gaussian distributions (see Appendix 8.1 and Equation S170), with

(S182) {εrN(μεr,σεr2),εfN(μεf,σεf2),1/εrN(μ1/εr,σ1/εr2),1/εfN(μ1/εf,σ1/εf2).

Here, σεr, σεf,σ1/εr, σ1/εf, and μεr, μεf,μ1/εr, μ1/εf are the standard deviations and mean values of εr, εf, 1/εr and 1/εf, respectively. Thus,

(S183) { μεr= εr ,μεf= εf ,

where the angle bracket ‘’ represents the average over the cell population, and εr and εf are the population-averaged values of εr and εf, respectively, which are both measurable in experiments. From the derivations shown in Appendix 8.1, we approximately have

(S184) {μ1/εr=1/μεr=1/εr,μ1/εf=1/μεf=1/εf.

Here, we use χεr, χεf, χ1/εr and χ1/εf to represent the CVs of εr, εf, 1/εr and 1/εf, respectively, with

(S185) {χεr=σεr/μεr,χεf=σεf/μεf,χ1/εr=σ1/εr/μ1/εf,χ1/εf=σ1/εf/μ1/εf

Similar to Equation S151, the CVs of 1/εr and 1/εf are roughly equal to those of εr and εf, respectively. Thus,

(S186) {χ1/εrχεr,χ1/εfχεf.

Combining Equations S181 and S182, and using the properties of Gaussian distributions, Δ¯ follows a Gaussian distribution:

(S187) Δ¯N(μΔ¯,σΔ¯2),

where μΔ¯ and σΔ¯ are the mean and standard deviation of Δ¯, respectively. Evidently, we have

(S188) {μΔ¯=μ1/εrμ1/εf,σΔ¯2=σ1/εr2+σ1/εf2.

Then, we proceed to calculate the relation between Prf and Δ¯ using Equation S187, and hence we obtain:

(S189) Prf=0+1σΔ¯2πe12(xμΔ¯σΔ¯)2dx=12[erf(μΔ¯2σΔ¯)+1].

Combining Equations S180, S183-S185, and S188-S189, we have:

(S190) Jf(E)Jf(E)+Jr(E)=12[erf(1εr/εf2χεr2+χεf2(εr/εf)2)+1].

Note that the normalized energy fluxes Jr(E) and Jf(E) are proportional to the measured ATP fluxes generated in respiration and fermentation, respectively. Hence, Equation S190 can be directly compared to experimental data. For yeast and tumor cells, due to a higher level of heterogeneity, the CVs of εr and εf, i.e., χεr and χεf, could be significantly higher than the corresponding values in E. coli, though their exact values are unknown. Consequently, we plot theoretical results with the values of χεr and χεf chosen as 0.25, 0.40, and 0.58 to compare with the experimental data for yeast and in vivo mouse tumors (Bartman et al., 2023; Shen et al., 2024). In Figure 5A, B, we observe that the theoretical results using χεr=χεf=0.58 agree well with the experimental data (Bartman et al., 2023; Shen et al., 2024), both on a log scale and linear scale. This demonstrates that our model has the potential to quantitatively illustrate the Crabtree effect in yeast and the Warburg effect in tumors.

Appendix 11

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 noted in the footnote of Appendix 1—table 2. With this calibration, the Jacetate(M)-λ data for E. coli (Basan et al., 2015) in Appendix 1—table 2 align with the curve shown in Figure 1C, which includes experimental data for E. coli from other sources. The second calibration applies to the data shown in Figures 3F and 1C (chemostat data for E. coli). The unit in the original reference (Holms, 1996) is mmol/(dry mass)g/h. To convert this to the unit mM/OD600/h, used 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 when calibrating with other experimental results. Therefore, there is a calibration factor of 0.6, and the conversion factor changes to 0.3.

Data from the inducible strains

Some of the experimental data in the original references (Basan et al., 2015; Hui et al., 2015) were obtained using E. coli strains with titratable systems (e.g. titratable ptsG, LacY). The Jacetate(M)-λ relation of these inducible strains generally aligns with the same curve as that of wild-type E. coli (Figure 1C). Since evolutionary treatment was not applied to the inducible strains, we approximate titration perturbation as a technique that mimics culturing the strains in a less efficient Group A carbon source.

Experimental data sources

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

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

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

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

The batch culture data for yeast shown in Figure 5 were derived from the source data of the reference’s extended figure 4c-d (Shen et al., 2024). The chemostat data for yeast shown in Figure 5 were derived from the source data of the reference’s figure 3d-e (Shen et al., 2024), where glucose is the limiting nutrient. We excluded the reference’s data for I. orientalis under condition C2, where the ATP flux was abnormally small. The mouse tumor in vivo data shown in Figure 5 were derived from the source data of the reference’s figure 4e-g (Shen et al., 2024), which were originally reported by Bartman et al., 2023, the same research group as Shen et al., 2024. We did not include the cancer cell line data shown in figure 4a-c of Shen et al., 2024 since it appears that the proteomic data and flux data were obtained from two different references with inconsistent culturing conditions.

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

Data availability

All study data are included in the manuscript and supporting files. All model results were generated using analytical formulas, with the relevant formulas and parameters specified in the manuscript and appendices. Source data files have been provided for Figures 1–5 and Appendix 1—figures 2–4.

References

  1. Book
    1. Kampen NG
    (1992)
    Stochastic Processes in Physics and Chemistry
    Elsevier.
  2. Book
    1. Neidhardt FC
    2. Ingraham JL
    3. Schaechter M
    (1990)
    Physiology of the Bacterial Cell
    Sinauer Associates, Inc.
  3. Book
    1. Neidhardt FC
    (1996)
    Escherichia coli and Salmonella: Cellular and Molecular Biology
    ASM Press.
  4. Book
    1. Nelson DL
    2. Cox MM
    (2008)
    Lehninger Principles of Biochemistry
    Macmillan.

Article and author information

Author details

  1. Xin Wang

    School of Physics, Sun Yat-sen University, Guangzhou, China
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    wangxin36@mail.sysu.edu.cn
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6479-395X

Funding

National Natural Science Foundation of China (12004443)

  • Xin Wang

National Natural Science Foundation of China (12474207)

  • Xin Wang

Guangzhou Municipal Science and Technology Bureau (202102020284)

  • Xin Wang

Sun Yat-sen University (The Hundred Talents Program)

  • Xin Wang

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

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

Version history

  1. Sent for peer review:
  2. Preprint posted:
  3. Reviewed Preprint version 1:
  4. Reviewed Preprint version 2:
  5. Reviewed Preprint version 3:
  6. Version of Record published:

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.94586. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2024, Wang

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,242
    views
  • 101
    downloads
  • 0
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Xin Wang
(2025)
Overflow metabolism originates from growth optimization and cell heterogeneity
eLife 13:RP94586.
https://doi.org/10.7554/eLife.94586.4

Share this article

https://doi.org/10.7554/eLife.94586