Introduction

In diverse ecological communities, species interact with other species in the community temporally (Li and Chesson 2016) and/or spatially (Hart et al. 2017). Understanding the mechanisms behind maintenance of species diversity has been one of the central goals of ecological research for decades (Levine et al. 2017). While our primary understanding of species coexistence empirically and theoretically has mostly come from study of species pairs (Chesson 2000; Barabás et al. 2012; Kraft et al. 2015; Hart et al. 2016), understanding the dynamics and coexistence of multiple species from the viewpoint of species pairs becomes unfeasible as species richness increases (Barabás et al. 2016; Bairey et al. 2016).

In competition models, the underlying processes that facilitate species coexistence require parameter trade-offs in competitive interactions to stabilize multispecies coexistence (Barabás et al. 2016). Moreover, competitive interactions between species in a diverse community are difficult to structure using generalizations from simple coexistence rules. In species-rich communities, underlying processes that cannot be captured by pairwise species interactions can emerge, and have been suggested to promote species coexistence (Abrams 1983). Such underlying processes, for example can be intransitive or “rock-paper-scissors” interactions (Laird and Schamp 2006; Gallien et al. 2017; Saavedra et al. 2017). Intransitive interactions are inherently pairwise in nature but they form interaction chains that favour species coexistence. In such three species systems, there is also the possibility of interactions among species to be not pairwise in nature. For instance, instead of pairwise interactions, it is also possible that a species could modulate pairwise interactions between two other species leading to interactions that are of higher order in nature (Levine et al. 2017). Advancing our understanding of species coexistence in diverse communities requires an understanding of mechanisms that could inherently be high-dimensional.

When diverse species communities are modelled explicitly by considering only pairwise interactions, multiple species coexistence is not always possible (Bairey et al. 2016; Levine et al. 2017). Recent studies have suggested that in order for large number of species to coexist, interactions between species should not be constrained to species pairs but should include higher order interactions (HOIs) (Grilli et al. 2017). Although such HOIs have been suggested to stabilize large ecological communities, their underlying role in maintaining species coexistence, from the purview of modern coexistence theory (MCT), is unknown. This is primarily due to the difficulty of integrating MCT for more than two species simultaneously (Saavedra et al. 2017). MCT states that coexistence is possible when fitness differences between species are smaller than their niche differences (Chesson 2000). In MCT, coexistence of species can be understood from a mutual invasibility criteria, where the invasion growth rate of a species is analytically decomposed into stabilizing niche differences and fitness differences (Kremer and Klausmeier 2013; Gallien et al. 2017). Stabilizing niche differences increases the probability for species coexistence, while fitness differences increase the probability of competitive exclusion. Importantly, fitness and niche differences can be quantified from the terms of the Lotka-Volterra pairwise competition model as (Chesson 2000, 2018; Saavedra et al. 2017):

$$ \underset{\mathrm{Niche}\ \mathrm{Overla}{\mathrm{p}}^{-1}}{\underset{\hbox{---} }{\sqrt{\frac{\alpha_{11}{\alpha}_{22}}{\alpha_{12}{\alpha}_{21}}}}}>\underset{\mathrm{Fitness}\ \mathrm{difference}}{\underset{\kern1.25em \hbox{---} }{\frac{r_1}{r_2}\sqrt{\frac{\alpha_{22}{\alpha}_{21}}{\alpha_{11}{\alpha}_{12}}}}}>\underset{\mathrm{Niche}\ \mathrm{Overla}\mathrm{p}}{\underset{\hbox{---} }{\sqrt{\frac{\alpha_{12}{\alpha}_{21}}{\alpha_{11}{\alpha}_{22}}}}} $$
(1)

Here αij is the pairwise competition coefficients between species i and j derived from a generalized Lotka-Volterra model, and ri is the intrinsic growth rate of a species i. If niche overlap is greater than the fitness difference between two species, then coexistence is not possible. Under certain simplifying assumptions, one can integrate results from HOIs into the traditional framework of pairwise species coexistence. Such an integration would make the relevance and the understanding of HOIs more complete.

Here, using a three-species Lotka-Volterra model, we demonstrate the importance of HOIs in maintaining and disrupting species coexistence. Specifically, using the invasibility criterion, we modified pairwise-interspecific and pairwise-intraspecific coefficients of the Lotka-Volterra model in a way that allowed us to create a range of fitness differences ranging from low to high. We then show how negative three-way HOIs (HOIs that intensify pairwise competitions) and positive three-way HOIs (HOIs that alleviate pairwise competition) can stabilize species coexistence in fitness regions, where species coexistence is impossible if only pairwise interactions are considered. In addition, we also evaluate how four-way HOIs could impact species coexistence under regions of extreme fitness difference. We then extend our three-species HOIs case, to a multispecies competitive community, and show that the conditions under which HOIs stabilize species coexistence in the three-species case still holds in a multispecies community. We thus highlight the possible mechanisms by which HOIs could promote coexistence in species-rich communities.

Methods

Higher order interactions

Lotka-Volterra model for three-way HOIs (Wilson 1992; Billick and Case 1994; Levine et al. 2017) can be written as:

$$ \frac{d{N}_i}{dt}={N}_i{r}_i\left(1-{\sum}_j^n{\alpha}_{ij}{N}_j+{\sum}_j^n{\sum}_k^n{\beta}_{ij k}{N}_j{N}_k\right) $$
(2)

and four-way HOIs can be written as:

$$ \frac{d{N}_i}{dt}={N}_i{r}_i\left(1-{\sum}_j^n{\alpha}_{ij}{N}_j+{\sum}_j^n{\sum}_k^n{\beta}_{ij k}{N}_j{N}_k+\kern0.5em {\sum}_j^n{\sum}_k^n{\sum}_l^n{\gamma}_{ij k l}{N}_j{N}_k{N}_l\right) $$
(2.1)

