1. Computational and Systems Biology
Download icon

Calibration and analysis of genome-based models for microbial ecology

  1. Stilianos Louca Is a corresponding author
  2. Michael Doebeli
  1. University of British Columbia, Canada
Tools and Resources
Cited
12
Views
1,358
Comments
0
Cite as: eLife 2015;4:e08208 doi: 10.7554/eLife.08208

Abstract

Microbial ecosystem modeling is complicated by the large number of unknown parameters and the lack of appropriate calibration tools. Here we present a novel computational framework for modeling microbial ecosystems, which combines genome-based model construction with statistical analysis and calibration to experimental data. Using this framework, we examined the dynamics of a community of Escherichia coli strains that emerged in laboratory evolution experiments, during which an ancestral strain diversified into two coexisting ecotypes. We constructed a microbial community model comprising the ancestral and the evolved strains, which we calibrated using separate monoculture experiments. Simulations reproduced the successional dynamics in the evolution experiments, and pathway activation patterns observed in microarray transcript profiles. Our approach yielded detailed insights into the metabolic processes that drove bacterial diversification, involving acetate cross-feeding and competition for organic carbon and oxygen. Our framework provides a missing link towards a data-driven mechanistic microbial ecology.

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

eLife digest

Microbes like bacteria and yeast play important roles in the environment, human health and even some industrial processes. However, it is difficult to understand the roles of microbes in these situations because many different types of microbes often live together in complex communities. Some of the microbes may compete with each other for resources like oxygen or sugar. Others may rely on one another for survival. For example, one microbe may feed on molecules that are released as waste from another microbe.

To better understand these microbial communities, we first need to understand the processes by which each microbe uses nutrients and releases waste molecules that influence other microbes. Researchers have used a technique called ‘genome sequencing’ to reconstruct the networks of genes and chemical reactions that are involved in these processes, and to build computer models of microbial communities in different environments.

However, the existing models can be labor intensive and do not allow researchers to easily use statistics to analyse them. To address this problem, Louca and Doebeli created a new computer model with built-in statistical tools that accurately predicts the interactions in communities that contain multiple strains of a bacterium called Escherichia coli. First, Louca and Doebeli grew a single strain of E. coli in the laboratory for many generations, which led to the evolution of the bacteria so that two new strains emerged. One of the new strains was more efficient at using sugar as a food source than the other and sometimes released a molecule called acetate. The other new strain became more efficient at using this acetate.

Next, Louca and Doebeli used data that had been collected for each individual strain, to test whether the model could recreate the way that the new strains had evolved together. The model accurately predicted that the two new strains would gradually replace the original strain. The strain that was more efficient at using sugar emerged first, which led to extra acetate being available for the other new strain that became more efficient at using acetate.

Louca and Doebeli's findings demonstrate for the first time that data collected for individual microbes can be used to explain the dynamics and evolution of small communities of microbes using mathematical models. The next step is to test this approach on larger communities in the environment.

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

Introduction

Metabolic interactions are an emergent property of microbial communities (Morris et al., 2013; Chiu et al., 2014). Even the simplest life forms can only be understood in terms of biological consortia characterized by shared metabolic pathways and distributed biosynthetic capacities (Klitgord and Segrè, 2010; McCutcheon and Moran, 2012; Husnik et al., 2013). For example, glucose catabolism to carbon dioxide or methane is a multi-step process often involving several organisms that indirectly exchange intermediate products through their environment (Stams, 1994). Microbial communities are thus complex systems comprising several interacting components that cannot be fully understood in isolation. In fact, metabolic interdependencies between organisms are at least partially responsible for our current inability to culture the great majority of prokaryotes (Schink and Stams, 2006). Understanding the emergent dynamics of microbial communities is crucial to harnessing these multicomponent assemblages and using synthetic ecology for medical, environmental and industrial purposes (Brenner et al., 2008).

Genome sequencing has enabled the reconstruction of full-scale cell-metabolic networks (Henry et al., 2010), which have provided a firm basis for understanding individual cell metabolism (Varma and Palsson, 1994; Duarte et al., 2004; Klitgord and Segrè, 2010). Recent work indicates that multiple cell models can be combined to understand microbial community metabolism and population dynamics (Stolyar et al., 2007; Klitgord and Segrè, 2010; Zengler and Palsson, 2012; Chiu et al., 2014; Harcombe et al., 2014). These approaches assume knowledge of all model parameters such as stoichiometric coefficients, maintenance energy requirements or extracellular transport kinetics, a requirement that is rarely met in practice (Feist et al., 2008; Harcombe et al., 2014). Experiments and monitoring of environmental samples could provide valuable data to calibrate microbial community models, for example, via statistical parameter estimation, but appropriate tools are lacking. So far, the standard approach has been to obtain each parameter through laborious specific measurements or from the available literature, or to manually adjust parameters to match observations (Mahadevan et al., 2002; Chiu et al., 2014; Harcombe et al., 2014). Furthermore, statistical model evaluation and sensitivity analysis is typically performed using ad hoc code, thus increasing the effort required for the construction of any new model. Consequently, the experimental validation of genome-based microbial community models and their application to biological questions are rare (Meadows et al., 2010; Harcombe et al., 2014).

We have developed MCM (Microbial Community Modeler), a mathematical framework and computational tool that unifies model construction with statistical evaluation, sensitivity analysis and parameter calibration. MCM is designed for modeling multi-species microbial communities, in which the metabolism and growth of individual cell species is predicted using genome-based metabolic models. Cells in the community interact in a dynamical environment in which metabolite concentrations and other environmental variables influence, and are influenced by, microbial metabolism. Unknown model parameters can be automatically calibrated (fitted) using experimental data such as cell densities, nutrient concentrations or rate measurements. To demonstrate the potential of MCM, we modeled a bacterial community that has emerged from in vitro evolution experiments, during which an ancestral strain repeatedly diversified into two distinct ecotypes. Experiments with microbes have an established tradition as model systems for understanding ecological and evolutionary processes (Elena and Lenski, 2003; Kassen and Rainey, 2004). We show that the predictions derived from MCM are in very good agreement with the outcomes of several monoculture and coculture experiments. While the experimental results described below have been found over the course of several years (Friesen et al., 2004; Spencer et al., 2007; Le Gac et al., 2008; Herron and Doebeli, 2013), it is only now that a mechanistic model has managed to unify them in a clear, unambiguous and synergistic manner. The analysis presented here thus provides strong credence to a large body of experimental work that was done in our lab over the course of roughly a decade.

Model

In MCM, a microbial community model is a set of differential equations for the population densities of the cell species comprising the community and of the ambient concentrations of utilized nutrients (metabolites), coupled to optimization problems for the cell-specific rates of reactions involving these metabolites. Each cell is characterized by its metabolic potential, that is, the genetically determined subset of reactions it can catalyze, as well as any available metabolite transport mechanisms. The reaction rates and metabolite exchange rates (i.e. the metabolism) of each cell are assumed to depend on its metabolic potential as well as on the current environmental conditions, such as metabolite concentrations. Through their metabolism, in turn, cells act as sinks and sources of metabolites in the environment. Additional metabolite fluxes, such as oxygen diffusion from the atmosphere into the growth medium of a modeled bacterial culture, can be included in the model.

At any point in time, individual cell metabolism is determined using flux balance analysis (FBA) (Orth et al., 2010), a widely used framework in cell-metabolic modeling (Varma and Palsson, 1994; Duarte et al., 2004; Klitgord and Segrè, 2010; Freilich et al., 2011; Chiu et al., 2014). In FBA, cell metabolism is assumed to be regulated in such a way that the rate of biosynthesis is maximized (Varma and Palsson, 1994; Feist and Palsson, 2010). The chemical state of cells is assumed to be steady, leading to stoichiometric constraints that need to be satisfied for any particular combination of intracellular reaction rates. Reaction rates, on the other hand, are limited due to finite enzyme capacities. Metabolite uptake/export rates can also be limited due to finite diffusion rates or limited transmembrane transporter efficiency. For example, uptake rates can be Monod-like functions of substrate concentrations (Mahadevan et al., 2002; Harcombe et al., 2014). Taken together, cell-metabolic potential, stoichiometric consistency, reaction rate limits and transport rate limits define the constraints of a linear optimization problem for each cell species at each point in time. The optimized biosynthesis rate is translated into a cell production rate by dividing by the cell's mass, thus defining the species' population growth (Figure 1).

Framework used by MCM.

(A) Conceptual framework used by MCM. Cells (colored shapes) optimize their metabolism for maximal growth and influence their environment via metabolite exchange (small colored arrows). Additional external fluxes can also affect the environment (large grey arrows). The environment, in turn, influences each cell's metabolism. (B) Computational framework used by MCM. Each iteration consists of four steps: flux balance analysis (FBA) is used to translate cell-metabolic potentials and environmental conditions (1) into a linear optimization problem for the growth rate of each cell species (2). The set of possible reaction rates corresponds to a polytope in high-dimensional space. Solving the optimization problems (3) yields predictions on microbial metabolite exchange rates (4). Metabolic fluxes and cell growth rates are used to predict metabolite and cell concentrations in the next iteration (1).

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

The central assumption of individual cells maximizing biosynthesis, subject to environmental and physiological constraints, is rooted in the idea that evolution has shaped regulatory mechanisms of unicellular organisms in such a way that they strive for maximum growth whenever possible. Biosynthesis has been experimentally verified as an objective for Saccharomyces cerevisiae and E. coli (Burgard and Maranas, 2003; Gianchandani et al., 2008; Harcombe et al., 2013). The assumption of maximized biosynthesis is less valid for genetically engineered organisms or those exposed to environments that are radically different from the environments that shaped their evolution (Segrè et al., 2002). Despite its limitations, FBA has greatly contributed to the understanding of several genome-scale metabolic networks and metabolic interactions between cells (Stolyar et al., 2007; Klitgord and Segrè, 2010; Orth et al., 2010; Freilich et al., 2011; Chiu et al., 2014; Harcombe et al., 2014). One advantage of FBA models over full biochemical cell models is their independence of intracellular kinetics and gene regulation, which limits the number of required parameters to stoichiometric coefficients and uptake kinetics.

