Spatial structure favors microbial coexistence except when slower mediator diffusion weakens interactions

Microbes often exist in spatially structured environments and many of their interactions are mediated through diffusible metabolites. How does such a context affect microbial coexistence? To address this question, we use a model in which the spatial distributions of species and diffusible interaction mediators are explicitly included. We simulate the enrichment process, examining how microbial species spatially reorganize and how eventually a subset of them coexist. In our model, we find that slower motility of cells promotes coexistence by allowing species to co-localize with their facilitators and avoid their inhibitors. We additionally find that a spatially structured environment is more influential when species mostly facilitate each other, rather than when they are mostly competing. More coexistence is observed when species produce many mediators and consume some (not many or few) mediators, and when overall consumption and production rates are balanced. Interestingly, coexistence appears to be disfavored when mediators are diffusing slowly because that leads to weaker interaction strengths. Overall, our results offer new insights into how production, consumption, motility, and diffusion intersect to determine microbial coexistence in a spatially structured environment.


23
Microbes are rarely found in isolation in nature. Instead, they are found coexisting with one 24 another in complex networks of interactions [1]. Given the differences among taxa and the 25 competitive forces that act between them, a fundamental question in microbial community 26 ecology is how this coexistence is maintained [2][3][4]. And because many important industrial, 27 environmental, and health-related processes rely on microbial communities to function (e.g. 28 anaerobic granules, microbial mats, and gut microbiota, respectively), understanding the 29 conditions that favor microbial coexistence is critical to sustaining these systems. 30 Spatial structure and organization may shape coexistence via numerous mechanisms [5][6][7][8][9][10][11], often 31 by modulating the interactions among individuals. For example, in a spatially structured 32 environment where progeny is more likely to be in the vicinity of parents, intensified 33 intrapopulation competition can give less competitive species a chance to survive [12]. In other 34 conditions, spatial isolation can allow organisms with conflicting abiotic needs to flourish in 35 appropriate environments [13,14]. The interplay between dispersal and competition can also 36 allow coexistence between species that are more competitive growers and species that are better 37 at dispersing and colonizing [15]. 38 Spatial heterogeneity has been invoked as a mechanism for microbial coexistence since the 39 pioneering work by Gause [16]. And although general concepts of coexistence are expected to 40 apply equally to microbes, microbial communities may be affected by spatial structure in unique 41 ways because of the scale and multiplicity of microbial interactions. An important and ubiquitous 42 example of this are microbial interactions that are mediated via diffusible metabolites-including 43 resources and metabolic byproducts. Spatial structure can stabilize these interactions and support 44 coexistence, for example, by allowing cheaters to be excluded from beneficial interactions 45 [17,18], or by permitting facilitative chemical interactions while preventing the inhibitory effects 46 of an interacting organism's physical presence [19]. And while it is clear that the outcomes of 47 interactions via diffusible mediators in structured environments may depend on mediator 48 diffusion rates [20,21] and the larger network of antagonistic and cooperative interactions [22], 49 how such factors translate into community-level consequences is not well understood. 50 Prior reports that address coexistence of metabolically interacting microbes in a spatially 51 structured environment are scarce. In an implicit model, Murrell and Law have shown in a 52 modified Lotka-Volterra model that when interspecific competition operates over shorter 53 distances than intraspecific competition a spatially structured environment can lead to species 54 coexistence by allowing for aggregation [23]. And in recent work with explicit modeling of 55 space, Weiner et al examined coexistence in territorial populations interacting through diffusible 56 mediators and found that metabolic tradeoffs allow for the coexistence of more species than the 57 number of nutrients [24].

58
Our model is distinct from previous work in that we allow overlap and dispersal of populations 59 through the shared space. Our motivation is to capture situations in which microbes can disperse 60 inside a matrix that defines the spatial structure. An example of this is the mucosal layer of the 61 digestive or respiratory tract, in which stratification is possible, yet the distribution of different 62 species populations can overlap. Another example is in yogurt or cheese, where spatial structure 63 exists, but populations are not territorial. We modify a previously developed mediator-explicit 64 model [25] to account for spatial structure and the dispersion of species in the same space. Here, 65 we limit our study to one-dimensional (1D) spatial structure as a starting point. We examine in 66 our model conditions under which coexistence is favored. We should emphasize that even 67 though we choose our parameters within the range of typical values observed among microbial 68 communities, the purpose of this work is not to recapture a specific community. Instead, by 69 examining a range of values for parameters such as metabolite diffusion and species dispersal, 70 we hope to gain a better understanding of how rates of these processes can affect species 71 coexistence.

