To search for ecological patterning rules in microbial communities, we built a three-dimensional fitness model that ignored molecular details and instead focused on the fitness effects of interactions (‘The fitness model’ in ‘Materials and methods’). Specifically, two populations of cells, marked as red and green, were initially randomly distributed on a surface (Figure 1A). Cells grew horizontally until sufficiently confined, at which point they grew upward to accommodate the birth of new cells (Figure 1—figure supplement 1, Video 1). Thus, no active cell motility was present during community growth. The growth rate of a focal cell was determined by its basal fitness as a single cell and by its interactions with other cells in a defined interaction neighborhood. To reflect the negative fitness effects of intra- and inter-population competition for shared resources ([∼ ∼]), the fitness of the focal cell was decreased as the total population size in the interaction neighborhood increased. In addition, the focal and partner cells affected each other's fitness positively (↑), negatively (↓), or neutrally (∼). The magnitude of fitness effect is quantitative. Thus, to obtain qualitative ecological patterning rules, we focused on strong interactions in which ↑ and ↓ exert fitness effects large enough to be comparable to the recipient's basal fitness. In most simulations using the fitness model, the basal fitnesses of both partners were non-zero, and therefore ↑ represented strong facultative interactions: a participant could survive on its own at its basal fitness but fared much better in the presence of its partner.

Video 1

Download asset

We first analyzed the population composition of communities formed in the fitness model. Simulations (Figure 1—figure supplement 2 and Figure 1—source data 1) and analytical calculations (‘Requirements for steady-state ratios in the six types of communities’ in ‘Materials and methods’) show that for interactions benefiting at least one partner, including cooperation ([↑ ↑]), commensalism ([∼ ↑]), and exploitation ([↓ ↑]), different initial partner ratios can potentially converge over time (Figure 1B). Ratio convergence requires ‘balanced’ fitness. For example, for A[∼ ↑]B to achieve ratio convergence, the basal fitness of A must be higher than that of B and after gaining the fitness benefit from commensalism, B must be able to grow at least as fast as A. This way, a situation incompatible to ratio convergence, that is, one partner always fitter than the other, does not occur. Using engineered competitive, obligatory commensal, and obligatory cooperative yeast communities (see below), we confirmed that population ratios converged in the commensal and cooperative but not competitive communities (Figure 2—figure supplement 1). This convergence of population ratios reflects a balance between supply and consumption when growth is limited by supply (Shou et al., 2007): if, for instance, the supplier population suddenly increased in relative abundance, then each individual in the consumer population would receive more benefit and grow faster, which would return the ratio to its original value.

We next examined patterns in vertical cross-sections of communities simulated by the fitness model (Figure 1A, parameters in Figure 1—source data 1), because patterning along the x and y directions can depend on the initial spatial distribution of cells whereas patterning along the vertical z direction results from growth under the fitness influences of ecological interactions. In [∼ ∼], [∼ ↓], and [↓ ↓], red and green populations primarily formed columns that are spatially segregated from each other (Figure 1C–E). In [↑ ↓] and [∼ ↑], frequently one of the populations (green) either formed a column or became covered by the partner population (red) (Figure 1F–G). Only in cooperation ([↑ ↑]) conferring large fitness benefits to both partners, the two partner populations appeared to be extensively ‘intermixed’, manifested as the two different cell types successively piling on top of each other (Figure 1H).

To compare levels of intermixing in different communities, we defined an ‘intermixing index’ as the average number of cell type transitions spanning community height (‘Spatial analysis’ in ‘Materials and methods’). Since the intermixing index can be a function of community height, we compared intermixing indexes of simulated communities at equivalent heights. Statistically significant differences were observed between strong cooperation versus other types of interactions (Figure 1I). Thus, we predicted that partner intermixing would distinguish strong cooperation from other ecological interactions.

To test the prediction that strong cooperation is the only ecological interaction capable of driving partner intermixing, we engineered yeast communities engaged in competitive ([∼ ∼]), obligatory commensal ([∼ ↑]), and obligatory cooperative ([↑ ↑]) metabolic interactions. Competitive communities represented the baseline intra- and inter-population competition common to ‘all’ communities, while commensalism served as the most stringent control to be discriminated against.

All engineered yeast communities consisted of two non-mating S. cerevisiae strains, a G strain expressing GFP or YFP and an R strain expressing DsRed (see ‘Materials and methods’). In competitive communities, prototrophic R and G competed for shared nutrients in agarose and for space (Figure 2—figure supplement 1A). Depending on their genetic backgrounds, the two strains engaged in either equal-fitness or unequal-fitness competition. In obligatory commensal communities, $R_{\to A}^{\leftarrow L}$ took in lysine from the media and overproduced adenine to feed the adenine-requiring $G^{\leftarrow A}$ (Figure 2—figure supplement 1B). In obligatory cooperative communities (previously described as Cooperation that is Synthetic and Mutually Obligatory, or “CoSMO”), $R_{\to A}^{\leftarrow L}$ overproduced adenine and required lysine while $G_{\to L}^{\leftarrow A}$ overproduced lysine and required adenine (Shou et al., 2007) (Figure 2—figure supplement 1C). When mixed, the two cooperative strains could form a viable community, growing from low to high densities in synthetic minimal medium (SD) lacking adenine and lysine (Shou et al., 2007).

To predict and extrapolate experimental results of yeast communities, we developed a three-dimensional model based on the consumption, release, and diffusion of metabolites in the yeast communities (‘the diffusion model’ in ‘Materials and methods’). Specifically, the diffusion model assumed that metabolites diffused in the community and agarose (Figure 2—figure supplement 2) and that cell growth depended on the local concentration of its limiting metabolite according to Monod's equation (Monod, 1949). Most parameters in the diffusion model were measured experimentally (Figure 2—source data 1). We reasoned that if predictions from the diffusion model were consistent with experimental observations, we could use this model to simulate experiments that would be technically difficult to perform.

