The proteome-dependent kinetic yeast model is based on a system of ordinary differential equations (ODEs) that describe the mass balances of all intra- and extra-cellular metabolites. This system of ODEs is solved with the ode15s function in MATLAB 2020b, for which the absolute tolerance is set to 1e-4, and all variables are constrained to be higher than zero with the “nonnegative” setting. A detailed description of the final proteome-dependent kinetic yeast model used is given in Supplementary Material S1.

To predict which proteome composition is the most competitive for defined environmental conditions, a Monte Carlo approach is used. The metabolic behavior of 1,000 random proteomes, generated around a seed proteome, is compared based on an objective function. Under steady-state conditions, the minimization of the residual substrate concentration was used as an objective function. Under dynamic conditions, the minimization of a time-weighted average substrate concentration was used, to promote fast consumption of the available substrate, therefore selecting for competitive proteomes:   ∫ 0 t c y c l e c s ⋅ t   d t ∫ 0 t c y c l e t d t .

Subsequently, it is determined whether the solution is balanced. If the objective function is optimized and the solution is balanced, the objective function and the seed proteome are updated. In the next iteration, the proteome is then generated around this new seed proteome, with a maximum deviation of 25% per sector.

All proteins from experimental datasets are sorted in the same nine protein sectors that are used in the model, to allow for a direct comparison of the experimental proteomes and the optimized proteomes. The proteins are categorized per sector based on either the protein name or the description in the KEGG database (Goffeau et al., 1996; Kanehisa et al., 2016) (Supplementary Material S7). The whole dataset is sorted with the MATLAB 2020b functions “strcmp” and “contains,” which are used to search the dataset for specific names or keywords to group the proteins.

The proteome cost parameters are estimated by optimization with the MATLAB 2020b function fmincon. For all parameter optimizations, a multi-start approach is used. This approach minimizes the risk of reaching a local minimum in the solution space by starting the optimization from different initial guesses. The tolerance of the function is set to 1e-12 for all optimizations. For the estimation of the k_cat parameters, the difference between the experimental and simulated fluxes is minimized. Additional weight in the objective function was applied for the growth rate, as k_cat parameters have to be rejected if the maximum growth rate is not reached.

The amount of overcapacity in the yeast proteome is determined by introducing a 10th protein sector. This new protein sector does not have a function for the cells, and hence, only takes up space in the proteome. Therefore, the fraction of the proteome that can be allocated into the extra sector without altering the metabolic fluxes is defined as overcapacity. The overcapacity is estimated for each sector separately, to minimize the changes in each step. The sectors are sorted in a decreasing order and then optimized for overcapacity in this order. The amount of overcapacity in each sector is determined in a step-wise approach. Per iteration, one percent of the specific protein sector is removed and allocated into the extra sector. Subsequently, the fluxes of the adapted proteome are compared to the reference fluxes, and only if the change in the fluxes remains within the boundaries, the seed proteome is updated. This new seed proteome is then used for the next iteration, in which the sector size is again decreased by 1%. By decreasing the sector size by 1% of the current size, the step size is reduced with each iteration. If the flux profile deviates more than the threshold value, the adapted proteome allocation is rejected. The fluxes are evaluated based on the following criterion: the average value of the uptake and growth fluxes should not deviate more than 1% from the reference flux, to ensure that the same substrate uptake and growth rates are achieved.

We wanted to construct a proteome-dependent kinetic model, which was small, but still able to reproduce the main phenotypes observed for S. cerevisiae. Furthermore, it should be calibrated with available experimental data. We constructed the model based on the kinetic model of yeast glycolysis (Teusink et al., 2000) which we extended with reactions for the trehalose cycle and respiration pathway as well as a growth reaction (see Figure 1). Each (lumped) reaction has been associated with a proteome fraction resulting in a proteome-dependent kinetic model of yeast central carbon metabolism and growth.

The Embden–Meyerhof glycolytic pathway has been implemented as three lumped reactions (uptake, upper, and lower glycolysis) with three intermediates: G6P, FBP, and pyruvate. The stoichiometry of the growth reaction was based on Suarez-Mendez et al. (2016). The NADPH requirement was assumed to be met by using the pentose phosphate pathway, which in sum (together with PGI) converts one G6P to six CO₂ and 12 NADPH. The required NADPH flux was balanced by a respective consumption of G6P. The ATP demand for growth has been derived from Della-Bianca et al. (2014) taking into account that the demand was expressed as catabolized glucose amounts. Furthermore, the trehalose cycle was included as two lumped reactions, based on an existing kinetic model of the trehalose cycle (Smallbone et al., 2011) (see Supplementary Material S1 for details).

