A compact model of Escherichia coli core and biosynthetic metabolism

To assemble the iCH360 model, we started from the core metabolic reactions present in ECC and extended them with the pathways required for the biosynthesis of the main biomass building blocks, including the twenty amino acids, the four nucleotides, and both saturated and unsaturated fatty acids (Fig. 1, Table 1, and Supplementary Figures 1-4). In addition, while not performing a comprehensive review of iML1515, we applied a small number of corrections to the original reactions based on literature knowledge (Supplementary Information A.1).

The final assembled model (Figure 1) comprises $304$ compartment-specific metabolites and $323$ metabolic reactions mapped to $360$ genes, thus qualifying as a medium-scale model ranging in between ECC and iML1515 (Supplementary Figure 5). While similar in scale, our model and ECC2 present a fundamental structural difference. ECC2 was built by systematically removing reactions from its genome-scale parent (iJO1366) [hadicke_ecolicore2_2017]. Thus, its metabolic space spans the production of all compounds consumed in the iJO1336 biomass reaction. In contrast, the metabolic space of iCH360 only reaches biomass building blocks, without explicitly considering peripheral pathways, such as assembly of cell-membrane components, de novo synthesis of cofactors, and active transport of ions in the cell.

To make the model comparable to its parent model iML1515, we computed an equivalent biomass reaction, wherein all biomass requirements not present in our model are converted into an equivalent metabolic cost in terms of precursors, based on manually curated bioproduction routes (Table 2, Supplementary Information A.2).

To check how its significantly smaller size and complexity affects the solution space of our model, we first looked at the maximally achievable biomass production flux under a range of growth conditions (Supplementary Figure 6). Across most conditions considered, the model can achieve biomass fluxes comparable to its genome-scale parent. The main differences exist in anaerobic growth on fumarate, alpha-ketoglutarate (AKG) and glycerol, where our model predicts zero growth, while the GEM achieves some (albeit small) biomass production rate. In practice, fermentation in these scenarios is in fact biologically unrealistic (e.g. [boecker_enabling_2022]). The reduced model is thus a good basis for metabolic simulation frameworks with few constraints such as FBA, since the reduction procedure, at least in this test, limits the original solution space of the GEM to more physiologically relevant regions.

To investigate the metabolic capabilities of our model even further, we looked at production envelopes (projections of the model’s solution space onto a smaller set of dimensions, also known as phenotypic phase planes), describing production rates for biomass and a range of metabolites. Figure 2 shows the computed envelopes for aerobic growth on glucose. The reduced model has a production envelope comparable with iML1515, where the main apparent difference is in the production of acetate. These differences can be ascribed to the absence of a series of degradation pathways in our model. These pathways allow the genome-scale model to channel additional carbon flux towards acetyl-CoA and, as a result, to reach higher theoretical acetate production yields of up to $\approx 3$ mole acetate per mole of glucose. Since the physiological relevance of such high acetate yields, in the absence of heterologously expressed pathways, has already been questioned [hadicke_ecolicore2_2017], we chose not to include in the model those additional reactions that would allow for higher acetate yields. A similar pattern of agreement with iML1515 was found when comparing production envelopes under different growth conditions (Supplementary Figure 7).

Metabolic models often contain annotations that connect model elements to entries in biological databases, such as EcoCyc [keseler_ecocyc_2017], KEGG [kanehisa_kegg_2017] and MetaNetX [moretti_metanetxmnxref_2021]. However, even for the subset of reactions included in iCH360, the annotations present in iML1515 were incomplete and in part outdated. To fill these gaps, we extended and corrected the original annotations through a mixture of automated querying and manual curation (Figure 3A). Notably, the annotations to EcoCyc, the main knowledgebase for E. coli, are complete so that every metabolic reaction in iCH360 is uniquely mapped to a corresponding EcoCyc ID (with the exception of three reactions, for which a match on the database could only be found for a different choice of redox cofactor, namely NADH instead of NADPH). Additionally, we found $134$ deprecated annotations pointing to the MetaNetX namespace, which were consequently updated with the most up-to-date IDs.

Leveraging the complete mapping with the EcoCyc database, we parsed, assembled, and manually curated a data structure that enhances the stoichiometric model with a detailed layer of information about enzymes, polypeptides, and genes related to the network (Methods). This data structure takes the form of a weighted graph, wherein nodes represent biological entities (reactions, proteins, genes or compounds) and edges represent (potentially quantitative) functional relationships between them (Methods, Figure 3B, Supplementary Tables 1 and 2), such as catalysis, regulation, protein modification, and protein subunit composition. This graph contains the gathered information in a unified form, allowing users to perform a number of tasks related to metabolic modelling applications. For example, by explicitly mapping reactions to their catalysing enzymes rather than to the associated genes (which may additionally include protein activators or cofactors), it simplifies the definition of meaningful enzyme capacity constraints in the model. Similarly, since the graph topology implicitly defines associations between reactions and genes, boolean gene-protein-reaction (GPR) rules needed for in silico knockout studies can be generated accounting for both catalytic and non-catalytic requirements for each reaction (Supplementary Information A.3). Crucially, these GPR rules can be re-generated as needed whenever the graph is updated with new nodes or edges. Finally, we employed it to estimate the abundance of protein complexes included in the model from measured polypeptide abundances (Supplementary Information A.4) using an automated procedure (see Result section 2.5).

