Structured Interactions Drive Abrupt Transitions in the Spatial Organization of Microbial Communities

TL;DR

This work addresses how multispecies microbial communities organize spatially when interspecies motility regulation is local and structured. The authors develop a minimal run-and-tumble model with a sigmoidal motility response controlled by an interspecies interaction matrix $A_{SS'}$, and they analyze both a two-species case and large communities. They find a sharp, nucleation-driven transition from a well-mixed motile state to large stationary clusters as the density of motility-suppressing links increases; in large communities, transition propensity can be dictated by network topology, with modular and heterogeneous matrices promoting clustering at lower link densities. The results imply that structured ecological interactions, not just abundances or overall interaction density, shape microbial patchiness and can inform inference of interspecies networks from spatial patterns and biofilm development.

Abstract

Bacteria possess diverse mechanisms to regulate their motility in response to environmental and physiological signals, enabling them to navigate complex habitats and adapt their behavior. Among these mechanisms, interspecies recognition enables cells to modulate their movement based on the ecological identity of neighboring species. Here, we introduce a model in which we assume bacterial species recognizes each other and interact via local signals that either enhance or suppress the motility of neighboring cells. Through large-scale simulations and a coarse-grained stochastic model, we demonstrate the emergence of a sharp transition driven by nucleation processes: increasing the density of motility-suppressing interactions drives the system from a fully mixed, motile phase to a state characterized by large, stationary bacterial clusters. Remarkably, in systems with a large number of interacting species, this transition can be triggered solely by altering the structure of the motility-regulation interaction matrix while maintaining species and interaction densities constant. In particular, we find that heterogeneous and modular interactions promote the transition more readily than homogeneous random ones. These results contribute to the ongoing effort to understand microbial interactions, suggesting that structured, non-random ones may be key to reproducing commonly observed spatial patterns in microbial communities.

Figures (10)

• Figure 1: (A) Mean number of non-motile bacteria at equilibrium as a function of the fraction of bacteria belonging to species "a”. Blue squares indicate the mean values from 100 runs of the microscopic simulation, with error bars representing standard deviations. The orange line shows results obtained from the mean of 100 realizations of the theoretical stochastic model simulated using the Gillespie algorithm. The "critical” region, where the transition may stochastically occur or not, is highlighted in slight red. The inset shows the distribution of the results for each of the 100 simulations for $n_a/n=0.23$. All the simulations here are performed with a total number of bacteria $n=16000$, an interaction radius of $R=0.025$ and a maximum velocity $v_0=0.05$. (B) Top: temporal dynamics of the fraction of motile bacteria from a representative run of the microscopic simulation at $n_a/n = 0.25$. Bottom: graphical representations of the system at three different time steps, highlighting the nucleation process. Bacteria of species "a" and "b" are depicted in red and black respectively.
• Figure 2: Steady-state spatial configurations of $n = 100000$ bacteria partitioned into $N = 100$ different species, shown for two distinct values of the number of motility suppression links. When the number of suppression links is insufficient, the system remains in a fully motile, well-mixed phase. Upon reaching a critical value, the system may undergo a transition to a stationary clustered state in which all bacteria cease motion. Simulations are performed with $R = 0.01$ and $v_0 = 0.05$ within a two-dimensional box of size $[0,1] \times [0,1]$.
• Figure 3: Mean fraction of non-motile bacteria at equilibrium $n^S/n$ as a function of the density of species-specific motility suppression interactions $2L_-/N(N-1)$. The orange and the violet curves correspond to a homogeneous interactions matrix generated by Erdős–Rényi (ER) model and a heterogeneous one provided by the Barabási–Albert (BA) model. Dots indicate averages over 100 independent microscopic simulations, with error bars representing standard deviations. The yellow region highlights the bistable domain where the system’s phase depends solely on the structure of the interactions matrix. All simulations are performed within a two-dimensional box of size $[0,1] \times [0,1]$ and with $n = 100000$ bacteria, $N = 100$ species, $R = 0.01$, and $v_0 = 0.05$.
• Figure 4: Heatmap showing the average number of non-motile bacteria at equilibrium for modular motility suppression matrices generated using a Stochastic Block Model. The x-axis represents $p_{\text{in}}$, the probability of forming interactions within groups, while the y-axis corresponds to $p_{\text{out}}$, the probability of establishing interactions between different groups. White dashed lines denote curves of constant interactions density, indicating configurations with identical numbers of interactions but different combinations of $p_{\text{in}}$ and $p_{\text{out}}$. The red dotted line represents the scenario of fully homogeneous interactions, where $p_{\text{in}} = p_{\text{out}}$. A step size of 0.02 was used for both $p_{\text{in}}$ and $p_{\text{out}}$, and each grid cell displays the mean value computed from 100 independent microscopic simulations. All simulations were conducted within a two-dimensional domain of size $[0,1] \times [0,1]$, with $n = 100000$ bacteria, $N = 100$ species, $R = 0.01$, and $v_0 = 0.05$.
• Figure 5: On the left, the complementary cumulative distribution function (CCDF) for the number of bacteria in a cell of size $R \times R$ ($R = 0.01$). The different solid curves correspond to various densities of suppression interactions, spanning from the critical to the super-critical regime (see Fig. 3 in the main text). Dashed lines indicate the best fits using a discrete power-law probability distribution. On the right, stationary spatial patterns are displayed for $C = 0.32, 0.34, 0.36, 0.38$. Each color represents a distinct species and the homogeneous interactions matrix is generated according to the Erdős-Renyi model. All simulations are performed in a $[0,1] \times [0,1]$ domain with $n = 100000$, $N = 100$ species, interaction radius $R = 0.01$, and velocity $v_0 = 0.05$.