where αij and βijk are pairwise interactions and three-way HOIs, respectively, and γijkl represents four-way HOIs. Here, higher order terms could broadly be defined as non-additive effects on per capita fitness of a species. In case of three-way HOIs, e.g. in βijk, the per capita effect of species j on i depends on the density of the third species (k)—this represents a non-additive effect, because the fitness of i cannot simply be predicted by the densities and per capita interaction strengths of the other two species. HOIs could intensify or alleviate the pairwise competition between two species depending on the sign of βijk as negative or positive, respectively. In the case of four-way HOIs, γijkl captures the non-additive effect of species l on a three-way interaction of species k and j on species i. We, first, evaluate the effect of three-way and four-way HOIs on a three species (n = 3) community, and later on a larger multispecies (n = 50) communities (see “Positive three-way and four-way higher order interactions on species coexistence”). Note that in Eqs. 2 and 2.1, we have a positive sign in front of three-way and four-way HOIs, while Letten and Stouffer (2019) and Kleinhesselink et al. (2019) have a negative sign; however, the interpretations are similar. The study by Kleinhesselink et al. (2019) specifically focuses on competitive HOIs (negative HOIs), and touches partly on positive HOIs. However, in our study, we focus on HOIs that could be either positive or negative in nature, such that positive HOIs (βijk, γijkl) will have positive values and negative HOIs (−βijk, −γijkl) will have negative values. In addition to this, Letten and Stouffer’s (2019) definition of HOIs differs from our formulation. In their interpretation of HOIs, for instance, βijk would mean a non-additive effect of both the pairwise effect of species k on species i and species j on species i. We, however, in our study follow the earlier definition of HOIs from Wilson (1992), Billick and Case (1994), and Levine et al. (2017), where HOIs are defined as “a change in the interaction of two species caused by a third species” (Billick and Case 1994). In addition, in our study, “the third” species that modifies the pairwise interaction could also be one of the two species involved in that pairwise interaction.

In our three species model, we make the following assumptions while analytically deriving the invasion growth rate:

  1. 1)

    There is interspecific competitive interaction between species 1 and 2, but not with species 3. This means that the matrix of competitive interactions will be as follows:

    $$ \left(\begin{array}{ccc}{\alpha}_{11}& {\alpha}_{12}& {\alpha}_{13}\\ {}{\alpha}_{21}& {\alpha}_{22}& {\alpha}_{23}\\ {}{\alpha}_{31}& {\alpha}_{32}& {\alpha}_{33}\end{array}\right)=\left(\begin{array}{ccc}{\alpha}_{11}& {\alpha}_{12}& 0\\ {}{\alpha}_{21}& {\alpha}_{22}& 0\\ {}0& 0& {\alpha}_{33}\end{array}\right) $$
    (3)
  2. 2)

    Only interspecific HOIs are taken into account. This means that terms such as βiii = γiiii = 0, where i = 1, 2, 3.

  3. 3)

    Species 3 influences species 1, and species 2 only through HOIs. Also, species 3 does not get influenced by species 1 or 2 through HOIs, i.e. β3jk = 0; where j, k ϵ [1, 2] and γ3jkl = 0, where j, k, l ϵ [1, 2, 3].

These assumptions were made to ensure that the number of terms in the HOI model is tractable for analysis (see Appendix 1). Using these assumptions, we can expand the three-way HOI model in Eq. (2) as:

$$ \frac{d{N}_i}{dt}={N}_i{r}_i\left(1-{\alpha}_{i1}{N}_1-{\alpha}_{i2}{N}_2+{\beta}_{i11}{N}_1^2+{\beta}_{i12}{N}_1{N}_2+{\beta}_{i13}{N}_1{N}_3+{\beta}_{i21}{N}_2{N}_1+{\beta}_{i22}{N}_2^2+{\beta}_{i23}{N}_2{N}_3\right) $$

where i = 1, 2 (for species 3, see Appendix 1). We can then write the HOI matrix for each of the species as:

$$ \left(\begin{array}{ccc}{\beta}_{i11}& {\beta}_{i12}& {\beta}_{i13}\\ {}{\beta}_{i21}& {\beta}_{i22}& {\beta}_{i23}\\ {}{\beta}_{i31}& {\beta}_{i32}& {\beta}_{i33}\end{array}\right) $$
(4)

where i = 1, 2; βi31 = βi32 = βi33 = 0; and βiii = 0. For instance, for species 1, the three-way HOI matrix would be:

$$ \left(\begin{array}{ccc}{\beta}_{i11}& {\beta}_{i12}& {\beta}_{i13}\\ {}{\beta}_{i21}& {\beta}_{i22}& {\beta}_{i23}\\ {}{\beta}_{i31}& {\beta}_{i32}& {\beta}_{i33}\end{array}\right)=\left(\begin{array}{ccc}0& {\beta}_{112}& {\beta}_{113}\\ {}{\beta}_{121}& {\beta}_{122}& {\beta}_{123}\\ {}0& 0& 0\end{array}\right) $$

Similarly, we can expand on the four-way HOI model in Eq. 2.1 as (for details, see Appendix 2):

$$ \frac{d{N}_i}{dt}={N}_i{r}_i\left(1-{\alpha}_{i1}{N}_1-{\alpha}_{i2}{N}_2+{\beta}_{i11}{N}_1^2+{\beta}_{i12}{N}_1{N}_2+{\beta}_{i13}{N}_1{N}_3+{\beta}_{i21}{N}_2{N}_1+{\beta}_{i22}{N}_2^2+{\beta}_{i23}{N}_2{N}_3+{\gamma}_{i121}{N}_1^2{N}_2+{\gamma}_{i122}{N}_1{N}_2^2+{\gamma}_{i123}{N}_1{N}_2{N}_3+{\gamma}_{i131}{N}_1^2{N}_3+{\gamma}_{i132}{N}_1{N}_3{N}_2+{\gamma}_{i133}{N}_1{N}_3^2+{\gamma}_{i211}{N}_2{N}_1{N}_1+{\gamma}_{i212}{N}_2^2{N}_1+{\gamma}_{i213}{N}_2{N}_1{N}_3+{\gamma}_{i221}{N}_2^2{N}_1+{\gamma}_{i222}{N}_2^3+{\gamma}_{i223}{N}_2^2{N}_3+{\gamma}_{i231}{N}_2{N}_3{N}_1+{\gamma}_{i232}{N}_2^2{N}_3+{\gamma}_{i233}{N}_2{N}_3^2\right) $$