Through this annotation graph, $318$ metabolic reactions in the model are linked to $289$ catalysing enzymes, with more than $25\%$ of the reactions being catalysed by multiple enzymes (isozymes). Such enzymatic redundancy plays an important role, for example, when designing metabolic engineering strategies to prevent flux through a pathway. However, the different isoenzymes of a reaction need not all be equally important, and treating them as completely equivalent may generate inaccuracies in some phenotypic predictions. For example, while phosphofructokinase activity in E. coli (reaction PFK) is known to be carried by two isozymes, pfkA and pfkB, the latter is known to be responsible only for minor activity under normal physiological conditions [babul_phosphofructokinases_1978]. If these differences are not taken into account, metabolic models can overpredict redundancy in the network, for instance when assessing phenotypic predictions in response to gene knockouts [monk_iml1515_2017].

To address this issue and make the model usable in a wider range of modelling scenarios, we qualified catalytic edges in the graph as either primary or secondary catalysis (Methods, Figure 3C). A catalytic relationship between a reaction and an enzyme was annotated as secondary whenever the enzyme, according to experiments, accounts only for negligible activity for the reaction in the wild-type strain. As experimental evidence we considered in vitro and in vivo complementation studies. Through this curation process, a total of $72$ catalytic edges were functionally annotated as secondary.

To test how this functional annotation can support phenotypic predictions, we investigated the disruption of each essential reaction in the model resulting from the knockout of each one of its associated genes (Methods). We classified the outcome of each possible knockout as: (i) complete disruption, if a reaction loses all of its catalytic edges, (ii) full primary disruption, if a reaction loses all of its primary catalysis edges, but secondary ones are left, (iii) partial primary disruption, if the reaction loses some, but not all, of its primary catalysis edges, or (iv) secondary disruption, if some or all secondary edges are disrupted, but none of the primary ones. We classified the effect of in silico knockouts for aerobic growth on $9$ different carbon sources and compared the results against a large dataset of mutant fitness data, obtained via competitive fitness assays, from [price_mutant_2018].

Based on this analysis, we found that disruptions of primary edges as a result of a knockout were significantly associated with greater fitness losses than disruptions of secondary edges (Wilcoxon rank-sum test, $p<10^{-6}$, Figure 3D). In addition, primary disruptions that were only partial were associated with more contained fitness losses. Finally, minor fitness gains could be identified when mutants were associated with secondary and partial primary disruptions, but not when complete primary disruptions occurred. We did not find as significant differences between complete disruptions and full primary disruptions (Supplementary Figure 8), supporting the idea that catalytic relationships annotated as secondary are unlikely to be strong enough to stand-in for disrupted primary ones under normal physiological conditions.

To extend the scope and applicability of iCH360, we constructed a version of the model that allows for enzyme-constrained flux simulations, which we denote as EC-iCH360. We constructed EC-iCH360 in the sMOMENT [bekiaris_automatic_2020] format (Methods). The sMOMENT framework is inherently simple and generates the same solution space as more complex model formats, such as GECKO, unless reaction-specific capacity constraints are specified in the latter [sanchez_improving_2017, bekiaris_automatic_2020]. Since it requires unique reaction-enzyme mappings, we used the knowledge graph to remove all secondary catalytic relationships in the model (Section 2.4).

We first parametrized the model by defining a flux-specific enzyme cost for each reaction (in units of grams of enzyme per unit flux), using as values the estimated in vivo turnover number from [heckmann_kinetic_2020]. Using this enzyme-constrained model, we then predicted enzyme abundances and compared them to experimental enzyme abundances, estimated from a dataset of proteomic measurements for aerobic growth on eight different carbon sources [schmidt_quantitative_2016](Methods). This analysis led to generally good predictions, with root mean squared error (RMSE, computed for log₁₀-transformed enzyme abundances) ranging from $0.53$ to $0.62$ (Supplementary Figure 9). In order to assess the nature of residuals between measurements and predictions, we investigated the geometric mean of enzyme abundances across all eight conditions (Supplementary Figure 9, bottom right panel). Notably, averaging predicted and measured enzyme abundances, respectively, across conditions did not significantly improve the RMSE, indicating that the abundances of individual enzymes were systematically over- or under-predicted across conditions.