To examine vertical patterning in yeast communities, we first used top-view time-lapse images to infer patterns which were subsequently verified by cryosectioning. This is because confocal and two-photon microscopy cannot penetrate deep into yeast communities (Váchová et al., 2009). In competitive communities, whether in the diffusion model or experiments, time-lapse top-views suggested population segregation. For equal-fitness competitive communities (Figure 2A, left; Video 2), individual cells (i) grew into microcolonies (ii) which continued to grow and expand until they reached neighboring microcolonies (iii). After this stage, cells were unable to push other cells to the side, and further cell divisions mainly occurred in the vertical z direction (Figure 1—figure supplement 1). Consequently, columns of primarily a single cell type formed, and top-views of the community remained static (compare iv and v). Competitive communities composed of populations with different fitnesses developed similarly (Figure 2A, right; Video 3), except that the fitter population expanded in the top view during growth (compare the green population in iv and v).

Figure 2 with 4 supplements see all

Download asset Open asset

Figure 2—source data 1

Definitions and values of parameters used in the diffusion model.

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

Download elife-00230-fig2-data1-v1.xls

Video 2

Download asset

Video 3

Download asset

In obligatory commensal communities, time-lapse top-views suggested a ‘layered’ pattern with one population covering the other (Figure 2B; Video 4). We found that in both the diffusion model and experiments, $R_{\to A}^{\leftarrow L}$ initially grew rapidly by consuming lysine in the agarose (compare i, ii, and iii). In contrast, $G^{\leftarrow A}$ initially grew slowly, presumably limited by the low level of adenine released by $R_{\to A}^{\leftarrow L}$ before $R_{\to A}^{\leftarrow L}$ became abundant (compare i, ii, and iii). During the growth of $R_{\to A}^{\leftarrow L}$, some $G^{\leftarrow A}$ cells had been pushed to the top layer of the community (iii). As the $R_{\to A}^{\leftarrow L}$ population continued to expand and release adenine, $G^{\leftarrow A}$ started to grow rapidly (compare iii and iv). Eventually, $R_{\to A}^{\leftarrow L}$ stopped growing after lysine in the agarose had been depleted. Since the amount of adenine released during the lifetime of every $R_{\to A}^{\leftarrow L}$ cell could support the birth of multiple $G^{\leftarrow A}$ cells (Shou et al., 2007), $G^{\leftarrow A}$ population outnumbered and covered $R_{\to A}^{\leftarrow L}$ (v).

Video 4

Download asset

In contrast, time-lapse top-views of obligatory cooperative communities suggested population intermixing (Figure 2C; Video 5). After plating on agarose lacking adenine and lysine (i), cells in both populations divided once or twice by utilizing metabolites stored in their vacuoles (Shou et al., 2007). $R_{\to A}^{\leftarrow L}$ released adenine and entered the death phase while $G_{\to L}^{\leftarrow A}$ continued to grow by utilizing the released adenine (ii) (Shou et al., 2007). $G_{\to L}^{\leftarrow A}$ entered the death phase and released lysine after a significant delay due to their better starvation tolerance compared to $R_{\to A}^{\leftarrow L}$ (Shou et al., 2007). Released lysine supported the growth of surviving $R_{\to A}^{\leftarrow L}$ cells into microcolonies (iii). $R_{\to A}^{\leftarrow L}$ in turn released adenine which promoted growth of nearby $G_{\to L}^{\leftarrow A}$ cells and led to partial covering up of $R_{\to A}^{\leftarrow L}$ microcolonies by rapidly growing $G_{\to L}^{\leftarrow A}$ (iv and v, insets). Abundant $G_{\to L}^{\leftarrow A}$ cells provided enough lysine for rapid growth of local $R_{\to A}^{\leftarrow L}$ cells, which gave rise to patches of $R_{\to A}^{\leftarrow L}$ cells on top of the community (vi, inset). These community growth kinetics were consistent with previous measurements in liquid cultures (Shou et al., 2007).

Video 5

Download asset

To confirm that partners intermixed in cooperative but not non-cooperative communities, we obtained vertical cross-sections of competitive, commensal, and cooperative communities at their maximal community heights (‘Cryosectioning’ for experimental communities in ‘Materials and methods’). We found that in both the diffusion model and experiments, equal-fitness competitive communities formed segregated columns (Figure 2D, left). Unequal-fitness (Figure 2D, right) and obligatory commensal communities (Figure 2E) displayed a layered pattern in which the top portion of a community was dominated by one partner. In contrast to competitive and obligatory commensal communities, obligatory cooperative communities appeared to show population intermixing with patches of red and green cells emerging on top of each other (Figure 2F and Figure 2—figure supplement 4A and B). Cooperative intermixing appeared to be less in the diffusion model than in yeast communities. This is presumably because the diffusion model assumed that a confined cell would bud strictly upward whereas in yeast communities, cell divisions could stray to the side (Figure 3—figure supplement 1F).

To compare levels of intermixing in different communities, we quantified the intermixing index of each community (Figure 2G, symbols for experiments and lines for results from the diffusion model). As expected, equal-fitness competition resulted in small intermixing indexes at all community heights (Figure 2G, grey). In unequal-fitness competition (Figure 2G, black) and obligatory commensalism (Figure 2G, brown), the formation of a layered pattern caused a small increase in the intermixing index that subsequently leveled off. In contrast, the intermixing index in obligatory cooperative communities increased proportionally to community height in both experiments and the diffusion model (Figure 2G, magenta). This proportionality suggested the existence of a characteristic patch size (‘Spatial analysis’ in ‘Materials and methods’), denoted ${\lambda }_{\text{z}}^{*}$, of 10–20 µm. The characteristic patch size was independent of initial conditions (Figure 2—figure supplement 3A). Indeed, a calculation of the patch size based on experimentally-determined release, diffusion, and consumption of exchanged metabolites yielded comparable results (‘Calculating the characteristic patch size in cooperative yeast communities’ in ‘Materials and methods’).