73
A spatial mediator-explicit model of microbial communities 74 In our mediator-explicit model, species interact through metabolites that they produce and/or 75 consume (Fig 1A-B) [25]. Each species can produce a subset of metabolites and consume a 76 subset of metabolites. Each of the metabolites in the shared environment can in turn influence 77 any of the species by increasing or decreasing their growth rate (i.e. facilitation or inhibition, 78 respectively) compared to how each species grows in the absence of interactions [25]. We also 79 assume that different interaction mediators additively influence the overall growth rate of each 80 species (see "Model description" in Methods). 81 We assume a 1D spatial structure which preserves the spatial context but allows the diffusion of 82 metabolites and dispersal of species. Multiple metabolites and species can be present in a single 83 location. Both metabolite diffusion and species dispersal are modeled as random walk processes, 84 characterized with a diffusion coefficient and a dispersal coefficient, respectively. In a typical 85 simulation, we start from an initial distribution in which populations occupy adjacent, non-86 overlapping spatial locations at low initial density. This choice is made to impose a reproducible 87 initial condition that emphasizes the role of space. Each simulation starts with a network of 88 interactions in which interaction strengths, production and consumption links, and production 89 and consumption rates are assigned randomly. The initial pool typically contains 10 species and 90 5 interaction mediators. We simulate community enrichment through rounds of growth and 91 dilution [25,26] for 100 generations, and assess the richness of each resulting community (i.e., the 92 number of species stably persisting in the community). We have chosen 100 generations of 93 growth, because we have observed that often this is enough to reliably decide which species 94 stably persist in the community (Fig 1-figure supplement 1). At each dilution step, we assume 95 that the overall spatial distribution of the community is preserved and all populations at all 96 locations are diluted with the same factor. We recognize that this assumption is not universally 97 true; however, we adopt it as an approximation, in the absence of additional information about a 98 particular community. Such a dilution preserves some of the spatial structure of the community 99 in the next round of growth and could represent a biofilm getting partially washed away by rain 100 or in a microfluidic device, gut microbiota after a defecation event, or a broken-off portion of a 101 granule initiating a new granule. We use a well-mixed version [25]-devoid of any spatial 102 context-with the same set of parameters for species properties and interactions (i.e. dynamics during the course of enrichment. In a simple example, we show that interactions and 106 subsequently the population dynamics are affected by growing in a well-mixed versus spatial 107 environment (Fig 1-figure supplement 3). We explored the impact of the overall spatial extent of 108 the community and found that within an order of magnitude of change, the outcomes remained 109 the same (Fig 1-figure supplement 4).

110
The shift from interspecies competition to intraspecies competition can favor coexistence in a 111 spatially structured environment. To assess this impact, we imposed a cap on total cell number at 112 each location in space. As this cap became more restrictive, it suppressed the most competitive 113 species and led to higher coexistence (Fig 1-figure supplement 5). Since our focus in this 114 manuscript is the impact of interspecies interactions, in the rest of this manuscript we pick the 115 total cell number cap at a level (k Y = 10 9 cells/ml) that minimizes the impact of imposed 116 intrapopulation competition.

117
A spatial environment favors coexistence more when facilitation among species 118 is prevalent 119 We first examined how the prevalence of facilitative versus inhibitory interactions impacted 120 coexistence in spatial communities. In our simulations, we dictated the ratio of facilitative and 121 inhibitory interactions in the initial pool of species. Our results show that, similar to a well-122 mixed environment, more facilitative interactions lead to higher richness in communities that 123 emerge from enrichment ( Fig 1C, along the x-axis). Additionally, we observe that spatial 124 communities show more coexistence than well-mixed communities when facilitation among 125 species is prevalent (Fig 1C, spatial versus well-mixed). The same pattern, although less 126 pronounced, was present when instead of richness we used the Shannon index to assess the 127 diversity of resulting communities (Fig 1-figure supplement 2). Our explanation is that species interactions and dampens inhibitory interactions, leading to more coexistence. This is supported 149 by our data which shows that the position of specific species with respect to other species that 150 facilitate or inhibit it can impact the population dynamics (Fig 1-figure supplement 8). Because 155 At low species dispersal, self-organization is one of the mechanisms that can lead to a difference 156 between spatial and well-mixed communities (Fig 1-figure supplement 3). In a simplified 157 interpretation, self-organization can be in the form of co-localization driven by facilitation or 158 segregation driven by inhibition. In our simulations, we observed that co-localization had a 159 stronger effect on coexistence. The positive influence was reinforced by more growth in the 160 vicinity of the partner, leading to a stronger representation of facilitation in spatial communities.

161
In contrast, segregation only had a modest effect on weakening the impact of inhibition. As a 162 result, there is more similarity between well-mixed and spatial communities in the absence of 163 strong facilitative interactions ( Fig 1C). when each mediator influenced too many or too few species on average. We note that these 174 trends were largely the same between spatial and well-mixed communities.

175
Our explanation is that a larger range of production offers more opportunities for interaction, 176 which through the enrichment process lead to the selection of facilitative subsets that coexist 177 [25]. A low mediator influence range works in the opposite direction, reduces opportunities for 178 interactions and results in lower coexistence. Very high mediator influence range potentially 179 leads to more self-facilitation (i.e. producing a metabolite that is beneficial to the producer 180 species), which our data suggests can lead to take-over by a single species and a lower 181 coexistence as a result (Fig 2-figure supplement 1).

182
Coexistence is higher when there is balance between production and 183 consumption of mediators 184 We next asked how the rates of production and consumption of mediators would influence 185 coexistence. To address this question, we surveyed a range of average rates of production and 186 consumption. We observed that the highest levels of coexistence occurred when there was a species stably present at the end of a simulation) was calculated for 500 simulated instances and marked 193 on the color bar. Each simulation started with 10 species and 5 mediators and ran for 100 generations. The 194 x-axis represents the average number of species influenced by a mediator and the y-axis represents the 195 average number of mediators produced by each species. Other simulation parameters are listed in Table 1. 196 197 balance between consumption and production rates among species, with slightly higher 198 production than consumption (Fig 3).

