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

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

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

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.

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 (O₂), ammonium (NH4+), nitrate (NO3-), and nitrite (NO2-). The model-predicted gene abundances were compared directly to observations from qPCR (gene copies L⁻¹, 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⁻¹ s⁻¹ vs. mRNA copies L⁻¹, 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.

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.

The approach is readily integrated into existing models based on concentration state variables (Reed et al., 2014).

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.

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

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.

Another gene-level ABM was developed by Hellweger (2013) to explore the mechanisms underlying adaptation of E. coli to tetracycline resistance.

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.

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.

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

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

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

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.

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

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.

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.

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

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 N₂ to NH4+, 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 CO₂ is divided over two cell types: one that converts glucose to acetate while the other converts acetate to CO₂ (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).

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

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

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.