What caused partner intermixing during cooperation? Using the diffusion model, we found that if cooperative benefits were not spatially localized because of instant distribution of benefits throughout the community or because of excessive supply levels, intermixing was diminished (Figure 2—figure supplement 3B). Based on this result and based on patterns observed in top-views and vertical cross-sections of cooperative communities, we propose that local deviations from the steady-state ratio coupled with localized large cooperative benefits cause partners to ‘take turns’ to grow, which leads to population intermixing (Figure 2H). In summary, in communities engaging in strong cooperation, but not in communities governed by other types of ecological interactions, the intermixing index increases proportionally as a function of community height. This proportionality is due to a fixed patch size determined by localized nutrient supply and consumption.

Under what conditions can we observe partner intermixing in cooperation? First, we experimentally tested partner intermixing in obligatory yeast cooperative communities initiated under different population ratios and densities. Next, we used the diffusion model to examine communities in which obligatory cooperative partners interacted with different dynamics. Finally, using the fitness model, we tested intermixing in facultative cooperation by varying the relative magnitude of cooperative benefits compared to the basal fitness of the two interacting populations. We utilized yeast communities when experimentally possible and otherwise took advantage of the diffusion and the fitness models.

We found that intermixing was insensitive to initial conditions in the yeast obligatory cooperative communities. The initial partner ratio did not significantly affect the level of intermixing (Figure 3—figure supplement 1D). At very high initial cell densities, we observed significant intermixing in communities directly above the inoculation area even in the absence of cooperation (Figure 3—figure supplement 1E, yellow shading). This is because high cell densities put different cell types in close proximity to one another, and cell divisions that were not perfectly straight upward (Figure 3—figure supplement 1F) caused intermixing. However, we reasoned that communities beyond the inoculation area might reveal patterns indicative of the underlying interactions, because these regions are formed by cell growth under the fitness influences of interactions. To test this possibility, we spotted cell mixtures at high densities on agarose and allowed the community to expand to new territories beyond the inoculation area. Even though communities directly above the inoculum always appeared highly intermixed (Figure 3A, ‘Center’ and Figure 2—figure supplement 4C), in edge sections, significant intermixing was only observed in cooperative communities (Figure 3A, ‘Edge’ and Figure 2—figure supplement 4D).

Figure 3 with 1 supplement see all

Download asset Open asset

Figure 3—source data 1

Parameter values used in the fitness model in Figure 3B.

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

Download elife-00230-fig3-data1-v1.xls

Intermixing is insensitive to interaction dynamics, so long as large cooperative benefits remain sufficiently localized for both partners. As described above, in the diffusion model, excessive supply amounts or instant distribution of benefits throughout the community diminished intermixing (Figure 2—figure supplement 3B). The former can occur if the local availability of cooperative benefits is not growth-limiting. However, cooperative benefits are unlikely to be available in large excess because of the potential fitness cost of generating benefits and because of competition for these benefits in the partner population. As far as diffusion is considered, the diffusion model showed significant population intermixing in obligatory cooperative communities with diffusion constants varying more than 10-fold (Figure 3—figure supplement 1A). Additional calculations showed that intermixing was largely insensitive to diffusion kinetics, because the characteristic patch size was related to the diffusion constant by 1/5 power (‘Calculating the characteristic patch size in cooperative yeast communities’ in ‘Materials and methods’). The diffusion model also showed that intermixing was present with or without a delay in $G_{\to L}^{\leftarrow A}$ supplying lysine benefits to the partner (Figure 3—figure supplement 1B). Furthermore, even though asymmetry between properties of partners (Figure 2C) would seem to facilitate intermixing (Figure 2H), asymmetry in partner properties was not required to generate intermixed patterns: partners that grew, died, and released and consumed metabolites at identical rates intermixed (Figure 3—figure supplement 1C).

Using the fitness model, we found that in facultative cooperation, small fitness benefits generated less intermixing than large fitness benefits (Figure 3B). This result is intuitive: facultative cooperation should resemble obligatory cooperation if fitness benefits are large for both partners. When only one partner receives a large fitness benefit, facultative cooperation should resemble obligatory commensalism. Finally, small fitness benefits for both partners will make facultative cooperation resemble competition. In facultative cooperation with smaller fitness benefits, intermixing would be apparent if communities could grow to greater heights. Further experiments are required to test these predictions. Together, these results suggest that intermixing relies on spatial localization of cooperative benefits that are sufficiently large to both partners, and is otherwise insensitive to initial conditions or the detailed kinetics of interactions.

To test whether cooperative patterning applies to other biological systems, we examined redox-coupling in a two-species methane-producing biofilm consisting of the bacterium Desulfovibrio vulgaris and the archaeon Methanococcus maripaludis. In the absence of sulfate and hydrogen, the two species engage in obligatory mutualism: D. vulgaris ferments lactate and promotes the growth of M. maripaludis by supplying the electron donor H₂, while M. maripaludis promotes the growth of D. vulgaris by consuming the H₂ gas which can be inhibitory to D. vulgaris under these conditions (Figure 4A). Similar types of syntrophic interactions leading to methane production typically occur in microbial consortia that digest organic compounds in freshwater sediments, sewage treatment plants, and the guts of ruminants (Conrad et al., 1985; Schink, 1997).

Figure 4

Download asset Open asset

Vertical cross-sections of independent D. vulgaris–M. maripaludis biofilms indeed exhibited increasing intermixing as a function of community height (Figure 4B). Thus, naturally mutualistic microbes cooperating through a coupling mechanism different from metabolite exchange also exhibited a significant level of intermixing. Other known cooperative communities, including those degrading herbicide pollutants (Breugelmans et al., 2008) and those colonizing teeth (Palmer et al., 2001) also seemed to display intermixed patterns, although we do not know whether the intermixing indexes of these communities increased linearly with height.