To increase the predictive ability of the model, we thus adjusted the turnover numbers by fitting them to experimental measurements through a custom heuristic (Methods, Supplementary Information A.5). By simultaneously fitting all available conditions, we ensured that our adjustment procedure is robust to condition-specific biases. Further, by introducing regularisation within our adjustment scheme, we penalised large deviations of parameters from the original dataset, increasing the robustness of the procedure to overfitting. Our adjusted parameter set shows a mean absolute deviation (computed for log₁₀-transformed turnover values) from the original parameter set of $\approx 0.22$ (Supplementary Figure 10) and results in significantly better agreement with experimental measurements across conditions (Figure 4A). Further, mean enzyme abundances across conditions are very well predicted with the adjusted parameter set (Figure 4B), implying that residuals between measurements and predictions are now to be attributed to variability in enzyme cost across conditions (which simple frameworks such as enzyme-constrained FBA cannot account for), rather than systematic over- or under-predictions. Thus, each adjusted turnover parameter can be thought of as a ”typical” apparent $k_{\rm cat}$ value, incorporating average saturation trends for an enzyme across growth conditions.

Despite the small size of iCH360, we found the explicit enumeration of its elementary flux modes (EFMs) to be intractable. This is not necessarily surprising, since the topology of a metabolic network is a much better determinant of EFM count than its sheer size. Metabolic networks can possess different types of redundancies, such as the presence of alternative pathways for the production of the same metabolite, the use of alternative cofactors for the same catalytic step in a pathway, or the presence of alternative transporters for the uptake/excretion of a compound. While knowledge about these redundancies is often valuable, including them in the model can increase the number of EFMs exponentially, hampering or even preventing EFM-based analyses.

To address this issue, we identified and removed a small set of such redundancies in iCH360, using available information from literature, whenever possible, to ensure the most physiologically relevant alternative was maintained (Supplementary Table 3). This resulted a metabolic submodel of iCH360, a model variant that we denote by iCH360red. iCH360red contains $305$ metabolic reactions ($18$ less than iCH360) and shares the same production envelopes with its parent model for a number of metabolites of interest (Supplementary Figure 12). While the number of EFMs in iCH360red is still relatively large ($\approx 13.5$ millions for aerobic growth on glucose, see Supplementary Table 4), it is not prohibitive for most types of EFM-based analysis, and their explicit enumeration does not require high-performance computing (Methods).

We used the obtained EFMs to study the possible combinations of achievable growth rates and yields in the network [wortel_metabolic_2018]. To this end, we considered growth on glucose as a scenario and computed, for each EFM from iCH360red, its yield and its achievable cell growth rate, which we estimate based on the enzyme costs defined for the enzyme-constrained model (Methods). Based on this analysis, we identified a front of Pareto-optimal EFMs, along which any increase in the growth rate will necessarily lead to a reduction in yield (Figure 5). Along the Pareto front, we observe a transition from a purely respiratory mode at maximum yield (Supplementary Figure 13) to a mixed respiratory-fermentative mode at maximum growth (Supplementary Figure 14). Quantitatively, the extent of this tradeoff was rather modest: the EFM with maximal yield reaches almost the maximal growth rate, so only minor gains in growth rate can be achieved by using other, fermentative modes along the Pareto front. However, it is worth noting that this analysis was performed using a simple capacity-based enzyme cost function, which summarises different determinants of enzyme cost such as driving force and substrate saturation [noor_protein_2016] in one single value, assumed to be constant. Hence, repeating the analysis with a more complete enzyme-cost function [noor_protein_2016] could help elucidate the nature of this tradeoff. To demonstrate that the shape and size of the Pareto front depend strongly on growth conditions, we also simulated an environment with very low oxygen levels. We implemented this by increasing the flux-specific enzyme cost of the oxygen-dependent reaction in the respiratory chain, equivalent to assuming a lower enzyme efficiency due to a lower oxygen level (Methods). Results (Supplementary Figure 15) show a much broader front of Pareto optimal EFMs, indicating that the nature of the observed tradeoff is indeed condition-dependent. Notably, a similar dependence of the Pareto front on extracellular oxygen availability has been previously observed in a small-scale model of E. coli core metabolism [wortel_metabolic_2018].

In order to study the effect of external conditions on optimal metabolic strategies in more detail, we used another framework that does not require an enumeration of EFMs and that allows for additional flux bounds, for example, a lower bound on the ATP consumption for cell maintenance. The saturation FBA (satFBA) framework [muller_resource_2015] is a variant of enzyme-constrained FBA, wherein a fixed enzyme cost per flux is assumed for all metabolic reactions in a model, except for the substrate transporter, for which a complete kinetic law is used (Methods). Since the external substrate concentration is a simple parameter, screening this concentration is equivalent to screening the values of the transporter efficiency. Here we used satFBA to simulate how the growth-maximising solution of the network varies in response to changing extracellular glucose concentration. By solving the satFBA problem for a range of glucose concentrations, we predicted the dependence of the cell’s growth rate on substrate concentration, resulting in the typically observed Monod curve (Figure 6A). If the problem does not contain any further flux bounds (so that the magnitude of fluxes at the optimum is solely limited by the maximum enzyme availability) the solution of satFBA problems will be an elementary flux mode [muller_resource_2015]. Hence, in this case we can use satFBA to explore how a cell should switch across elementary modes as a function of the growth environment.

