High-order species interactions shape ecosystem diversity

1A major challenge of theoretical ecology is explaining the stable coexistence of multi-species communities. While species interactions are important for maintaining natural diversity, their random nature destabilizes large communities. Even a community that satisfies the exclusion principle, or is otherwise inherently stabilized, can have its stability jeopardized by random interactions among species. This challenge was embodied in the seminal work by May, showing that the number of species that can stably coexist in randomly interacting communities is inversely proportional to the strength of their effective pairwise interactions. In ecological communities with strictly random pairwise interactions, this result predicts an upper threshold on diversity, as observed in random pairwise interaction models. Constrained properties of interactions, such as trophic level, intervality, broad degree distributions, positive and negative reciprocity, distribution skewness and connectance between differentially self-regulating species have been suggested as ways to explain the coexistence of large communities. Less is known about the feasibility, and possibly even the necessity, of diversity within the strict context of unconstrained random interactions.

2Species can exhibit interactions that inherently involve multiple species. While modelling in ecology often assumes species communities with pairwise interactions (Fig. 1a), species can also interact in higher-order combinations, that is, the interactions between two species can be modulated by other species. For instance, while a microbial species can produce an antibiotic that inhibits the growth of a competing species, this pairwise inhibitory effect can be attenuated by a third species that produces an enzyme that degrades the antibiotic. This third species thus modifies the interaction between the antibiotic-producing and sensitive species, without having a direct effect on any of them in isolation (Fig. 1b). The expression or activity of the antibiotic-degrading enzyme may in turn be inhibited by compounds produced by a fourth species, thereby generating a four-way interaction (Fig. 1c). Another general class of high-order interactions may arise when species exhibit adaptive behaviour, such as a predator switching its prey when a more preferred prey becomes available, or a prey that reacts to the presence of a predator by decreasing its foraging activity directed at a third species. Such behavioural effects can quickly lead to interactions at even higher order; for example, the third species may subsequently increase its activity against a fourth species, generating an effect that is inherently four-way. While the paradigm of pairwise interactions can capture the density-mediated effect of a predator on a species devoured by its prey that arises from the predator's impact on the prey's density, it misses the often larger trait-mediated effect of the predator on this third species through its impact on the prey's behaviour.

3While the importance of high-order interactions has been recognized, their general impact on critical community size has not been studied. One class of three-way interactions, whereby one species attenuates the negative interactions between two others, can stabilize well-mixed communities in specific network configurations. Three-way interactions have further been shown to promote diversity in patchy environments by increasing sensitivity to initial conditions. While increasing the order of interactions weakens their destabilizing effect and can increase the variance of species abundances at stationary states, it is unclear how high-order interactions affect the critical threshold on community size, and whether the classical result that random interactions impose a cap on community diversity still applies.

4Here, we combine simulations and theoretical analysis of random communities to examine the relation between diversity and stability in ecosystems with both pairwise and higher-order interactions. We find that high-order interactions impose a lower bound on the number of species, in addition to the upper bound imposed by pairwise interactions. Altogether, interactions of different orders therefore define an optimally stable range of diversities.

5We use a replicator dynamics model to simulate communities with interactions of different orders. Our model assumes N species with abundances xi, and whose per-capita growth rates are determined by the sum of pairwise interactions Aijxj (where Aij determines the effect of species j on species i), of three-way interactions Bijkxjxk (where species j and k have a joint effect on species i) and so on. This model therefore represents the Taylor expansion of a general system that includes trait-mediated effects as well as pairwise interactions:

6The model assumes random interactions perturbing an inherently stable community. The stability is provided by the first term (−xi), which reflects the assumption that species are self-limiting in high concentrations. The pairwise interactions matrix is given by, where is a random matrix whose elements are drawn from a Gaussian distribution with mean 0 and variance 1, and α determines the strength (variance) of the pairwise interactions. The three-way interactions are set by the tensor, where is a random three-dimensional (N × N × N) tensor whose elements are drawn from a Gaussian distribution with mean 0 and variance 1, and β represents the strength of three-way interactions; and similarly with four-way interactions whose strength scales with a parameter γ. As the total biomass is often determined by external resources, such as nutrients, energy and space, we keep the sum of all species abundances constant by subtracting from the individual growth rate of each species the average growth rate of the species in the community.

