Self-organized biodiversity and species abundance distribution patterns in ecosystems with higher-order interactions

TL;DR

The paper extends the Generalized Lotka-Volterra framework to include higher-order interactions (HOIs) and a small dispersal rate, showing HOIs stabilize coexistence and generate diverse dynamics such as equilibria, oscillations, and chaos. Through analytical stability assessments (Jacobian eigenvalues and a Lyapunov-based global stability function) and extensive simulations, the authors demonstrate that HOIs prevent ecosystem collapse and reproduce universal species abundance distributions observed across multiple ecological communities. The findings indicate HOIs are a generic, mechanistic driver of self-organized biodiversity, capable of reproducing empirical rank-abundance patterns and explaining variability in community structure beyond pairwise interactions. The work provides a general, applicable framework for modeling complex ecological systems and highlights the importance of incorporating HOIs in biodiversity theory and ecological forecasting.

Abstract

Explaining the emergence of self-organized biodiversity and species abundance distribution patterns remians a fundamental challenge in ecology. While classical frameworks, such as neutral theory and models based on pairwise species interactions, have provided valuable insights, they often neglect higher-order interactions (HOIs), whose role in stabilizing ecological communities is increasingly recognized. Here, we extend the Generalized Lotka-Volterra framework to incorporate HOIs and demonstrate that these interactions can enhance ecosystem stability and prevent collapse. Our model exhibits a diverse range of emergent dynamics, including self-sustained oscillations, quasi-periodic (torus) trajectories, and intermittent chaos. Remarkably, it also reproduces empirical species abundance distributions observed across diverse natural communities. These results underscore the critical role of HOIs in structuring biodiversity and offer a broadly applicable theoretical framework for capturing complexity in ecological systems

Figures (18)

• Figure 1: Higher-order interactions prevent ecosystems collapse. (a) Schematic representation of the generalized model incorporating both pairwise and higher-order interactions among $S$ consumer species. (b) Species collapse occurs when dynamics are governed solely by pairwise interactions. (c) The introduction of higher-order interactions stabilizes the system, enabling persistent coexistence through self-organized dynamics, even under the same initial conditions as in (b). (d-g) Representative time series of species abundances from simulations with $S = 32$ species, illustrating the emergence of diverse dynamical regimes. For full simulation details, see SM Sec. \ref{['Simulation details']}.
• Figure 2: Emergence of chaotic dynamics induced by pairwise and higher-order interactions. (a-c) The blue and green dots represent the simulation results for initial conditions of $N_{25}(0)=0.2$ and $N_{25}(0)=0.20001$, respectively, with all other initial values set to $N_{i}(0)=0.2$ ($i=1,\cdots,32$). (a) Bifurcation diagram of the system (\ref{['GLV']}) as a function of the mean higher-order interaction strength $\mu_2$. (b-c) Representative chaotic dynamics at $\mu_2=-0.45$, consistent with the bifurcation structure in panel (a), shown via a two-dimensional phase space projection and the corresponding Poincaré map. (d) Sensitivity analysis of system (\ref{['GLV']}) under the same parameters as in panels (b-c). The divergence $\Delta N_{25}$ quantifies the sensitivity to initial conditions, indicative of chaos. All simulations were performed with $S = 32$ species. See SM Sec. \ref{['Simulation details']} for simulation details.
• Figure 3: Higher-order interactions shape species distribution patterns across ecological communities. (a) Comparison of species abundance distributions across communities. Hollow markers represent simulations with pairwise interactions only; solid markers incorporate both pairwise and HOIs. Observed data are from published empirical studies MichaelRoswell2021MichaelSCrossley2020Bird23EstradaA2001. (b) HOIs reproduce the characteristic S-shaped rank-abundance curves observed in empirical ecosystems HubbellSP2001FuhrmanJA2008MichaelSCrossley2020. Solid markers denote empirical data; hollow markers show corresponding simulation outcomes. (c) Direct visual comparison of species abundance distributions in a bird community. Empirical observations are sourced from Cornell Lab: Birds of the World Bird6Bird7Bird8Bird9Bird11Bird15Bird16Bird17Bird18Bird23Bird30Bird34Bird36Bird37Bird38Bird51. (d) Comparison of Shannon diversity indices and Kolmogorov-Smirnov (K-S) test p-values quantifies the similarity between observed and simulated distributions in (c). At the 0.05 significance threshold, none of the p-values indicate statistically significant differences. (a-d) All simulations were evaluated at time $t = 1.0 \times 10^{5}$. See SM Sec.\ref{['Simulation details']} for full details.
• Figure S1: The stability of the coexistent equilibrium. (a-b, d-e) Time courses and phase diagrams in the scenario involving pairwise and higher-order interactions, respectively. (c, f) The eigenvalues of the Jacobian of this community when the scenario involving pairwise and higher-order interactions is considered. (c) The coexistent equilibrium is stable when all the eigenvalues have a negative real part corresponding to \ref{['stability']}(a-b). (f) The red dotted displays the eigenvalues with a positive real part, indicating instabilitycorresponding to \ref{['stability']}(d-e). The simulations involve $S = 5$ species. In (a-f): $r_{i}=0.4$, $d_{i}=0.001$, $\alpha_{ij}=\mathcal{N}(-0.25, 0.1)$$\beta_{ijk}=\mathcal{N}(\mu_{2}, 0.3)$, $i, j , k= 1,\dots, 5$; In (a-c): $\mu_{2}=-0.25$; In (d-f): $\mu_{2}=-0.15$.
• Figure S2: The influence of stochasticity on species coexistence and stability. (a-b) Phase diagrams in the scenario involving pairwise and higher-order interactions, respectively. The species’ stability fraction in each pixel was calculated from 60 random repeats. (c-d) The coexisting fraction of species is introduced by higher-order interactions. The fraction of species that coexist under different mean strengths of higher-order interactions. In (a-b) the parameter values are provided in Fig. \ref{['pairwise higher-order']}. In (c-d): $r_{i}=0.3$, $d_{i}=0.001$, $\beta_{ijk}=\mathcal{N}(\mu_{2}, 0.3)$, $i, j , k= 1,\dots, 32$; In (c): $\alpha_{ij}=\mathcal{N}(\mu_{1}, 0.3)$; In (d): $\alpha_{ij}=\mathcal{N}(\mu_{1}, 0.5)$.