Introduction {#s1}
============

Biological interactions drive pattern formation at different levels of organization ([@bib34]), ranging from developmental patterning within multicellular organisms and biofilms ([@bib43]; [@bib29]; [@bib47]; [@bib7]), to ecological patterning within multi-species communities ([@bib27]; [@bib39]; [@bib31]). Patterning, reflecting the relative spatial positioning of individuals with respect to each other, can be critical for the proper functioning of a community. Consider microbial communities: in a synthetic community, three bacterial species, each contributing an essential benefit while simultaneously competing for these benefits, can only grow when they are separated by an intermediate distance ([@bib20]); different types of patterning are correlated with different levels of biofilm growth ([@bib6]; [@bib4]); branching colony morphology allows more effective spreading across a nutrient-poor surface ([@bib28]); and in waste treatment granules, the layered pattern of bacteria and archaea is thought to facilitate the sequential degradation of substrates ([@bib40]).

Despite the wide-ranging importance of microbial communities in, for example, human health and the biogeochemical cycling of elements, it is still unclear how cell--cell and cell--environment interactions govern the patterning of communities ([@bib13]). Understanding the mechanistic basis of pattern formation from observations of natural communities is stymied by the multitude of cell--cell and cell--environment interactions, as well as environmental variations within and across communities. Thus, it is not uncommon to observe qualitatively different patterns in samples of essentially the same type of community ([@bib6]; [@bib50]; [@bib11]).

To circumvent the lack of control in natural communities, we employed mathematical and experimental systems to systematically investigate how different types of ecological interactions might lead to distinct community patterning. Interactions can be classified into different ecological types based on their fitness effects on the interacting partners. We focused on the fitness effects rather than the molecular mechanisms of interactions, because diverse molecular mechanisms, ranging from physical associations in cell coaggregates and biofilms ([@bib22]) to chemical interactions such as quorum sensing ([@bib37]), toxin warfare ([@bib46]), and metabolite supply ([@bib6]), all have fitness consequences which can be positive, neutral, or negative. Among different ecological interactions, we have placed a special emphasis on strong cooperation, interactions with large positive fitness effects on both partners including obligatory cooperation. This is because 1) it is important in a wide variety of microbial communities ranging from syntrophic systems critical for nutrient cycling ([@bib42]; [@bib14]; [@bib30]) to pathogenic biofilms ([@bib19]; [@bib22]; [@bib13]); and 2) the codependence between cooperative partners poses special challenges for isolating and culturing cells ([@bib42]; [@bib30]).

