Abstract
The ubiquitous existence of microbial communities marks the importance of understanding how species interact within the community to coexist and their spatial organization. We study a two-species mutualistic cross-feeding model through a stochastic cellular automaton on a square lattice using kinetic Monte Carlo simulation. Our model encapsulates the essential dynamic processes such as cell growth, and nutrient excretion, diffusion and uptake. Focusing on the interplay among nutrient diffusion and individual cell division, we discover three general classes of colony morphology: co-existing sectors, co-existing spirals, and engulfment. When the cross-feeding nutrient is widely available, either through high excretion or fast diffusion, a stable circular colony with alternating species sector emerges. When the consumer cells rely on being spatially close to the producers, we observe a stable spiral. We also see one species being engulfed by the other when species interfaces merge due to stochastic fluctuation. By tuning the diffusion rate and the growth rate, we are able to gain quantitative insights into the structures of the sectors and the spirals.
Citation: Lin J, Sun H, Dong J (2022) Emergence of sector and spiral patterns from a two-species mutualistic cross-feeding model. PLoS ONE 17(10): e0276268. https://doi.org/10.1371/journal.pone.0276268
Editor: Attila Csikász-Nagy, Pázmány Péter Catholic University: Pazmany Peter Katolikus Egyetem, HUNGARY
Received: May 17, 2022; Accepted: October 4, 2022; Published: October 19, 2022
Copyright: © 2022 Lin et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper and its Supporting information files.
Funding: JJD, JL: DMR-1702321, MCB-2029480 (National Science Foundation, US); HS: MCB-2029580 (National Science Foundation, US). https://nsf.gov/ The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
1 Introduction
Understanding the emergence of complex patterns from simple dynamical rules and unraveling the rich behaviors is a continuous pursuit in physics. Microbial communities exist ubiquitously in nature, where multi-species occupy spatial spaces and interact with one another. Within microbial communities, how they survive, expand, and maintain species diversity provides much fodder for theoretical and experimental explorations. To make robust progress towards understanding a complex colony, it is important to identify the essential components within the colony [1], the spatial structure of the colony [2], and the stability of the colony diversity [3] with various interactions within species and between species and the environment.
One way through which species interact with one another is by exchanging metabolites, namely metabolite cross-feeding. In a simple two-species system, there are several possible cross-feeding mechanisms depending on whether one or both species produce cross-feeding metabolite, and whether one or both species benefit from the process. Among all possible species interactions, commensalism, syntrophy and mutualism [4–7] are often studied. In a system of two species, commensalism is defined where one species, the “producer”, produces metabolic byproduct that is essential for the growth of the other species, the “consumer”. The metabolite has a neutral effect on the growth of the producer. A common example of commensal interaction is that typically there are billions of the bacterium Staphylococcus epidermidis feeding on dead human skin cells while we, the “producer”, do not react to them [8]. Syntrophy is similar to commensalism in that the producer excretes metabolites that benefit the consumer. However, the accumulation of the metabolite impacts the producer negatively by, for instance, changing the environment pH or access to oxygen. Thus the presence of the consumer helps reduce the environmental stress and thus the survival of the producer. Syntrophic interactions are closely associated with microbes under anoxic conditions and energy constraints [5]. Mutualism is when each species excretes a distinct metabolite that is essential for the other species. The accumulation of each metabolite can impact the producer neutrally or negatively, further underscoring the significance of the presence of consumer. Since population range expansion tends to create local homogeneity due to the “founder effect” [9] while mutualistic interactions promotes species intermixing (see e.g. [10]), it is of particular interest to us to investigate the interplay between cross-feeding and spatial structure of the colony.
In this study, we focus on a two-species symmetric mutualistic cross-feeding model in two-dimensional space, where one species excrete metabolic nutrients, to be taken up by the other species, and vice versa. While both species divide and grow, they also compete for space to expand. Therefore, the colony spatial pattern depends crucially on the availability of the cross-feeding metabolites, which diffuse in open environments, and the local cross-feeding dynamics. We focus on the condition under which we can find stable coexistence of both species and the corresponding colonial spatial structure. The stability of microbial communities [3] has many important implications, including the production of pharmaceutical products [11], indication of environmental condition [12], reflection on climate changes, and so on. We expect that through a simple growth model, we can begin to unravel the complex behaviors of ecological systems.
Taking a similar approach as discrete cellular automata systems [13–15], we devise a stochastic cellular automaton on a square lattice as a minimal model to investigate a symmetric mutualistic cross-feeding between two species. The dynamical evolution of the system is simulated using the kinetic Monte Carlo method (KMC) [16–18]. This method is effective in describing variety of phenomena [19], including reaction-diffusion systems [20, 21], structures and properties of materials [22], and equilibrium and non-equilibrium chemistry [23]. We use KMC to analyze a range of parameters that give rise to distinct colony morphology.
In Section 2, we provide details of our model and simulation algorithm. We characterize the distinct colony patterns in Section 3 and analyze the interplay between cell growth and cross-feeding. Finally, we discuss in Section 4 the generality of our findings. We suggest that the approach presented in the article provides preliminary predictions and powerful guidance on further studies, both in theory and in experiments, on multi-species colony pattern formation.
2 Model specification
In this study, we focus on a symmetric set of mutualistic cross-feeding mechanism: Either species (type 1 or 2) has an individual growth rate λ1,2 which depends on the availability of the metabolite, molecule B (or A), produced by the other species and some other generic nutrient that is abundant in the environment. The metabolite excretion rates are γA,B and the excreted molecule, if not taken up, diffuses with a diffusion coefficient of DA,B. The yield YA,B is a parameter that converts the nutrient molecules taken up by a species into cell biomass, typically measured as mass of nutrient per dry mass of cell or concentration of nutrient per unit optical density [24]. The dynamics of the system can be described in the following set of equations:(1)(2)(3)
Note that the cell concentration ρ1,2 and the nutrient concentration nA,B depend on both space and time. The parameters and reactions of our model are summarized in a schematic shown in Fig 1.
[Figure omitted. See PDF.]
As species 1 cells metabolize nutrient molecules A (produced by species 2) at a growth rate λ1, they excrete molecules B with rate γB as a metabolic by-product. Molecule B is the essential nutrient for species 2 cells’ growth and is taken up with yield YB. The same for the growth kinetics for cells of species 2.
The overarching question in this study is how spatiotemporal dynamics governed by cell growth and cross-feeding gives rise to distinct colony morphology and species diversity. We construct our simulation model on a 2-dimensional (2D) square lattice and incorporate cell doubling (growth), nutrient excretion and uptake, and nutrient diffusion. Compared to the nutrient molecules, the cell diffusion is much slower and is not included in this model. In this system, cells of both species are identical in size, occupying a single lattice site and cannot overlap. Each lattice size is the same length of the cell. Nutrient molecules have negligible spatial content and thus infinitely many of them can occupy the same lattice site. When a cell at lattice site is chosen to divide, the individual growth rate λ1,2 is determined by nA,B, the total nutrient concentration at and its four nearest-neighbor (n.n.) lattice sites:(4)
Here indicates the maximal single cell growth rate. The instantaneous growth rate is modified by a Monod factor KA,B [25] to reflect the nutrient dependence, as shown in Eq 4. The daughter cell is placed in one of the four n.n. lattice sites, chosen at random provided that it is empty. This is similar to the Eden model in surface growth [26, 27]. If a nutrient excretion process is chosen, then a randomly chosen producer cell produces a nutrient molecule which is placed at one of the n.n. sites of the producer cell. If an uptake process is chosen to occur, a cell is picked at random. We compute the nutrient concentration in the 4 n.n. sites and the cell picks up one nutrient molecule chosen at random. When nutrient molecules diffuse, each diffusion step puts them in one of the four n.n. sites as well. In our simulation scheme, cells that do not divide due to spatial limitation still consumes nutrients and excrete metabolic by-product in an attempt to mimic a basal level of metabolic maintenance. This, however, can be revised into other scenarios to capture different metabolic processes which we will not delve into here.
The system is typically initialized with an inoculum patch of L × L lattice sites with ρ0L2 cells, half of each species randomly distributed within the patch. L2 is much smaller than the final colony size and ρ0 is 1/4 unless otherwise specified. We implement the kinetic Monte Carlo algorithm [17, 18] using the direct method to simulate the dynamics of the system outlined in Fig 1. The individual rates for each process are summarized in Table 1. Specifically, ρ1,2 is the total number of cell type 1 (or 2), while nA,B is the total number of nutrient molecule A (or B) in the entire system at a given time t. The propensity kj of each reaction, which scales with reaction probability, is determined by the respective reaction rate (time-independent) and the amount of reagents in the system at time t, and k0 is the sum of all propensities.
[Figure omitted. See PDF.]
Let r1, r2 ∈ (0, 1) be two independent, uniformly distributed random numbers. At time t, the j-th reaction (j ∈ [1, 8]) takes place if:Afterwards t is advanced by Δt where
The iteration continues until the system reaches a final time Tmax, typically after an obvious pattern is established. The system is large enough so that the colony can keep growing without running into the boundaries. When we refer to growth rate in the remainder of the study, we almost exclusively mean the maximal single cell growth rate. For ease of notation we will therefore drop the asterisk and denote it as λ1,2. Throughout this study, the parameters are symmetric for both species, namely they will have the same individual growth rate, nutrient excretion and uptake rates, as well as the nutrient diffusion rate. Each species produce one type of nutrient and are simultaneously consuming the nutrient produced by the other species.
3 Results
3.1 Sectors with high nutrient diffusion
When the essential nutrients for both species’ growth are readily available, it has been shown that an initially mixed populations of two species will segregate into sector-like domains as a result of random fluctuations at the expanding colony frontier experimentally [28, 29] and theoretically [30–32]. The domain boundaries perform super-diffusive random walks with respect to the radial growth. In the presence of mutualistic cross-feeding, when sufficient nutrients are produced by the cells, the system should behave in the same way as the nutrient-rich scenario in the absence of cross-feeding. We set out to first verify our model in this limit by comparing the following cases: 1) the environment provides nutrients for both species without cross-feeding; 2) cross-feeding with high nutrient excretion rates γA,B and high diffusion coefficients DA,B.
In the first case, the 2D system is initially seeded with an equal amount of species 1 and 2 cells, 50 of each type randomly distributed in a 20 × 20 lattice. Every lattice site, namely the “environment”, starts with 15 of A and 15 of B molecules. Since nutrients are widely available, both species grow without relying on cross-feeding. However, they are subject to the spatial constraint and only cells adjacent to at least one empty lattice site can divide. Fig 2 shows the population of each species ρ1,2(t) over time. In the early stage (t ≲ 20), the colony quickly expands with little spatial constraint. It then transitions to a constant radial expansion. In our system, only cells at the colony frontier have access to empty space to grow and thus contributing to colony expansion. Assuming an approximate circular colony of radius r, we have the total colony size ρ ≈ πr2 and the number of cells at the expanding frontier 2πr. The radial expansion of the colony is approximately:(5)This gives a constant colony expansion speed that depends on the individual growth rate λ and the total population within the colony grows quadratic with time.
[Figure omitted. See PDF.]
A) Colony population ρ1,2(t) and B) colony radius r(t). Equal number of species 1 and 2 are seeded in a 20 × 20 patch with initial density ρ0 = 1/4. λ1,2 = 1, DA,B = 1, YA,B = 1, with nutrients placed throughout the lattice.
In Fig 2, we show the population ρ1,2(t) and the colony radius r(t) for a system with nutrients available from the environment. In our simulation, we determine the colony radius r(t) by measuring the distances between the center and all cells adjacent to an empty lattice site, and then taking the average among all the distances. The colony expands radially with both cell populations grow at t2.11. The expansion speed dr/dt, measured to be 1.17, is consistent with the individual growth rate λ1,2 = 1.
The two species within the colony are separated by meandering interfaces. A colony snapshot at t = 600 with interface marked in red is shown in Fig 3. As expected, the initially mixed two populations segregate into sectors as the individual cells divide and the colony expands. Because the fluctuation of the interfaces is intimately related to the colony morphology, we analyze the inter-species interfaces to see whether the fluctuations resemble a standard random walk. We adopt a similar quantification method as discussed in Ref. [28] with details included in S1 File.
[Figure omitted. See PDF.]
Colony snapshot at t = 600. Equal number of species 1 and 2 are seeded in a 20 × 20 patch with initial density ρ0 = 1/4. λ1,2 = 1, DA,B = 1, YA,B = 1 with nutrients placed throughout the lattice. Scale bar indicates the width of 100 cells.
For an interface that fluctuates like a standard random walk, we expect the mean-square displacement to scale linearly with time, . In our case, as shown in Fig 4 indicates a super-diffusive behavior which is consistent with the findings in Ref. [28]. In this scenario, the adjacent interfaces extends linearly with time while their fluctuation grows sub-linearly. This establishes a long-time stability of the sector morphology in 2-dimensional space.
[Figure omitted. See PDF.]
Interface mean square displacement as a function of the sliding-window size x with nutrients from the environment. Symbols (blue, color online): 4 individual interface trajectories in Fig 3; Black line: average of the 4 trajectories. Red line: power-law fit with exponent 1.33 ± 0.18. Equal number of species 1 and 2 are seeded in a 20 × 20 patch with initial density ρ0 = 1/4. λ1,2 = 1, DA,B = 1, YA,B = 1 with nutrients placed throughout the lattice.
Now we turn to the second case to see how nutrient excretion rate γA,B under high nutrient diffusion DA,B impacts the colony morphology. Intuitively, when cross-feeding produces abundant nutrients through high excretion and/or high yield, while fast diffusion brings nutrients to the growing frontier, the colony dynamics should resemble the first case where nutrients are readily available from the environment and the two species do not depend on each other spatially.
When we set both excretion rates and diffusion to be high, we indeed see the colony pattern develops similarly as the case where nutrients are pre-seeded in the environment. In Fig 5, we show a colony snapshot at t = 140 with cross-feeding at a high nutrient excretion rate γA,B = 10 and high nutrient diffusion rate DA,B = 500. Similar to Fig 3, the colony again separates into sectors of same-type species. The interface mean-square displacement , same as in Fig 4. Fig 6 shows the cell population ρ1,2(t) and colony radius r(t). The overall population displays a quadratic growth while r2 increases faster than ρ. This is a result of deviation from circular colony growth: We measure the colony circumference roughness as defined in [27] for the two growth conditions shown in Figs 3 and 5, and indeed find an increase in roughness in the latter case.
[Figure omitted. See PDF.]
A) Colony snapshot at t = 140 for cross-feeding with high nutrient excretion. Inoculum condition is the same as in Fig 3. λ1,2 = 1, DA,B = 500, γA,B = 10, YA,B = 1. Scale bar indicates the width of 100 cells. B) Interface mean square displacement as a function of the sliding-window size x with high nutrient excretion. Symbols (blue, color online): 4 individual interface trajectories. Black line: average of the 4 individual ones. Red line: power-law fit with exponent 1.33 ± 0.40.
[Figure omitted. See PDF.]
A) Colony population ρ1,2(t) and B) colony radius r(t). Equal number of species 1 and 2 are seeded in a 20 × 20 patch with initial density ρ0 = 1/4. λ1,2 = 1, DA,B = 500, γA,B = 10, YA,B = 1.
The species’ reliance on cross-feeding develops as we dial down the nutrient excretion rate while holding nutrient diffusion constant. In Fig 7, we keep DA,B = 500 as in Fig 5 while reducing γA,B to 1. The nutrient concentration at the expanding front is no longer keeping up with diffusion in the same radial direction. As a result, the interface fluctuation increases over time, . This means at least some of the interfaces will collide and annihilate in a two-dimensional system, as seen in Fig 7 at later times.
[Figure omitted. See PDF.]
A)- C) Colony snapshots at t = 230, 300, and 420. D) Interface mean square displacement as a function of the sliding-window size x. Symbols (blue, color online): 6 individual interface trajectories. Black line: average of the 4 individual ones. Red line: power-law fit with exponent 2.51. Inoculum condition is the same as in Fig 3. λ1,2 = 1, DA,B = 500, γA,B = 1, YA,B = 1. Scale bar indicates the width of 100 cells.
The coalesce of interfaces clearly points to the possibility of one species being crowded out, or “engulfed”, by the other when all of the interfaces eventually merge. The specific time at which engulfment occurs will depend on the dynamics of the system. In this case, the colony expands radially with constant speed. However, the mean square displacement of the interface grows faster than t2, shown in Fig 7(D). This means that as time increases, the two interfaces will eventually coalesce.
3.2 Spirals emerge when cross-feeding limits growth
In the previous section, we discuss the emergence of sector patterns under high nutrient availability and the merging of sectors when nutrient concentration in the radial direction drops below the consumption. The stochasticity at the growing tip of the interface gives rise to the sector patterns of the colony. Here we hone in on the case where globally there is a scarcity of nutrients due to low excretion and low diffusion, thus the proliferation of either species depends on the availability of the cross-feeding nutrients.
Starting the system with a square patch of side-length 20 and initial density ρ0 = 1/4 as before, we set the excretion rate γA,B and the diffusion rate DA,B to be 1. Here with γA,B = YA,B = 1. This means each producer produces one molecule and the consumer needs one molecule to be at half-maximal growth rate.
We are surprised to see a stable spiral pattern of the colony emerge after a transient period. In Fig 8, we show the nascent colony patterns at t = 100, 150, 200 and 250. In the early stage of the colony growth (t ≲ 100), the small domains of the two species quickly merge due to the short separation distance and the resulting interfaces start to bend from the radial direction. Since cell type 1 needs nA (produced by cell type 2) to grow, only the ones near a type 2 cell are growing when diffusion is slow. The same is true for cell type 2. Both types of cells grow along the interface between them, giving rise to multiple branches of spirals with cell growth in the tip of the branches and in the radial direction. Unlike the quadratic population growth (ρ ∼ t2)—equivalently constant speed radial growth—in the sector pattern case discussed in the previous section, the overall population in the spiral pattern case grows more slowly due to both nutrient and spatial limitation. In Fig 9, we see that both populations grow ∼ t3/2, and the colony radius r expands sub-linearly in time, indicating a decreasing radial advancing speed.
[Figure omitted. See PDF.]
A)-D) Colony snapshots from t = 100 to 250. Equal number of species 1 and 2 are seeded in a 20 × 20 patch with initial density ρ0 = 1/4. λ1,2 = 1, DA,B = 1, γA,B = 1, YA,B = 1. Scale bar indicates the width of 100 cells.
[Figure omitted. See PDF.]
A) Colony population ρ1,2(t) and B) colony radius r(t) (right). Equal number of species 1 and 2 are seeded in a 20 × 20 patch with initial density ρ0 = 1/4. λ1,2 = 1, DA,B = 1, γA,B = 1, YA,B = 1.
To further investigate the colony radial growth, we examine the nutrient profiles together with the colony morphology, shown in Fig 10. We observe that, with the slow diffusion DA,B = 1, the regions with metabolites closely mirror the spatial distribution of the producer cells. The nutrient depletion regions also mirror the spatial distribution of the consumers. In addition, the bounds of the colony does not exceed that of the nutrient profiles. Collectively, this indicates that the nutrient diffusion limits the colony expansion. The tip of the spiral has a low metabolite concentration because of the newly founded producers. The nutrient profiles shown in Fig 10B and 10C are the net result of nutrient production and consumption. Within the cross-section of a branch of producers, the overall nutrient production per unit time is constant due to the fixed number of producers confined in space. With regard to consumption, the diffusion of nutrient determines how far the consumers can grow. As visualized in Fig 10B and 10C, there is a build-up of nutrient in the center of the cross-section while the concentration drops to zero in the radial directions due to complete consumption of nutrients that diffuse into the region. We will see next that this indeed impacts the colony morphology when individual growth rates are increased.
[Figure omitted. See PDF.]
A) Snapshots of the colony morphology at t = 300, 400, 500 and 1000. B) and C) show the corresponding nutrient profiles. λ1,2 = 1, DA,B = 1, γA,B = 1, YA,B = 1. Scale bar indicates a width of 100 cells.
The characteristics and evolution of the spiral patterns depend on the single cell growth rate λ1,2 when we keep the other parameters the same. In Fig 11, snapshots of the colony are taken at t = 1000 for growth rates λ1,2 = 2, 5 and 10 along with λ1,2 = 1 for comparison. With larger growth rates, the overall colony is larger in radius after growing for the same amount of time. In Fig 12, we can see a slight growth rate dependence in the colony radial expansion speed dr(t)/dt, and in all cases the radial expansion is slowing down over time. Unlike the sector patterns discussed earlier, the radial expansion in the spiral patterns comes mainly from the branch wrapping around. We can see from Fig 11 that, with a higher growth rate, the colony has a “thinner” wrapping branch. As alluded to previously, the steady state of nutrient concentration is established within a spiral branch. Since consumers can only get metabolites through diffusion DA,B and the time scale of significance is the growth rate λ1,2, the characteristic distance for the metabolite diffusion, and thus the characteristic width of the branch, can be estimated using dimensional analysis as at a given time t. Therefore, as λ1,2 increases, we see a “thinner” wrapping branch emerging.
[Figure omitted. See PDF.]
Colony snapshots at t = 1000 with λ1,2 = A) 1, B) 2, C) 5, D) 10. In all cases, DA,B = 1, γA,B = 1, YA,B = 1. Scale bar indicates a width of 100 cells.
[Figure omitted. See PDF.]
λ1,2 = 1, 2, 5, 10. In all cases, DA,B = 1, γA,B = 1, YA,B = 1. The power-law fit is performed to times after a clear spiral is established.
Moreover, the width of the growing branches show interesting dynamics. In our simulation, we measure the spiral branch width at a given time in the following way: We take the stretch of the branch still in contact with the open space, calculate the distance between the inner (in contact with the existing colony) and the outer (in contact with the open space) parts, and average over the entire stretch. As the colony continues to evolve, we see that w(t)∼t1/2 for λ = 1 in Fig 13. The dynamics of the branch width w(t) also shows a slight growth rate dependence.
[Figure omitted. See PDF.]
λ1,2 = 1, 2, 5, 10. In all cases, DA,B = 1, γA,B = 1, YA,B = 1. The power-law fit is performed to times after a clear spiral is established.
In addition to characterizing the growing colony in terms of the colony radius r(t) and wrapping branch width w(t), we examine the angular position of the advancing branch tips. Setting a reference time t0 when a clear spiral pattern emerges, we record the angular position of the advancing tips as θi(t0), i = 1, 2 in the case with 2 advancing spiral tips. We then measure the subsequent angular positions θi(t) with respect to θi(t0). For example, after a transient of multiple interfaces merging around t = 250 in Fig 8, we set the angular position of the tips of the branches to be the reference angular position.
In a stable spiral, we observe at least two advancing branches, one of each species. It is possible to have more as we will see later. This means both species can stably co-exist. We also note that as interfaces stochastically coalesce, the remaining growing tips of the two species settle into a stable angular separation around π, the maximum value in the case where there are two growing branches. At first, the growing tips are at random angular positions of the colony, depending on where the interfaces coalesce. When t = 250 in Fig 8, for instance, the angular separation between the two species’ advancing front is between π/2 and π. In this case, the consumer farther behind the producer (green species shown in Fig 10) temporarily advance faster until the phase difference between the two advancing tips reaches the maximal value of π. The phase separation between the two advancing fronts then remains at π given the symmetric choices of the parameters. The long-time dynamics of the tip angular position θ(t) also shows a slight growth rate dependence, as shown in Fig 14. We show only one advancing front because the behavior of the other one is the same.
[Figure omitted. See PDF.]
λ1,2 = 1, 2, 5, 10. In all cases, DA,B = 1, γA,B = 1, YA,B = 1. θ(t = 0) = 0 is chosen when the first angular position measurement is taken. The power law fit is performed between t = 103 and 5 × t3.
The branch width w(t), the advancing tip angular position θ(t) and the colony radius r(t) are inter-connected: In the steady state of two growing spirals, the number of newborn cells in a given time interval dt equal to:(6)
The first term is the growth in the advancing tip and the second refers to the radial growth along the wrapping branch. As discussed before, the former contributes most of the colony growth. Our simulation results shown in Figs 12–14 give self-consistent scaling behaviors across the colony characteristics.
The above discussions illustrate the spiral morphology when growth is limited by metabolite diffusion. It is natural to also look into the case when DA,B is increased for a specific growth rate. This time, we are surprised to see distinct branch patterns emerging for very large DA,B. Using λA,B = 10 for faster simulation results, we increased DA,B from 10 to 500, shown in Fig 15. As DA,B increases, the characteristic width of the branches increases as discussed above. Additionally, the outer edge of the colony develops individual branches, reminiscent of diffusion limited aggregations seen in other computational models, for example in Refs. [33, 34]. With higher nutrient diffusion, the characteristic width of the region where the cross-feeding metabolites exist is greater than that of the average width of the consumers. Therefore stochastically some consumer cells grow more radially than others, giving rise to the island structure. It is also worth noting that the spiral colony could have more than 2 growing branches, see for example the DA,B = 100 case in Fig 15 where 4 branches are present. And again, the phase difference between neighboring advancing tips is maximized at 2π/4. Our simulations indicate that the number of branches is a stochastic result set in during the early stage of the colony when interfaces merge. For symmetry reasons, there is always an even number of branches, with two-branch spiral the most likely scenario.
[Figure omitted. See PDF.]
Colony snapshots at t = 400 for different nutrient diffusion DA,B = 10, 100, 250, 500. In all cases, λ1,2 = 10, γA,B = 1, YA,B = 1.
4 Summary and discussions
Revealing and understanding the mechanisms of multi-species co-existence is a central question in ecological pursuits. The species coexistence emerges from spatial structure of their living community. Motivated by a theoretical cross-feeding model, our observations presented in this study suggest that the varying spatial availability of the cross-feeding metabolites can lead to drastically different, stable, and stunning, co-existence colony patterns.
Adopting a cellular automaton model on a square lattice using kinetic Monte Carlo simulation algorithm, we set out to study the spatial features of a two-species mutualistic cross-feeding mechanism. With symmetric parameters for both species, we discover non-trivial colony patterns with stable co-existence of both species as a result of cross-feeding nutrient availability through production and diffusion for a wide range of individual cell growth rates and nutrient diffusion rates. We focus on the effects of λ1,2 because they set the intrinsic timescale of the system. Together with the underlying lattice spacing and nutrient diffusion DA,B, they also set the length scale of the system. Qualitatively, the colony morphology evolves from a stable co-existing expanding sector s when nutrients are plentiful to a stable co-existing spiral with species grow along the interface. In the former regime, neither species need to rely on spatial proximity to the other. The initially mixed inoculum self-separates and grows into radially expanding sectors separated by super-diffusive interfaces. In the latter regime, nutrient availability is localized near the producer cells. Within the growing colony, instead of radially expanding, the interface between the two species bends and each species grow s along the other. Transitioning between the radially extending interfaces and the spiral interfaces, there is also the possibility of interfaces coalesce, resulting in one species being engulfed by the other. Beyond analyzing the macroscopic colony patterns, we also investigated the colony characteristics including colony radius, and the spiral branch width.
Spiral patterns and spiral waves are commonplace in nature, occurring both in biological systems and chemical reactions. What we revealed here is a general mechanism for spiral patterns to emerge in a generic cross-feeding two-species system. The lattice-based model is versatile and readily adaptable to explore the vast parameter space outlined in this work. Here we focus on the symmetric mutualistic interaction between the two species. We have also performed preliminary studies on cases with asymmetric growth rates and diffusion coefficients, discussed in S2 File. We see the preservation of the spiral pattern with small differences between either the growth rates, or the diffusion coefficients. This promises a wider parameter space where we can observe spiral patterns as a result of mutualistic cross-feeding. Beyond our study presented here, we see multiple promising directions for further investigation: On the theoretical front, we can adapt the system to studying other types of cross-feeding such as commensalism and syntrophy. It is also possible to incorporate negative effects of the accumulation of the excess metabolites, as suggested by recent experiments [35], that impact local species diversity. At present, the metabolic rule is such that for cells that do not divide due to limitation of space, they continue to take up and excrete metabolites as a form of cell maintenance. This condition can be revised for different growth limitations and to incorporate cell deaths. The simple construction of this model also serves as an ideal primer to quickly scan the parameter space before extending into more sophisticated, yet computationally more intense, three-dimensional models such as Ref. [36]. In this study, we focus on the diffusion of the cross-feeding metabolites. We will direct our future work in including cell motility inter-cell mechanical interactions to see how that impacts spatial structures of the colony.
Supporting information
S1 File. Measurement of interface fluctuations.
https://doi.org/10.1371/journal.pone.0276268.s001
S2 File. Colony patterns with asymmetric mutualism.
https://doi.org/10.1371/journal.pone.0276268.s002
S3 File. All codes used to generate the data are published on Github: https://github.com/mr7Jacky/mutualistic-crossfeeding.
https://doi.org/10.1371/journal.pone.0276268.s003
Acknowledgments
The authors acknowledge the financial support from the National Science Foundation through grants DMR-1702321, MCB-2029480, and MCB-2029580. This work benefited from fruitful discussions with T. Hwa and B. Li. JJD is grateful for the hospitality of T. Hwa at UCSD during her sabbatical visit, where this project initiated.
Citation: Lin J, Sun H, Dong J (2022) Emergence of sector and spiral patterns from a two-species mutualistic cross-feeding model. PLoS ONE 17(10): e0276268. https://doi.org/10.1371/journal.pone.0276268
About the Authors:
Jiaqi Lin
Roles: Data curation
Affiliation: Department of Computer Science, Bucknell University, Lewisburg, Pennsylvania, United States of America
https://orcid.org/0000-0003-1383-6723
Hui Sun
Roles: Data curation, Software, Supervision, Writing – review & editing
Affiliation: Department of Mathematics, California State University, Long Beach, California, United States of America
JiaJia Dong
Roles: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing
E-mail: [email protected]
Affiliation: Department of Physics & Astronomy, Bucknell University, Lewisburg, Pennsylvania, United States of America
https://orcid.org/0000-0001-9908-2699
1. Konopka A. What is microbial community ecology? The ISME journal. 2009;3(11):1223–1230. pmid:19657372
2. Davey ME, O’toole GA. Microbial biofilms: from ecology to molecular genetics. Microbiology and molecular biology reviews. 2000;64(4):847–867. pmid:11104821
3. Tilman D. The ecological consequences of changes in biodiversity: a search for general principles. Ecology. 1999;80(5):1455–1474.
4. Stams AJ, Plugge CM. Electron transfer in syntrophic communities of anaerobic bacteria and archaea. Nature Reviews Microbiology. 2009;7(8):568–577. pmid:19609258
5. Morris BE, Henneberger R, Huber H, Moissl-Eichinger C. Microbial syntrophy: interaction for the common good. FEMS microbiology reviews. 2013;37(3):384–406. pmid:23480449
6. Blanchard AE, Lu T. Bacterial social interactions drive the emergence of differential spatial colony structures. BMC systems biology. 2015;9(1):1–13. pmid:26377684
7. Smith NW, Shorten PR, Altermann E, Roy NC, McNabb WC. The classification and evolution of bacterial cross-feeding. Frontiers in Ecology and Evolution. 2019;7:153.
8. Ghosh AR. Appraisal of microbial evolution to commensalism and pathogenicity in humans. Clinical Medicine Insights: Gastroenterology. 2013;6:CGast–S11858. pmid:24833938
9. Excoffier L, Foll M, Petit RJ. Genetic consequences of range expansions. Annual Review of Ecology, Evolution, and Systematics. 2009;40:481–501.
10. Momeni B, Brileya KA, Fields MW, Shou W. Strong inter-population cooperation leads to partner intermixing in microbial communities. elife. 2013;2:e00230. pmid:23359860
11. Dao H, Lakhani P, Police A, Kallakunta V, Ajjarapu SS, Wu KW, et al. Microbial stability of pharmaceutical and cosmetic products. Aaps Pharmscitech. 2018;19(1):60–78. pmid:29019083
12. Tobor-Kapłon MA, Bloem J, Römkens PF, Ruiter Pd. Functional stability of microbial communities in contaminated soils. Oikos. 2005;111(1):119–129.
13. Wolfram S. Statistical mechanics of cellular automata. Rev Mod Phys. 1983;55(3):601.
14. Wolfram S. Cellular automata as models of complexity. Nature. 1984;311(5985):419–424.
15. Ermentrout GB, Edelstein-Keshet L. Cellular automata approaches to biological modeling. Journal of theoretical Biology. 1993;160(1):97–133. pmid:8474249
16. Voter AF. Introduction to the kinetic Monte Carlo method. In: Radiation effects in solids. Springer; 2007. p. 1–23.
17. Gillespie DT. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. Journal of Computational Physics. 1976;22(4):403–434.
18. Gillespie DT. Exact stochastic simulation of coupled chemical reactions. The Journal of Physical Chemistry. 1977;81(25):2340–2361.
19. Andersen M, Panosetti C, Reuter K. A practical guide to surface kinetic Monte Carlo simulations. Frontiers in chemistry. 2019;7:202. pmid:31024891
20. Donev A, Bulatov VV, Oppelstrup T, Gilmer GH, Sadigh B, Kalos MH. A first-passage kinetic Monte Carlo algorithm for complex diffusion–reaction systems. Journal of Computational Physics. 2010;229(9):3214–3236.
21. Katsoulakis MA, Vlachos DG. Coarse-grained stochastic processes and kinetic Monte Carlo simulators for the diffusion of interacting particles. The Journal of chemical physics. 2003;119(18):9412–9427.
22. Piana S, Gale J. Three-dimensional kinetic Monte Carlo simulation of crystal growth from solution. Journal of crystal growth. 2006;294(1):46–52.
23. Reuter K, Scheffler M. First-principles kinetic Monte Carlo simulations for heterogeneous catalysis: Application to the CO oxidation at Ru O 2 (110). Physical Review B. 2006;73(4):045433.
24. Neidhardt FC, Ingraham JL, Schaechter M. Physiology of the bacterial cell; a molecular approach. 589.901 N397. Sinauer associates; 1990.
25. Monod J. The growth of bacterial cultures. Annual review of microbiology. 1949;3(1):371–394.
26. Eden M, et al. A two-dimensional growth process. In: Proceedings of the fourth Berkeley symposium on mathematical statistics and probability. vol. 4. University of California Press Berkeley; 1961. p. 223–239.
27. Barabási AL, Stanley HE. Fractal concepts in surface growth. Cambridge university press; 1995.
28. Hallatschek O, Hersen P, Ramanathan S, Nelson DR. Genetic drift at expanding frontiers promotes gene segregation. PNAS. 2007;104(50):19926–19930. pmid:18056799
29. Müller MJ, Neugeboren BI, Nelson DR, Murray AW. Genetic drift opposes mutualism during spatial population expansion. Proceedings of the National Academy of Sciences. 2014;111(3):1037–1042. pmid:24395776
30. Korolev K, Nelson DR. Competition and cooperation in one-dimensional stepping-stone models. Physical Review Letters. 2011;107(8):088103. pmid:21929209
31. Lavrentovich MO, Nelson DR. Asymmetric mutualism in two-and three-dimensional range expansions. Physical review letters. 2014;112(13):138102. pmid:24745459
32. Menon R, Korolev KS. Public good diffusion limits microbial mutualism. Physical review letters. 2015;114(16):168102. pmid:25955075
33. Nadell CD, Foster KR, Xavier JB. Emergence of spatial structure in cell groups and the evolution of cooperation. PLoS Comput Biol. 2010;6(3):e1000716. pmid:20333237
34. Tronnolone H, Tam A, Szenczi Z, Green J, Balasuriya S, Tek EL, et al. Diffusion-limited growth of microbial colonies. Scientific Reports. 2018;8(1):1–11. pmid:29662092
35. Goldschmidt F, Regoes RR, Johnson DR. Metabolite toxicity slows local diversity loss during expansion of a microbial cross-feeding community. The ISME journal. 2018;12(1):136–144. pmid:28914879
36. Warren MR, Sun H, Yan Y, Cremer J, Li B, Hwa T. Spatiotemporal establishment of dense bacterial colonies growing on hard agar. Elife. 2019;8:e41093. pmid:30855227
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer
© 2022 Lin et al. This is an open access article distributed under the terms of the Creative Commons Attribution License: http://creativecommons.org/licenses/by/4.0/ (the “License”), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.
Abstract
The ubiquitous existence of microbial communities marks the importance of understanding how species interact within the community to coexist and their spatial organization. We study a two-species mutualistic cross-feeding model through a stochastic cellular automaton on a square lattice using kinetic Monte Carlo simulation. Our model encapsulates the essential dynamic processes such as cell growth, and nutrient excretion, diffusion and uptake. Focusing on the interplay among nutrient diffusion and individual cell division, we discover three general classes of colony morphology: co-existing sectors, co-existing spirals, and engulfment. When the cross-feeding nutrient is widely available, either through high excretion or fast diffusion, a stable circular colony with alternating species sector emerges. When the consumer cells rely on being spatially close to the producers, we observe a stable spiral. We also see one species being engulfed by the other when species interfaces merge due to stochastic fluctuation. By tuning the diffusion rate and the growth rate, we are able to gain quantitative insights into the structures of the sectors and the spirals.
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer
Details
; Sun, Hui; Dong, JiaJia