199
Our justification for the observed pattern is that in one extreme where production is too high 200 (lower right corner of each plot), mediators will build up in the environment. This will put the 201 community in a regime in which consumption is not enough to create a feedback, i.e. "reusable 202 mediators" as discussed in [25], which leads to lower coexistence. In the other extreme, when 203 consumption is too high (upper left corner of each plot), metabolites that mediate the interactions 204 will be depleted from the environment, leading to an effectively weaker interaction and thus 205 lower coexistence. However, when production is slightly higher than consumption, metabolite 206 quantities are sufficient to create strong interactions and facilitation feedbacks, leading to higher 207 coexistence. While coexistence is slightly higher in the spatial communities compared to well-208 mixed ones, the production-consumption trends apply equally to spatial and well-mixed 209 communities, as expected. simulation) is calculated for 500 simulated instances. Each simulation started with 10 species and 5 216 mediators and ran for 100 generations. Color bar represents the average richness. The x-axis represents 217 the average production rate of mediators and the y-axis represents the average consumption rate of 218 mediators. Other simulation parameters are listed in Table 1. 219 Limited species dispersal in a spatial environment allows more coexistence, 220 especially when facilitation is common 221 Because species dispersal is a major factor in preserving community spatial structure, we 222 examined how the dispersal coefficient affected coexistence outcomes. For this, we kept the 223 diffusion coefficient of the mediators fixed and surveyed mean richness among many instances 224 of communities randomly assembled (n = 500). When the diffusion coefficient for species 225 approaches zero and cells remain in their original spatial location, we observe higher levels of 226 coexistence (Fig 4). We also observed that the impact of lower dispersal is stronger in 227 communities in which most interactions are facilitative rather than inhibitory. Our explanation is 228 that lower dispersal rates mean that species grow best in spatial locations that are more 229 supportive for their growth, which is in the vicinity of their beneficial partners and away from 230 competitors or inhibitors. As discussed in Fig 1, such self-organization effectively amplifies 231 facilitative interactions and de-emphasizes inhibitory interactions, leading to a higher 232 coexistence. This is also consistent with the observation that the effect of dispersal rate is 233 strongest when the proportion of facilitative interactions is highest. As the dispersal coefficient 234 increases, the self-organization gets washed away by dispersal of cells to less than ideal locations 235 for their growth and its benefit for coexistence diminishes.