We consider cases of both positive and negative HOIs, such that in positive HOIs, all β and γ are positive, and in negative HOIs, β and γ are negative, while calculating the invasion growth rates. We would like to state that these equations are representative of only competitive communities as the pairwise competition coefficients are purely negative for both species 1 and species 2.

Invasion growth rate and coexistence theory

The invasion growth rate, ri, is the per capita rate of increase in a species’ abundance—when it is rare—in presence of the other species, which is at equilibrium in the community. This means \( \frac{1}{N_i}\frac{d{N}_i}{dt}>0 \) for species i in the community. We here consider three species, i.e. i = 1, 2, and 3. Invasion growth rates of species 1 in the presence of three-way HOIs in the community can be written as (see Appendix 1):

$$ {r}_1^{\ast }=1-\left(\frac{\alpha_{12}\left({\alpha}_{22}-{\beta}_{223}{N}_3^{\ast}\right)-{\beta}_{122}-{\beta}_{123}\left({\alpha}_{22}-{\beta}_{223}{N}_3^{\ast}\right)}{{\left(\ {\alpha}_{22}-{\beta}_{223}{N}_3^{\ast}\right)}^2}\right) $$
(5)

where \( {N}_3^{\ast }=\frac{1}{\alpha_{33}} \) (Appendix 1).

Similarly, invasion growth rate of species 1 in the presence of four-way HOIs can be written as (see Appendix 2):

$$ {r}_1^{\ast }=1-\kern0.5em {\alpha}_{12}{N}_2+{\beta}_{122}{N}_2^2+{\beta}_{123}{N}_2{N}_3+{\gamma}_{1222}{N}_2^3+{\gamma}_{1223}{N}_2^2{N}_3+{\gamma}_{1232}{N}_2^2{N}_3+{\gamma}_{1233}{N}_2{N}_3^2. $$

We evaluated the effect of 3-way and 4-way HOIs in promoting coexistence by comparing invasion growth rates of species 1 in the presence and absence of HOIs. Following Gallien et al. (2017), we created scenarios where pairwise competitive matrix (3) varied from purely symmetric pairwise interactions to asymmetric interactions with gradually increasing pairwise fitness differences. As pairwise fitness difference (calculated from Eq. 1) increases, niche difference (calculated from Eq. 1) should also increase correspondingly, to stabilize species coexistence. This was done by modifying the pairwise interaction matrix to (Gallien et al. 2017):

$$ \left(\begin{array}{ccc}{\alpha}_{11}& {\alpha}_{12}& {\alpha}_{13}\\ {}{\alpha}_{21}& {\alpha}_{22}& {\alpha}_{23}\\ {}{\alpha}_{31}& {\alpha}_{32}& {\alpha}_{33}\end{array}\right)=x\ \left(\begin{array}{ccc}2\sqrt{\theta }& \theta & 0\\ {}1& 2\sqrt{\theta }& 0\\ {}0& 0& {\alpha}_{33}\end{array}\right) $$
(6)

where x = 0.01. The above modified matrix ensures that niche overlap (calculated using Eq. 1) between species 1 and species 2 is fixed at 0.5 even when fitness differences, controlled by θ, increases linearly. As θ is varied from 0 to 7 (θ = [0, 7]), fitness difference between species 2 and 1 increases linearly. Note that such increases in θ affected only pairwise competition coefficients, i.e. αij. Also to be noted is that the fitness difference between species 2 and species 1 is given from Eq. 1 as \( \sqrt{\frac{\alpha_{11}{\alpha}_{12}}{\alpha_{22}{\alpha}_{21}}} \), and for species 1 over species 2 is given as\( \sqrt{\frac{\alpha_{22}{\alpha}_{21}}{\alpha_{11}{\alpha}_{12}}} \). For certain θ values (θ = [5, 7]), fitness difference exceeds niche difference. Consequently, following the pairwise coexistence rule, species 1 can never invade species 2. This means that pairwise coexistence is impossible for these θ values. Here, θ = 5 defines the boundary between stable pairwise coexistence and competitive exclusion, when there are no HOIs. Note that species 3 remains unaffected by this modification and only participates in HOIs. Also note that Eq. 1 defines the conditions for species coexistence given by pairwise fitness difference and niche overlap for pairwise interactions (i.e. in the absence of any HOIs). However, fitness difference and niche overlap in the presence of HOIs can be calculated differently (see Appendix 3).

Negative three-way and four-way higher order interactions on species coexistence

Next, to evaluate whether the presence of negative three-way HOI stabilizes pairwise coexistence in scenarios where fitness differences are extreme, we estimated invasion growth rates of species 1 (\( {r}_1^{\ast } \)) when species 2 is present at equilibrium numerically and in the presence of HOIs, and compared with the invasion growth rate of species 1 (\( {r}_1^{\ast } \)) when HOIs are zero, i.e. in pairwise interaction scenario. This was done by varying fitness difference (Eq. 6) by increasing θ values from 0 to 7, such that fitness difference increased from 0 to 3. For each θ value, we calculated invasion growth rate for species 1 \( \Big({r}_1^{\ast } \)) and species 2 \( \Big({r}_2^{\ast } \)) in the presence of HOI, and in the absence of HOI (pairwise interaction case). When HOIs are absent, it is expected that given niche overlap is at 0.5, and pairwise coexistence becomes impossible as fitness differences increase. The importance of HOIs will be evident if species 1 (or species 2) could increase its invasion growth rates and invade even when fitness differences are large, i.e. they could increase in abundance even when differences between the two species in terms of fitness are large. For sensitivity analysis of invasion growth rate to HOIs, see Appendix 4.