At low glucose concentrations, the glucose transporter operates at low saturation and glucose uptake is enzymatically expensive, leading to a high-yield, purely respiratory metabolic mode at the optimum (Figure 6B and C). As substrate availability is increased, the cost of substrate uptake decreases and higher growth rates are achieved by switching to lower yield, acetate-secreting modes, in line with typical observations of overflow metabolism in E. coli [basan_overflow_2015, chen_energy_2019]. However, the satFBA formalism can also be used with additional flux bounds, thus going beyond EFM-based analysis. For example, if a positive lower bound on ATP hydrolysis is added as a maintenance requirement, optimal solutions to the satFBA problem will no longer be elementary modes, and the yield of the optimal solution no longer follows a piecewise constant profile (Supplementary Figure 16)

Living systems operate out of thermodynamic equilibrium, and thermodynamics places strong constraints on the operation of metabolic systems. In any metabolic state, the flux directions must follow the signs of thermodynamic forces, which depend on metabolite concentrations and equilibrium constants. To provide these constants as parts of our model, we used the component contribution framework [noor_consistent_2013] to estimate the standard Gibbs free energy ($\Delta_{r}G^{\prime\circ}$) of each reaction (Methods). These estimates account for compartment-specific chemical environment parameters, such as pH, pMg, and ionic strength, and were corrected to account for protons and charge translocation in multi-compartment reactions.

The resulting parameter set covers the vast majority of the model (over 97 % of metabolic reactions, with the remaining being not covered by the component contribution database used here) and accounts for the uncertainty in the estimates through a multivariate covariance matrix. Accounting for correlations between the different $\Delta_{r}G^{\prime\circ}$ values becomes important when imposing thermodynamic constraints on the model. For example, the fatty acid biosynthesis subsystem in the model consists of repeated elongation cycles, where a short sequence of chemical transformations is repeatedly performed on a growing carbon chain. As a result, even if the $\Delta_{r}G^{\prime\circ}$ for each reaction in the pathway is known with some uncertainty, this uncertainty is tightly correlated across the reactions, which constrains the set of achievable thermodynamic states in the network.

Using this set of thermodynamic constants, we first tested whether some typical flux distributions obtained from the model are thermodynamically feasible. To this end, we considered flux distributions generated by parsimonious Flux Balance Analysis (pFBA) across $12$ growth conditions and computed their max-min driving force (MDF) [noor_pathway_2014], accounting for uncertainty in the estimates (Methods). We found a positive MDF for each of the flux distributions, indicating that all pFBA solutions tested are thermodynamically feasible. Notably, we found that the computed MDF values cluster very clearly in three groups, corresponding to aerobic growth on glycolytic substrates, aerobic growth on gluconeogenic substrates, and anaerobic conditions (Figure 7A).

Having confirmed that our reference FBA-derived flux distributions are thermodynamically feasible, we then employed an alternative flux prediction method that ensures thermodynamically feasibilty by construction. The probabilistic metabolic optimisation (PMO) framework [gollub_probabilistic_2021] uses a mixed-integer quadratic programming approach to compute a set of fluxes, metabolite concentrations, and reaction driving forces which is probabilistically most in agreement with experimentally measured metabolite concentrations (Methods). Using this framework, we computed a maximum-likelihood thermodynamic state for the model, simulating aerobic growth on glucose. In the thermodynamic state computed by PMO, we found that all metabolites lie within physiologically reasonable concentrations (1 $\mu$M – 1 mM). In addition, all ”anomalous concentrations” identified by the framework (metabolite concentrations lying more than one standard deviation away from the mean experimental value) had been identified and explained previously [gollub_probabilistic_2021], and can potentially be addressed by lumping together those reactions for which substrate channelling is known to happen in literature.