In communities with more than two species, indirect interactions can obscure direct interactions (Wootton, 2002). For example, if A promotes B which inhibits C, it will appear that A inhibits C. This is akin to indirect interactions between gene products in a cell. To test whether intermixing between cooperators was affected by other members of a community, we used the fitness model to simulate communities composed of five interacting species (Figure 5 and Figure 5—figure supplement 1). We randomly assigned one of the six possible ecological interactions between each pair of populations, and consequently, each network had 10 pair-wise interactions. The fitness effects of ↑ and ↓ were sufficiently large to be comparable to the recipient's basal fitness. A total of 240 pair-wise interactions in 24 independent communities were examined in the fitness model. In most cases (26 out of 31), cooperative pairs intermixed (Figure 5A,D). Occasionally, commensal pairs (Figure 5B, ② [∼ ↑] ③) showed substantial intermixing (3 out of 58 commensal pairs, Figure 5D), or cooperative pairs (e.g., Figure 5C, ① [↑ ↑] ③) showed little intermixing (5 out of 31, Figure 5D), presumably due to interactions through other community members (Figure 5B, ③ indirectly promoted ② through ④; Figure 5C, ③ promoted ⑤ which inhibited ①). These deviations are consistent with the notion that the intermixing index reflects the ‘overall’ interaction between two partners, integrating any additional fitness effects of indirect interactions through other community members. Thus, strongly cooperating partners intermix while deviations from this expectation reflect the presence of strong indirect interactions.

Figure 5 with 1 supplement see all

Download asset Open asset

Patterning is driven by cell–cell and cell–environment interactions. Here, we examined how the net fitness effects of cell–cell interactions could influence patterning. Using defined mathematical and biological systems under controlled spatial environments, we have established two expectations for ecological patterning between two interacting partners.

The first expectation is that interactions conferring large fitness benefits to at least one partner ([↑ ↑], [∼ ↑], and [↓ ↑]) can potentially result in ratio convergence (Figure 1—figure supplement 2 and Figure 2—figure supplement 1). Ratio convergence confers stability to community composition. The composition of a community is defined to be locally (or globally) stable if after small (or any-size) deviations from it, either due to initial conditions different from it or perturbations, the community eventually returns to this composition. An analytical calculation based on the fitness model showed that in [↑ ↑] and [∼ ↑], global stability can be achieved whereas in [↓ ↑], local stability can be achieved (Figure 1—figure supplement 2; ‘Requirements for steady-state ratios in the six types of communities’ in ‘Materials and methods’). Experimentally, obligatory cooperation and obligatory commensalism led disparate initial ratios to converge (Figure 2—figure supplement 1). In contrast, in competition between two populations of equal fitness, population ratio is determined by the initial value until a perturbation resets the value which will remain unchanged until the next perturbation. In competition between two populations of unequal fitness, the fitter population should increase in frequency and therefore, population ratio is unstable. In amensalism and mutual antagonism, population ratios at best have an unstable fixed point (Figure 1—figure supplement 2B). This means that even if the population ratio may seem ‘fixed’, slight deviations will send it farther and farther away from the original fixed value. The reason for this instability, using mutual antagonism as an example, is because if A inhibits B and B inhibits A (A[↓ ↓]B), then an increase in population A will facilitate inhibition of B and therefore make B less able to inhibit A, further increasing the relative abundance of A.

The second expectation is that strong cooperation results in partner intermixing in the direction vertical to the surface of initial colonization. Cooperative intermixing requires spatially localized benefits that are sufficiently large to both partners (Figure 2—figure supplement 3 and Figure 3B), and otherwise appears to be robust against variations in initial conditions (Figure 2—figure supplement 3A, Figure 3A, and Figure 3—figure supplement 1D) or reaction dynamics (Figure 3—figure supplement 1A–C). We have observed cooperative intermixing in the fitness model simulating strong facultative cooperation (Figure 1) and in obligatory experimental systems including engineered yeast communities (Figure 2) and syntrophic biofilms (Figure 4). The intermixing indexes in these communities increased linearly as a function of community height, because these communities exhibited fixed-size patches associated with the spatial localization of benefits (Figure 2—figure supplement 3; ‘Calculating the characteristic patch size in cooperative yeast communities’ in ‘Materials and methods’). This is in stark contrast with other ecological interactions which lead to segregated or layered patterns (Figures 1 and 2) in which the intermixing index remained low or increased transiently before leveling off. In simulated multi-species communities, strongly cooperative pairs intermixed in most cases; cases where cooperative pairs did not intermix or non-cooperative pairs intermixed were likely caused by strong indirect interactions through other partners. In a theoretical study using a one-dimensional stepping-stone model, strongly cooperating partners also appeared much more mixed than competing partners (Korolev and Nelson, 2011).

Our work defines the expected pattern created by pair-wise interactions if the fitness effects of interactions are the main driving force of patterning. The generality of these ecological patterning conclusions awaits further experimental validations. Many simple communities grown in laboratory environments conform to these findings. For instance, competing bacterial species form a columnar or layered pattern (Palmer et al., 2001; Kreth et al., 2005; An et al., 2006; Hallatschek et al., 2007; Bernstein et al., 2012). Burkholderia sp. LB400 and Pseudomonas sp. B13(FR1) show spatial segregation as competitors when grown on citrate, but when grown on chlorobiphenyl they engage in metabolic commensalism and exhibit a layered pattern (Nielsen et al., 2000). Similarly, when grown on a non-selective carbon source, Comamonas testosteroni WDL7 outcompeted and covered Variovorax sp. WDL1, but when fed with an aromatic compound as the sole carbon source, competition switched to cooperation and the two populations intermixed (Breugelmans et al., 2008). Cooperative oral bacteria species intermixed when grown on saliva (Palmer et al., 2001; Periasamy and Kolenbrander, 2009). If cooperative intermixing is suspected, then examining whether the intermixing index increases linearly as a function of community height will be informative to exclude transient increases.

