Emergent Subpopulation Behavior Uncovered with a Community Dynamic Metabolic Model of Escherichia coli Diauxic Growth

A single GEM was used to model the average E. coli metabolic state, and we compared the simulation results with the original data from Varma and Palsson (18) (see Fig. S1 in the supplemental material). The parameters for the simulations are reported in Table 1 and 2, and the only flux constraints that we calibrated to the data were the oxygen uptake rate and the maximal acetate secretion rate. A fixed cell death rate (Table 1) was introduced as previously described, using a value from the literature (19). In these simulations, a lower absolute level of flux variation at each simulation time step was observed with the MOMA implementation (Fig. S2). We used the same GEM to reproduce the results from Enjalbert et al. (22), changing only the initial values for biomass, glucose, and acetate (Fig. 1; see also Fig. S5). Although the pFBA simulation (Fig. 1a) showed a brief shift to growth on acetate at the time of glucose exhaustion (GE) (∼4 h), the MOMA simulation predicted complete growth arrest already occurring at that point, with a minimal acetate consumption to satisfy the ATP maintenance requirement implemented in the GEM (Fig. 1b). Both simulations well captured glucose consumption and acetate secretion, but neither was able to reproduce the slow acetate consumption observed experimentally. Even after fine-tuning the constraint on acetate uptake to achieve a perfect match of the acetate consumption data from Varma and Palsson (18), the model could not reproduce the acetate concentration dynamics of the corresponding data from Enjalbert et al. (22) (data not shown). Therefore, we decided to avoid fine-tuning of the acetate uptake (Table 1). Both the pFBA and MOMA simulations showed an abrupt change in the flux distribution upon shifting from glucose consumption to acetate consumption (Fig. S3). We evaluated the agreement between the experiment and the simulation with the R² distance between in vivo and in silico data for biomass (pFBA R² = 0.989; MOMA R² = 0.982), glucose (pFBA R² = 0.993; MOMA R² = 0.993), and acetate (pFBA R² = 0.277; MOMA R² = 0.409). In Fig. 2, we compare the flux distributions of our simulation results to the experimental results reported by Enjalbert et al. (22) for overexpression/underexpression of key genes associated with glucose and acetate metabolism (represented graphically in the top panels). First, we computed the flux solutions for E. coli growing on either glucose or acetate exponentially (data not shown) and compared the fluxes through the relevant reactions in E. coli growing on acetate to those in E. coli growing on glucose. Fig. 2a shows the absolute values for the flux results in the two simulations, normalized to values between 0 and 1 for direct comparison with the qualitative representation of the gene expression data (with a value of 0 for nonexpressed genes and a value of 1 for expressed genes). The simulation results were consistent with the results of the experiments, with active reactions (dark green) related to acetate consumption and anabolism (ACKr, PPCK, FBP, ICL, MALS) and inactive reactions (white) related to glycolysis (PFK and PYK) during growth on acetate and vice versa during growth on glucose. PPS did not carry flux in either simulation. We then used the simulation results presented in Fig. 1 to compare the metabolic fluxes before and after glucose exhaustion (GE), i.e., before and after the single E. coli model shifted from growth on glucose to growth on acetate. Enjalbert et al. (22) compared gene expression levels between samples taken at time GE plus 30 min and at time GE minus 115 min. However, Fig. 1 shows that according to the simulation, growth had already stopped after 30 min from the GE point. Indeed, comparing the absolute values of fluxes taken at time GE plus 30 min and at time GE minus 115 min, we found that both the pFBA and MOMA simulations qualitatively captured the downregulation trends, whereas neither simulation reproduced the observed upregulation (data not shown). Fig. 2b shows the difference in absolute values of fluxes taken at time GE plus 18 min and at time GE minus 115 min, time points where growth is still observed in pFBA simulations. In this case, both simulations qualitatively captured most of the upregulation/downregulation trends. Figure S4 shows the metabolic network (modified from the map for the E. coli core model constructed by the use of Escher [28]), with the data from the reactions performed as described for Fig. 2b highlighted and color-coded according to the gene expression data. Finally, we reproduced the other experimental scenarios from Enjalbert et al. (22) with the uniform population model, adjusting only the initial values of biomass, glucose, and acetate. We observed that a uniform shift was able to reproduce the biomass profile well when high acetate concentrations were present in the medium (Fig. S5c), while this was not the case when only low acetate concentrations were available (Fig. S5b).