236
At intermediate levels of dispersal, the trend reversed and well-mixed communities showed more 237 coexistence compared to spatial communities. This is interesting because at the limit of 238 extremely rapid diffusion (shown with a '∞' sign in Fig 4) when we kept the species distribution 239 uniform across the spatial extent, coexistence outcomes matched the well-mixed case, as 240 expected. We found two factors that contributed to this trend. The first contribution came from and 5 mediators and ran for 100 generations. Other simulation parameters are listed in Table 1. The error  264 bars are 95% confidence intervals generated by bootstrapping 100 samples. 265 species will be more confined in space and some of the metabolite will leak out to other species; 266 in the other extreme, in the very high dispersal rates all distributions will become uniform and 267 the distinction between self and others diminishes. To test this, we tested weaker self-facilitation 268 links in our simulations and observed that this change led to higher coexistence in communities 269 with intermediate dispersal coefficients but not in well-mixed communities or communities with 270 low dispersal coefficients (Fig 4-figure supplement 2). It is a matter of debate how prevalent self-  Coexistence is disrupted when the diffusion of mediators is too slow 276 The rate of diffusion of metabolites also has the potential to affect coexistence. We investigated 277 coexistence over a range of mediator diffusion coefficients. We still typically observe a higher 278 mean richness for spatial communities compared with the well-mixed communities (Fig 5).  Table 1. The error bars are 95% confidence intervals generated by bootstrapping 100 samples. 293 consumed by nearby species and do not travel long enough to reach other members of the 294 community, the interaction-driven mechanism of coexistence is disrupted.

296
Our results dispel the common presumption that a spatially structured environment will 297 universally lead to more coexistence. We find that, compared to a well-mixed environment, a 298 spatial environment can favor or disfavor coexistence depending on the balance between species 299 dispersal and the diffusion of interaction mediators. Interestingly, a lower species dispersal rate 300 favors coexistence, but this effect can be diminished or even reversed if accompanied by low 301 mediator diffusion rates. Coexistence is favored when species have a broad range of 302 consumption and an intermediate range of production of interaction mediators. Additionally, we 303 predict more coexistence when there is a balance between overall production and consumption 304 rates for mediators.

305
The spatial structure of microbial communities has been extensively studied for example in 306 simulating the development of biofilms [27][28][29][30], for specific interactions among species 307 [17,31,32], or for modeling game-theory dynamics [33][34][35][36]. However, as there is often a tradeoff 308 between the incorporation of detailed mechanisms and generality of conclusions [37,38], we 309 chose in this work to explore a simple, general model of chemically-mediated microbial 310 interactions. We assumed, for example, that mediators affected species by additively influencing 311 their growth rates. Although it is possible (and even probable) that mediator effects could be 312 multiplicative, nonlinear, or otherwise context-dependent and that they may impact other model 313 parameters, we chose here to present what we felt to be the simplest case. Exploration of 314 alternative implementations of mediator effects would make a fascinating follow-up to this work. 315 We have made assumptions in our model to simplify the configuration and make the analyses 316 and interpretations easier. We asked if making these assumptions more realistic would affect our 317 conclusions. For example, we have assumed no carrying capacity limit for the growth of our 318 populations. We explored the effect of imposing a total population limit, enforced at each spatial 319 location, and found that it did alter our conclusions (Fig 1-figure supplement 5). However, 320 because the relationship between carrying capacity and coexistence has been explored 321 extensively elsewhere, we chose parameters to minimize this impact, allowing us to focus on 322 other interspecies interactions (beyond competition) and relative rates of diffusion and dispersal. 323 We also tested the impact of the spatial extent of the community (Z), and observed that our 324 results were largely unaffected if the community's spatial extent was changed by an order of 325 magnitude (Fig 1-figure supplement 4). The effect of larger changes in the spatial extent can be 326 examined by scaling the diffusion and dispersal coefficients accordingly.  If spatial organization of cells matters, we also expect that the initial spatial position of species in 338 the community impacts coexistence. To test this, we started from 100 simulations instances and 339 in each case, we tried 100 rearrangements, each obtained by shuffling the spatial position of 340 species, while keeping the species properties and interactions intact. Interestingly, in many cases 341 coexistence was affected (Fig 4-figure supplement 3), indicating that the adjacency to partners is 342 an important determinant of spatial coexistence (as also suggested by Fig 1-figure supplement 3   343 and Fig 1-figure supplement 7). When we examined the effect of the fac:inh ratio on these 344 outcomes, we observed that larger changes in richness when facilitation interactions were more 345 prevalent in the community (Fig 4-figure supplement 4), which aligns with many of our other 346 results showing that facilitation amplifies the positive effect of spatial structure on coexistence.

