- 1Centre for Computational Biology, Institute for Microbiology and Infection, School of Biosciences, University of Birmingham, Birmingham, United Kingdom
- 2Laboratory of Microbiology, Wageningen University and Research, Wageningen, Netherlands
- 3Department of Physics, School of Agricultural Engineering of Barcelona, Universitat Politècnica de Catalunya–BarcelonaTech, Castelldefels, Spain
- 4Department of Plant Pathology, University of California, Davis, Davis, CA, United States
- 5Laboratory of Synthetic Microbiology, Key Laboratory of Systems Bioengineering (Ministry of Education), School of Chemical Engineering and Technology, Tianjin University, Tianjin, China
- 6Civil and Environmental Engineering Department, Marine and Environmental Sciences Department, Bioengineering Department, Northeastern University, Boston, MA, United States
Models are important tools in microbial ecology. They can be used to advance understanding by helping to interpret observations and test hypotheses, and to predict the effects of ecosystem management actions or a different climate. Over the past decades, biological knowledge and ecosystem observations have advanced to the molecular and in particular gene level. However, microbial ecology models have changed less and a current challenge is to make them utilize the knowledge and observations at the genetic level. We review published models that explicitly consider genes and make predictions at the population or ecosystem level. The models can be grouped into three general approaches, i.e., metabolic flux, gene-centric and agent-based. We describe and contrast these approaches by applying them to a hypothetical ecosystem and discuss their strengths and weaknesses. An important distinguishing feature is how variation between individual cells (individuality) is handled. In microbial ecosystems, individual heterogeneity is generated by a number of mechanisms including stochastic interactions of molecules (e.g., gene expression), stochastic and deterministic cell division asymmetry, small-scale environmental heterogeneity, and differential transport in a heterogeneous environment. This heterogeneity can then be amplified and transferred to other cell properties by several mechanisms, including nutrient uptake, metabolism and growth, cell cycle asynchronicity and the effects of age and damage. For example, stochastic gene expression may lead to heterogeneity in nutrient uptake enzyme levels, which in turn results in heterogeneity in intracellular nutrient levels. Individuality can have important ecological consequences, including division of labor, bet hedging, aging and sub-optimality. Understanding the importance of individuality and the mechanism(s) underlying it for the specific microbial system and question investigated is essential for selecting the optimal modeling strategy.
Introduction
Microbes are important drivers of biogeochemical cycles in all ecosystems and impact their environments in a plethora of ways. For example, in lakes, the harmful cyanobacterium Microcystis aeruginosa can bloom and produce toxins that make the water unsafe to drink (Paerl et al., 2011). The common gut bacterium Bacteroides fragilis produces a chemical that helps the host develop its immune system (Atarashi et al., 2013).
Models are important tools for understanding and managing ecosystems. They can be used to advance scientific understanding by interpreting field observations and aid in hypothesis testing. For example, Jöhnk et al. (2008) used a model to quantify the roles of temperature range and buoyancy regulation in the fitness of the toxic cyanobacterium Microcystis during heat waves. Buffie et al. (2015) applied the model of Stein et al. (2013) to infer an antagonistic interaction in the gut between the pathogen Clostridium difficile and another species of that genus, Clostridium scindens. For ecosystem management, models can be used to answer “what if” questions and make predictions about the effects of future environmental conditions. For example, Blumberg and Di Toro (1990) used a model to predict the effects of climate warming on phytoplankton and dissolved oxygen in a lake. Bucci et al. (2016) predicted the composition of the mouse gut microbiota following infection with C. difficile.
In the past decades, microbiology has experienced rapid advances in observational and experimental technologies, resulting in substantial progress in the understanding of microbes at the molecular level. For example, nitrogen (N) fixation by the cyanobacterium Anabaena involves a division of labor among N-fixing heterocysts and photosynthesizing vegetative cells. The nitrogen-containing β-aspartyl-arginine is produced by cyanophycinase in heterocysts, transferred intercellularly to vegetative cells where it is converted to aspartate and arginine by isoaspartyl dipeptidase (Burnat et al., 2014). Another example involves transcription of genes to messenger RNA (mRNA) and translation to proteins, which is performed by RNA polymerase (RNAP) and the ribosome complex, respectively. In bacteria, those can form a single transcribing and translating “expressome” complex, with known 3D structure and functional consequences on transcriptional pausing, backtracking and termination (Kohler et al., 2017). Characterization of ecosystems is following the same trend. For example, lakes used to be characterized using bulk measures, like Chlorophyll a and total phosphorus concentrations, but observations are now often at the molecular level, including gene expression (transcript levels) (Vila-Costa et al., 2013). Animal and human microbiota are now routinely characterized using multiple omics technologies, such as community characterization using bacterial 16S ribosomal RNA (rRNA) polymerase chain reaction (Costello et al., 2009), and increasingly meta-genomics, transcriptomics and proteomics (Wang et al., 2015).
The development of models is lagging behind as most models still do not make use of molecular level understanding or observations. It is recognized that there is a substantial gap between our microbial ecology models and current microbiology knowledge and environmental observations (Fuhrman et al., 2013; Trivedi et al., 2013; Hellweger, 2015; Dick, 2017; Stec et al., 2017). For example, lake phytoplankton models still simulate phytoplankton biomass concentrations (e.g., μg Chlorophyll a L−1) and the effect of a nutrient on the growth rate using an equation developed 75 years ago (Monod model). Likewise, most models of the gut aggregate species into functional groups based on metabolic pathways (Kettle et al., 2015). Models are now being developed that explicitly resolve genes and make predictions at the population and ecosystem level.
This paper has two parts. First, we review existing modeling approaches. Here, we focus on mechanistic models that explicitly include genes and simulate population-level properties (e.g., microbe concentration, nutrient uptake) rather than empirical models. One aspect in which the existing approaches differ is their representation of microbial individuality. The second part of our review will use examples to explain why including individuality is important.
Part 1: Review of Existing Modeling Approaches
In this section we describe three modeling approaches that have been used to bridge the gap between genes and ecosystems, including metabolic flux, gene-centric and agent-based modeling (ABM). We illustrate each approach using a hypothetical ecosystem, where two microbial species grow and interact via three metabolites (Figure 1). We then discuss a number of examples from the literature, focusing mostly on the modeling aspects of the studies. Then we highlight the weaknesses and strengths of each approach. Finally, we characterize the models along a number of dimensions, including space, time, function, heterogeneity, species diversity and genes.
Figure 1. Hypothetical ecosystem used to illustrate different modeling approaches. Species 1 takes up metabolites A and B, produces metabolites C, D, E, and F and excretes metabolite C. Species 2 takes up metabolites A and C and produces metabolites D, G, and H.
Literature Selection Criteria
The review is focused on the use of gene-level models for advancing understanding and making predictions of microbial ecosystems. To keep the scope of the review manageable, we included only quantitative models, which are required for predictions, although qualitative models may be sufficient to advance understanding. We applied the following selection criteria: (1) model uses a mechanistic (vs. empirical) approach, (2) model explicitly considers at least one actual gene or protein; (3) model includes some form of direct or indirect interaction among microbes; (4) model includes multiple microbial species (or strains) or phenotypes in different locations; and (5) model makes predictions at the population level. We therefore exclude empirical models that correlate observed gene distributions to environmental factors and function (e.g., carbon export in the ocean, Guidi et al., 2016), models that use hypothetical genes or digital genomes describing behavioral traits (e.g., Lenski et al., 1999; Clark et al., 2011), scale up single-cell models using multiple independent simulations where the cells do not interact (e.g., Emonet and Cluzel, 2008; Labhsetwar et al., 2013) and studies that infer interactions from comparison of metabolic networks and do not make quantitative predictions (e.g., Levy and Borenstein, 2013; Zelezniak et al., 2015).
Metabolic Flux Modeling
Definition
This approach builds on the genome-scale, constraint-based modeling approach most commonly applied to single species (Feist et al., 2008). In this approach, the genome sequence is used to derive a network of potential metabolic reactions by a combination of automated and manual (curation) steps. Then, a flux distribution is predicted, typically by optimizing the flux distribution to maximize an objective function, like maximization of biomass production (Schuster et al., 2008). The extension of this approach to multiple species builds on efforts to extend it to multiple compartments of higher eukaryotic organisms. There are three approaches to multi-species metabolic flux modeling, which we will refer to as environmentally coupled, directly linked and aggregated approaches. The environmentally coupled approach builds on the dynamic flux balance analysis (FBA) approach (Varma and Palsson, 1994), where the microbes and extracellular metabolites are represented using concentration state variables. The growth rate and metabolite fluxes are computed from FBA assuming a common pool for extracellular metabolites and that the system is in a steady state during each time step. The directly linked approach explicitly links the metabolic networks of the species using exchange reactions. This is conceptually the same way in which multi-compartment organisms are modeled. The aggregated approach (also referred to as pooled, supra-organism or enzyme soup approach) involves constructing one network by combining the individual networks and removing duplicates. This ignores cellular boundaries and is most applicable to metagenomic datasets. Box 1 illustrates these three approaches for the hypothetical ecosystem shown in Figure 1. This approach has also been referred to as Ecosystems Biology (Klitgord and Segrè, 2011) or Community Systems Biology (Zengler and Palsson, 2012).
Box 1. From genes to ecosystems using metabolic flux modeling.
Single species
The starting point for a metabolic flux model is a set of mass balance equations:
where x (mmol gDW−1, i.e., per gram biomass dry weight) is a vector of metabolite concentrations, S is the stoichiometric matrix and v (mmol gDW−1 h−1) is a vector of reaction rates for uptake, excretion, internal metabolism and growth. Typically, a steady-state is assumed so the derivatives are zero. The stoichiometric matrix (S) for species 1 of the hypothetical ecosystem is presented in Table B1, where columns are reactions and rows are metabolites. Lower and upper bounds for the reaction rates, determined based on thermodynamics, enzyme kinetics or measurements, can be included in the optimization procedure.
There are typically infinitely many solutions that satisfy the equation. For example, in species 1 (Figure B1.1), biomass (metabolite F) can be produced by any combination of two pathways (AEF or ADF). Computational methods are available that decompose the stoichiometric matrix into unique sets of functional units (pathways) such as elementary modes or extreme pathways (Papin et al., 2004). A more common approach, flux-balance analysis (FBA), involves optimizing reaction rates to maximize the value of some objective function using linear programing (LP). Several objective functions have been used, such as minimizing ATP production and maximizing production of some metabolite, but maximizing biomass production yield or rate is often considered to be the most appropriate in an ecological context. When biomass production is maximized, it is assumed that the cell regulates fluxes through its metabolic network in a way that maximizes biomass production. The corresponding objective function for the species 1 of our example system is to maximize the production of metabolite F (VMetE1 + VPrdF1 or VGrowth1 in Table B1). This is relatively simple and real models typically use a more complex biomass growth function, e.g., a genome-scale model may include various precursors (e.g., G6P, F6P) and cofactors (e.g., ATP, NADH). Algorithms that integrate gene expression data are also available (Becker and Palsson, 2008). FBA is fundamentally a steady-state approach, but a pseudo-time-variable model can be constructed (Varma and Palsson, 1994; Mahadevan et al., 2002).
Multiple species—environmentally coupled models
Figure B1.1 illustrates the dynamic, multi-species metabolic flux modeling methodology. The model includes state variables for microbial biomass (X) and extracellular metabolites (C). The microbes grow according to a growth rate (μ) and consume/produce metabolites according to specific consumption/production rates (V). Those values are calculated from the metabolic flux models, which are optimized to maximize the growth yield subject to a number of constraints, including a maximum consumption rate for each metabolite based on its concentration. A simulation will proceed in a step-wise manner: (1) Calculate the constraints based on all C. (2) Optimize the metabolic model of each species, which yields μ and V. (3) Calculate the new X for both species based on μ. (4) Calculate the new C for both metabolites based on V from both species. Repeat. When the metabolic model does not lead to a viable solution, a simple death routine can be invoked (Zhuang et al., 2011). It is conceptually straightforward to include other reactions (e.g., between extracellular compounds) and transport (Scheibe et al., 2009).
Figure B1.1. Multi-species metabolic flux modeling—environmentally coupled models. After Figure 2 in Zhuang et al. (2011). X (gDW L−1) = microbial biomass concentration, C (mmol L−1) = extracellular metabolite concentrations, μ (h−1) = specific growth rates, V (mmol gDW−1 h−1) = specific flux velocities.
Multiple species—directly linked model
Figure B1.2 illustrates the multi-species metabolic flux modeling approach developed by Stolyar et al. (2007). The metabolic models for each species (Figure B1.1) are combined into a single model. Exchange of metabolites among the species occurs by directly linking their reactions, which constrains them to be the same. This is equivalent to assuming there is no change in the extracellular metabolite concentrations. The model is optimized to maximize a weighted combination of the biomasses.
Multiple species—aggregated model
Figure B1.3 illustrates the multi-species metabolic flux modeling approach developed by Taffs et al. (2009). The reactions and metabolites for the two species (as shown in Figure B1.1) are merged into a single model and a single objective function is used to determine the flux distribution.
Figure B1.2. Multi-species metabolic flux modeling—directly linked model. After Figure 2 in Stolyar et al. (2007). The metabolic models for each species (Figure B1.1) are combined into one model, with exchange reactions linking their metabolisms.
Figure B1.3. Multi-species metabolic flux modeling—aggregated model. After Figure S2 in Taffs et al. (2009). The metabolic models for each species (as shown in Figure B1.2) are merged into one model.
Examples
There have been several applications of metabolic flux models to communities of microbes. For recent reviews see Zengler and Palsson (2012), Biggs et al. (2015), Tan et al. (2015), Zomorrodi and Segrè (2016), Perez-Garcia et al. (2016), and Gottstein et al. (2016).
Environmentally coupled models
Scheibe et al. (2009) applied FBA to learn about the growth of Geobacter and uranium bioremediation in a contaminated groundwater site where Geobacter dominates the community. They coupled a genome-scale FBA model to a two-dimensional reactive transport model. The FBA model computes growth rate and fluxes based on ambient acetate, Fe(III) and ammonia concentrations in each grid element. Those growth rates and fluxes are then used by the reactive transport model to compute the Geobacter biomass, acetate, Fe(III) and ammonia concentrations, as well as other processes like U(VI) reduction. The new ambient concentrations are then again used by the FBA model to compute the growth rate and fluxes at the next time step and so on. Due to computational constraints, the FBA calculations were done a priori for 1,000 combinations of metabolite concentrations and stored in a look-up table, rather than a dynamic coupling between the models. One of the main advantages of the FBA-based approach is that it allows for variable substrate utilization and growth yields, which is not supported by conventional models. The model was able to make predictions of similar quality as the previous reactive transport model (i.e., without FBA component), but it did so without the need to calibrate rate parameters (Figure 2).
Figure 2. Comparison of observations (symbols), traditional model (solid lines) and FBA-based model (dashed lines). Reproduced from Scheibe et al. (2009) with permission. The figure shows acetate and U(VI) concentrations at a groundwater bioremediation site. Concentration time series are presented at 3.7, 7.3, and 14.6 m distance from the acetate injection gallery. Acetate increases at progressively later times as the distance from the injection gallery increases. Consistent with this, U(VI) decreases at progressively later times. Colors identify single wells.
Tzamali et al. (2009) and Tzamali et al. (2011) used the dynamic FBA approach to simulate the interaction among various E. coli strains, including wild type and single gene knockouts. For various substrates, they identified potential communities of co-existing strains. For example, growth on pyruvate supported communities with up to 6 strains. The most efficient community of 4 mutants produced 2.2% more biomass than a pure culture of the wild type.
Zhuang et al. (2011) developed a dynamic, genome-scale FBA model of two species in competition in a uranium-contaminated aquifer. Rhodoferax and Geobacter both oxidize acetate and reduce Fe(III), but only Geobacter can reduce U(VI), rendering it less soluble and therefore contributing to the clean-up of the site. The FBA models of the two species calculate growth and metabolite production/consumption rates, which are used to integrate biomass and metabolite concentration state variables. The model predicted that, under low-ammonia conditions, Rhodoferax is outcompeted by Geobacter, which can fix nitrogen, and that this promotes respiration (vs. biomass production) and associated U(VI) reduction, which are patterns consistent with observations.
Zhuang et al. (2012) expanded the model by Zhuang et al. (2011) and applied it to design remediation scenarios. In particular, they used two separate FBA models for attached and planktonic Geobacter to differentiate their functions: planktonic cells reduce U(VI) and attached cells reduce Fe(III). Attachment and detachment rates were used to transfer biomass among these two fractions. This illustrates one approach by which heterogeneity can be simulated in these types of models.
Harcombe et al. (2014) developed dynamic FBA models of two and three species on a two-dimensional grid, where biomass grows and dies, extracellular metabolites are consumed and produced, and biomass and metabolites move by diffusion. Cole et al. (2015) extended the dynamic FBA approach further to three dimensions and used it to simulate growth of E. coli in colonies on agar. The model was able to simulate the small-scale environmental heterogeneity in dissolved oxygen and nutrient concentrations, and the resulting phenotypic differentiation of the bacteria (i.e., fermenting cells in the interior). Other multi-species, environmentally coupled metabolic flux models were presented by Salimi et al. (2010), Hanly and Henson (2011), Hanly and Henson (2013), Biggs and Papin (2013), Chiu et al. (2014) and Louca and Doebeli (2015). Zomorrodi et al. (2014) presented a dynamic version of the multi-level optimization routine presented previously (Zomorrodi and Maranas, 2012, see below).
Directly linked models
Stolyar et al. (2007) developed an FBA model of two microbes that are mutualistic in the absence of sulfate, Desulfovibrio vulgaris and Methanococcus maripaludis. In the scenario evaluated, D. vulgaris grows on lactate, producing H2, formate, CO2 and acetate, which support the growth of M. maripaludis. The model consists of three compartments, representing the metabolism of the two species and the exchange between them. The metabolite fluxes in the central metabolism of each species and exchange reactions are represented using 89 and 82 equations, respectively. The third compartment represents the exchange flux of H2, formate, CO2 and acetate, where H2 and formate were not allowed to accumulate in the medium, so that their rates of production by D. vulgaris and consumption by M. maripaludis are the same. The combined model was optimized to maximize biomass production of both species, with a larger weight for D. vulgaris, based on observations. However, the biomass ratio of the two species is constrained by the exchange reaction, so it was relatively invariant to the weights used. The model suggested that the H2 was essential, but that formate could be eliminated.
Wintermute and Silver (2010) applied the FBA modeling approach at the genome scale to 46 E. coli mutants, each incapable to synthesize an essential metabolite. Growth experiments were conducted with 1,035 binary strain combinations. A joint FBA of each pair allowing for exchange of all shared metabolites between the strains was developed. The models were optimized to minimize the difference between the flux distributions of the wildtype and mutant (minimization of metabolic adjustment, MOMA, Segrè et al., 2002). The idea behind this objective function is that the regulatory system is still based on the wildtype and has not yet adjusted to the mutation. The joint FBA models were consistent with the finding that pairings of mutants blocked in the same biosynthetic pathway rarely show synergistic growth (4% of the cases) while pairings of mutants in separate pathways did so in 18% of cases. The model correctly predicted that strains grow best when they require small amounts of metabolites that are cheap to produce by the other strain. The ability of simple stoichiometric models to predict fitness costs and benefits of metabolic cross-feeding is encouraging.
Klitgord and Segrè (2010) applied the FBA modeling approach to binary pairs of seven species and identified the media composition that would support symbiosis. They developed genome-scale FBA models of all possible binary pairs and did a systematic search for media compositions that would support growth of the pair but not the individual species.
Huthmacher et al. (2010) generated an FBA model of the metabolism of the malaria causing Plasmodium falciparum and its host, the erythrocyte (red blood cell). By constraining the metabolic network with gene expression data of P. falciparum, they were able to predict metabolic fluxes for different life cycle stages of the pathogen.
Zomorrodi and Maranas (2012) developed a community FBA modeling framework and applied it to a number of systems, including those of Stolyar et al. (2007) and Taffs et al. (2009). A novel aspect in this work is the consideration of multiple objective functions, including maximization of growth of each species as well as biomass production at the community level, which can be used to explore tradeoffs between selfish and altruistic driving forces.
Other multi-species, directly linked metabolic flux models were produced by Taffs et al. (2009), Bizukojc et al. (2010), Bordbar et al. (2010), Freilich et al. (2011), Khandelwal et al. (2013), Nagarajan et al. (2013), Shoaie et al. (2013), Ye et al. (2014), El-Semman et al. (2014), Merino et al. (2015) and Heinken and Thiele (2015).
Aggregated models
Taffs et al. (2009) applied different approaches to model three species (oxygenic phototrophs, filamentous anoxygenic phototrophs and sulfate-reducing bacteria) in the thermophilic, phototrophic mat communities from Octopus and Mushroom Springs in Yellowstone National Park (USA). One of their approaches does not consider compartments, but lumps all reactions into one species (see Box 1). This approach ignores compartmentalization and the fact that intermediate intracellular metabolites from one species may not be available to another. However, it does not require assigning individual enzymes or reactions to species, functional groups or guilds and is well suited for data from metagenomics. A unique aspect of this study is the use of elementary mode analysis (EMA), which is an alternative to FBA and characterizes the set of all possible flux distributions, rather than just the optimal one.
Tobalina et al. (2015) applied the aggregated approach to naphthalene-contaminated soil communities. An interesting aspect of that study was that the model was based on metaproteomics data, which implicitly accounts for regulation.
Cerqueda-García and Falcón (2016) applied the aggregated approach to study the metabolism of communities in microbial mats and microbialites (living carbonate rock structures similar to corals and stromatolites). Starting with metagenomic datasets, they reconstructed a metabolic network, and then used EMA to identify feasible pathways through this network for C and N assimilation. They identified a number of alternative CO2 fixation pathways, which were not identified for these systems previously.
Strengths
• The FBA approach can directly utilize molecular data, genomics, transcriptomics, proteomics and metabolomics, from pure laboratory cultures and the environment (e.g., metagenomics) as long as annotations are available, which is increasingly the case.
• The approach is comprehensive in terms of functions and metabolites. This is likely to be increasingly useful, as recent observations from a number of environments suggest that bacteria have a high substrate specificity (Kindaichi et al., 2013; Salcher et al., 2013). For example, when a freshwater community was presented with 14 radiolabeled low-molecular weight (LMW) organic substrates, the two most abundant microbes belonging to the Actinobacteria ac1 and Alphaproteobacteria LD12 tribes had no overlap in their substrate acquisition spectra. The concept of dissolved organic carbon (DOC) as a common currency for heterotrophic microbes is too simplistic. One of the main applications of FBA has been to understand complex substrate uptake patterns.
Weaknesses
• The directly linked and aggregated approaches assume the system to be in a steady-state. The environmentally coupled approach also assumes steady-state flux distributions during each time step, but flux distributions can change from time to time. For many cases this assumption will be sensible, but for others not. For example, planktonic bacteria experience a very heterogeneous nutrient regime and may experience nutrient patches with short durations (~60 s, Taylor and Stocker, 2012), comparable to the time required for gene expression, protein translation and maturation. Genome-scale models are being developed that go beyond steady-state metabolite fluxes (e.g., include dynamic transcript, protein and metabolite pools, Karr et al., 2012) and this technology will eventually be applied at the ecosystem scale.
• It is not always clear what objective function should be used to optimize the flux distribution (Schuster et al., 2008). Maximization of biomass production seems like a good choice from a biotechnological perspective. However, there are cases where it is advantageous to divert production away from biomass, including to storage products, toxins or EPS (Merino et al., 2015), which may conflict with the biomass objective. Moreover, in a well-mixed, stable environment, specific growth rate will likely be maximized by natural selection while in a spatially structured environment such as a biofilm, the biomass yield is likely to be maximized by natural selection (Kreft, 2004).
• The approach typically entails specifying a biomass composition, and commonly this is applied across different conditions. However, the biomass composition is known to change (Benyamini et al., 2010).
• Growth dilution of metabolites, other than the ones used in the growth equation (see above), is typically ignored (Benyamini et al., 2010). Specifically, there should be a “–μ x” on the right-hand side of Equation B1.1. Accounting for growth dilution is conceptually straightforward but it requires specifying the metabolite concentrations, which are not typically available at the genome scale. Metabolomics data can help to fill this gap, but this would be difficult for all metabolites, times and locations in the model and impossible for prediction simulations. Another hurdle is the computational cost. The metabolite dilution FBA (MD-FBA) model of Benyamini et al. (2010) uses mixed-integer linear programming (MILP, vs. LP used by FBA), which is computationally more demanding than LP. This limitation may be especially important for applications that require solutions for multiple species, times and locations.
• The approach does not account for individual heterogeneity (see Part 2).
Gene-Centric Modeling
Definition
In the gene-centric or functional gene approach, the model is built based on genetic information, as in metabolic flux modeling, but focused on capturing the dynamic behavior of specific genes or gene activities in the system. Thus, the biogeochemical fluxes are based on the genetic composition of the microbial community. Microbes are grouped based on specific functional or proxy genes and tracked using corresponding concentration state variables. This is similar in spirit to modeling functional groups (e.g., N-fixers, lactate producers, Le Quéré et al., 2005; Kettle et al., 2015). The concentrations of genes (e.g., number of gene copies per liter) are simulated using mass balance differential equations, which is how typical microbial ecology models simulate species. The rate of gene production (or growth) can be tied to the Gibbs free energy released by the reaction catalyzed by the corresponding enzyme. The approach is illustrated in Box 2.
Box 2. From genes to ecosystems using gene centric modeling.
The following description is based on Reed et al. (2014), but adapted to the hypothetical ecosystem considered here and some of the nomenclature is altered to facilitate comparison with the other approaches. The first step in the development of the model is to identify the functional genes. For the hypothetical example, we will use prdF and metC and ignore the others (Figure B2). Thus, state variables prdFg and metCg (no. L−1) represent species 1 and 2, respectively. The reactions mediated by these gene products are assumed to be the limiting reactions along the pathway, but exchange with extracellular metabolites requires accounting for the input and output of the entire pathway. For prdF, the overall reaction is:
The production rate of metC genes as a result of metabolism associated with the metC gene (RmetC, no. L−1 d−1) is:
metCg (no. L−1) is the metC gene concentration. FT,metC is a thermodynamic potential factor, which accounts for the chemical energy available to drive the metabolism, and can be estimated from the energy yield of the associated reaction. μmax,metC (d−1) is the maximum specific growth rate. CC (molC L−1) and Km,metC,C (molC L−1) are the extracellular concentration and half-saturation constant for metabolite C. The production rate of prdF genes (RprdF) is:
The last fraction accounts for inhibition by metabolite C. The mass balance equation for extracellular metabolite B is (transport and other reactions are omitted for clarity):
γprdF,B and γprdF,A are stoichiometric coefficients. Here, 2 mol B are consumed for every 1 mol A, so γprdF,B = 2 and γprdF,A = 1. qprdF (no. gDW−1, i.e., per gram biomass dry weight) is the intracellular concentration of prdF genes, which depends on the number of gene copies in the genome. YprdF (gDW molA−1) is the yield, which depends on the energy yield of the associated reaction. For metabolite A, we have to consider consumption by both species:
Here, A and C are consumed in equal amounts, so γmetC,A = γsmetC,C = 1. Metabolite C is consumed by the reaction associated with metC and produced by the reaction associated with prdF:
Here, 2 mol C are produced for every 1 mol A consumed, so γprdF,C = 2 and γprdF,A = 1. The mass balance for gene prdF is:
where kd (d−1) is the death rate. The first term in the gene mass balance equation accounts for the production of the gene due to the growth associated with its reaction and the second term accounts for mortality.
Examples
Reed et al. (2014) presented the gene-centric approach and applied it to study nitrogen cycling in the Arabian Sea oxygen minimum zone (OMZ). The model includes eight functional genes, including those for denitrification (nitrate reductase, narG, nitrite reductase, nirK), aerobic ammonia oxidation (ammonia monooxygenase, amoA) and anaerobic ammonium oxidation (anammox, hydrazine oxidoreductase, hzo), as well as relevant metabolites, including dissolved oxygen (O2), ammonium (), nitrate (), and nitrite (). The model-predicted gene abundances were compared directly to observations from qPCR (gene copies L−1, Figure 3). The authors also compared model-predicted changes in gene abundances over time to observed mRNA concentrations in a qualitative manner (gene copies L−1 s−1 vs. mRNA copies L−1, Figure 3). An interesting problem addressed by this model is the dual role of nitrogen as an energy source and biomass component, where the latter is not considered by the gene-centric approach. This was handled by calculating the total biomass increase/decrease and removing/adding corresponding amounts of N from/to the extracellular metabolite pools. The model was used to show that denitrification is the dominant nitrogen loss process in this area, which is different from many other OMZs, where anammox dominates.
Figure 3. Comparison of gene-centric model predictions (solid lines) to observations (dotted line in top left and symbols in other panels) of metabolites, gene and mRNA levels in the Arabian Sea oxygen minimum zone. From Reed et al. (2014). Copyright © (2014) by the National Academy of Sciences. The figure shows an oxygen low from about 400 to 900 m and a coincident low in ammonia monooxygenase (amoA) DNA and mRNA, and an increase in anaerobic ammonium oxidation (anammox, hydrazine oxidoreductase, hzo) DNA and mRNA.
Reed et al. (2015) applied the gene-centric approach to simulate functional genes for sulfur, nitrite, ammonia, methane and hydrogen oxidation and associated metabolites in a submarine hydrothermal vent plume. The pathways for oxidation of a number of reduced sulfur species (e.g., hydrogen sulfide, thiosulphate) co-occur in one species (SUP05), so the pathways were combined into one functional gene in the model. The authors compared their model-predicted relative abundance (%) of functional genes to observations. The hydrogen concentration is relatively low in this system and their model predicted no significant increase in hydrogenase gene abundance due to aerobic hydrogen oxidation. However, substantial quantities of hydrogenase genes were observed suggesting that they may be produced because they co-occur with a gene that does experience substantial growth as observed (SUP05 may have genes for sulfur and hydrogen oxidation, Anantharaman et al., 2013). When this coupling is included in the model, it was able to reproduce the observations.
Louca et al. (2016) presented a gene-centric model of six functional genes and eight metabolites for a number of dissimilatory redox pathways involved in nitrogen and sulfur cycling in a seasonally anoxic fjord (Saanich Inlet, Vancouver Island, Canada). That model extends the gene-centric modeling approach by explicitly simulating mRNA and proteins, assuming their production rates are proportional to the corresponding reaction rates and subjecting them to transport and decay processes. Model predictions for mRNA and proteins were compared to observations on a qualitative basis. The model was used to gain insights into the sulfur and nitrogen pathways in this system. For example, the model predicted incomplete denitrification by the SUP05 clade, which results in leakage of nitrite that supports anammox and loss of nitrogen.
Strengths
• The approach is readily integrated into existing models based on concentration state variables (Reed et al., 2014).
• The approach makes quantitative predictions of gene levels that can be compared directly to observations.
Weaknesses
• While this modeling approach is readily applied to chemotrophs where there is a direct link between the rate of the reaction and the growth rate of the microbe, it is less clear how to apply it to a phytoplankton species that may be limited by nitrate, but uses the energy derived from sunlight to reduce it to ammonia for incorporation into amino acids. Heterotrophs growing on a complex mixture of dissolved organic matter (DOM) may also be difficult to model with this approach.
• The extension of the method to mRNA and proteins (Louca et al., 2016) includes simulating them as independent concentration variables. This does not account for their natural co-existence in the cell and may lead to some odd effects, like mRNA and protein appearing in locations where there are no corresponding genes.
• The method supports multiple co-occurring genes (see Reed et al., 2014 for equations), but that is based on constant fractions within a community, which may change dynamically and spatially in a natural community (e.g., species succession in phytoplankton). This method is also more difficult to implement. The reader is invited to rework the example in Box 2 for a model that uses metA and metC (which co-occur in species 1) as functional genes.
• The approach does not account for individual heterogeneity (see Part 2).
Agent-Based Modeling (ABM)
Definition
ABM or individual-based modeling (IBM) involves simulating microbes as individuals. This is in contrast to the traditional population-level approach, where microbes are simulated as concentration state variables. ABM is already an established modeling technology in microbial ecology (Hellweger et al., 2016a). Microbial ABMs increasingly resolve intracellular mechanisms and the extension to genes is a natural progression and already well on the way. The approach is illustrated in Box 3. This approach has also been referred to as Systems BioEcology (Hellweger, 2009).
Box 3. From genes to ecosystems using agent-based modeling.
Here the approach for dynamic, molecular-level, mechanistic modeling is illustrated by application to the hypothetical ecosystem (Figure 1). The model explicitly resolves genes, transcripts, proteins and metabolites (Figure B3). Following the central dogma of biology, genes are transcribed by the RNA polymerase to yield transcripts (mRNA), which are translated by the ribosomes to yield proteins, which then carry out various functions. Once a biomass (in the case of strain 1 the metabolite E, QE) threshold is reached, the DNA polymerase is induced, which synthesizes DNA. Once that is complete, cell division is induced and the cell divides, which involves division of all intracellular components. The approach entails explicitly modeling genes, transcripts and proteins. However, typically only a handful of representative genes are simulated using a coarse-grained approach (Castellanos et al., 2004; Hellweger, 2013).
Figure B3. Agent-based modeling of microbes. Gene/protein: rpoMH/RNAP, RNA polymerase; rptMH/RPT, ribosome; ftsMH/Fts, cell division; polMH/Pol, DNA polymerase; dumMH, dummy (accounts for genes not explicitly considered); uptA/UptA, uptake A; metA/MetA, metabolism A; excC/ExcC, excretion C. Substrates and metabolites (extracellular, concentration C; intracellular, quota Q): CA, substrate A; QH, metabolite H, etc. After Figure 1 of Hellweger (2009).
To illustrate the approach, we present the equations for uptB1 transcription, UptB1 synthesis, UptB1 rate and QB mass balance. Here, intracellular concentrations are defined on a per biomass dry weight (DW) basis, but carbon and volume can also be used. The uptB1 transcript (uptB1t, mol mRNA gDW−1, i.e., per gram biomass dry weight) mass balance equation is:
where kS,T (bp RNAP−1 s−1) is the transcription rate, LDNA (bp) is the total DNA length, RNAP1 (mol protein gDW−1) is the RNA polymerase level, yuptB1 is the uptB1 expression level (may depend on various factors), uptB1g is the number of uptB1 gene copies, kd,T (s−1) is the mRNA decay rate and μg (s−1) is the specific growth rate. The UptB1 (mol gDW−1) protein mass balance is:
where kS,P (nt RPT−1 s−1) is the translation rate, TxL (mol mRNA gDW−1) is the total mRNA, RPT1 (nmol protein gDW−1) is the ribosome level and kd,P,UptB1 (s−1) is the UptB1 decay rate. The UptB1 rate (VUptB1, mol gDW−1 s−1) is:
where kUptB1 (molB molUptB1−1 s−1) is the UptB1 catalytic rate constant, Km,UptB1 (molB L−1) is the half-saturation constant and Ki,UptB1 (molB gDW−1) is the inhibition constant. The intracellular metabolite B (QB, mol gDW−1) mass balance is:
where VPrdF1 (molF gDW−1 s−1) is the PrdF1 reaction rate (the factor 2 accounts for 2 mol B per 1 mol F, see Figure 1).
Examples
ABM was used by Hellweger (2009) to explore the role of photosynthesis genes (psbA, hli) carried by viruses that infect the marine cyanobacterium Prochlorococcus. The idea is that these genes help to maintain the host photosynthesis apparatus during the latent period, increasing energy to support the replication of the virus. The model simulates individual viruses and host cells and explicitly resolves mechanisms of gene expression, protein synthesis, photosynthesis and events associated with infection at the molecular level. The model was calibrated to observations of virus and host gene transcript and protein levels and then used to simulate population dynamics in the water column of the Sargasso Sea. Modeled populations were diverse, including multiple virus types (different combinations of psbA and hli copies) and cells with different light histories, cell cycle phases and infection stages. Using competition experiments between virus strains that have different combinations of psbA and hli, and evolution experiments (i.e., gene packaging error), the model predicted an optimal gene content that matched that of the wild-type.
An ABM of the cyanobacterium Synechococcus and its circadian clock was constructed by Hellweger (2010). The model structure is similar to the Prochlorococcus model described above. A new feature was the explicit simulation of the concentration of proteins with different phosphorylation states and their interaction. The modeled population includes cells at different phases in their cell and circadian cycles and gene expression levels (psbAI luminescence) were compared to observations at the individual level.
Mina et al. (2013) used an ABM of genetically-engineered quorum sensing E. coli cells in a three-dimensional microfluidics chamber. The model simulates a heterogeneous population of individual, motile cells, each with a number of genes (luxI, aiiA, and yemGFP) and associated proteins, which communicate via a diffusible substance. They showed that autoinducer oscillations on the population level do not follow simply from synchronizing single cell oscillations. Single cells can switch between a state of constant signal concentration and oscillations, depending on the parameters of the positive and negative feedback loops in the gene regulatory network. Yet in a population of these cells, only the oscillatory state is stable—once cell density exceeds a threshold.
Hellweger et al. (2014) built an ABM of the yeast Saccharomyces cerevisiae and used it to explore the fitness effect of age-correlated stress resistance. The model explicitly simulates the regulation of the proteins Tsl1 and Tps3, which synthesize the stress protectant trehalose. Their expression is modeled using constant, age-dependent and stochastic terms. The population is diverse consisting of cells in different phases of their cell cycles as well as different ages, damage and Tsl1/Tsl3 expression levels. The modeled heterogeneity was compared to observations obtained using flow cytometry. Comparison of the various expression strategies showed that age-correlated stress resistance can be beneficial under some conditions.
A model of Anabaena—nitrogen interactions was developed by Hellweger et al. (2016b). This model simulates the uptake of various forms of nitrogen and early intracellular assimilation pathways. Uptake and intracellular reactions are mediated by enzymes (e.g., GlnA) and their expression is controlled by a number of regulatory proteins (e.g., NtcA). A novel feature of this model is the explicit simulation of cell differentiation and division of labor. When fixed nitrogen is depleted, the cells become nitrogen-stressed and some differentiate into heterocysts, which are anoxic cells that fix nitrogen and pass the fixed nitrogen to their neighboring vegetative cells (Figure 4C). The model was informed by observations from 269 laboratory experiments from 55 papers published from 1942 to 2014, including transcript levels and enzyme activities (Figures 4A,B). Hellweger et al. (2016b) also applied the model to a hypothetical lake, but validation by comparison to field observations was not performed.
Figure 4. Anabaena—nitrogen interaction model. (A,B) Comparison to observations of transcript levels in filaments (fil., i.e., all cells), vegetative cells (veg.) and heterocysts (het.) (A) and enzyme activities (B) of cells grown under different conditions (data are from Martin-Figueroa et al., 2000). (C) N fluxes for growth on N2. Red, heterocysts; Green, vegetative cells. Numbers are fluxes in pmol N cell−1 d−1. From Hellweger et al. (2016b). The figure shows that heterocysts have higher levels of nitrogenase (nifH) transcripts, higher levels of glutamine synthetase (glnA, GS) transcripts and enzyme levels and lower levels of glutamate synthase (glsF, GOGAT) transcripts and enzyme activities. These observations support the model where N2 is fixed in heterocysts and combined with glutamate (GLU) that is imported from adjacent vegetative cells to yield glutamine (GLU), which is then exported to vegetative cells and further incorporated into labile nitrogen (LN) and structural nitrogen (SN) pools.
Another gene-level ABM was developed by Hellweger (2013) to explore the mechanisms underlying adaptation of E. coli to tetracycline resistance.
Strengths
• The main advantage of ABM is the ability to resolve intra-population heterogeneity. We will discuss the importance of individuality in the second part of this review.
• In ABM, the description of the system is very flexible and not constrained by having to use one specific mathematical formalism (e.g., the stoichiometric matrix of the FBA approach). For example, on/off control of a gene by light can be modeled simply using an “if” statement (if light is on, then turn on gene, otherwise turn it off). It is much more difficult to incorporate this into a stoichiometric matrix or differential mass balance equation.
Weaknesses
• This approach is relatively complex and difficult to apply. Although it can theoretically be extended to the whole-genome scale, past models have focused on a handful of genes, transcripts, proteins and metabolites. This is due to the limited availability of rate formulations and parameters, and the difficulty of calibrating a model with numerous non-linear feedbacks.
Summary of Examples
The models reviewed above are characterized along six different dimensions, including space, time, function, heterogeneity, species diversity, and genes (Figure 5). These dimensions were selected as they highlight differences between the reviewed models, nevertheless, the list is not exhaustive and other dimensions can be used, like types of interactions. The figure illustrates that, together, the population of past models covers the entire space. However, no single model or approach has covered the entire space by itself. The gene centric approach is amenable to space, time and species diversity and those dimensions have been explored in past models. Function and genes dimensions are linked in this approach, and limited because each species is typically associated with only one function. Simulating individual heterogeneity is difficult with the population-level gene centric approach and has not been explored in past models. Metabolic flux models are routinely genome-scale, and past models have pushed the boundaries along other dimensions, including space, time (although only quasi-time-variable), and species diversity. Function is often limited to metabolism and heterogeneity to phenotypes. The agent-based approach is flexible along most dimensions, but limited in terms of gene coverage and models with more than a handful of genes have yet to be developed.
Figure 5. Comparison of modeling approaches along various dimensions. Based on metabolic flux (N = 32), gene-centric (N = 3) and agent-based (N = 6) models included in this review. Dimensions: Space: 0 (i.e., well-mixed reactor), 1, 2 or 3 dimensional; Time: steady-state, quasi-time-variable (e.g., dynamic FBA), time-variable; Function: None, metabolism, + regulation, + division, + additional functions; Heterogeneity (individuality): None, types (i.e., phenotypes), individuals; Species diversity: None, one, two, three or more; Genes: None, a handful, core/central metabolism, whole genome. Symbols represent averages and “error bars” the range between minimum (i.e., dimensions covered by all models) and maximum. For example, agent-based models have been developed with zero to three spatial dimensions, and the average is 0.67.
Part 2: The Importance of Individuality
A key distinction between the modeling approaches reviewed above is their consideration of individuality and heterogeneity. It is now well established that microbial populations in the environment and laboratory exhibit substantial heterogeneity in properties and behaviors. There are probably cases where this individuality averages out and is of no consequence to the ecology or biogeochemistry of the system (e.g., steady-state growth of a single species on a single nutrient). However, there are also cases where individuality has been shown to critically affect the fitness of a population. A thorough understanding of individual heterogeneity and its potential ecological consequences is critical for selecting the most appropriate modeling strategy.
In this part, we review the importance of individuality. There are a number of mechanisms that produce heterogeneity in a population, which we refer to as sources. Once heterogeneity is introduced, it can be maintained and amplified in a number of ways. Finally, there are a number of important ecological consequences of heterogeneity (Figure 6).
From a modeling perspective, the distinction of sources, amplifiers and consequences is important. Specifically, sources of heterogeneity are included in the design of the model. In other words, there are equations or parts of the model code that produce heterogeneity. For example, stochastic cell division asymmetry can be included in an agent-based model by randomly varying the daughter biomass from the perfect 50/50 split after division. Amplifiers are mechanisms that operate at the individual level and change the cell's properties. The resulting additional heterogeneity is not prescribed in the model design, but it emerges from running it (i.e., it is a model output). For example, heterogeneity in birth sizes may lead to heterogeneity in generation times without any added equation or code. Consequences are also not included in the design of a model, but they emerge as population- or ecosystem-level properties rather than individual-level properties.
Sources of Individuality
There are many mechanisms that can produce and maintain individual heterogeneity. Here we consider a mechanism that would lead a colony growing up from a single cell to become heterogeneous to be a “source of heterogeneity.” Of course these mechanisms can also operate and produce/modify heterogeneity in other scenarios.
Stochastic Interactions of Molecules
Intracellular “concentrations” of transcription factors and macromolecules (DNA, mRNA, proteins) are often low. For example, natural populations have on average less than one transcript per gene (Cottrell and Kirchman, 2016). That means the continuum assumption underlying deterministic chemical reaction kinetics is not met and corresponding regulatory or signaling networks can exhibit substantial stochasticity. This leads to intra-population heterogeneity and when coupled with positive feedbacks can lead to bi-stability and phenotypic differentiation (Veening et al., 2008). For example, the expression level of heterodisulfide reductase subunit A (central for respiration in sulfate reducers) in Desulfovibrio vulgaris cells varied by as much as 50-fold in a sample of 30 individual cells (Qi et al., 2014). Another example is stochasticity in the chemotaxis regulatory network. Low concentrations of signaling molecules, specifically phosphorylated CheY, lead to behavioral variability of individuals, and this can be reduced by increasing the concentration of this element in the network (Korobkova et al., 2004). However, stochastic gene expression may be more the exception than the rule as the expression of most genes in E. coli does not show any bursts (Silander et al., 2012). Metabolic pathways are usually considered to be unaffected by stochasticity because of the higher numbers of metabolic enzymes and metabolites in the cell, but the stochastic expression of catabolic enzymes has been found to lead to fluctuations in growth rate that can perturb the expression of other enzymes (Kiviet et al., 2014). Heterogeneity may also arise from changes to the DNA, including mutation, recombination and methylation (Avery, 2006; van der Woude, 2011).
Stochastic Cell Division Asymmetry
Another source of cell-to-cell variability is stochastic partitioning of cellular components during cell division. This may be due to low copy numbers of mRNAs, proteins, plasmids and genomes (Huh and Paulsson, 2011; Jahn et al., 2015), or imperfections in the cell division machinery leading to unequal daughter cell sizes and consequently asymmetry in all cellular components. Bacteria can control this heterogeneity by molecular mechanisms that are increasingly understood. For example, interactions of MinCDE proteins with themselves and the polar membranes set up a spatial gradient inside the cell that favors assembly of the FtsZ cell division ring in the middle of the cell (Kieser and Rubin, 2014). Missing components of this regulatory system have been implicated in the higher division asymmetry observed for Mycobacterium smegmatis compared to E. coli (Aldridge et al., 2012).
Deterministic Cell Division Asymmetry
In addition to stochastic processes, there are deterministic mechanisms that lead to cell division asymmetry. Replication in budding bacteria and yeast obviously produces two different individuals and population heterogeneity. In Saccharomyces cerevisiae, the mother cell is larger and accumulates damage, including bud scars, extrachromosomal DNA circles (ERCs) and carbonylated proteins, which are retained preferentially by the mother cell during division by binding to special cellular compartments (Unruh et al., 2013; Figure 7). Also, for cells that divide by apparently symmetric binary fission, one “daughter” inherits the old pole and one the new pole. The old pole may have accumulated more damage or other properties over its longer lifetime. For example, division in E. coli is associated with asymmetric segregation of damaged protein aggregates (Lindner et al., 2008). The aggregates diffuse by stochastic Brownian motion but they are too large to enter the nucleoid region and therefore get trapped at the poles (Coquel et al., 2013). The asymmetry goes beyond damage. For example, since the outer membrane (OM) is synthesized mostly along the cylindrical part of the cell, old poles have older OM and may include proteins that were previously expressed and now repressed (Ursell et al., 2012). Bacteria with flagella at one pole will also generate two different daughter cells at division (Christen et al., 2010).
Figure 7. Deterministic cell division asymmetry in the budding yeast Saccharomyces cerevisiae leads to heterogeneity in size and bud scars (stained bright). Scale bar is 5 μm. Reproduced from van Deventer et al. (2015) with permission.
Small-Scale Environmental Heterogeneity
Spatially structured microenvironments constitute another driver of heterogeneity. For example, microbes in colonies, biofilms or granular sludge flocs experience gradients in dissolved oxygen and nutrient concentrations (Wimpenny and Coombs, 1983; Matsumoto et al., 2010). The response of the microbes to these different conditions, including growth and acclimation, leads to a microbial population with heterogeneous properties. For example, the growth rate of Pseudomonas putida cells in biofilms was monitored using a reporter consisting of the growth rate-regulated rrnBP1 promoter and unstable GFP (Sternberg et al., 1999). Cells along the periphery of the biofilm were observed to grow rapidly, whereas those on the inside grew slower or not at all.
Differential Transport and Environmental Heterogeneity
Transport, whether passive (e.g., with water or air) or active (e.g., chemotaxis), can act differently on individuals within a homogenous population. Water flow velocities tend to be larger near the center of conduits (e.g., pores, pipes, rivers). Even for a uniform flow field, advection is generally associated with diffusion causing individuals from one location to be transported to different locations. On leaf surfaces, cells can be dispersed to different locations and grow into microcolonies, followed by detachment and colonization elsewhere (van der Wal et al., 2013). Also, transport by active mechanisms, like chemotaxis, entails stochastic variability and can lead to different paths of individuals (due to stochastic signaling, see above) (Korobkova et al., 2004). Often the environment exhibits substantial heterogeneity at this scale. For example, the nutrient concentrations in surface waters are highly heterogeneous, with microscale patches created by lysing cells, phytoplankton exudates or marine snow (Stocker, 2012; Taylor and Stocker, 2012; Zehr et al., 2017). Similarly, nutrient availability on plant leaf surfaces varies greatly at a micrometer scale and often correlates with local topography (Remus-Emsermann et al., 2011). Even environments that are designed to be homogenous, like strongly agitated small-scale fermentors, can be heterogeneous (Dunlop and Ye, 1990). When differential transport occurs in a heterogeneous environment, it can lead to intra-population heterogeneity as microbes respond to their individual conditions (i.e., by gene expression, nutrient uptake, and growth). By the same argument, the population at any given location may be comprised of individuals with vastly different life histories (Bucci et al., 2012). For example, depending on how they were formed (by a mechanism of staying-together or coming-together), aggregates of bacteria on leaf surfaces may consist of cells that are either clonal with similar life histories or represent a variety of previous leaf surface experiences (Tecon and Leveau, 2012).
Amplifiers of Individuality
The heterogeneity produced by the above sources can manifest itself in a number of properties and behaviors, which can then feed forward and produce heterogeneity in other properties and behaviors, effectively amplifying the overall heterogeneity.
Nutrient Uptake, Metabolism, and Growth
Stochastic gene expression, or any of the other primary sources of heterogeneity discussed above, may result in different levels of some functional enzyme and behavior, such as nutrient uptake. For example, the assimilation of nitrate and urea is very heterogeneous when a cultured population of nitrate-acclimated, marine dinoflagellate Prorocentrum minimum cells is exposed to a sudden input of urea (a preferred N source) (Figure 8; Matantseva et al., 2016). This may in turn affect nutrient metabolism and growth. For example, single-cell observations for Methylobacterium extorquens AM1 showed high variability in cell size at division, division time (2.5-fold range) and growth rate (Strovas et al., 2007). Consequently, even some bulk housekeeping functions, like metabolism and growth, which are generally considered to be relatively homogenous, can be very heterogeneous, even in cultured populations.
Figure 8. Heterogeneity in urea uptake by P. minimum at the single-cell level. (A) P. minimum cells in UV light. (B) 15N-urea uptake by P. minimum cells depicted as 12C15N−/12C14N− ratio. Scale bar is 5 μm. Reproduced from Matantseva et al. (2016) with permission.
Cell Cycle Asynchronicity
Asymmetric division can lead to an asynchronous population (i.e., where cells are in different phases in the cell cycle) because size can be a major checkpoint for various cell cycle phases. Since the cells perform different tasks at different phases in the cell cycle, this translates into a population with heterogeneous behavior. For example, in Saccharomyces cerevisiae, 800 genes are cell-cycle regulated (Spellman et al., 1998) and in Caulobacter crescentus, over 500 genes (Laub et al., 2000). In photosynthetic microorganisms, such as microalgae and cyanobacteria, gene expression is also tied to the light-dark cycle, often via a circadian clock (Ito et al., 2009). Unless the population grows at a generation time of 1 day, it will consist of cells with various phase differences between their cell and diel cycles. This effectively adds another dimension of variation and increases the number of phenotypes and population heterogeneity.
Age and Damage
Asymmetric segregation of damage during cell division produces younger and older cells and therefore an age distribution in the population. This affects the growth rates and other behaviors of cells. For instance, damaged protein aggregates are partitioned asymmetrically in E. coli and new-pole cells with less damage have a 4% higher specific growth rate (Lindner et al., 2008). In S. cerevisiae, older cells also grow slower and they synthesize more of the stress protectant trehalose (Levy et al., 2012).
Ecological Consequences of Individuality
In many cases the heterogeneity may simply average out and be of little consequence to the fitness of the population. However, there are a number of cases where heterogeneity has been shown to have important ecological consequences.
Division of Labor
Phenotypic differentiation forms the basis for a division of labor, where different cells carry out complementary tasks that benefit the population. For example, oxygenic photosynthesis and nitrogen fixation are incompatible processes. Specifically, the enzyme nitrogenase, encoded by genes nifH, nifD and nifK and responsible for reducing N2 to , breaks down in the presence of oxygen. To overcome this problem, the filamentous cyanobacterium Anabaena can differentiate into two types: photosynthesizing vegetative cells and nitrogen fixing heterocysts (Flores and Herrero, 2010). Another example includes evolved populations of E. coli where the labor of converting glucose to CO2 is divided over two cell types: one that converts glucose to acetate while the other converts acetate to CO2 (Harvey et al., 2014). There is also altruistic division of labor, which is the tasked sacrifice of some members of the group to benefit others. For example, Salmonella cells that invade the gut tissue get killed by the host immune system, but not before triggering a host response that kills other bacteria in the gut lumen but not the subpopulation of Salmonella cells that stayed behind and now have a competitive advantage in the gut lumen (Ackermann et al., 2008). The basis for this strategy lies in the stochastic expression of genes coding for a Type III Secretion System (T3SS) within the clonal Salmonella population: only a subset of cells within this population express a T3SS and it is this subset of cells that is capable of invading the gut tissue. Another example is the split into motile and immotile subpopulations of Pseudomonas aeruginosa that only together can generate mushroom-shaped biofilms (Ghanbari et al., 2016). Division of labor does not have to involve a direct effect of one phenotype on the other, but it may simply involve growing on different nutrients, like nitrate and urea, allowing the population to maximize uptake (Matantseva et al., 2016).
Bet Hedging
The future is uncertain and may bring unpredictable changes in stresses or any other environmental factors. If cells could react instantly to changes in their environment, a good strategy may be to rely on sensing and responding, but if the response is too slow, it is better to maintain a diversity of phenotypes (Kussell and Leibler, 2005), which is referred to as bet hedging. In the context of stress resistance, bet hedging is when a population contains some cells that are ill-adapted to the current environment but better adapted to potential future stresses. An important example are persister cells that are produced spontaneously, make up a small fraction of the population, are inhibited in growth (dormant) but can survive antibiotics (Balaban et al., 2004, 2013; Figure 9). Dormancy is widespread in microbes (Lennon and Jones, 2011), and when it involves a fraction of the cells and is not purely responsive to environmental conditions (i.e., at least partially spontaneous) it can also be considered a bet hedging strategy. For example, the cyanobacterium Anabaena forms akinetes that sink to the sediment bed and can serve as a seedbank for future blooms. This can be considered a bet hedging strategy to protect the population from being wiped out by some factor (washout, grazing), but it may also be a case of division of labor, as the migration to the sediment bed allows the cells to accumulate nutrients (Hellweger et al., 2008). Bet-hedging strategies do not have to involve an all-or-nothing differentiation, but can be gradual. For example, populations of the budding yeast have a gradual (and age-correlated) distribution of the stress-protectant trehalose (Levy et al., 2012). Finally, bet hedging is not restricted to stress resistance, but it may involve nutrient acquisition. For example, diversifying chemotactic behaviors in clonal populations could be an adaptation to foraging in variable environments (Frankel et al., 2014).
Figure 9. Survival of persister cells under antibiotic treatment. Growth of a hipA7 mutant, which produces a larger fraction of persisters, in microfluidic channels. Times are in hours:min. Bacteria are exposed to three phases, including growth medium (GM1), ampicillin (A) and then washing and again growth medium (GM2). Persister cells are marked with a red arrow. Reproduced from Balaban et al. (2004) with permission.
Aging
Aging is a strategy for eliminating damage from a population by concentrating it in a few cells that will eventually be discarded (i.e., die of old age). The alternative is to repair or eliminate the damage in some way. The evidence for aging to provide a significant ecological benefit in microbes is elusive, probably because the extent of damage segregation varies between species and environmental conditions. For example, for E. coli, one study showed asymmetric partitioning of damaged protein aggregates and decreased growth rates of older cells (Lindner et al., 2008), but in another study growth rates were not observed to decrease over many generations (Wang et al., 2010). Thus, there is an ongoing debate about the ecological benefits of aging in bacteria (Clegg et al., 2014; Koleva and Hellweger, 2015). Several recent studies suggest that aging does not increase fitness or does not occur under benign conditions but instead is a stress response at the population level (Coelho et al., 2013; Clegg et al., 2014; Iyer-Biswas et al., 2014; Vedel et al., 2016).
Sub-Optimality
In the absence of conditions that make heterogeneity advantageous (division of labor, bet hedging or aging), it is disadvantageous or sub-optimal. For any given set of (constant) conditions, there is only one optimal behavior that maximizes fitness (for one species). This has been explored in the context of nutrient assimilation and the effect on growth. Nutrient quotas of phytoplankton can be quite heterogeneous. This heterogeneity leads to a reduction in growth rate, compared to a hypothetical population with uniform quotas, due to the non-linearity of the underlying process (Bucci et al., 2012; Fredrick et al., 2013). When the cell's environment (and thus the heterogeneity) is controlled using microfluidic culturing technology, the growth rate increases compared to flask cultures (Dusny et al., 2012). Another case of sub-optimality stems from a mismatch between cell and environmental cycles. When a population of cells in different phases of their cell cycle grows in a cyclically varying environment (i.e., diel cycle in light or temperature), the cells have different alignments between these two cycles. It is reasonable to expect that some of those alignments may be more optimal than others, leading to sub-optimality.
Consideration of Individuality in Models
From the above review, it is clear that individual heterogeneity can have important effects on the ecology of microbes and the ecosystems harboring them. Any model that is to capture these effects, whether for advancing understanding or making predictions, has to be able to simulate the production and amplification of this heterogeneity. Therefore, when selecting a modeling strategy it is important to understand upfront the role of heterogeneity in the system, and how it is produced and amplified. Then, a modeling approach can be selected. For example, the ecology of infectious bacteria is in many ways controlled by bet hedging, which builds on individual heterogeneity. In some cases, such as persisters that can survive antimicrobial or other stresses, there are only two phenotypes. This type of heterogeneity has been modeled with the metabolic flux approach (Zhuang et al., 2012). Other cases involve a more gradual differentiation, like age-correlated stress resistance. Resolving this type of heterogeneity has been modeled using the agent-based approach (Hellweger et al., 2014).
Author Contributions
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
This review grew out of a recent @ASM conference entitled “The individual microbe: Single-cell analysis and agent-based modeling.” We would like to thank the American Society of Microbiology (ASM) for support and the conference participants for the stimulating discussions. Two reviewers provided constructive criticism. J-UK thanks the UK National Centre for Replacement, Refinement and Reduction of Animals in Research (NC3Rs) for funding the development of eGUT (grant NC/K000683/1). JL acknowledges funding by USDA-NIFA-AFRI through grant #2013-02075. WZ acknowledges funding by the National Science Foundation of China through grants No. 21621004 and No. 31370115. FH acknowledges funding by NSF through grants DEB-1240894 and ENG/ECCS-1404163.
References
Ackermann, M., Stecher, B., Freed, N. E., Songhet, P., Hardt, W.-D., and Doebeli, M. (2008). Self-destructive cooperation mediated by phenotypic noise. Nature 454, 987–990. doi: 10.1038/nature07067
Aldridge, B. B., Fernandez-Suarez, M., Heller, D., Ambravaneswaran, V., Irimia, D., Toner, M., et al. (2012). Asymmetry and aging of mycobacterial cells lead to variable growth and antibiotic susceptibility. Science 335, 100–104. doi: 10.1126/science.1216166
Anantharaman, K., Breier, J. A., Sheik, C. S., and Dick, G. J. (2013). Evidence for hydrogen oxidation and metabolic plasticity in widespread deep-sea sulfur-oxidizing bacteria. Proc. Natl. Acad. Sci. U.S.A. 110, 330–335. doi: 10.1073/pnas.1215340110
Atarashi, K., Tanoue, T., Oshima, K., Suda, W., Nagano, Y., Nishikawa, H., et al. (2013). Treg induction by a rationally selected mixture of Clostridia strains from the human microbiota. Nature 500, 232–236. doi: 10.1038/nature12331
Avery, S. V. (2006). Microbial cell individuality and the underlying sources of heterogeneity. Nat. Rev. Microbiol. 4, 577–587. doi: 10.1038/nrmicro1460
Balaban, N. Q., Gerdes, K., Lewis, K., and McKinney, J. D. (2013). A problem of persistence: still more questions than answers? Nat. Rev. Microbiol. 11, 587–591. doi: 10.1038/nrmicro3076
Balaban, N. Q., Merrin, J., Chait, R., Kowalik, L., and Leibler, S. (2004). Bacterial persistence as a phenotypic switch. Science 305, 1622–1625. doi: 10.1126/science.1099390
Becker, S. A., and Palsson, B. O. (2008). Context-specific metabolic networks are consistent with experiments. PLoS Comput. Biol. 4:e1000082. doi: 10.1371/journal.pcbi.1000082
Benyamini, T., Folger, O., Ruppin, E., and Shlomi, T. (2010). Flux balance analysis accounting for metabolite dilution. Genome Biol. 11:R43. doi: 10.1186/gb-2010-11-4-r43
Biggs, M. B., Medlock, G. L., Kolling, G. L., and Papin, J. A. (2015). Metabolic network modeling of microbial communities. Wiley Interdiscip. Rev. Syst. Biol. Med. 7, 317–334. doi: 10.1002/wsbm.1308
Biggs, M. B., and Papin, J. A. (2013). Novel multiscale modeling tool applied to Pseudomonas aeruginosa biofilm formation. PLoS ONE 8:e78011. doi: 10.1371/journal.pone.0078011
Bizukojc, M., Dietz, D., Sun, J., and Zeng, A.-P. (2010). Metabolic modelling of syntrophic-like growth of a 1,3-propanediol producer, Clostridium butyricum, and a methanogenic archeon, Methanosarcina mazei, under anaerobic conditions. Bioprocess Biosyst. Eng. 33, 507–523. doi: 10.1007/s00449-009-0359-0
Blumberg, A. F., and Di Toro, D. M. (1990). Effects of climate warming on dissolved oxygen concentrations in lake Erie. Trans. Am. Fish. Soc. 119, 210–223.
Bordbar, A., Lewis, N. E., Schellenberger, J., Palsson, B. Ø., and Jamshidi, N. (2010). Insight into human alveolar macrophage and M. tuberculosis interactions via metabolic reconstructions. Mol. Syst. Biol. 6:422. doi: 10.1038/msb.2010.68
Bucci, V., Nunez-Milland, D., Twining, B., and Hellweger, F. (2012). Microscale patchiness leads to large and important intraspecific internal nutrient heterogeneity in phytoplankton. Aquatic Ecol. 46, 101–118. doi: 10.1007/s10452-011-9384-6
Bucci, V., Tzen, B., Li, N., Simmons, M., Tanoue, T., Bogart, E., et al. (2016). MDSINE: microbial dynamical systems INference engine for microbiome time-series analyses. Genome Biol. 17:121. doi: 10.1186/s13059-016-0980-6
Buffie, C. G., Bucci, V., Stein, R. R., McKenney, P. T., Ling, L., Gobourne, A., et al. (2015). Precision microbiome reconstitution restores bile acid mediated resistance to Clostridium difficile. Nature 517, 205–208. doi: 10.1038/nature13828
Burnat, M., Herrero, A., and Flores, E. (2014). Compartmentalized cyanophycin metabolism in the diazotrophic filaments of a heterocyst-forming cyanobacterium. Proc. Natl. Acad. Sci. U.S.A. 111, 3823–3828. doi: 10.1073/pnas.1318564111
Castellanos, M., Wilson, D. B., and Shuler, M. L. (2004). A modular minimal cell model: purine and pyrimidine transport and metabolism. Proc. Natl. Acad. Sci. U.S.A. 101, 6681–6686. doi: 10.1073/pnas.0400962101
Cerqueda-García, D., and Falcón, L. I. (2016). Metabolic potential of microbial mats and microbialites: autotrophic capabilities described by an in silico stoichiometric approach from shared genomic resources. J. Bioinform. Comput. Biol. 14:1650020. doi: 10.1142/S0219720016500207
Chiu, H.-C., Levy, R., and Borenstein, E. (2014). Emergent biosynthetic capacity in simple microbial communities. PLoS Comput. Biol. 10:e1003695. doi: 10.1371/journal.pcbi.1003695
Christen, M., Kulasekara, H. D., Christen, B., Kulasekara, B. R., Hoffman, L. R., and Miller, S. I. (2010). Asymmetrical distribution of the second messenger c-di-GMP upon bacterial cell division. Science 328, 1295–1297. doi: 10.1126/science.1188658
Clark, J. R., Daines, S. J., Lenton, T. M., Watson, A. J., and Williams, H. T. (2011). Individual-based modelling of adaptation in marine microbial populations using genetically defined physiological parameters. Ecol. Modell. 222, 3823–3837. doi: 10.1016/j.ecolmodel.2011.10.001
Clegg, R. J., Dyson, R. J., and Kreft, J.-U. (2014). Repair rather than segregation of damage is the optimal unicellular aging strategy. BMC Biol. 12:52. doi: 10.1186/s12915-014-0052-x
Coelho, M., Dereli, A., Haese, A., Kühn, S., Malinovska, L., DeSantis, M. E., et al. (2013). Fission yeast does not age under favorable conditions, but does so after stress. Curr. Biol. 23. 1844–1852. doi: 10.1016/j.cub.2013.07.084
Cole, J. A., Kohler, L., Hedhli, J., and Luthey-Schulten, Z. (2015). Spatially-resolved metabolic cooperativity within dense bacterial colonies. BMC Syst. Biol. 9:15. doi: 10.1186/s12918-015-0155-1
Coquel, A.-S., Jacob, J.-P., Primet, M., Demarez, A., Dimiccoli, M., Julou, T., et al. (2013). Localization of protein aggregation in Escherichia coli is governed by diffusion and nucleoid macromolecular crowding effect. PLoS Comput. Biol. 9:e1003038. doi: 10.1371/journal.pcbi.1003038
Costello, E. K., Lauber, C. L., Hamady, M., Fierer, N., Gordon, J. I., and Knight, R. (2009). Bacterial community variation in human body habitats across space and time. Science 326, 1694–1697. doi: 10.1126/science.1177486
Cottrell, M. T., and Kirchman, D. L. (2016). Transcriptional control in marine copiotrophic and oligotrophic bacteria with streamlined genomes. Appl. Environ. Microbiol. 82, 6010–6018. doi: 10.1128/AEM.01299-16
Dick, G. J. (2017). Embracing the mantra of modellers and synthesizing omics, experiments and models. Environ. Microbiol. Rep. 9, 18–20. doi: 10.1111/1758-2229.12491
Dunlop, E. H., and Ye, S. J. (1990). Micromixing in fermentors: metabolic changes in Saccharomyces cerevisiae and their relationship to fluid turbulence. Biotechnol. Bioeng. 36, 854–864. doi: 10.1002/bit.260360816
Dusny, C., Fritzsch, F. S. O., Frick, O., and Schmid, A. (2012). Isolated microbial single cells and resulting micropopulations grow faster in controlled environments. Appl. Environ. Microbiol. 78, 7132–7136. doi: 10.1128/AEM.01624-12
El-Semman, I. E., Karlsson, F. H., Shoaie, S., Nookaew, I., Soliman, T. H., and Nielsen, J. (2014). Genome-scale metabolic reconstructions of Bifidobacterium adolescentis L2-32 and Faecalibacterium prausnitzii A2-165 and their interaction. BMC Syst. Biol. 8:41. doi: 10.1186/1752-0509-8-41
Emonet, T., and Cluzel, P. (2008). Relationship between cellular response and behavioral variability in bacterial chemotaxis. Proc. Natl. Acad. Sci. U.S.A. 105, 3304–3309. doi: 10.1073/pnas.0705463105
Feist, A. M., Herrgård, M. J., Thiele, I., Reed, J. L., and Palsson, B. Ø. (2008). Reconstruction of biochemical networks in microorganisms. Nat. Rev. Microbiol. 7, 129–143. doi: 10.1038/nrmicro1949
Flores, E., and Herrero, A. (2010). Compartmentalized function through cell differentiation in filamentous cyanobacteria. Nat. Rev. Microbiol. 8, 39–50. doi: 10.1038/nrmicro2242
Frankel, N. W., Pontius, W., Dufour, Y. S., Long, J., Hernandez-Nunez, L., and Emonet, T. (2014). Adaptability of non-genetic diversity in bacterial chemotaxis. Elife 3:e03526. doi: 10.7554/eLife.03526
Fredrick, N. D., Berges, J. A., Twining, B. S., Nuñez-Milland, D., and Hellweger, F. L. (2013). Use of agent-based modeling to explore the mechanisms of intracellular phosphorus heterogeneity in cultured phytoplankton. Appl. Environ. Microbiol. 79, 4359–4368. doi: 10.1128/AEM.00487-13
Freilich, S., Zarecki, R., Eilam, O., Segal, E. S., Henry, C. S., Kupiec, M., et al. (2011). Competitive and cooperative metabolic interactions in bacterial communities. Nat. Commun. 2:589. doi: 10.1038/ncomms1597
Fuhrman, J., Follows, M., and Forde, S. (2013). Applying “-omics” data in marine microbial oceanography. Eos Trans. Am. Geophys. Union 94, 241–241. doi: 10.1002/2013EO270006
Ghanbari, A., Dehghany, J., Schwebs, T., Müsken, M., Häussler, S., and Meyer-Hermann, M. (2016). Inoculation density and nutrient level determine the formation of mushroom-shaped structures in Pseudomonas aeruginosa biofilms. Sci. Rep. 6:32097. doi: 10.1038/srep32097
Gottstein, W., Olivier, B. G., Bruggeman, F. J., and Teusink, B. (2016). Constraint-based stoichiometric modelling from single organisms to microbial communities. J. R. Soc. Interface 13:20160627. doi: 10.1098/rsif.2016.0627
Guidi, L., Chaffron, S., Bittner, L., Eveillard, D., Larhlimi, A., Roux, S., et al. (2016). Plankton networks driving carbon export in the oligotrophic ocean. Nature 532, 465–470. doi: 10.1038/nature16942
Hanly, T. J., and Henson, M. A. (2011). Dynamic flux balance modeling of microbial co-cultures for efficient batch fermentation of glucose and xylose mixtures. Biotechnol. Bioeng. 108, 376–385. doi: 10.1002/bit.22954
Hanly, T. J., and Henson, M. A. (2013). Dynamic metabolic modeling of a microaerobic yeast co-culture: predicting and optimizing ethanol production from glucose/xylose mixtures. Biotechnol. Biofuels 6:44. doi: 10.1186/1754-6834-6-44
Harcombe, W. R., Riehl W. J., Dukovski, I., Granger, B. R., Betts, A., Lang A. H., et al. (2014). Metabolic resource allocation in individual microbes determines ecosystem interactions and spatial dynamics. Cell Rep. 7, 1104–1115. doi: 10.1016/j.celrep.2014.03.070
Harvey, E., Heys, J., and Gedeon, T. (2014). Quantifying the effects of the division of labor in metabolic pathways. J. Theor. Biol. 360, 222–242. doi: 10.1016/j.jtbi.2014.07.011
Heinken, A., and Thiele, I. (2015). Anoxic conditions promote species-specific mutualism between gut microbes in silico. Appl. Environ. Microbiol. 81, 4049–4061. doi: 10.1128/AEM.00101-15
Hellweger, F., Fredrick, N., and Berges, J. (2014). Age-correlated stress resistance improves fitness of yeast: support from agent-based simulations. BMC Syst. Biol. 8:18. doi: 10.1186/1752-0509-8-18
Hellweger, F. L. (2009). Carrying photosynthesis genes increases ecological fitness of cyanophage in silico. Environ. Microbiol. 11, 1386–1394. doi: 10.1111/j.1462-2920.2009.01866.x
Hellweger, F. L. (2010). Resonating circadian clocks enhance fitness in cyanobacteria in silico. Ecol. Modell. 221, 1620–1629. doi: 10.1016/j.ecolmodel.2010.03.015
Hellweger, F. L. (2013). Escherichia coli adapts to tetracycline resistance plasmid (pBR322) by mutating endogenous potassium transport: in silico hypothesis testing. FEMS Microbiol. Ecol. 83, 622–631. doi: 10.1111/1574-6941.12019
Hellweger, F. L. (2015). 100 Years since Streeter and Phelps: it is time to update the biology in our water quality models. Environ. Sci. Technol. 49, 6372–6373. doi: 10.1021/acs.est.5b02130
Hellweger, F. L., Clegg, R. J., Clark, J. R., Plugge, C. M., and Kreft, J.-U. (2016a). Advancing microbial sciences by individual-based modelling. Nat. Rev. Microbiol. 14, 461–471. doi: 10.1038/nrmicro.2016.62
Hellweger, F. L., Fredrick, N. D., McCarthy, M. J., Gardner, W. S., Wilhelm, S. W., and Paerl, H. W. (2016b). Dynamic, mechanistic, molecular-level modelling of cyanobacteria: Anabaena and nitrogen interaction. Environ. Microbiol. 18, 2721–2731. doi: 10.1111/1462-2920.13299
Hellweger, F. L., Kravchuk, E. S., Novotny, V., and Gladyshev, M. I. (2008). Agent-based modeling of the complex life cycle of a cyanobacterium (Anabaena) in a shallow reservoir. Limnol. Oceanogr. 53, 1227–1241. doi: 10.4319/lo.2008.53.4.1227
Huh, D., and Paulsson, J. (2011). Non-genetic heterogeneity from stochastic partitioning at cell division. Nat. Genet. 43, 95–100. doi: 10.1038/ng.729
Huthmacher, C., Hoppe, A., Bulik, S., and Holzhütter, H.-G. (2010). Antimalarial drug targets in Plasmodium falciparum predicted by stage-specific metabolic network analysis. BMC Syst. Biol. 4:120. doi: 10.1186/1752-0509-4-120
Ito, H., Mutsuda, M., Murayama, Y., Tomita, J., Hosokawa, N., Terauchi, K., et al. (2009). Cyanobacterial daily life with Kai-based circadian and diurnal genome-wide transcriptional control in Synechococcus elongatus. Proc. Natl. Acad. Sci. U.S.A. 106, 14168–14173. doi: 10.1073/pnas.0902587106
Iyer-Biswas, S., Wright, C. S., Henry, J. T., Lo, K., Burov, S., Lin, Y., et al. (2014). Scaling laws governing stochastic growth and division of single bacterial cells. Proc. Natl. Acad. Sci. U.S.A. 111, 15912–15917. doi: 10.1073/pnas.1403232111
Jahn, M., Günther, S., and Müller, S. (2015). Non-random distribution of macromolecules as driving forces for phenotypic variation. Curr. Opin. Microbiol. 25, 49–55. doi: 10.1016/j.mib.2015.04.005
Jöhnk, K. D., Huisman, J. E. F., Sharples, J., Sommeijer, B. E. N., Visser, P. M., and Stroom, J. M. (2008). Summer heatwaves promote blooms of harmful cyanobacteria. Glob. Chang. Biol. 14, 495–512. doi: 10.1111/j.1365-2486.2007.01510.x
Karr, J. R., Sanghvi J. C., Macklin D. N., Gutschow M. V., Jacobs J. M., Bolival, B., et al. (2012). A whole-cell computational model predicts phenotype from genotype. Cell 150, 389–401. doi: 10.1016/j.cell.2012.05.044
Kettle, H., Louis, P., Holtrop, G., Duncan, S. H., and Flint, H. J. (2015). Modelling the emergent dynamics and major metabolites of the human colonic microbiota. Environ. Microbiol. 17, 1615–1630. doi: 10.1111/1462-2920.12599
Khandelwal, R. A., Olivier, B. G., Röling, W. F. M., Teusink, B., and Bruggeman, F. J. (2013). Community flux balance analysis for microbial consortia at balanced growth. PLoS ONE 8:e64567. doi: 10.1371/journal.pone.0064567
Kieser, K. J., and Rubin, E. J. (2014). How sisters grow apart: mycobacterial growth and division. Nat. Rev. Microbiol. 12, 550–562. doi: 10.1038/nrmicro3299
Kindaichi, T., Nierychlo, M., Kragelund, C., Nielsen, J. L., and Nielsen, P. H. (2013). High and stable substrate specificities of microorganisms in enhanced biological phosphorus removal plants. Environ. Microbiol. 15, 1821–1831. doi: 10.1111/1462-2920.12074
Kiviet, D. J., Nghe, P., Walker, N., Boulineau, S., Sunderlikova, V., and Tans, S. J. (2014). Stochasticity of metabolism and growth at the single-cell level. Nature 514, 376–379. doi: 10.1038/nature13582
Klitgord, N., and Segrè, D. (2010). Environments that induce synthetic microbial ecosystems. PLoS Comput. Biol. 6:e1001002. doi: 10.1371/journal.pcbi.1001002
Klitgord, N., and Segrè, D. (2011). Ecosystems biology of microbial metabolism. Curr. Opin. Biotechnol. 22, 541–546. doi: 10.1016/j.copbio.2011.04.018
Kohler, R., Mooney, R. A., Mills, D. J., Landick, R., and Cramer, P. (2017). Architecture of a transcribing-translating expressome. Science 356, 194–197. doi: 10.1126/science.aal3059
Koleva, K. Z., and Hellweger, F. L. (2015). From protein damage to cell aging to population fitness in E. coli: insights from a multi-level agent-based model. Ecol. Model. 301, 62–71. doi: 10.1016/j.ecolmodel.2015.01.024
Korobkova, E., Emonet, T., Vilar, J. M. G., Shimizu, T. S., and Cluzel, P. (2004). From molecular noise to behavioural variability in a single bacterium. Nature 428, 574–578. doi: 10.1038/nature02404
Kreft, J.-U. (2004). Biofilms promote altruism. Microbiology 150, 2751–2760. doi: 10.1099/mic.0.26829-0
Kussell, E., and Leibler, S. (2005). Phenotypic diversity, population growth, and information in fluctuating environments. Science 309, 2075–2078. doi: 10.1126/science.1114383
Labhsetwar, P., Cole, J. A., Roberts, E., Price, N. D., and Luthey-Schulten, Z. A. (2013). Heterogeneity in protein expression induces metabolic variability in a modeled Escherichia coli population. Proc. Natl. Acad. Sci. U.S.A. 110, 14006–14011. doi: 10.1073/pnas.1222569110
Laub, M. T., McAdams, H. H., Feldblyum, T., Fraser, C. M., and Shapiro, L. (2000). Global analysis of the genetic network controlling a bacterial cell cycle. Science 290, 2144–2148. doi: 10.1126/science.290.5499.2144
Lennon, J. T., and Jones, S. E. (2011). Microbial seed banks: the ecological and evolutionary implications of dormancy. Nat. Rev. Microbiol. 9, 119–130. doi: 10.1038/nrmicro2504
Lenski, R. E., Ofria, C., Collier, T. C., and Adami, C. (1999). Genome complexity, robustness and genetic interactions in digital organisms. Nature 400, 661–664. doi: 10.1038/23245
Le Quéré, C., Harrison, S. P., Colin Prentice, I., Buitenhuis, E. T., Aumont, O., Bopp, L., et al. (2005). Ecosystem dynamics based on plankton functional types for global ocean biogeochemistry models. Glob. Chang. Biol. 11, 2016–2040. doi: 10.1111/j.1365-2486.2005.1004.x
Levy, R., and Borenstein, E. (2013). Metabolic modeling of species interaction in the human microbiome elucidates community-level assembly rules. Proc. Natl. Acad. Sci. U.S.A. 110, 12804–12809. doi: 10.1073/pnas.1300926110
Levy, S. F., Ziv, N., and Siegal, M. L. (2012). Bet hedging in yeast by heterogeneous, age-correlated expression of a stress protectant. PLoS Biol. 10:e1001325. doi: 10.1371/journal.pbio.1001325
Lindner, A. B., Madden, R., Demarez, A., Stewart, E. J., and Taddei, F. (2008). Asymmetric segregation of protein aggregates is associated with cellular aging and rejuvenation. Proc. Natl. Acad. Sci. U.S.A. 105, 3076–3081. doi: 10.1073/pnas.0708931105
Louca, S., and Doebeli, M. (2015). Calibration and analysis of genome-based models for microbial ecology. Elife 4:e08208. doi: 10.7554/eLife.08208
Louca, S., Hawley, A. K., Katsev, S., Torres-Beltran, M., Bhatia, M. P., Kheirandish, S., et al. (2016). Integrating biogeochemistry with multiomic sequence information in a model oxygen minimum zone. Proc. Natl. Acad. Sci. U.S.A. 113, E5925–E5933. doi: 10.1073/pnas.1602897113
Mahadevan, R., Edwards, J. S., and Doyle, F. J. (2002). Dynamic flux balance analysis of diauxic growth in Escherichia coli. Biophys. J. 83, 1331–1340. doi: 10.1016/S0006-3495(02)73903-9
Martin-Figueroa, E., Navarro, F., and Florencio, F. J. (2000). The GS-GOGAT pathway is not operative in the heterocysts. Cloning and expression of glsF gene from the cyanobacterium Anabaena sp. PCC 7120. FEBS Lett. 476, 282–286. doi: 10.1016/S0014-5793(00)01722-1
Matantseva, O., Skarlato, S., Vogts, A., Pozdnyakov, I., Liskow, I., Schubert, H., et al. (2016). Superposition of individual activities: urea-mediated suppression of nitrate uptake in the dinoflagellate Prorocentrum minimum revealed at the population and single-cell levels. Front. Microbiol. 7:1310. doi: 10.3389/fmicb.2016.01310
Matsumoto, S., Katoku, M., Saeki, G., Terada, A., Aoi, Y., Tsuneda, S., et al. (2010). Microbial community structure in autotrophic nitrifying granules characterized by experimental and simulation analyses. Environ. Microbiol. 12, 192–206. doi: 10.1111/j.1462-2920.2009.02060.x
Merino, M. P., Andrews, B. A., and Asenjo, J. A. (2015). Stoichiometric model and flux balance analysis for a mixed culture of Leptospirillum ferriphilum and Ferroplasma acidiphilum. Biotechnol. Prog. 31, 307–315. doi: 10.1002/btpr.2028
Mina, P., di Bernardo, M., Savery, N. J., and Tsaneva-Atanasova, K. (2013). Modelling emergence of oscillations in communicating bacteria: a structured approach from one to many cells. J. R. Soc. Interface 10:20120612. doi: 10.1098/rsif.2012.0612
Nagarajan, H., Embree, M., Rotaru, A.-E., Shrestha, P. M., Feist, A. M., Palsson, B. Ø., et al. (2013). Characterization and modelling of interspecies electron transfer mechanisms and microbial community dynamics of a syntrophic association. Nat. Commun. 4:2809. doi: 10.1038/ncomms3809
Paerl, H. W., Hall, N. S., and Calandrino, E. S. (2011). Controlling harmful cyanobacterial blooms in a world experiencing anthropogenic and climatic-induced change. Sci. Total Environ. 409, 1739–1745. doi: 10.1016/j.scitotenv.2011.02.001
Papin, J. A., Stelling, J., Price, N. D., Klamt, S., Schuster, S., and Palsson, B. O. (2004). Comparison of network-based pathway analysis methods. Trends Biotechnol. 22, 400–405. doi: 10.1016/j.tibtech.2004.06.010
Perez-Garcia, O., Lear, G., and Singhal, N. (2016). Metabolic network modeling of microbial interactions in natural and engineered environmental systems. Front. Microbiol. 7:673. doi: 10.3389/fmicb.2016.00673
Qi, Z., Pei, G., Chen, L., and Zhang, W. (2014). Single-cell analysis reveals gene-expression heterogeneity in syntrophic dual-culture of Desulfovibrio vulgaris with Methanosarcina barkeri. Sci. Rep. 4:7478. doi: 10.1038/srep07478
Reed, D. C., Algar, C. K., Huber, J. A., and Dick, G. J. (2014). Gene-centric approach to integrating environmental genomics and biogeochemical models. Proc. Natl. Acad. Sci. U.S.A. 111, 1879–1884. doi: 10.1073/pnas.1313713111
Reed, D. C., Breier, J. A., Jiang, H., Anantharaman, K., Klausmeier, C. A., Toner, B. M., et al. (2015). Predicting the response of the deep-ocean microbiome to geochemical perturbations by hydrothermal vents. ISME J. 9, 1857–1869. doi: 10.1038/ismej.2015.4
Remus-Emsermann, M., de Oliveira, S., Schreiber, L., and Leveau, J. (2011). Quantification of lateral heterogeneity in carbohydrate permeability of isolated plant leaf cuticles. Front. Microbiol. 2:197. doi: 10.3389/fmicb.2011.00197
Salcher, M. M., Posch, T., and Pernthaler, J. (2013). In situ substrate preferences of abundant bacterioplankton populations in a prealpine freshwater lake. ISME J. 7, 896–907. doi: 10.1038/ismej.2012.162
Salimi, F., Zhuang, K., and Mahadevan, R. (2010). Genome-scale metabolic modeling of a clostridial co-culture for consolidated bioprocessing. Biotechnol. J. 5, 726–738. doi: 10.1002/biot.201000159
Scheibe, T. D., Mahadevan, R., Fang, Y., Garg, S., Long, P. E., and Lovley, D. R. (2009). Coupling a genome-scale metabolic model with a reactive transport model to describe in situ uranium bioremediation. Microb. Biotechnol. 2, 274–286. doi: 10.1111/j.1751-7915.2009.00087.x
Schuster, S., Pfeiffer, T., and Fell, D. A. (2008). Is maximization of molar yield in metabolic networks favoured by evolution? J. Theor. Biol. 252, 497–504. doi: 10.1016/j.jtbi.2007.12.008
Segrè, D., Vitkup, D., and Church, G. M. (2002). Analysis of optimality in natural and perturbed metabolic networks. Proc. Natl. Acad. Sci. U.S.A. 99, 15112–15117. doi: 10.1073/pnas.232349399
Shoaie, S., Karlsson, F., Mardinoglu, A., Nookaew, I., Bordel, S., and Nielsen, J. (2013). Understanding the interactions between bacteria in the human gut through metabolic modeling. Sci. Rep. 3:2532. doi: 10.1038/srep02532
Silander, O. K., Nikolic, N., Zaslaver, A., Bren, A., Kikoin, I., Alon, U., et al. (2012). A genome-wide analysis of promoter-mediated phenotypic noise in Escherichia coli. PLoS Genet. 8:e1002443. doi: 10.1371/journal.pgen.1002443
Spellman, P. T., Sherlock, G., Zhang, M. Q., Iyer, V. R., Anders, K., Eisen, M. B., et al. (1998). Comprehensive identification of cell cycle–regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. Mol. Biol. Cell 9, 3273–3297. doi: 10.1091/mbc.9.12.3273
Stec, K. F., Caputi, L., Buttigieg, P. L., D'Alelio, D., Ibarbalz, F. M., Sullivan, M. B., et al. (2017). Modelling plankton ecosystems in the meta-omics era. Are we ready? Mar. Genomics 32, 1–17. doi: 10.1016/j.margen.2017.02.006
Stein, R. R., Bucci, V., Toussaint, N. C., Buffie, C. G., Rätsch, G., Pamer, E. G., et al. (2013). Ecological modeling from time-series inference: insight into dynamics and stability of intestinal microbiota. PLoS Comput. Biol. 9:e1003388. doi: 10.1371/journal.pcbi.1003388
Sternberg, C., Christensen, B. B., Johansen, T., Toftgaard Nielsen, A., Andersen, J. B., Givskov, M., et al. (1999). Distribution of bacterial growth activity in flow-chamber biofilms. Appl. Environ. Microbiol. 65, 4108–4117.
Stocker, R. (2012). Marine microbes see a sea of gradients. Science 338, 628–633. doi: 10.1126/science.1208929
Stolyar, S., Van Dien, S., Hillesland, K. L., Pinel, N., Lie, T. J., Leigh, J. A., et al. (2007). Metabolic modeling of a mutualistic microbial community. Mol. Syst. Biol. 3:92. doi: 10.1038/msb4100131
Strovas, T. J., Sauter, L. M., Guo, X., and Lidstrom, M. E. (2007). Cell-to-cell heterogeneity in growth rate and gene expression in Methylobacterium extorquens AM1. J. Bacteriol. 189, 7127–7133. doi: 10.1128/JB.00746-07
Taffs, R., Aston, J. E., Brileya, K., Jay, Z., Klatt, C. G., McGlynn, S., et al. (2009). In silico approaches to study mass and energy flows in microbial consortia: a syntrophic case study. BMC Syst. Biol. 3:114. doi: 10.1186/1752-0509-3-114
Tan, J., Zuniga, C., and Zengler, K. (2015). Unraveling interactions in microbial communities - from co-cultures to microbiomes. J. Microbiol. 53, 295–305. doi: 10.1007/s12275-015-5060-1
Taylor, J. R., and Stocker, R. (2012). Trade-offs of chemotactic foraging in turbulent water. Science 338, 675–679. doi: 10.1126/science.1219417
Tecon, R., and Leveau, J. H. J. (2012). The mechanics of bacterial cluster formation on plant leaf surfaces as revealed by bioreporter technology. Environ. Microbiol. 14, 1325–1332. doi: 10.1111/j.1462-2920.2012.02715.x
Tobalina, L., Bargiela, R., Pey, J., Herbst, F.-A., Lores, I., Rojo, D., et al. (2015). Context-specific metabolic network reconstruction of a naphthalene-degrading bacterial community guided by metaproteomic data. Bioinformatics 31, 1771–1779. doi: 10.1093/bioinformatics/btv036
Trivedi, P., Anderson, I. C., and Singh, B. K. (2013). Microbial modulators of soil carbon storage: integrating genomic and metabolic knowledge for global prediction. Trends Microbiol. 21, 641–651. doi: 10.1016/j.tim.2013.09.005
Tzamali, E., Poirazi, P., Tollis, G., and Reczko, M. (2009). Computational identification of bacterial communities. Int. J. Biol. Life Sci. 1, 185–191.
Tzamali, E., Poirazi, P., Tollis, I. G., and Reczko, M. (2011). A computational exploration of bacterial metabolic diversity identifying metabolic interactions and growth-efficient strain communities. BMC Syst. Biol. 5:167. doi: 10.1186/1752-0509-5-167
Unruh, J. R., Slaughter, B. D., and Li, R. (2013). Quality control: putting protein aggregates in a bind. Curr. Biol. 23, R74–R76. doi: 10.1016/j.cub.2012.12.005
Ursell, T. S., Trepagnier, E. H., Huang, K. C., and Theriot, J. A. (2012). Analysis of surface protein expression reveals the growth pattern of the Gram-negative outer membrane. PLoS Comput. Biol. 8:e1002680. doi: 10.1371/journal.pcbi.1002680
van der Wal, A., Tecon, R., Kreft, J.-U., Mooij, W. M., and Leveau, J. H. J. (2013). Explaining bacterial dispersion on leaf surfaces with an individual-based model (PHYLLOSIM). PLoS ONE 8:e75633. doi: 10.1371/journal.pone.0075633
van der Woude, M. W. (2011). Phase variation: how to create and coordinate population diversity. Curr. Opin. Microbiol. 14, 205–211. doi: 10.1016/j.mib.2011.01.002
van Deventer, S., Menendez-Benito, V., van Leeuwen, F., and Neefjes, J. (2015). N-terminal acetylation and replicative age affect proteasome localization and cell fitness during aging. J. Cell Sci. 128, 109–117. doi: 10.1242/jcs.157354
Varma, A., and Palsson, B. O. (1994). Stoichiometric flux balance models quantitatively predict growth and metabolic by-product secretion in wild-type Escherichia coli W3110. Appl. Environ. Microbiol. 60, 3724–3731.
Vedel, S., Nunns, H., Košmrlj, A., Semsey, S., and Trusina, A. (2016). Asymmetric damage segregation constitutes an emergent population-level stress response. Cell Systems 3, 187–198. doi: 10.1016/j.cels.2016.06.008
Veening, J. W., Smits, W. K., and Kuipers, O. P. (2008). Bistability, epigenetics, and bet-hedging in bacteria. Annu. Rev. Microbiol. 62, 193–210. doi: 10.1146/annurev.micro.62.081307.163002
Vila-Costa, M., Sharma, S., Moran, M. A., and Casamayor, E. O. (2013). Diel gene expression profiles of a phosphorus limited mountain lake using metatranscriptomics. Environ. Microbiol. 15, 1190–1203. doi: 10.1111/1462-2920.12033
Wang, P., Robert, L., Pelletier, J., Dang, W. L., Taddei, F., Wright, A., et al. (2010). Robust growth of Escherichia coli. Curr. Biol. 20, 1099–1103. doi: 10.1016/j.cub.2010.04.045
Wang, W.-L., Xu, S.-Y., Ren, Z.-G., Tao, L., Jiang, J.-W., and Zheng, S.-S. (2015). Application of metagenomics in the human gut microbiome. World J. Gastroenterol. 21, 803–814. doi: 10.3748/wjg.v21.i3.803
Wimpenny, J. W. T., and Coombs, J. P. (1983). Penetration of oxygen into bacterial colonies. Microbiology 129, 1239–1242. doi: 10.1099/00221287-129-4-1239
Wintermute, E. H., and Silver, P. A. (2010). Emergent cooperation in microbial metabolism. Molecular Systems Biology 6:407. doi: 10.1038/msb.2010.66
Ye, C., Zou, W., Xu, N., and Liu, L. (2014). Metabolic model reconstruction and analysis of an artificial microbial ecosystem for vitamin C production. J. Biotechnol. 182, 61–67. doi: 10.1016/j.jbiotec.2014.04.027
Zehr, J. P., Weitz, J. S., and Joint, I. (2017). How microbes survive in the open ocean. Science 357, 646–647. doi: 10.1126/science.aan5764
Zelezniak, A., Andrejev, S., Ponomarova, O., Mende, D. R., Bork, P., and Patil, K. R. (2015). Metabolic dependencies drive species co-occurrence in diverse microbial communities. Proc. Natl. Acad. Sci. U.S.A. 112, 6449–6454. doi: 10.1073/pnas.1421834112
Zengler, K., and Palsson, B. O. (2012). A road map for the development of community systems (CoSy) biology. Nat. Rev. Microbiol. 10, 366–372. doi: 10.1038/nrmicro2763
Zhuang, K., Izallalen, M., Mouser, P., Richter, H., Risso, C., Mahadevan, R., et al. (2011). Genome-scale dynamic modeling of the competition between Rhodoferax and Geobacter in anoxic subsurface environments. ISME J. 5, 305–316. doi: 10.1038/ismej.2010.117
Zhuang, K., Ma, E., Lovley, D. R., and Mahadevan, R. (2012). The design of long-term effective uranium bioremediation strategy using a community metabolic model. Biotechnol. Bioeng. 109, 2475–2483. doi: 10.1002/bit.24528
Zomorrodi, A. R., Islam, M. M., and Maranas, C. D. (2014). d-OptCom: dynamic multi-level and multi-objective metabolic modeling of microbial communities. ACS Synth. Biol. 3, 247–257. doi: 10.1021/sb4001307
Zomorrodi, A. R., and Maranas, C. D. (2012). OptCom: a multi-level optimization framework for the metabolic modeling and analysis of microbial communities. PLoS Comput. Biol. 8:e1002363. doi: 10.1371/journal.pcbi.1002363
Keywords: microbial ecology, gene-centric modeling, metabolic flux modeling, agent-based modeling, individuality, heterogeneity, single cell
Citation: Kreft JU, Plugge CM, Prats C, Leveau JHJ, Zhang W and Hellweger FL (2017) From Genes to Ecosystems in Microbiology: Modeling Approaches and the Importance of Individuality. Front. Microbiol. 8:2299. doi: 10.3389/fmicb.2017.02299
Received: 25 July 2017; Accepted: 07 November 2017;
Published: 27 November 2017.
Edited by:
Steve Lindemann, Purdue University, United StatesReviewed by:
Christopher Scott Henry, Argonne National Laboratory (DOE), United StatesMatthew B. Sullivan, Ohio State University, United States
Copyright © 2017 Kreft, Plugge, Prats, Leveau, Zhang and Hellweger. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ferdi L. Hellweger, ferdi@coe.neu.edu