The combination of FBA with a varying environmental metabolite pool, as implemented by MCM, is known as dynamic flux balance analysis (DFBA) (Mahadevan et al., 2002; Chiu et al., 2014; Harcombe et al., 2014). In contrast to conventional FBA, DFBA models are dynamical because cell densities and environmental metabolite concentrations both change with time, and the rate of change of each cell density and metabolite concentration depends on the current cell densities and metabolite concentrations (Mahadevan et al., 2002; Harcombe et al., 2014). Because metabolites can be depleted or produced by several cell species, the environmental metabolite pool mediates the metabolic interactions between cells (Schink and Stams, 2006). For example, oxygen uptake rates might depend on environmental oxygen concentrations, which in turn are reduced by cellular respiration. Similarly, cells might excrete acetate as a byproduct of glucose catabolism, which then becomes available to other cells. The metabolic optimization of individual cells striving for maximal growth, while modifying their environment, leads to non-trivial community dynamics that can include competition, cooperation and exploitation. The cell-centric nature of DFBA differs fundamentally from other flux balance analyses of microbial communities that assume an optimization of a community-wide objective such as total biomass synthesis (Stolyar et al., 2007; Klitgord and Segrè, 2011; Zomorrodi and Maranas, 2012). Such an assumption is at least questionable from an evolutionary perspective and likely not appropriate for communities comprising several species (Mitri and Foster, 2013).

Recent work suggests that DFBA is a promising approach to microbial ecological modeling (Meadows et al., 2010; Chiu et al., 2014; Harcombe et al., 2014). For example, Harcombe et al. (2014) designed a computational tool (COMETS) based on DFBA, which was able to accurately predict equilibrium compositions of mixed bacterial cultures grown on petri dishes. However, COMETS offers limited model versatility in terms of uptake and reaction kinetics and only has few environmental feedback mechanisms (namely, varying extracellular metabolite concentrations). Furthermore, it assumes complete knowledge of all required model parameters and provides no generic statistical model analysis. Hence, while COMETS sets an important precedent, considerable work is still needed to make DFBA a practical approach in microbial ecosystem modeling. MCM extends Harcombe et al.'s framework to more versatile microbial ecological models that include arbitrary reaction kinetics (e.g., subject to product-inhibition) as well as dynamical environmental variables (e.g., pH) that influence, and are influenced by, microbial metabolism. In addition, MCM supports cell models in which internal molecules act as dynamical constraints that further restrict the FBA solution space, for example to account for post-transcriptional regulation or delays in enzyme synthesis (Blazier and Papin, 2012). These so called regulatory FBA models have been shown to improve the fidelity of conventional FBA models for E. coli and S. cerevisiae (Covert et al., 2001; Covert and Palsson, 2002; Covert et al., 2004; Herrgård et al., 2006), however their application to microbial communities remains untested. MCM can statistically evaluate models against data, analyze their sensitivity to varying parameters (Cariboni et al., 2007), and estimate the uncertainty of model predictions in the face of stochasticity (Hammersley and Handscomb, 1964). Perhaps most importantly, MCM can automatically calibrate unknown model parameters to data, for example obtained from monoculture experiments (as demonstrated below), from bioreactor experiments involving multiple species (Louca and Doebeli, 2015) or from environmental samples of unculturable communities (Figure 2; see the ‘Materials and methods’ and the Supplement for details). MCM can thus be used to understand the dynamics of realistic microbial ecosystems, ranging from the soil or groundwater to mixed laboratory cultures and bioreactors.

Overview of MCM's working principle and functionalities: A microbial community model is specified using human-readable configuration files in terms of metabolites, reactions, the metabolic potential of cell species and any additional environmental variables.

Models with multiple ecosystem compartments are also possible. A script with MCM commands controls the analysis of the model and, if needed, its calibration using experimental data. The calibrated model can also be used to create new, more complex models (as exemplified in this article).

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

Results and discussion

Successional dynamics of a microbial community

In a series of laboratory evolution experiments with E. coli (strain B REL606; Yoon et al., 2012) in glucose-acetate supplemented medium, two metabolically distinct strains consistently evolved from the ancestral (A) strain (Le Gac et al., 2008; Spencer et al., 2008; Herron and Doebeli, 2013). When grown in monoculture with the same medium composition, all three strains exhibit diauxic growth curves with a fast glucose-driven growth phase followed by slower growth on acetate. However, the three strains differ in their efficiencies to catabolize glucose and acetate: Strain SS (slow switcher) is a better glucose utilizer when compared to strain A, and the depletion of glucose only leads to a slow switch to acetate consumption. On the other hand, the FS (fast switcher) strain has evolved to be a better acetate utilizer, initiating acetate consumption at higher remnant glucose concentrations than strains A and SS. This acetate specialization is based on a tradeoff in the citric acid cycle and comes at the cost of being a less competitive glucose consumer.

Replicated serial dilution experiments starting with strain A monocultures have shown a consistent phenotypic diversification, involving an initial invasion of the SS phenotype and a subsequent invasion of the FS phenotype, leading to the eventual extinction or near-extinction of the ancestor and the stable coexistence of the SS and FS phenotypes (Figure 3) (Le Gac et al., 2008; Spencer et al., 2008; Tyerman et al., 2008; Herron and Doebeli, 2013). Genome sequencing revealed that this metabolic diversification can be attributed to point-mutations in genes linked to glucose and acetate uptake kinetics and metabolism (Herron and Doebeli, 2013). The successional dynamics of the three phenotypes are thus likely driven by adaptations to a changing metabolic niche space, defined by fluctuating glucose, acetate and, potentially, oxygen availabilities (Le Gac et al., 2008; Tyerman et al., 2008; Herron and Doebeli, 2013). An understanding of the underlying ecological processes would shed light on the ecology and evolution of natural microbial communities with shared catabolic pathways.

Estimated relative cell densities of the A, SS and FS types during three replicated evolution experiments by Herron and Doebeli (2013, Figure 2), starting with the same ancestral E. coli strain.

Within each of the three experiments (AC), the illustrated SS or FS lineage comprises several strains with varyingly pronounced SS or FS phenotypes, respectively. Cell generations were translated to days by assuming an average of 6.7 generations per day (Herron and Doebeli, 2013).

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

To mechanistically explain the observed community dynamics, we used MCM to construct a model comprising the ancestral and the two evolved E. coli types. By keeping track of pathway activation, cell densities, metabolic fluxes and nutrient concentrations, we gained detailed insight into the processes driving the successional dynamics of metabolic diversification.

Experimental calibration

Based on a published cell-metabolic template for the ancestral E. coli strain comprising over 2000 reactions (Yoon et al., 2012), we first constructed three separate cell models for the phenotypes A, SS and FS, respectively. In these preliminary models, cells grew on a substrate pool that resembled previous batch-fed monoculture experiments with glucose-acetate supplemented minimal medium (Le Gac et al., 2008). Cell-specific oxygen, acetate and glucose uptake rate limits were Monod-like functions of substrate concentrations (Emerson and Hedges, 2008; Millero, 2013). We calibrated several physiological parameters for each cell type to measured chemical concentration and cell density profiles, using least squares fitting (Figure 4). MCM automatically calibrates free parameters to data through an optimization algorithm that involves step-wise exploration of parameter space and repeated simulations (see ‘Materials and methods’ and Supplementary Material).

Calibration of E. coli cell models.

Continuous curves: Time course of cell densities, glucose concentration, acetate concentration and oxygen concentration (columns 1–4, respectively) predicted by MCM for monocultures of strain A, SS and FS (rows 1–3, respectively) grown on glucose-acetate medium. Points are data used for model calibration, and were obtained from analogous monoculture growth experiments (Le Gac et al., 2008). Oxygen data were not available for strain A.

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

We then constructed the microbial community (MC) model by combining the three calibrated cell models into a community growing in a common substrate pool. The environmental context resembles Herron & Doebeli's evolution experiments (Herron and Doebeli, 2013). In particular, the model includes realistic oxygen depletion-repletion dynamics (Gupta and Rao, 2003), glucose and acetate depletion by microbial consumption, as well as daily dilutions into fresh glucose-acetate supplemented medium at a factor 1:100. The microbial community initially consists mostly of type A (1010 cells/l), while both SS as well as FS cells are assumed to be rare (1 cell/l). Because the model is deterministic, the invasion or extinction of each type only depends on its growth rate in a possibly changing environment, but not on random mutation events, nor on demographic stochastic fluctuations.

Predicting microbial community dynamics

Simulations of the MC model reproduced the successional dynamics observed in Herron & Doebeli's experiments: An initial replacement of the ancestor by the SS type is followed by an invasion of the FS type, leading to the eventual coexistence of the SS and FS types and extinction of the ancestral strain (Figure 5A). Interestingly, FS can also invade in the absence of SS, however invasion occurs much slower and FS reaches lower densities than in the presence of SS (Figure 5—figure supplement 1). This is consistent with an early presence of the FS lineage at low densities in the evolution experiments (Figure 3), indicating that some of the first FS mutations already confer a slight advantage over the ancestor when FS is rare (Herron and Doebeli, 2013).

Figure 5 with 2 supplements see all
Dynamics of the E. coli microbial community model.

(A) Relative cell densities of the A, SS and FS types over time. (B) Acetate concentration over time. (C), (D) and (E): SS and FS cell densities, relative cell densities and growth rates over time, respectively, during stable coexistence. (F), (G) and (H): Cell-specific glucose, acetate and oxygen uptake rates over time, respectively. Negative values correspond to export. (I), (J) and (K): Glucose, acetate and oxygen concentrations over time, respectively. Diurnal fluctuations in all figures are due to daily dilutions into fresh medium. Tics on the time axes in (CK) mark points of dilution.

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

Time series of acetate concentrations (Figure 5B) link the observed successional dynamics of the three types to a gradually changing metabolic niche space: The replacement of type A by the more efficient glucose specialist SS leads to an accumulation of acetate and facilitates the invasion of the FS type. The specialization of the SS and FS types on glucose and acetate, respectively (Figure 6A), enables their long-term coexistence on glucose-acetate enriched medium through frequency dependent competition (Friesen et al., 2004; Le Gac et al., 2008; Herron and Doebeli, 2013). In fact, cell-specific acetate exchange rates reveal that the SS type temporarily excretes acetate during short intervals, which is concurrently and subsequently consumed by the FS type (Figure 5G). This periodic acetate cross-feeding is an evolutionarily emergent property of the microbial community (Treves et al., 1998). The temporary production of acetate by the SS type is consistent with previous SS-FS coculture experiments, which revealed slightly increased acetate concentrations towards the end of the SS exponential growth phase (Spencer et al., 2007). An evolved increase of acetate excretion by E. coli in glucose minimal medium has also been reported by Harcombe et al. (2013).