Due to a lack of kinetic models of yeast TCA cycle and oxidative phosphorylation, the two respiratory reactions (from cytosolic NADH and pyruvate, v_NDE, and v_TCA, respectively) have been implemented using general Michaelis–Menten kinetics. However, the two reactions are interdependent—both connect to the electron transport chain—and consequently, the rate is determined by the same proteome fraction. A maximum value for the rate of the two reactions combined is defined, reflecting the capacity in the electron transport chain, limited by the provided proteome sector size (see Supplementary Material S1).

The biomass reaction contains many complex reactions, and the kinetics of the full process currently cannot be derived from basic principles. Therefore, a holistic approach based on experimental observations was chosen, i.e., the growth rate has been found to correlate with the energy charge (Boer et al., 2010). Here, the growth rate is described by a sigmoid function that is the most sensitive within the range of an energy charge between 0.7 and 0.9 as observed for growing cells (Boer et al., 2010).

The specific activity for the defined pathways has a major impact on model predictions. To obtain realistic values, the specific enzyme activities (k_cat) were estimated from experimental omics datasets. In the proposed model, the k_cat,i for each reaction i is defined as the maximum reaction rate per fraction of proteome (mol/Cmol_x/h), where 100% proteome reflects 500 mg protein per g_X (Ertugay and Hamamci, 1997). Hence, the maximum rate of the reaction i ( V m a x , i ) with a given sector fraction φ i   is: V max , i =   φ i ⋅ k cat , i . From this, c →   enzymatic rate V i is calculated by multiplying the V max , i with the function f i ( c →    ) describing the effects on the enzymatic rate due to substrate and product concentrations as well as effects by allosteric activators and inhibitors (see Supplementary Material S1 for specification of f i ( c →   ) for each reaction): V i =   V max , i ⋅ f i ( c →   ) .

The k_cat parameters have been estimated by combining the proteome and fluxome measurements under batch conditions. The proteome fractions were taken from de Godoy et al. (2008) using S. cerevisiae grown under batch conditions with a defined glucose minimal medium and aligned according to the protein classification in the KEGG database. Specifically, grouping all proteins with the KEGG BRITE label “Genetic Information Processing” and all proteins with the “Metabolism” label that were not classified as “Central Carbon Metabolism” or “Energy Metabolism” being assigned to the “growth protein sector,” assuming that their size is growth rate-dependent in the minimal medium. Furthermore, for the calculations, it was assumed that the whole proteome sector of cells grown under excess substrate at the maximal growth rate was used.

The corresponding flux distribution, i.e., under batch conditions was obtained from Heyland et al. (2009) with the exception of fluxes for the trehalose cycle—these were obtained from the feast/famine experiments conducted by Suarez-Mendez et al. (2017). For both trehalose synthesis and degradation, the maximum value of the flux reached in one feast/famine cycle was used, which was 5.10⋅10⁻³ mol/Cmol_X/h for trehalose synthesis and 4.09·10⁻³ mol/Cmol_X/h for the degradation of trehalose. The k_cat value for maintenance was set to 0.0155 mol/Cmol_X/h, which is the maintenance requirement measured at near-zero growth rates (Vos et al., 2016).

To obtain the k_cat parameters, parameter optimization was performed, estimating the parameters which produced the smallest deviation between the simulated and experimental fluxes (Heyland et al., 2009), using the batch proteome composition taken from de Godoy et al. (2008) (see Supplementary Material S2 for details). Using this approach, the proteome-dependent kinetic model was able to largely reproduce the experimental flux distribution (Table 1), and this k_cat calibration was used in all further calculations.