We used two GEMs (and the same parameter values as before) to model E. coli monocultures as a mixture of two populations, one adapted to grow on glucose and one adapted to grow on acetate. The two models, the E. coli glucose (EC_Gl) model and the E. coli acetate (EC_Ac) model, were hence constrained to exclusively take up the corresponding carbon source. Two transition functions, dependent on acetate or glucose concentrations, were introduced to model cellular differentiation and cellular shift from one population to the other (see Materials and Methods for details). We ran simulations to compare the different scenarios investigated experimentally by Enjalbert et al. (22). The initial values for biomass, glucose, and acetate were adjusted to the corresponding data sets. The transition rates, as well as the initial population ratios, were chosen following the assumption, supported by a simple mathematical model, that populations in constant environments will converge to a constant ratio (see Text S1 in the supplemental material for details). Data in Fig. 3 show simulations performed under the same conditions as those described for Fig. 1a, with the same absolute initial biomass values, distributed in this case as 95% EC_Gl and 5% EC_Ac. This initial ratio was chosen by considering the range of steady-state values for the population ratio (reported in Table S1 in the supplemental material) as well as considering that it is reasonable to assume that a higher number of cells would be adapted to grow on glucose, which is the carbon source on which laboratory cultures are usually maintained. Data in Fig. 3a show the simulation results for a scenario without transitions between the two states, whereas the results of Fig. 3b were obtained with active transition functions, defined here by constant transition rates as reported in Table 2. Although both panels a and b of Fig. 3 capture well the biomass (R² = 0.987 and R² = 0.990, respectively) and glucose concentrations (R² = 0.996 and R² = 0.997, respectively), only the simulation that included the population transition realistically reproduced the acetate consumption levels (R² = 0.336 and R² = 0.951 respectively) as well as a lag phase before culture crash. Neither of the simulations captured the eventual recovery of growth hinted at by the last data points. We reproduced two other results (where only biomass measurements were available) from Enjalbert et al. (22), again using the same GEMs and changing only the initial conditions (biomass quantity and distribution among EC_Gl and EC_Ac) and the experimental setup accordingly. By modeling the population transition with the same constant rate, we were able to explain the biomass profile in the case where E. coli was grown on 15 mM glucose and, after glucose exhaustion, the acetate concentration was maintained at around 4 mM (Fig. S5e) (R² = 0.986), but not in the case where E. coli was grown on 15 mM glucose and 32 mM acetate and, after glucose exhaustion, the acetate concentration was maintained at the same high level (Fig. S5f) (R² = 0.727). We therefore introduced a dependency of the transition functions on the substrate concentration (see Materials and Methods for details) that well captures all the experimental scenarios with a unique set of parameters (Fig. S5g, h, and i). Data in Fig. 4a show that an E. coli population starting with 95% EC_Gl and 5% EC_Ac describes well the biomass dynamics (R² = 0.985) and matches the glucose exhaustion point, observed after around 4 h when acetate was maintained at 4 mM. Again, without fine-tuning the GEM simulation parameters, Fig. 4b shows that an E. coli population starting with 75% EC_Gl and 25% EC_Ac reproduced the biomass measurements (R² = 0.940) and the glucose exhaustion point after around 4 h also in the experimental setup with acetate maintained at 32 mM. The effect of adjusting the initial biomass ratios in the different experimental conditions is shown in Fig. S6. Overall, the simulations starting with 95% EC_Gl and 5% EC_Ac or 75% EC_Gl and 25% EC_Ac did not show strong differences, but further reducing the percentage of EC_Gl (and leaving the range of steady-state values of Table S1) resulted in drastic changes to the shape of the growth curves. The initial condition of a 75% EC_Gl and 25% EC_Ac population distribution for Fig. 4b is also justified by a difference in the initial experimental values for the biomass quantity (see Fig. S7).