Metabolic differentiation of the A, SS and FS types.

(A) Predicted cell-specific net metabolite uptake rates in coculture. (B) Predicted cell-specific reaction rates in coculture, for acs (acetyl-CoA synthesis), ack (acetate synthesis), pta (acetyl phosphate synthesis), ppc (oxaloacetate synthesis from phosphoenolpyruvate), pdh (decarboxylation of pyruvate to acetyl-CoA) and pyk (pyruvate synthesis from phosphoenolpyruvate). Rates in (A) and (B) are averaged over all time points within the first 100 days of evolution, with reversed reactions or net metabolite export represented by negative rates. (C) and (D): Simplified model subset of E. coli acetate and glucose metabolism, showing pathway activations in type SS (C) and FS (D) relative to type A during exponential growth in monoculture. Non-bracketed numeric values are ratios of predicted fluxes in the evolved types over fluxes in type A. Bracketed values are ratios of mRNA harvested from monoculture experiments by Le Gac et al. (2008), for comparison. A ratio of 0/0 indicates zero flux in both the evolved and ancestral type, a ratio of 1 corresponds to an unchanged flux or mRNA, a ratio of 0 corresponds to complete deactivation in the evolved type. Darker arrows indicate increased predicted fluxes in the evolved type. Flux predictions correspond to the time points of mRNA measurements, that is, 3.5 hr after dilution for SS and 4 hr after dilution for A and SS (Le Gac et al., 2008).

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

It should be noted that cell metabolism depends on substrate concentrations and is subject to strong temporal variation. In particular, acetate excretion by SS cells correlates strongly with oxygen limitation (Figure 5G,K). The excretion of acetate by E. coli as a byproduct of oxygen-limited glucose catabolism has been observed experimentally and explained using flux balance analysis (Mahadevan et al., 2002). In the absence of oxygen limitation, complete aerobic glucose catabolism to carbon dioxide is preferred over incomplete glucose catabolism with acetate excretion. On the other hand, oxygen limitation leads to an energetic tradeoff between complete glucose catabolism and efficient oxygen utilization, resulting in the excretion of acetate.

Furthermore, the depletion of oxygen during cell growth makes oxygen a temporary limiting resource for all cells (Figure 5K). Shortly after dilution into fresh medium, the exponential growth of the SS type on glucose leads to a rapid drop of oxygen to nanomolar concentrations. Despite oxygen diffusion into the medium, oxygen remains at sub-saturation levels for several more hours because the slow-growing acetate-consuming FS cells still consume oxygen after the growth of SS cells has halted. Differences in SS and FS growth rates (Figure 5C,E) thus mitigate competition for oxygen through temporal niche separation. Hence, oxygen likely plays an important role in the metabolic diversification, as previously hypothesized by Le Gac et al. (2008). This shows that the splitting of metabolic pathways across specialists can be caused by the composite effects of competition for electron donors and electron acceptors.

Consistent with differential substrate usage, average cell-specific reaction rates (Figure 6B) reveal differences in pathway activation: The transformation of acetate into acetyl-CoA by acetyl-CoA synthetase (acs) is decreased in type SS and increased in type FS, when compared to the ancestral type. Furthermore, the conversion of phosphoenolpyruvate to oxaloacetate (ppc), the conversion of phosphoenolpyruvate to pyruvate (pyk) and the decarboxylation of pyruvate to acetyl-CoA (pdh), linking the glycolysis pathway to the citric acid cycle, are all upregulated in the SS type when compared to the FS type. Similar differences in pathway activation also exist during early exponential growth in monoculture (Figure 6C,D), because FS grows partly on acetate and SS excretes acetate (Figure 4F,J). Previous microarray profiles of mRNA concentrations during exponential growth in monocultures (Le Gac et al., 2008) found an upregulation of acetate consumption genes in FS and acetate excretion genes in SS compared to A, qualitatively confirming our predictions (Figure 6C,D). Interestingly, our simulations suggest a significant downregulation of glucose catabolism (pyk, pdh and ppc) in FS compared to A, which contradicts the transcript profiles (Figure 6D). However, mRNA was harvested from well-aerated flasks, while the monoculture experiments (Figure 4) and evolution experiments (Figure 3) were performed in test tubes where oxygen can become limiting (Andersen and von Meyenburg, 1980). Oxygen becomes particularly scarce in the FS tubes (Figure 4K) and temporarily limits glucose catabolism, which would explain the strong downregulation not reflected in the transcript profiles (Le Gac et al., 2008). Furthermore, while broad pathway activation patterns could be qualitatively reproduced in our system, this might be harder in other cases due to post-transcriptional regulation or post-translational modifications (Blazier and Papin, 2012).

The periodic (seasonal) changes in glucose and acetate concentrations in batch culture have previously been shown to promote coexistence of the SS and FS types, in analogy to the maintenance of phytoplankton diversity via fluctuations of resource availability (Sommer, 1984; Spencer et al., 2007). Experiments with SS-FS batch cocultures revealed that the SS type quickly dominates over the FS type, when restricted to the first glucose-rich season through frequent dilution into fresh growth medium. Reciprocally, when SS and FS are grown in solution resembling the second glucose-depleted acetate-rich season, the FS type quickly dominates over the SS type (Spencer et al., 2007). Accordingly, in a full batch cycle the relative SS cell density has been shown to culminate within 4–8 hr and to gradually decrease afterwards (Friesen et al., 2004, Figure 6B), in consistence with our simulations (Figure 5D). Simulations of the SS and FS batch coculture restricted to the first or second season, analogous to Spencer et al.'s experiments, reproduce these observations and verify the role of periodic variation of glucose and acetate concentrations in maintaining the coexistence of both types (Figure 7, see the ‘Materials and methods’ for details).

Figure 7 with 1 supplement see all
Predicted relative cell densities of the SS and FS types in batch coculture when restricted to either the first glucose-rich (A) or second glucose-depleted (B) season.

In (A), restriction to the first season was achieved by shorter dilution periods which prevented the complete depletion of glucose. In (B), restriction to the second season was achieved by using the glucose-depleted acetate-rich solution, produced by the full-batch coculture, as growth medium (see the Methods for details). See Figure 7—figure supplement 1 for results from analogous experiments by Spencer et al. (2007).

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

Conclusions

The models presented here make detailed predictions about the microbial dynamics in the considered experiments. First, after calibration the cell models largely explain the data from the monoculture experiments (Figure 4). Second, the predictions for pathway activation in the three strains (Figure 6) are roughly consistent with transcription profiles. Third, simulations of the microbial community consisting of all three strains (Figure 5) reproduce the successional dynamics of diversification observed in the evolution experiments (Figure 3). Fourth, simulations of the SS-FS cocultures restricted to either the glucose-rich or glucose-depleted season reproduce the dominance of the SS or FS type (Figure 7), respectively, in consistence with previous co-culture experiments. It is important to note that only data from monoculture experiments were used to calibrate the cell models for the three strains (A, SS and FS). In particular, no information from co-culture experiments was used in the setup of the microbial community model, and thus there was no a priori knowledge about what the emergent community dynamics would be. Hence, our work conceptually produced non-trivial predictions that could be compared to experimental observations, although all experiments had already been performed.

Our work sheds light on the fundamental problem of metabolic diversification and the emergence of shared catabolic pathways. In particular, our model allowed quantitative predictions for the metabolic fluxes for each strain in coculture, revealing temporary cross-feeding as an emergent property of the evolved community (Treves et al., 1998). Cross-feeding, conventionally seen as a beneficial interaction (Morris et al., 2013), thus emerged as a form of niche segregation driven by competition for organic carbon and oxygen. Because both evolved types prefer glucose whenever available at high concentrations, but exchange acetate under oxygen limitation, the community constantly switches between competitive and beneficial interactions. Natural microbial populations might thus also oscillate between negative and positive interactions, for example depending on oxygen levels. Our findings also support previous suggestions that microbial evolution can be driven by deterministic ecological processes (Wood et al., 2005; Oxman et al., 2008; Herron and Doebeli, 2013). In this case, the observed diversification is due to competition for limiting resources whose use is constrained by basic metabolic tradeoffs. Other instances of ecological diversification in microbial evolution experiments, for example as reported by Plucain et al. (2014), might be explained using a similar approach.

More generally, we have demonstrated how MCM can be used to experimentally calibrate and combine genome-based cell models to predict the emergent dynamics of microbial communities. Our framework thus provides a starting point for designing microbial communities with particular metabolic properties, such as optimized catabolic performance. While MCM is designed for genome-based metabolic models, it can also accommodate conventional functional group models. In these models, different ecological functions such as photosynthesis, heterotrophy or nitrification are performed by distinct populations whose metabolic activity is determined, for example, by Michaelis–Menten kinetics and whose growth is described by simple substrate-biomass yield factors (Hood et al., 2006; Reed et al., 2014). Hence, natural microbial communities could be modeled even if annotated genomes are not available for each member species. While functional group models general require fewer parameters, their calibration remains a challenge (Panikov and Sizova, 1996). In MCM, model calibration becomes analogous to coefficient estimation in conventional multivariate regression, and can be used to estimate poorly known parameters such as stoichiometric coefficients, growth kinetics or extracellular transport coefficients (MCM user manual, Supplementary file 1, section 12). To our knowledge, no existing comparable framework offers the flexibility combined with the statistical functionality of MCM. In view of the increasing availability of genome-scale metabolic models (Feist et al., 2008), our work provides a missing link to a predictive and synthetic microbial ecology.

Materials and methods

MCM overview

MCM is a mathematical and computational framework for the construction, simulation, statistical analysis and calibration of microbial community models (Figure 2). Models are specified in special files that define all metabolites, reactions, cell species and environmental variables. MCM is controlled through custom scripts, that is, text files containing a sequence of special commands, such as for running simulations or fitting parameters. MCM includes tools for the conversion of conventional genome-scale FBA models, such as generated by the Model SEED pipeline (Henry et al., 2010) based on sequenced genomes, into a draft MC model.