S. cerevisiae is a Crabtree-positive yeast, and thus fermentation is observed next to oxidative phosphorylation at substrate uptake rates above an observed “critical” rate (Barford and Hall, 1979). The ability of the model to reproduce the Crabtree effect is assessed by optimizing proteomes for dilution rates in the range from 0.05 h⁻¹ to 0.4 h⁻¹. The proteome optimization was started at the dilution rate of 0.4 h⁻¹ using the experimental batch proteome as a starting value. The most competitive proteome out of 1,000 randomly generated proteome allocations was selected using minimization of the residual substrate concentration as an objective function. Subsequently, this procedure was repeated for the next lower dilution rate. The optimal proteome allocation of the previous dilution rate was used as a starting value. To validate the model, the predicted fluxes and metabolite concentrations were compared with a flux and metabolome dataset (Suarez-Mendez et al., 2016) at different dilution rates under chemostat conditions. This comparison of predicted and measured fluxes and metabolite concentrations can be found in Supplementary Material S3 and in Supplementary Figures S2, S3, respectively.

The experimental data for ethanol production and oxygen consumption in Figure 2 show that the ethanol production starts at a dilution rate of 0.28 h⁻¹ (Rieger et al., 1983; Van Hoek et al., 1998). Above this critical dilution rate, the oxygen consumption rate decreases, while ethanol production keeps increasing. Ethanol production is first predicted by the model for a dilution rate of 0.25 h⁻¹, which is a lower rate than the experimental data. Furthermore, there is no decrease in the oxygen consumption rate above a dilution rate of 0.28 h⁻¹ for the optimized proteomes, which was observed in experimental studies (Van Hoek et al., 1998). From the model, this can be explained by the proteome-specific ATP production “cost”: Respiration has a high yield compared to fermentation (Table 2.). Hence, reducing the size of the respiration proteome sector will not be predicted by the model as it is not beneficial. The predicted plateau originates from a constraint that was introduced manually (12% of the proteome for respiration) to reflect the maximum oxygen consumption rate measured by Rieger et al. (1983) after long-term evolution. The continuous increase in the ethanol production rate can then be explained by the increasing need for ATP with an increasing growth rate while respiration is at its maximum.

This result conflicts with the discussed dataset of Van Hoek et al. (1998) as well as the model predictions of Nilsson and Nielsen (2016), which was partly based on this experimental dataset. This mismatch and conclusions will be discussed in more detail later. Notably, there is also experimental evidence from previous studies that the predicted plateau is reasonable. It was shown that the respiratory repression observed by Van Hoek et al. (1998) could be negated upon long-term adaptation (Barford and Hall, 1979; Rieger et al., 1983; Postma et al., 1989), and a stable maximum oxygen uptake rate above a dilution rate of 0.28 h⁻¹ was found.

To test the hypothesis of short- vs. long-term evolution, the proteome optimization approach was performed with a reduced number of generated proteomes and compared to the experimental data of Van Hoek et al. (1998) (Figure 2, red line). With a high number of generated proteomes for the optimization, the experimental findings of long-term chemostats could be reproduced. From these predictions, we hypothesize that cells not exposed to long-term glucose-limited conditions did not yet reach the “optimal” proteome allocation and respective metabolic phenotype. This set number of 1,000 simulations was chosen because only very limited further optimization of the objective function was observed after this number of simulations. As such, 1,000 simulations were concluded as sufficient to reach the optimum. Work on adapted glucose-grown cultures, at which point glucose repression on respiration disappears, is cultivated for at least 50 generations at the same dilution rate (Barford and Hall, 1979). A work by Van Hoek et al. (1998) describes the Crabtree effect with its typical glucose repression of respiration, by cultivating cultures at the same dilution rate for seven generations. Therefore, a set number of 100 simulations was chosen to reflect this state of limited adaptation of the proteome from batch growth conditions.

Looking into the global trends in the fully evolved proteome allocation at different dilution rates (Figure 3, see Supplementary Material S4 for sensitivity analysis), an increase in the dilution rate can be seen for nearly all sectors leading to the unused space (in the following called overcapacity sector, last panel). The overcapacity sector accounts for the fraction of the proteome which remains unused within the optimized proteomes. Before discussing specific trends, the high dilution rates will be highlighted. Even close to the maximal growth rate, the model predicts a small overcapacity sector. Nevertheless, please note that the batch and very high dilution rate might still have different optimization criteria; here, in the model, minimal substrate concentration was applied as the objective function. Because of the optimization approach, some robustness is required that was not further tuned as the fraction is rather small (7%) and does not change trends. Additionally, the algorithm samples from an enumerated number of randomly generated proteomes and therefore requires some buffer for robustness.

