Fingerprints of High-Dimensional Coexistence in Complex Ecosystems

The coexistence of many competing species in an ecological community is a long-standing theoretical and empirical puzzle. Classic approaches in ecology assume that species fitness and interactions in a given environment are mainly driven by a few essential species traits, and coexistence can be explained by trade-offs between these traits. The apparent diversity of species is then summarized by their positions (“ecological niches”) in a low-dimensional trait space. Yet, in a complex community, any particular set of traits and trade-offs is unlikely to encompass the full organization of the community. A diametrically opposite approach assumes that species interactions are disordered, i.e., essentially random, as might arise when many species traits combine in complex ways. This approach is appealing theoretically, and can lead to novel emergent phenomena, fundamentally different from the picture painted by low-dimensional theories. Nonetheless, fully disordered interactions are incompatible with many-species coexistence, and neither disorder nor its dynamical consequences have received direct empirical support so far. Here we ask what happens when random species interactions are minimally constrained by coexistence. We show theoretically that this leads to testable predictions. Species interactions remain highly disordered, yet with a “diffuse” statistical structure: interaction strengths are biased so that successful competitors subtly favor each other, and correlated so that competitors partition their impacts on other species. We provide strong empirical evidence for this pattern, in data from grassland biodiversity experiments that match our predictions quantitatively. This is a first-of-a-kind test of disorder on empirically measured interactions, and unique evidence that species interactions and coexistence emerge from an underlying high-dimensional space of ecological traits. Our findings provide a new null model for inferring interaction networks with minimal prior information and a set of empirical fingerprints that support a statistical physics-inspired approach of complex ecosystems.

Popular Summary

How biodiversity persists over many generations, without any single species consuming or outcompeting all the others, is an enduring puzzle in ecology. This puzzle suffers from a surfeit of solutions: Over a century, dozens of explanations have been proposed for how two or more species can coexist. Many of these mechanisms have then been backed up by experiments and observations in various ecological settings. Can we tell whether the coexistence of many species in an ecosystem is ensured by a few mechanisms or by many acting together? Here, we propose a way of demonstrating the latter.

When the coexistence of species is ruled by many interacting mechanisms—a scenario we call high-dimensional coexistence—then species interactions should appear as random as possible while still allowing every species to persist. We show that these minimal assumptions lead to a statistical structure in which species that are more successful within their ecosystem are less likely to compete strongly with each other. This structure comes with measurable statistical fingerprints, which we successfully identify in experimental plant communities.

While traditional approaches often suggest that each ecosystem may only be understood separately by characterizing its specific dominant mechanisms, high-dimensional coexistence offers a different outlook. We expect that the behavior of such ecosystems—their diversity, their dynamics, and their long-term persistence—will be more universal, in the way that macroscopic physical phenomena emerge from, and are insensitive to most details of, the underlying microscopic world.

Figure 1

Dimensionality of interaction patterns. The network of pairwise interactions between $S$ species can be represented by a $S\times S$ matrix (square box) where each element ${\beta }_{ij}$ denotes the effect of species $j$ on ${\eta }_i$, the relative yield of species $i$ defined in Eq. (1). We assume that interactions $\beta$ are constant in time, and entirely determined by underlying species traits and limiting factors $\theta$, for instance representing resources, pathogens, or behaviors. In a low-dimensional trait space, each species is characterized by a few traits; therefore the interaction matrix can be ordered by trait values and displays a simple organization. In a high-dimensional trait space, each interaction is a complex combination of factors, potentially unique and independent of other interactions. This can lead to a matrix without conspicuous order, often modeled as a random matrix [11, 12, 13, 14].

Figure 2