MCM can accommodate microbial communities comprising genome-based cell models with arbitrary environmental variables, metabolite exchange kinetics and regulatory mechanisms. For example, environmental variables may be stochastic processes (e.g., representing climate), or specified using measured data (e.g., redox potential in bioreactor experiments), or depend on metabolite concentrations (e.g., pH determined by acetate concentration) or even be dynamical (e.g., temperature increasing at a rate proportional to biomass production rates). This versatility allows for the incorporation of complex environmental feedbacks, such as host immune responses in gut microbiota (Karlsson et al., 2011). Metabolite uptake and export rate limits can be arbitrary functions of metabolite concentrations or environmental variables. Similar interdependencies are possible for reaction rate limits, thus allowing the inclusion of inhibitory or regulatory mechanisms (Covert et al., 2008). Metabolite concentrations can be explicitly specified, for example, using measured time series, or depend dynamically on microbial export and other external fluxes. Effects of phage predation (Jensen et al., 2006), reaction energetics (Reed et al., 2014) or stochastic environments can also be incorporated.

MCM keeps track of a multitude of output variables such as cell densities, reaction rates, metabolite concentrations and metabolite exchange rates. Because each reaction can be formally associated with a particular enzyme, in turn encoded by a particular gene, MCM also makes predictions about gene densities as a product of cell densities and gene copy numbers per cell. Metabolic activity statistics (e.g., Figure 6A,B) facilitate the identification of metabolic interactions such as cross-feeding (Morris et al., 2013). The predicted time courses of output variables can be statistically evaluated against time series ranging from chemical concentrations, rate measurements to cell densities and metagenomics.

MC models can include arbitrary abstract (symbolic) numeric parameters with a predefined value range or probability distribution. Symbolic parameters can represent, for example, stoichiometric coefficients, gene copy numbers, cell life expectancies, half-saturation constants or environmental variables. The inclusion of symbolic parameters enables a high-level analysis of microbial communities: For example, MCM can automatically calibrate (fit) unknown symbolic parameters to time series using maximum–likelihood parameter estimation (Eliason, 1993). The likelihood of the data, given a particular parameter choice, is calculated by assuming a mixed deterministic-stochastic model in which the deterministic part is given by the model predictions, and the stochastic part is given by normally distributed errors. The likelihood is minimized using an iterative optimization algorithm involving step-wise parameter adjustments and repeated simulations. Other fitting algorithms are also available, such as maximization of the average coefficient of determination (R2), which is equivalent to weighted least-squares fitting. Because MCM can calibrate unknown measurement units, raw uncalibrated data (e.g., optical cell densities with no calibration to colony forming units, Figure 4A) can also be used.

In this paper single-cell models were calibrated to monoculture experiments, however models can also be calibrated using data from experimental or natural communities that include unculturable species (MCM user manual, Supplementary file 1, sections 7 and 12; Louca and Doebeli, 2015). In general, fitted parameters need not be directly connected to the data used for calibration, as long as a change in the parameters influences the predictions that are being compared to the data. While this is a general principle of parameter estimation (Tarantola, 2005), in practice the uncertainty of calibrated parameters (e.g., in terms of confidence intervals) increases when their influence on the ‘goodness of fit’ is weaker. Moreover, alternative parameter combinations can sometimes yield a comparable match to the data, especially if multiple parameters influence the same variables (inverse problem degeneracy). Local fitting optima can be detected through repeated randomly seeded calibrations (see next section), and overfitting can be partially avoided by keeping the number of free parameters at a bare minimum. Nevertheless, in certain cases good knowledge of the system or previous literature may be required to identify the most plausible calibrations. Finally, we emphasize that MCM is, after all, merely a framework enabling the construction, calibration and analysis of microbial community models. MCM models are thus limited by the same caveats and assumptions as other constraint-based metabolic models (Blazier and Papin, 2012; Antoniewicz, 2013) and any predictions made by MCM should be subject to similar scrutiny.

Calibration of E. coli cell models

E. coli strains were obtained from an evolution experiment performed in a batch culture environment with daily dilutions into glucose-acetate supplemented Davis minimal medium (Spencer et al., 2008; Tyerman et al., 2008). For each phenotype, three clones were isolated from population 20 after 150 days and used for three independent monoculture growth experiments. Optical densities, as well as glucose, acetate and oxygen concentration data from these experiments were used to calibrate the individual cell-metabolic models for the A, SS and FS phenotypes. Oxygen measurements were not available for type A. Experimental details and results are described by Le Gac et al. (2008).

In the models, the limiting nutrients are assumed to be oxygen, glucose and acetate; all other nutrients can be taken up at an arbitrary rate. Oxygen, glucose and acetate uptake rate limits were described by Monod-like kinetics. The maximum cell-specific oxygen uptake rate was set to 1.008×1013 mol/(dcell), according to Varma and Palsson (1994). The oxygen half-saturation constant was set to 1.21×10−7 M according to Stolper et al. (2010). Oxygen was assumed to be initially at atmospheric saturation levels (0.217 mM at 37 C) and repleted at a rate proportional to its deviation from saturation (Gupta and Rao, 2003).

The fitted parameters for each cell type were the maximum cell-specific uptake rates and half-saturation constants for glucose and acetate, as well as initial cell densities and non-growth associated ATP maintenance energy requirements. The initial glucose and acetate concentrations were set to the average value measured at the earliest sampling point (1 hr after incubation) for each type. The oxygen mass transfer coefficient (M/day per M deviation) was initially fitted individually for each type together with all other parameters, and then fixed to the average of all three initial fits. All other parameters were then again fitted individually for each type. Parameter fitting was done by maximizing the average coefficient of determination (R2) using the MCM command fitMCM. A total of 237 data points were used to fit 19 parameters (Table in Supplementary file 2). To reduce the possibility of only reaching a local maximum, fitting was repeated 100 times for each strain starting at random initial parameter values and the best fit among all 100 runs was used. While some fitting runs reached alternative local maxima, the best overall fit was reached in most cases.

Cell densities were directly compared to optical density (OD) measurements. The appropriate calibrations were estimated by MCM and ranged within 8.2×10111.3×1012 cells/(LOD). These estimates are consistent with previous experimental calibrations (Lawrence and Maier, 1977) yielding 0.26 g dry weight/(LOD), which corresponds to 1.4×1012 cells/(LOD) (assuming a cell dry weight of 1.8×10−13 g in the stationary phase; Fagerbakke et al., 1996).

Simulation of the microbial community model

The microbial community model was simulated using the MCM command runMCM. Initial glucose and acetate concentrations were set to the average of all values measured at the earliest sampling point of the monoculture incubations. Cell death was not explicitly included, because of lack of appropriate data for calibration and because daily dilutions by far exceeded cell death as a factor of cell population reduction.

Robustness of the SS-FS coexistence

To verify the robustness of the stable SS-FS coexistence in coculture, we randomly varied each fitted model parameter uniformly within an interval spanning 10% above and 10% below its calibrated value. Both types coexisted in 50 out of 50 random simulations (Figure 5—figure supplement 2).

Seasonal restriction of the SS-FS cocultures

Simulations of the SS-FS cocultures restricted to the first glucose-rich or second glucose-depleted season, as opposed to the full batch cycle, were performed in analogy to the experiments by Spencer et al. (2007). More precisely, to model the first season experiment we changed the dilution rate to 1/32 every 5 hr, so that at the end of each batch cycle glucose was not yet completely depleted. Similarly, for the second season experiment we changed the dilution rate to 1/32 every 19 hr, and adjusted the growth medium to resemble the glucose-depleted acetate-rich solution reported by Spencer et al. (no glucose, 3.59 mM acetate). Initial cell densities were set to 1010 cells/l for both types. All other model parameters were kept unchanged. The original experiments by Spencer et al. (2007) were performed at higher dilution rates (4 and 15 hr for the first and second season experiment, respectively), however in our simulations neither the FS nor SS type could persist at these high dilution rates. We note that the strains used in our work (Le Gac et al., 2008) had evolved in separate evolution experiments using a different growth medium than those by Spencer et al. (2007).

Obtaining MCM

MCM is open source and available at http://www.zoology.ubc.ca/MCM.

References

  1. 1
    Are growth rates of Escherichia coli in batch cultures limited by respiration?
    1. KB Andersen
    2. K von Meyenburg
    (1980)
    Journal of Bacteriology 144:114–123.
  2. 2
  3. 3
    Integration of expression data in genome-scale metabolic network reconstructions
    1. AS Blazier
    2. JA Papin
    (2012)
    Frontiers in Physiology, 3, 10.3389/fphys.2012.00299.
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
    Integrating metabolic, transcriptional regulatory and signal transduction models in Escherichia coli
    1. MW Covert
    2. N Xiao
    3. TJ Chen
    4. JR Karr
    (2008)
    Bioinformatics 24:2044–2050.
  12. 12
  13. 13
  14. 14
    Maximum likelihood estimation: logic and practice
    1. SR Eliason
    (1993)
    Newbury Park, CA: SAGE Publications.
  15. 15
    Chemical oceanography and the marine carbon cycle
    1. S Emerson
    2. J Hedges
    (2008)
    Cambridge, UK: Cambridge University Press.
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
    Experimental evidence for sympatric ecological diversification due to frequency-dependent competition in Escherichia coli
    1. ML Friesen
    2. G Saxer
    3. M Travisano
    4. M Doebeli
    (2004)
    Evolution 58:245–260.
  21. 21
  22. 22
  23. 23
    Monte carlo methods
    1. JM Hammersley
    2. DC Handscomb
    (1964)
    Methuen.
  24. 24
  25. 25
  26. 26
  27. 27
    Integrated analysis of regulatory and metabolic networks reveals novel regulatory mechanisms in Saccharomyces cerevisiae
    1. MJ Herrgård
    2. BS Lee
    3. V Portnoy
    4. BØ Palsson
    (2006)
    Genome Research 16:627–635.
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
    Correction for the inherent error in optical density readings
    1. J Lawrence
    2. S Maier
    (1977)
    Applied and Environmental Microbiology 33:482–484.
  37. 37
    Metabolic changes associated with adaptive diversification in Escherichia coli
    1. M Le Gac
    2. MD Brazas
    3. M Bertrand
    4. JG Tyerman
    5. CC Spencer
    6. REW Hancock
    7. M Doebeli
    (2008)
    Genetics 178:1049–1060.
  38. 38
    Transient dynamics of competitive exclusion in microbial communities
    1. S Louca
    2. M Doebeli
    (2015)
    Environmental Microbiology, 10.1111/1462-2920.13058.
  39. 39
    Dynamic flux balance analysis of diauxic growth in Escherichia coli
    1. R Mahadevan
    2. JS Edwards
    3. FJ Doyle III
    (2002)
    Biophysical Journal 83:1331–1340.
  40. 40
  41. 41
  42. 42
    Chemical oceanography. Marine science series (4th edn)
    1. FJ Millero
    (2013)
    Bosa Roca: Taylor & Francis.
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
    Syntrophism among prokaryotes
    1. B Schink
    2. A Stams
    (2006)
    In: M Dworkin, S Falkow, E Rosenberg, KH Schleifer, E Stackebrandt, editors. The prokaryotes. New York: Springer. pp. 309–335.
  51. 51
    Analysis of optimality in natural and perturbed metabolic networks
    1. D Segrè
    2. D Vitkup
    3. GM Church
    (2002)
    Proceedings of the National Academy of Sciences of USA 99:15112–15117.
  52. 52
  53. 53
    Seasonal resource oscillations maintain diversity in bacterial microcosms
    1. CC Spencer
    2. G Saxer
    3. M Travisano
    4. M Doebeli
    (2007)
    Evolutionary Ecology Research 9:775.
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
    Inverse problem theory and methods for model parameter estimation
    1. A Tarantola
    (2005)
    Philadelphia: Society for Industrial and Applied Mathematics.
  59. 59
  60. 60
  61. 61
    Stoichiometric flux balance models quantitatively predict growth and metabolic by-product secretion in wild-type Escherichia coli W3110
    1. A Varma
    2. BO Palsson
    (1994)
    Applied and Environmental Microbiology 60:3724–3731.
  62. 62
  63. 63
  64. 64
  65. 65