We have described the expected patterning between two cell populations if the fitness effects of interactions are the major driving force. These expectations are abstract in the sense that in reality, no interactions can exist in the absence of molecular mechanisms or an abiotic environment. Assigning expected patterns to different types of ecological interactions will allow us to identify deviations from expectations. Determining the causes of such deviations will help us better understand the biology of a community.

Fitness effects of ecological interactions rely on molecular mechanisms to manifest themselves. Different types of molecular mechanisms can conceivably alter patterning. For instance, in the fitness and the diffusion models and in the S. cerevisiae communities, cells divide upward when sufficiently confined horizontally. This type of cell rearrangement has also been observed in bacterial colonies (Kreft et al., 1998) and was adopted to model biofilm growth (Xavier and Foster, 2007). However, it is conceivable that if at least one population actively moves (Houry et al., 2012) or grows hyphae to penetrate the entire community, two populations might appear intermixed even if they do not cooperate. In biofilms of Pseudomonas aeruginosa, two populations differing only in fluorescent colors (Klausen et al., 2003) showed modest intermixing even though the expectation for equal-fitness competition is a columnar pattern with an intermixing index close to zero (Figure 2D, left). This modest intermixing was caused by P. aeruginosa differentiating into non-motile ‘stalk’ cells that anchored to the surface and motile cells that climbed up to form the mushroom-like caps (Klausen et al., 2003).

Environmental influences can also alter ecological patterning. For instance, if two cooperating populations have very different preferences for oxygen, then the two populations are likely not to intermix and instead form layers in which the aerobic population is exposed to oxygen while the anaerobic population is protected from oxygen.

Variability in patterning has been observed within and between communities (Christensen et al., 2002; Wilmes et al., 2008; Dekas et al., 2009), even if they were grown in laboratory-controlled environments (Christensen et al., 2002). Stochastic events such as environmental variability, mutations, or fluctuations in initial conditions can all lead to variable patterns. For instance, if the straight-columnar pattern expected for equal-fitness competition is observed for the majority of community but a layered pattern is observed in occasional locations, then fitter mutants may have arisen in these locations, initiating unequal-fitness competition (Hallatschek et al., 2007; Korolev et al., 2012). For mutually antagonistic interactions, the fitness model showed that population ratios can quickly diverge. Thus, which population eventually dominates depends on the initial population ratio (Figure 1—figure supplement 3). In this case, stochastic variations in initial conditions can result in dramatically different patterns, giving rise to a phenomenon equivalent to ‘survival of the first’. Indeed, different patterns have been observed for communities formed by two antagonistic bacteria species (Kuramitsu et al., 2007).

In summary, our work is conceptually analogous to that of the competitive exclusion principle (Gause, 1934). The competitive exclusion principle, also known as Gause's law, states that two species competing for the exact same resources cannot stably coexist. Analogous to how the competitive exclusion principle has created a framework to examine forces that cause species coexistence, our work on ecological patterning will hopefully lay the ground for examining mechanisms that shape patterning in microbial communities. We encourage comments, especially those pertaining to the generality of our conclusions, to be posted to the eLife website.

In competitive communities, equal-fitness G and R strains were respectively WS931 (MATa ste3::kanMX4 trp1-289::pRS404(TRP)-ADHp-venus-YFP) and WS937 (MATa ste3::kanMX4 trp1-289::pRS404(TRP)-ADHp-DsRed.T4), both from the S288C background. For unequal-fitness competitive communities, G was replaced by WS1246 (MATa ho::loxP AMN1-BY Supercontig17(27163-27164)::TDH3p-yEGFP-loxP-kanMX-loxP) from the RM11 background with a 20% fitness advantage over the S288C background. In commensal communities, $R_{\to A}^{\leftarrow L}$ and $G^{\leftarrow A}$ strains were respectively WS950 (MATa ste3::kanMX4 lys2Δ0 ade4::ADE4(PUR6) trp1-289::pRS404(TRP)-ADHp-DsRed.T4) and WS932 (MATa ste3::kanMX4 ade8Δ0 trp1-289::pRS404(TRP)-ADHp-venus-YFP), both from the S288C background. In cooperative communities, $R_{\to A}^{\leftarrow L}$ and $G_{\to L}^{\leftarrow A}$strains were respectively WS950 and WS954 (MATa ste3::kanMX4 ade8Δ0 lys21::LYS21(fbr) trp1-289::pRS404(TRP)-ADHp-venus-YFP) from the S288C background.

For yeast communities, agarose columns were prepared by pouring SD minimal medium (Sherman, 2002) with 2% low melting temperature agarose in 1.3-ml deep-well plates (U96 DeepWell from Nunc, Penfield, NY), and after agarose had solidified, adding drops of melted SD agarose to make the surface flat. For commensal yeast communities, SD was supplemented with lysine (80 μM final concentration). Batch cultures of yeast strains were grown to exponential phase at 30°C in SD with supplements when necessary. Cells were washed free of supplements if any, mixed, and filtered on top of MF membrane (HAWP04700 from Millipore, Billerica, MA). Disks were cut from the membrane using a 6-mm-diameter Harris Uni-Core puncher and transferred to the top of agarose columns, unless otherwise stated.