We used the PMO-derived thermodynamic state to identify candidate bottlenecks, in terms of enzyme demand, across the network. To this end, we first computed the flux-force efficacy of each reaction in the thermodynamic state [noor_pathway_2014] (Methods). The flux-force efficacy is a unitless quantity, ranging from $0$ to $1$, denoting the ratio between net flux (forward minus backward flux) and total flux (forward plus backward flux) of a reaction. Reactions operating at low flux-force efficacy have a lower net flux due to two reasons: (1) the forward flux is lower in absolute terms, and (2) the higher backward flux is counter-productive and subtracts from the forward flux. Therefore, to achieve a given required net flux, the cell has to invest more resources in maintaining a higher enzyme level. To identify potential thermodynamic bottlenecks leading to high enzyme costs in specific reactions, we thus screened all reactions for their predicted flux-force efficacy and flux (Figure 7B) and identified those predicted to operate at low efficacy, while carrying significant flux (Figure 7B, labelled points). To assess the extent to which the thermodynamic operating states predicted by PTA are predictive of enzyme investment, we split reactions in two groups based on whether their predicted flux-force efficacy sits above or below $50\%$, and compared the distributions of measured enzyme abundances between the two groups (Methods, Supplementary Figure 17). While there are many determinants of enzyme abundance beyond thermodynamics, including enzyme turnover, affinity, or regulation, we observed a significant difference between the two distributions ($p<0.001$, two-sided Wilcox sum-rank test), with the low efficacy having approximately a 3-fold higher median enzyme abundance than the high efficacy group.

All the relevant pathways included in the model were assembled and curated based on information available in the EcoCyc [keseler_ecocyc_2017] and KEGG [kanehisa_kegg_2017] databases. The respective reactions were then extracted from iML1515 and parsed into a new model. In order to compute an equivalent biomass reaction, we first collected all the pathways required for the production of any component present in the iML1515 biomass reaction, but not in our model. These additional pathways (available in the repository supporting this manuscript) were manually curated based on available literature and database annotation to ensure they represent the biologically most relevant bioproduction route for each biomass component. By adding these pathways to iCH360, we obtained an extended model that was able to predict growth rates directly through the original iML1515 biomass reaction. The equivalent biomass function was then computed based on a reference flux distribution computed on this extended model, as explained in Supplementary Information A.2. Model assembly, manipulation and validation were performed using the COBRA Toolbox [ebrahim_cobrapy_2013]. Extension of database cross-annotation for the model reactions was performed through a mixture of automated database query and manual curation.

iCH360 can be visualised through a series of custom-built maps using the metabolic visualisation tool Escher [king_escher_2015] (see Figure 1). There are three main ways to visualise the model or solutions thereof. First, a complete map of the model, including all of its reactions, can be used (Figure 1). In order to provide a more compact representation of the network, a compressed second variant of the same map was constructed (Supplementary Figure 18). Here, long biosynthetic linear pathways were lumped into single pseudo-reactions, which only show the net production or consumption of metabolites by the pathway while omitting intermediates. Finally, individual maps for each of the main subsystems in the model are provided.