Decision letter

  1. Wenying Shou
    Reviewing Editor; Fred Hutchinson Cancer Research Center,, United States

eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see review process). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.

Thank you for submitting your work entitled “Predicting the dynamics of microbial communities using genome-based metabolic models” for peer review at eLife. Your submission has been favorably evaluated by Naama Barkai (Senior Editor) and three reviewers, one of whom is a member of our Board of Reviewing Editors.

The reviewers have discussed the reviews with one another and the Reviewing editor has drafted this decision to help you prepare a revised submission.

This manuscript presents (a) a computational pipeline for calibrating genome-scale models of metabolism through fitting to experimental data and (b) the application of calibrated E. coli models to the study of a previously published evolutionary experiment, in which two mutants were found to coexist, after outcompeting the ancestral strain. What have we learned? From the point of view of this biological system, this paper provided yet another piece of data to the many already existing that the niche for FS emerges after the SS rise in frequency, but does not go to fixation. This has been known for quite a while, but it is nice to see the model agree with this finding. The calibration of individual genome scale models prior to studying communities seems a nice addition to existing frameworks, in fact one that should become a standard approach, whenever possible, for ecosystem-level modeling. The modeling and calibration software developed by the authors seems very well documented.

Shared major concerns:

1) Your paper did not explicitly compare model predictions with experiments. We would like to see a comparison of predicted metabolic and population dynamics with experimental coculture dynamics. There are no new experiments to test any predicted effect of perturbations. To at least test that FS can't invade alone is totally easy. Presumably one could also find an acetate concentration that would have been high enough to obviate the need for SS to come first (or a shorter duration or larger dilution that would have left enough acetate by the time of the transfer to have had the same effect). Another major prediction – that was pretty cool – is that acetate is excreted during a short window of time. This would be an excellent experiment to validate the model's predictions. Finally, you have to actually show us the comparison of predicted fluxes and mRNA that is supposed to be such a great fit. If it really is that good, we'd be rather impressed, for there is a huge literature about how this is not the case.

2) MCM seems to be sold too hard and does not responsibly acknowledge other DFBA platforms for communities like COMETS that do nearly everything mentioned here except the statistical fitting of parameters. We rather like that extension, but it seems to be advertised as more than it is.

3) Given that the authors propose their approach as broadly applicable to studying microbial communities, we think it may be important for them to comment on the realistic applicability of the method beyond E. coli wild type and mutants, or other well-annotated models. Would their method be useful, for example, for improving models that lack precision not just at the level of uptake kinetics, but at the level of the stoichiometry itself?

Specific major concerns:

Reviewer 1:

1) I don't understand Figure 3. If SS released acetate under oxygen limitation (∼1/4 day after daily dilution, Figure 3 i,j), then why should acetate be accumulated that quickly – immediately after daily dilution (Figure 3f)?

2) The Introduction should clearly discuss where the field stands currently and how MCM pushes the field forward. I am not sure how, for example, COMETS of the Segre group could have predicted the spatial-temporal dynamics in microbial communities if they had not somehow calibrated their model. Only after I talked to an expert did I realize that usually parameter choice is done through looking for values in literature, which leaves a few tunable parameters. Then, parameters are adjusted manually to fit experimental dynamics. This is because some of these parameters may not be directly measured experimentally. A thoughtful discussion on this will be helpful.

3) The end of conclusion mentioned that parameter estimation does not necessarily require monoculture measurements. This is a critical point, and should be formally demonstrated (rather than hidden in a supplementary file). For example, the authors could model the three-member community with parameters derived from cocultures of two members starting at arbitrary initial compositions. This is to mimic cases (e.g. soil microbes) where many species are not individually culturable.

4) The flow of the paper is suboptimal, especially to an outside reader. For example, you can move “MCM overview” to immediately after the Introduction. You may also want to add concrete examples to your figures to demonstrate how MCM works in reality.

Reviewer 2:

1) The novelty of this work is mainly in the combination of different approaches and data, rather than in the approaches and data themselves. I find fascinating that the model recapitulates the observations, and that model-predicted fluxes are consistent with previously measured gene expression data. However, it is not clear to me what aspects of the insight provided by the model were not known or suspected before, given that extensive work was done on this system.

2) As mentioned above, I like the calibration approach. However, I think that some important information is missing. First, there should be a table (fine as supplementary) detailing the values of fitted parameters, and any available comparative values from the literature (only a few examples are provided in the Methods section). Also: is the fitted ATP maintenance the value for the non-growth associated maintenance, or the coefficients of growth-associated maintenance? Second, it would important for readers to know whether the solutions found by the fitting algorithm are unique, and how sensitive the result are to parameter precision. Sensitivity analyses are discussed in the user manual, and seem to have been applied to a different microbial community. I think it would be particularly important to know whether the main result of stable coexistence is sensitive to the choice of fitted parameters.

Reviewer 3:

1) The text reads that “these observations are in exact agreement with microarray transcriptional profiles”. Given the strength of “exact agreement”, I was very surprised that there was no display of experiment versus model. This was a major finding, but no figures or statistical analysis of what “exact” means.

A related point: it is claimed that “MCM also makes predictions about gene densities” because each flux is associated with an enzyme. This is actually a fairly ludicrous claim. There have been paper after paper showing that, not only can there be a lack of correlation between flux and enzyme activity, they can even be negatively correlated. This is a central tenet of Metabolic Control Analysis, and has been well documented and commented upon by folks such as Dan Fraenkel (2003, Current Opinion in Microbiology) and Hans Westerhoff (Rossell et al., 2006. PNAS). This possibility must be addressed to make the reader aware that there should never be the assumption that flux is proportionally related to enzyme levels for all enzymes. If that were so, every enzyme would have a control coefficient of 1, which is absolutely impossible because the sum of control coefficients from the entire cell is 1.

Returning to the point above, then, if there really is a good quantitative correlation between the array data from the Le Gac work and the DFBA model here – for matching timepoints as the mRNA was harvested – it would be quite a nice finding. This analysis is absolutely essential to the paper.

2) An advance stated is MCM itself. I am not sold upon exactly how this is really a major advance beyond a variety of other DFBA approaches, including prior work on communities. For example, a very large number of MCM features that are described in a way that comes across as them being novel also exist in COMETS (2014, Cell Reports). At the very least, acknowledge where this is so, and use it as an opportunity to more precisely say where this goes above and beyond. For example, I like the inherent fitting of data and perhaps more could be made of exactly how this works. COMETS does not use fitting to each substrate, and frankly an objective means to parameterize when canonical parameters fail to work well is a nice step.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled “Calibration and analysis of genome-based models for microbial ecology” for further consideration at eLife. Your revised article has been favorably evaluated by Naama Barkai (Senior Editor) and three reviewers, one of whom is a member of our Board of Reviewing Editors.

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

Reviewer 1:

For this to be a useful Tools and Resources article, I'd still like to see a discussion on MCM's limitations and requirements in the main text. For example, what is enough data for calibrating a model? Under what situations will MCM give you several solutions and what do you do in this case? What should be avoided or advocated when using MCM? You probably discussed these in the manual, but it will be helpful to summarize the main points in the text.

Reviewer 2:

The authors have overall addressed my concerns, partially by answering the issues raised, partially by transforming their work into a Tools and Resources article.

There are still two points that need some revision:

1) Even if biological insight is no more the main focus of the paper, I still think that the authors could explain a little bit in the Introduction why people may care about the specific example they study. For example, I like the following rebuttal of the authors to one of the questions: “While most of the results have been found experimentally over the course of several years, it is only now that a mechanistic model has managed to unify many of them in such a clear, unambiguous and synergistic manner. These modeling results provide very strong credence to a large body of experimental work that was done in our lab over the course of roughly a decade.” I think that a slightly expanded version of this text would help orient the broad readership towards understanding why their example is interesting (in addition to the fact that it works).