7As expected from the analysis by May and subsequent results, without high-order interactions we find that coexistence is lost when the pairwise interaction strength α is increased beyond a critical threshold that decreases with the number of species. We ran simulations with strictly pairwise interactions (α>0 with β=0, γ=0) for varying interaction strengths α and different numbers of species N. When the pairwise interaction strength α is small the simulations reach a steady state where all species coexist (Supplementary Fig. 1a), but when α increases simulations start exhibiting species extinctions (Supplementary Fig. 1b). Therefore, for each number of species N, we can define a critical threshold αc as the strength of pairwise interactions at which community feasibility (positive abundance for all species at equilibrium) is lost for a given fraction of the simulations (5%, Supplementary Fig. 1c,d). When the number of species N is increased, the critical threshold αc decreases as 1/N, in accordance with the result by May (Fig. 2a). Thus, as communities become more diverse they are more sensitive to random pairwise interactions (Fig. 2b).

8Strikingly, this classic relation between diversity and stability is inverted when high-order interactions are considered. We repeated the analysis above with strictly three-way interactions (only β>0) and strictly four-way interactions (only γ>0). Similar to the case of pairwise interactions, when the interaction strengths β or γ are increased beyond a certain level, βc or γc, feasibility is lost and extinctions abound. However, in contrast to pairwise interactions, the critical thresholds for high-order interactions do not decrease with the number of species. The critical threshold of three-way interactions βc remains unchanged, while the critical threshold of four-way interactions γc increases in proportion to the number of species (Fig. 2a). These four-way interaction communities can therefore withstand stronger interactions the more species they have. For a fixed strength of four-way interactions, these communities are stabilized, rather than destabilized, as the number species increases (Fig. 2c; increased stability with diversity can also appear in reciprocity models). Thus, while pairwise interactions put a higher bound on diversity, it is a lower bound that is created by high-order interactions. Provided that the total abundance is constant (does not change with the number of species), these results are robust to various changes in the underlying model: they hold for different distributions of the random coefficients (Supplementary Fig. 2), for connectance values smaller than 1 (Supplementary Fig. 3) and when high-order diagonal terms are set to zero (Cijk≠0 only when i≠j≠k≠i; Supplementary Fig. 4). Similar scalings of the critical strengths of the different interaction orders with the number of species also occur when the dynamic equations are replaced by a Lotka–Volterra model (Supplementary Fig. 5).