Similarly, for the effect of negative four-way HOI on species coexistence, we estimated the invasion growth rates of species 1 and species 2 by varying pairwise fitness difference (Eq. 6) done by increasing θ values such that fitness difference increased from 0 to 3. We then compared the invasion growth rate of species 1 (and species 2) in the presence of four-way HOI with the pairwise interaction case (i.e. without the presence of any HOI).

Positive three-way and four-way higher order interactions on species coexistence

Next, we evaluated whether the presence of positive three-way and positive four-way HOIs stabilized pairwise species coexistence by calculating the invasion growth rates of species 1 and species 2 for a range of fitness differences. We created a range of pairwise fitness difference by varying θ in Eq. 6 (similar to above) and evaluated whether pairwise species coexistence is stabilized when positive three-way and positive four-way HOIs are introduced. We then compared the invasion growth rates of species 1 and species 2 in the presence of positive three-way and positive four-way HOIs and in the absence of HOI (i.e. pairwise interaction scenario).

Multispecies coexistence and higher order interactions

Multispecies generalization of two-species pairwise coexistence rule, however, is complicated. In multispecies communities, while all species pairs must satisfy the pairwise coexistence rule, this does not, however, guarantee stability (Barabás et al. 2016). For multispecies competitive communities, the stability and feasibility of coexistence can be evaluated from Weyl’s inequality (Fulton 2000), as described below. A pairwise competitive community can be written as:

$$ \frac{d{N}_i}{dt}={N}_i{r}_i\left(1-\sum \limits_j^n{\alpha}_{ij}{N}_j\right) $$

where, αij represents pairwise competitive interactions. αij is the element in the ith row and jth column of a matrix of pairwise competitive interactions. This pairwise matrix of competitive interactions can be denoted, by, say A. Here, we consider A as a symmetric matrix, such that, αij = αji.

Now, A can be decomposed into inter- and intraspecific matrices, say B and C, respectively, where B is a matrix of only interspecific competitive interactions while C is a matrix of intraspecific competitive interactions. Now, since C contains only intraspecific coefficients, C is a diagonal matrix, whereas B has zeros in the diagonal. The off-diagonal entries of B capture the symmetric pairwise-interspecific competitive interactions. Hence, we can write, A = B + C. Since, A is symmetric, B and C are also symmetric, and hence, all their eigenvalues are real. For an S species community, B and C will have S eigenvalues and these eigenvalues can be ordered as b1 ≥ b2 ≥ b3 ≥ … ≥ bs and c1 ≥ c2 ≥ c3 ≥ … ≥ cs. The Weyl’s inequality now states that the necessary and sufficient condition for multispecies coexistence to be stable is:

$$ {b}_1+{c}_s<0,\kern1em {b}_1+{c}_1<0, $$
(7)

When all intraspecific coefficients (or the diagonal elements of C) are equal, all the eigenvalues have the same value, i.e. cs = c1 = c, and the necessary and sufficient condition for stable species coexistence then becomes (Barabás et al. 2016):

$$ {b}_1+c<0. $$
(7a)

Weyl’s inequality, as defined, is a necessary and sufficient condition for stabilization of coexistence of the all the species in a community. However, it should be noted that Weyl’s inequality and stabilization of coexistence is only valid for pairwise competition matrices that are symmetric in nature. When A is not symmetric, decomposition of A into interspecific and intraspecific parts is not useful. Furthermore, whether multispecies coexistence is stabilized will be dependent on the eigenvalue distribution evaluated from the Jacobian matrix.

Effect of three-way HOIs on multispecies coexistence: when Weyl’s inequality was satisfied

We structured our simulations in the following way—we took a 50-species competitive community where pairwise intraspecific competition coefficients, i.e. diagonal elements of C, were kept the same at a value of 0.4. The pairwise-interspecific competitive interactions of B were drawn from a random uniform distribution in a way that it satisfied the Weyl’s inequality. Typically, Weyl’s inequality was fulfilled when intraspecific effects were much larger than all the interspecific effects (Barabás et al. 2016) and pairwise species coexistence was stabilized. Hence, for this particular case, we fixed intraspecific pairwise competition coefficient at 0.4, whereas interspecific pairwise competition was randomly sampled from a U (0.01, 0.04), where “U” stands for uniform random distribution. In all of our multispecies simulations, we ensured that species went extinct when densities dropped below 10−8 or were termed infeasible when it increased to 1012.

Effect of negative three-way HOIs on species coexistence: when Weyl’s inequality was satisfied

Next, we then evaluated how negative intra- and interspecific HOIs can destabilize species coexistence even when Weyl’s inequality was satisfied. For that, as a first-case example scenario, we ensured negative intraspecific HOIs to be stronger in magnitude (note the magnitude not the direction), drawn randomly from – U (0.08, 0.1), than negative interspecific HOIs, drawn randomly from – U (0.01, 0.04), and evaluated whether species coexistence was further stabilized. Note that here “−”in front of U (0.01, 0.04) or U (0.08, 0.1) means a negative sign. For a second case example scenario, we took negative intraspecific HOI to be lower in magnitude (drawn randomly from – U (0.01, 0.04)), than interspecific HOI (drawn randomly from – U (0.08, 0.1)) and evaluated whether it stabilized or destabilized species coexistence when Weyl’s inequality was satisfied. We used this range of HOI strengths following Bairey et al.’s (2016) study, where they showed that the critical strength of three-way HOIs that leads to community feasibility is of the order of −0.07 ± 0.04. However, we go beyond this threshold and explore across different HOI parameter values (see next paragraph).