Information about the enzymes and proteins behind the stoichiometric model was collected in a knowledge graph. To build this graph, all available data on reaction-protein association and subunit composition were retrieved by automatic querying of the BioCyc database through its REST-based data retrieval API (https://biocyc.org/web-services.shtml). This information was then extended and curated based on a comparison with existing iML1515 GPR annotations. The resulting data were used to generate a directed graph wherein nodes represent biological entities (such as reactions, proteins, genes, and metabolites) and edges represent functional dependencies across them, including catalysis, subunit composition, post-translational modifications, and others. A complete list of node and edge types in the graph are provided in Supplementary Tables 1 and 2, respectively. All polypeptide nodes were annotated with their molecular mass (parsed from EcoCyc), enabling recursive computation of the molecular mass of any protein node in the graph (see Supplementary Information A.3). All manipulation and analyses of the graph data structure were performed using the NetworkX Python package [hagberg_exploring_2008], and the final data structure provided to the user in both Cytoscape and GML formats.

Catalytic edges in the graph, i.e. edges connecting reaction and enzyme nodes, were manually annotated as either primary or secondary, based on available evidence for the activity of each enzyme in the model with respect to its associated reactions. More specifically, a catalytic edge between a reaction-enzyme pair was labelled as secondary whenever the enzyme had been shown in literature to account for only minor catalytic activity for a reaction when compared to another isozyme. In this case, the references to the relevant literature were included as metadata for that edge. Whenever sufficient information was not available, all isozyme edges for a reaction were conservatively treated as primary.

For the catalytic disruption test, we identified condition-specific essential reactions in the following way: in the model, for each condition considered, we determined the reactions whose knockout led to inability to produce biomass. A small number of known false-positives, which are essential in iCH360 only due to its lack of certain reactions or pathways, were excluded from the analysis. For each essential reaction, we considered all associated genes and investigated the removal of each of them, individually, from the graph. The effect of the simulated knockout was propagated across the graph by removing all nodes from the graph for which the gene is required according to the nodes’ boolean GPR rules (see Supplementary Information A.3). Finally, the outcome of each simulated knockout was catalogued based on the reaction-level disruption it caused (see main text). Whenever multiple reactions were disrupted by the knockout of a gene, the strongest disruption among them was assigned to the knockout, in the following precedence order: complete disruption, full primary disruption, partial primary disruption, secondary disruption. The analysis was repeated for a total of $9$ growth conditions, and each condition-gene pair, labelled with the assigned disruption type, was compared with experimentally measured mutant relative fitness (averaged across available replicates for that condition-gene combination), using the dataset from [price_mutant_2018].

In the sMOMENT formalism, positive and negative fluxes in each reaction are formally described, respectively, as positive fluxes in separate “forward” and “backward” versions of the reaction. To construct the enzyme-constrained model in sMOMENT format, reversible reactions in the model were therefore duplicated in the model to separately represent fluxes in forward or backward direction. Direction-specific turnover number estimates were parsed from [heckmann_kinetic_2020] and a default value of $65~{}\rm{s}^{-1}$ was used for transporters as in [heckmann_machine_2018]. To account for the fact that [heckmann_kinetic_2020] reports values as turnover numbers per polypeptide, the values were multiplied by the number of polypeptide subunits in each enzyme. For each reaction-enzyme pair, an enzyme cost per unit flux (in units of ${\rm g\cdot h\cdot mmol^{-1}}$) was then defined as:

where $a_{i}$ is the enzyme cost of reaction-enzyme pair $i$, $k_{\text{cat},i}$ is the turnover rate estimate for the pair (here, in units of $\rm h^{-1}$), $M_{i}$ is the molecular mass of the enzyme involved (in ${\rm kDa}$), and $\sigma$ is a unitless condition-specific scaling factor (typically interpreted as an average enzyme saturation value). Then, a unique enzyme was assigned to each reaction. To this end, secondary catalytic relationships were first discarded and, for reactions with multiple annotated primary isoenzymes, the enzyme with the highest measured abundance in the integrated PAX Database [huang_paxdb_2023] was heuristically chosen. Based on these costs per unit flux, an enzyme capacity constraint was introduced as:

where $v_{i}$ denotes the flux in reaction $i$ (in ${\rm mmol~{}h^{-1}~{}gDW^{-1}}$) and $E_{\text{tot}}$ is a parameter denoting the total amount of enzyme mass (in g/gDW) that can be allocated to the flux mode. The constraint is enforced by augmenting the stoichiometric matrix of the model with an additional enzyme supply pseudoreaction (upper-bounded by $E_{\text{tot}}$) and an enzyme pool pseudometabolite consumed in each reaction with stoichiometry $e_{i}$ [bekiaris_automatic_2020].

Turnover numbers from [heckmann_kinetic_2020] were adjusted based on the nonlinear programming (NLP) formulation detailed in Supplementary Information A.5). The reference flux distributions required by the procedure were obtained using the original (unadjusted) set of turnover numbers and the bounds on allowable adjustments ($\mathbf{l}$ and $\mathbf{u}$ in Eq. 29) set to $\pm 2$ (corresponding to a maximal 100-fold increase or reduction of each parameter from the original value). The linear program used to obtain reference flux distributions was formulated and solved with GUROBI [gurobi], while the nonlinear program used for turnover adjustment was formulated and solved with the open-source optimisation package CasADi [Andersson2019]. To investigate the effect of the ridge regularisation hyperparameter ($\rho$ in equation 29), we solved the adjustment problem for a broad range of values of this parameter and, each time, we computed the RMSE between measurements and predictions of enzyme abundance, as well as the mean absolute deviation of between original and adjusted turnover values, both computed for log₁₀-transformed data (Supplementary Figure 11). Based on this information, a value of $\rho=1$, after which any further decrease in the amount of regularisation results in marginal reduction of the RMSE, was chosen to compute the final set of adjusted turnover numbers. Upon reparametrisation into apparent turnover numbers, (see Supplementary Information A.5.3), this set of adjusted parameters was used to parametrise the enzyme-constrained model variant EC-iCH360.

To validate enzyme allocation predictions against experimental values, we first retrieved measured polypeptide abundances for each growth condition from [schmidt_quantitative_2016] and imputed missing values using, whenever available, abundance values from the PAX Database [huang_paxdb_2023]. Then we estimated enzyme counts across conditions from polypeptide counts using non-negative least-squares estimation (Supplementary Information A.4). Next, we converted them into mass abundances (in units of ${\rm g~{}gDW^{-1}}$) based on the molecular mass of each complex (see Supplementary Information A.3) and assuming a cell mass of $2.8\times 10^{-13}\rm{g}$ (BIONUMBER ID 103904 [milo_bionumbers--database_2010]). Enzyme allocation predictions for each condition were then computed via EC-iCH360 by fixing the growth rate to the experimentally measured one from [schmidt_quantitative_2016] and minimising the total enzyme cost. To compute predictions using the turnover number estimates from [heckmann_kinetic_2020], the average saturation coefficient $\sigma$ was estimated, for each condition as:

where $\mathcal{M}$ denotes the index set of model enzymes for which measurements are available and $\bar{E_{\rm{tot}}}$ is the total measured model enzyme abundance for that condition. This value of $\sigma$ was then used to scale predicted enzyme abundances before comparing them to the experimentally measured ones. To compute predictions with the adjusted parameter set, the scaling factors for each condition were obtained as part of the fitting procedure (see Supplementary Information A.5). All reported root mean squared errors were computed on log₁₀ transformed enzyme abundances, excluding zero-predictions from the computation.

Elementary Flux Modes (EFMs) enumeriation for the submodel iCH360red was performed, for each growth condition, using EFMtools [terzer_large-scale_2008]. Filtered modes (Supplementary Table 4) were defined as those supporting nonzero biomass flux and, for aerobic modes, nonzero oxygen uptake. In addition, in the aerobic case, filtered modes exclude those carrying flux in either of three reactions – Pyruvate-Formate Lyase (PFL), Fumarate Reductase (FRD2), and the menaquinone-dependent Dihydroorotate dehydrogenase (DHORD5) – known to be only physiologically active under anaerobic conditions [zhang_inactivation_2001, cecchini_succinate_2002, andrews_anaerobic_1977]

To analyse trade-offs between growth rate and yield, we computed the yield of each EFM as the ratio of its biomass flux $v_{\rm BM}$ to its glucose uptake flux. In FBA models, fluxes are given in units of mmol/gDW/h, except for the biomass flux, whose conventional unit is $\rm{h}^{-1}$ because the biomass flux in these models is scaled to describe the biomass production rate per amount of biomass present. This means that our yields are given in units of gDW/mmol. To determine the cell growth rate allowed by a flux distribution, we first computed its absolute enzyme cost

where $c_{\rm enz}$ is the enzyme cost of the flux distribution (an enzyme mass, measured in g/gDW), $v_{i}$ is the flux of reaction $i$, and $a_{i}$ is the enzyme cost per unit flux in reaction $i$, computed as per equation 1 using the set of adjusted turnover numbers. Clearly, both $v_{\rm BM}$ and $c_{\rm enz}$ depend on the particular choice of scaling of the mode (though their ratio, $v_{\rm BM}/c_{\rm enz}$, doesn’t). In order to obtain an estimate for the growth rate that is scale-invariant, we thus rescaled all modes to the same enzyme cost $f$ (in g/gDW), and looked at the resulting flux through the biomass reaction. More formally, the achievable growth rate $\mu$ (in $\rm h^{-1}$) for an (arbitrarily scaled) EFM with biomass flux $v_{\rm BM}$ and enzyme cost $c_{\rm enz}$ was computed as:

For simplicity, we approximate $f$ by a constant value of $f=0.285$ g/gDW, which we obtained by taking the minimum enzyme investment required by the enzyme-constrained model to support the experimentally measured growth rate reported in the proteomic dataset [schmidt_quantitative_2016] for the condition of interest in this analysis (aerobic growth on glucose).

To simulate low-oxygen conditions, the cost per unit flux of all oxygen-consuming reactions (CYTBO3_4pp, CYTBDpp, CYTBD2pp) was increased by a $1000$-fold, mimicking the physiological state in which these reaction operate at low substrate saturation.

