Research Article | Molecular Biology and Physiology

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

Antonella Succurro, Daniel Segrè, Oliver Ebenhöh
Sarah Glaven, Editor
DOI: 10.1128/mSystems.00230-18

ABSTRACT

Microbes have adapted to greatly variable environments in order to survive both short-term perturbations and permanent changes. A classical and yet still actively studied example of adaptation to dynamic environments is the diauxic shift of Escherichia coli, in which cells grow on glucose until its exhaustion and then transition to using previously secreted acetate. Here we tested different hypotheses concerning the nature of this transition by using dynamic metabolic modeling. To reach this goal, we developed an open source modeling framework integrating dynamic models (ordinary differential equation systems) with structural models (metabolic networks) which can take into account the behavior of multiple subpopulations and smooth flux transitions between time points. We used this framework to model the diauxic shift, first with a single E. coli model whose metabolic state represents the overall population average and then with a community of two subpopulations, each growing exclusively on one carbon source (glucose or acetate). After introduction of an environment-dependent transition function that determined the balance between subpopulations, our model generated predictions that are in strong agreement with published data. Our results thus support recent experimental evidence that diauxie, rather than a coordinated metabolic shift, would be the emergent pattern of individual cells differentiating for optimal growth on different substrates. This work offers a new perspective on the use of dynamic metabolic modeling to investigate population heterogeneity dynamics. The proposed approach can easily be applied to other biological systems composed of metabolically distinct, interconverting subpopulations and could be extended to include single-cell-level stochasticity.

IMPORTANCE Escherichia coli diauxie is a fundamental example of metabolic adaptation, a phenomenon that is not yet completely understood. Further insight into this process can be achieved by integrating experimental and computational modeling methods. We present a dynamic metabolic modeling approach that captures diauxie as an emergent property of subpopulation dynamics in E. coli monocultures. Without fine-tuning the parameters of the E. coli core metabolic model, we achieved good agreement with published data. Our results suggest that single-organism metabolic models can only approximate the average metabolic state of a population, therefore offering a new perspective on the use of such modeling approaches. The open source modeling framework that we provide can be applied to model general subpopulation systems in more-complex environments and can be extended to include single-cell-level stochasticity.

INTRODUCTION

In natural environments, microorganisms are exposed to high fluctuations of nutrient and micronutrient availability and have therefore evolved adaptation strategies, both short term (to respond to temporary perturbations) and long term (to increase evolutionary fitness) (1). We still lack a sound theoretical understanding of the mechanisms driving such strategies, but the recent technological advances in high-throughput experimental techniques pave the way to novel approaches that integrate experimental and theoretical biology (2). Theoretical ecology describes ecosystems in mathematical terms as dynamic organism-environment interactions (3). As in statistical physics, individual behaviors in an ensemble result in observable emergent patterns that can be modeled with mathematical equations (4). This is the case for the earliest models of population dynamics developed by Verhulst (5), Lotka (6), and Volterra (7) and for the pioneering work of Jacques Monod in modeling microbial growth (8). With the rising academic and industrial interest in the “microbiome,” systems biology approaches are becoming a new standard (9) and more methods for the mathematical modeling of microbial communities are being developed (10, 11).

In constraint-based stoichiometric modeling, the metabolic network model of an organism is reconstructed from its annotated genome and described mathematically as a stoichiometric matrix (S). After imposition of the steady-state assumption and introduction of thermodynamic and biological boundaries for the metabolic fluxes (v→ ), flux balance analysis (FBA) (12, 13) defines an optimization problem in order to identify one particular flux distribution in the solution space. As long as the objective function (which imposes further biological assumptions on the system) is linear in the fluxes, the optimization problem can be solved by linear programming (LP). FBA returns a unique solution for the objective function, but the metabolic flux distribution is generally not unique, especially in genome-scale metabolic network models (GEMs). On the basis of the hypothesis that metabolism has evolved to make efficient use of resources and minimize waste, two specific methods were developed to extend FBA: parsimonious FBA (pFBA) (14) and minimization of metabolic adjustment (MOMA) (15). In pFBA, a second LP is defined such that the value of the objective function is set to the FBA solution and the new objective is the minimization of the overall fluxes. MOMA was developed to simulate the response to the perturbation introduced by gene deletion and is based on the principle that the organism would readjust its metabolism to a minimally different configuration with respect to the wild-type optimum. Another extension of FBA, dynamic FBA (dFBA) (16), allows partial recovery of the dynamic information lost under the conditions of the steady-state assumption. In the static optimization approach (SOA) that underlies dFBA, time is divided into discrete intervals and a new FBA problem is solved at time ti after updating of the external conditions according to the FBA solution at time ti−1. Approaches to modeling of microbial communities with GEMs have been recently reviewed by Succurro and Ebenhöh (17).