A major difference between this model and earlier studies (Nilsson and Nielsen, 2016) is that the proteome space limit is not reached at the critical growth rate (D = 0.28 h⁻¹). At the critical dilution rate (D = 0.28 h⁻¹), the overcapacity sector still has a significant fraction (21%). As briefly discussed earlier, Nilsson and Nielsen (2016) postulated that the Crabtree effect could be explained by the catalytic efficiency of the fermentation and respiration pathways expressed as ATP per amount of protein used in the pathway (Table 2). To estimate these catalytic efficiencies, Nilsson and Nielsen (2016) used the fluxes and specific enzyme activities for fermentation and respiration, under the assumption that all enzymes operate at half of their maximum specific activity , whereas in this model, the estimation of the catalytic efficiency is based upon the proteome and fluxome dataset, using dynamic saturation of enzymes. The estimation proposed by Nilsson and Nielsen (2016) subsequently produced a proteome composition in which the mass of all respiration proteins is 19 times larger than the protein mass of all glycolysis enzymes , while from proteome measurements it was observed that the mass of all respiration proteins is 0.3 times the size of the mass of all glycolysis proteins (de Godoy et al., 2008; Elsemman et al., 2022). This large difference in proteome allocation between glycolysis and respiration causes the catalytic efficiency of fermentation to be overestimated. The conclusion that the proteome is fully allocated after the critical growth rate is reached leads to the prediction that the “optimal” endpoint of proteome allocation is reached, which cannot explain datasets by Barford and Hall (1979); Rieger et al. (1983). Additional modeling studies by Elsemman et al. (2022) suggest that the decrease in oxygen consumption at higher growth rates observed by Van Hoek et al. (1998) is not caused by a limitation in proteome capacity but rather by a maximum rate of mitochondria biogenesis, in which long-term adaptation could overcome the described glucose repression of respiration.

The proteome compositions, especially at low dilution rates were characterized by a significant overcapacity sector. The kinetic proteome allocation approach could not yet answer why the cells maintained such an excess proteome. As discussed earlier, the hypothesis for a proteome overcapacity is competitiveness and robustness including dynamic environmental conditions. Overcapacity could enable faster substrate uptake rates and enable a competitive advantage and outcompete slower consuming microbes (Jannasch, 1967). Furthermore, excess capacity could enable a robust, balanced functioning of pathways such as glycolysis (van Heerden et al., 2014) under dynamic substrate conditions.

To test these hypotheses, we studied the predicted proteome allocation under different repetitive substrate-feeding regimes using the proteome-dependent kinetic model, using the minimization of the time-weighted residual substrate concentration as the objective function. With this approach, we were able to select competitive proteomes with fast substrate uptake rates. As a reference dynamic condition, an experimentally explored feeding regime was chosen, i.e., a cycle length of 400 s of which 20 s was used to feed the culture (D = 2 h⁻¹), leading to the average dilution rate of D = 0.1 h⁻¹ over the complete cycle (Suarez-Mendez et al., 2014).

Proteome allocations and respective metabolic phenotypes were then compared to the steady-state at the same (average) growth rate. First, we studied the maximum, minimum, and average enzyme saturation (V/V_max) under dynamic conditions compared to the enzyme saturation under chemostat conditions (Table 3). Under dynamic conditions, the maximal enzyme saturation is much higher (up to 92% for the respiration reaction) than that under chemostat conditions (77% for respiration). Nevertheless, the average enzyme saturation over the whole cycle is actually lower than that under the reference chemostat state (for respiration, 25% compared to 77% at steady-state). This indicates that the proteome optimization to some extent focuses on the “peak” flux, especially for the large sectors of respiration and growth, indicating high usage of the available flux capacity while on average leaving a large overcapacity over the whole cycle. This enables a rapid consumption of the substrate as soon as it becomes available, which was the optimization criteria.

We were now curious to see how the perturbation strength would influence the proteome allocation. Therefore, the length of the feeding period was varied at the same average dilution rate, resulting in different substrate perturbation intensities. The respective predicted proteome allocations were calculated and compared (Figure 4) for the different ratios of feeding time over cycle time (TF/TC). TF/TC values were chosen as log2 increments from the experimentally used TF/TC value of 1/20 (Suarez-Mendez et al., 2014).