347
Although these results are tantalizing, a detailed examination of the spatial organization of  Overall, we believe this work revisits how spatial structure-and spatial self-organization-358 affects community assembly and coexistence. In our model, which emphasizes the contributions 359 of interspecies interactions, we find that the impact of spatial structure on coexistence largely 360 arises from two processes: (1) spatial self-organization, which can improve coexistence by 361 favoring facilitation over inhibition, and (2) localization of interactions, which can promote 362 coexistence in association with self-organization or hamper coexistence by slowing down and 363 weakening species interactions.

365
Model description 366 Our model is an extension of a model introduced earlier [25] in which a set of species interact 367 with each other through diffusible mediators. Each mediator is produced by a subset of species, 368 consumed by a subset of species, and has a positive or negative influence on the growth rate of 369 some species (Fig 1).

409
To assess coexistence, we use a criterion similar to [25]. In short, any species whose density 410 drops below a pre-specified extinction threshold (ExtTh) is considered extinct. Among species 411 that persist throughout the simulation, only those are considered to coexist whose relative 412 frequency does not drop by more than 10% in the last 20 generations of the simulation. We 413 consider these species to be 'stably present' in the community. Species whose relative frequency 414 declines faster than this threshold are assumed to go extinct later and are not considered to be 415 part of coexisting communities. The only exception to this criterion is the data in Fig 4-  of growth (one round of growth between dilutions). There was minimal change (less than 1% difference) 547 after 100 generations of growth in both well-mixed and spatial contexts. This observation was consistent 548 regardless of the initial ratio of facilitation to inhibition interactions, at fac:inh = 10:90 (top), 50:50 549 (middle), and 90:10 (bottom). In all cases, D Cell = 5×10 -9 cm 2 /hr and D Med = 1.8×10 -2 cm 2 /hr. fac:inh ratios for both spatial communities (blue squares) and well-mixed communities (red circles). Each 555 ratio was run 500 times with the richness (number of species stably surviving at the end of a simulation) 556 averaged over all the simulations. Each simulation started with 10 species and 5 mediators and ran for 100 557 generations. The error bars are 95% confidence intervals generated by bootstrapping 100 samples. Here, 558 the species dispersal coefficient is × cm 2 /hr. 559 560 Fig 1-figure supplement 3. Comparing the spatial distribution of species at different dispersal rates 562 illustrates the impact of dispersal on coexistence. (A) Species are initially stacked over the spatial 563 extent of the community. We examine a simple case with 5 species and 3 mediators (interaction network 564 shown in the inset). Different progressions are observed at high (B, D Cell = 5×10 -6 cm 2 /hr) versus low (C, 565 D Cell = 5×10 -9 cm 2 /hr) dispersal rates, leading to different coexistence outcome. In both cases, D Med = 566 1.8×10 -2 cm 2 /hr. 567 569 Fig 1-figure supplement 4. Species interactions and dynamics are different in spatial versus well-570 mixed environments, leading to different coexistence outcomes. For simplicity, we consider an initial 571 pool of 5 species with 3 mediators. We kept the same interaction network (A) in both spatial and well-572 mixed conditions and simulated the species dynamics. In the interaction network, hollow arrows indicate 573 the flow of mediators (either production or consumption). The thickness of the lower end of each arrow 574 shows the production rate and the thickness of the upper end of the arrow shows the strength of the 575 mediator influence on the corresponding species. The basal growth rate of each species is shown as 576 different shades of gray, with darker shades indicating larger basal growth rates. For spatial communities, 577 the population density represents the sum of all densities across different spatial extents. We ran the 578 simulation in both cases through 10 rounds of growth and dilution. Different population dynamics and 579 ultimately different coexistence outcomes are observed in these cases after simulating 100 generations. 580 Note that even though species 4 and 5 in the spatial community and species 1, 4, and 5 in the well-mixed 581 community are still present in the final communities, because they exhibit a trend of decline in relative 582 frequency, we consider them not to coexist (in accordance with our coexistence criterion, see Methods).