FBA and dFBA have been applied to the study of one of the most basic examples of metabolic transitions: diauxie (16, 18, 19). Discovered in the model organism Escherichia coli in 1941 by Monod (8, 20), diauxie remains a topic of active research (21–23). Under aerobic conditions with glucose as the sole carbon source (and generally the preferred one), E. coli secretes acetate during growth, which it then consumes once the glucose is exhausted. The molecular mechanisms driving this transition are still not completely understood, but over the last few years the fundamental roles of stochasticity and population heterogeneity have been demonstrated experimentally (24), often with the support of mathematical models. Indeed, in unpredictable natural environments fluctuating conditions of nutrient availability and variable fitness landscapes, homogeneous populations are more likely to face extinction, and bet hedging provides a selective advantage (25). Single-cell studies have suggested that the observed biphasic growth possibly represents the effect of stochastic gene expression (21), eventually coregulated by memory mechanisms (26). Kotte et al. (27) systematically investigated bistability in a clonal E. coli population. After ruling out responsive switching as a homogeneous adaptation, their results strongly suggested that the heterogeneous adaptation that results in two coexisting phenotypes was driven by responsive diversification (where a single phenotype diversifies in response to environmental changes) rather than stochastic switching (where the two phenotypes would coexist from the beginning). Although stochastic mathematical models have been proposed to support those findings, metabolic modeling approaches are only considered suitable to describe homogeneous systems, with single-organism GEMs representing the average population metabolic state.

Varma and Palsson (18) performed the first dFBA of E. coli, with a single GEM growing aerobically first on glucose and then on the secreted acetate. Here we present a study of E. coli diauxic growth on these two carbon sources, with the bacterial population modeled either as having an average, unique metabolic state (standard FBA and dFBA approach) or as being composed of two E. coli populations adapted to one of the two carbon sources. We used a modeling approach that integrates ordinary differential equation (ODE) models with dFBA, extending methods typically applied to study the dynamics of multispecies communities to the investigation of emergent patterns from individual behavior in monocultures. We implemented three approaches: (i) we modeled a homogeneous and yet smooth shift, with a single E. coli GEM, by adapting the MOMA algorithm; (ii) we introduced the hypothesis of subpopulations growing on specific carbon sources and model population transition as a purely stochastic mechanism; and (iii) we introduced an environment-driven response. Our results suggest that diauxie, rather than being modeled as a coordinated metabolic shift, can be modeled as the emergent pattern resulting from subpopulations optimizing growth on different substrates in response to environmental changes. This is much in agreement with experimental evidence from, e.g., Kotte et al. (27) and offers a new perspective on the use of dynamic metabolic modeling to investigate population dynamics. The proposed approach can easily be transferred to studies of generic subpopulations or communities and ultimately can be expanded to investigate single-cell dynamics.

RESULTS

We ran simulations with an open source modeling framework developed to model ecosystem dynamics. The models are ODE systems solved with integrating routines that at each integration step solve an FBA problem. We first validated the E. coli GEM on the data from Varma and Palsson (18) (who reported the first dFBA of the glucose-acetate shift) and then used the calibrated model to reproduce the independent sets of experiments described by Enjalbert et al. (22) (who experimentally analyzed E. coli grown in aerobic batch systems with different concentrations of glucose and acetate). In the standard dFBA approach, a population is modeled with a unique GEM and fluxes instantaneously change to adapt to new environmental conditions. In reality, however, transcriptional changes and flux rerouting may cause delays which are not captured by existing algorithms. Furthermore, dFBA might predict metabolic states in which more carbon sources are simultaneously utilized, and it is not obvious that such an approach would correctly capture the complexity of a population diversifying into metabolically distinct subpopulations. Therefore, we modified the dFBA algorithm by taking advantage of optimization strategies previously developed for different biological issues and implemented novel concepts as well. In particular, we used either pFBA (14) or an adaptation of MOMA (15) to solve the FBA problem at each time step, replicating the standard dFBA approach or implementing a homogeneous and yet smooth shift, respectively. The MOMA algorithm was integrated into the dFBA routine by imposing the constraint that the solution of the FBA problem at time ti had to be minimally different from the solution at time ti−1. We tested three different hypotheses: (i) homogeneous, smooth population shift; (ii) stochastic population shift; and (iii) environment-driven subpopulation differentiation. We observed that dFBA performed with both pFBA and MOMA predicted abrupt transitions from acetate catabolism to acetate anabolism and that condition-specific parameterizations were necessary to reproduce the different data. We then modeled two E. coli subpopulations growing exclusively on glucose or acetate. For this, we extended the standard dFBA approach to include the process of population shifts. We tested whether purely stochastic switches (ii) or, rather, responsive diversification (iii) could capture the diauxic behavior by modeling the population transitions either with constant rates (ii) or with a heuristic function dependent on carbon source concentrations (iii). We observed that only model iii could reproduce data from different experiments with a unique set of parameters. We did not find significant improvements using MOMA rather than pFBA within the same metabolic state, so the simpler pFBA implementation was used in the subpopulation simulations where each model was fixed into one metabolic configuration. Further details of the modeling approach are provided in Materials and Methods.

E. coli diauxie modeled with a uniform population.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 R2 distance between in vivo and in silico data for biomass (pFBA R2 = 0.989; MOMA R2 = 0.982), glucose (pFBA R2 = 0.993; MOMA R2 = 0.993), and acetate (pFBA R2 = 0.277; MOMA R2 = 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).

View this table:
TABLE 1

Fixed parameters for all simulationsa

View this table:
TABLE 2

Parameters of the simulationsa

FIG 1
FIG 1