The growth sector fraction increased with the perturbation intensity, suggesting that this strategy was the most effective measure to survive the higher substrate concentration variations (from faster feeding) and consequently high flux dynamics. The growth reaction seemed to act as an efficient and fast sink for substrate and ATP. However, in reality, the growth sector does not consist of a single reaction and may not be able to provide a rapid response upon glucose influx. For this reason, two other scenarios were additionally evaluated: 1) the regulation of the trehalose cycle upon repeated substrate pulses and 2) the regulation of the ratio between upper and lower glycolysis (see Supplementary Material S6).

The trehalose cycle has been described to function as a “safety valve” upon large changes in the glycolytic flux (Thevelein and Hohmann, 1995; Blomberg, 2000; van Heerden et al., 2014; Vicente et al., 2018). Under dynamic conditions in yeast, it was found that a significant amount of imported glucose was recycled through the trehalose cycle, especially during periods of high flux changes (Suarez-Mendez et al., 2017). To evaluate the effect of storage metabolism activity under dynamic conditions, the reference condition [D = 0.1 h-1, TF/TC = 0.05, (Suarez-Mendez et al., 2014)], was further analyzed. We varied the trehalose sector size between 0 and 1% (Figure 5) and compared the response of metabolism using FBP and Pi as indicators. A balanced metabolic response will lead to repetitive cycles in FBP and Pi. Such repetitive response was observed for proteomes with a trehalose sector larger than 0.1%. Increasing the trehalose sector above 0.1% leads to reduced fluctuations in G6P/FBP and Pi, suggesting a more robust metabolic response. Simulated changes in FBP and Pi are in line with results from previous work by van Heerden et al. (2014).

To evaluate the prediction accuracy and trends of the predicted proteomes under dynamic conditions, the simulated proteome adaptation from chemostat to feast/famine conditions was compared with the experimentally measured proteome fold changes between chemostat and feast/famine conditions (Verhagen et al., 2022) (Figure 6). Proteins of trehalose/glycogen storage, ribosomes, and oxidative phosphorylation were used as proxies for the storage, growth, and respiration sectors, respectively (proteins categorized in the same way as calibration approach, see Methods).

The model predicted the experimentally observed changes in upper and lower glycolysis (Figure 6). The enzyme TDH catalyzes the glyceraldehyde dehydrogenase reaction (TDH), which forms ATP using Pi. However, if upper and lower glycolysis reactions are imbalanced during high fluxes, this reaction becomes a bottleneck, leading to the accumulation of FBP and subsequently to an imbalanced metabolism. Therefore, it was expected that TDH had to be upregulated under substrate-fluctuating conditions to facilitate balanced intermediates, which was reflected in both the model predictions as well as in the experimental dataset. The predicted change in lower glycolysis is larger than that in the experimental data. This is likely caused by the fact that simulated proteomes for chemostat conditions contain no overcapacity in the lower glycolysis sector, whereas experimental proteomes under chemostat conditions appear to contain more overcapacity in this proteome sector. As such, the fold change between measured and simulated values is higher. Furthermore, the model reproduced the average change observed for the uptake sector, although it should be noted that effects of individual iso-enzymes (especially with regard to HXK/GLK, which catalyzes the first step of glycolysis) were not taken into account in the current model.

Significant deviations between experimental and predicted fractions were observed for the storage sector. This was significantly decreased experimentally, while the resource-dependent kinetic model predicted an increase. Experimentally, a decrease of 28%, from 0.25 to 0.2% of the proteome, was observed, while an increase to 0.2% of the proteome was predicted in the model. Possible reasons for this difference in sector size could be: 1) the synthesis of trehalose has additional functions in the cell which are not represented in the model—it is described that trehalose plays an important role in different stress responses, including severe substrate limitation at low dilution rates (see also Supplementary Figure S3). 2) The measured and predicted proteomes do not include neither post-translational modifications, which are known to significantly affect the k_cat’s of enzymes in the trehalose cycle (Sengupta et al., 2011), nor changes that could occur during cell-cycle progression.

Furthermore, there could also be a bias from the experimental setup–the differences in the trehalose sector, combined with the observed increase of the lower glycolysis sector, suggest that the experimental chemostat proteome is potentially already primed for dynamic environments, compared to experimental conditions, and as such is more robust than the predicted optimized chemostat proteomes.