We investigated patterning in three-dimensional communities grown from two fluorescently-marked populations of cells initially randomly distributed on top of a surface ([Figure 1A](#fig1){ref-type="fig"}). Starting with a generalized model based on fitness effects of ecological interactions between two populations (A and B) occurring at a local scale ('fitness model'), we predicted: 1) interactions benefiting at least one partner could potentially allow initially disparate partner ratios to converge over time, and 2) unlike other types of ecological interactions that caused partner segregation or layering of one population over the other (A over B or B over A), strongly cooperating partners intermixed by forming patches that successively accumulated on top of each other (A over B over A over B, etc). We tested these predictions experimentally in obligatory cooperative systems including engineered yeast communities and syntrophic methanogenic biofilms. Finally, we used the fitness model to show 'strongly cooperating partners intermix' could be generalized to communities consisting of multiple species.10.7554/eLife.00230.003Figure 1.The fitness model generates two ecological patterning predictions.(**A**) In all simulated and experimental communities (see 'Materials and methods'), two populations of cells, marked in red and green, were initially randomly distributed on a surface unless otherwise stated. The two populations engaged in one of the six ecological interactions. Population ratios for the entire community and patterns of vertical cross sections were examined. (**B**) The fitness model predicts that strong interactions beneficial to at least one partner can potentially lead to the convergence of initially disparate population ratios ([Figure 1---figure supplement 2](#fig1s2){ref-type="fig"}). (**C**)--(**H**) Representative vertical cross-sections of simulated communities, each engaging in one of the six types of ecological interactions, are presented. The fitness effects of ↑and ↓are large compared to the non-zero basal fitness of the recipient ([Figure 1---source data 1](#SD1-data){ref-type="supplementary-material"}), and therefore \[↑ ↑\] is strong facultative cooperation. (**I**) Vertical cross-sections of single-cell thickness from cooperative communities show significantly more intermixing than those from other communities (*n* = 28 sections; p\<0.01, Mann--Whitney *U* test). An intermixing index of 6 (red dotted line) or above separates cooperative from non-cooperative communities in our simulations. To reduce the correlation of sections sampled from the same community, nearest sections were separated by at least seven sections.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.003](10.7554/eLife.00230.003)10.7554/eLife.00230.004Figure 1---source data 1.Parameter values used in the fitness model.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.004](10.7554/eLife.00230.004)10.7554/eLife.00230.005Figure 1---figure supplement 1.Cell rearrangement in simulations follows experimental observations on cells that are not actively motile.(**A**) To incorporate realistic assumptions about cell rearrangement in three-dimensional communities, we monitored the growth of a single fluorescent yeast cell into a microcolony on top of solid agarose (Video 1). Starting from one budding cell (a, hour 0.5), cells initially divided and pushed other cells to the side such that all cells remained in the same plane (b, hour 9). When the group of cells reached approximately 5-cell radius (c, hour 14), cells in the center could no longer push other cells to the side and instead budded upward, as indicated by higher intensities in the fluorescence images. Continued growth of microcolony forced more cells in the middle to send their progeny to upper layers, while cells close to the edge could still push other cells to the side and remained on the agarose surface (d, hour 17.5). Scale bars are 20 μm. (**B**) In simulations of the diffusion and fitness models, a focal budding cell (star) was assumed to bud toward the nearest empty space within its coplanar (x,y) confinement neighborhood of n-cell radius (left panel, magenta circle). If the nearest empty space was not unique, a random choice was made. Cells along the path (yellow rim) to the nearest empty space shifted positions (left panel, cyan arrows) to accommodate the new cell. During rearrangement, if a cell moved to a new position with no cells immediately below it, the cell was lowered until it landed on another cell or the agarose surface. If the confinement neighborhood was filled, the cell budded upward and pushed up all cells above it (right panel).**DOI:** [http://dx.doi.org/10.7554/eLife.00230.005](10.7554/eLife.00230.005)10.7554/eLife.00230.006Figure 1---figure supplement 2.The fitness model predicts that convergence of population ratios is possible when an interaction benefits at least one partner.(**A**) The fitness model demonstrated that interactions conferring large benefit to at least one partner could lead to ratio convergence. Simulation parameters are listed in [Figure 1---source data 1](#SD1-data){ref-type="supplementary-material"}. (**B**) We explored conditions for community ratio convergence ('Requirements for steady-state ratios in the six types of communities' in 'Materials and methods') by looking for the existence of stable steady states in which the two partners in an interaction neighborhood would reach identical fitness ${\hat{r}}_{G} - {\hat{r}}_{R} = 0$ (solid circles). Open and solid circles denote unstable and stable steady states, respectively. Here, $\phi_{G*}$ is the steady state fraction occupancy of G so that $\phi_{G*}$:$\phi_{R*}$ is the steady state ratio.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.006](10.7554/eLife.00230.006)10.7554/eLife.00230.007Figure 1---figure supplement 3.Strong mutual antagonism can lead to rapid divergence of population ratios.The fitness model shows that if two partner populations inhibit each other sufficiently strongly compared to their basal fitness, population ratios may diverge rapidly.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.007](10.7554/eLife.00230.007)

Results {#s2}
=======

All microbial communities exhibit intra- and inter-population competition as cells compete for shared resources, including space. Therefore, for simplicity, we use \[∼ ∼\] to denote the 'baseline' competition. The fitness effect of baseline competition can start at zero if shared resources are in excess, but will eventually become negative as shared resources become limited. In 'addition' to competition, interactions between two partners can exert positive, negative, or no fitness effects on one or both partners. The 'net' fitness effects of all interactions between two partners, including competition, on the two partners can be represented as two symbols in a square bracket. There are six possibilities:\[∼ ∼\] (no fitness effects other than those from baseline competition), \[∼ ↑\] (commensalism, in which one partner experiences nothing more than competition whereas the other enjoys an overall fitness benefit even when competition has been taken into consideration), \[∼ ↓\] (amensalism), \[↓ ↓\] (mutual antagonism), \[↓ ↑\] (victim-exploiter), and \[↑ ↑\] (cooperation). The identities of partners may be added to the notation such that A\[↓ ↑\]B would mean that the overall interaction (including competition) inhibits A and promotes B. Under this notation, inter-population toxin-warfare ([@bib46]) would be \[↓ ↓\] while inter-population cooperation based on the exchange of distinct net benefits would be \[↑ ↑\]. Some interactions may at the first sight seem cooperative, but the net interaction may turn out not to be cooperative after considering the negative effects of competition ([@bib20]). In the game theory definition, cooperative acts incur fitness costs to the performers. Here we have taken a more liberal definition of cooperation to include mutually-beneficial interactions that may or may not involve fitness costs.

Predicting ecological patterning in simple communities using a generalized fitness model {#s2-1}
----------------------------------------------------------------------------------------

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](#fig1){ref-type="fig"}). Cells grew horizontally until sufficiently confined, at which point they grew upward to accommodate the birth of new cells ([Figure 1---figure supplement 1](#fig1s1){ref-type="fig"}, [Video 1](#video1){ref-type="other"}). 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.Yeast cells bud to the sides when there is available space and bud upward when sufficiently confined (corresponding to [Figure 1---figure supplement 1A](#fig1s1){ref-type="fig"}). To infer the process of cell rearrangement in three-dimensional communities, we monitored how single YFP-fluorescent yeast cells grew into microcolonies on top of solid agarose. Initially, dividing cells pushed other cells to the side such that all cells remained in the same plane. When a cell was sufficiently confined from the sides by other cells (approximately within a 5-cell radius), it could no longer bud to the side and instead budded upward, as indicated by higher intensities in the fluorescence images. Continued growth of microcolony forced more cells in the middle to send their progeny to upper layers, while cells close to the edge could still push other cells to the side and remain on the agarose surface. All images are taken with the same exposure time. Scale bar is 20 μm.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.008](10.7554/eLife.00230.008)10.7554/eLife.00230.008

We first analyzed the population composition of communities formed in the fitness model. Simulations ([Figure 1---figure supplement 2](#fig1s2){ref-type="fig"} and [Figure 1---source data 1](#SD1-data){ref-type="supplementary-material"}) 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](#fig1){ref-type="fig"}). 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](#fig2s1){ref-type="fig"}). This convergence of population ratios reflects a balance between supply and consumption when growth is limited by supply ([@bib45]): 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](#fig1){ref-type="fig"}, parameters in [Figure 1---source data 1](#SD1-data){ref-type="supplementary-material"}), 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](#fig1){ref-type="fig"}). In \[↑ ↓\] and \[∼ ↑\], frequently one of the populations (green) either formed a column or became covered by the partner population (red) ([Figure 1F--G](#fig1){ref-type="fig"}). 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](#fig1){ref-type="fig"}).

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](#fig1){ref-type="fig"}). Thus, we predicted that partner intermixing would distinguish strong cooperation from other ecological interactions.

Partner intermixing in engineered obligatory cooperative yeast communities and a corresponding diffusion model {#s2-2}
--------------------------------------------------------------------------------------------------------------

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](#fig2s1){ref-type="fig"}). Depending on their genetic backgrounds, the two strains engaged in either equal-fitness or unequal-fitness competition. In obligatory commensal communities, $R_{\rightarrow 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](#fig2s1){ref-type="fig"}). In obligatory cooperative communities (previously described as Cooperation that is Synthetic and Mutually Obligatory, or "CoSMO"), $R_{\rightarrow A}^{\leftarrow L}$ overproduced adenine and required lysine while $G_{\rightarrow L}^{\leftarrow A}$ overproduced lysine and required adenine ([@bib45]) ([Figure 2---figure supplement 1C](#fig2s1){ref-type="fig"}). 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 ([@bib45]).

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](#fig2s2){ref-type="fig"}) and that cell growth depended on the local concentration of its limiting metabolite according to Monod\'s equation ([@bib33]). Most parameters in the diffusion model were measured experimentally ([Figure 2---source data 1](#SD2-data){ref-type="supplementary-material"}). 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 ([@bib61]). In competitive communities, whether in the diffusion model or experiments, time-lapse top-views suggested population segregation. For equal-fitness competitive communities ([Figure 2A](#fig2){ref-type="fig"}, left; [Video 2](#video2){ref-type="other"}), 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](#fig1s1){ref-type="fig"}). 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](#fig2){ref-type="fig"}, right; [Video 3](#video3){ref-type="other"}), except that the fitter population expanded in the top view during growth (compare the green population in iv and v).10.7554/eLife.00230.009Figure 2.Obligatory cooperation, but not competition or obligatory commensalism, results in substantial partner intermixing in engineered yeast communities and in communities simulated using the diffusion model.Competitive communities of strains with equal fitness (equal-fitness competition, abbreviated as 'Eq-fitness Comp.') showed population segregation as suggested by static late-stage top-views (**A**, left) and columnar cross-section patterns (**D**, left). When competing strains had different fitnesses (unequal-fitness competition, abbreviated as 'Uneq-fitness Comp.'), the fitter population *G* expanded during growth, as evident in top views (**A**, right) and vertical cross-sections (**D**, right). In obligatory commensal communities, since one supplier could support the birth of multiple consumers, consumers eventually overgrew and covered suppliers (top-views in **B** and vertical cross-sections in **E**). Obligatory cooperative communities showed substantial population intermixing as suggested by alternating cell types in top-views (**C**, 6× magnification insets in experiments) and patchy cross-section patterns (**F**). Top views of communities from the diffusion model integrate intensity and color over height such that brighter colors represent higher cell numbers and yellowness indicates the simultaneous presence of green and red. Scale bar: 100 μm. All communities started from total 500 cells/mm^2^ and *R*:*G* = 1:1. (**G**) Quantification of intermixing in experimental (symbols) and diffusion-model (lines) communities showed that while the intermixing index remained low for commensal (brown) and competitive (grey and black) communities, it increased linearly with community height in obligatory cooperative (magenta) communities. Results from the diffusion model underestimated intermixing indexes because a confined cell was modeled to divide strictly vertically upward ([Figure 1---figure supplement 1](#fig1s1){ref-type="fig"}) without allowing cell movements in horizontal directions ([Figure 3---figure supplement 1F](#fig3s1){ref-type="fig"}). (**H**) A conceptual model illustrates the development of intermixing over time in a strongly cooperative community with 1:1 steady-state population ratio. Local deviations from the steady-state ratio result from asymmetric partner properties and/or stochastic fluctuations in cell numbers (i). The under-abundant population (red) grows faster than its neighboring over-abundant partner (green). Due to the spatial localization of large cooperative benefits, cells near population borders grow faster than those farther away. Consequently, cells from the initially under-abundant red population at the border divide the fastest. Progeny that pile on the green partner have more access to cooperative benefits than those who do not (ii), which favors intermixing. When the previously over-abundant partner becomes under-abundant (iii), piling-up in the opposite direction occurs (iv, green on red).**DOI:** [http://dx.doi.org/10.7554/eLife.00230.009](10.7554/eLife.00230.009)10.7554/eLife.00230.010Figure 2---source data 1.Definitions and values of parameters used in the diffusion model.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.010](10.7554/eLife.00230.010)10.7554/eLife.00230.011Figure 2---figure supplement 1.In engineered yeast communities, obligatory cooperation and obligatory commensalism allow initially different partner ratios to converge over time.Engineered yeast communities consisted of two non-mating fluorescent strains of yeast, *G* and *R*, engaged in a metabolic interaction. (**A**) In competitive communities, two prototrophic strains engaged in equal-fitness or unequal-fitness competition for shared resources. (**B**) In obligatory commensal communities, $R_{\rightarrow A}^{\leftarrow L}$ took in lysine from the media and overproduced adenine to feed $G^{\leftarrow A}$. (**C**) In obligatory cooperative communities, $R_{\rightarrow A}^{\leftarrow L}$ overproduced adenine and required lysine while $G_{\rightarrow L}^{\leftarrow A}$ overproduced lysine and required adenine. (**D**)--(**H**) Symbols: experiments; lines: results from the diffusion model. (**D**) In competitive communities, the ratio between two equally-fit populations remained steady (grey), whereas the ratio between two populations with unequal fitness (*G* 20% fitter than *R*) changed monotonically (black) in favor of the fitter population. Population ratios of obligatory commensal (**E**) and obligatory cooperative (**F**) communities converged regardless of their initial values (different colors representing different initial conditions). (**G--H**) Ratio fluctuations observed in experimental cooperative communities were better recapitulated in diffusion models that incorporated additional specifics of the yeast system. (**G**) Incorporating a slower diffusion constant in the community (∼20 μm^2^/s) than that in agarose (360 μm^2^/s) allowed ratio convergence to more closely resemble experimental data. (**H**) Incorporating experimental observations, such as better starvation tolerance of adenine-requiring cells compared to lysine-requiring cells and the consequent delay in lysine release compared to adenine release ([@bib45]), allowed ratio convergence to more closely resemble experimental data. All communities were initiated at 3000 cells/mm^2^ on top of a 24-mm-high column of agarose. Population ratios in experiments were measured using flow cytometry (see 'Materials and methods'). As predicted by the fitness model, population ratios converged as time progressed in obligatory cooperative and obligatory commensal but not competitive communities both in the experiments and the diffusion model.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.011](10.7554/eLife.00230.011)10.7554/eLife.00230.012Figure 2---figure supplement 2.Basic assumptions in the diffusion model.In the diffusion model, community grew on the top surface of an agarose column (light blue). Cell grids (black borders) represented cells at different states: live R (light red), live G (light green), dead R (dark red), dead G (dark green), and no cells (white). (**A**) When the diffusion constants of nutrients in the yeast community were assumed to equal those in agarose (360 μm^2^/s), nutrient grids (orange) were used to simulate the distribution of nutrients. (**B**) When the diffusion constants of nutrients in the yeast community were modeled to be that of the fluorescent dye Sulforhodamine 101 (∼20 μm^2^/s; [Figure 2---figure supplement 1G](#fig2s1){ref-type="fig"}), we implemented community diffusion and agarose diffusion grids of 3- and 12-cell grid width, respectively, to accommodate two diffusion constants.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.012](10.7554/eLife.00230.012)10.7554/eLife.00230.013Figure 2---figure supplement 3.Cooperative communities exhibit a characteristic patch size associated with the spatial localization of benefits.(**A**) Starting from densities and ratios spanning orders of magnitude, experimental obligatory cooperative yeast communities developed a consistent characteristic vertical patch size $\text{λ}_{\text{z}}^{*}$ of ∼10 to 20 μm. (**B**) The diffusion model shows that when both released nutrients were instantly and uniformly distributed throughout a community ('Coop. Instant Distr.'), intermixing was diminished compared to when nutrients diffused in the community (360 µm^2^/s, 'Coop.'). Excessive release of nutrients (200× compared to typical values in [Figure 2---source data 1](#SD2-data){ref-type="supplementary-material"}) also significantly reduced population intermixing ('Coop. Excessive Release'), since the benefits exceeded the uptake capacity of partner cells and therefore no longer remained localized. All simulations started with 2000 total cells/mm^2^ and 1:1 population ratio. Box plots show the 25th to 75th percentile range, with the median marked with a line and whiskers extending to the 95% confidence interval.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.013](10.7554/eLife.00230.013)10.7554/eLife.00230.014Figure 2---figure supplement 4.Obligatory cooperative yeast partners intermix.Vertical cross-sections in (**A**) and (**B**) were from different growth stages of obligatory cooperative yeast communities starting at total 500 cells/mm^2^ and 1:1 population ratio. Vertical cross-section in (**C**) and (**D**) were respectively from the center and the edge of an obligatory cooperative yeast community started from a high-density inoculum and allowed to expand to new territories. Areas enclosed in the yellow boxes were 3× magnified in the right panels to show details of intermixing. The top-view schematic insets in (**C**) and (**D**) show the relative position of the corresponding cross-section with regard to the initial inoculum (dashed line). All scale bars are 100 µm.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.014](10.7554/eLife.00230.014)Video 2.Top views of an equal-fitness competitive community suggest population segregation (corresponding to [Figure 2A](#fig2){ref-type="fig"}, left). Competitive communities of strains with equal fitness showed population segregation as suggested by static late-stage top-views. The community started from a uniform distribution of total 500 cells/mm^2^ and R:G = 1:1. Intensities of both fluorescent channels in images at different times were normalized to the same maximum value for better representation of patterns throughout growth. Blanks in the video were due to removal of the dish to sample other replicate communities for flow-cytometry or sectioning. Scale bar is 100 μm.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.015](10.7554/eLife.00230.015)10.7554/eLife.00230.015Video 3.In unequal-fitness competition, the fitter population gradually covers the less fit population (corresponding to [Figure 2A](#fig2){ref-type="fig"}, right). Here, G is fitter than R. The community started from a uniform distribution of total 500 cells/mm^2^ and R:G = 1:1. Intensities of both fluorescent channels in images at different times were normalized to the same maximum value for better representation of patterns throughout growth. Scale bar is 100 μm.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.016](10.7554/eLife.00230.016)10.7554/eLife.00230.016

In obligatory commensal communities, time-lapse top-views suggested a 'layered' pattern with one population covering the other ([Figure 2B](#fig2){ref-type="fig"}; [Video 4](#video4){ref-type="other"}). We found that in both the diffusion model and experiments, $R_{\rightarrow 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_{\rightarrow A}^{\leftarrow L}$ before $R_{\rightarrow A}^{\leftarrow L}$ became abundant (compare i, ii, and iii). During the growth of $R_{\rightarrow A}^{\leftarrow L}$, some $G^{\leftarrow A}$ cells had been pushed to the top layer of the community (iii). As the $R_{\rightarrow A}^{\leftarrow L}$ population continued to expand and release adenine, $G^{\leftarrow A}$ started to grow rapidly (compare iii and iv). Eventually, $R_{\rightarrow 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_{\rightarrow A}^{\leftarrow L}$ cell could support the birth of multiple $G^{\leftarrow A}$ cells ([@bib45]), $G^{\leftarrow A}$ population outnumbered and covered $R_{\rightarrow A}^{\leftarrow L}$ (v).Video 4.Top views of an obligatory commensal community suggest population layering (corresponding to [Figure 2B](#fig2){ref-type="fig"}). For a detailed explanation of the growth kinetics of the community ($R_{\rightarrow A}^{\leftarrow L}\left\lbrack \operatorname{\sim\uparrow} \right\rbrack G^{\leftarrow A}$), please refer to the main text. The community started from a uniform distribution of total 500 cells/mm^2^ and $R_{\rightarrow A}^{\leftarrow L}$:$G^{\leftarrow A} = 1:1$. Intensities of both fluorescent channels in images at different times were normalized to the same maximum value for better representation of patterns throughout growth. Blanks in the video were due to removal of the dish to sample other replicate communities for flow-cytometry or sectioning. Scale bar is 100 μm.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.017](10.7554/eLife.00230.017)10.7554/eLife.00230.017

In contrast, time-lapse top-views of obligatory cooperative communities suggested population intermixing ([Figure 2C](#fig2){ref-type="fig"}; [Video 5](#video5){ref-type="other"}). After plating on agarose lacking adenine and lysine (i), cells in both populations divided once or twice by utilizing metabolites stored in their vacuoles ([@bib45]). $R_{\rightarrow A}^{\leftarrow L}$ released adenine and entered the death phase while $G_{\rightarrow L}^{\leftarrow A}$ continued to grow by utilizing the released adenine (ii) ([@bib45]). $G_{\rightarrow L}^{\leftarrow A}$ entered the death phase and released lysine after a significant delay due to their better starvation tolerance compared to $R_{\rightarrow A}^{\leftarrow L}$ ([@bib45]). Released lysine supported the growth of surviving $R_{\rightarrow A}^{\leftarrow L}$ cells into microcolonies (iii). $R_{\rightarrow A}^{\leftarrow L}$ in turn released adenine which promoted growth of nearby $G_{\rightarrow L}^{\leftarrow A}$ cells and led to partial covering up of $R_{\rightarrow A}^{\leftarrow L}$ microcolonies by rapidly growing $G_{\rightarrow L}^{\leftarrow A}$ (iv and v, insets). Abundant $G_{\rightarrow L}^{\leftarrow A}$ cells provided enough lysine for rapid growth of local $R_{\rightarrow A}^{\leftarrow L}$ cells, which gave rise to patches of $R_{\rightarrow A}^{\leftarrow L}$ cells on top of the community (vi, inset). These community growth kinetics were consistent with previous measurements in liquid cultures ([@bib45]).Video 5.Top views of an obligatory cooperative community suggest populations intermixing (corresponding to [Figure 2C](#fig2){ref-type="fig"}). For a detailed explanation of the growth kinetics of the community ($R_{\rightarrow A}^{\leftarrow L}\left\lbrack \operatorname{\uparrow\ \uparrow} \right\rbrack G_{\rightarrow L}^{\leftarrow A}$), please refer to the main text. The community started from a uniform distribution of total 500 cells/mm^2^ and $R_{\rightarrow A}^{\leftarrow L}$:$G_{\rightarrow L}^{\leftarrow A} = 1:1$. Intensities of both fluorescent channels in images at different times were normalized to the same maximum value for better representation of patterns throughout growth. Blanks in the video were due to removal of the dish to sample other replicate communities for flow-cytometry or sectioning. The video started at ∼300 hr, after the formation of $R_{\rightarrow A}^{\leftarrow L}$ microcolonies. Scale bar is 100 μm.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.018](10.7554/eLife.00230.018)10.7554/eLife.00230.018

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](#fig2){ref-type="fig"}, left). Unequal-fitness ([Figure 2D](#fig2){ref-type="fig"}, right) and obligatory commensal communities ([Figure 2E](#fig2){ref-type="fig"}) 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](#fig2){ref-type="fig"} and [Figure 2---figure supplement 4A and B](#fig2s4){ref-type="fig"}). 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](#fig3s1){ref-type="fig"}).

To compare levels of intermixing in different communities, we quantified the intermixing index of each community ([Figure 2G](#fig2){ref-type="fig"}, 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](#fig2){ref-type="fig"}, grey). In unequal-fitness competition ([Figure 2G](#fig2){ref-type="fig"}, black) and obligatory commensalism ([Figure 2G](#fig2){ref-type="fig"}, 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](#fig2){ref-type="fig"}, 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](#fig2s3){ref-type="fig"}). 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](#fig2s3){ref-type="fig"}). 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](#fig2){ref-type="fig"}). 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.

Strongly cooperative partners intermix under a wide range of conditions {#s2-3}
-----------------------------------------------------------------------

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](#fig3s1){ref-type="fig"}). 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](#fig3s1){ref-type="fig"}, 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](#fig3s1){ref-type="fig"}) 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](#fig3){ref-type="fig"}, 'Center' and [Figure 2---figure supplement 4C](#fig2s4){ref-type="fig"}), in edge sections, significant intermixing was only observed in cooperative communities ([Figure 3A](#fig3){ref-type="fig"}, 'Edge' and [Figure 2---figure supplement 4D](#fig2s4){ref-type="fig"}).10.7554/eLife.00230.019Figure 3.Strongly cooperating populations intermix under a wide range of conditions.(**A**) In engineered yeast communities, even though both obligatory cooperative and non-cooperative communities directly above the high-density inoculation spot showed high population intermixing ('Center'), edge sections ('Edge') of obligatory cooperative communities showed significantly more intermixing than those from non-cooperative communities (Mann--Whitney *U* test, *p*\<5 × 10^−5^). A total of 10^6^ *R* and *G* cells at a 1:1 ratio were deposited in an inoculation spot of ∼2 mm^2^, corresponding to 10 cell layers. Communities were allowed to grow and expand beyond the inoculation spot on a 4-mm-tall agarose pad of 500 mm^2^ area. Vertical sections from the edges were taken at a height approximately half of the maximum community height at the center. Box plots show the 25th to 75th percentile range, with the median marked with a line and whiskers extending to the 95% confidence interval. (**B**) In simulations using the fitness model, facultative cooperation conferring smaller fitness benefits required greater community heights to exhibit a significant level of intermixing. The strength of facultative cooperation is shown as the ratio of fitness benefit received from each cooperative partner cell in the interaction neighborhood to the basal fitness of the focal cell. Simulation parameters can be found in [Figure 3---source data 1](#SD3-data){ref-type="supplementary-material"}. Error bars indicate 95% confidence interval.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.019](10.7554/eLife.00230.019)10.7554/eLife.00230.020Figure 3---source data 1.Parameter values used in the fitness model in [Figure 3B](#fig3){ref-type="fig"}.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.020](10.7554/eLife.00230.020)10.7554/eLife.00230.021Figure 3---figure supplement 1.Intermixing is observed in obligatory cooperative communities over a wide range of conditions.(**A**) Partner intermixing is insensitive to diffusion constants. In the diffusion model, obligatory cooperative communities with diffusion constants ranging from 20 to 360 µm^2^/s (corresponding to the diffusion constants of Sulforhodamine 101 in yeast community and agarose, respectively) both showed significant partner intermixing during growth. (**B**) Partner intermixing is insensitive to the kinetics of interactions. With or without incorporating details such as significantly better starvation tolerance in adenine-requiring cells compared to lysine-requiring cells and the consequent delay in lysine release compared to adenine release ([@bib45]), the diffusion model (lines) yielded substantial population intermixing. Diffusion constant within the community was assumed to be 20 μm^2^/s. In both (**A**) and (**B**), communities started from 3000 total cells/mm^2^ at 1:1 population ratio. (**C**) Intermixing does not depend on asymmetry between cooperating partners. We examined symmetric obligatory cooperation by assuming that for the two partners, all live cells released nutrients at a constant and identical rate without delay, and that other properties such as rates of growth, death, and nutrient consumption were also identical. Partner intermixing in this symmetric obligatory cooperative system (Sym. coop) was similar to that in asymmetric obligatory yeast cooperative system (Coop). Simulated communities started from total 500 cells/mm^2^ at 1:1 population ratio. (**D**) Cooperative intermixing is insensitive to initial ratios. Engineered yeast communities were initiated at a total cell density of 3000 cells/cm2 and grown to near the carrying capacity of the agarose column. Box plots show 25th to 75th percentile values of the intermixing index, with the median marked with a line and whiskers extending to 95% confidence interval. To facilitate comparisons, intermixing indexes in communities with different heights were linearly normalized to correspond to a height of ∼40 cell layers (200 µm). (**E**) Experimental obligatory cooperative communities (magenta) showed more intermixing than obligatory commensal (brown), equal-fitness competitive (grey), and unequal-fitness competitive (black) communities, except when the initial cell densities exceeded 10% confluence (yellow shading). Large intermixing indexes at high initial cell densities were due to the fact that the budding direction of a cell under spatial confinement was not exactly vertical (**F**). To facilitate comparisons, intermixing indexes in communities with different heights were linearly normalized to correspond to a height of ∼30 cell layers (150 µm). Communities were initiated at R:G = 1:1 and grown until near the carrying capacity of the agarose. Box plots show the 25th to 75th percentile range, with the median marked with a line and whiskers representing 95% confidence intervals. (**F**) An equal-fitness competitive community was initiated at confluent density (50000 cells/mm^2^) and R:G = 100:1 to ensure that single green cells were laterally confined by red cells. Cross-sections showed that the early progeny of a single cell was distributed in a division zone of ∼3-cell rather than 1-cell width. Lateral budding by late progeny further widens the zone to ∼5-cell diameter in a mature community. Scale bar is 20 µm.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.021](10.7554/eLife.00230.021)

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](#fig2s3){ref-type="fig"}). 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](#fig3s1){ref-type="fig"}). 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_{\rightarrow L}^{\leftarrow A}$ supplying lysine benefits to the partner ([Figure 3---figure supplement 1B](#fig3s1){ref-type="fig"}). Furthermore, even though asymmetry between properties of partners ([Figure 2C](#fig2){ref-type="fig"}) would seem to facilitate intermixing ([Figure 2H](#fig2){ref-type="fig"}), 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](#fig3s1){ref-type="fig"}).

Using the fitness model, we found that in facultative cooperation, small fitness benefits generated less intermixing than large fitness benefits ([Figure 3B](#fig3){ref-type="fig"}). 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.

Communities of naturally mutualistic microbes exhibit intermixing {#s2-4}
-----------------------------------------------------------------

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~2~, while *M. maripaludis* promotes the growth of *D. vulgaris* by consuming the H~2~ gas which can be inhibitory to *D. vulgaris* under these conditions ([Figure 4A](#fig4){ref-type="fig"}). 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 ([@bib8]; [@bib41]).10.7554/eLife.00230.022Figure 4.Obligatory cooperation through redox-coupling leads to partner intermixing.(**A**) In the absence of sulfate and hydrogen, the bacterium *Desulfovibrio vulgaris* (Dv) and the archaeon *Methanococcus maripaludis* (Mm) cooperate through redox coupling. Dv ferments lactate and produces mainly acetate, CO~2~, and H~2~. However, this reaction is not thermodynamically favorable unless H~2~ is kept at very low concentrations. H~2~ is used by Mm to reduce CO~2~ to methane. (**B**) In cooperative biofilms of Dv (green) and Mm (red), the intermixing index increased with height. Cell identification relied on FISH (see 'Materials and methods'). Scale bar: 20 µm.**DOI:** [http://dx.doi.org/10.7554/eLife.00230.022](10.7554/eLife.00230.022)

Vertical cross-sections of independent *D. vulgaris*--*M. maripaludis* biofilms indeed exhibited increasing intermixing as a function of community height ([Figure 4B](#fig4){ref-type="fig"}). 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 ([@bib5]) and those colonizing teeth ([@bib36]) also seemed to display intermixed patterns, although we do not know whether the intermixing indexes of these communities increased linearly with height.

Most of the strongly-cooperative pairs intermix in simulated multi-species communities {#s2-5}
--------------------------------------------------------------------------------------

In communities with more than two species, indirect interactions can obscure direct interactions ([@bib51]). 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](#fig5){ref-type="fig"} and [Figure 5---figure supplement 1](#fig5s1){ref-type="fig"}). 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](#fig5){ref-type="fig"}). Occasionally, commensal pairs ([Figure 5B](#fig5){ref-type="fig"}, ② \[∼ ↑\] ③) showed substantial intermixing (3 out of 58 commensal pairs, [Figure 5D](#fig5){ref-type="fig"}), or cooperative pairs (e.g., [Figure 5C](#fig5){ref-type="fig"}, ① \[↑ ↑\] ③) showed little intermixing (5 out of 31, [Figure 5D](#fig5){ref-type="fig"}), presumably due to interactions through other community members ([Figure 5B](#fig5){ref-type="fig"}, ③ indirectly promoted ② through ④; [Figure 5C](#fig5){ref-type="fig"}, ③ 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.10.7554/eLife.00230.023Figure 5.Most of the strongly-cooperative pairs intermixed in simulated five-species communities.(**A-C**) Examples of networks in which cooperative pairs intermixed (**A**), non-cooperative pairs intermixed (**B**), or cooperative pairs did not intermix (**C**) are shown. In the schematic network diagrams, line termini of →, ⊣, and --- represent ↑, ↓, and **∼**, respectively; cooperative pairs are highlighted in magenta. Simulations were performed using the fitness model. The basal fitness for each population was chosen randomly from a range spanning 0.03--0.05/hr. The fitness effect from each partner cell in the interaction neighborhood was either 0 for ∼, or otherwise randomly chosen to be approximately 2--3% of basal fitness to achieve strong interaction. Initial population ratios were randomly assigned such that no population was initially lower than 5% of the total population. Simulations were run for 10 generations, and vertical cross-sections of the final communities were examined for intermixing. We considered intermixing index exceeding a threshold of 6 ([Figure 1I](#fig1){ref-type="fig"}) as significant (red dotted lines). The remaining nine panels in [Figure 5B and C](#fig5){ref-type="fig"} are provided in [Figure 5---figure supplement 1](#fig5s1){ref-type="fig"}. (**D**) Quantifying intermixing in a total of 240 interactions from 24 independent communities showed that most of the cooperative pairs intermixed (magenta). Indirect interactions through other community members could make cooperative pairs not intermixed (grey) or non-cooperative (commensal) pairs intermixed (brown).**DOI:** [http://dx.doi.org/10.7554/eLife.00230.023](10.7554/eLife.00230.023)10.7554/eLife.00230.024Figure 5---figure supplement 1.The complete results of [Figure 5B and C](#fig5){ref-type="fig"} (panel A and B, respectively).**DOI:** [http://dx.doi.org/10.7554/eLife.00230.024](10.7554/eLife.00230.024)

Discussion {#s3}
==========

Using controlled systems to identify patterns driven by ecological interactions {#s3-1}
-------------------------------------------------------------------------------

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](#fig1s2){ref-type="fig"} and [Figure 2---figure supplement 1](#fig2s1){ref-type="fig"}). 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](#fig1s2){ref-type="fig"}; '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](#fig2s1){ref-type="fig"}). 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](#fig1s1){ref-type="fig"}). 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](#fig2s3){ref-type="fig"} and [Figure 3B](#fig3){ref-type="fig"}), and otherwise appears to be robust against variations in initial conditions ([Figure 2---figure supplement 3A](#fig2s3){ref-type="fig"}, [Figure 3A](#fig3){ref-type="fig"}, and [Figure 3---figure supplement 1D](#fig3s1){ref-type="fig"}) or reaction dynamics ([Figure 3---figure supplement 1A--C](#fig3s1){ref-type="fig"}). We have observed cooperative intermixing in the fitness model simulating strong facultative cooperation ([Figure 1](#fig1){ref-type="fig"}) and in obligatory experimental systems including engineered yeast communities ([Figure 2](#fig2){ref-type="fig"}) and syntrophic biofilms ([Figure 4](#fig4){ref-type="fig"}). 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](#fig2s3){ref-type="fig"}; '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](#fig1 fig2){ref-type="fig"}) 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 ([@bib60]).

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 ([@bib36]; [@bib25]; [@bib2]; [@bib16]; [@bib3]). *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 ([@bib35]). 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 ([@bib5]). Cooperative oral bacteria species intermixed when grown on saliva ([@bib36]; [@bib201]). 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.

Deviations from the expected ecological patterns suggest other major patterning forces {#s3-2}
--------------------------------------------------------------------------------------

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 ([@bib24]) and was adopted to model biofilm growth ([@bib52]). However, it is conceivable that if at least one population actively moves ([@bib17]) 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 ([@bib21]) showed modest intermixing even though the expectation for equal-fitness competition is a columnar pattern with an intermixing index close to zero ([Figure 2D](#fig2){ref-type="fig"}, 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 ([@bib21]).

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.

What can cause variability in patterning? {#s3-3}
-----------------------------------------

Variability in patterning has been observed within and between communities ([@bib6]; [@bib50]; [@bib11]), even if they were grown in laboratory-controlled environments ([@bib6]). 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 ([@bib16]; [@bib23]). 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](#fig1s3){ref-type="fig"}). 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 ([@bib26]).

In summary, our work is conceptually analogous to that of the competitive exclusion principle ([@bib15]). 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.

Materials and methods {#s4}
=====================

Engineering yeast strains {#s4-1}
-------------------------

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_{\rightarrow 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_{\rightarrow A}^{\leftarrow L}$ and $G_{\rightarrow 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.

Culturing communities {#s4-2}
---------------------

For yeast communities, agarose columns were prepared by pouring SD minimal medium ([@bib44]) 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 ([@bib48]), a modified basal salt medium without choline chloride. Headspace was sparged with anoxic 80% N~2~:20% CO~2~, 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 {#s4-3}
--------------

Flow cytometry was performed in a modified 4-laser FACS Calibur machine (DxP; CyTek Development, Fremont, CA). Cells in communities were diluted to ∼10^6^ cells/ml in H~2~O. Each sample (90 µl) was supplemented with 10 µl of fluorescent bead stock (Thermo Scientific Fluoro-Max Cat\# R0300 at ∼8 × 10^6^/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^4^--10^5^ 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.

Fluorescence imaging {#s4-4}
--------------------

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

Cryosectioning {#s4-5}
--------------

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 ([@bib38]) 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 ([@bib32]).

Fluorescence in situ hybridization (FISH) {#s4-6}
-----------------------------------------

Biofilms of *D. vulgaris* and *M. maripaludis* were fixed in 4% paraformaldehyde for 4 hr, and then embedded in polyacrylamide ([@bib10]). 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 ([@bib62]) 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~2~O and quickly dried with compressed air ([@bib1]). 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.

Spatial analysis {#s4-7}
----------------

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$$IM = {\left( {\sum\limits_{x_{i}}{c\left( x_{\mathbf{i}} \right)\, h\left( x_{\mathbf{i}} \right)}} \right)/\left( {\sum\limits_{x_{i}}{h\left( x_{\mathbf{i}} \right)}} \right)}.$$

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$$\lambda_{\text{z}}^{*} = {\left( {\sum\limits_{x_{i}}{\left\lbrack h\left( x_{i} \right) \right\rbrack^{2}/\left\lbrack 1 + c\left( x_{i} \right) \right\rbrack}} \right)/\left( {\sum\limits_{x_{i}}{h\left( x_{i} \right)}} \right)}.$$

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.

Calculating the characteristic patch size in cooperative yeast communities {#s4-8}
--------------------------------------------------------------------------

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*~ = (2*l/c*)^3^, 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*^2^/2*D*~1~, where *D*~1~ = 360 μm^2^/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$$\beta_{L} = N_{u}vt_{c} \leq \left( \frac{2l}{c} \right)^{3}\frac{\alpha_{L}}{T}\frac{l^{2}}{2D_{1}};$$thus,$$l \geq \left\lbrack {\frac{\beta_{L}}{\alpha_{L}}\frac{D_{1}Tc^{3}}{4}} \right\rbrack^{1/5} \approx 50\,\mu\text{m}.$$

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](#fig3s1){ref-type="fig"}). 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^2^/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 ([@bib45]), further reducing the estimated interaction length scale.

The fitness model {#s4-9}
-----------------

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 ([@bib12]) defined by *l*-cell-width to the left, right, front, back, above, and below. Let $\phi_{i}$ and $\phi_{j}$ be the fraction occupancy of *i* and *j* in the interaction neighborhood, respectively. The growth rate of the focal cell is $r_{i} = \left\lbrack r_{i0} + r_{ij}\phi_{j}\left( 1 - \phi_{i} \right) \right\rbrack\left\lbrack 1 - \chi\left( \phi_{i} + \phi_{j} \right) \right\rbrack$. $r_{i0}$ is the basal fitness (growth rate of *i* without any interactions); $r_{ij}\phi_{j}\left( 1 - \phi_{i} \right)$ 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; $\left\lbrack 1 - \chi\left( \phi_{i} + \phi_{j} \right) \right\rbrack$ 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](#SD1-data){ref-type="supplementary-material"}). Parameters used are listed in [Figure 1---source data 1](#SD1-data){ref-type="supplementary-material"}. χ = 0.8, *l* = 3, and *n* = 2 in all cases. See [Source code 1](#SD4-data){ref-type="supplementary-material"} for an example (and [Source code 4](#SD4-data){ref-type="supplementary-material"} for the MATLAB function).

The diffusion model {#s4-10}
-------------------

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](#fig2s2){ref-type="fig"}). 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](#SD2-data){ref-type="supplementary-material"}.

Nutrient concentrations as a function of space and time are based on the diffusion equation$$\frac{\partial S}{\partial t} = \nabla \cdot \left( D\nabla S \right) - U + Q,$$with$$U = v_{m}\frac{S}{S + K_{MM}}n_{u}$$$$Q = \rho\, n_{q}.$$

[Equation (3)](#equ3){ref-type="disp-formula"} 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 ([@bib49]) by cells (*U*), and iii) in cooperative communities, release of nutrients by the partner population (*Q*). In [equations (4)](#equ4){ref-type="disp-formula"} and [(5)](#equ5){ref-type="disp-formula"}, *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 ([@bib9]) 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$$r\left( S \right) = r_{m}\frac{S}{S + K_{M}},$$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](#fig1s1){ref-type="fig"}). 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](#fig2s2){ref-type="fig"}). 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](#fig2){ref-type="fig"}).

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*^3^ cells in each diffusion grid and$$S\left( nc \right)^{3} > {{n^{3}v_{m}}\frac{S}{S + K_{MM}}{dt_{u},}}$$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)](#equ7){ref-type="disp-formula"}, we obtain$$dt_{u} < \frac{S + K_{MM}}{v_{m}}c^{3}$$

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 ([@bib18])$$dt_{u} < \frac{1}{2}\frac{\left( nc \right)^{2}}{D},$$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](#fig2s2){ref-type="fig"}). 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](#SD5-data SD6-data){ref-type="supplementary-material"} for examples (and [Source code 4](#SD4-data){ref-type="supplementary-material"} for the MATLAB function).

Requirements for steady-state ratios in the six types of communities {#s4-11}
--------------------------------------------------------------------

Let $\phi_{R}$ and $\phi_{G}$ be the fraction occupancy of *R* and *G* in an interaction neighborhood, respectively. $\phi_{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} = \left\lbrack r_{i0} + r_{ij}\phi_{j}\left( 1 - \phi_{i} \right) \right\rbrack\left\lbrack 1 - \chi\left( \phi_{i} + \phi_{j} \right) \right\rbrack$ where *i* = *R* or *G* and *j* = *G* or *R*. For simplicity, we assume that $\phi_{R} + \phi_{G} = 1$, which leads to ${\hat{r}}_{i} = {r_{i}/\left( 1 - \chi \right)} = r_{i0} + r_{ij}\phi_{j}$.

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

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

10.7554/eLife.00230.025

###### 

Summary of steady-state occupancy, conditions to achieve steady-state, and the stability of steady-state for six types of ecological interactions.

**DOI:** [http://dx.doi.org/10.7554/eLife.00230.025](10.7554/eLife.00230.025)

  Interaction     Steady-state occupancy                                                              Steady-state condition                             Stability
  --------------- ----------------------------------------------------------------------------------- -------------------------------------------------- --------------------------------------------------------------------------------------------------------------
  *G*\[↑ ↑\]*R*   $\phi_{G*} = \frac{r_{G0} - r_{R0} + r_{int}}{2r_{int}}$                            $r_{int} > \left| r_{G0} - r_{R0} \right|$         Globally stable
  *G*\[∼↓\]*R*    $\phi_{G*} = \sqrt{\frac{r_{R0} - r_{G0}}{r_{int}}}$                                $r_{int} > {r_{R0} - r_{G0}} > 0$                  Unstable
  *G*\[∼↑\]*R*    $\phi_{G*} = \sqrt{\frac{r_{G0} - r_{R0}}{r_{int}}}$                                $r_{int} > {r_{G0} - r_{R0}} > 0$                  Globally stable
  *G*\[↓↓\]*R*    $\phi_{G*} = \frac{r_{R0} - r_{G0} + r_{int}}{2r_{int}}$                            $r_{int} > \left| r_{G0} - r_{R0} \right|$         Unstable
  *G*\[↓↑\]*R*    $\phi_{G*} = \frac{1}{2} + \sqrt{\frac{r_{G0} - r_{R0}}{2r_{int}} - \frac{1}{4}}$   $\begin{array}{l}                                  Locally stable (when $\phi_{G} > \phi_{G,c}$; see [Figure 1---figure supplement 2](#fig1s2){ref-type="fig"})
                                                                                                      {r_{int} > {r_{G0} - r_{R0}} > 0} \\               
                                                                                                      {r_{int} < {2\left( r_{G0} - r_{R0} \right)}} \\   
                                                                                                      \end{array}$                                       

Funding Information
===================

This paper was supported by the following grants:

-   Gordon and Betty Moore Foundation to Babak Momeni.

-   Life Science Research Foundation to Babak Momeni.

-   National Science Foundation to Kristen A Brileya.

-   United States Department of Energy DE-AC02-05CH11231 to Kristen A Brileya, Matthew W Fields.

-   National Institutes of Health DP2 OD006498-01 to Babak Momeni, Wenying Shou.

-   W. M. Keck Foundation to Wenying Shou.

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.

Additional information
======================

The authors have declared that no competing interests exist.

BM, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article.

KAB, Conception and design, Acquisition of data, Drafting or revising the article.

MWF, Conception and design, Drafting or revising the article.

WS, Conception and design, Analysis and interpretation of data, Drafting or revising the article.

Additional files
================

10.7554/eLife.00230.026

###### MATLAB code for simulating cooperative communities in the fitness model.

Other ecological interactions can be simulated by changing the values of the interaction parameters *r*~*12*~ and *r*~*21*~.

**DOI:** [http://dx.doi.org/10.7554/eLife.00230.026](10.7554/eLife.00230.026)

10.7554/eLife.00230.027

###### MATLAB code for simulating cooperative communities in the diffusion model.

In this version, the diffusion constant inside the community is assumed to be similar to the diffusion constant of agarose ([Figure 2---figure supplement 2A](#fig2s2){ref-type="fig"}).

**DOI:** [http://dx.doi.org/10.7554/eLife.00230.027](10.7554/eLife.00230.027)

10.7554/eLife.00230.028

###### MATLAB code for simulating competitive communities in the diffusion model.

In this version, the diffusion constant inside the community is assumed to be smaller than the diffusion constant of agarose ([Figure 2---figure supplement 2B](#fig2s2){ref-type="fig"}).

**DOI:** [http://dx.doi.org/10.7554/eLife.00230.028](10.7554/eLife.00230.028)

10.7554/eLife.00230.029

###### MATLAB function for in-plane rearrangement of cells in both the fitness and the diffusion models.

**DOI:** [http://dx.doi.org/10.7554/eLife.00230.029](10.7554/eLife.00230.029)

10.7554/eLife.00230.030

Decision Letter

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

Thank you for choosing to send your work entitled "Cooperation leads to partner intermixing in microbial communities" for consideration at *eLife*. Your article has been evaluated by a Senior Editor and 3 reviewers, one of whom is a member of our Board of Reviewing Editors.

The Reviewing Editor and the other reviewers discussed their comments before we reached this decision, and the Reviewing Editor has assembled the following comments based on the reviewers' reports.

The study deals with a very interesting mixture of experimental, theoretical, and simulation approaches to arrive at basic rules of spatial patterning between obligatory cooperative partners in single-celled systems. The studies are based on a previously developed system of two yeast strains that can enter cooperative interactions. The present paper deals with the development of a spatial patterning model in conjunction with a series of experiments that use different starting conditions to infer the stability of the observed patterns. The results of these experiments are then compared with the respective simulation results to validate the theoretical approach. This in turn is taken to infer further, more complex interactions in multi-partner systems. An additional experimental validation comes from the use of two bacterial species that interact on biochemically very different principles.

The referees agree that the paper is very interesting, technically sound, and the whole work rests on a very rich database. The main research question addressed is certainly significant, and the authors elegantly investigate it by using a combination of experimental and theoretical approaches. There are no severe methodological problems that need to be addressed. However, there are some concerns about the generality of the findings. The current version reads in places as if general rules were derived; but, strictly speaking, they are only tested for systems with obligate cooperation. Also, the possibility of antagonistic interactions is not addressed at all. Still, it appears that the paper is a very strong starting point for future work into extended directions. Hence, we will by happy to offer publication in *eLife*, provided the following points are addressed in a revised version.

A\) The authors should be very explicit about the limitations of their current study, as well as the strengths. This could be presented in an extended discussion.

B\) One of the reviewers suggests additional experiments in this context, but after further discussions, we consider them to be optional:

"The authors investigate a very special case of cooperation, namely an obligate interaction. Both their synthetic yeast system and the Desulfovibrio--Methanococcus interaction represent a case of an obligate interaction: that is, none of the two partners can grow without the respective other. In this light, the finding that cooperative interactions mix more than other ecological interactions seems much less surprising. It may simply be the consequence of the properties that were previously built into the system by the experimenter (i.e., the obligate mutual dependency) rather than a pattern that generally emerges in cooperative interactions (which are not necessarily obligate in nature). Support the interpretation that the partner intermixing observed by the authors is limited to obligatory cooperation.

Another very prominent example of cooperation in a microbial community is the cooperative production of biofilms (e.g., Rainey & Rainey 2003 *Nature*). Here we would also see spatial patterns in a biofilm, but I doubt that one would find the enhanced population mixing (given that two biofilm-formers were differentially labeled) as shown by the authors.

A way to manipulate the degree of mutual dependency experimentally would be to add lysine and/or adenine to the medium that is used to grow the yeast strains (as has been done to mimic the commensalism). My expectations for such an experiment would be that with more 'public goods' added, the less tight the interaction gets and the less pronounced the population mixing will be. Such an outcome would further support the interpretation that the conclusions drawn do not apply to cooperative interactions within microbial communities in general, but be rather limited to obligate cooperative interactions."

Please respond to this comment either by doing the experiment or by clearly discussing these points in the manuscript.

10.7554/eLife.00230.031

Author Response

The current version reads in places as if general rules were derived; but, strictly speaking, they are only tested for systems with obligate cooperation. Also, the possibility of antagonistic interactions is not addressed at all.

When we originally wrote "general rules," we meant general rules on patterning driven by the fitness effects of ecological interactions. We have modified "general rules" to "ecological patterning." We have also eliminated the usage of "rules." We tested facultative cooperation in the fitness model, and, unfortunately, experiments will require a significant body of work. Antagonistic interactions were tested in the fitness model. We have added a discussion about antagonistic interactions.

*A) The authors should be very explicit about the limitations of their current study, as well as the strengths. This could be presented in an extended discussion*.

We are not sure whether the "limitations" refer to the lack of experimental data for interactions, such as facultative cooperation, or the existence of other patterning forces that can cause deviations from the ecological patterning rules.

The former point is discussed below. We do not believe that the latter can be regarded as limitations of our study. There are different forces that shape patterns of communities, one of which is the fitness effects of ecological interactions between populations. We show that if such fitness effects are responsible for community patterns, strong cooperation is the only type of interaction that leads to population intermixing. To make an analogy, species coexistence is not considered a limitation of the competitive exclusion principle. There are other forces that make species not compete for exactly the same resources, and the competitive exclusion principle provides a framework of expectations that allows us to recognize and study those forces. We do agree that the original discussions on forces that cause deviations from the ecological patterning expectations could benefit from being more explicit and organized. Thus, we have expanded and reorganized discussion into three subsections.

*B) One of the reviewers suggests additional experiments in this context, but after further discussions, we consider them to be optional*:

*"The authors investigate a very special case of cooperation, namely an obligate interaction. Both their synthetic yeast system and the Desulfovibrio--Methanococcus interaction represent a case of an obligate interaction: that is, none of the two partners can grow without the respective other. In this light, the finding that cooperative interactions mix more than other ecological interactions seems much less surprising. It may simply be the consequence of the properties that were previously built into the system by the experimenter (i.e., the obligate mutual dependency) rather than a pattern that generally emerges in cooperative interactions (which are not necessarily obligate in nature). Seems to support the interpretation that the partner intermixing observed by the authors is limited to obligatory cooperation*.

The reviewer is correct in pointing out that all our experiments were performed on obligatory cooperative systems. We only considered facultative cooperation where we showed that the larger the fitness benefits are compared to the basal fitnesses of the two partners, the more quickly intermixing index increases as a function of community height. To make this point more explicit, we have added to the main text:

"Using the fitness model, we found that in facultative cooperation, small fitness benefits generated less intermixing than large fitness benefits. This result is intuitive: patterns of facultative cooperation should resemble those of obligatory cooperation (or competition) when benefits to both partners are large (or small) compared to the basal fitnesses of cooperative partners. 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."

In we showed that excessive benefits diminish the degree of intermixing. Excessive benefits are unrelated to whether cooperation is facultative or not. For example, if **A** completely depends on **B** for a metabolite, **B** can still supply this metabolite in excess.

*Another very prominent example of cooperation in a microbial community is the cooperative production of biofilms (e.g., Rainey & Rainey 2003 Nature). Here we would also see spatial patterns in a biofilm, but I doubt that one would find the enhanced population mixing (given that two biofilm-formers were differentially labeled) as shown by the authors*.

The cooperative production of biofilm (Rainey & Rainey, 2003) relied on cooperation among cells of the *same* population. This type of intra-population cooperation is different from the inter-population cooperation **A\[**↑ ↑**\]B** described here, in which the benefit produced by **A** is different from that produced by **B**. We have now made this point more explicit in the main text:

"..., while inter-population cooperation based on the exchange of distinct benefits would be **\[**↑ ↑**\]**."

*A way to manipulate the degree of mutual dependency experimentally would be to add lysine and/or adenine to the medium that is used to grow the yeast strains (as has been done to mimic the commensalism). My expectations for such an experiment would be that with more 'public goods' added, the less tight the interaction gets and the less pronounced the population mixing will be. Such an outcome would further support the interpretation that the conclusions drawn do not apply to cooperative interactions within microbial communities in general, but be rather limited to obligate cooperative interactions.*"

We have considered the experimental scheme proposed by the reviewer. However, we realized that if we supplement adenine and lysine to the agarose, the cells would simply exhaust the supplements first (competition) and then the interaction would switch to obligatory cooperation in addition to competition for other shared resources. One could argue that a steady supply of a very low concentration of adenine and lysine would create a facultative cooperative system, but there is no easy way of evenly supplying all cells in the community. In addition, any flow system would also wash out some of the released metabolites. To create *bona fide* facultative cooperation, we would need to either screen for natural microbes for facultative cooperation or create reduction-of-function (rather than loss-of-function) alleles in enzymes mediating adenine and lysine biosynthesis, and ideally these mutations should also be non-revertible. This would call for two genetic screens for appropriate mutations, which is beyond the scope of our study. Our fitness model predicts that in facultative cooperation, if for both partners the fitness benefits are large compared to the basal fitness, intermixing will increase rapidly as a function of community height.