Diauxic growth of E. coli modeled as a uniform population under batch conditions. Simulation data (lines) are compared to data from Enjalbert et al. (22) (squares) as a function of time. Biomass data (blue, top subplots) and glucose and acetate data (red and yellow, bottom subplots) are shown. The flux distribution at each time step was obtained with pFBA (a) or MOMA (b). gDW, grams dry weight; C, concentration; Q, quantity.

FIG 2
FIG 2

Comparisons of experimental information on gene expression levels with simulated flux distributions. The top plots qualitatively represent the gene expression data from Enjalbert et al. (22). Flux solutions in the simulations for the reactions associated with the reported key genes are compared between two independent simulations with E. coli exponentially growing either on acetate or on glucose (a) and within the same simulation (b) (growth on glucose simulated with MOMA or pFBA) (Fig. 1) before and after the point of glucose exhaustion. Genes: acs, acetyl coenzyme A synthetase; pck, phosphoenolpyruvate carboxykinase; pps, phosphoenolpyruvate synthetase; fbp, fructose-1,6-biphosphatase; icl, isocitrate lyase; mls, malate synthase; pfkA, phosphofructokinase; pykF, pyruvate kinase; ppc, phosphoenolpyruvate carboxylase; icd, isocitrate dehydrogenase. Reactions (BiGG identifiers): ACKr, acetate kinase; PPCK, phosphoenolpyruvate carboxykinase; PPS, phosphoenolpyruvate synthase; FBP, fructose-bisphosphatase; ICL, isocitrate lyase; MALS, malate synthase; PFK, phosphofructokinase; PYK, pyruvate kinase; PPC, phosphoenolpyruvate carboxylase; ICDHyr, isocitrate dehydrogenase (NADP).

FIG S1

Diauxic growth of E. coli modeled as a uniform population in batch (a and e) and fedbatch (b and f) conditions. Simulation (connected dots; each dot represents a simulation step) data are compared to data (squares) from Varma and Palsson (18) as a function of time. Biomass (blue), glucose (red), and acetate (cyan) data are shown. The flux distribution at each time step was obtained with pFBA (a and b) or MOMA (e and f). Below each plot, left to right, are reported data corresponding to agreement between the biomass, glucose, and acetate experimental data with simulation results. Download FIG S1, PDF file, 0.4 MB.
Copyright © 2019 Succurro et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

FIG S2

(Left plots) Euclidean distance between the flux solution vectors at simulation steps i and i + 1 plotted against simulation time (lower x axis) and simulation steps (upper x axis, in red). (Right plots) Absolute distance between the flux data at simulation steps i and i + 1 for the reactions yielding the higher change, plotted against simulation time (lower x axis) and simulation steps (upper x axis, in red). Notice the nonlinear relation between simulation time and simulation step data, originating from the variable integrator time step choice. Download FIG S2, PDF file, 0.3 MB.
Copyright © 2019 Succurro et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

FIG S3

(Left plots) Euclidean distance between the flux solution vectors at simulation steps i and i + 1 plotted against simulation time (lower x axis) and simulation steps (upper x axis, in red). (Right plots) Absolute distance between the flux data at simulation steps i and i + 1 for the reactions yielding the higher change, plotted against simulation time (lower x axis) and simulation steps (upper x axis, in red). Notice the nonlinear relation between simulation time and simulation step data, originating from the variable integrator time step choice. Download FIG S3, PDF file, 0.1 MB.
Copyright © 2019 Succurro et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

FIG S4

Escher (28) map of the metabolic network of the E. coli core. The pathways highlighted in different colors are those that were found to differ between the time before and the time after glucose exhaustion (GE at −115 min and GE at +18 min, respectively) in the simulations. In red are highlighted the reactions found to carry higher fluxes after the switch to growth on acetate, in blue those carrying lower fluxes, and in gray those that did not change. The names of the reactions are color-coded in the same way but in accordance with the experimental results from Enjalbert et al. (22). Download FIG S4, EPS file, 2.2 MB.
Copyright © 2019 Succurro et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

FIG S5

Simulations of the three experiments from Enjalbert et al. (22). (Left column) Growth on 15 mM glucose. (Middle column) Growth on 15 mM glucose and fed-batch 4 mM acetate. (Right column) Growth on 15 mM glucose and fed-batch 32 mM acetate. (Top row) Uniform population model. (Middle row) Subpopulation model with stochastic shift (fixed-rate transition). (Bottom row) Subpopulation model with responsive shift (Hill function transition). Download FIG S5, PDF file, 0.3 MB.
Copyright © 2019 Succurro et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