For the coculture of D. vulgaris and M. maripaludis, a CDC reactor (BioSurface Technologies Corp., Bozeman, MT) was used for anaerobic biofilm growth. Biofilm coupon holders were modified to hold glass microscope slides (Fisher Scientific #12-544-1, Fisher Scientific, Pittsburgh, PA) cut to 7.6 × 1.8 cm. Cocultures were grown in CCM (Walker et al., 2009), a modified basal salt medium without choline chloride. Headspace was sparged with anoxic 80% N₂:20% CO₂, and the reactor was maintained at 30°C with stirring (150 rpm). The reactors were inoculated with planktonic coculture, and after cell attachment to the glass slides (48 hr), biofilms were allowed to develop and grow over time in the presence of planktonic cells, or initially drained of planktonic cells. Patterns in the two experimental regimes were similar.

Flow cytometry was performed in a modified 4-laser FACS Calibur machine (DxP; CyTek Development, Fremont, CA). Cells in communities were diluted to ∼10⁶ cells/ml in H₂O. Each sample (90 µl) was supplemented with 10 µl of fluorescent bead stock (Thermo Scientific Fluoro-Max Cat# R0300 at ∼8 × 10⁶/ml, counted using a Z2 Coulter counter and a hemacytometer) as a reference to determine total cell density and 3 µl of 1 µM ToPro3 (Invitrogen, Grand Island, NY) to mark dead cells. The laser and filter configurations for different fluorophores were 50 mW 488 nm laser with 530/30 filter for YFP, 75 mW 561 nm laser with 575/26 filter for DsRed, and 25 mW 639 nm laser with 660/20 filter for ToPro3. Using an automatic micro-sampling system (DxP; CyTek Development, Fremont, CA) samples in 96-well plates were processed for 60 s at a flow rate of 0.5–1 µl/s, recording 10⁴–10⁵ events. FlowJo software (Tree Star, Ashland, OR) was used to measure the ratio of different populations of fluorescent cell against the bead standard in order to calculate cell densities.

A Nikon TE2000 inverted microscope equipped with a Prior stage controller, a Sutter Lambda XL fluorescent lamp, and a Photometric SnapHQ CCD camera was controlled by custom LabVIEW software to auto-focus and acquire images. All images were taken at 10× magnification using a Nikon long working distance CFI Plan Fl objective (NA 0.30, WD 16). For imaging YFP- and DsRed-tagged strains, ET500/20×-ET535/30m-T515LP and ET545/30×-ET620/60m-T570LP filter cubes were used, respectively. Timelapse imaging took place in a 30°C chamber (In Vivo Scientific microscope incubator).

To obtain vertical cross-sections of communities through cryosectioning, we froze communities in liquid nitrogen for 15 s and fixed them in methanol at −20°C. After 20 min, the communities were transferred to a pre-cooled empty plate at −20°C to allow methanol to evaporate, which typically took 4 hr. The communities were embedded in optimal cutting temperature (OCT) compound for 10 min at room temperature and subsequently frozen over dry ice and kept at −20°C for sectioning. We also froze down communities without fixing by directly embedding them in OCT (Piccirillo et al., 2010) and immediately freezing them on dry ice. Results obtained from the two procedures yielded similar conclusions.

For sectioning, communities embedded in OCT blocks were mounted on a cryotome. The blade was adjusted to ensure cross-sections of the community were vertical. For each community, typically thirty to fifty 14-μm (3-cell-thick) sections were cut and transferred to glass slides. Cross-sections were imaged using the fluorescence microscope as described above. Only images of cross-sections that were minimally perturbed by the fixing and sectioning processes were included in the analysis. More details of the cryosectioning method can be found in (Momeni and Shou, 2012).

Biofilms of D. vulgaris and M. maripaludis were fixed in 4% paraformaldehyde for 4 hr, and then embedded in polyacrylamide (Daims et al., 2006). The embedded biofilms were dehydrated and hybridized in 1 ml buffer solution (0.9 M NaCl, 20 mM Tris–HCl (pH 8), 0.01% SDS, and 35% deionized formamide) with 3 ng each of probes EUB338 (GCT GCC TCC CGT AGG AGT) 5′ and 3′-labeled with Cy3 and ARCH915 (GTG CTC CCC CGC CAA TTC CT) 5′ and 3′-labeled with Cy5 (Stoecker et al., 2010) for 5 or 8 hr at 46°C in a humid chamber. Next, samples were washed in 50 ml prewarmed washing buffer (70 mM NaCl, 20 mM Tris–HCl pH 8, and 5 mM EDTA) at 47°C for 20 min, then dipped in ice cold ddH₂O and quickly dried with compressed air (Amann et al., 1995). Finally, each sample was mounted with Citifluor AF1 antifadent (Citifluor Ltd., Leicester, United Kingdom) and viewed using a Leica TCS SP5 II inverted confocal laser scanning microscope with 488, 561 and 633 nm lasers and appropriate filter sets for Cy3 and Cy5. Confocal voxel size was typically 0.24 × 0.24 × 0.49 μm.

For both simulated and experimental sections, our unit of analysis was the size of one CCD camera field of view under a 10× objective which has a width of 0.7 mm. In any unit of analysis, community height is the value such that 90% of height values are below it. This choice is made to exclude artifacts such as height spikes in simulations. Images of experimental community sections were rotated in ImageJ such that the x–z axes of the frame matched those of the community. We further digitized these images into f(x,z) with assigned values of +1, −1, and 0 for pixels identified as population 1, population 2, and no-signal, respectively. No-signal pixels were defined as having a fluorescence intensity per unit exposure time less than 10–20% above the background in fluorescence channels. For the remaining pixels, the intensity values of 90th percentiles were found in each fluorescence channel. These values were used to normalize the corresponding green- or red-fluorescence intensity for each pixel whose identity was then assigned to be the color with the higher normalized value. In simulated communities with more than two populations, analysis was performed on two focal colors at a time, and in each analysis, pixels identified as other colors were treated as no-signal.

To compare levels of intermixing in different communities, we estimated the number of cell type transitions spanning community height. We define the intermixing index as the average number of color changes along community height: let h(x_i) and c(x_i) respectively be the local height and the number of color changes along z at the lateral position x_i. The intermixing index IM can be calculated as

Intermixing is small for segregated patterns with few color changes along height, and increases when patches of different cell types successively appear on top of each other. Note that this choice of intermixing index yields small values when one population is very rare compared to the other. We weighted c(x_i) by the local community height h(x_i), thus giving more emphasis to taller regions. This is because taller regions have gone through more growth regulated by the fitness effects of interactions.

To estimate the vertical patch size ${\lambda }_{\text{z}}^{*}$, the ratio of height h(x_i) to number of cell-type layers, 1 + c(x_i), was averaged across the community cross-section

For engineered yeast communities, typically 10–20 frames from different locations of each community were included in the analysis to ensure an unbiased representation of community patterns. For biofilms of D. vulgaris and M. maripaludis, vertical cross-sections ∼2.4 μm apart were sampled from confocal z-stack images of biofilms.

The characteristic patch size can also be calculated as the effective length of interaction between two partners. We calculate how far inside the community released lysine can diffuse before being consumed; similar discussion applies to adenine. Assume one cell has released β_L fmole of lysine that diffuses at most a distance l before being consumed. We define the sphere of radius l as the diffusion domain of the release event. The number of consuming cells within the diffusion domain can be estimated as N_u = (2l/c)³, where c is the diameter of a cell. The average time nutrient can diffuse in the diffusion domain before being consumed is t_c = l²/2D₁, where D₁ = 360 μm²/s is the diffusion constant within the community. Assuming that each consuming cell takes up lysine with a rate v ≤ α_L/T, with α_L being the amount of lysine required for a cell division and T being the minimum doubling time, we have

thus,

Within radius l, the average number of release events within time t_c is ${\left(2l/c\right)}^3{t}_c{d}_G\approx 0.2$. Thus, the probability of a second release event occurring in the same diffusion domain before nutrient from the first release event has been consumed is low. Therefore, l defines the interaction length scale. Note that diffusion constant contributes by fifth root to l, and therefore, different diffusion constants in the community should not considerably alter patterns, as observed in simulations (Figure 3—figure supplement 1A). The calculated l of 50 μm is larger than the experimental patch size ${\lambda }_{\text{z}}^{*}$ of ∼10 to 20 μm. What could account for this discrepancy? Experimentally measured diffusion constant of Sulforhodamine 101 in community (20 µm²/s) reduces l to ∼30 µm. In addition, cells may take up more nutrients than what is required for producing one daughter and may store extra nutrients in vacuoles (Shou et al., 2007), further reducing the estimated interaction length scale.

The individual-based fitness model followed cell growth in a three-dimensional simulation grid of 100 × 100 × 300 cells with periodic boundary conditions along the x and y directions. Consider population i (i = R or G) interacting with population j (j = G or R). Without loss of generality, consider a focal cell from population i. The growth of the focal cell is influenced by cells in its cubic three-axial interaction neighborhood (Deutsch and Dormann, 2004) defined by l-cell-width to the left, right, front, back, above, and below. Let ${\varphi }_i$ and ${\varphi }_j$ be the fraction occupancy of i and j in the interaction neighborhood, respectively. The growth rate of the focal cell is $r_i=[r_{i0}+r_{ij}{\varphi }_j(1-{\varphi }_i\left)\right][1-\chi ({\varphi }_i+{\varphi }_j\left)\right]$. $r_{i0}$ is the basal fitness (growth rate of i without any interactions); $r_{ij}{\varphi }_j(1-{\varphi }_i)$ represents the fitness effect on i by j, which increases with partner abundance and decreases with recipient abundance due to intra-population competition for partner; $[1-\chi ({\varphi }_i+{\varphi }_j\left)\right]$ reflects intra- and inter-population competition for shared resources with fitness decreasing as the neighborhood becomes more occupied. Cells were inoculated in the bottom surface of the simulation grid. In each simulation time step Δt, the probability of cell division is r_i Δt. A cell would divide either to the side if there was space within its (x, y) planar confinement neighborhood of n-cell radius or upward otherwise (Figure 1—figure supplement 1). Parameters used are listed in Figure 1—source data 1. χ = 0.8, l = 3, and n = 2 in all cases. See Source code 1 for an example (and Source code 4 for the MATLAB function).

The individual-based diffusion model followed actions of cells (nutrient uptake, cell division, cell death, and possibly release of nutrients) and the distribution of nutrients in a three-dimensional simulation grid. Since cell division and death occur at a time-scale much longer than diffusion and nutrient uptake, we used a multi-grid scheme in both space and time. In this model, a three-dimensional simulation domain consisted of cell grids representing individual cells and nutrient grids representing nutrient concentrations (Figure 2—figure supplement 2). Cell actions and nutrient distributions were updated at discrete time steps over the simulation domain. Simulations were typically performed over an agarose domain of 0.75-mm length × 0.75-mm width × 24-mm depth and a community domain of 0.75-mm length × 0.75-mm width × 0.3-mm height with parameters listed in Figure 2—source data 1.

Nutrient concentrations as a function of space and time are based on the diffusion equation

with

Equation (3) states that S, the amount of limiting nutrient in a diffusion grid, depends on three processes: i) diffusion of nutrient with diffusion constant D, ii) uptake of nutrients (Walther et al., 2005) by cells (U), and iii) in cooperative communities, release of nutrients by the partner population (Q). In equations (4) and (5), n_u and n_q are the number of consuming and releasing cells within the diffusion grid, respectively, K_MM is the Michaelis-Menten constant for uptake, v_m is the maximum uptake rate per cell, and ρ is the release rate per cell. To solve this diffusion equation, we used a finite difference time–domain method (Crank, 1980) with no-flow ($\partial S/\partial z=0$) boundary conditions applied to the top and bottom surfaces of the simulation domain and periodic boundary conditions applied to the four vertical sides of the domain.