Next, to evaluate whether the case examples are robust across a range of HOI parameter values, we did a series of replicate simulations for the case when Weyl’s inequality was satisfied. For each replicate simulation, negative interspecific HOIs were sampled from random – U (0.01, 0.09) and negative intraspecific HOIs were sampled from random – U (0.01 + d, 0.09 + d), where d increased from − 0.15 to 0.15 in steps of 0.05, such that when d = − 0.15, interspecific HOIs sampled are strictly of greater magnitude than intraspecific HOI, and when d = 0.15, magnitude of intraspecific HOIs is strictly greater than interspecific HOIs. For each d, we did 30 replicate simulations and noted the number of species that coexisted at the end of each replicate simulation.

Effect of positive three-way HOIs on species coexistence: when Weyl’s inequality was satisfied

We first show some example cases similar to “Effect of negative three-way HOIs on species coexistence: when Weyl’s inequality was satisfied”, where we first draw positive interspecific HOIs randomly from U (0.0003, 0.0005) and intraspecific HOIs from U (0.0007, 0.0009) and in the second case interspecific HOIs from U (0.0007, 0.0009) and intraspecific HOIs from U (0.0003, 0.0005) to evaluate whether strength of positive intra- vs. interspecific HOIs matters in stabilizing or destabilizing species coexistence when Weyl’s inequality was satisfied. Note the difference in the order of magnitudes used for parameterization of positive HOIs, which we discuss in the “Results” and in the “Discussion” sections.

Further, to evaluate whether the case examples are robust across a range of HOI parameter values, we did a series of replicate simulations for the case when Weyl’s inequality was satisfied. For each replicate simulation, positive interspecific HOIs were sampled from random U (0.0005, 0.0007) and positive intraspecific HOIs were sampled from random U (0.0005 + d, 0.0007 + d), where d increased from − 0.0002 to 0.0002 in steps of 0.0001, such that when d = − 0.0002, interspecific HOIs sampled are strictly of greater magnitude than intraspecific HOI, and when d = 0.0002 magnitude of intraspecific HOIs are strictly greater than interspecific HOIs. For each d, we did 30 replicate simulations and noted the number of species coexistence at the end of each simulation.

Effect of three-way HOIs on multispecies coexistence: when Weyl’s inequality was not satisfied

Typically, Weyl’s inequality was not fulfilled when intraspecific effects were smaller or similar in magnitude than all the interspecific effects (Barabás et al. 2016) and pairwise species coexistence was thus then destabilized. Hence, for this particular case, we fixed intraspecific pairwise competition coefficient at 0.01, whereas interspecific pairwise competition was randomly sampled from a U (0.01, 0.04). We thus wanted to evaluate whether HOIs (positive or negative) can promote species coexistence when pairwise species coexistence is destabilized.

Effect of negative three-way HOIs on species coexistence: Weyl’s inequality not satisfied

We first show some example cases similar to “Effect of negative three-way HOIs on species coexistence: when Weyl’s inequality was satisfied”, where we first draw intraspecific HOIs randomly from – U (0.08, 0.1) and interspecific HOIs from – U (0.01, 0.04) and in the second case intraspecific HOIs from – U (0.01, 0.04) and interspecific HOIs from – U (0.08, 0.1) to evaluate whether strength of intra vs. interspecific HOIs matters in stabilizing species coexistence when Weyl’s inequality was not satisfied.

For a wider parameter search, similar to “Effect of negative three-way HOIs on species coexistence: when Weyl’s inequality was satisfied”, negative interspecific HOIs were sampled from random – U (0.01, 0.09) and negative intraspecific HOIs were sampled from random – U (0.01 + d, 0.09 + d), where d increased from − 0.15 to 0.25 in steps of 0.05, such that when d = − 0.15, interspecific HOIs sampled are strictly of greater magnitude than intraspecific HOI, and when d = 0.25, magnitude of intraspecific HOIs is strictly greater than interspecific HOIs. For each d, we did 30 replicate simulations and noted the number of species coexistence at the end of each simulation.

Effect of positive three-way HOIs on species coexistence: Weyl’s inequality not satisfied

We first show some example cases similar to the above sections, where we first draw positive intraspecific HOIs randomly from U (0.0007, 0.0009) and interspecific HOIs from U (0.0003, 0.0005) and in the second case intraspecific HOIs from U (0.0003, 0.0005) and interspecific HOIs from (0.0007, 0.0009) to evaluate whether strength of intra- vs. interspecific HOIs matters in stabilizing species coexistence when Weyl’s inequality was not satisfied.

For a wider parameter search, for each replicate simulation, positive interspecific HOIs were sampled from random U (0.0005, 0.0007) and positive intraspecific HOIs were sampled from random U (0.0005 + d, 0.0007 + d), where d increased from − 0.0002 to 0.0002 in steps of 0.0001, such that when d = − 0.0002, interspecific HOIs sampled are strictly of greater magnitude than intraspecific HOI, and when d = 0.0002, magnitude of intraspecific HOIs is strictly greater than interspecific HOIs. For each d, we did 30 replicate simulations, and noted the number of species coexistence at the end of each simulation.

Results

Negative higher order interactions

In the pairwise interaction case, given a niche overlap of 0.5, coexistence between species 1 and 2 was only possible when fitness differences between them ranged from 0 to 2. Beyond a fitness difference of 2, coexistence was not possible, as the invasion growth rate of species 1 became negative (Fig. 1). Similarly, for species 2 also, invasion growth rate in the pairwise interaction case was negative once fitness difference increased beyond 2 (Fig. S1).

Fig. 1
figure 1