E. coli diauxie modeled with a mixed population.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 (ECGl) model and the E. coli acetate (ECAc) 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% ECGl and 5% ECAc. 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 (R2 = 0.987 and R2 = 0.990, respectively) and glucose concentrations (R2 = 0.996 and R2 = 0.997, respectively), only the simulation that included the population transition realistically reproduced the acetate consumption levels (R2 = 0.336 and R2 = 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 ECGl and ECAc) 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) (R2 = 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) (R2 = 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% ECGl and 5% ECAc describes well the biomass dynamics (R2 = 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% ECGl and 25% ECAc reproduced the biomass measurements (R2 = 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% ECGl and 5% ECAc or 75% ECGl and 25% ECAc did not show strong differences, but further reducing the percentage of ECGl (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% ECGl and 25% ECAc population distribution for Fig. 4b is also justified by a difference in the initial experimental values for the biomass quantity (see Fig. S7).

FIG 3
FIG 3

Diauxic growth of E. coli modeled as a mixture of two E. coli populations, ECGl and ECAc, growing exclusively on glucose and acetate, respectively, without (a) or with (b) the possibility of shifting from one population to the other. Simulation data (lines) are compared to data from Enjalbert et al. (22) (squares) as a function of time. The upper plots show simulation results (obtained using pFBA) for ECGl and ECAc biomass data (light blue and aqua, respectively) and the observable E. coli biomass data (black line simulation, blue dots). The bottom plots show glucose and acetate (red and yellow, respectively).

FIG 4
FIG 4

Diauxic growth of E. coli modeled as a mixture of two E. coli populations, ECGl and ECAc, growing exclusively on glucose and acetate, respectively, with the possibility of shifting from one population to the other. (a) E. coli grows on 15 mM glucose; after the glucose was exhausted, the acetate concentration was kept at about 4 mM. (b) E. coli grows on 15 mM glucose and 32 mM acetate; after the glucose was exhausted, the acetate concentration was maintained at around the same concentration. The upper plots show simulation results (obtained using pFBA) for ECGl and ECAc biomass data (light blue and aqua lines, respectively) and the observable E. coli biomass data (black line simulation, blue dots) (data from Enjalbert et al. [22]). The bottom plots show simulation results for glucose and acetate (red and yellow lines, respectively).

TEXT S1

Simple mathematical model of two coexisting subpopulations in constant environments. Download Text S1, PDF file, 0.1 MB.
Copyright © 2019 Succurro et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

FIG S6

Simulation of two subpopulations (ECGlc and ECAc) growing on glucose only (a) or on glucose and (b) a low level of acetate or (c) a high level of acetate. Experimental data from Enjalbert et al. (22) are shown as blue circles. The different curves show simulation results with same parameters except for initial biomass ratios. Download FIG S6, PDF file, 0.2 MB.
Copyright © 2019 Succurro et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

FIG S7

Results from simulating ECGlc and ECAc growing on glucose and a high level of acetate (see main text). The initial biomass value was taken from time point “GE—4hr” (a) or “GE—5hr” (b) in the data from Enjalbert et al. (22), with the initial population ratio as reported on the plot. All other simulation parameters were the same and were as described in the manuscript. Download FIG S7, PDF file, 0.07 MB.
Copyright © 2019 Succurro et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

TABLE S1

Scan of the allowed values for ϕ and ψ under different conditions and respective values of log Γ. Under low-acetate conditions (minimal ψ), any value of ϕ is allowed, but we report a subset here. Download Table S1, PDF file, 0.02 MB.
Copyright © 2019 Succurro et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

Lag time for growth on acetate explained with population distribution.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 ECGl and ECAc 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% ECGl and 25% ECAc 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 ECGl and ECAc 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.

View this table:
TABLE 3

Percentage of ECGl biomass under M9G and M9GA conditions at indicated time points relative to glucose exhaustion

FIG 5
FIG 5

Simulation data (dark and light gray points) and experimental data (orange and yellow points) (data from Enjalbert et al. [22]) for the delay in the growth of daughter cultures before they reached maximal growth after the medium switch. Mother cultures are grown on either 15 mM glucose (M9G) (a) or 15 mM glucose and 32 mM acetate (M9GA) (b). Daughter cultures are reinoculated into fresh media with either 15 mM glucose (M9G; square and plus markers) or 45 mM acetate (M9A; diamond and cross markers). The simulation error bars are obtained by adjusting the initial population ratios (obtained by sampling the simulated mother cultures) by ±15%.

FIG S8

(a to d) Phase space of the lag time for typical daughter cultures from M9G mothers (a and b) and M9GA mothers (c and d) grown under M9G (a and c) or in M9A (b and d) conditions. When the glucose (acetate) value is zero, the transition ϕ (ψ) is in principle reduced to the noise component ϕ0 (ψ)0. (e to h) Effects of adjusting the transition parameters VMψ and VMϕ to values of 0.2, 0.5, and 0.8 and Kψ to values of 10 (e and f) and 30 (g and h) for the population ratio (top panel in the plots) and for the lag time (bottom panel in the plots) for daughter cultures sampled at different times with respect to the time of glucose exhaustion. Each point in panels e and g represents results from a daughter culture from the M9G mother culture, while each point in panels f and h represents results from the M9GA mother culture. Download FIG S8, PDF file, 0.2 MB.
Copyright © 2019 Succurro et al.

This content is distributed under the terms of the Creative Commons Attribution 4.0 International license.

DISCUSSION

We have investigated a fundamental example of metabolic adaptation, namely, the diauxic growth of E. coli on glucose and acetate, aiming to test whether a dynamic metabolic modeling approach can capture diauxie in monocultures of E. coli as the observable emergent result of individual (subpopulation) behavior. To this end, we first developed a modeling framework to integrate dynamic models (ODE systems) with structural models (metabolic networks) and then performed simulations to reproduce published experimental results in silico.

Avoiding fine-tuning of model 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.

Standard dFBA allows for abrupt metabolic readjustments.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 ti, the flux solution differs minimally from that at time ti-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.

Monocultures can be modeled as multisubpopulation systems to capture individual heterogeneity.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.

The lag phase corresponding to growth on different substrates can be explained by population distributions.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.

Subpopulations in the dynamic metabolic modeling approach.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.

Outlook.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.

MATERIALS AND METHODS

FBA methods.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 sij being the stoichiometric coefficients of metabolite i taking part in reaction j. FBA defines and solves the following LP problem:maximize z⃗ (1)subject to:  Sv⃗=0 (2)  l.b.j≤vj≤u.b.j (3)with l.b.j and u.b.j representing, respectively, a lower and upper bound on flux vj.

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 ti-1 is used to compute the MOMA solution at time ti 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.

Modeling framework integrating ODE and FBA.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 v⃗ as follows:dqidt=F(p⃗;v⃗,q⃗,t) (4)

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.

E. coli uniform population model.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 ECany 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:l.b.EX_Gl_eECany=−VMGl[Gl][Gl]+KMGl (5)l.b.EX_Ac_eECany=−VMAc[Ac][Ac]+KMAc (6)

The ODE system is defined as follows:dBMECanydt=vμECany·BMECany−δBMECany (7)dGldt=vEX_Gl_eECany·BMECany+ξfed-batch (8)dAcdt=vEX_Ac_eECany·BMECany (9)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.

E. coli mixed-population model.Two E. coli core models are loaded and defined as either a glucose consumer (ECGl) model or an acetate consumer (ECAc) model by switching off uptake of the other carbon source as follows:l.b.EX_Gl_eECGl=−VMGl[Gl][Gl]+KMGl (10)l.b.EX_Gl_eECAc=0 (11)l.b.EX_Ac_eECGl=0 (12)l.b.EX_Ac_eECAc=−VMAc[Ac][Ac]+KMAc (13)

The ODE system is defined as follows:dBMECGldt=(vμECGl−ψ−δ)·BMECGl+εϕBMECAc (14)dBMECAcdt=(vμECAc−ϕ−δ)·BMECAc+εψBMECGl (15)dGldt=vEX_Gl_eECGl·BMECGl+vEX_Gl_eECAc·BMECAc+ξfed-batch (16)dAcdt=vEX_Ac_eECGl·BMECGl+[vEX_Gl_eECAc+ζ·H(t−tx)]·BMECAc (17)where ζ · H(t − tx) is a heaviside function activated at time tx of glucose exhaustion in order to keep the acetate level constant, ψ and ϕ are functions that model the cellular shift from ECGl to ECAc and from ECAc to ECGl, 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:ψ([Ac])=ψ0+VMψ[Ac]n[Ac]n+KMψn (18)ϕ([Gl])=ϕ0+VMϕ[Gl]n[Gl]n+KMϕn (19) where they are constant transition rates for VMϕ=VMψ=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.

Switch experiment simulations.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 ECGl and ECAc used as the initial condition for the reinoculation simulations. The percentage of ECGl 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):tlag=t1−ln⁡(X1/X0)μmax (20)where X0 is the total initial E. coli biomass, X1 is the total E. coli biomass value at time t1 (1.5 h as described previously [22]), and μmax values are used as described by Enjalbert et al. (22).

Published experimental data.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).

Data availability.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).

ACKNOWLEDGMENTS

A.S., D.S., and O.E. initiated the project; A.S. designed the final project, developed the modeling framework, performed simulations, and wrote the initial manuscript; all of us contributed to the interpretation of the results and to the final manuscript. The authors are grateful to Brice Enjalbert for kindly providing the original experimental data used to compare with the simulation results and acknowledge Richard M. Twyman for assistance with manuscript editing.

A.S. and O.E. are supported by the Deutsche Forschungsgemeinschaft, Cluster of Excellence on Plant Sciences CEPLAS (EXC 1028). A.S. was supported also by funding from the European Commission Seventh Framework Marie Curie Initial Training Network project AccliPhot (grant agreement PITN-GA-2012-316427). D.S. acknowledges funding from the U.S. Department of Energy (DE-SC0012627), the NIH (5R01DE024468 and R01GM121950), the National Science Foundation (grants 1457695 and NSFOCE-BSF 1635070), MURI grant W911NF-12-1-0390, the Human Frontiers Science Program (RGP0020/2016), and the Boston University Interdisciplinary Biomedical Research Office. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

FOOTNOTES

    • Received October 16, 2018.
    • Accepted November 27, 2018.

This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license.

REFERENCES

  1. 1.↵
    1. van Boxtel C,
    2. van Heerden JH,
    3. Nordholt N,
    4. Schmidt P,
    5. Bruggeman FJ
    . 2017. Taking chances and making mistakes: non-genetic phenotypic heterogeneity and its consequences for surviving in dynamic environments. J R Soc Interface 14:20170141. doi:10.1098/rsif.2017.0141.
    OpenUrlCrossRefPubMed
  2. 2.↵
    1. Succurro A,
    2. Moejes FW,
    3. Ebenhöh O
    . 2017. A diverse community to study communities: integration of experiments and mathematical models to study microbial consortia. J Bacteriol 199:e00865-16. doi:10.1128/JB.00865-16.
    OpenUrlAbstract/FREE Full Text
  3. 3.↵
    1. Murray JD
    (ed). 2004. Mathematical biology of interdisciplinary applied mathematics, vol 17. Springer, New York, NY. http://link.springer.com/10.1007/b98868.
  4. 4.↵
    1. Hagstrom GI,
    2. Levin SA
    . 2017. Marine ecosystems as complex adaptive systems: emergent patterns, critical transitions, and public goods. Ecosystems 20:458–476. doi:10.1007/s10021-017-0114-3.
    OpenUrlCrossRef
  5. 5.↵
    1. Verhulst PF
    . 1838. Notice sur la loi que la population poursuit dans son accroissement. Corresp Mathématique Phys 10:113–121.
    OpenUrl
  6. 6.↵
    1. Lotka A
    . 1925. Elements of physical biology. Williams and Wilkins, Baltimore, MD.
  7. 7.↵
    1. Volterra V
    . 1926. Fluctuations in the abundance of a species considered mathematically. Nat 118:558–560. doi:10.1038/118558a0.
    OpenUrlCrossRefWeb of Science
  8. 8.↵
    1. Monod J
    . 1949. The growth of bacterial cultures. Annu Rev Microbiol 3:371–394. doi:10.1146/annurev.mi.03.100149.002103.
    OpenUrlCrossRefWeb of Science
  9. 9.↵
    1. Widder S,
    2. Allen RJ,
    3. Pfeiffer T,
    4. Curtis TP,
    5. Wiuf C,
    6. Sloan WT,
    7. Cordero OX,
    8. Brown SP,
    9. Momeni B,
    10. Shou W,
    11. Kettle H,
    12. Flint HJ,
    13. Haas AF,
    14. Laroche B,
    15. Kreft JU,
    16. Rainey PB,
    17. Freilich S,
    18. Schuster S,
    19. Milferstedt K,
    20. van der Meer JR,
    21. Groβkopf T,
    22. Huisman J,
    23. Free A,
    24. Picioreanu C,
    25. Quince C,
    26. Klapper I,
    27. Labarthe S,
    28. Smets BF,
    29. Wang H,
    30. Soyer OS
    . 2016. Challenges in microbial ecology: building predictive understanding of community function and dynamics. ISME J 10:2557–2568. doi:10.1038/ismej.2016.45.
    OpenUrlCrossRefPubMed
  10. 10.↵
    1. Song HS,
    2. Cannon W,
    3. Beliaev A,
    4. Konopka A
    . 2014. Mathematical modeling of microbial community dynamics: a methodological review. Process 2:711–752. doi:10.3390/pr2040711.
    OpenUrlCrossRef
  11. 11.↵
    1. Zomorrodi AR,
    2. Segrè D
    . 2016. Synthetic ecology of microbes: mathematical models and applications. J Mol Biol 428:837–861. doi:10.1016/j.jmb.2015.10.019.
    OpenUrlCrossRefPubMed
  12. 12.↵
    1. Fell DA,
    2. Small JR
    . 1986. Fat synthesis in adipose tissue. An examination of stoichiometric constraints. Biochem J 238:781–786. doi:10.1042/bj2380781.
    OpenUrlAbstract/FREE Full Text
  13. 13.↵
    1. Varma A,
    2. Palsson BO
    . 1993. Metabolic capabilities of Escherichia coli II. Optimal growth patterns. J Theor Biol 165:503–522. doi:10.1006/jtbi.1993.1203.
    OpenUrlCrossRef
  14. 14.↵
    1. Lewis NE,
    2. Hixson KK,
    3. Conrad TM,
    4. Lerman JA,
    5. Charusanti P,
    6. Polpitiya AD,
    7. Adkins JN,
    8. Schramm G,
    9. Purvine SO,
    10. Lopez-Ferrer D,
    11. Weitz KK,
    12. Eils R,
    13. König R,
    14. Smith RD,
    15. Palsson BØ
    . 2010. Omic data from evolved E. coli are consistent with computed optimal growth from genome-scale models. Mol Syst Biol 6:390. https://doi.org/10.1038/msb.2010.47.
    OpenUrlAbstract/FREE Full Text
  15. 15.↵
    1. Segrè D,
    2. Vitkup D,
    3. Church GM
    . 2002. Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci U S A 99:15112–15117. doi:10.1073/pnas.232349399.
    OpenUrlAbstract/FREE Full Text
  16. 16.↵
    1. Mahadevan R,
    2. Edwards JS,
    3. Doyle FJ
    . 2002. Dynamic flux balance analysis of diauxic growth in Escherichia coli. Biophys J 83:1331–1340. doi:10.1016/S0006-3495(02)73903-9.
    OpenUrlCrossRefPubMedWeb of Science
  17. 17.↵
    1. Succurro A,
    2. Ebenhöh O
    . 2018. Review and perspective on mathematical modeling of microbial ecosystems. Biochem Soc Trans 46:BST20170265. https://www.ncbi.nlm.nih.gov/pubmed/29540507.
    OpenUrl
  18. 18.↵
    1. Varma A,
    2. Palsson BO
    . 1994. Stoichiometric flux balance models quantitatively predict growth and metabolic by-product secretion in wild-type Escherichia coli W3110. Appl Environ Microbiol 60:3724–3731. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC201879/.
    OpenUrlAbstract/FREE Full Text
  19. 19.↵
    1. Harcombe WR,
    2. Delaney NF,
    3. Leiby N,
    4. Klitgord N,
    5. Marx CJ
    . 2013. The ability of flux balance analysis to predict evolution of central metabolism scales with the initial distance to the optimum. PLoS Comput Biol 9:e1003091. doi:10.1371/journal.pcbi.1003091.
    OpenUrlCrossRefPubMed
  20. 20.↵
    1. Monod J
    . 1942. Recherches sur la croissance des cultures bactèriennes. Hermann & cie, Paris, France.
  21. 21.↵
    1. Boulineau S,
    2. Tostevin F,
    3. Kiviet DJ,
    4. ten Wolde PR,
    5. Nghe P,
    6. Tans SJ
    . 2013. Single-cell dynamics reveals sustained growth during diauxic shifts. PLoS One 8:e61686. doi:10.1371/journal.pone.0061686.
    OpenUrlCrossRefPubMed
  22. 22.↵
    1. Enjalbert B,
    2. Cocaign-Bousquet M,
    3. Portais JC,
    4. Letisse F
    . 2015. Acetate exposure determines the diauxic behavior of Escherichia coli during the glucose-acetate transition. J Bacteriol 197:3173–3181. doi:10.1128/JB.00128-15.
    OpenUrlAbstract/FREE Full Text
  23. 23.↵
    1. Sandberg TE,
    2. Lloyd CJ,
    3. Palsson BO,
    4. Feist AM
    . 2017. Laboratory evolution to alternating substrate environments yields distinct phenotypic and genetic adaptive strategies. Appl Environ Microbiol 83:e00410-17. https://doi.org/10.1128/AEM.00410-17.
  24. 24.↵
    1. Siegal ML
    . 2015. Shifting sugars and shifting paradigms. PLoS Biol 13:e1002068. doi:10.1371/journal.pbio.1002068.
    OpenUrlCrossRef
  25. 25.↵
    1. Solopova A,
    2. van Gestel J,
    3. Weissing FJ,
    4. Bachmann H,
    5. Teusink B,
    6. Kok J,
    7. Kuipers OP
    . 2014. Bet-hedging during bacterial diauxic shift. Proc Natl Acad Sci U S A 111:7427–7432. doi:10.1073/pnas.1320063111.
    OpenUrlAbstract/FREE Full Text
  26. 26.↵
    1. Lambert G,
    2. Kussel E
    . 2014. Memory and fitness optimization of bacteria under fluctuating environments. PLoS Genet 10:e1004556. doi:10.1371/journal.pgen.1004556.
    OpenUrlCrossRefPubMed
  27. 27.↵
    1. Kotte O,
    2. Volkmer B,
    3. Radzikowski JL,
    4. Heinemann M
    . 2014. Phenotypic bistability in Escherichia coli’s central carbon metabolism. Mol Syst Biol 10:736–736. doi:10.15252/msb.20135022.
    OpenUrlAbstract/FREE Full Text
  28. 28.↵
    1. King ZA,
    2. Dräger A,
    3. Ebrahim A,
    4. Sonnenschein N,
    5. Lewis NE,
    6. Palsson BO
    . 2015. Escher: a Web application for building, sharing, and embedding data-rich visualizations of biological pathways. PLoS Comput Biol 11:e1004321. doi:10.1371/journal.pcbi.1004321.
    OpenUrlCrossRefPubMed
  29. 29.↵
    1. Beg QK,
    2. Vazquez A,
    3. Ernst J,
    4. de Menezes MA,
    5. Bar-Joseph Z,
    6. Barabasi A-L,
    7. Oltvai ZN
    . 2007. Intracellular crowding defines the mode and sequence of substrate uptake by Escherichia coli and constrains its metabolic activity. Proc Natl Acad Sci U S A 104:12663–12668. doi:10.1073/pnas.0609845104.
    OpenUrlAbstract/FREE Full Text
  30. 30.↵
    1. van Heerden JH,
    2. Wortel MT,
    3. Bruggeman FJ,
    4. Heijnen JJ,
    5. Bollen YJM,
    6. Planque R,
    7. Hulshof J,
    8. O'Toole TG,
    9. Wahl SA,
    10. Teusink B
    . 2014. Lost in transition: start-up of glycolysis yields subpopulations of nongrowing cells. Sci 343:1245114. doi:10.1126/science.1245114.
    OpenUrlAbstract/FREE Full Text
  31. 31.↵
    1. Nath K,
    2. Koch A
    . 1971. Protein degradation in Escherichia coli. II. Strain differences in degradation of protein and nucleic acid resulting from starvation. J Biol Chem 246:6956–6967. http://www.jbc.org/content/246/22/6956.abstract?ijkey=1c8784a9a1a98d492436344cd3dc11ad63e5de75&keytype2=tf_ipsecsha.
    OpenUrlAbstract/FREE Full Text
  32. 32.↵
    1. Borek E,
    2. Ponticorvo L,
    3. Rittenberg D
    . 1958. Protein turnover in micro-organisms. Proc Natl Acad Sci U S A 44:369–374. http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC335432. doi:10.1073/pnas.44.5.369.
    OpenUrlFREE Full Text
  33. 33.↵
    1. Pine MJ
    . 1965. Heterogeneity of protein turnover in Escherichia coli. Biochim Biophys Acta 104:439–456. doi:10.1016/0304-4165(65)90349-1.
    OpenUrlCrossRefPubMedWeb of Science
  34. 34.↵
    1. DeGennaro CM,
    2. Savir Y,
    3. Springer M
    . 2016. Identifying metabolic subpopulations from population level mass spectrometry. PLoS One 11:e0151659. doi:10.1371/journal.pone.0151659.
    OpenUrlCrossRef
  35. 35.↵
    1. Chu D
    . 2017. Limited by sensing - a minimal stochastic model of the lag-phase during diauxic growth. J Theor Biol 414:137–146. doi:10.1016/j.jtbi.2016.10.019.
    OpenUrlCrossRef
  36. 36.↵
    1. Gomez JA,
    2. Höffner K,
    3. Barton PI
    . 2014. DFBAlab: a fast and reliable MATLAB code for dynamic flux balance analysis. BMC Bioinformatics 15:409. doi:10.1186/s12859-014-0409-8.
    OpenUrlCrossRefPubMed
  37. 37.↵
    1. Zomorrodi AR,
    2. Islam MM,
    3. Maranas CD
    . 2014. d-OptCom: dynamic multi-level and multi-objective metabolic modeling of microbial communities. ACS Synth Biol 3:247–257. doi:10.1021/sb4001307.
    OpenUrlCrossRefPubMedWeb of Science
  38. 38.↵
    1. Ebrahim A,
    2. Lerman JA,
    3. Palsson BO,
    4. Hyduke DR
    . 2013. COBRApy: COnstraints-Based Reconstruction and Analysis for Python. BMC Syst Biol 7:74. doi:10.1186/1752-0509-7-74.
    OpenUrlCrossRefPubMed
  39. 39.↵
    1. van der Walt S,
    2. Colbert SC,
    3. Varoquaux G
    . 2011. The NumPy array: a structure for efficient numerical computation. Comput Sci Eng 13:22–30. http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=5725236. doi:10.1109/MCSE.2011.37.
    OpenUrlCrossRef
  40. 40.↵
    1. Jones E,
    2. Oliphant T,
    3. Peterson P
    . 2001. SciPy: Open source scientific tools for Python. http://www.scipy.org/.
  41. 41.↵
    1. Hunter JD
    . 2007. Matplotlib: a 2D graphics environment. Comput Sci Eng 9:90–95. http://scitation.aip.org/content/aip/journal/cise/9/3/10.1109/MCSE.2007.55. doi:10.1109/MCSE.2007.55.
    OpenUrlCrossRef
  42. 42.↵
    1. McKinney W
    . 2011. pandas: a foundational python library for data analysis and statistics. http://pandas.sourceforge.net/.
  43. 43.↵
    1. Orth JD,
    2. Palsson BØ,
    3. Fleming RMT
    . 2010. Reconstruction and use of microbial metabolic networks: the core Escherichia coli metabolic model as an educational guide. EcoSal Plus 4. http://www.asmscience.org/content/journal/ecosalplus/10.1128/ecosalplus.10.2.1. doi:10.1128/ecosalplus.10.2.1.
    OpenUrlCrossRef
  44. 44.↵
    1. Rohatgi A
    . 2017. WebPlotDigitizer v4.0. https://automeris.io/WebPlotDigitizer. Accessed 10 January 2018.
  45. 45.
    1. Saint-Ruf C,
    2. Taddei F,
    3. Matic I
    . 2004. Stress and survival of aging Escherichia coli rpoS colonies. Genetics 168:541–546. http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC1448099. doi:10.1534/genetics.104.028704.
    OpenUrlAbstract/FREE Full Text
  46. 46.
    1. Gosset G
    . 2005. Improvement of Escherichia coli production strains by modification of the phosphoenolpyruvate:sugar phosphotransferase system. Microb Cell Fact 4:14. http://www.microbialcellfactories.com/content/4/1/14. doi:10.1186/1475-2859-4-14.
    OpenUrlCrossRefPubMed
  47. 47.
    1. Harcombe WR,
    2. Riehl WJ,
    3. Dukovski I,
    4. Granger BR,
    5. Betts A,
    6. Lang AH,
    7. Bonilla G,
    8. Kar A,
    9. Leiby N,
    10. Mehta P,
    11. Marx CJ,
    12. Segrè D
    . 2014. Metabolic resource allocation in individual microbes determines ecosystem interactions and spatial dynamics. Cell Rep 7:1104–1115. http://www.cell.com/article/S2211124714002800/fulltext. doi:10.1016/j.celrep.2014.03.070.
    OpenUrlCrossRefPubMed
View Abstract
Back to top
Download PDF
Citation Tools
Print
Alerts
Email

KEYWORDS

diauxic growth
metabolic network modeling
microbial communities
population heterogeneity

Related Articles

Cited By...