Saturation FBA calculations were performed by optimising biomass production in the enzyme-constrained model, setting the saturation coefficient ($\sigma$ in Eq. (1) to the value fitted as part of the turnover number adjusting procedure (Section 4.7) for all reactions except the glucose transporter (GLCptspp). The saturation of the glucose transporter, $\sigma_{\rm up}$, was computed as a function of external glucose concentration assuming irreversible Michaelis Mentens kinetics, so that:

where $\rm[Glc]$ is the external glucose concentration (in mM) and a value of $0.116$ mM was used for the Michaelis constant $K_{\rm m}$ [wortel_metabolic_2018]

Estimates of the free energies of reactions, and their uncertainties, were obtained using the component contribution framework previously described in [noor_consistent_2013]. Several reactions in the model involve protein side groups, such as the acyl-carrying protein (ACP), or cofactors, such as glutaredoxin, for which a decomposition in terms of chemical groups is not available. As a result, thermodynamic constants for these reactions cannot be directly estimated through database-based implementations of the component contribution method, such as eQuilibrator [beber_equilibrator_2022]. However, since these non-decomposable protein groups are conserved in all reactions within our model, their net contribution to the reaction thermodynamics is, at least from a group contribution perspective, null. Hence, we first augmented the group incidence matrix of the eQuilibrator database (see Section S1.1 in [beber_equilibrator_2022]) by treating protein groups as non-decomposable units and manually curated the group decomposition of the metabolites containing them. The mean and covariance of the Gibbs standard free energy of reaction ($\Delta_{r}G^{\circ}$) estimates for all reactions in the model were then computed using this augmented matrix as described in [beber_equilibrator_2022].

In order to compute transformed Gibbs free energies of reactions ($\Delta_{r}G^{\prime\circ}$), for which an exact chemical definition of each metabolite is required, protein side groups were replaced by an appropriate chemical moiety best approximating the chemical environment of the metabolite. Specifically, ACP- groups were replaced by a phosphopantetheine group, the natural prosthetic group of acyl-carrier proteins, with a methyl group at the attachment site to the protein scaffold. Similarly, glutaredoxin was replaced by its Cys-Pro-Tyr-Cys active site, with the two cysteines being either free (for the reduced form of the cofactor) or linked by a disulfide bridge (for the oxidised state).

Corrections for reactions occurring across different compartments were computed as described in [beber_equilibrator_2022], using compartment-specific pH, pMg, ionic strength and potentials from [gollub_enkie_2023].

The max-min driving force (MDF) of each reference flux distribution was computed by extending the original formulation described in [noor_pathway_2014] to account for correlated uncertainty in the estimates of the thermodynamic constants. Let $\Delta G_{r}^{\prime\circ}\in\mathbb{R}^{N}$ be a random vector representing the (uncertain) standard Gibbs free energy of reaction for the reactions in the network. Importantly, this vector includes only balanced metabolic reactions and excludes pseudoreactions such as exchange reactions. To describe our uncertain knowledge, we assume that the vector $\Delta G_{r}^{\prime\circ}$ follows a multivariate normal distribution:

where $\bar{\Delta G_{r}^{\prime\circ}}$ and $\Sigma$ are the mean and covariance of the estimates obtained through the component contribution methods. This random vector can equivalently be expressed as:

where $\mathbf{z}$ is a standard normal random vector in $\mathbb{R}^{q}$, with $q=\text{rank}(\Sigma)$, and $\mathbf{Q}\in\mathbb{R}^{N\times q}$ is a square root of the covariance matrix, i.e. it satisfies:

In order to integrate this probabilistic description within a typical constrained-optimisation formulation, we define the set $\mathcal{D}_{\alpha}$ as the $\alpha$-level confidence region around the mean of $\Delta G_{r}^{\prime\circ}$. Using Eq. (7), and noting that the squared norm of a standard normal random variable is known to be Chi-squared distributed, we can represent this set as:

where $\mathbf{m}\in\mathbb{R}^{q}$ is a free parameter vector and $\chi_{q}(.)$ denotes the quantile function (inverse cumulative distribution function) of a chi-squared distribution with $q$ degrees of freedom. We can thus account for the uncertainty in thermodynamic estimates by treating the free energies of reaction as decision variables, rather than known parameters, and constraining their value to belong to ($\mathcal{D}_{\alpha}$). For a given reference flux distribution $\mathbf{v}$, this results in the following quadratically-constrained program (QCP):

where $b$ is the min driving force (or the MDF after optimisation), $\mathbf{c}\in\mathbb{R}^{m}$ is a vector of log-metabolite concentrations, $\mathbf{l_{c}};\mathbf{u_{c}}\in\mathbb{R}^{m}$ are (log) lower and upper bounds on these concentrations, $\mathbf{S}\in\mathbb{R}^{m\times N}$ is the stoichiometric matrix of the model, ${\rm R}$ is the Boltzmann gas constant, and $\rm T$ is the temperature used for the computation of the free energy estimates. For our MDF calculations, we used a confidence level of $90\%$ and set the bounds on the concentration of all metabolites to the physiologically plausible range of ($1\upmu\text{M}$, $10\text{mM}$). The above quadratically constrained program was formulated and solved the the GUROBI package [gurobi].

Probabilistic metabolic optimisation (PMO) of the model was performed using the PTA python package [gollub_probabilistic_2021], providing the software with the curated thermodynamic estimates generated for this model. The default values available through the package for growth on M9 medium with glucose (which include measurements from [gerosa_pseudo-transition_2015]) were used as priors for the concentration of metabolites, and the growth rate was bounded below by the reported value in [gerosa_pseudo-transition_2015]. Further, for this analysis, which requires every flux in the solution to have a well-defined directionality, two transhydrogenase reactions (NADH17pp, and THD2pp) were allowed to operate reversibly [gollub_probabilistic_2021].

To analyse the thermodynamic state computed by PTA, we compute the flux-force efficacy $\eta$ for each reaction as (see Figure 8):

where $\Delta_{r}G^{\prime}$ is the Gibbs free energy of reaction in the thermodynamic state, ${\rm R}$ is the universal gas constant and ${\rm T}$ is the temperature considered for the analysis ($310.15~{}K$). To compare the PTA-predicted flux force efficacies with experimental measurements of enzyme abundance, we selected reactions with carrying at least $2.5\%$ of the glucose uptake flux, and pooled these into two groups, corresponding to $\eta>0.5$ ($72$ reactions) and $\eta<0.5$ ($31$ reactions). $\eta=0.5$ corresponds to $\Delta_{r}G^{\prime}\approx-1.1~{}{\rm RT}=-2.8$ kJ/mol.