Effect of negative three-way HOIs on invasion growth rate (y-axis) of species 1 for pairwise species competition (dashed grey lines) and negative three-way HOIs (grey solid lines) for a range of fitness difference (x-axis). The red-dashed line marks the y-intercept at zero. Each panel of the plot compares invasion growth rate of species 1 under pairwise competition and under negative three-way HOIs. The panel label refers to the values of the HOI terms. For instance, the panel “symmetric” would mean that all the HOI terms in the HOI matrix have the exact same magnitude of − 0.01, and the panel β223 = − 0.01 (top row, species 1) would mean all the elements of HOI matrix are zero except β223 which is at − 0.01; β122 = − 0.01 (i.e. more negative and hence more increase in strength of intraspecific competition of species 2 α22). Panel β123 = − 0.01, would mean all terms of HOI matrix are zero except β123. Note that invasion growth rate of species 1 was negative across the range of fitness difference in negative three-way HOIs when interspecific competition was intensified more than intraspecific competition (panels: β122 =  − 0.01 and β123 =  − 0.01 for species 1)

Full size image

We found that, when three-way HOIs were negative, invasion growth rate of species 1 was positive across the range of fitness differences only when species 3 intensified intraspecific competition of species 2 (β223) more than it intensified interspecific competition (β123) (Fig. 1, β223 = − 0.01). However, if all negative three-way HOIs had the same magnitude, species coexistence was impossible even with low fitness differences (Fig. 1, symmetric case). Similarly, for species 2 (see Fig. S1) only when species 3 intensified pairwise-intraspecific competition of species 1, i.e., β113, invasion growth rate was positive across a wide range of fitness differences.

Negative four-way HOIs could also promote coexistence, even when fitness differences between two species were high (Fig. 3), provided their strength was an order of magnitude lower in comparison with interspecific and intraspecific pairwise and three-way higher order interactions (Appendix Figs. S3, S6, S7).

Positive higher order interactions

Generally, positive three-way HOIs led to species coexistence in fitness regions where species coexistence is impossible if only pairwise interactions are considered, even when there was substantial fitness difference between the two species (Fig. 2). For example, when species 3 alleviated intraspecific competition of species 2 (β223) more than it alleviated interspecific competition (of species 2 or 1), invasion growth rate of species 1 increased non-linearly, as fitness difference increased. This particular result could be understood by looking at the invasion growth rate of species 1 in the presence of non-zero β223 HOI (i.e. effect of species 3 on intraspecific interaction of species 2) while the rest of the HOI terms are zero, which becomes:

$$ {r}_1^{\ast }=1-\frac{\alpha_{12}}{\left({\alpha}_{22}-{\beta}_{223}{N}_3^{\ast}\right)} $$
(8)
Fig. 2
figure 2

Effect of positive three-way HOIs on invasion growth rate (y-axis) of species 1 for pairwise species competition (dashed grey lines) and positive three-way HOIs (grey solid lines in the figure) shown for a range of fitness difference (x-axis). The red-dashed line marks the y-axis at zero invasion growth rate. The panel label refers to the values of the terms in the HOI matrix. For instance, the panel “symmetric” would mean that all the HOI terms in the HOI matrix have the exact same magnitude of 0.01, and the panel β223 = 0.01 (top row, species 1) would mean all the elements of HOI matrix are zero except for β223 which is at 0.01 (i.e. more positive and hence more decrease in strength of intraspecific competition of species 2, α22). Panel β213 = − 0.01 would mean all terms of HOI matrix are zero except β213. Note that invasion growth rate of species 1 was positive across the range of fitness difference in all positive three-way HOIs

Full size image

We specifically wanted to evaluate whether HOIs could stabilize (or destabilize) pairwise species coexistence in regions where pairwise species coexistence is impossible, for instance in regions of extreme fitness difference with considerable niche overlap. We thus varied α12, α22 in a way given by Eq. 6 such that as fitness difference increased (due to θ varying from 0 to 7 in Eq. 6) interspecific effect of species 2 on species 1, i.e. α12 increased more rapidly than intraspecific competition α22 of species 2. For instance, as θ increases in Eq. 6, α12 increases at the same rate as θ, while α22 in increases by 2\( \sqrt{\theta } \). So, at high θ values, fitness difference between the two species becomes extreme and pairwise coexistence becomes impossible. But, with \( {N}_3^{\ast }=\frac{1}{\alpha_{33}}=10 \) and with the presence of only one HOI term β223 of 0.01 causes the invasion growth rate of species 1 (Eq. 8) to increase rapidly as fitness difference increases. In all the cases, positive HOIs lead to very high invasion growth rates (Fig. 3). We observed similar results for species 2 too (see Fig. S2).

Fig. 3
figure 3

Invasion growth rate of species 1 (y-axis), for pairwise species competition (dashed grey lines) and negative four-way HOIs (grey lines) shown for a range of fitness differences (x-axis). The red-dashed line marks the y-intercept at zero. Each panel of the plot compares invasion growth rate of species 1 in pairwise competition (dashed grey lines) and in negative four-way HOIs (solid grey continuous lines). The leftmost top panel in the first row of the figure (“symmetric”) denotes the case where all HOI terms are negative and have the same strength of − 0.001. In this panel, invasion growth rate of species 1 in pure pairwise interaction with species 2 (dashed line) is plotted against invasion growth rate of species 1 in four-way HOIs (grey continuous line) where all the terms of HOIs have same strength. In the other panels, invasion growth rate of species 1 in pure pairwise competition (dashed line) is plotted against invasion growth rate of species 1 in HOIs, when only one term from HOIs is perturbed, denoted by the panel labelled as γ1222, γ1223, γ1232, γ1233 with the rest of the terms of the four-way HOI matrix being 0. For example, γ1223 = − 0.001 would mean all the HOI terms are zero except γ1223 which is at − 0.001. Here N2 = N3 = 5

Full size image