Cell growth rate r in the model is dictated by Monod's equation

in which r_m is the maximum growth rate when nutrients are abundant, S is the concentration of the limiting nutrient, and K_M is the S at which half maximal growth rate is achieved. We assume that individual cells take up nutrients with K_MM = K_M, and once they have accumulated the required amount of the limiting nutrient, cell division occurs.

To incorporate realistic assumptions about cell rearrangement upon division in the community, we monitored single yeast cells growing into microcolonies on solid media (Figure 1—figure supplement 1A). Initially, each cell budded in the same plane and pushed others in its immediate neighborhood to the side. Once a cell was completely surrounded on each side by roughly five cells, it budded upward. The same process was implemented in the diffusion model (Figure 2—figure supplement 2). It should be noted that by forcing the confined cells to bud strictly upward, the model underestimates intermixing. As a result, communities show more vertical features in simulations than experiments (Figure 2F).

Temporally, the update time-step for nutrient diffusion and uptake (∼1 s) was smaller than that for cell division and death (∼360 s). Spatially, cell actions took place on a cell grid with single-cell resolution (5 μm), while nutrient distributions were followed on a diffusion grid at lower spatial resolution (∼15 to 60 μm). These values were selected considering the trade-off between simulation time and accuracy, while ensuring the stability of simulations. For instance, at each time-step within a diffusion grid, the total amount of nutrients consumed should be considerably smaller than available nutrients. In other words, assuming that the length of a diffusion grid is nc where c is the length of a cell grid which is equivalent to the size of a cell, there are at most n³ cells in each diffusion grid and