Enjalbert et al. (22) showed different trends in the lag times of E. coli cultures required to achieve maximal growth after GE. In their switch experiments, they sampled at different time points “mother cultures” of E. coli cells growing under batch conditions on 15 mM glucose alone (“M9G” condition) or on 15 mM glucose and 32 mM acetate (“M9GA” condition) and reinoculated the sampled cells as “daughter cultures” into fresh medium exclusively containing glucose (M9G condition) or acetate (M9A condition). We replicated this experiment in silico by running first simulations under the M9G and M9GA conditions. For the M9G mother culture, we used the simulation of the mixed EC_Gl and EC_Ac population shown in Fig. 3b, because the experimental conditions were the same. We did not have an experimental reference data set for the M9GA mother culture, and we simulated a new scenario similar to that shown in Fig. 4b, with the same initial population composed of 75% EC_Gl and 25% EC_Ac but without the feeding of additional acetate. The GE time points were about 4.6 h for M9G and 3.9 h for M9GA, consistent with the observations of Enjalbert et al. (22) (data not shown). The in silico mother cultures were sampled at regular time intervals to obtain the initial biomass distribution of EC_Gl and EC_Ac for the daughter cultures (reported in Table 3), and the lag time was computed for each daughter culture (see Materials and Methods for details). Data in Fig. 5 show the simulation results compared with the experimental data from Enjalbert et al. (22). The error bars were obtained for the simulated lag time data by adjusting the initial biomass ratio of the daughter cultures by ±15%. A quantitative agreement between simulation and experimental results was achieved only in the M9G-M9G switch experiment (Fig. 5a) with the correct prediction of almost zero lag time for the daughter cells, but the trend for the delay to reach maximal growth was in general qualitatively reproduced also for the other scenarios. According to the simulations, cultures switched from M9G to M9A (Fig. 5a) need about 1.5 h before reaching maximal growth, which is more than twice the duration observed experimentally. For cultures pregrown in M9GA (Fig. 5b), we observed both in simulations and in experiments a decreasing lag time for daughter cultures sampled after GE for the M9GA-M9A switch and an increasing lag time for the M9GA-M9G switch. Additional studies are represented in Fig. S8. In particular, panels a to d of Fig. S8 show the dependence of the lag time in the daughter cultures on the maximal transition values and panels e to h of Fig. S8 show the same dependence, including the distribution of the biomass ratio in the mother cultures, for a limited set of parameters.

One recurrent criticism of stoichiometric and constraint-based modeling approaches, such as FBA, is that they can easily be adjusted to reproduce experimental results by ad hoc changes of flux constraints. Indeed, we observed that a condition-specific fine-tuning of the constraint on acetate uptake could reproduce fairly well the growth dynamics of the different experiments (data not shown). However, the change of such a constraint from one experimental condition to another is not biologically justified. Although some extensions of the FBA approach such as FBA with molecular crowding (FBAwMC [29]) provide reasonable ways to constrain the metabolic fluxes and were shown to reproduce carbon consumption hierarchies, they also require extensive parameterization. We therefore chose to use the basic FBA approach, limiting the number of constraints imposed and with parameters mostly from experimental measurements (Table 1). In the case of oxygen uptake and acetate secretion, we calibrated the constraints using data from Varma and Palsson (18), where an E. coli diauxic shift from glucose to acetate was first simulated using a genome-scale model. The FBA parameters were left unchanged to reproduce the independent experiments reported by Enjalbert et al. (22). The use of an independent set of data to calibrate the FBA model parameters is a possible way to improve the confidence in subsequent results. Further model parameters of the ODE system were chosen according to reasonable hypotheses and were adjusted slightly to achieve fair agreement with the experimental results in a manner that was consistent among all the simulations. The initial conditions were specific to the experiments that we aimed to reproduce.