Generally, positive four-way HOIs led to species coexistence despite extreme fitness differences. As strength of positive four-way HOIs increased, invasion growth rate of species 1 when species 2 and species 3 are present also increased (Appendix Figs. S4, S5, S8, S9).

Effect of negative HOIs on multispecies coexistence when Weyl’s inequality was satisfied

When Weyl’s inequality was satisfied, a 50 species community was feasible and stable in the absence of HOIs (Fig. 4c). Interestingly, in the presence of negative HOIs, even when Weyl’s inequality was satisfied, there were parameter regions where species coexistence could be destabilized. This was particularly evident when the magnitude intraspecific HOIs were strictly less strong in comparison with interspecific HOIs (Fig. 5a).

Fig. 4
figure 4

Multispecies coexistence in the presence and absence of negative HOIs that either satisfy Weyl’s inequality (b1 + cs < 0) or do not (b1 + cs > 0). When 50 species in a competitive community compete in pairwise manner and satisfy Weyl’s inequality (c), species coexistence was always stabilized. However, in the presence of negative HOIs where magnitude of interspecific HOIs was strictly greater than intraspecific HOIs, i.e. βiik < βijk (a), species coexistence was destabilized (d). When the opposite happens: (b), i.e. βiik > βijk which suggests magnitude of intraspecific HOIs is greater than interspecific HOIs, species coexistence was again stabilized. When Weyl’s inequality was not satisfied (b1 + cs > 0) (f), pairwise coexistence was impossible. In the presence of negative HOIs, species coexistence was not possible if magnitude of interspecific HOIs was greater than intraspecific HOIs (g), while species coexistence was stabilized if magnitude of intraspecific HOIs was greater than interspecific HOIs (h)

Full size image
Fig. 5
figure 5

Each black dot represents replicate simulations (n = 30) of negative HOI model for each d values that capture the difference in the range of sampling of negative three-way interspecific HOIs and negative three-way intraspecific HOIs. Pairwise competition matrix, i.e. α matrix is symmetric. Interspecific negative HOIs are randomly sampled for each replicate simulation from U (0.01, 0.09) but intraspecific negative HOIs were sampled from U (0.01 + d, 0.09 + d). Thus, negative d would mean the overall magnitude of negative intraspecific HOIs are of lower strength than the magnitude of negative interspecific HOIs, and positive d would mean vice versa. a Weyl’s inequality is satisfied: pairwise species coexistence destabilized by the parameterization of the HOIs when d is negative, but further stabilized by HOIs when d is positive. This indicates that the magnitude of negative intraspecific HOIs sampled should be higher than the magnitude of negative interspecific HOIs. b Weyl’s inequality not satisfied: pairwise coexistence will always be destabilized. But HOIs can still promote and stabilize coexistence provided d is positive and large, indicating that intraspecific HOIs should be strictly greater in magnitude then interspecific negative HOIs

Full size image

Effect of positive HOIs on species coexistence when Weyl’s inequality was satisfied

From our analytical results, we found that positive HOIs could lead to extremely high invasion growth rates as high as 250 (Fig. 2) and are extremely sensitive to changes in parameter values (Fig. S10). With such high invasion growth rates for species in a high-dimensional system such as a 50-species community, species coexistence becomes unfeasible with higher strength in positive HOIs. Hence, in our multispecies simulations, we used the parameter positive HOI values that lead to feasible densities (i.e. below 1012). Our results also demonstrated that when Weyl’s inequality is satisfied, species coexistence is possible irrespective of differences in intra- and interspecific positive HOIs (Figs. S14 and S16).

Effect of negative HOIs on species coexistence when Weyl’s inequality was not satisfied

Failing to fulfil the Weyl’s inequality led to disruption of pairwise species coexistence (Fig. 4g). But even when Weyl’s inequality was not fulfilled, and hence pairwise species coexistence was impossible, negative HOIs could still promote and stabilize species coexistence provided intraspecific HOIs were strictly greater in magnitude than interspecific HOIs (Fig. 5b).

Effect of positive HOIs on species coexistence when Weyl’s inequality was not satisfied

When Weyl’s inequality was not satisfied, irrespective of differences in magnitude of inter- and intraspecific positive HOIs, stable coexistence of all species was not possible (see Figs. S14, S15, and S16).

Discussion

While the assumption that interactions between species is pairwise is pervasive in coexistence theory (Valladares et al. 2015; Levine et al. 2017; Gallien et al. 2017), the role of HOIs in stabilizing or destabilizing coexistence has been relatively understudied both empirically and theoretically (Abrams 1983; Wilson 1992; Grilli et al. 2017; Letten and Stouffer 2019). Recently, empirical understanding of species coexistence has been sought through MCT where invading potential of a species in the presence of an established competitor species is explored (Grainger et al. 2019). In this context, the role of HOIs needs to be elucidated to understand multispecies coexistence, and its underlying mechanisms (Levine et al. 2017). Our study fills this gap by using concepts from MCT (Chesson 2000), to evaluate the effect of HOIs on species coexistence in a simple three-species Lotka-Volterra model and in a more complex multispecies community. This study shows that HOIs can, in general, promote species coexistence in parameter spaces where pairwise coexistence is unstable, provided certain conditions are fulfilled.

Using a three-species Lotka-Volterra model, our results showed that positive three-way HOI’s stabilize coexistence across a wide range of fitness differences irrespective of differences in strength of inter- and intraspecific competition, while negative three-way HOI’s can stabilize coexistence only if intraspecific competition was strengthened more than interspecific competition. If, however, in the case of three-way negative HOIs, interspecific competition was strengthened more in relation to intraspecific competition, species coexistence was impossible even when fitness differences were negligible.

