Microbiome association diversity reflects proximity to the edge of instability

TL;DR

This work tackles how microbiome macroecological laws coexist with interspecific interactions and health/disease states. It introduces the Interacting Stochastic Logistic Model (ISLM), a balanced Gaussian interaction ensemble that preserves single-species statistics while enabling realistic interspecific covariation, analyzed via dynamical mean-field theory and random-matrix stability. A key contribution is the distance-to-instability parameter $g$, inferred from short-time covariances, which places communities on a May-like phase diagram and explains broadening of association patterns as systems approach criticality. Applying this framework to synthetic data, environmental microbiomes, and human gut cohorts shows that real communities operate near the edge of instability and that healthy guts sit closer to this edge (broader associations) than dysbiotic states, providing a fast, data-efficient dynamical marker of microbiome health with potential diagnostic utility.

Abstract

Recent advances in metagenomics have revealed macroecological patterns or "laws" describing robust statistical regularities across microbial communities. Stochastic logistic models (SLMs), which treat species as independent -- akin to ideal gases in physics -- and incorporate environmental noise, reproduce many single-species patterns but cannot account for the pairwise covariation observed in microbiome data. Here we introduce an interacting stochastic logistic model (ISLM) that minimally extends the SLM by sampling an ensemble of random interaction networks chosen to preserve these single-species laws. Using dynamical mean-field theory, we map the model's phase diagram -- stable, chaotic, and unbounded-growth regimes -- where the transition from stable fixed-point to chaos is controlled by network sparsity and interaction heterogeneity via a May-like instability line. Going beyond mean-field theory to account for finite communities, we derive an estimator of an effective stability parameter that quantifies distance to the edge of instability and can be inferred from the width of the distribution of pairwise covariances in empirical species-abundance data. Applying this framework to synthetic data, environmental microbiomes, and human gut cohorts indicates that these communities tend to operate near the edge of instability. Moreover, gut communities from healthy individuals cluster closer to this edge and exhibit broader, more heterogeneous associations, whereas dysbiosis-associated states shift toward more stable regimes -- enabling discrimination across conditions such as Crohn's disease, inflammatory bowel syndrome, and colorectal cancer. Together, our results connect macroecological laws, interaction-network ensembles, and May's stability theory, suggesting that complex communities may benefit from operating near a dynamical phase transition.

Figures (5)