The flux distributions obtained from FBA solutions represent an average picture of the metabolic state of a population, which is in general modeled using a single genome-scale model. Therefore, standard dFBA implementations, in which the FBA constraints evolve according to the updated external conditions, reproduce the average change in metabolic state of the population in response to external variations. This is equivalent to assuming that a population undergoes a coordinated, uniform metabolic shift under changing environmental conditions. Furthermore, such transitions are generally abrupt with dFBA models. We therefore tested two alternative approaches to simulate the diauxic shift in uniform E. coli monocultures, solving the FBA problems either with pFBA (mostly equivalent to the usual dFBA implementations) or with an adaptation of the MOMA algorithm. In the latter case, instead of minimizing the difference in flux distribution between a “wild-type” GEM and a modified one (original MOMA implementation), we used the same concept to integrate the dFBA system while also imposing the following condition: at initial time t_i, the flux solution differs minimally from that at time t_i-1, where the time steps are set by the integration routine. In contrast to our expectations, however, this approach did not achieve smoother metabolic adjustments in the system in response to the changing external conditions. Instead, both implementations resulted in abrupt changes in the flux distributions following the shift from glucose to acetate metabolism (see Fig. S3 in the supplemental material). More-sophisticated implementations of a dynamic MOMA model (e.g., computing the minimal adjustment based on a subset of biologically relevant variables) might succeed in achieving smooth metabolic transitions but would require the introduction of additional parameters and ad hoc biological hypotheses. In a similar way, biologically justified extensions of FBA such as FBAwMC (29) might provide better descriptions of an average and uniform population-level metabolism but typically need the empirical determination of large numbers of organism-specific parameters.

With the introduction of two basic assumptions (first, that there are two distinct metabolic states consuming either glucose or acetate; second, that transition from one state to the other is driven by glucose and acetate concentrations), we were able to capture all the experimental trends published by Enjalbert et al. (22) with the same computational model. The transitions between the two states were modeled as Hill functions of the corresponding substrate concentrations with a noise offset representing a constant, small noise component in cell regulation. Although other transition laws could have been chosen, Hill functions conveniently model concentration-dependent shifts between two states. For example, when acetate is highly abundant, more cells in the glucose consumption state shift to the acetate consumption state in response to the change in environment. Finally, the introduction of a transition efficiency term was motivated by the observation that cells can get “lost in transition,” an effect that was estimated to account for the death of ∼7% of yeast cells, which cannot initiate glycolysis following a shift to high glucose levels (30). Using a simple mathematical model (see Text S1 in the supplemental material), we identified ranges for the parameters of the transition functions and selected reasonable values that would return good agreement between simulations and experiments. Both the values for the constant transition rate (4% h^–1) and the values for the maximal transition rate (20% h^–1) were in good agreement with measured average protein turnover rates in E. coli cultures from the literature (31–33). Simulation results were mostly in very good agreement with the experimental data, and our results strongly further support the idea, suggested over the last few years by results from independent studies of different organisms (21, 25, 34), that monocultures represent an ensemble of subpopulations in different metabolic states that are partially regulated by the environmental conditions. When the acetate concentrations were too low to support growth, it was sufficient to model the transition as a constant random process. In contrast, in order to reproduce the data under conditions with high acetate concentrations, we needed to introduce an active transition rate dependent on substrate concentrations. Interestingly, this assumption alone was sufficient to model the experimentally observed growth rate, without further fine-tuning of model parameters. The introduction of substrate-dependent transition functions is also consistent with the experimental observations of Kotte et al. (27), supporting the hypothesis that a monoculture undergoes diversification in response to environmental changes.