Three approaches to coexistence. (a) An equilibrium coexistence mechanism can be expressed as a specific pattern in the network of species interactions. To allow coexistence, i.e., an equilibrium where all species have positive relative yields $\eta$ as defined in Eq. (1), the traits $\theta$ which determine species interactions $\beta (\theta )$ must follow particular relations or constraints. We show here the classic low-dimensional example of a trade-off between two traits (competitive ability and colonization rate [9, 10]; see Supplemental Material [24]). (b) Our inferential approach reverses the direction of reasoning: given that we observe species coexisting at equilibrium in nature, we can infer the most likely distribution of interactions that could have caused this equilibrium, i.e., the posterior distribution $P(\mathit{\beta }|\overrightarrow{\eta })$ given a prior $P(\mathit{\beta })$ over all possible interaction networks. If we choose a high-dimensional prior (independent interactions), communities drawn from the posterior distribution retain random-looking but correlated interactions, with simple and predictable statistical features: a “diffuse clique” pattern described below. (c) Community assembly (see Sec. 3) involves a larger pool of $S_{\mathrm{pool}}>S$ species, here with independent interactions. Through ecological dynamics, some of these species go extinct or cannot invade, leading to a smaller persistent group of $S$ species whose interactions $\beta$ have been dynamically selected to allow coexistence. Dotted arrows: The assembled community in (c) precisely follows the statistical properties of communities drawn from the posterior distribution in (b), if we assume independent interactions for both the inference prior and pool distribution [25]. Low-dimensional structures, such as illustrated here in (a), can be seen as possible but very improbable matrices in this distribution, and can exhibit different patterns from our diffuse clique prediction. We can interpret the posterior distribution $P(\mathit{\beta }|\overrightarrow{\eta })$ as a null model for coexistence in the presence of many biological mechanisms. The resulting diffuse clique pattern represents a form of collective organization, where coexistence arises, not from particular species traits, but from statistical biases distributed over all interactions.

Figure 3

The diffuse clique structure is characterized by two statistical patterns in how successful species compete. Any typical interaction network drawn from the posterior distribution $P(\mathit{\beta }|\overrightarrow{\eta })$ (see Fig. 2) will have similar statistical features. (a) Trend in the expectation $E[{\beta }_{ij}|\overrightarrow{\eta }]$ of competition strength Eq. (4). If all species competed with equal strength $\overline{\beta }$, we would expect any given species $i$ to achieve the relative yield ${\eta }_i={\eta }^{*}$ [Eq. (3)]. Given how much ${\eta }_i$ differs from this baseline, we can infer how interactions most likely deviate from $\overline{\beta }$. We show this deviation $\mathrm{\Delta }({\eta }_i,{\eta }_j)$ for simulated data, in matrix form: each element shows the bias in the interaction effect of species $j$ (column) on species $i$ (row). Left to right: An unsuccessful species (low ${\eta }_j$) competes indiscriminately against others (white, $\mathrm{\Delta }=0$), whereas a successful species (high ${\eta }_j$) is a biased competitor. Top to bottom: Species with ${\eta }_i<{\eta }^{*}$ experience stronger competition from successful species (red, $\mathrm{\Delta }>0$), whereas species with ${\eta }_i>{\eta }^{*}$ experience weaker competition from them (blue, $\mathrm{\Delta }<0$). Together, these biases indicate the existence of a “clique” of species that compete less against each other, and more against all others, thus achieving higher relative yield than the baseline ${\eta }^{*}$. (b) Correlation structure (5) between columns of the interaction matrix: the effects of two successful species are anticorrelated, less likely to target the same species, whereas competition from unsuccessful species is indiscriminate.

Figure 4

Validation of the Lotka-Volterra model. The consistency of the equilibrium description (2) was validated through a series of tests. For the Wageningen grassland experiment (with $S=8$ species), we find accurate predictions for the full fitted Lotka-Volterra model (dashed lines show a $1:1$ relationship between prediction and measurement). (a) We compare the equilibrium values of ${\eta }_i=B_i/K_i$ in the full-diversity plot ($S=8$) to the prediction ${\eta }_i=1-{\sum }_j{\beta }_{ij}{\eta }_j$, where interactions ${\beta }_{ij}$ were inferred from all partial compositions (plots with $S<8$). (b) Carrying capacities $K_i$ inferred as the intercept of the multilinear regression from polyculture plots ($S>1$) agree with direct measurements in monocultures ($S=1$). (c),(d) Even when species are too numerous to allow a precise inference of all pairwise interactions, we can predict statistics of ${\eta }_i$ using only the mean and variance of interactions ${\beta }_{ij}$ in tested species pairs, following the methodology in Ref. [23]. Across the Wageningen, Cedar Creek Big Bio (16 species), and Jena (60 species) experiments, we find good agreement in predicted versus measured values for the mean in (c) and coefficient of variation in (d) of relative yields in the full-diversity plots.