2) Most importantly, I think that in talking about their software, as well as previous tools, the authors should be very careful to clearly distinguish between features that ere possible (currently or in principle) in the different tools, vs. features that are actually presented in detail and tested in the manuscript. I am referring in particular to the last paragraph in the subsection “Model”, which is quite problematic in a number of ways. For example, regulation (both allosteric and transcriptional) in FBA models is notoriously a very tricky, overall unresolved problem. MCM has the potentially useful feature of allowing users to set rules to limit fluxes as a function of other parameters. This feature, described in one page in the user manual, is oversold in the main text as a capacity of MCM to include regulation. Similarly superficial is the description of the inclusion of phages in the model. Again, this is described shortly in the user manual with no justification or testing. All these features are unnecessarily used as a way of contrasting MCM with previous software, leading to unjustified conclusions, e.g. that COMETS “offers limited model versatility in terms of uptake kinetics” (COMETS does allow different parameters for different molecules and organisms, even if parameters were chosen to be equal in the specific simulations presented in the COMETS paper), that COMETS “seems limited to Petri dishes”, whereas “MCM can be used to understand the dynamics of realistic microbial communities, ranging from the soil or groundwater to artificial communities and bioreactors”. Both COMETS and MCM, after all, were based on the same underlying modeling framework, and tested on laboratory systems. I think the authors should definitely mention (with added accuracy) – perhaps in the Discussion – the additional features of their software that are not described and tested in detail in the main text. At the same time, I think they should limit their claims of major innovation to the components that are actually tested and presented in detail. The capacity to perform calibration on individual organisms, followed by proof of predictive capacity of the global dynamics on a highly interesting system is elegant and brilliant, and I think it is unnecessarily weakened by these other dubious claims.

Reviewer 3:

First off, I strongly agree with Reviewer 2's comments above. They have done much to improve the paper to show more and fix much of the language. As noted, however, saying why this matters in the paper as well as they did in the rebuttal would be great, and they still need to back down from hyping a large series of bells and whistles that may be great features but are not demonstrated here. There are plenty of nice advances in their method and they should stick to these.

As for my own concerns, the direct comparison to data and clarification that expression and flux should not be expected to perfectly correlate is a massive relief. The qualitative agreement is not bad in Figure 6c and d, and can rest on its own laurels. The figures are also particularly improved in their clarity.

I still find the concept of an “average flux” over such a dynamic experiment to be strange, but at least it is now explained. Using the same time point as the RNA data is far better.

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

Author response

Shared major concerns:

1) Your paper did not explicitly compare model predictions with experiments. We would like to see a comparison of predicted metabolic and population dynamics with experimental coculture dynamics. There are no new experiments to test any predicted effect of perturbations. To at least test that FS can't invade alone is totally easy. Presumably one could also find an acetate concentration that would have been high enough to obviate the need for SS to come first (or a shorter duration or larger dilution that would have left enough acetate by the time of the transfer to have had the same effect).

To address this concern, we would like to first clarify our methodology: to calibrate our microbial community model we have only used data from monoculture experiments, and we have completely ignored any information from co-culture experiments. This calibration yields the parameters that we used to characterize each of the three strains (ancestor, SS and FS). We then put these three calibrated strains into an artificial MCM community, which predicted emergent dynamics of the microbial community. Because the model strains were calibrated in monoculture, there was no guarantee as to what these emergent dynamics would be. Our main result is that in fact, the emergent dynamics very closely correspond to those that had previously been found in a series of co-culture experiments. Thus, conceptually our model yielded predictions that were borne out by experiments, albeit experiments that had already been performed. Nevertheless, the editors' and reviewers' suggestion for a more extensive experimental validation is well taken, and we have now expanded that part of the manuscript.

In summary:

1) Our MCM simulations predict a clear succession of the three phenotypes and an eventual coexistence of the SS and FS phenotypes (Figure 5a). This has been shown repeatedly in previous evolution experiments.

2). Our simulations of a single batch of FS and SS co-culture predict a culmination of the SS relative density within 4-8 hours, and a subsequent decline due to the growth of the FS type (Figure 5d). This has been shown experimentally in Friesen et al. (2004, Figure 6B therein). We now discuss this observation in the manuscript (last paragraph of the “Predicting microbial community dynamics” section).

3) Our simulations of the FS and SS co-culture predict a small increase in acetate concentration, shortly after depletion of the glucose within each batch cycle (Figure 5b and 5j). This has been shown experimentally by Spencer et al. (2007, p. 779). We have added an appropriate comment to the manuscript (second paragraph of the “Predicting microbial community dynamics” section).

4) Our simulations predict differential pathway activation between the SS and FS types in co-culture, particularly of the acs, ack, pta, ppc, pdh and pyk genes, which are important determinants of the acetate and glucose pathways. These patterns are qualitatively consistent with microarray transcript profiles by Le Gac et al. (2008). Also see our response to related questions by the reviewers.

5) Additional MCM simulations of the FS and SS co-culture, restricted to either the glucose-rich season or the glucose-depleted season, predict a rapid extinction of the FS or SS type, respectively (Figure 7). This behavior was found in the experiments by Spencer et al. (2007) and we now discuss this in the article (last paragraph of the “Predicting microbial community dynamics” section).

Again, it is important to point out that all the experimental findings mentioned in points 1.-5. above were obtained in experimental co-cultures, whereas our calibration involved only monocultures, so that the subsequent MCM simulations with the calibrated strains in co-culture yield dynamics that are not part of the model assumptions, but instead is an emergent property of the MCM system. The fact that this dynamics closely corresponds to various salient experimental findings in co-culture constitutes a validation of the MCM predictions.

Regarding the invasion of the FS type in the absence of SS: Simulations predict that FS also invades in the absence of type SS, however the invasion speed is much slower and FS reaches lower densities than in the presence of the SS type. This is consistent with an early presence of the FS lineage at low densities in the evolution experiments, as reported in Herron and Doebeli (2013), indicating that some of the first FS mutations already confer a slight advantage over the ancestor when FS is rare. We now discuss these observations in the article (first paragraph of the “Predicting microbial community dynamics” section), and we have added the appropriate figures (Figure 5-figure supplement 1 and Figure 3 in the main article).

Performing additional co-culture experiments is, unfortunately, easier said than done because these experiments were performed several years ago and the appropriately skilled staff has long left our lab. Our goal is to present a novel mathematical and computational framework, and to use our wealth of pre-existing experimental data to validate its practicality and accuracy. The fact that the experiments were done beforehand, we argue, does not influence information flow nor does it compromise the validity of our arguments.

Another major prediction – that was pretty cool – is that acetate is excreted during a short window of time. This would be an excellent experiment to validate the model's predictions.

Previous SS-FS coculture experiments (Spencer et al., 2007) have indeed shown a temporary increase of acetate concentration towards the end of the SS exponential growth phase, validating our prediction. We now mention this in the article (second paragraph of the “Predicting microbial community dynamics” section).

Finally, you have to actually show us the comparison of predicted fluxes and mRNA that is supposed to be such a great fit. If it really is that good, we'd be rather impressed, for there is a huge literature about how this is not the case.

Our predicted pathway activation patterns are qualitatively consistent with the mRNA measurements by Le Gac et al. (2008). That is, the upregulation of acetate consumption genes in FS and acetate production genes in SS, when compared to the ancestor, is reflected in the transcript profiles. We now elaborate more on this comparison in the manuscript (second-last paragraph of the “Predicting microbial community dynamics” section). We also added a figure comparing the predictions to the experimental data (Figure 6c,d). Please also see our responses to related questions by the other reviewers.

2) MCM seems to be sold too hard and does not responsibly acknowledge other DFBA platforms for communities like COMETS that do nearly everything mentioned here except the statistical fitting of parameters. We rather like that extension, but it seems to be advertised as more than it is.

We now modified our manuscript to include a better description of existing techniques (including COMETS; Harcombe et al. 2014). We also better highlight what we think are the novelties of our work. We encourage the reviewers to look at MCM's User Manual (Supplementary file 2) to see the versatility and wealth of functionalities introduced by MCM.

For example, COMETS offers limited model versatility in terms of uptake kinetics, no gene regulation and few environmental feedback mechanisms (namely, only dynamic extracellular metabolite concentrations). Furthermore, it assumes complete knowledge of all required model parameters and provides no generic statistical model analysis. Hence, while COMETS sets an important precedent, it is of limited applicability to microbial ecological questions outside of petri dishes.

MCM extends Harcombe et al.'s approach to more versatile microbial ecological models that include, for example, arbitrary dynamical environmental variables (e.g. pH and temperature), stochastic processes, phage predation, product inhibition, arbitrary uptake kinetics, reaction energetics and gene regulation. Furthermore, MCM can perform statistical model evaluation, sensitivity analysis, Monte Carlo-based uncertainty analysis as well as automatic parameter fitting to data. We now point out these extensions in the Model section.

3) Given that the authors propose their approach as broadly applicable to studying microbial communities, we think it may be important for them to comment on the realistic applicability of the method beyond E. coli wild type and mutants, or other well-annotated models.

While MCM is designed for genome-based metabolic models, it can also accommodate conventional functional models in which the growth of each functional group is determined, for example, by simple Michaelis-Menten kinetics (Hood et al., 2006) or the thermodynamics of biocatalyzed redox reactions (Reed et al., 2014). Hence, natural microbial communities could be modeled even if annotated genomes are not available for each member species. While functional models require fewer parameters, their calibration remains a challenge (Panikov and Sizova, 1996). MCM would thus also be useful to such conventional modeling efforts. We now discuss this in the Conclusions.

Would their method be useful, for example, for improving models that lack precision not just at the level of uptake kinetics, but at the level of the stoichiometry itself?

Thanks to MCM's design, any numerical model parameter can in principle be calibrated, provided that suitable data is available. That includes stoichiometric coefficients, uptake kinetics, extracellular transport coefficients (e.g. for O2, as in our model), cell death rates, product inhibition kinetics or environmental pH levels. We now make this more clear in the Conclusions and the Methods.

Specific major concerns:

Reviewer 1:

1) I don't understand Figure 3. If SS released acetate under oxygen limitation (∼1/4 day after daily dilution, Figure 3 i,j), then why should acetate be accumulated that quickly – immediately after daily dilution (Figure 3f)?

Acetate is included in the growth medium with every dilution (see the Methods for details). In fact, the acetate produced by the SS type is only a small fraction of the total acetate present. We now make this more clear in section “Experimental calibration”, by emphasizing that the dilution is performed “into fresh glucose-acetate supplemented medium”. Please note that Figure 3 is now renamed to Figure 5.

2) The Introduction should clearly discuss where the field stands currently and how MCM pushes the field forward.

We now added the following to the Introduction:

“So far, the standard approach has been to obtain each parameter through laborious specific measurements or from the available literature, or to manually adjust parameters to match observations (Mahadevan et al., 2002; Chiu et al., 2014; Harcombe et al., 2014). Furthermore, statistical model evaluation and sensitivity analysis is typically performed using ad-hoc code, thus increasing the effort required for the construction of any new model.”