9To better understand the relation between the diversity of a community and its sensitivity to interactions of different orders, we examine a modified pairwise interactions matrix that captures the high-order interactions around a fixed point. Considering an approximation for how the interactions of different orders contribute to the effective pairwise interactions considered in the analysis by May (namely, the Jacobian), we derive a scaling relationship for the dependence of the feasible range of N on the different orders of interactions (see the ‘Methods' section, Supplementary Fig. 6):

10Intuitively, this condition tests whether the s.d. of the sum of all random interactions affecting a species near a fixed point is smaller than the self-stabilizing term (assuming fixed total abundance, see the ‘Methods' section).

11This theoretical analysis explains the different scaling of each of the interaction orders with N in accordance with our simulation results (Fig. 2a). When high-order interactions are turned off (only α>0), the critical strength of pairwise interactions αc decreases as 1/N. The critical threshold of exclusively three-way interactions βc does not depend on the number of species, while that of strictly four-way interactions γc is proportional to the number of species. Thus, interactions of any given order destabilize communities, and the effects of different orders are additive, but they scale differently with diversity. Small differences between these theoretical slopes and the numerical slopes of Fig. 2a might result from the small species numbers of the simulations, for which the analytic analysis holds only approximately.

12Since communities with a small number of species are sensitive to high-order interactions, while communities with a large number of species become sensitive to pairwise interactions, we asked whether the combination of pairwise and high-order interactions might dictate an intermediate range of species numbers that optimizes feasibility. To answer this question, we ran simulations with both pairwise and four-way interactions for three different community sizes N (for simplicity we set β=0; the effect of β>0 is discussed below). For each value of N, we scanned α and γ and identified the domain where communities are feasible (Fig. 3a). Since the critical value of α decreases with N, while the critical value of γ increases, there are regions where an intermediate-sized community would be feasible, while either smaller or larger communities would not (area within the N=8 domain but not within the N=5 or the N=18 domain). To understand how the feasibility of specific communities depends on the number of species, we ran simulations where we gradually add species with random interactions to communities with given values of interaction strengths α and γ: pairwise-dominated, four-way dominated and mixed. When pairwise interactions are dominant, feasibility decreases when N is increased beyond a threshold (Fig. 3b). In contrast, when four-way interactions are dominant, the community becomes sensitive to the removal of species (Fig. 3d). In the mixed case, when pairwise and high-order interactions are combined, feasibility is maximized for an intermediate number of species; both the removal and the addition of species destabilize the community (Fig. 3c).

13The combination of pairwise and high-order interactions thus defines a range of stable diversities. Assuming that a given natural ecosystem can be characterized by its average strengths of interactions of different orders (α, β, γ…), which are characteristics of its species and the environment, the upper and lower bounds on its diversity N can be derived as the roots of equation 2 (see the ‘Methods' section; Supplementary Fig. 7). As the total strength of interactions α+β+γ increases, these upper and lower bounds grow closer and the range of allowed diversities narrows around a defined number of species N* (Fig. 4 for β=0; using β=γ gives qualitatively similar results, Supplementary Fig. 8). This value of optimally stable diversity is determined by the relative strengths of the four-way and pairwise interactions (, see the ‘Methods' section). High-order dominated and strong interactions can thus require communities to maintain a large number of species to remain feasible.

14The combination of pairwise and high-order interaction strengths determines both upper and lower bounds on the diversity of multi-species communities. While very little is known about the abundance and strength of high-order interactions in nature, these results bring a new perspective to understand much empirical evidence showing that communities are sensitive not only to the addition of species but also to the opposite perturbation, namely the removal of species. Considering the perplexing effects of high-order interactions may allow a better understating of why high diversity is often a necessity rather than an option for natural ecosystems.

15To find the probability that a random community made up of N species with specified strengths of pairwise and high-order interactions is feasible, we numerically integrate equation 1 with different randomizations of the interaction matrices. For a given N, we generate a set of random pairwise, three- and four-way interaction matrices, and then for a given choice of α, β, γ we substitute into equation 1:, and. We then follow the species abundances as they change from a uniform initial value (xi=1/N; in the stable regime, the equations should not be sensitive to initial conditions) using the MATLAB ode45 solver with default integration properties. Simulations were ended when they reach a persistent state (steady-state or bounded oscillations). We repeat this process R times (∼300) with different sets of random pairwise, three- and four-way interaction matrices.

16We concluded that a community had reached a persistent state by demanding that the entropy (−Σixi log(xi)) remain constant or bounded in its fluctuations. We evolve the equations for a fixed number of time units (Δt=10), and after each of those intervals, we calculate the ratio between the minimal value taken by the species distribution entropy in two different periods of time: the last tenth and the last three-tenths of the simulation time length. If the ratio between those two values, as well as the similar ratio between the entropy maxima in those two periods, falls within a range of of 1, we concluded the simulation. If this condition was not met after a long time (104 time units) we concluded the simulation as well, and to account for the possibility that a species would have become extinct after this time threshold, we counted those simulations as unfeasible in the calculation of the lower error-bars.

17We defined a community as feasible if all the species existed at the end of the simulation, with an abundance above a defined threshold (set as 10−5/N). The probability that a community would be feasible was then defined as the fraction of stable communities out of the R random communities simulated. To find the critical strength of interactions for a given number of species, we used a fixed collection of R sets of pairwise, three- and four-way interaction matrices and increased the interaction parameters until 5% of the communities exhibited extinctions. In Fig. 2 (as well as Supplementary Figs 2–5) this was done by increasing the strength of interactions of a given order, while the rest are kept at 0. In Fig. 3a this was done along radial lines in log (α), log (γ) space: α=0.001·1.2r·cos θ, γ=0.05·1.3r·sin θ, such that for each value of θ, we increased r until stability was lost.

18The Jacobian of the system at a given point is:

19Where, is the Kronecker delta. At a fixed point, the first term vanishes. We may also neglect the terms in the last line which are derivatives of the normalization, as we will justify below. Assuming stability of all species coexisting, the value of species abundances at this fixed point must scale as xi≈1/N (because the total abundance is fixed to 1). Equation 1 can therefore be approximated as pairwise dynamics:

20with

21the fixed point will be stable if the eigenvalues of the Jacobian all have negative real parts. Within our assumption xi≈1/N, we may focus on the eigenvalues of Aeff, which differs from the Jacobian by a constant factor N that does not affect the signs of the eigenvalues (as previously reported, the actual Jacobian in a steady-state deviates somewhat from this simple analysis due to deviation of the fixed-point xi from 1/N; Supplementary Fig. 9).

22The entries of the random component of Aeff are sums of independent identically distributed Gaussians with 0 mean, and the variance of such a sum is additive. For simplicity, we assume that only elements with j>k>l are non-zero (such as in Supplementary Fig. 4), so the three-way interactions term averages over N-1 such elements and its variance is therefore ∼N times smaller, and similarly with the higher-order terms. The variance of the random component of Aeff is thus given by:

23by Girko's circular law, as N→∞, the eigenvalues of a random matrix M of order N × N will form a disk of radius around the origin on the complex plane. In a large enough community, and taking into account the stabilizing term, the eigenvalues of Aeff will therefore be approximately distributed in a disc around (−1, 0) (Supplementary Fig. 6), so that the community will be stable with high probability if the radius of the disc is smaller than 1, or:

24If we remove the assumption that only elements with j>k>l are non-zero, the high-order interaction strengths in the last equation are multiplied by constant factors, but the scaling is unaffected.

25This analysis does not depend on the specific distribution of the matrix elements, as long as their variance is fixed, and our simulations indeed show it holds for different distributions (Supplementary Fig. 2). Similar scaling results can be obtained relying on Gershgorin's circle theorem (see ref. for the case of pairwise interactions; generalization to high-order case follows from the same reasoning applied above). The normalization term should not have a significant impact on the instantaneous stability analysis, since by the same reasoning, its contribution to the variance of the Jacobian would be N times smaller than the interaction term, as it is itself an average of N interaction terms. However, this analysis does depend on the sum of species abundances being constant (and xi thereby scaling as 1/N); a different scaling of critical interaction strengths will appear when the total abundance scales with the number of species.

26Note that local stability analysis models exclusively pairwise interactions, since the equation is linearized around a point and the high-order interactions are therefore embedded in the coefficients of the effective pairwise interactions obtained in the linearization. Our analysis therefore extends the result by May in the sense that it suggests how the interaction strength in the Jacobian May analyses should scale with the number of species with respect to the strengths of interactions of different orders.

27The same scaling of diversity with interactions of different orders can also be obtained by demanding that the expected variance of the random interactions be weaker than the self-stabilizing term. The n-species interaction term is given by:

28with an n-dimensional interaction tensor whose variance is 1, and α(n) the strength of the n-species interactions. Since the value of species abundances near a fixed point with all species coexisting scales as xi≈1/N, and from the additivity of the variance of independent random variables, the standard deviation of such a term scales as. Comparing this value with the stabilizing term |−xi|≈1/N, we get the condition α(n)<Nn−3, in agreement with our result (α<1/N for pairwise interactions, β<1 for three-species interactions and so forth). The additivity of the variance allows us to extend this reasoning to obtain equation 2 when interactions of different orders are combined.

29Restricting interactions up to the fourth order, communities are stable if they satisfy the equation:

30demanding an equality and solving for N gives the lower and upper bounds for number of species:

31these two bounds narrow around a defined optimal number of species when the discriminant vanishes:, which gives.

32MATLAB code for the replicator and Lotka–Volterra models is available from the corresponding author on request.

33The simulation results are available from the corresponding author on request.

34How to cite this article: Bairey, E. et al. High-order species interactions shape ecosystem diversity. Nat. Commun. 7:12285 doi: 10.1038/ncomms12285 (2016).