With standard dFBA simulations, the metabolic transition during the shift from one carbon source to another is abrupt, and no lag phase is observable. Such an outcome is rarely the case, and, most remarkably, the duration of the lag phase between the exhaustion of the favored carbon source and the resumption of optimal growth on the alternate carbon source is highly variable under different environmental conditions. This observation can easily be explained as an emergent property of subpopulation dynamics. Our simulations are consistent with the explanation that the delay in the resumption of full growth actually depends on the relative abundances of the two subpopulations. Although the simulation results did not reproduce the experimental data quantitatively, all qualitative trends were fully explained. Several factors may explain these discrepancies. For example, the lack of experimental data concerning the mother cultures (in terms of biomass, glucose, and acetate dynamics) made it impossible to calibrate the initial model population. This could introduce a significant bias in the later sampling and determination of the subpopulation ratio, thus strongly influencing the quantification of the lag time, which is highly correlated with the population distribution (Fig. S8). Solopova et al. (25) showed that the density of a Lactococcus lactis population (translating in practice to the rate at which the primary carbon source was consumed) played a significant role in determining the proportion of cells successfully transitioning to growth on the secondary carbon source. The connection between lag time and subpopulation distribution could in principle be exploited to estimate initial population distributions from lag time measurements. However, it is difficult to assess the robustness and reliability of such predictions with the currently available data, and further investigation, including experiments devoted to determination of initial conditions, is therefore required. An additional source of the discrepancies between our quantitative results and the experimental measurements could have been the experimental procedure itself. For example, abrupt changes in conditions, such as the reinoculation of daughter cultures into a different medium in the switch experiments, might select for additional adaptation strategies. Interestingly, we observed a dramatic improvement in the quantitative agreement between experiment and simulation by relaxing the condition imposing no growth for populations inoculated on the “wrong” carbon source (data not shown). By allowing the glucose-consuming population sampled from glucose mother cultures to growth more slowly on acetate, we mimicked a situation in which cells store resources and are able to survive a bit longer. On the other hand, allowing reduced growth on acetate (glucose) for the glucose consumer (acetate consumer) population that was exposed to both carbon sources in the mixed mother cultures could be a proxy for a memory effect. Bacterial cells do show memory effects in response to changes in environmental conditions (26), but more-systematic experiments would be necessary to carefully and reproducibly determine the lag times as functions of external parameters to explore this potential explanation further. Finally, data obtained using a recent stochastic model of the regulatory network of diauxic growth in E. coli suggest that the limitations of biological sensors are responsible for the lag phase (35). From these results, we can infer that the transition functions, which currently depend on the absolute concentration of one carbon source at a time, might not be able to capture the fine details of population shifts in our model. A possible extension would be to introduce more-complex transition mechanisms dependent on the relative concentrations of primary and secondary carbon sources, a process whose elucidation would need dedicated experiments for the construction and validation of the new transition functions.

We developed a modeling framework to perform FBA simulations using embedment in a system of ODEs. Building on previous methods and approaches (19, 36), we further extended the standard dFBA implementation and introduced novel concepts. In particular, standard dFBA approaches assume that fluxes can instantaneously change to adapt to new environmental conditions, and flux solutions at subsequent time steps might differ significantly. This is an obvious limitation when aiming to capture diauxic shift, where lag phases, highly dependent on the environmental conditions, are typically observed. We implemented the MOMA algorithm (originally developed to model the response to genetic perturbations in static FBA) in dFBA to minimize the metabolic adjustments between different time points. Furthermore, we integrated dynamic mechanisms into dFBA that cannot be included in metabolic models, such as population transitions. Indeed, although the use of dFBA to model subpopulations bears some similarities to the use of other platforms for the simulation of microbial communities, a notable difference in our formulation is the capacity of subpopulations to interconvert. The current study relied on the a priori knowledge that only two carbon sources would be available to E. coli, thus motivating the development of a two-subpopulation community, but in principle, an arbitrary number of subpopulations can be defined and more generic transition functions introduced. Further experiments, in particular, single-cell studies, could be designed to define and parameterize these transition functions. Thanks to the object-oriented (OO) design of the framework, it is relatively easy to introduce other functions regulating the constraints on specific reaction fluxes in the FBA problem. In this way, different hypotheses can be extensively tested to improve understanding of how to capture regulatory dynamics in dFBA. Notably, the methods developed in this framework to study population heterogeneity could then be transferred to other platforms that are more specific for microbial community modeling where different features are implemented (e.g., spatial structure [19] or community-level objectives [37]). Finally, the framework could also be developed further to include stochastic mechanisms, such as mutations that would alter the function of metabolic genes. Indeed, our implementation of the dFBA algorithm is able to call different methods at each time step, e.g., to update the flux rates, and a regulatory function with random components could in principle be defined.