583
(B-C) In this example, proximity of Species 1 to a beneficial Species 2 in the spatial community allows a 584 strong boost in the growth of Species 1, leading to its coexistence with Species 2 and 3. In contrast, 585 Species 1 in the well-mixed community receives a smaller portion of C 1 , not enough to allow Species 1 to 586 keep up with the other species. In C, only coexisting species and mediators related to them are shown. 587 mixed environments, even when inhibition is prevalent. For simplicity, we consider an initial pool of 5 590 species with 3 mediators, but rather than fac:inh = 90:10 in Fig 1-figure supplement 4, here we initially 591 assign fac:inh = 10:90. We kept the same interaction network (A) in both spatial and well-mixed 592 conditions and simulated the species dynamics. In the interaction network, hollow arrows indicate the 593 flow of mediators (either production or consumption). The thickness of the lower end of each arrow 594 shows the production rate and the thickness of the upper end of the arrow shows the strength of the 595 mediator influence on the corresponding species. The basal growth rate of each species is shown as 596 different shades of gray, with darker shades indicating larger basal growth rates. For spatial communities, 597 the population density represents the sum of all densities across different spatial extents. We ran the 598 simulation in both cases through 10 rounds of growth and dilution. Different population dynamics are 599 observed in these cases. Note that even though species 4 in the spatial community is still present in the 600 final communities, because they exhibit a trend of decline in relative frequency, we consider them not to 601 coexist (in accordance with our coexistence criterion, see Methods). (B-C) In this example, proximity of 602 Species 1 to an inhibitory Species 2 in the spatial community leads to rapid exclusion of Species 1. Early 603 extinction of species 1 removes its inhibition of Species 4 and allows Species 4 to sustain longer in the 604 spatial context compared to the well-mixed community. In C, only coexisting species and mediators 605 related to them are shown. 606 607 Fig 1-figure supplement 6. Spatial distance between species can modulate the strength of their 608 interaction. For simplicity, here we consider only two species engaged in commensalism (Species 1 609 providing a benefit to Species 2) through a single mediator (see insets). We examined a no-interaction 610 control (A), versus when the species were initially close (B) or far (C) from each other within the spatial 611 extent of the communities. All parameters are similar to Table 1, except the mediator diffusion coefficient 612 which is chosen at D Med = 1.8×10 -3 cm 2 /hr to exaggerate the impact of mediator diffusion in space. The 613 panel on the left shows the initial distribution of the two populations in space and the panel on the right 614 shows the population dynamics. Notably, the ratio of Species 1 to Species 2 drops by 10-fold, when 615 Species 2 is farther away from Species 1 (C versus B). The basal growth rate of Species 1 and 2 are 0.1/hr 616 and 0.08/hr, respectively. We ran the simulation for 5 rounds of growth and dilution. 617 does not have a large impact on spatial coexistence. All parameters other than the community's spatial 620 extent are kept fixed. Here (A) D Cell = 5×10 -8 cm 2 /hr and (B) D Cell = 5×10 -9 cm 2 /hr. Legend shows the 621 values of fac:inh assigned in the initial pool. In these simulations, the spatial resolution dz is kept at 0.05 622 mm in all cases and the initial extent of each species is 1/10 th of the total spatial extent at the beginning of 623 each simulation. 624 625 Fig 1-figure supplement 8. Imposing a local carrying capacity favors species coexistence. All 627 parameters other than the local carrying capacity k Y are kept fixed. Imposing a local carrying capacity on 628 the total cell number (see the equations described in Methods) lowers the competitive advantage of the 629 most successful species and allows more coexistence. This is similar to the idea of strengthened 630 intrapopulation competition in a spatially structured environment that has been discussed in depth in the 631 past as one of the reasons that spatial environments can support more coexistence. To focus on the effect 632 of interspecies interaction, we assume k Y = 10 9 cells/ml in our other simulations throughout this 633 manuscript to minimize the impact of k Y (i.e. a forced intrapopulation competition) on our results. 634