Figure 5

Diffuse structure in the Wageningen grassland experiment. (a) Using many species combinations (56 sets of $S<8$ species), we can fit interaction coefficients ${\beta }_{ij}$ in Eq. (2) by multilinear regression and compute their mean $E[\beta ]\approx 0.4$. (b) Using the relative yields ${\eta }_i$ in the full community ($S=8$ species), we predict the theoretical pattern of expectations (4) for each interaction [Fig. 3; here we see only negative biases in blue, indicating that all interactions are less competitive than the prior mean $\overline{\beta }$]. Abbreviations for species names are given in Sec. 4. (c)–(e) In our diffuse pattern, individual coefficients are expected to exhibit a large spread around their expectation $E[{\beta }_{ij}|{\eta }_i,{\eta }_j]$, and must be binned to compute empirical statistics. (c) We group coefficients ${\beta }_{ij}$ by their associated ${\eta }_i$, ${\eta }_j$ to compute the running average $A[{\beta }_{ij}|{\eta }_i,{\eta }_j]$ (gray curve, $\pm 1$ standard deviation in shaded area). Its proximity to the dashed $1:1$ line is one of multiple metrics of theory-data agreement. (d) We show the empirical running average $A[{\beta }_{ij}|{\eta }_i,{\eta }_j]$ in matrix form (median value for each species pair $i$, $j$), to be compared with the theoretical prediction $E[{\beta }_{ij}|{\eta }_i,{\eta }_j]$ in (b). This is a particular way of smoothing the empirical matrix ${\beta }_{ij}$ shown in (a); see Sec. 4. (e) In another test of the theory, we group coefficients by species, i.e., either by row or by column in the matrix ${\beta }_{ij}$. We then compute how coefficients in row $i$ (column $j$) vary with ${\eta }_j$ (${\eta }_j$). Each solid point is the linear regression slope within a row or column. The empirical values and their spread around the trend are both in good agreement with our theoretical predictions (mean trend in solid line, spread in dashed lines, $\pm 1$ SD). The linear slope of the points matches quantitatively the theoretical trend in the mean, which is approximately linear as expected for large $S$.

Figure 6

Cross-experiment validation of the diffuse clique structure. Multiple quantitative tests are needed to confirm agreement between data and theoretical predictions while ruling out counterexamples. Therefore, we compute a series of metrics (defined in the Supplemental Material [24]) of agreement for conditional expectations $E[{\beta }_{ij}|{\eta }_i,{\eta }_j]$, correlations $\mathrm{corr}({\beta }_{ij},{\beta }_{ik}|{\eta }_i,{\eta }_j,{\eta }_k)$, and diffuseness (whether interactions appear random). The first two metrics test the pattern shown in Fig. 3. We do not show each metric directly, because its predicted range is specific to each experiment. Instead, we assign a similarity score $s=+1$ if predicted and observed values differ by less than 2 standard deviations, $-1$ if they differ by more. We also compute a discrimination score $D$ between 0 and 1 measuring how well the metric differentiates between our theory and a counterexample (see Sec. 4). The product $s\times D$ is shown here. Other metrics were also tested, and the last column reports the number (or range across treatments or simulation runs) of positive results out of all 13 tests shown in the Supplemental Material [24]. The first and last row display simulated communities that constitute either counterexamples (low-dimensional trade-off patterns) or examples of our theory.

Figure 7

Spectrum of conditioned $\mathit{\beta }$. Crosses: eigenvalues of a matrix $\mathit{\beta }$ conditioned with the pattern defined by Eqs. (4) and (5), with $S=50$. Dashed line: circle within which ${\mathit{\beta }}^{\mathrm{cond}}$ should have eigenvalues due to $\sigma \mathbf{A}$. Small circles: predicted positions of additional two eigenvalues. Other parameters: ${\eta }_{1\dots 10}=3$, ${\eta }_{11..50}$ are i.i.d variables, uniformly sampled in [0, 1.5]; $\overline{\beta }=0.05$; $\sigma =0.5$.