Emergence of sector and spiral patterns
from a two-species mutualistic cross-feeding model

Jiaqi.Lin, Hui. Sun, JiaJia Dong^∗ Department of Computer Science, Bucknell University Department of Mathematics, Cal State, Long Beach Department of Physics & Astronomy, Bucknell University [email protected]

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.

I 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, 5, 6, 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, 14, 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, 17, 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 II, we provide details of our model and simulation algorithm. We characterize the distinct colony patterns in Section III and analyze the interplay between cell growth and cross-feeding. Finally, we discuss in Section IV 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.

II 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 $\lambda_{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 $\gamma_{A,B}$ and the excreted molecule, if not taken up, diffuses with a diffusion coefficient of $D_{A,B}$ . The yield $Y_{A,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:

$\displaystyle\frac{\partial{\rho}_{1}}{\partial t}$	$\displaystyle=\lambda_{1}(n_{A})\rho_{1};~{}~{}\frac{\partial{\rho}_{2}}{\partial t}=\lambda_{2}(n_{B})\rho_{2};$	(1)
$\displaystyle\frac{\partial n_{A}}{\partial t}$	$\displaystyle=\gamma_{A}\rho_{2}-\lambda_{1}(n_{A})\rho_{1}/Y_{A}+D_{A}\cdot\nabla^{2}n_{A};$	(2)
$\displaystyle\frac{\partial n_{B}}{\partial t}$	$\displaystyle=\gamma_{B}\rho_{1}-\lambda_{2}(n_{B})\rho_{2}/Y_{B}+D_{B}\cdot\nabla^{2}n_{B}.$	(3)

Note that the cell concentration $\rho_{1,2}$ and the nutrient concentration $n_{A,B}$ depend on both space and time. The parameters and reactions of our model are summarized in a schematic shown in Fig.1.

Refer to caption — Figure 1: Mutualistic cross-feeding scheme. As species 1 cells metabolize nutrient molecules $A$ (produced by species 2) at a growth rate $\lambda_{1}$ , they excrete molecules $B$ with rate $\gamma_{B}$ as a metabolic by-product. Molecule $B$ is the essential nutrient for species 2 cells’ growth and is taken up with yield $Y_{B}$ . 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 $\vec{r}(x,y)$ is chosen to divide, the individual growth rate $\lambda_{1,2}$ is determined by $n_{A,B}$ , the total nutrient concentration at $\vec{r}$ and its four nearest-neighbor (n.n.) lattice sites:

\lambda_{1,2}=\lambda_{1,2}^{\ast}\dfrac{\sum_{\rm{n.n.}}n_{A,B}}{\sum_{\rm{n.n.}}n_{A,B}+K_{A,B}}

(4)

Here $\lambda_{1,2}^{\ast}$ indicates the maximal single cell growth rate. The instantaneous growth rate is modified by a Monod factor $K_{A,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\times L$ lattice sites with $\rho_{0}L^{2}$ cells, half of each species randomly distributed within the patch. $L^{2}$ is much smaller than the final colony size and $\rho_{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, $\rho_{1,2}$ is the total number of cell type $1$ (or $2$ ), while $n_{A,B}$ is the total number of nutrient molecule $A$ (or $B$ ) in the entire system at a given time $t$ . The propensity $k_{j}$ 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 $k_{0}$ is the sum of all propensities.

#	reaction	propensity	notes
1	type 1 cell divide	$k_{1}=\lambda_{1}\rho_{1}$	$\lambda_{1}=\lambda_{1}^{\ast}\dfrac{n_{A}}{n_{A}+K_{A}}$ ;
2	type 1 cell excrete $n_{B}$	$k_{2}=\gamma_{B}\rho_{1}$
3	type 1 cell take up $n_{A}$	$k_{3}=\dfrac{\lambda_{1}}{Y_{A}}\rho_{1}$
4	type 2 cell divide	$k_{4}=\lambda_{2}\rho_{2}$	$\lambda_{2}=\lambda_{2}^{\ast}\dfrac{n_{B}}{n_{B}+K_{B}}$ ;
5	type 2 cell excrete $n_{A}$	$k_{5}=\gamma_{A}\rho_{2}$
6	type 2 cell take up $n_{B}$	$k_{6}=\dfrac{\lambda_{2}}{Y_{B}}\rho_{2}$
7	nutr. $B$ diffuses	$k_{7}=D_{B}n_{B}$
8	nutr. $A$ diffuses	$k_{8}=D_{A}n_{A}$
		$k_{0}=\sum_{j=1}^{8}k_{j}$

Table 1: Processes involved in the system and the corresponding rates.

Let ${\rm r}_{1},{\rm r}_{2}\in(0,1)$ be two independent, uniformly distributed random numbers. At time $t$ , the $j$ -th reaction ( $j\in[1,8]$ ) takes place if:

\sum_{i=1}^{j-1}k_{i}<k_{0}\cdot{\rm r}_{1}\leq\sum_{i=1}^{j}k_{i}.

Afterwards $t$ is advanced by $\Delta t$ where

\Delta t=-\dfrac{\ln{\rm r}_{2}}{k_{0}}.

The iteration continues until the system reaches a final time $T_{\rm max}$ , 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 $\lambda_{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.

III Results

III.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 both experimentally [28, 29] and theoretically [30, 31, 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 $\gamma_{A,B}$ and high diffusion coefficients $D_{A,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\times 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 $\rho_{1,2}(t)$ over time. In the early stage ( $t\lesssim 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 $\rho\approx\pi r^{2}$ and the number of cells at the expanding frontier $2\pi r$ . The radial expansion of the colony is approximately:

\dfrac{d}{dt}\left(\pi r^{2}\right)=\lambda\cdot 2\pi r\rightarrow\dfrac{dr}{dt}=\lambda,\rho\sim~{}t^{2}

(5)

This gives a constant colony expansion speed that depends on the individual growth rate $\lambda$ and the total population within the colony grows quadratic with time.

In Fig.2, we show the population $\rho_{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 $t^{2.11}$ . The expansion speed $dr/dt$ , measured to be 1.17, is consistent with the individual growth rate $\lambda_{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.

For an interface that fluctuates like a standard random walk, we expect the mean-square displacement to scale linearly with time, $\overline{y^{2}}\sim t$ . In our case, $\overline{y^{2}}\sim t^{4/3}$ 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.

Now we turn to the second case to see how nutrient excretion rate $\gamma_{A,B}$ under high nutrient diffusion $D_{A,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 $\gamma_{A,B}=10$ and high nutrient diffusion rate $D_{A,B}=500$ . Similar to Fig.3, the colony again separates into sectors of same-type species. The interface mean-square displacement $\overline{y^{2}}\sim t^{1.33}$ , same as in Fig.4. Fig.6 shows the cell population $\rho_{1,2}(t)$ and colony radius $r(t)$ . The overall population displays a quadratic growth while $r^{2}$ increases faster than $\rho$ . 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.

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 $D_{A,B}=500$ as in Fig.5 while reducing $\gamma_{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, $\overline{y^{2}}\sim r^{2}\sim t^{2.5}$ . 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.

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 $\overline{y^{2}}$ grows faster than $t^{2}$ , shown in Fig.7(D). This means that as time increases, the two interfaces will eventually coalesce.

III.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 $\rho_{0}=1/4$ as before, we set the excretion rate $\gamma_{A,B}$ and the diffusion rate $D_{A,B}$ to be 1. Here with $\gamma_{A,B}=Y_{A,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\lesssim 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 $n_{A}$ (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 ( $\rho\sim t^{2}$ ) – 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 $\sim t^{3/2}$ , and the colony radius $r$ expands sub-linearly in time, indicating a decreasing radial advancing speed.

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 $D_{A,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 C 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 C, 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.

The characteristics and evolution of the spiral patterns depend on the single cell growth rate $\lambda_{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 $\lambda_{1,2}=2,5$ and $10$ along with $\lambda_{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 $D_{A,B}$ and the time scale of significance is the growth rate $\lambda_{1,2}$ , the characteristic distance for the metabolite diffusion, and thus the characteristic width of the branch, can be estimated using dimensional analysis as $w(t)\sim\sqrt{2D_{A,B}/\lambda_{1,2}}$ at a given time $t$ . Therefore, as $\lambda_{1,2}$ increases, we see a “thinner” wrapping branch emerging.

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)\sim t^{1/2}$ for $\lambda=1$ in Fig.13. The dynamics of the branch width $w(t)$ also shows a slight growth rate dependence.

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 $t_{0}$ when a clear spiral pattern emerges, we record the angular position of the advancing tips as $\theta_{i}(t_{0})$ , $i=1,2$ in the case with 2 advancing spiral tips. We then measure the subsequent angular positions $\theta_{i}(t)$ with respect to $\theta_{i}(t_{0})$ . 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 $\pi$ , 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 $\pi/2$ and $\pi$ . 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 $\pi$ . The phase separation between the two advancing fronts then remains at $\pi$ given the symmetric choices of the parameters. The long-time dynamics of the tip angular position $\theta(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.

The branch width $w(t)$ , the advancing tip angular position $\theta(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:

d\rho\approx 2\times\left(w(t)\cdot r(t)\cdot d\theta+\pi r(t)\cdot dw\right)

(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, 13, and 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 $D_{A,B}$ is increased for a specific growth rate. This time, we are surprised to see distinct branch patterns emerging for very large $D_{A,B}$ . Using $\lambda_{A,B}=10$ for faster simulation results, we increased $D_{A,B}$ from 10 to 500, shown in Fig.15. As $D_{A,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 $D_{A,B}=100$ case in Fig.15 where 4 branches are present. And again, the phase difference between neighboring advancing tips is maximized at $2\pi/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.

IV 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 $\lambda_{1,2}$ because they set the intrinsic timescale of the system. Together with the underlying lattice spacing and nutrient diffusion $D_{A,B}$ , they also set the length scale of the system. Qualitatively, the colony morphology evolves from a stable co-existing expanding sector 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 grows 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. 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 Measurement of interface fluctuations.

Measurement of interface fluctuations. We develop the algorithm in MATLAB code. To characterize the fluctuating boundaries separating the two species, we define the lattice site $p_{i}$ along the interface as the following: Each lattice site can take on values 0 (empty), or 1 or 2 depending on the occupant cell type. $p_{i}$ is on the interface if it is occupied by a cell and there are only two unique non-zero values among itself and its four nearest-neighbor sites. The collection of points $\{p_{i}\}$ defines the interface and its length is $L$ .

For each interface, we quantify the fluctuation along the interface as a function of a sliding window $x\in[1,L]$ . Given $x$ , We first find the best linear fit $y(x^{\prime})$ for the segment connecting from $p_{0}$ to $p_{x-1}$ . We then compute the mean square displacement $y_{0}^{2}$ between the interface and $y(x^{\prime})$ . This process is repeated by moving the segment along the interface one point at a time. Then the average of $y_{i}^{2}$ gives us the mean square displacement $\overline{y^{2}}(x)$ for a given window size $x$ .

S2 Colony patterns with asymmetric mutualism.

Colony patterns with asymmetric mutualism. Our work focuses on the mutualistic cross-feeding case where all parameters are the same for both species. However, the spiral patterns are robust even when the parameters are asymmetric. Below we show a few scenarios where the nutrient diffusion rates and the growth rates are different for the two species. Given the large parameter space, we will reserve a comprehensive study on all possible asymmetric configurations in another study.

As one of the nutrient diffusion rate, $D_{B}$ (molecules excreted by species 1 to be taken up by species 2), is reduced while keeping $D_{A}$ the same, we see the system continue to emerge into a stable spiral pattern, shown in Fig.16. When $D_{B}\ll D_{A}$ with the same excretion rate, the slower-diffusing nutrient leads to local concentration lower than the Monod constant $K_{B}$ , thus reducing the effective growth rate for species 2 (blue cells). In this case, species 2 lose out in competing for available growth space and is engulfed by the faster growing species 1, as shown in Fig.16(A). The decrease in one of the nutrient diffusion rates also led to the decrease in overall colony size due to the coupling of individual growth rates $\lambda_{1,2}$ and local nutrient concentration $n_{A,B}$ .

When the maximal individual cell growth rate $\lambda_{1,2}$ is varied, we again observe a stable spiral pattern when $\lambda_{1}$ and $\lambda_{2}$ are comparable. In Fig.17(A-B), species 2 (blue cells) have a much smaller maximal growth rate. At the nascent stage of the colony, the limited number of empty lattice sites are taken up by species 1 (green cells), resulting in the engulfment pattern. When $\lambda_{2}$ is increased, we see the emergence of a stable spiral pattern again as shown in Fig.17(C-D). It is worth noting that when the two species have different growth rates, the slower growing cells (blue in Fig.17) develop a wider wrapping branch, which is consistent with what we discussed in the main article.

S3 All codes used to generate the data are published onGithub

All codes used to generate the data are published on 365 Github: https://github.com/mr7Jacky/mutualistic-crossfeeding

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.

References

[1] Konopka A. What is microbial community ecology? The ISME journal. 2009;3(11):1223–1230.
[2] Davey ME, O’toole GA. Microbial biofilms: from ecology to molecular genetics. Microbiology and molecular biology reviews. 2000;64(4):847–867.
[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.
[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.
[6] Blanchard AE, Lu T. Bacterial social interactions drive the emergence of differential spatial colony structures. BMC systems biology. 2015;9(1):1–13.
[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.
[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.
[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.
[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.
[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.
[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. doi:10.1073/pnas.0710150104.
[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.
[30] Korolev K, Nelson DR. Competition and cooperation in one-dimensional stepping-stone models. Physical Review Letters. 2011;107(8):088103.
[31] Lavrentovich MO, Nelson DR. Asymmetric mutualism in two-and three-dimensional range expansions. Physical review letters. 2014;112(13):138102.
[32] Menon R, Korolev KS. Public good diffusion limits microbial mutualism. Physical review letters. 2015;114(16):168102.
[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.
[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.
[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.
[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.

Emergence of sector and spiral patterns from a two-species mutualistic cross-feeding model