There is extensive experimental evidence that bacteria differentiate into subpopulations as a result of survival strategies (25, 27). Simulations based on standard dFBA model the dynamics of cells by predicting the putative average behavior of a whole population. For example, if a population of cells globally utilizes a combination of two carbon sources, dFBA would predict metabolic states in which the two carbon sources are utilized simultaneously. Our model assumes that cells are either in the glucose-consuming state or the acetate-consuming state, with an instantaneous transition between these two subpopulations that follows a simplistic rule which cannot capture intermediate states. This simplification is both practical and plausible for observing population dynamics as the emergent properties of individual behavior, and it works well in dynamically changing environments with a continuous transition. However, rather than having a well-defined metabolic state, especially during the transition between states, cells might exhibit a mixed state, which could be described as a superposition of “pure” states, analogous to the state vectors in quantum physics. Furthermore, our approach suggests a fundamental difference in the strategies used to account for metabolic fluxes in heterogeneous populations, because the average fluxes in a uniform population might differ from the cumulative average fluxes of subpopulations. Further investigations of this novel concept of superimposed metabolic states will provide a promising new approach to study the principles of metabolic regulation.

In stoichiometric models, the stoichiometric matrix S(m×n) is defined with the row and column dimensions corresponding to the numbers of metabolites m and reactions n respectively, the elements s_ij being the stoichiometric coefficients of metabolite i taking part in reaction j. FBA defines and solves the following LP problem:subject to:with l.b._j and u.b._j representing, respectively, a lower and upper bound on flux v_j.

The steady-state assumption (equation 2) gives a system of equations that is underdetermined and has an infinite number of solutions. Constraints on the fluxes (equation 3) allow us to restrict the solutions to a convex solution space but still result in an infinite number of solutions. The definition of an objective (equation 1) selects one solution, but this is still generally not unique for large (genome-scale) metabolic networks.

We consider two extensions to the FBA problem definition, namely, pFBA (14) and MOMA (15). We then use these two methods to solve the FBA problem in an approach similar to dFBA (16). Assuming that metabolism evolves toward the efficient utilization of resources, pFBA finds the minimal flux distribution that returns the same objective as that defined by the FBA problem. We use the pFBA implementation from COBRApy (38) with maximal flux through the biomass reaction as the objective function. Considering that metabolism must respond quickly to perturbations, MOMA implements a quadratic algorithm to find the FBA solution after gene deletion that is most similar to the optimal wild-type configuration. In our case, we do not introduce modifications to the metabolic network but rather require that the MOMA solution obtained at time t_i-1 is used to compute the MOMA solution at time t_i as the minimally different solution that satisfies the objective function. Also in this case, the objective function is maximal flux through the biomass reaction. We use the MOMA implementation from COBRApy (38) in the linear approximation, with a slight modification to allow the LP problem to be reset in an iterative manner, which is necessary to run MOMA within the dFBA approach.

In the SOA of dFBA, the boundary conditions in equation 3 are updated at discrete time steps according to the solution of the FBA problem in the previous time interval. Assuming quasi-steady-state conditions, i.e., that metabolism readjustments are faster than external environmental changes, dFBA can approximate the dynamic response of a GEM to a changing environment. Our approach is an extension of dFBA. The model is built as a system of ODEs, whose dimension depends on the dynamics to be modeled. Each ODE describes the variation in time of biomass, metabolites, or other regulatory/dynamic processes. The biomasses and the metabolites can be but are not necessarily linked to the corresponding variables in a GEM. Their ODEs vary according to a function that can then depend on the flux solutions $\vec{v}$ as follows:

The ODE system is then solved using integration routines with an automated choice of time step. Each integration step solves the FBA problem (or pFBA or MOMA problem) to obtain the current reaction rates for equation 4, updates the metabolite quantities according to the FBA solution, recomputes the flux boundaries of equation 3 according to specific reaction kinetics (typically Michaelis-Menten enzyme kinetics), and redefines the FBA problems with the new boundaries and/or other regulatory mechanisms defined by the user.

The modeling framework is written in Python (Python Software Foundation, https://www.python.org/) following the object-oriented programming (OOP) paradigm for efficiency and flexibility. The framework uses functionality from the following third-party packages: numpy (39), scipy (40), matplotlib (41), COBRApy (38), and pandas (42). In particular, we use COBRApy methods to solve the FBA problems and Python integrators from the scipy.integrate method ode to solve the system of ODEs.

We used a previously reported core version of E. coli GEM (43) downloaded from http://bigg.ucsd.edu/models/e_coli_core. The E. coli EC_any model is constrained with respect to the consumption of “any” carbon source (i.e., glucose [Gl] and acetate [Ac]) solely by the environmental conditions, and the lower bound of the exchange reactions (EX_Gl_e and EX_Ac_e, respectively) follows two simple Michaelis-Menten kinetics equations:

The ODE system is defined as follows:where v_μ is the reaction rate of the biomass function (proxy for growth rate) in the FBA model, δ is the cell death rate, and ξ_fed-batch is a positive rate under fed-batch conditions and zero under batch conditions. Parameters and initial conditions are summarized in Table 2. Either pFBA or MOMA can be used to solve the FBA problem.

Two E. coli core models are loaded and defined as either a glucose consumer (EC_Gl) model or an acetate consumer (EC_Ac) model by switching off uptake of the other carbon source as follows:

The ODE system is defined as follows:where ζ · H(t − t_x) is a heaviside function activated at time t_x of glucose exhaustion in order to keep the acetate level constant, ψ and ϕ are functions that model the cellular shift from EC_Gl to EC_Ac and from EC_Ac to EC_Gl, respectively, and 0<ε<1 is a positive factor representing the transition efficiency. The ψ and ϕ functions are modeled as Hill functions with a noise offset as follows:where they are constant transition rates for $V_M^{\varphi }=V_M^{\psi }=0$ . For the simulations presented here, we used a Hill coefficient value of n = 5. Indeed, the simulations seemed to work best for a transition function with a high degree of cooperativity, and the results are robust with respect to small deviations relative to this value. The other parameters and initial conditions, specific to the different simulations, are summarized in Table 2. For mixed-population simulations, pFBA is used to solve the FBA problem.

Two E. coli mixed-population model simulations are run as “mother cultures” as shown in Table 2 for M9G and M9GA conditions (glucose and glucose plus acetate, respectively). From each mother culture, we sample 11 time points between –1 and +1.5 h from the corresponding GE time (4.6 h for M9G and 3.9 h for M9GA) to obtain the biomass ratio between EC_Gl and EC_Ac used as the initial condition for the reinoculation simulations. The percentage of EC_Gl biomass at these time points is shown in Table 3. The daughter cultures are then grown under M9G glucose-only or M9A acetate-only conditions (see Table 2), yielding 44 simplified simulations, corresponding to 11 for each of the following 4 switch experiments: M9G to M9G, M9G to M9A, M9GA to M9G, and M9GA to M9A. For each simulation, the lag time is computed according to Enjalbert et al. (22):where X₀ is the total initial E. coli biomass, X₁ is the total E. coli biomass value at time t₁ (1.5 h as described previously [22]), and μ_max values are used as described by Enjalbert et al. (22).

Experimental data (values with standard deviations, when available) from Enjalbert et al. (22) were kindly provided by B. Enjalbert. The data from Varma and Palsson (18) were extracted from the original publication using WebPlotDigitizer (44).

The version of the modeling framework used to obtain the results presented in this manuscript (v1.1) is publicly available with instructions to install and run simulations at https://github.com/QTB-HHU/daphne_ecoli-diauxie. The development version is hosted on https://gitlab.com/asuccurro/dfba-ode-framework, and people interested in contributing can request access by contacting A. Succurro (corresponding author).