Three-way positive HOIs in our model leads to a decrease in the per capita strength of pairwise competition between two species, while three-way negative HOIs lead to an increase in the per capita strength of competition. Moreover, three-way positive HOIs could lead to disproportionately high invasion growth rates, particularly when they alleviated intraspecific competition (Fig. 2). Such high invasion growth rate could affect stability and lead to infeasible species densities (Terry et al. 2017, 2018). In a multispecies community, having a higher strength in positive three-way HOIs could lead to disproportionately high growth rates for all species, which could lead to community infeasibility. However, when fitness differences between two species were extremely high, positive HOIs could lead to species coexistence by decreasing interspecific competition more than intraspecific competition. For example, for extreme fitness difference between species 1 and species 2, say fitness difference of 3, positive HOIs (β123 = 0.01), which decreases interspecific competition of species 2 on species 1, will lead to species coexistence and proportionately low but positive invasion growth rate (Fig. 2). An earlier study had reported that HOIs could positively influence fitness of species by alleviating the dominating effect of neighbouring species (Bairey et al. 2016; Mayfield and Stouffer 2017).

Invasion growth rates were generally sensitive to changes in the strength of HOIs, for both positive and negative HOIs (Appendix Fig. S11), which suggests that slight parameter changes in HOIs have the potential to affect species coexistence. Hence, under restricted parameter space, HOIs could stabilize species coexistence (AlAdwani and Saavedra 2019). Four-way HOIs could lead to species coexistence across a range of fitness differences, if their direction was positive (Appendix Figs. S5–6; S8–9). Negative four-way HOIs could also promote coexistence, provided their strength was an order of magnitude lower in comparison with pairwise interactions (Appendix Figs. S3–4, S7–8). It is possible that four-way interactions could be prevalent in species communities, but empirically parameterizing such four-way interactions would be a difficult task (Mayfield and Stouffer 2017). In addition, incorporating four-way HOIs in models parameterized from empirical data might not always provide additional explanatory power (AlAdwani and Saavedra 2019).

In a multispecies community, pairwise coexistence rule does not hold (Barabás et al. 2016). However, multispecies coexistence could be understood by analyzing the Weyl inequality (Fulton 2000). When Weyl’s inequality was not satisfied, pairwise species coexistence was impossible. Generally, for Weyl’s inequality to be satisfied, intraspecific competition of all the species in a community should be significantly larger than their interspecific competitive coefficients. In our models, positive HOIs decrease pairwise competition in a multispecies community, which consequently led to species coexistence provided self-regulation was strong (Mayfield and Stouffer 2017) and Weyl’s inequality was satisfied (Fig. 5; appendix Fig. S15). When Weyl’s equality was not satisfied, species coexistence becomes impossible both with pairwise and with positive HOIs. Pairwise coexistence in such a scenario was disrupted particularly because strength of intraspecific competition was on average lower than interspecific competition between species. However, surprisingly, even in such a scenario, negative HOIs could stabilize species coexistence (Fig. 4h). Weyl’s inequality is valid only for symmetric matrices, and in our multispecies simulations, we specifically focused on symmetric pairwise competition matrices. However, even when we used asymmetric pairwise competition matrices, the effect of three-way HOIs on multispecies coexistence were similar to the results obtained from symmetric pairwise competition matrices (Fig. S15).

If intraspecific HOIs are strictly greater than interspecific HOIs, species coexistence in a large competitive community was stabilized, irrespective of whether the pairwise competition matrices were symmetric or asymmetric. This is analogous to the two-species coexistence rule that species must limit themselves more than they limit competitors. In general, the simplest way to generalize multispecies coexistence in the presence of three-way negative HOIs was that—when pairwise coexistence for multispecies community was impossible (Weyl’s inequality not satisfied), intraspecific competition should be strengthened more than interspecific competition by three-way HOIs. If even a single species in the multispecies community violated this rule, coexistence in the multispecies community was disrupted, even in the presence of HOIs (Fig. S12). On the contrary, when pairwise multispecies coexistence was possible (Weyl’s inequality was satisfied and pairwise-intraspecific effects were substantially larger than interspecific effects), having intraspecific HOIs that have similar strength as interspecific HOIs could still stabilize species coexistence (Appendix Fig. S13; Figs. 4 and 5).

The Lotka-Volterra models of competition have been used extensively to understand mechanisms that could promote species coexistence through pairwise interactions (Barabás et al. 2016) and through HOIs (Abrams 1983; Wilson 1992; Bairey et al. 2016; Mayfield and Stouffer 2017; Grilli et al. 2017; Letten and Stouffer 2019). Using the Lotka-Volterra models, our results add to the emerging body of knowledge of HOIs and the prospect of extending simple tractable dynamics to understanding complex multispecies dynamics. It is straightforward, although challenging, to estimate pairwise competition coefficients from experiments that involve manipulating competitor densities, or from long-term observational field data. However, there are obstacles in estimating higher order coefficients and to singularly understand their effects on species coexistence (Mayfield and Stouffer 2017). To implement a higher order model empirically with just three competitors would require no less than 27 parameters at the very least. To collect such amount of empirical data, whether observational or experimental is an enormous task. Nonetheless, invasion growth rates could be estimated from empirical data that are used to parameterize community models or mechanistic competition models by explicitly incorporating HOIs (Tilman 1994; Letten and Stouffer 2019).

Our three-species analytical model might seem biologically unrealistic, but was specifically built to make the analysis simple and tractable. For instance, with all the three species interacting both in pairwise and three-way HOIs, the number of HOI terms increases to 27, and including four-way HOIs, the number of terms goes up to 54. With this number of HOI terms, assessing and analyzing the effect of HOIs on species coexistence could be manifold. Although we have much to understand about the effects of HOIs empirically, it is clear that effects of HOIs on species coexistence were dependent on the strength as well as the direction of the HOIs. That being said, HOI terms in ecological models could increase the number of equilibrium points exponentially, though such equilibrium points can be ecologically feasible only under a restricted set of parameter space (AlAdwani and Saavedra 2019).