Feasibility and Single Parameter Scaling of Extinctions in Multispecies Lotka-Volterra Ecosystems

TL;DR

The paper develops a random-matrix-theory framework to analyze feasibility and extinctions in multispecies Lotka-Volterra ecosystems in the weak-interaction limit. By deriving the ensemble-averaged abundance $\overline N$, its variance $\Sigma_N$, and higher moments via Neumann-series expansions, it shows that fixed-point abundances are Gaussian and that feasibility breaks before stability as the species number grows. It further connects the distribution of negative fixed-point abundances to extinction events under LV dynamics and uncovers a single-parameter scaling law governed by $\eta = \overline N/\sqrt{2\Sigma_N}$ for the number of extinctions. The results hold across general cross-diagonal correlation $\gamma$ and offer predictive, statistically grounded insights into ecosystem coexistence and collapse, with strong numerical validation. Key contributions include analytic extinction probabilities and demonstrated data-collapse across diverse parameter regimes.

Abstract

Multispecies ecosystems modelled by generalized Lotka-Volterra equations exhibit stationary population abundances, where large number of species often coexist. Understanding the precise conditions under which this is at all feasible and what triggers species extinctions is a key, outstanding problem in theoretical ecology. Using standard methods of random matrix theory, I show that distributions of species abundances are Gaussian at equilibrium, in the weakly interacting regime. One consequence is that feasibility is generically broken before stability, for large enough number of species. I further derive an analytic expression for the probability that $n=0,1,2,...$ species go extinct and conjecture that a single-parameter scaling law governs species extinctions. These results are corroborated by numerical simulations in a wide range of system parameters.

Figures (6)

• Figure 1: Average (left panel) and variance (right) of the abundance distribution at stable equilibria for Eq. \ref{['lv']}, $\mu=0$, and a random, homogeneous distribution $k_i \in [0.5,3.5]$, i.e. with $\kappa=2$ and $\chi=0.75$, as a function of $\sigma$ and for $\gamma=-0.5$ (green), 0 (red) and 0.5 (black). Solid lines give the RMT predictions of Eqs. \ref{['avgfins']} and \ref{['sigmu']}. Data correspond to averages over all species when the equilibrium of Eq. \ref{['lv']} has been reached for 500 realizations of the random interaction matrix $\mathbb A$ with $S=37$ ($\times$) and 157 (circles).
• Figure 2: Normal (a) and semilog (b) plots of species abundance distributions for $\mu=0$, $k_i \equiv1$, and $\sigma=0.2$, $\gamma=-0.5$ (black circles), $\sigma=0.4$, $\gamma=-0.5$ (red), $\sigma=0.5$, $\gamma=0$ (blue). Solid lines give Gaussian distributions with average ($\overline N = 0.98$, 0.93 and 1.01) and variances ($\Sigma_N=0.038$, 0.138, 0.31) given by Eqs. \ref{['avgfins']} and \ref{['sigmu']}. Distributions are calculated from fixed point solutions to Eq. \ref{['lv']} (dots) and from solutions to Eq. \ref{['equi2']} ($\times$), over 1000 realizations of the interaction matrix. The blue (red) dot at $N_i^*=0$ correspond to $\simeq 6.5$ (0.8) extinctions per realization, in agreement with Eq. \ref{['probNe']}. Inset: Abundances $N_i$ at a fixed point solution to Eq. \ref{['lv']} (black) and for Eq. \ref{['equi2']} (red) for the same realization of $\mathbb A$, for $\mu=0$, $k_i \equiv1$, $\sigma=0.5$ and $\gamma=0$. There are six species with negative abundances at the fixed point and six species with abundances below 10$^{-20}$ and still going down after 2 $\times 10^{7}$ Runge-Kutta iterations of Eq. \ref{['equi2']}.
• Figure 3: Distribution of number $N_e$ of extinct species, from the species with negative abundancies in Eq. \ref{['equi2']} (+) and from Eq. \ref{['probNe']} (dots), with $k_i \equiv1$. (a) $S=607$, $\mu=0$, $\sigma=0.35$ and $\gamma=-0.5$ (blue), 0 (black) and 0.5 (red); (b) $S=607$, $\sigma=0.35$, $\gamma=0$ and $\mu=0$ (black), 2 (red) and 4 (blue); (c) $S=607$, $\mu=0$, $\gamma=0$ and $\sigma=0.23$ (blue), 0.29 (red) and 0.35 (black); (d) $\mu=0$, $\gamma=0$, $\sigma=0.35$ and $S=57$ (green), 157 (blue), 307 (red) and 607 (black). All distributions are calculated over 1000 different realizations of the interaction matrix $\mathbb A$.
• Figure 4: Normal (a) and semilog (b) plots of the distribution of number $N_e$ of extinct species, for $S=157$, $\mu=0$, $k_i \equiv1$, $\gamma=0$, $\sigma=0.35$ (red symbols) and 0.4 (black). Distributions are calculated from fixed point solutions obtained by time-evolving Eq. \ref{['lv']} (circles) and from solutions to Eq. \ref{['equi2']} (+), over 1000 different realizations of the interaction matrix $\mathbb A$. (c) Scaling of the ratio of the average number of extinct species vs. ratio of the average over standard deviation of the abundance distribution, for 24 different sets of parameters, $\sigma \in [0,0.6]$, $\mu \in [0,1]$, $\gamma \in [-1,1]$, $S \in [100,300]$, and $k_i \equiv 1$ (circles) as well as $k_i \in [0.5,3.5]$ (crosses). Averages are calculated over 100 to 500 realizations of the random matrix $\mathbb A$ for each set of parameter, by time-evolving Eq. \ref{['lv']} until a stationary solution is reached. Symbols of the same color correspond to sets varying either $\sigma$ or $\gamma$, all other parameters being fixed. When time-evolving Eq. \ref{['lv']}, extinctions are defined at a threshold $N_c=10^{-20}$.
• Figure S1: Convergence of the series expansion for $\mu=0$ and $k_i \equiv \kappa = 1$. The solid curves give Eq. \ref{['avg1fin']} truncated at the zeroth (red), second (green), fourth (violet) and sixth (black) order in $\sigma$.