• Figure 1: From ideal-gas-like descriptions to interacting ensembles. Non-interacting models such as the SLM capture marginal statistical patterns of empirical data (Grilli's laws) but fail to reproduce interspecific associations, while Bayesian inference can fit both at the cost of losing interpretability. A. Analogy between physics and ecology: the ideal gas corresponds to the uncoupled SLM, whereas the Van der Waals gas plays the role of the interacting SLM (ISLM) that explicitly includes interactions among species. B. Three macroecological laws observed in empirical microbiomes: (i) for each species, temporal abundances follow a Gamma distribution (abundance-fluctuation distribution, AFD); (ii) across species, mean abundances are approximately lognormally distributed (mean-abundance distribution, MAD); and (iii) variance scales as the square of the mean (Taylor's law with exponent $\approx 2$). These patterns concern marginal statistics and ignore correlations between species’ relative abundances. C. Pairwise associations: two species can be positively correlated ($i$ and $j$) or negatively correlated ($i$ and $k$). Empirical microbiomes display a broad histogram of pairwise correlations, whose width we denote by $\Delta$, much wider than predicted by uncoupled models. D. Statistics of the synthetic ISLM. First row: schematic visualization of the interaction network. Second row: histogram of the non-zero entries of $G_{ij}$ (see Materials and Methods). Third row: distribution of carrying capacities, set as in \ref{['eq: ensemble 2']}. Fourth row: link sparsity; a fraction $1-C$ of links is zero, while a fraction $C$ is non-zero.
• Figure 2: The synthetic model exhibits fixed-point, chaotic, and unbounded dynamical phases.A. Heat map of the largest Lyapunov exponent (LLE) as a function of $\sigma$ and $C$ (negative values are plotted as zero for clarity). Three regimes can be distinguished: a fixed-point phase (FP), a chaotic phase (C), and an unbounded phase (U) in which trajectories blow up. The white solid line and the black dashed line delimit the FP--C and C--U boundaries, respectively. Simulations were performed for communities with $S = 10^3$ species (to suppress finite-size effects), without environmental noise ($\sigma_0 = 0$, the effect of noise is further explored in the SI--Section 5.B), and averaged over $M = 50$ independent realizations of the quenched interactions. B--D. Representative time series illustrating the dynamical behavior in each regime shown in panel A.
• Figure 3: Grilli's laws are preserved in the interacting ISLM model, while approaching the edge of instability broadens the correlation histogram and enhances association diversity.B.1--B.3 Spectral density of the effective interaction matrix in \ref{['eq: Aeff']} for three values of the interaction strength: B.1 $\sigma = 0$ (pure SLM), B.2 $\sigma = 1/\sqrt{8}$, and B.3 $\sigma = 1/\sqrt{2}$. The connectivity is fixed at $C = 1/2$, so panel B.3 lies exactly at the edge of instability. We use $\mu_{\mathrm{LN}} = 0$ and $\sigma_{\mathrm{LN}} = 1$ for the log-normal parameters. C.1--C.3 Abundance fluctuation distribution (AFD) for the same values of $\sigma$. Points show simulation results for a community with $S=200$ species, integrated up to $t_{\max} = 100$ with time step $h = 0.01$, noise variance $\sigma_0 = 1$, and growth time scale $\tau = 1$, averaged over $M = 500$ independent realizations. The black line is the best-fit Gamma distribution. D.1--D.3 Mean abundance distribution (MAD). Green dots are simulation results; the black curve is the log-normal distribution from which the values of $-1/D_{ii}$ are drawn. E.1--E.3 Taylor's law. Points show the empirical relationship between mean and standard deviation from simulations; the black line is the theoretical scaling $\sigma_i \propto \langle x_i\rangle^2$ (plotted as $y = x^2$). H.1--H.3 Association diversity for the three values of $\sigma$. For small and intermediate interaction strengths (H.1--H.2) the correlation histograms are narrow, while at the edge of instability (H.3) they become much broader. Light points correspond to different realizations of the quenched disorder, and the dark line shows their average, highlighting the multitude of possible correlation patterns compatible with the same single-species macroecological laws.
• Figure 4: Empirical biomes are found more likely at the edge of instability. A. Distance to the edge of instability, measured via the estimator of \ref{['eq: g estimator']}. The color map shows the value of $g$. The black line shows the boundary of instability ($\sigma\sqrt{C}=1$). Below this boundary, the fixed point solution is stable and, above it, it becomes linearly unstable. Four different lines are shown in violet, grey, blue and green, indicating each of the four microbiomes analyzed: Seawater (Panel B), River (Panel C), Glacier (Panel D) and Lake (Panel E), respectively. Lines depict the possible combinations of $(\sigma, C)$ giving rise to the reported values of $g$, measured for each biome. B-E Histograms of correlations for the four microbiomes, each with their own inferred value of $g$. Points show empirical values; shades show the best fit found by simulating the dynamics for $S=200$ and $M=10^3$ different realizations of the quenched disorder at the value of $\sigma$ needed to reproduce the empirical $g$ (with fixed $C=0.5$).
• Figure 5: Dysbiosis-associated conditions are farther from the critical point than healthy cohorts. (A--C) Distributions of pairwise species--species correlations $\rho(C_{ij})$ for Healthy controls (A), Ulcerative Colitis (B), and Crohn's disease (C) gut microbiomes. For each cohort we generated $M=1000$ independent subsamples; each subsample contains $T_2=200$ randomly selected samples (without replacement within a subsample). An occupancy filter of $0.2$ was applied so that only taxa present in at least $20\%$ of the selected samples were retained, yielding $\sim$100 taxa per dataset. (D) Distributions of the estimated distance to criticality ${g}$ for the four groups (two healthy cohorts (group 1 and group 2); two disease cohorts (UC, CD). Brackets report the AUC p-value between (i) the two disease cohorts (U--U), (ii) the two healthy cohorts (H--H), and (iii) pooled Unhealthy versus pooled Healthy (U--H). Small $p$-values indicate evidence for distributional differences; the U--H contrast shows clear separation (at least for the available data), whereas H--H is not significant and U--U shows a modest shift. As an interpretable effect size, the AUC equals the probability that a randomly chosen value from one group exceeds a value from the other; in this setting it quantifies separability of both distributions. The inset (Reshuffled) shows $g$ after random reshuffling that destroys the association between taxa and samples; the distinctive cohort structure disappears, confirming that the observed separations arise from genuine ecological organization rather than sampling artifacts.