where S is the concentration of the limiting nutrient in the spatial grid, K_MM is the Michaelis-Menten coefficient for nutrient uptake, v_m is the maximum uptake rate, and dt_u is the uptake time-step. After simplifying equation (7), we obtain

In the worst case scenario of S being much smaller than K_MM, using parameters of our engineered yeast strains, we find dt_u ∼ 0.5 s. Thus, we chose dt_u = 0.5–1 s as the time-step for updating the nutrient uptake and diffusion equations. To ensure stability of the finite-difference equations for diffusion, the diffusion grid-size (nc) and the time-step (dt_u) for the diffusion equation have to satisfy (Iserles, 2009)

where D is the diffusion constant in the region of interest. From this relation, we chose diffusion grid size nc = 50 μm, and consequently, each diffusion grid contains 10 × 10 × 10 cells (Figure 2—figure supplement 2A). Since diffusion is fast within each grid (∼1 s), this choice of grid size is unlikely to introduce a notable error in our calculations. See Source codes 2 and 3 for examples (and Source code 4 for the MATLAB function).

Let ${\varphi }_R$ and ${\varphi }_G$ be the fraction occupancy of R and G in an interaction neighborhood, respectively. ${\varphi }_{G*}$is the fraction occupancy of G that leads to equal fitness of the two populations and thereby result in a steady state ratio within the interaction neighborhood. Following the assumptions of the fitness model, the growth rate of each population is $r_i=[r_{i0}+r_{ij}{\varphi }_j(1-{\varphi }_i\left)\right][1-\chi ({\varphi }_i+{\varphi }_j\left)\right]$ where i = R or G and j = G or R. For simplicity, we assume that ${\varphi }_R+{\varphi }_G=1$, which leads to ${\widehat{r}}_i=r_i/(1-\chi )=r_{i0}+r_{ij}{\varphi }_j$.

In two-population cooperative communities, using the simplifying assumption of $r_{RG}=r_{GR}=r_{\mathrm{int}}>0$, we have ${\widehat{r}}_G-{\widehat{r}}_R=(r_{G0}-r_{R0}+r_{\mathrm{int}})-2r_{\mathrm{int}}{\varphi }_G$ (Figure 1—figure supplement 2B). Setting ${\widehat{r}}_G-{\widehat{r}}_R$ to 0, the community can achieve a steady-state value of ${\varphi }_{G*}=(r_{G0}-r_{R0}+r_{\mathrm{int}})/2r_{\mathrm{int}}$. To satisfy $0<{\varphi }_{G*}<1$, $r_{\mathrm{int}}>|r_{G0}-r_{R0}|$ which means that the interaction term has to be strong enough to overcome the difference in the basal fitness of the two populations. Here, ${\varphi }_{G*}$ is stable, since ${\widehat{r}}_G-{\widehat{r}}_R$ is positive (favoring G) when ${\varphi }_G<{\varphi }_{G*}$ and negative (favoring R) when ${\varphi }_G>{\varphi }_{G*}$. At ${\varphi }_{G*}$, two populations grow at the same rate and population ratio is fixed at R:G = (1 − ${\varphi }_{G*}$):${\varphi }_{G*}$.

Similar analysis shows the existence of a stable partner ratio under commensalism G[∼↑]R. A steady-state ratio ${\varphi }_{G*}$ can exist under exploitation. However, initial ratios below the critical value ${\varphi }_{G,c}$ (Figure 1—figure supplement 2B) will not converge to ${\varphi }_{G*}$. Other interactions ([∼ ∼], [∼ ↓], and [↓ ↓]) do not converge to a stable ratio. Similar conclusions on ratio convergence hold if we assume ${\widehat{r}}_i=r_{i0}+r_{ij}f\left({\varphi }_j\right)$ for any continuous function f that monotonically increases with ${\varphi }_j$ (proof not shown). Results are summarized in Table 1.

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

We would like to thank current and past Shou group members (Adam Waite, Chi-Chun Chen, David Skelding, Aric Capel, and Justin Burton), our Hutch colleagues (Harmit Malik, Daniel Gottschling, Suzannah Rutherford, Katie Peichel, Maulik Patel, Nitin Phadnis, Sarah Holte, William Hazelton, and Steve Andrews), our UW colleagues (Ben Kerr, David Stahl, Kristina Hillesland, Josephine Chandler, Ajai Dendekar, E. Peter Greenberg), Sri Ram, and Wallace Marshall for discussions. We are grateful to Kristina Hillesland for connecting the Shou and Fields groups.

© 2013, Momeni et al.

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