We now also moved the Model section so that it immediately follows the Introduction. There we explain the underlying concepts and clarify our contributions to the field of microbial community modeling (last paragraph).

I am not sure how, for example, COMETS of the Segre group could have predicted the spatial-temporal dynamics in microbial communities if they had not somehow calibrated their model.

Segre's group chose uptake kinetic parameters, as well as diffusion coefficients, to be equal for all substrates and based on typical values from the literature. Not doubting the quality of their work, we are unsure about the biological realism of such assumptions, but we prefer not to discuss this in the paper. Certainly, in many other cases a differentiation of kinetics between different metabolites and a calibration of parameters will be needed. Currently, no systematic automated method exists for estimating arbitrary model parameters.

Only after I talked to an expert did I realize that usually parameter choice is done through looking for values in literature, which leaves a few tunable parameters. Then, parameters are adjusted manually to fit experimental dynamics. This is because some of these parameters may not be directly measured experimentally. A thoughtful discussion on this will be helpful.

We now discuss this in the Introduction and the Conclusions. In the Model section (last paragraph) we further elaborate on the advancements of MCM compared to COMETS.

3) The end of conclusion mentioned that parameter estimation does not necessarily require monoculture measurements. This is a critical point, and should be formally demonstrated (rather than hidden in a supplementary file). For example, the authors could model the three-member community with parameters derived from cocultures of two members starting at arbitrary initial compositions. This is to mimic cases (e.g. soil microbes) where many species are not individually culturable.

We designed MCM so that it is completely irrelevant whether the fitted parameters belong to a one- or multiple-species model. Including an additional example in the main article demonstrating this would, in our opinion, not add anything new conceptually, neither is it required for the completeness of our arguments. In fact, as we explain above, our goal was to demonstrate that independent data from monoculture experiments (i.e. without knowledge of the dynamics emerging in co-culture) can be used to calibrate a model, which would then predict novel emergent dynamics of a microbial community.

We do realize that some readers might be skeptical about MCM's ability to calibrate multi-species models. Hence, we now specify the section in the supplement where the interested reader can see an elaborate example using real data. We are reluctant to include that example in the main article, because it would confuse our main example of the E. coli evolution experiment. We would be happy to create a separate supplementary section focused on such an example, if the editors or reviewers believe that is appropriate.

We have recently completed a separate project where we use MCM to study the stability of microbial communities in bioreactors. That work is based on calibration of a multi-species model to real data, and has recently been submitted for review in another journal.

4) The flow of the paper is suboptimal, especially to an outside reader. For example, you can move “MCM overview” to immediately after the Introduction.

We have now restructured the manuscript in line with the reviewer's suggestions. We now introduce the model framework (in particular FBA and DFBA) and MCM's principles right after the Introduction. We then introduce the microbial evolution experiments. We still keep a more technical “MCM overview” section in the Methods.

You may also want to add concrete examples to your figures to demonstrate how MCM works in reality.

We have now added appropriate comments to the MCM overview figure and its legend (Figure 2 in the revised manuscript). We hope that the Methods section and our extensive MCM User Manual (Supplementary file 2) will provide sufficient information to the interested reader.

Reviewer 2:

1) The novelty of this work is mainly in the combination of different approaches and data, rather than in the approaches and data themselves. I find fascinating that the model recapitulates the observations, and that model-predicted fluxes are consistent with previously measured gene expression data. However, it is not clear to me what aspects of the insight provided by the model were not known or suspected before, given that extensive work was done on this system.

Our main contribution is the computational framework, MCM, which combines microbial community model construction with statistical model evaluation, sensitivity analysis and automated, statistically rigorous calibration to data. The versatility and functionalities of MCM are discussed and demonstrated throughout our paper. We now elaborate more on the novelties of MCM compared to existing techniques (e.g. in the Introduction and the Model sections).

Furthermore, to the best of our knowledge, this is the first time that monocultures of multiple strains were used to calibrate individual FBA models, which were then combined into a microbial community to predict the dynamics emerging in co-culture.

We go through some lengths to show that this approach, and MCM for that matter, is capable of correctly predicting the dynamics of microbial communities, both in terms of population dynamics as well as pathway activation patterns and metabolic fluxes (e.g. periodic acetate production by SS). In view of the conventional approach for obtaining model parameters (see Introduction), we expect that systematic model calibration as demonstrated in our paper will find widespread use in microbial systems biology.

In terms of the results on the E. coli system, we agree that our paper is mainly a synthesis. While most of the results have been found experimentally over the course of several years, it is only now that a mechanistic model has managed to unify many of them in such a clear, unambiguous and synergistic manner. These modeling results provide very strong credence to a large body of experimental work that was done in our lab over the course of roughly a decade.

To conclude, as we focus mainly on methodological novelties and synthesis, our paper is better suited as a Tools and Resources paper. As suggested, we have adjusted our title accordingly.

2) As mentioned above, I like the calibration approach. However, I think that some important information is missing. First, there should be a table (fine as supplementary) detailing the values of fitted parameters, and any available comparative values from the literature (only a few examples are provided in the Methods section). Also: is the fitted ATP maintenance the value for the non-growth associated maintenance, or the coefficients of growth-associated maintenance?

The fitted ATP maintenance is the non-growth associated maintenance, in terms of ATP molecules per unit time. The appropriate “loss reaction” is part of the E. coli FBA model template, published by Yoon et al. (2012). Growth-associated maintenance is incorporated as inefficiencies of biomass synthesis from precursors biomolecules. We now make that more clear in the Methods section. We have also included a table listing the fitted parameters (Supplementary file 3), as suggested by the reviewer.

The value of non-growth associated maintenance energy was calibrated, along with the other model parameters, using the available monoculture data. For clarification, in general fitted parameters need not be directly connected to the data used for calibration, as long as a change of the parameters influences the variables that are being compared to the data. This is a general principle of parameter estimation, which is further exemplified in the MCM user manual (e.g. section 2.6). However, the uncertainty of calibrated parameters increases when their influence on the “goodness of fit” is weaker.

Second, it would important for readers to know whether the solutions found by the fitting algorithm are unique, and how sensitive the result are to parameter precision. Sensitivity analyses are discussed in the user manual, and seem to have been applied to a different microbial community. I think it would be particularly important to know whether the main result of stable coexistence is sensitive to the choice of fitted parameters.

To reduce the possibility of only reaching a local maximum, fitting was repeated 100 times for each strain starting at random initial parameter values and the best fit among all 100 runs was used (MCM handles this automatically). While some fitting runs reached alternative local maxima, the best overall fit was reached in most cases. We have now added this information to the Methods section.

We also analyzed the sensitivity of stable SS-FS coexistence in co-culture by randomly choosing each fitted parameter within an interval spanning 10% below and 10% above its calibrated value. In 50 out of 50 random simulations both types coexisted. We now mention this in the Methods section and have added an appropriate figure to the supplement (Figure 5-figure supplement 2).

Reviewer 3:

1) The text reads that “these observations are in exact agreement with microarray transcriptional profiles”. Given the strength of “exact agreement”, I was very surprised that there was no display of experiment versus model. This was a major finding, but no figures or statistical analysis of what “exact” means.

FBA makes predictions about reaction rates, and a correlation with mRNA levels is expected to be at most qualitative, not quantitative.

Our choice of wording was unfortunate, and we now corrected it to emphasize the qualitative nature of the agreement. That is, the upregulation of acetate consumption genes in FS and acetate production genes in SS, when compared to the ancestor, is reflected in the mRNA profiles. We now elaborate this comparison in the manuscript (at the end of the “Predicting microbial community dynamics” section). We also added a figure comparing the predictions to the experimental data (Figure 6c,d).

A related point: it is claimed that “MCM also makes predictions about gene densities” because each flux is associated with an enzyme. This is actually a fairly ludicrous claim. There have been paper after paper showing that, not only can there be a lack of correlation between flux and enzyme activity, they can even be negatively correlated. This is a central tenet of Metabolic Control Analysis, and has been well documented and commented upon by folks such as Dan Fraenkel (2003, Current Opinion in Microbiology) and Hans Westerhoff (Rossell et al., 2006. PNAS). This possibility must be addressed to make the reader aware that there should never be the assumption that flux is proportionally related to enzyme levels for all enzymes. If that were so, every enzyme would have a control coefficient of 1, which is absolutely impossible because the sum of control coefficients from the entire cell is 1.

We appreciate the reviewer's (correct) observation. To clarify, we do not assume that flux is correlated with enzyme activity in any way. Predictions about gene densities (i.e. gene copies per liter) made by MCM are based on predicted cell densities and the number of gene copies per cell. Hence, for example, the density of the amo gene (associated with aerobic ammonium oxidation) would be the sum of cell densities for all cell species capable of producing the amo enzyme (i.e., in the context of FBA, capable of potentially performing the amo reaction) multiplied by the number of amo gene copies per cell. Gene copies can be different across cell species. This is explained in detail in the MCM user manual. We now also make that more clear in the manuscript (“MCM overview” section in the Methods).

Of course, the very assumption that each reaction corresponds to one gene via one enzyme can be questioned, e.g. when protein complexes are encoded by several genes. Hence, the final interpretation or use of MCM's predictions about gene densities should be case-dependent and subject to the user's judgment.

Returning to the point above, then, if there really is a good quantitative correlation between the array data from the Le Gac work and the DFBA model here – for matching timepoints as the mRNA was harvested – it would be quite a nice finding. This analysis is absolutely essential to the paper.

We now elaborate on that comparison in more detail. We also added a figure, which directly compares our predictions with Le Gac's results (Figure 6c,d). To clarify, we do not claim a quantitative correlation exists between reaction rates and mRNA concentrations. Rather, our predictions on the upregulation of acetate production genes in SS and acetate consumption genes in FS, when compared to the ancestor, are consistent with analogous changes in mRNA levels. We now make that more clear.

2) An advance stated is MCM itself. I am not sold upon exactly how this is really a major advance beyond a variety of other DFBA approaches, including prior work on communities. For example, a very large number of MCM features that are described in a way that comes across as them being novel also exist in COMETS (2014, Cell Reports). At the very least, acknowledge where this is so, and use it as an opportunity to more precisely say where this goes above and beyond. For example, I like the inherent fitting of data and perhaps more could be made of exactly how this works. COMETS does not use fitting to each substrate, and frankly an objective means to parameterize when canonical parameters fail to work well is a nice step.

As mentioned in our responses to Reviewer 1, in contrast to COMETS, MCM can accommodate more versatile microbial ecological models that include, for example, dynamical environmental variables (such as pH or temperature, that can change in response to microbial metabolism), stochastic processes, phage predation, product inhibition, arbitrary uptake kinetics, reaction energetics and gene regulation. Furthermore, MCM can perform statistical model evaluation and sensitivity analysis, as well as automatic parameter fitting to data.

In the Model section (last two paragraphs) we now compare our work to previous DFBA approaches, and in particular COMETS, and we further emphasize the importance of MCM's fitting capacities in the Conclusions section. We now also elaborate more on the fitting algorithms used by MCM in the Methods section.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Reviewer 1:

For this to be a useful Tools and Resources article, I'd still like to see a discussion on MCM's limitations and requirements in the main text. For example, what is enough data for calibrating a model? Under what situations will MCM give you several solutions and what do you do in this case? What should be avoided or advocated when using MCM? You probably discussed these in the manual, but it will be helpful to summarize the main points in the text.

The amount of required data is highly case-specific and indeed a non-trivial question that cannot be answered in general. MCM calibrations will always be subject to the same fundamental limitations as other non-linear inverse problems. We now added a discussion of the points raised by the reviewer, in the “MCM overview” section: “In this paper single-cell models were calibrated to monoculture experiments […] MCM models are thus limited by the same caveats and assumptions as other constraint-based metabolic models (Blazier and Papin, 2012; Antoniewics, 2013) and any predictions made by MCM should be subject to similar scrutiny.”

Reviewer 2:

The authors have overall addressed my concerns, partially by answering the issues raised, partially by transforming their work into a Tools and Resources.

There are still two points that need some revision:

1) Even if biological insight is no more the main focus of the paper, I still think that the authors could explain a little bit in the Introduction why people may care about the specific example they study. For example, I like the following rebuttal of the authors to one of the questions: “While most of the results have been found experimentally over the course of several years, it is only now that a mechanistic model has managed to unify many of them in such a clear, unambiguous and synergistic manner. These modeling results provide very strong credence to a large body of experimental work that was done in our lab over the course of roughly a decade.” I think that a slightly expanded version of this text would help orient the broad readership towards understanding why their example is interesting (in addition to the fact that it works).

We now finish the Introduction with the following passage: “To demonstrate the potential of MCM, we modeled a bacterial community that has emerged from in-vitro evolution experiments […] a large body of experimental work that was done in our lab over the course of roughly a decade.”

2) Most importantly, I think that in talking about their software, as well as previous tools, the authors should be very careful to clearly distinguish between features that ere possible (currently or in principle) in the different tools, vs. features that are actually presented in detail and tested in the manuscript. I am referring in particular to the last paragraph in the subsection “Model”, which is quite problematic in a number of ways. For example, regulation (both allosteric and transcriptional) in FBA models is notoriously a very tricky, overall unresolved problem. MCM has the potentially useful feature of allowing users to set rules to limit fluxes as a function of other parameters. This feature, described in one page in the user manual, is oversold in the main text as a capacity of MCM to include regulation. Similarly superficial is the description of the inclusion of phages in the model. Again, this is described shortly in the user manual with no justification or testing. All these features are unnecessarily used as a way of contrasting MCM with previous software, leading to unjustified conclusions, e.g. that COMETS “offers limited model versatility in terms of uptake kinetics” (COMETS does allow different parameters for different molecules and organisms, even if parameters were chosen to be equal in the specific simulations presented in the COMETS paper), that COMETS “seems limited to Petri dishes”, whereas “MCM can be used to understand the dynamics of realistic microbial communities, ranging from the soil or groundwater to artificial communities and bioreactors”. Both COMETS and MCM, after all, were based on the same underlying modeling framework, and tested on laboratory systems. I think the authors should definitely mention (with added accuracy) – perhaps in the Discussion – the additional features of their software that are not described and tested in detail in the main text. At the same time, I think they should limit their claims of major innovation to the components that are actually tested and presented in detail. The capacity to perform calibration on individual organisms, followed by proof of predictive capacity of the global dynamics on a highly interesting system is elegant and brilliant, and I think it is unnecessarily weakened by these other dubious claims.

Regarding gene regulation:

Incorporating post-transcriptional or post-translational gene regulation in FBA models is indeed tricky for fundamental reasons, and we likely have overstated MCM's potential to that end. The incorporation of gene regulation into MCM models follows published regulatory FBA approaches (Covert et al., 2001; Covert and Palsson, 2002; Covert et al., 2004; Herrgard et al., 2006) that have been shown to improve the fidelity of E. coli and S. cerevisiae FBA models, but are yet to be tested at the microbial community level. Hence, the underlying theory is definitely not our innovation, and MCM merely provides a convenient tool for constructing and calibrating such models. We now clarify this in the last paragraph of the “Model” section. We also further explain regulatory FBA and non-stationary cell models in the MCM User Manual.

Regarding COMETS:

We would like to clarify that while COMETS does indeed allow different parameters for different molecules and organisms, the assumed kinetics are all the same, namely single-substrate Michaelis-Menten uptake kinetics (with an optional Hill coefficient) and constant reaction rate bounds. More complicated kinetics, or multi-substrate dependencies, or environmental influences (such as effects of temperature or pH) are simply not possible. We have now rephrased our conclusions on COMETS to the following: “Hence, while COMETS sets an important precedent, considerable work is still needed to make DFBA a practical approach in microbial ecosystem modeling.” We have also adjusted the remaining text in that paragraph (last paragraph of the “Model” section), where we explain the main distinguishing features of MCM.

Regarding MCM's untested features:

At the end of the day MCM is a framework facilitating the construction, analysis and calibration of microbial ecosystem models. MCM's versatility should thus be understood in terms of possible model structures and possible couplings between components (e.g. environment-cell coupling). The actual model will always be subject to the user's choice. Hence, for example, saying that MCM “supports phage predation” is merely a claim that MCM supports model structures that have conventionally been used to model phage predation (at the population level; for an example see Jensen et al., 2006). Similarly, when we say “MCM supports reaction energetics” this means that MCM allows for the calculation of free energies and their incorporation into models (e.g. influencing growth factors), regardless of how this is used in any particular case (for an example see Reed et al., 2014).

We now make this clear in the “MCM overview” section, saying: “Finally, we emphasize that MCM is, after all, merely a framework enabling the construction, calibration and analysis of microbial community models. MCM models are thus limited by the same caveats and assumptions as other constraint-based metabolic models (Blazier and Papin, 2012; Antoniewicz, 2013) and any predictions made by MCM should be subject to similar scrutiny.”

Furthermore, to avoid confusion, we omit several of MCM's features from the main text. For the ones kept, we refer to existing literature that demonstrates their potential use.

Reviewer 3:

First off, I strongly agree with Reviewer 2's comments above. They have done much to improve the paper to show more and fix much of the language. As noted, however, saying why this matters in the paper as well as they did in the rebuttal would be great, and they still need to back down from hyping a large series of bells and whistles that may be great features but are not demonstrated here. There are plenty of nice advances in their method and they should stick to these.

We have now further emphasized the significance of our experimental system in the Introduction. We also rephrased the comparison of MCM to previous tools in a more conservative way. Please see our responses to Reviewer 2 for details.

As for my own concerns, the direct comparison to data and clarification that expression and flux should not be expected to perfectly correlate is a massive relief. The qualitative agreement is not bad in Figure 6c and d, and can rest on its own laurels. The figures are also particularly improved in their clarity.

I still find the concept of an “average flux” over such a dynamic experiment to be strange, but at least it is now explained. Using the same time point as the RNA data is far better.

For the comparison of predicted pathway activations to transcript microarray profiles (Figure 6c and d) we do indeed use the same time points as the RNA data. We now make that more clear by ending the figure caption as follows: “Flux predictions correspond to the time points of mRNA measurements, i.e. 3.5 hours after dilution for SS and 4 hours after dilution for A and SS (Le Gac et al., 2008).”

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

Article and author information

Author details

  1. Stilianos Louca

    Institute of Applied Mathematics, University of British Columbia, Vancouver, Canada
    Contribution
    SL, Conceived and developed MCM. Designed the E. coli community model and performed the calibration and simulations. Analyzed and interpreted the simulations. Wrote the article
    For correspondence
    stilianos.louca@gmail.com
    Competing interests
    No competing interests declared.
  2. Michael Doebeli

    Department of Zoology, University of British Columbia, Vancouver, Canada
    Contribution
    MD, Interpreted the simulations. Contributed to writing the article. Supervised the project
    Competing interests
    MD: Reviewing editor, eLife.

Funding

Natural Sciences and Engineering Research Council of Canada (NSERC) (Discovery Grant, 21990)

  • Michael Doebeli

Pacific Institute for Mathematical Sciences (PIMS) (International Graduate Training Centre (IGTC) in Mathematical Biology)

  • Stilianos Louca

University of British Columbia (UBC) (Graduate Student Fellowship)

  • Stilianos Louca

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

Acknowledgements

SL acknowledges the financial support of the PIMS IGTC for Mathematical Biology as well as the Department of Mathematics, UBC. SL and MD acknowledge the support of NSERC. We thank Evan Durno and Matthew Osmond for comments. We thank Sung Ho Yoon (Korea Research Institute of Bioscience & Biotechnology) for providing us with a copy of the REL606 cell-metabolic model (Yoon et al., 2012). We thank Mickaël Le Gac (Laboratoire d'Ecologie Pélagique, France) for providing us with the raw data of his monoculture incubation experiments (Le Gac et al., 2008).

Reviewing Editor

  1. Wenying Shou, Reviewing Editor, Fred Hutchinson Cancer Research Center,, United States

Publication history

  1. Received: April 18, 2015
  2. Accepted: September 17, 2015
  3. Version of Record published: October 16, 2015 (version 1)

Copyright

© 2015, Louca and Doebeli

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

Metrics

  • 1,358
    Page views
  • 329
    Downloads
  • 12
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Comments

Download links

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

Downloads (link to download the article as PDF)

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

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

Further reading

    1. Neuroscience
    Sebastian A Vasquez-Lopez et al.
    Short Report Updated