Intraspecific predator interference promotes biodiversity in ecosystems

eLife assessment

This manuscript is an important contribution, assessing the role of intraspecific consumer interference in maintaining diversity using a mathematical model. Consistent with long-standing ecological theory, the authors convincingly show that predator interference allows for the coexistence of multiple species on a single resource, beyond the competitive exclusion principle. Notably, the model matches observed rank-abundance curves in several natural ecosystems.

https://doi.org/10.7554/eLife.93115.3.sa0

Significance of the findings:

Important: Findings that have theoretical or practical implications beyond a single subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Convincing: Appropriate and validated methodology in line with current state-of-the-art

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
eLife digest
Introduction
Results
Discussion
Materials and methods
Appendix 1
Appendix 2
Appendix 3
Appendix 4
Appendix 5
Appendix 6
Appendix 7
Appendix 8
Appendix 9
Data availability
References
Article and author information
Metrics

Abstract

Explaining biodiversity is a fundamental issue in ecology. A long-standing puzzle lies in the paradox of the plankton: many species of plankton feeding on a limited variety of resources coexist, apparently flouting the competitive exclusion principle (CEP), which holds that the number of predator (consumer) species cannot exceed that of the resources at a steady state. Here, we present a mechanistic model and demonstrate that intraspecific interference among the consumers enables a plethora of consumer species to coexist at constant population densities with only one or a handful of resource species. This facilitated biodiversity is resistant to stochasticity, either with the stochastic simulation algorithm or individual-based modeling. Our model naturally explains the classical experiments that invalidate the CEP, quantitatively illustrates the universal S-shaped pattern of the rank-abundance curves across a wide range of ecological communities, and can be broadly used to resolve the mystery of biodiversity in many natural ecosystems.

eLife digest

The surface waters of the ocean are teeming with microscopic creatures known as plankton, which get carried across vast distances by the currents. In a single ecosystem, thousands of plankton species may coexist, all competing for very few types of food sources. According to the principle of competitive exclusion, this should not be the case. Indeed, this theory states that the population levels of two species competing for the same resource cannot remain steady over time – or more generally, that the number of consumer species in an ecosystem cannot be higher than the number of resource types on which they rely. And yet, the Earth abounds with examples where a limited variety of resources supports a large number of competing yet coexisting consumer species. This is known as the paradox of the plankton.

Many models have been proposed to explain how the limitations set by the competitive exclusion principle can be overcome, yet it is still unknown how to resolve the paradox of the plankton in a steady environment. In response, Kang et al. set out to test whether a phenomenon known as predator interference, which emerges when two or more individuals of the same species compete for the same resources, could help address the paradox of the plankton.

To test this idea, Kang et al. developed a mathematical model of predator interference for multiple species of plankton feeding on a limited variety of food sources. The model put predators of the same species into encountering pairs to simulate predator interference. In this scenario, a wide range of predator species were able to live alongside each other with the numbers of each type of predator remaining steady over time.

These results can be understood as follows: as a species becomes more successful at extracting resources from its environment, its population grows – and with it, the number of individuals engaged in intraspecific interference. Locked in interference, these species become less effective at getting food, creating a negative feedback loop that slows down the expansion of the species, allowing others to occupy the same niche.

These findings may benefit ecologists and conservationists by offering insights into how to maintain biodiversity in ecosystems and protect endangered species. Further work is needed to test how well the model applies to other ecosystems.

Introduction

The most prominent feature of life on Earth is its remarkable species diversity: countless macro- and micro-species fill every corner on land and in the water (Pennisi, 2005; Hoorn et al., 2010; de Vargas et al., 2015; Daniel, 2005). In tropical forests, thousands of plant and vertebrate species coexist (Hoorn et al., 2010). Within a gram of soil, the number of microbial species is estimated to be 2000–18,000 (Daniel, 2005). In the photic zone of the world ocean, there are roughly 150,000 eukaryotic plankton species (de Vargas et al., 2015). Explaining this astonishing biodiversity is a major focus in ecology (Pennisi, 2005). A great challenge stems from the well-known competitive exclusion principle (CEP): two species competing for a single type of resources cannot coexist at constant population densities (Gause, 1934; Hardin, 1960), or generically, in the framework of consumer-resource models, the number of consumer species cannot exceed that of resources at a steady state (MacArthur and Levins, 1964; Levin, 1970; McGehee, 1977). On the contrary, in the paradox of plankton, a limited variety of resources supports hundreds or more coexisting species of phytoplankton (Hutchinson, 1961). Then, how can plankton and many other organisms somehow liberate the constraint of CEP?

Ever since MacArthur and Levin proposed the classical mathematical proof for CEP in the 1960s (MacArthur and Levins, 1964), various mechanisms have been put forward to overcome the limits set by CEP (Chesson, 2000). Some suggest that the system never approaches a steady state where the CEP applies, due to temporal variations (Hutchinson, 1961; Levins, 1979), spatial heterogeneity (Levin, 1974), or species’ self-organized dynamics (Koch, 1974; Huisman and Weissing, 1999). Others consider factors such as toxins (Czárán et al., 2002), cross-feeding (Goyal and Maslov, 2018; Goldford et al., 2018; Niehaus et al., 2019), spatial circulation (Villa Martín et al., 2020; Gupta et al., 2021), ‘kill the winner’ (Thingstad, 2000), pack hunting (Wang and Liu, 2020), collective behavior (Dalziel et al., 2021), metabolic trade-offs (Posfai et al., 2017; Weiner et al., 2019), co-evolution (Xue and Goldenfeld, 2017), and other complex interactions among the species (Beddington, 1975; DeAngelis et al., 1975; Arditi and Ginzburg, 1989; Kelsic et al., 2015; Grilli et al., 2017; Ratzke et al., 2020). However, questions remain as to what determines species diversity in nature, especially for quasi-well-mixed systems such as that of plankton (Pennisi, 2005; Sunagawa et al., 2020).

Among the proposed mechanisms, predator interference, specifically the pairwise encounters among consumer individuals, emerges as a potential solution to this issue. Predator interference is commonly described by the classical Beddington-DeAngelis (B-D) phenomenological model (Beddington, 1975; DeAngelis et al., 1975). Through the application of the B-D model, several studies (Cantrell et al., 2004; Hsu et al., 2013) have shown that intraspecific predator interference can break CEP and facilitate species coexistence. However, from a mechanistic perspective, the functional response of the B-D model can be formally derived from a scenario solely involving chasing pairs, representing the consumption process between consumers and resources, without accounting for pairwise encounters among consumer individuals (Wang and Liu, 2020; Huisman and De Boer, 1997). Disturbingly, it has been established that the scenario involving only chasing pairs is subject to the constraint of CEP (Wang and Liu, 2020), raising doubt regarding the validity of applying the B-D model to overcome the CEP.

In this work, building upon MacArthur’s consumer-resource model framework (Arthur, 1969; MacArthur, 1970; Chesson, 1990), and drawing on concepts from chemical reaction kinetics (Ruxton et al., 1992; Huisman and De Boer, 1997; Wang and Liu, 2020), we present a mechanistic model of predator interference that extends the B-D phenomenological model (Beddington, 1975; DeAngelis et al., 1975) by providing a detailed consideration of pairwise encounters. The intraspecific interference among consumer individuals effectively constitutes a negative feedback loop, enabling a wide range of consumer species to coexist with only one or a few types of resources. The coexistence state is resistant to stochasticity and can hence be realized in practice. Our model is broadly applicable and can be used to explain biodiversity in many ecosystems. In particular, it naturally explains species coexistence in classical experiments that invalidate CEP (Ayala, 1969; Park, 1954) and quantitatively illustrates the S-shaped pattern of the rank-abundance curves in an extensive spectrum of ecological communities, ranging from the communities of ocean plankton worldwide (Fuhrman et al., 2008; Ser-Giacomi et al., 2018), tropical river fishes from Argentina (Cody and Smallwood, 1996), forest bats of Trinidad (Clarke et al., 2005), rainforest trees (Hubbell, 2001), birds (Terborgh et al., 1990; Martínez et al., 2023), butterflies (De Vries, 1997) in Amazonia, to those of desert bees (Hubbell, 2001) in Utah and lizards from South Australia (Cody and Smallwood, 1996).

Results

A generic model of pairwise encounters

Here we present a mechanistic model of pairwise encounters (see Figure 1A), where $S_{C}$ consumer species ${C_{1}, \dots, C_{S_{C}}}$ compete for $S_{R}$ resource species ${R_{1}, \dots, R_{S_{R}}}$ . The consumers are biotic, while the resources can be either biotic or abiotic. For simplicity, we assume that all species are motile and move at certain speeds, namely, $v_{C_{i}}$ for consumer species $C_{i}$ and $v_{R_{l}}$ for resource species $R_{l}$ . For abiotic resources, they cannot propel themselves but may passively drift due to environmental factors. Each consumer is free to feed on one or multiple types of resources, while consumers do not directly interact with one another except through pairwise encounters.

Figure 1

Download asset Open asset

A model of pairwise encounters may naturally break CEP.

(A) A generic model of pairwise encounters involving $S_{C}$ consumer species and $S_{R}$ resource species. (B) The well-mixed model of (A). (**C–E**) Time courses of two consumer species competing for one resource species. (**F–H**) Positive solutions to the steady-state equations (see Equations S38 and S65): $\dot{R} = 0$ (orange surface), ${\dot{C}}_{1} = 0$ (blue surface), ${\dot{C}}_{2} = 0$ (green surface), that is the zero-growth isoclines. The black/red dot represents the unstable/stable fixed point, while the dotted lines in (E) are the analytical solutions of the steady-state abundances (marked with superscript ‘(A)’). See Appendix 1—tables 1 and 2 for the definitions of symbols. See Appendix 9 for simulation details.

Then, we explicitly consider the population structure of consumers and resources: some wander around freely, undergoing Brownian motions; others encounter one another, forming ephemeral entangled pairs. Specifically, when a consumer individual $C_{i}$ and a resource $R_{l}$ come close within a distance of $r_{i l}^{(C)}$ (see Figure 1A), the consumer can chase the resource and form a chasing pair: $C_{i}^{(P)} ⋁ R_{l}^{(P)}$ (see Figure 1B), where the superscript ‘(P)’ represents ‘pair’. The resource can either escape at rate $d_{i l}$ or be caught and consumed by the consumer with rate $k_{i l}$ . Meanwhile, when a $C_{i}$ individual encounters another consumer $C_{j}$ within a distance of $r_{i j}^{(I)}$ (see Figure 1A), they can stare at, fight against, or play with each other, thus forming an interference pair: $C_{i}^{(P)} ⋁ C_{j}^{(P)}$ (see Figure 1B). This paired state is evanescent, with consumers separating at rate $d_{i j}^{'}$ . For simplicity, we assume that all $r_{i l}^{(C)}$ and $r_{i j}^{(I)}$ are identical, respectively, that is $\forall i, j, l$ , $r_{i l}^{(C)} = r^{(C)}$ and $r_{i j}^{(I)} = r^{(I)}$ .

In a well-mixed system of size $L^{2}$ , the encounter rates among individuals, $a_{i l}$ and $a_{i j}^{'}$ (see Figure 1B), can be derived using the mean-field approximation: $a_{i l} = 2 r^{(C)} L^{- 2} \sqrt{v_{C_{i}}^{2} + v_{R_{l}}^{2}}$ and $a_{i j}^{'} = 2 r^{(I)} L^{- 2} \sqrt{v_{C_{i}}^{2} + v_{C_{j}}^{2}}$ (see Materials and methods, and Appendix 1—figure 1). Then, we proceed to analyze scenarios involving different types of pairwise encounters (see Figure 1B). For the scenario involving only chasing pairs, the population dynamics can be described as follows:

\begin{array}{ll} {\begin{cases} {\dot{x}}_{i l} = a_{i l} C_{i}^{(F)} R_{l}^{(F)} - (d_{i l} + k_{i l}) x_{i l}, \\ {\dot{C}}_{i} = \sum_{l = 1}^{S_{R}} w_{i l} k_{i l} x_{i l} - D_{i} C_{i}, i = 1, \dots, S_{C}, \\ {\dot{R}}_{l} = g_{l} ({R_{l}}, {x_{i}}, {C_{i}}), l = 1, \dots, S_{R}, \end{cases} \end{array}

where $x_{i l} \equiv C_{i}^{(P)} ⋁ R_{l}^{(P)}$ , g_l is an unspecified function, the superscript ‘(F)’ represents the freely wandering population, $D_{i}$ denotes the mortality rate of $C_{i}$ , and $w_{i l}$ is the mass conversion ratio (Wang and Liu, 2020) from resource $R_{l}$ to consumer $C_{i}$ . With the integration of intraspecific predator interference, we combine Equation 1 and the following equation:

{\dot{y}}_{i} = a_{i}^{'} [C_{i}^{(F)}]^{2} - d_{i}^{'} y_{i},

where $a_{i}^{'} = a_{i i}^{'}$ , $d_{i}^{'} = d_{i i}^{'}$ , and $y_{i} \equiv C_{i}^{(P)} ⋁ C_{i}^{(P)}$ represents the intraspecific interference pair (see Figure 1B). For the scenario involving chasing pairs and interspecific interference, we combine Equation 1 with the following equation:

{\dot{z}}_{i j} = a_{i j}^{'} C_{i}^{(F)} C_{j}^{(F)} - d_{i j}^{'} z_{i j}, i \neq j,

where $z_{i j} \equiv C_{i}^{(P)} ⋁ C_{j}^{(P)}$ stands for the interspecific interference pair (see Figure 1B). In the scenario where chasing pairs and both intra- and interspecific interference are relevant, we combine Equations 1–3, and the populations of consumers and resources are given by $C_{i} = C_{i}^{(F)} + \sum_{l} x_{i l} + 2 y_{i} + \sum_{j \neq i} z_{i j}$ and $R_{l} = R_{l}^{(F)} + \sum_{i} x_{i l}$ , respectively.

Generically, the consumption and interference processes are much quicker compared to the birth and death processes. Thus, in the derivation of the functional response, $F (R_{l}, C_{i}) \equiv k_{i l} x_{i l} / C_{i}$ , the consumption and interference processes are supposed to be in fast equilibrium. In all scenarios involving different types of pairwise encounters, the functional response in the B-D model is a good approximation only for a special case with $d_{i l} \approx 0$ and $R_{l} ≫ \sum_{i = 1}^{S_{C}} C_{i}$ (see Appendix 1—figure 2 and Appendix 3 for details).

To facilitate further analysis, we assume that the population dynamics of the resources follows the same construction rule as that in MacArthur’s consumer-resource model (Arthur, 1969; MacArthur, 1970; Chesson, 1990). Then,

\begin{array}{ll} g_{l} ({R_{l}}, {x_{i}}, {C_{i}}) = {\begin{cases} η_{l} R_{l} (1 - R_{l} / κ_{l}) - \sum_{i = 1}^{S_{C}} k_{i l} x_{i l} (f o r b i o t i c r e s o u r c e s); \\ ζ_{l} (1 - R_{l} / κ_{l}) - \sum_{i = 1}^{S_{C}} k_{i l} x_{i l} (f o r a b i o t i c r e s o u r c e s) . \end{cases} \end{array}

In the absence of consumers, biotic resources exhibit logistic growth. Here, $η_{l}$ and $κ_{l}$ represent the intrinsic growth rate and the carrying capacity of species $R_{l}$ . For abiotic resources, $ζ_{l}$ stands for the external resource supply rate of $R_{l}$ , and $κ_{l}$ is the abundance of $R_{l}$ at a steady state without consumers. For simplicity, we focus our analysis on abiotic resources, although all results generally apply to biotic resources as well. By applying dimensional analysis, we render all parameters dimensionless (see Appendix 7). For convenience, we retain the same notations below, with all parameters considered dimensionless unless otherwise specified.

Intraspecific predator interference facilitates species coexistence and breaks CEP

To clarify the specific mechanisms that can facilitate species coexistence, we systematically investigate scenarios involving different forms of pairwise encounters in a simple case with $S_{C} = 2$ and $S_{R} = 1$ . To simplify the notations, we omit the subscript/superscript ‘ $l$ ’ since $S_{R} = 1$ . For clarity, we assign each consumer species of unique competitiveness by setting that the mortality rate $D_{i}$ is the only parameter that varies with the consumer species.

First, we conduct the analysis within a deterministic framework using ordinary differential equations (ODEs). In the scenario involving only chasing pairs, consumer species cannot coexist at a steady state except for special parameter settings (sets of measure zero) (Wang and Liu, 2020). In practice, if all species coexist, the steady-state equations of the consumer species ( ${\dot{C}}_{i} = 0$ , i.e. the zero-growth isolines) yield $f_{i} (R^{(F)}) = D_{i}$ ( $i = 1, 2$ ), where $f_{i} (R^{(F)})$ is defined as $f_{i} (R^{(F)}) \equiv R^{(F)} / (R^{(F)} + K_{i})$ and $K_{i} \equiv (d_{i} + k_{i}) / a_{i}$ . These equations form two parallel surfaces in the $(C_{1}, C_{2}, R)$ coordinates, making steady coexistence impossible (Wang and Liu, 2020; see Figure 1C and F and Appendix 1—figure 3A–C).

Meanwhile, in the scenario involving chasing pairs and interspecific interference, if all species coexist, the zero-growth isolines of the three species (see Equation S65) correspond to three non-parallel surfaces $Ω_{i}^{'} (R, C_{1}, C_{2}) = D_{i}$ ( $i = 1, 2$ ), $G^{'} (R, C_{1}, C_{2}) = 0$ (see Figure 1G and Appendix 1—figure 3D; refer to Appendix 5 for definitions of $Ω_{i}^{'}$ and $G^{'}$ ), which can intersect at a common point (fixed point). However, this fixed point is unstable (see Figure 1G and Appendix 1—figure 3D, F), and thus one of the consumer species is doomed to extinction (see Figure 1D).

Next, we turn to the scenario involving chasing pairs and intraspecific interference. Likewise, steady coexistence requires (see Equation S38) that three non-parallel surfaces $Ω_{i} (R, C_{1}, C_{2}) = D_{i}$ ( $i = 1, 2$ ), $G (R, C_{1}, C_{2}) = 0$ cross at a common point (see Figure 1H and Appendix 1—figure 3G; refer to Appendix 4 for definitions of $Ω_{i}$ and $G$ ). Indeed, this naturally happens, and encouragingly the fixed point can be stable. Therefore, two consumer species may stably coexist at a steady state with only one type of resources, which obviously breaks CEP (see Figure 1E and Appendix 1—figure 4A). In fact, the coexisting state is globally attractive (see Appendix 1—figure 4A), and there exists a non-zero volume of parameter space where the two consumer species stably coexist at constant population densities (see Appendix 1—figure 4B, C), demonstrating that the violation of CEP does not depend on special parameter settings. We further consider the scenario involving chasing pairs and both intra- and interspecific interference (see Appendix 1—figure 5). Much as expected, the species coexistence behavior is very similar to that without interspecific interference.

Intraspecific interference promotes biodiversity in the presence of stochasticity

Stochasticity is ubiquitous in nature. However, it is prone to jeopardize species coexistence (Xue and Goldenfeld, 2017). Influential mechanisms such as ‘kill the winner’ fail when stochasticity is incorporated (Xue and Goldenfeld, 2017). Consistent with this, we observe that two notable cases of oscillating coexistence (Koch, 1974; Huisman and Weissing, 1999) turn into species extinction when stochasticity is introduced (see Appendix 1—figure 6A, B), where we simulate the models with stochastic simulation algorithm (SSA; Gillespie, 2007) and adopt the same parameters as those in the original references (Koch, 1974; Huisman and Weissing, 1999).

Then, we proceed to investigate the impact of stochasticity on our model using SSA (Gillespie, 2007). In the scenario involving chasing pairs and intraspecific interference, species may coexist indefinitely in the SSA simulations (see Figure 2A and Appendix 1—figure 4D). In fact, the parameter region for species coexistence in this scenario is rather similar between the SSA and ODEs studies (see Appendix 1—figure 6C, D). Similarly, in the scenario involving chasing pairs and both inter- and intraspecific interference, all species may coexist indefinitely in company with stochasticity (see Appendix 1—figure 5D).

Figure 2

Download asset Open asset

Intraspecific predator interference facilitates species coexistence regardless of stochasticity.

(**A, B**) Time courses of the species abundances simulated with ODEs, SSA, or IBM. (C) Snapshots of the IBM simulations. (**D, E**) A model of intraspecific predator interference explains two classical laboratory experiments that invalidate CEP. (D) In Ayala’s experiment, two *Drosophila* species coexist with one type of resources within a laboratory bottle (Ayala, 1969). (E) In Park’s experiment, two *Tribolium* species coexist for 2 years with one type of food (flour) within a lab (Park, 1954). See Appendix 1—figure 7C, D for the comparison between model results and experimental data using Shannon entropies.

To further mimic a real ecosystem, we resort to individual-based modeling (IBM; Grimm, 2013; Vetsigian, 2017), an essentially stochastic simulation method. In the simple case of $S_{C} = 2$ and $S_{R} = 1$ , we simulate the time evolution of a 2-D square system in a size of $L^{2}$ with periodic boundary conditions (see Materials and methods for details). In the scenario involving chasing pairs and intraspecific interference, two consumer species coexist for long with only one type of resources in the IBM simulations (see Figure 2B and C). Together with the SSA simulation studies, it is obvious that intraspecific interference still robustly promotes species coexistence when stochasticity is considered.

Comparison with experimental studies that reject CEP

In practice, two classical studies (Ayala, 1969; Park, 1954) reported that, in their respective laboratory systems, two species of insects coexisted for roughly years or more with only one type of resources. Evidently, these two experiments (Ayala, 1969; Park, 1954) are incompatible with CEP, while factors such as temporal variations, spatial heterogeneity, cross-feeding, etc. are clearly not involved in such systems. As intraspecific fighting is prevalent among insects (Boomsma et al., 2005; Dankert et al., 2009; Chen et al., 2002), we apply the model involving chasing pairs and intraspecific interference to simulate the two systems. Overall, our SSA results show good consistency with those of the experiments (see Figure 2D and E and Appendix 1—figure 7). The fluctuations in experimental time series can be mainly accounted by stochasticity.

A handful of resource species can support a wide range of consumer species regardless of stochasticity

To resolve the puzzle stated in the paradox of the plankton, we analyze the generic case where $S_{C}$ consumer species compete for $S_{R}$ resource species (with $S_{C} > S_{R}$ ) within the scenario involving chasing pairs and intraspecific interference. The population dynamics is described by equations combining Equations 1, 2 and 4. As with the cases above, each consumer species is assigned a unique competitiveness through a distinctive $D_{i}$ .

Strikingly, a plethora of consumer species may coexist at a steady state with only one resource species ( $S_{C} ≫ S_{R}$ , $S_{R} = 1$ ) in the ODEs simulations, and crucially, the facilitated biodiversity can still be maintained in the SSA simulations. The long-term coexistence behavior is exemplified in Figure 3 and Appendix 1—figures 8–10, involving simulations with or without stochasticity. The number of consumer species in long-term coexistence can be up to hundreds or more (see Figure 3 and Appendix 1—figure 8). To mimic real ecosystems, we further analyze cases with more than one type of resources, such as systems with $S_{R} = 3$ ( $S_{C} ≫ S_{R}$ ). Just like the case of $S_{R} = 1$ ( $S_{C} ≫ S_{R}$ ), an extensive range of consumer species may coexist indefinitely regardless of stochasticity (see Figure 3 and Appendix 1—figures 11–14).

Figure 3

Download asset Open asset

Intraspecific interference enables a wide range of consumer species to coexist with only one or a handful of resource species.

(**A, B**) Representative time courses simulated with ODEs and SSA. (**C, D**) A model of intraspecific predator interference illustrates the S-shaped pattern of the species’ rank-abundance curves across different ecological communities. The solid icons represent the experimental data (marked with ‘Exp’) reported in existing studies (Fuhrman et al., 2008; Cody and Smallwood, 1996; Terborgh et al., 1990; Martínez et al., 2023; Clarke et al., 2005; De Vries, 1997; Hubbell, 2001), where the bird community data were collected longitudinally in 1982 and 2018 (Terborgh et al., 1990; Martínez et al., 2023). The ODEs and SSA results were constructed from timestamp $t = 1.0 \times 10^{5}$ in the time series. In the K-S test, the probabilities (p-values) that the simulation results and the corresponding experimental data come from the same distributions are: $p_{ODEs}^{bird(1982)} = 0.17$ , $p_{ODEs}^{bird(2018)} = 0.26$ , $p_{ODEs}^{butterfly} = 0.70$ , $p_{ODEs}^{fish} = 0.88$ , $p_{ODEs}^{bat} = 0.42$ , $p_{SSA}^{bat} = 0.48$ , $p_{ODEs}^{lizard} = 0.96$ , $p_{SSA}^{lizard} = 0.54$ , $p_{ODEs}^{plankton} = 0.20$ , $p_{SSA}^{plankton} = 0.06$ . See Appendix 9 for simulation details and the Shannon entropies.

We further analyze the scenario involving chasing pairs and both intra- and interspecific interference, where multiple consumer species compete for one resource species. Similar to the scenario involving chasing pairs and intraspecific interference, all species coexist indefinitely in either ODEs or SSA simulation studies (see Appendix 1—figure 5F–H for the cases of $S_{C} = 6, 20$ and $S_{R} = 1$ ).

Intuitive understanding: an underlying negative feedback loop

For the case with only one resource species ( $S_{R} = 1$ ), if the total population size of the resources is much larger than that of the consumers (i.e. $R ≫ \sum_{i = 1}^{S_{C}} C_{i}$ ), the functional response $F \equiv k_{i} x_{i} / C_{i}$ and the steady-state population of each consumer and resource species can be obtained analytically (see Appendix 4.B-C for details). In fact, the functional response of a consumer species (e.g. $C_{i}$ ) is negatively correlated with its own population size:

F (R, C_{i}) \approx \frac{2 R}{\sqrt{(R + K_{i})^{2} + 8 β_{i} K_{i}^{2} C_{i}} + R + K_{i}},

where $β_{i} \equiv a_{i}^{'} / d_{i}^{'}$ . The analytical steady-state solutions are highly consistent with the numerical results (see Figure 1E and Appendix 1—figure 3H, I) and can even quantitatively predict the coexistence region of the parameter space (see Appendix 1—figure 3I).

Intuitively, the mechanisms of how intraspecific interference facilitates species coexistence can be understood from the underlying negative feedback loop. Specifically, for consumer species of higher competitiveness (e.g. $C_{i}$ ) in an ecological community, as the population size of $C_{i}$ increases during competition, a larger portion of $C_{i}$ individuals are then engaged in intraspecific interference pairs which are temporarily absent from hunting (see Equation S59 and Appendix 1—figure 15A, B). Consequently, the fraction of $C_{i}$ individuals within chasing pairs decreases (see Equation S59 and Appendix 1—figure 15A, B) and thus form a self-inhibiting negative feedback loop through the functional response (see Equation 5 and Appendix 1—figure 15C). This negative feedback loop prevents further increases in $C_{i}$ populations, results in an overall balance among the consumer species, and thus promotes biodiversity (see Appendix 4.C for details).

The S shape pattern of the rank-abundance curves in a broad range of ecological communities

As mentioned above, a prominent feature of biodiversity is that the species’ rank-abundance curves follow a universal S-shaped pattern in the linear-log plot across a broad spectrum of ecological communities (Fuhrman et al., 2008; Ser-Giacomi et al., 2018; Cody and Smallwood, 1996; Terborgh et al., 1990; Martínez et al., 2023; Clarke et al., 2005; Hubbell, 2001; De Vries, 1997). Previously, this pattern was mostly explained by the neutral theory (Hubbell, 2001), which requires special parameter settings that all consumer species share identical fitness. To resolve this issue, we apply the model involving chasing pairs and intraspecific interference to simulate the ecological communities, where one or three types of resources support a large number of consumer species ( $S_{C} ≫ S_{R}$ ). In each model system, the mortality rates of consumer species follow a Gaussian distribution where the coefficient of variation was taken around 0.3 (Menon et al., 2003; see Appendix 9 for details). For a broad array of the ecological communities, the rank-abundance curves obtained from the long-term coexisting state of both the ODEs and SSA simulation studies agree quantitatively with those of experiments (see Figure 3C and D and Appendix 1—figures 8–14), sharing roughly equal Shannon entropies and mostly being regarded as identical distributions in the Kolmogorov-Smirnov (K-S) statistical test (with a significance threshold of 0.05). Still, there is a noticeable discrepancy between the experimental data and SSA studies in terms of the species’ absolute abundances (e.g. see Appendix 1—figure 8C): those with experimental abundances less than 10 tend to be extinct in the SSA simulations. This is due to the fact that the recorded individuals in an experimental sample are just a tiny portion of that in the real ecological system, whereas the species population size in a natural community is certainly much larger than 10.

Discussion

The conflict between the CEP and biodiversity, exemplified by the paradox of the plankton (Hutchinson, 1961), is a long-standing puzzle in ecology. Although many mechanisms have been proposed to overcome the limit set by CEP (Hutchinson, 1961; Chesson, 2000; Levins, 1979; Levin, 1974; Koch, 1974; Huisman and Weissing, 1999; Czárán et al., 2002; Goyal and Maslov, 2018; Goldford et al., 2018; Villa Martín et al., 2020; Gupta et al., 2021; Thingstad, 2000; Wang and Liu, 2020; Dalziel et al., 2021; Posfai et al., 2017; Weiner et al., 2019; Xue and Goldenfeld, 2017; Beddington, 1975; DeAngelis et al., 1975; Arditi and Ginzburg, 1989; Kelsic et al., 2015; Grilli et al., 2017; Ratzke et al., 2020), it is still unclear how plankton and many other organisms can flout CEP and maintain biodiversity in quasi-well-mixed natural ecosystems. To address this issue, we investigate a mechanistic model with detailed consideration of pairwise encounters. Using numerical simulations combined with mathematical analysis, we identify that the intraspecific interference among the consumer individuals can promote a wide range of consumer species to coexist indefinitely with only one or a handful of resource species through the underlying negative feedback loop. By applying the above analysis to real ecological systems, our model naturally explains two classical experiments that reject CEP (Ayala, 1969; Park, 1954), and quantitatively illustrates the universal S-shaped pattern of the rank-abundance curves for a broad range of ecological communities (Fuhrman et al., 2008; Ser-Giacomi et al., 2018; Cody and Smallwood, 1996; Terborgh et al., 1990; Martínez et al., 2023; Clarke et al., 2005; Hubbell, 2001; De Vries, 1997).

In fact, predator interference has been introduced long ago by the classical B-D phenomenological model (Beddington, 1975; DeAngelis et al., 1975). However, the functional response of the B-D model involving intraspecific interference can be formally derived from the scenario involving only chasing pairs without consideration of pairwise encounters between consumer individuals (Wang and Liu, 2020; Huisman and De Boer, 1997; see Equations S8a and S24a). Yet, it has been demonstrated that the scenario involving only chasing pairs is under the constraint of CEP (Wang and Liu, 2020; see Appendix 1—figure 3A–C). Therefore, it is questionable regarding the validity of applying the B-D model to break CEP (Cantrell et al., 2004; Hsu et al., 2013). From a mechanistic perspective, we resolve these issues and show that the B-D model corresponds to a special case of our mechanistic model yet without the escape rate (see Appendix 1—figure 2 and Appendix 3 for details).

Our model is broadly applicable to explain biodiversity in many ecosystems. In practice, many more factors are potentially involved, and special attention is required to disentangle confounding factors. In microbial systems, complex interactions are commonly involved (Goyal and Maslov, 2018; Goldford et al., 2018; Hu et al., 2022), and species’ preference for food is shaped by the evolutionary course and environmental history (Wang et al., 2019). It is still highly challenging to fully explain how organisms evolve and maintain biodiversity in diverse ecosystems.

Materials and methods

Derivation of the encounter rates with the mean-field approximation

Request a detailed protocol

In the model depicted in Figure 1A, consumers and resources move randomly in space, which can be regarded as Brownian motions. At moment $t$ , a consumer individual of species $C_{i}$ moves at speed $v_{C_{i}}$ with velocity $v_{C_{i}} (t)$ , while a resource individual of species $R_{l}$ moves at speed $v_{R_{l}}$ with velocity $v_{R_{l}} (t)$ . Here, $v_{C_{i}}$ and $v_{R_{l}}$ are two invariants, while the directions of $v_{C_{i}} (t)$ and $v_{R_{l}} (t)$ change constantly. The relative velocity between the two individuals is $u_{C_{i} - R_{l}} (t) \equiv v_{R_{l}} (t) - v_{C_{i}} (t)$ , with a relative speed of $u_{C_{i} - R_{l}} (t)$ . Then, $u_{C_{i} - R_{l}} (t)^{2} = v_{C_{i}}^{2} + v_{R_{l}}^{2} - 2 v_{C_{i}} \cdot v_{R_{l}} \cdot \cos θ_{C_{i} - R_{l}} (t)$ , where $θ_{C_{i} - R_{l}} (t)$ represents the angle between $v_{C_{i}} (t)$ and $v_{R_{l}} (t)$ . This system is homogeneous, thus, $\bar{\cos θ_{C_{i} - R_{l}}} = 0$ , where the overline stands for the temporal average. Then, we obtain the average relative speed between the $C_{i}$ and $R_{l}$ individuals: $\bar{u_{C_{i} - R_{l}}} = \sqrt{v_{C_{i}}^{2} + v_{R_{l}}^{2}}$ . Likewise, the average relative speed between the $C_{i}$ and $C_{j}$ individuals is $\bar{u_{C_{i} - C_{j}}} = \sqrt{v_{C_{i}}^{2} + v_{C_{j}}^{2}}$ . Evidently, $\bar{u_{C_{i} - C_{i}}} = \sqrt{2} v_{C_{i}}$ . Meanwhile, the concentrations of species $C_{i}$ and $R_{l}$ in a 2-D square system with a length of $L$ are $n_{C_{i}} = C_{i} / L^{2}$ and $n_{R_{l}} = R_{l} / L^{2}$ , while those of the freely wandering $C_{i}$ and $R_{l}$ are $n_{C_{i}^{(F)}} = C_{i}^{(F)} / L^{2}$ and $n_{R_{l}^{(F)}} = R_{l}^{(F)} / L^{2}$ .

Then, we use the mean-field approximation to calculate the encounter rates $a_{i l}$ and $a_{i j}^{'}$ in the well-mixed system. In particular, we estimate $a_{i l}$ by tracking a randomly chosen consumer individual from species $C_{i}$ and counting its encounter frequency with the freely wandering individuals from resource species $R_{l}$ (see Appendix 1—figure 1). At any moment, the consumer individual may form a chasing pair with a $R_{l}$ individual within a radius of $r_{i l}^{(C)}$ (see Figure 1A). Over a time interval of $Δ t$ , the number of encounters between the consumer individual and $R_{l}$ individuals can be estimated by the encounter area and the concentration $n_{R_{l}}$ , which takes the value of $2 r_{i l}^{(C)} n_{R^{(F)}} \bar{u_{C_{i} - R}} Δ t$ (see Appendix 1—figure 1). Combined with $n_{R_{l}^{(F)}} = R_{l}^{(F)} / L^{2}$ , for all freely wandering $C_{i}$ individuals, the number of their encounters with $R^{(F)}$ during interval $Δ t$ is $\frac{2 r_{i l}^{(C)} \bar{u_{C_{i} - R}} C_{i}^{(F)} R^{(F)}}{L^{2}} Δ t$ . Meanwhile, in the ODEs, this corresponds to $a_{i} C_{i}^{(F)} R^{(F)} Δ t$ . Comparing both terms above, for chasing pairs, we have $a_{i l} = 2 r_{i l}^{(C)} L^{- 2} \bar{u_{C_{i} - R_{l}}} = 2 r_{i l}^{(C)} L^{- 2} \sqrt{v_{C_{i}}^{2} + v_{R_{l}}^{2}}$ . Likewise, for interference pairs, we obtain $a_{i j}^{'} = 2 r_{i j}^{(I)} L^{- 2} \bar{u_{C_{i} - C_{j}}} = 2 r_{i j}^{(I)} L^{- 2} \sqrt{v_{C_{i}}^{2} + v_{C_{j}}^{2}}$ . In particular, $a_{i i}^{'} = 2 \sqrt{2} v_{C_{i}} r_{i i}^{(I)} L^{- 2}$ .

Stochastic simulations

To investigate the impact of stochasticity on species coexistence, we use the stochastic simulation algorithm (SSA; Gillespie, 2007) and individual-based modeling (IBM; Vetsigian, 2017; Grimm, 2013) in simulating the stochastic process. In the SSA studies, we follow the standard Gillespie algorithm and simulation procedures (Gillespie, 2007).

In the IBM studies, we consider a 2D square system with a length of $L$ and periodic boundary conditions. In the case of $S_{C} = 2$ and $S_{R} = 1$ , consumers of species $C_{i}$ move at speed $v_{C_{i}}$ , while the resources move at speed $v_{R}$ . The unit length is $Δ l = 1$ , and all individuals move probabilistically. Specifically, when $Δ t$ is small so that $v_{C_{i}} Δ t ≪ 1$ , $C_{i}$ individuals jump a unit length with the probability $v_{C_{i}} Δ t$ . In practice, we simulate the temporal evolution of the model system following the procedures below.

Initialization

Request a detailed protocol

We choose the initial position for each individual randomly from a uniform distribution in the square space, which rounds to the nearest integer point in the $x$ - $y$ coordinates.

Moving

Request a detailed protocol

We choose the destination of a movement randomly from four directions ( $x$ -positive, $x$ -negative, $y$ -positive, $y$ -negative) following a uniform distribution. The consumers and resources jump a unit length with probabilities $v_{C_{i}} Δ t$ and $v_{R} Δ t$ , respectively.

Forming pairs

Request a detailed protocol

When a $C_{i}$ individual and a resource individual get close in space within a distance of $r^{(C)}$ , they form a chasing pair. Meanwhile, when two consumer individuals $C_{i}$ and $C_{j}$ stand within a distance of $r^{(I)}$ , they form an interference pair.

Dissociation

Request a detailed protocol

We update the system with a small time step $Δ t$ so that $d_{i} Δ t, k_{i} Δ t, d_{i j}^{'} Δ t ≪ 1$ . In practice, a random number $ς$ is sampled from a uniform distribution between 0 and 1, that is $U (0, 1)$ . If $ς < d_{i} Δ t$ or $ς < d_{i j}^{'} Δ t$ , then the chasing pair or interference pair dissociates into two separated individuals. One occupies the original position, while the other individual moves just out of the encounter radius in a uniformly distributed random angle. For a chasing pair, if $d_{i} Δ t < ς < (d_{i} + k_{i}) Δ t$ , then, the consumer catches the resource, and the biomass of the resource flows into the consumer populations (updated according to the birth procedure), while the consumer individual occupies the original position. Finally, if $ς > (d_{i} + k_{i}) Δ t$ or $ς > d_{i j}^{'} Δ t$ , the chasing pair or interference pair maintains the current status.

Birth and death

Request a detailed protocol

For each species, we use two separate counters with decimal precision to record the contributions of the birth and death processes, both of which accumulate in each time step. The incremental integer part of the counter will trigger updates in this run. Specifically, a newborn would join the system following the initialization procedure in a birth action, while an unfortunate target would be randomly chosen from the living population in a death action.

Appendix 1

Appendix 1—table 1

Illustrations of symbols in our generic model of pairwise encounters.

Symbols	Illustrations
$C_{i}$	The total population of consumer species $C_{i}$ .
$R_{l}$	The total population of resource species $R_{l}$ .
$C_{i}^{(F)}$	The freely wandering population of consumer species $C_{i}$ .
$R_{l}^{(F)}$	The freely wandering population of resource species $R_{l}$ .
$x_{i l}$	Chasing pairs formed between individuals from species $C_{i}$ and $R_{l}$ , i.e. $C_{i}^{(P)} ⋁ R_{l}^{(P)}$ .
$y_{i}$	Intraspecific interference pairs formed between individuals from species $C_{i}$ , i.e. $C_{i}^{(P)} ⋁ C_{i}^{(P)}$ .
$z_{i j}$	Interspecific interference pairs formed between individuals from species $C_{i}$ and $C_{j}$ , i.e., $C_{i}^{(P)} ⋁ C_{j}^{(P)}$ .
$r_{i l}^{(C)}$	The upper distance criterion for forming a chasing pair.
$r_{i j}^{(I)}$	The upper distance criterion for forming an interference pair.
$v_{C_{i}}$	The motility speed of consumer species $C_{i}$ .
$v_{R_{l}}$	The motility speed of resource species $R_{l}$ .
$S_{C}$	The number of consumer species.
$S_{R}$	The number of resource species.
$a_{i l}$	The encounter rate between a consumers and a resource.
$d_{i l}$	The escape rate within a chasing pair.
$k_{i l}$	The capture rate within a chasing pair.
$a_{i j}^{'}$	The encounter rate among consumer individuals.
$d_{i j}^{'}$	The separation rate within an interference pair.
$w_{i l}$	The mass conversion ratio from resource $R_{l}$ to $C_{i}$ .
$D_{i}$	The mortality rate of species $C_{i}$ .
$κ_{l}$	The steady-state population abundance of resources species $R_{l}$ in the absence of consumers.
$ζ_{l}$	The external resource supply rate of species $R_{l}$ .
$η_{l}$	The intrinsic growth rate of species $R_{l}$ for biotic resources (unused in all analyses).
$g_{l}$	The function describing the population dynamics of resource species $R_{l}$ .
$L$	The length of the 2-D square system where species coexist.
$C_{i}^{(F)} (+)$	We count $C^{(F)} (+)$ as $C^{(F)}$ , where ”(+)” signifies gaining biomass from resources.
$v_{C_{i}}$	The velocity of an individual of species $C_{i}$ .
$v_{R_{l}}$	The velocity of an individual of species $R_{l}$ .
$θ_{C_{i} - R_{l}}$	The angle between $v_{C_{i}}$ and $v_{R_{l}}$ .
$u_{C_{i} - R_{l}}$	The relative velocity between a consumer and a resource.
$u_{C_{i} - R_{l}}$	The relative speed between a consumer and a resource.
$n_{C_{i}}$	The concentration of species $C_{i}$ .
$n_{R_{l}}$	The concentration of species $R_{l}$ .
$n_{C_{i}^{(F)}}$	The concentration of the freely wandering $C_{i}$ .
$n_{R_{l}^{(F)}}$	The concentration of the freely wandering $R_{l}$ .

For all the symbols in Appendix 1—tables 1 and 2, the subscript ‘ $l$ ’ is omitted if $S_{R} = 1$ , and the subscript ‘ $i$ ’ is omitted if $S_{C} = 1$ .

Appendix 1—table 2

Illustrations of other symbols used in our manuscript.

Symbols	Illustrations/Definitions
$F$	The functional response.
$Ξ$	The searching efficiency.
$ς$	A random number sampled from a uniform distribution.
$U$	The uniform distribution.
$N (μ, σ)$	A Gaussian distribution with a mean of μ and a standard deviation of σ.
$\equiv$	An equal sign for equations defining the symbol on the left-hand side.
$Δ$	The competitive difference between two consumer species, defined as $Δ \equiv (D_{1} - D_{2}) / D_{2}$ .
$\hat{Δ}$	The supremum of the competitive difference tolerated for species coexistence.
$Δ t$	A short time interval.
$P_{i}$	The probability that a consumer individual of the ecological community belongs to species $C_{i}$ .
$p_{O D E s}$ , $p_{S S A}$	The p-value assessing the similarity of simulation results and experimental data.
$H$	The Shannon entropy: $H = - \sum_{i = 1}^{S_{C}} P_{i} \log_{2} (P_{i})$ .
$\bar{X}$	The time average of an arbitrary quantity $X$ .
$δ X$	The standard deviation of an arbitrary quantity $X$ .
$⟨ X ⟩$	The expectation of a random variable $X$ .
$τ$	The parameter that sets the time scale of a system.
$K_{i l}$	$K_{i l} \equiv \frac{k_{i l} + d_{i l}}{a_{i l}}$ .
$α_{i l}$	$α_{i l} \equiv \frac{D_{i}}{w_{i l} k_{i l}}$ .
$β_{i}$	$β_{i} \equiv \frac{a_{i}^{'}}{d_{i}^{'}}$ .
$γ$	$γ \equiv \frac{a_{12}^{'}}{d_{12}^{'}}$ .
$f_{i} (R^{(F)})$	$f_{i} (R^{(F)}) \equiv R^{(F)} / (R^{(F)} + K_{i})$ .
$ϕ_{0}$ , $ϕ_{1}$ , $ϕ_{2}$	$ϕ_{0} \equiv - C R^{2}$ , $ϕ_{1} \equiv 2 C R + K R + R^{2}$ , $ϕ_{2} \equiv 2 β K^{2} - K - C - 2 R$ .
$Λ$	The discriminant of Equation S13.
$ψ$ , $φ$	$ψ \equiv ϕ_{1} - (ϕ_{2})^{2} / 3$ , $φ \equiv ϕ_{0} - ϕ_{1} ϕ_{2} / 3 + 2 (ϕ_{2})^{3} / 27$ .
$ω$ , $θ_{1}$ , $θ_{2}$	$ω \equiv - 1 / 2 + i \sqrt{3} / 2$ , $θ_{1} \equiv (- φ / 2 + \sqrt{- Λ / 108})^{1 / 3}$ , $θ_{2} \equiv (- φ / 2 - \sqrt{- Λ / 108})^{1 / 3}$ .
$ψ^{'}$ , $φ^{'}$	$ψ^{'} \equiv (- 4 ψ / 3)^{1 / 2}$ , $φ^{'} \equiv \arccos (- (- ψ / 3)^{- 3 / 2} φ / 2) / 3$ .
$t_{h}$	The handling time in the B-D model.
$t_{w}$	The wasting time in the B-D model.
$u_{i} (R, C_{1}, C_{2})$	Expression of $x_{i}$ using $C_{1}$ , $C_{2}$ , and $R$ in Equation S33 involving intraspecific interference, see Equation S37.
$Ω_{i} (R, C_{1}, C_{2})$	$Ω_{i} (R, C_{1}, C_{2}) \equiv \frac{w_{i} k_{i}}{C_{i}} u_{i} (R, C_{1}, C_{2})$ .
$G (R, C_{1}, C_{2})$	$G (R, C_{1}, C_{2}) \equiv g (R, u_{1} (R, C_{1}, C_{2}), u_{2} (R, C_{1}, C_{2}), C_{1}, C_{2})$ , see Equations S33 and S38.
$o_{1}$ , $o_{2}$	$o_{1} \equiv \frac{ζ}{κ} - \frac{k_{1}}{2 β_{1} K_{1}} - \frac{k_{2}}{2 β_{2} K_{2}}$ , $o_{2} \equiv \frac{k_{1} (1 - α_{1})}{2 β_{1} α_{1} (K_{1})^{2}} + \frac{k_{2} (1 - α_{2})}{2 β_{2} α_{2} (K_{2})^{2}}$ .
$ϖ$	$ϖ \equiv \frac{1}{2} (\frac{1}{κ} - \frac{k_{2}}{2 ζ β_{2} K_{2}}) + \frac{1}{2} \sqrt{(\frac{1}{κ} - \frac{k_{2}}{2 ζ β_{2} K_{2}})^{2} + 2 \frac{k_{2} (1 - α_{2})}{ζ β_{2} α_{2} (K_{2})^{2}}}$ .
$ι_{1}$ , $ι_{2}$	$ι_{1} \equiv \frac{ζ}{κ} - \sum_{i = 1}^{S_{C}} \frac{k_{i}}{2 β_{i} K_{i}}$ , $ι_{2} \equiv \sum_{i = 1}^{S_{C}} \frac{k_{i} (1 - α_{i})}{2 β_{i} α_{i} (K_{i})^{2}}$ .
$u_{i}^{'} (R, C_{1}, C_{2})$	Expression of $x_{i}$ using $C_{1}$ , $C_{2}$ , and $R$ in Equation S61 involving interspecific interference, see Equation S64.
$Ω_{i}^{'} (R, C_{1}, C_{2})$	$Ω_{i}^{'} (R, C_{1}, C_{2}) \equiv \frac{w_{i} k_{i}}{C_{i}} u_{i}^{'} (R, C_{1}, C_{2}) .$
$G^{'} (R, C_{1}, C_{2})$	$G^{'} (R, C_{1}, C_{2}) \equiv g (R, u_{1}^{'} (R, C_{1}, C_{2}), u_{2}^{'} (R, C_{1}, C_{2}), C_{1}, C_{2})$ , see Equation S61 and S65.
$ϱ_{1}$ , $ϱ_{2}$	$ϱ_{1} \equiv \frac{ζ}{κ} - \frac{k_{1}}{γ K_{1}} - \frac{k_{2}}{γ K_{2}}$ , $ϱ_{2} \equiv \frac{k_{1} (1 - α_{2})}{γ K_{1} K_{2} α_{2}} + \frac{k_{2} (1 - α_{1})}{γ K_{1} K_{2} α_{1}}$ .
$χ_{1}$	$χ_{1} \equiv \frac{k_{2} γ (α_{1} - 1)}{K_{1} K_{2} α_{1} (4 β_{1} β_{2} - γ^{2})} + \frac{k_{1} γ (α_{2} - 1)}{K_{1} K_{2} α_{2} (4 β_{1} β_{2} - γ^{2})} - \frac{k_{1} 2 β_{2} (α_{1} - 1)}{K_{1}^{2} α_{1} (4 β_{1} β_{2} - γ^{2})} - \frac{k_{2} 2 β_{1} (α_{2} - 1)}{K_{2}^{2} α_{2} (4 β_{1} β_{2} - γ^{2})}$ .
$χ_{2}$	$χ_{2} \equiv \frac{k_{1} (γ - 2 β_{2})}{K_{1} (4 β_{1} β_{2} - γ^{2})} + \frac{k_{2} (γ - 2 β_{1})}{K_{2} (4 β_{1} β_{2} - γ^{2})} + \frac{ζ}{κ}$ .

Appendix 1—figure 1

Download asset Open asset

Estimation of the encounter rates with the mean-field approximations.

To calculate $a_{i l}$ in the chasing pair, we suppose that all individuals of species $R_{l}$ stand still while a $C_{i}$ individual moves at the speed of $\bar{u_{C_{i} - R_{l}}}$ (the relative speed). Over a time interval of $Δ t$ , the length of zigzag trajectory of the $C_{i}$ individual is approximately $\bar{u_{C_{i} - R_{l}}} Δ t$ , while the encounter area (marked with dashed lines) is estimated to be $2 r_{i l}^{(C)} \bar{u_{C_{i} - R}} Δ t$ . Then, we can estimate the encounter rate $a_{i l}$ using the encounter area and concentrations of the species (see Materials and methods for details). Similarly, we can estimate $a_{i j}^{'}$ in the interference pair.

Appendix 1—figure 2

Download asset Open asset

Functional response in scenarios involving different types of pairwise encounter.

(**A–C**) In the scenario involving only chasing pair, the red surface/line corresponds to the B-D model (calculated from Equation S10b), while the green surface/line represents the exact solutions to our mechanistic model (calculated from Equation S7b). The magenta (calculated from Equation S8b) and blue (calculated from Equation S9b) surfaces/lines represent the approximate solutions to our model (see Appendix 3.B). (**D–F**) In the scenario involving chasing pairs and intraspecific interference, the red surface/line corresponds to the B-D model (calculated from Equation S24b), while the green/line surface represents the exact solutions to our mechanistic model (calculated from Equation S17b). The blue surface/line (calculated from Equation S19b) and the magenta surface/line (calculated from Equation S21b) represent the quasi-rigorous and the approximate solutions to our model, respectively (see Appendix 3.C). (**G–I**) In the scenario involving chasing pairs and interspecific interference, the red surface/line corresponds to the B-D model (calculated from Equation S31b), while the green surface/line represents the quasi-rigorous solutions to our mechanistic model (calculated from Equation S27d). The blue surface/line (calculated from Equation S29d) represents the approximate solutions to our model (see Appendix 3.D). In (**A–C**): $k = 0.1$ , $a = 0.25$ . In (A): $d = 0$ . In (C): $R = 10^{4}$ , $C = 10^{3}$ . In (**D–F**): $a = 0.1$ , $k = 0.1$ , $a^{'} = 0.12$ , $d^{'} = 0.1$ . In (D): $d = 0$ . In (F): $R = 10^{6}$ , $C = 10^{5}$ . In (**G–I**): $a_{1} = a_{2} = 0.1$ , $k_{1} = k_{2} = 0.1$ , $a_{12}^{'} = 0.12$ , $d_{12}^{'} = 0.1$ . In (G): $d_{i} = 0$ . In (I): $R = 10^{6}$ , $C_{i} = 10^{5}$ .

Appendix 1—figure 3

Download asset Open asset

Numerical solutions in scenarios involving chasing pairs and different types of predator interference.

Here, $S_{C} = 2$ and $S_{R} = 1$ . $D_{i} (i = 1, 2)$ is the only parameter varying with the consumer species (with $D_{1} > D_{2}$ ), and $Δ \equiv (D_{1} - D_{2}) / D_{2}$ represents the competitive difference between the two species. (**A–C**) Scenario involving only chasing pairs. (A) If all consumer species coexist at steady state, then $f_{i} (R^{(F)}) / D_{i} = 1$ , where $f_{i} (R^{(F)}) = R^{(F)} / (R^{(F)} + K_{i})$ and $K_{i} = (d_{i} + k_{i}) / a_{i}$ . This requires that the three lines $y = f_{i} (R) / D_{i}$ and $y = 1$ share a common point, which is generally impossible. (B) The blue plane is parallel to the green one, and hence they do not have a common point. (C) Time courses of the species abundances in the scenario involving only chasing pairs. The two consumer species cannot coexist at steady state. (**D–F**) Scenario involving chasing pairs and interspecific interference. (**G–I**) Scenario involving chasing pairs and intraspecific interference. (**D, G**) Positive solutions to the steady-state equations: $\dot{R} = 0$ (orange surface), ${\dot{C}}_{1} = 0$ (blue surface), ${\dot{C}}_{2} = 0$ (green surface). The intersection point marked by black/red dots is an unstable/stable fixed point. (**E, H**) Comparisons between the numerical results and analytical solutions of the species abundances at fixed points. Color bars are analytical solutions while hollow bars are numerical results. The analytical solutions in (E) and (H) (marked with superscript ‘(A)’) were calculated from Equations S68 and S70 and Equation S41, Equation S43, respectively. (F) In this scenario, there is no parameter region for stable coexistence. The region below the red surface and above $Δ = 0$ represents unstable fixed points. (I) Comparisons between the numerical results and analytical solutions of the coexistence region. Here, $\hat{Δ}$ represents the maximum competitive difference tolerated for species coexistence. The red and cyan surfaces represent the analytical solutions (calculated from Equation S46) and numerical results, respectively. The numerical results in (C), (**D–F**) and (**G–I**) were calculated from Equations 1 and 4, Equation S42 and S61 and Equation S33 and S42, respectively. In (C): $a_{i} = 0.1$ , $k_{i} = 0.1$ , $w_{i} = 0.1$ , $d_{i} = 0.5$ , $(i = 1, 2)$ , $D_{1} = 0.002$ , $D_{2} = 0.001$ , $κ = 5$ , $ζ = 0.05$ . In (D): $a_{i} = 0.05$ , $d_{i} = 0.05$ , $k_{i} = 0.02$ , $w_{i} = 0.08$ , $(i = 1, 2)$ , $D_{1} = 0.0011$ , $D_{2} = 0.001$ , $κ = 20$ , $ζ = 0.01$ , $a_{12}^{'} = 0.06$ , $d_{12}^{'} = 0.01$ . In (E): $a_{i} = 0.04$ , $d_{i} = 0.2$ , $k_{i} = 0.1$ , $w_{i} = 0.3$ , $(i = 1, 2)$ , $D_{2} = 0.0008$ , $κ = 10$ , $ζ = 0.2$ , $a_{12}^{'} = 0.048$ , $d_{12}^{'} = 0.001$ . In (F): $a_{i} = 0.1$ , $k_{i} = 0.1$ , $w_{i} = 0.1$ , $(i = 1, 2)$ , $D_{2} = 0.001$ , $κ = 100$ , $ζ = 0.05$ , $a_{12}^{'} = 0.12$ . In (G): $a_{i} = 0.5$ , $a_{i}^{'} = 0.625$ , $d_{i} = 0.5$ , $d_{i}^{'} = 0.5$ , $k_{i} = 0.4$ , $w_{i} = 0.5$ , $(i = 1, 2)$ , $D_{2} = 0.02$ , $D_{1} = 1.2 D_{2}$ , $κ = 10$ , $ζ = 0.1$ . In (H): $a_{i} = 0.1$ , $a_{i}^{'} = 0.12$ , $k_{i} = 0.12$ , $w_{i} = 0.3$ , $d_{i} = 0.5$ , $d_{i}^{'} = 0.05$ , $(i = 1, 2)$ , $D_{2} = 0.02$ , $κ = 100$ , $ζ = 0.8$ . In (I): $a_{i} = 0.5$ , $k_{i} = 0.2$ , $d_{i} = 0.8$ , $w_{i} = 0.2$ , $(i = 1, 2)$ , $D_{2} = 0.008$ , $κ = 60$ , $ζ = 0.8$ .

Appendix 1—figure 4

Download asset Open asset

Intraspecific predator interference facilitates species coexistence regardless of stochasticity.

Here, we consider the case of $S_{C} = 2$ , $S_{R} = 1$ . (A) A representative trajectory of species coexistence in the phase space simulated with ODEs. The fixed point (shown in red) is stable and globally attractive. (**B, C**) 3D phase diagrams in the ODEs studies. Here, $D_{i}$ is the only parameter that varies with the two consumer species, and $Δ \equiv (D_{1} - D_{2}) / D_{2}$ measures the competitive difference between the two species. The parameter region below the blue surface yet above the red surface represents stable coexistence. The region below the red surface and above $Δ = 0$ represents unstable fixed points (an empty set). (D) An exemplified transection corresponding to the $Δ = 0.3$ plane in (C). (E) Time courses of the species abundances simulated with ODEs or SSA. (F) Representative trajectories of species coexistence in the phase space simulated with SSA. The coexistence state is stable and globally attractive (see (E) for time courses, SSA results). (**A–F**) were calculated or simulated from Equations 1, 2 and 4. In (A): $a_{i} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i} = 0.1$ , $d_{i}^{'} = 0.05$ , $w_{i} = 0.1$ , $k_{i} = 0.1$ , $(i = 1, 2)$ , $D_{1} = 0.0035$ , $D_{2} = 0.0038$ , $κ = 100$ , $ζ = 0.3$ . In (B): $a_{i} = 0.1$ , $d_{i} = 0.1$ , $w_{i} = 0.1$ , $k_{i} = 0.1$ , $(i = 1, 2)$ , $D_{2} = 0.001$ , $Δ \equiv (D_{1} - D_{2}) / D_{2}$ , $κ = 100$ , $ζ = 0.1$ . In (**C, D**): $a_{i} = 0.05$ , $a_{i}^{'} = 0.065$ , $w_{i} = 0.1$ , $k_{i} = 0.1$ , $(i = 1, 2)$ , $D_{2} = 0.002$ , $Δ \equiv (D_{1} - D_{2}) / D_{2}$ , $κ = 10$ , $s = 0.1$ . In (D): $Δ = 0.3$ . In (**E, F**): $a_{i} = 0.02$ , $a_{i}^{'} = 0.025$ , $d_{i} = 0.7$ , $d_{i}^{'} = 0.7$ , $w_{i} = 0.4$ , $k_{i} = 0.05$ , $(i = 1, 2)$ , $D_{1} = 0.0160$ , $D_{2} = 0.0171$ , $κ = 2000$ , $ζ = 5.5$ .

Appendix 1—figure 5

Download asset Open asset

Outcomes of multiple consumers species competing for one resource species involving chasing pairs and intra- and inter-specific interference.

(**A–E**) The case involving two consumer species and one resource species ( $S_{C} = 2$ , $S_{R} = 1$ ). Here, $D_{i}$ is the only parameter that varies with the consumer species (with $D_{1} > D_{2}$ ), and $Δ \equiv (D_{1} - D_{2}) / D_{2}$ measures the competitive difference between the two species. (A) A 3D phase diagram. The parameter region below the blue surface yet above the red surface represents stable coexistence, while that below the red surface and above $Δ = 0$ represents unstable fixed points (an empty set). (B) Time courses of the species abundances. Two consumer species may coexist with one type of resources at steady state. (C) Representative trajectories of species coexistence in the phase space. The fixed point (shown in red) is stable and globally attractive. (D) Consumer species may coexist indefinitely with the resources regardless of stochasticity. (E) Comparisons between numerical results and analytical solutions of the species abundances at fixed points. Color bars are analytical solutions while hollow bars are numerical results. The analytical solutions (marked with superscript ‘(A)’) were calculated from Equation S74 and S75. (**F–H**) Time courses of species abundances in cases involving 6 or 20 consumer species and one type of resources ( $S_{C} = 6$ or 20, $S_{R} = 1$ ). All consumer species may coexist with one type of resource at a steady state, and this coexisting state is robust to stochasticity. (**A–C, E, G**) ODEs results. (**D, F**) ODEs and SSA results. (H) SSA results. The numerical results in (**A–H**) were calculated or simulated from Equations 1-4. In (**A–C**): $a_{i} = 0.1$ , $a_{i}^{'} = 0.12$ , $k_{i} = 0.1$ , $w_{i} = 0.1$ , $(i = 1, 2)$ , $D_{2} = 0.004$ , $κ = 100$ , $ζ = 0.8$ , $a_{12}^{'} = 0.12$ , $d_{12}^{'} = 0.5$ . In (**B–C**): $d_{i} = 0.3$ , $d_{i}^{'} = 0.5$ , $(i = 1, 2)$ . In (D): $a_{i} = 0.1$ , $a_{i}^{'} = 0.14$ , $k_{i} = 0.12$ , $w_{i} = 0.15$ , $d_{i} = 0.3$ , $d_{i}^{'} = 0.5$ , $(i = 1, 2)$ , $D_{1} = 0.0125$ , $D_{2} = 0.012$ , $κ = 300$ , $ζ = 5.5$ , $a_{12}^{'} = 0.14$ , $d_{12}^{'} = 5$ . In (E): $a_{i} = 0.05$ , $a_{i}^{'} = 0.06$ , $k_{i} = 0.1$ , $w_{i} = 0.2$ , $d_{i} = 0.5$ , $d_{i}^{'} = 0.15$ , $(i = 1, 2)$ , $D_{2} = 0.008$ , $κ = 100$ , $ζ = 0.8$ , $a_{12}^{'} = 0.06$ , $d_{12}^{'} = 0.0005$ . In (F): $a_{i} = 0.1$ , $a_{i}^{'} = 0.14$ , $k_{i} = 0.15$ , $w_{i} = 0.18$ , $d_{i} = 2.8$ , $d_{i}^{'} = 0.02$ , $a_{i j}^{'} = 0.14$ , $d_{i j}^{'} = 0.15$ , $(i, j = 1, . . ., 6$ , $i \neq j)$ , $κ = 600$ , $s = 100$ , $D_{1} = 0.0091$ , $D_{2} = 0.0084$ , $D_{3} = 0.0088$ , $D_{4} = 0.0096$ , $D_{5} = 0.0082$ , $D_{6} = 0.0093$ . In (**G–H**): $a_{i} = 0.1$ , $a_{i}^{'} = 0.14$ , $k_{i} = 0.2$ , $w_{i} = 0.18$ , $d_{i} = 2.8$ , $d_{i}^{'} = 0.02$ , $a_{i j}^{'} = 0.14$ , $d_{i j}^{'} = 0.8$ , $D_{i} = N (1, 0.1) \times 0.005$ , $(i, j = 1, . . ., 20$ , $i \neq j)$ , $κ = 1000$ , $s = 500$ .

Appendix 1—figure 6

Download asset Open asset

The influence of stochasticity on species coexistence.

(**A, B**) Stochasticity jeopardizes species coexistence. Koch’s model (Koch, 1974) and Huisman-Weissing model (Huisman and Weissing, 1999) were simulated with SSA using the same parameter settings as their deterministic model. Nevertheless, both cases of oscillating coexistence are vulnerable to stochasticity. See (Koch, 1974) and (Huisman and Weissing, 1999) for the parameters. (**C, D**) Phase diagrams in the scenario involving chasing pairs and intraspecific interference. Here, $S_{C} = 2$ and $S_{R} = 1$ . $D_{i} (i = 1, 2)$ is the only parameter varying with the consumer species (with $D_{1} > D_{2}$ ), and $Δ \equiv (D_{1} - D_{2}) / D_{2}$ represents the competitive difference between the two species. (C) The ODEs results. (D) The SSA results (with the same parameter region as (C)). The species’ coexisting fraction in each pixel was calculated from 16 random repeats. (C) and (D) were calculated from Equations 1, 2 and 4. In (**C, D**): $a_{i} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i} = 0.5$ , $w_{i} = 0.1$ , $k_{i} = 0.1$ , $(i = 1, 2)$ , $κ = 100$ , $ζ = 5$ , $D_{2} = 0.0014$ .

Appendix 1—figure 7

Download asset Open asset

A model of intraspecific predator interference explains two classical experiments that invalidate CEP.

(A) In Ayala’s experiment (Ayala, 1969), two *Drosophila* species (consumers) coexisted for 40 weeks with the same type of abiotic resources within a laboratory bottle. The time averages ( $\bar{C_{i}}$ ) and standard deviation ( $δ C_{i}$ ) of the species’ relative abundances for the experimental data or SSA results are: $\bar{^{(R)} C_{D . s e r r a t a_C H_G r p 1}^{E x p (S S A)}} = 0.77 (0.80)$ , $δ^{(R)} C_{D . s e r r a t a_C H_G r p 1}^{Exp (SSA)} = 0.17 (0.08)$ , $\bar{^{(R)} C_{D . s e r r a t a_C H_G r p 2}^{Exp (SSA)}} = 0.78 (0.73)$ , $δ^{(R)} C_{D . s e r r a t a_C H_G r p 2}^{Exp (SSA)} = 0.14 (0.07)$ , where the superscript ‘(R)’ represents relative abundances. (B) In Park’s experiment (Park, 1954), two *Tribolium* species coexisted for 750 days with the same food (flour). The time averages ( $\bar{C_{i}}$ ) and standard deviations ( $δ C_{i}$ ) of the species’ abundances are: $\bar{C_{T . c o n f u s u m {_29}^{\circ} C}^{Exp (SSA)}} = 33.4 (28.8)$ , $δ C_{T . c o n f u s u m {_29}^{\circ} C}^{E x p (SSA)} = 6.0 (5.4)$ , $\bar{C_{T . c a s t a n e u m a {_29}^{\circ} C}^{E x p (SSA)}} = 48.8 (47.7)$ , $δ C_{T . c a s t a n e u m a {_29}^{\circ} C}^{E x p (SSA)} = 11.9 (9.9)$ . (**A, B**) The solid icons represent the experimental data, which are connected by the dotted lines for the sake of visibility. The solid lines stand for the SSA simulation results. (**C, D**) The Shannon entropies of each time point for the experimental or model-simulated communities shown in (B) and Figure 2D and E. Here, we calculated the Shannon entropies with $H (t) = - \sum_{i = 1}^{S_{C}} P_{i} (t) \log_{2} (P_{i} (t))$ , where $P_{i} (t)$ is the probability that a consumer individual belongs to species $C_{i}$ at the time stamp of $t$ . The time averages ( $\bar{H}$ ) and standard deviations ( $δ H$ ) of the Shannon entropies are: ${\bar{H}}_{E x p (S S A)}^{D r o s o p h i l a_A R_G r p 1} = 0.95 (0.97)$ , $δ H_{E x p (S S A)}^{D r o s o p h i l a_A R_G r p 1} = 0.06 (0.04)$ , ${\bar{H}}_{E x p (S S A)}^{D r o s o p h i l a_A R_G r p 2} = 0.94 (0.92)$ , $δ H_{E x p (S S A)}^{D r o s o p h i l a_A R_G r p 2} = 0.07 (0.07)$ , ${\bar{H}}_{E x p (S S A)}^{T r i b o l i u m {_24}^{\circ} C} = 0.96 (0.92)$ , $δ H_{E x p (S S A)}^{T r i b o l i u m {_24}^{\circ} C} = 0.02 (0.05)$ , ${\bar{H}}_{E x p (S S A)}^{T r i b o l i u m {_29}^{\circ} C} = 0.97 (0.94)$ , $δ H_{E x p (S S A)}^{T r i b o l i u m {_29}^{\circ} C} = 0.02 (0.05)$ . (**E–H**) Time courses of the species abundances in the scenario involving chasing pairs and intraspecific interference. The time series in (**E–H**) correspond to the long-term version of that shown in Figure 2D, Appendix 1—figure 7A, Figure 2E, Appendix 1—figure 7B, respectively. (**A, B, F–I**) were simulated from Equations 1, 2 and 4. In (A): $a_{i} = 0.3$ , $a_{i}^{'} = 0.33$ , $w_{i} = 0.018$ , $k_{i} = 4.8$ , $d_{i}^{'} = 5$ , $d_{i} = 5.5$ , $(i = 1, 2)$ , $D_{1} = 0.0132$ , $D_{2} = 0.010$ , $ζ = 35$ , $κ = 10000$ . In (B): $a_{i} = 0.3$ , $a_{i}^{'} = 0.36$ , $w_{i} = 0.02$ , $k_{i} = 4.5$ , $d_{i}^{'} = 4$ , $d_{i} = 4.5$ , $(i = 1, 2)$ , $D_{1} = 0.0122$ , $D_{2} = 0.010$ , $ζ = 35$ , $κ = 10000$ . In(A-B): $τ = 0.4$ Day (see Appendix 7).

Appendix 1—figure 8

Download asset Open asset

A model of intraspecific interference semi-quantitatively illustrates the rank-abundance curve of a plankton community ( $S_{C} ≫ S_{R}$ ).

(**A, B**) Intraspecific interference enables a wide range of consumers species to coexist with one type of resources. (A) Time courses of the species abundances simulated with ODEs. (B) Time series of the species abundances simulated with SSA (with theh same as parameter settings as (A)). (C) The rank-abundance curve of a plankton community. The solid dots represent the experimental data (marked with ‘Exp’) reported in a recent study (Ser-Giacomi et al., 2018) (TARA_139.SUR.180.2000.DNA), while the hollow dots and those with ‘+’ center are the ODEs and SSA results constructed from timestamp $t = 5.0 \times 10^{5}$ in the time series (see (A) and (B)), respectively. The Shannon entropies of the experimental data and simulation results for the plankton community are: $H_{Exp(ODEs,SSA)}^{Plankton} = 2.85 (2.18, 2.00)$ . In the model settings, $S_{C} = 200$ and $S_{R} = 1$ . $D_{i} (i = 1, \dots, S_{C})$ is the only parameter that varies with the consumer species, which was randomly drawn from a Gaussian distribution $N (μ, σ)$ . Here, μ and σ are the mean and standard deviation of the distribution. The numerical results in (**A–C**) were simulated from Equations 1, 2 and 4. In (**A–C**): $a_{i} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i} = 0.2$ , $d_{i}^{'} = 0.5$ , $w_{i} = 0.2$ , $k_{i} = 0.1$ , $D_{i} = 0.008 \times N (1, 0.38)$ , $(i = 1, \dots, 200)$ , $κ = 10^{5}$ , $ζ = 150$ .

$a_{i} = 0.1, a_{i}^{'} = 0.125, d_{i}^{'} = 0.5, d_{i} = 0.2, w_{i} = 0.2, k_{i} = 0.1, D_{i} = N (1, 0.38) \times 0.008$

Appendix 1—figure 9

Download asset Open asset

A model of intraspecific interference illustrates the rank-abundance curves across different ecological communities ( $S_{C} ≫ S_{R}$ ).

The solid dots represent the experimental data (marked with ‘Exp’) reported in existing studies (Hubbell, 2001; Holmes et al., 1986; Cody and Smallwood, 1996), while the hollow dots and those with ‘+’ center are the ODEs and SSA results constructed from timestamp $t = 1.0 \times 10^{5}$ in the time series (see Appendix 1—figure 10A–G), respectively. In the model settings, $S_{R} = 1$ , $S_{C} = 20$ (in (A)), 35 (in (B)) or 45 (in (C)). $D_{i}$ is the only parameter varying with the consumer species, which was randomly drawn from a Gaussian distribution. The Shannon entropies of the experimental data and simulation results for each ecological community are: $H_{Exp(ODEs,SSA)}^{bird} = 2.98 (3.06, 2.98)$ , $H_{Exp(ODEs,SSA)}^{bee} = 4.04 (4.02, 4.02)$ , $H_{Exp(ODEs,SSA)}^{fish} = 3.78 (3.40, 3.41)$ . In the Kolmogorov-Smirnov (K-S) test, the probabilities (p-values) that the simulation results and the corresponding experimental data come from identical distributions are: $p_{ODEs}^{bird} = 0.89$ , $p_{SSA}^{bird} = 0.88$ , $p_{ODEs}^{bee} = 0.47$ , $p_{SSA}^{bee} = 0.75$ , $p_{SSA}^{fish} = 0.77$ . With a significance threshold of 0.05, none of the p-values suggest there exists a statistically significant difference. The numerical results in (**A–C**) were simulated from Equations 1, 2 and 4. In (A): $a_{i} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i} = 0.5$ , $d_{i}^{'} = 0.6$ , $w_{i} = 0.22$ , $k_{i} = 0.1$ , $D_{i} = 0.016 \times N (1, 0.35)$ , $(i = 1, \dots, 20)$ , $ζ = 350$ , $κ = 10^{6}$ . In (B): $a_{i} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i} = 0.5$ , $d_{i}^{'} = 0.6$ , $w_{i} = 0.22$ , $k_{i} = 0.1$ , $D_{i} = 0.012 \times N (1, 0.35)$ , $(i = 1, \dots, 35)$ , $ζ = 350$ , $κ = 10^{6}$ . In (C): $a_{i} = 0.1$ , $a_{i}^{'} = 0.14$ , $d_{i} = 0.5$ , $d_{i}^{'} = 0.5$ , $w_{i} = 0.2$ , $k_{i} = 0.1$ , $D_{i} = 0.015 \times N (1, 0.32)$ , $(i = 1, \dots, 45)$ , $ζ = 550$ , $κ = 10^{6}$ .

Appendix 1—figure 10

Download asset Open asset

Time courses of the species abundances in the scenario involving chasing pairs and intraspecific interference.

The time series in (**A, E**), (**B, F**), (**C, G**), (**D, H**), (**I, J**) and (K) correspond to that shown in Appendix 1—figure 9A–C, Figure 3D (bat), Figure 3D (lizard) and Figure 3C (butterfly), respectively.

Appendix 1—figure 11

Download asset Open asset

Appendix 1—figure 12

Download asset Open asset

A model of intraspecific interference illustrates the rank-abundance curves across different plankton communities ( $S_{C} ≫ S_{R}$ ).

The solid dots represent the experimental data (marked with ‘Exp’) reported in a recent study (Fuhrman et al., 2008), while the hollow dots and those with ‘+’ center are the ODEs and SSA results constructed from timestamp $t = 1.0 \times 10^{5}$ in the time series (see Appendix 1—figure 13), respectively. The plankton community data were obtained separately from the Norwegian Sea (NS) and Pacific Station (PS). In the model settings, $S_{R} = 1$ (in (**B, C**)), 3 (in (A)); $S_{C} = 50$ (in (**A, C**)), 150 (in (B)). $D_{i}$ is the only parameter varying with the consumer species, which was randomly drawn from a Gaussian distribution. The Shannon entropies of the experimental data and simulation results for each plankton community are: $H_{Exp(ODEs,SSA)}^{plankton(NS)} = 4.67 (4.85, 4.90)$ for $S_{R} = 3$ , $H_{Exp(ODEs,SSA)}^{plankton(NS)} = 4.67 (4.74, 4.64)$ for $S_{R} = 1$ . In the K-S test, the p-values that the simulation results and the corresponding experimental data come from identical distributions are: $p_{ODEs}^{plankton(NS)} = 0.31$ , $p_{SSA}^{plankton(NS)} = 0.14$ for $S_{R} = 3$ , $p_{ODEs}^{plankton(NS)} = 0.46$ , $p_{SSA}^{plankton(NS)} = 0.37$ for $S_{R} = 1$ . With a significance threshold of 0.05, none of the p-values suggest there exists a statistically significant difference. The numerical results in (**A–C**) were simulated from Equations 1, 2 and 4. In (A): $a_{i l} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i l} = 0.5$ , $d_{i}^{'} = 0.2$ , $w_{i l} = 0.3$ , $k_{i l} = 0.2$ , $κ_{1} = 8 \times 10^{4}$ , $κ_{2} = 5 \times 10^{4}$ , $κ_{3} = 3 \times 10^{4}$ , $ζ_{1} = 280$ , $ζ_{2} = 200$ , $ζ_{3} = 150$ , $D_{i} = 0.035 \times N (1, 0.25)$ , $(i = 1, \dots, 50$ , $l = 1, 2, 3)$ . In (B): $a_{i} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i} = 0.3$ , $d_{i}^{'} = 0.3$ , $w_{i} = 0.3$ , $k_{i} = 0.2$ , $D_{i} = 0.025 \times N (1, 0.25)$ , $(i = 1, \dots, 150$ , $l = 1, 2, 3)$ , $ζ = 350$ , $κ = 10^{4}$ . In (C): $a_{i} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i} = 0.3$ , $d_{i}^{'} = 0.3$ , $w_{i} = 0.3$ , $k_{i} = 0.2$ , $ζ = 350$ , $κ = 10^{4}$ , $D_{i} = 0.03 \times N (1, 0.27)$ , $(i = 1, \dots, S_{C})$ .

Appendix 1—figure 13

Download asset Open asset

Appendix 1—figure 14

Download asset Open asset

A model of intraspecific interference illustrates the rank-abundance curve of a bird community ( $S_{C} ≫ S_{R}$ ).

(A) The rank-abundance curve. The solid dots represent the bird community data (marked with ‘Exp’) collected longitudinally within the same Amazonian region in 1982 (blue) and 2018 (cyan) (Terborgh et al., 1990; Martínez et al., 2023). The hollow dots are the ODEs results constructed from timestamp $t = 1.0 \times 10^{5}$ in the time series (see (B)). (B) Time courses of the species abundances simulated with ODEs. In the model settings, $S_{R} = 3$ and $S_{C} = 250$ is the only parameter varying with the consumer species, which was randomly drawn from a Gaussian distribution. The Shannon entropies of the experimental data and simulation results for the bird community are $H_{Exp(ODEs)}^{bird(1982/2018)} = 5.67 / 6.63 (7.31)$ . In the K-S test, the p-values that the simulation results and the corresponding experimental data come from identical distributions are: $p_{ODEs}^{bird(1982)} = 0.28$ , $p_{ODEs}^{bird(2018)} = 0.46$ . With a significance threshold of 0.05, none of the p-values suggest there exists a statistically significant difference. (C) Time courses of the species abundances simulated with ODEs corresponding to Figure 3C (bird), and the simulation parameters are the same as Figure 3C (bird). The numerical results in (**A–C**) were simulated from Equations 1, 2 and 4. In (**A, B**): $a_{i l} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i l} = 0.5$ , $d_{i}^{'} = 0.6$ , $w_{i l} = 0.3$ , $k_{i l} = 0.2$ , $D_{i} = 0.032 \times N (1, 0.17)$ , $(i = 1, \dots, 250,$ $l = 1, 2, 3)$ , $κ_{1} = 5 \times 10^{4}$ , $κ_{2} = 3 \times 10^{4}$ , $κ_{3} = 10^{4}$ , $ζ_{1} = 100$ , $ζ_{2} = 70$ , $ζ_{3} = 40$ .

Appendix 1—figure 15

Download asset Open asset

Intraspecific interference results in an underlying negative feedback loop and thus promotes biodiversity.

(**A, B**) The fraction of consumer individuals engaged in pairwise encounter. Here, $S_{C} = 40$ and $S_{R} = 1$ . x_i represents $C_{i}^{(P)} \lor R^{(P)}$ , and y_i represents $C_{i}^{(P)} \lor C_{i}^{(P)}$ and $y_{i} / C_{i}$ stand for the fractions of consumer individuals within a chasing pair and an interference pair, respectively. The numerical results were calculated from Equation S55, while the analytical solutions (marked with superscript ‘(A)’) were calculated from Equation S59. The orange surface in (A) is an overlap of the red and yellow surfaces. (C) The formation of intraspecific interference results in a self-inhibiting negative feedback loop. In (**A, B**): $a_{i} = 0.0015, a_{i}^{'} = 0.0021, d_{i} = 0.1, d_{i}^{'} = 0.05, k_{i} = 5$ . In (B): $R = 2000$ .

Appendix 2

The classical proof of Competitive Exclusion Principle (CEP)

In the 1960s, MacArthur (MacArthur and Levins, 1964) and Levin (Levin, 1970) put forward the classical mathematical proof of CEP. We rephrase their idea in the simple case of $S_{C} = 2$ and $S_{R} = 1$ , that is two consumer species $C_{1}$ and $C_{2}$ competing for one resource species $R$ . In practice, this proof can be generalized into higher dimensions with several consumer and resource species. The population dynamics of the system can be described as follows:

{\begin{cases} {\dot{C}}_{i} = C_{i} (f_{i} (R) - D_{i}), i = 1, 2; \\ \dot{R} = g (R, C_{1}, C_{2}) . \end{cases}

Here, $C_{i}$ and $R$ represent the population abundances of consumers and resources, respectively, while the functional forms of $f_{i} (R)$ and $g (R, C_{1}, C_{2})$ are unspecific. $D_{i}$ stands for the mortality rate of the species $C_{i}$ . If all consumer species can coexist at steady state, then $f_{i} (R) / D_{i} = 1$ $(i = 1, 2)$ . In a 2-D representation, this requires that three lines $f_{i} (R) / D_{i} = 1$ and $y = 1$ share a common point, which is commonly impossible unless the model parameters satisfy special constraint (sets of Lebesgue measure zero). In a 3-D representation, the two planes corresponding to $f_{i} (R) / D_{i} = 1$ $(i = 1, 2)$ are parallel, and hence do not share a common point (see Wang and Liu, 2020 for details).

Appendix 3

Comparison of the functional response with Beddington-DeAngelis (B-D) model

A B-D model

In 1975, Beddington proposed a mathematical model (Beddington, 1975) to describe the influence of predator interference on the functional response with hand-waving derivations. In the same year, DeAngelis and his colleagues considered a related question and put forward a similar model (DeAngelis et al., 1975). Essentially, both models are phenomenological, and they were called B-D model in the subsequent studies. In practice, the B-D model can be extended into scenarios involving different types of pairwise encounters with Beddington’s modelling method. In this section, we systematically compare the functional response in B-D model with that of our mechanistic model in all the relevant scenarios.

Recalling Beddington’s analysis, the model (Beddington, 1975) consists of one consumer species $C$ and one resource species $R$ ( $S_{C} = 1$ , $S_{R} = 1$ ). In a well-mixed system, an individual consumer meets a resource with rate $a$ , while encounters another consumer with rate $a^{'}$ . There are two other phenomenological parameters in this model, namely, the handling time $t_{h}$ and the wasting time $t_{w}$ . Both can be determined by specifying the scenario and using statistical physics modeling analysis. In fact, Beddington analyzed the searching efficiency $Ξ_{B-D}$ rather than the functional response $F_{B-D}$ , yet both can be reciprocally derived with $Ξ_{B-D} \equiv F_{B-D} / R$ . Here $R$ stands for the population abundance of the resources, and the specific form of $Ξ_{B-D}$ is (Beddington, 1975):

Ξ_{B-D} (R, C) = \frac{a}{1 + a t_{h} R + a^{'} t_{w} C^{'}},

where $C^{'} = C - 1$ , and $C$ stands for the population abundance of the consumes. Generally, $C ≫ 1$ , and thus $C^{'} \approx C$ .

B Scenario involving only chasing pairs

Here, we consider the scenario involving only chasing pair for the simple case with one consumer species $C$ and one resource species $R$ ( $S_{C} = 1, S_{R} = 1$ ). When an individual consumer is chasing a resource, they form a chasing pair:

C^{(F)} + R^{(F)} ⇌_{d}^{a} C^{(P)} \lor R^{(P)} \overset{k}{\to} C^{(F)} (+),

where the superscript ‘(F)’ stands for populations that are freely wandering, and ‘(+)’ signifies gaining biomass (we count $C^{(F)} (+)$ as $C^{(F)}$ ). $C^{(P)} \lor R^{(P)}$ represents chasing pair (where ‘(P)’ signifies pair), denoted as $x$ . $a$ , $d$ and $k$ stand for encounter rate, escape rate and capture rate, respectively. Hence, the total number of consumers and resources are $C \equiv C^{(F)} + x$ and $R \equiv R^{(F)} + x$ . Then, the population dynamics of the system follows:

{\begin{cases} \dot{x} = a C^{(F)} R^{(F)} - (k + d) x, \\ \dot{C} = w k x - D C, \\ \dot{R} = g (R, x, C) . \end{cases}

Here, the functional form of $g (R, x, C)$ is unspecific, while $D$ and $w$ represent the mortality rate of the consumer species and biomass conversion ratio (Wang and Liu, 2020), respectively. Since consumption process is generically much faster than the birth/death process, in deriving the functional response, the consumption process is supposed to be in fast equilibrium (i.e. $\dot{x} = 0$ ). Then, we can solve for $x$ with:

x^{2} - (R + C + K) x + R C = 0,

where $K = \frac{k + d}{a}$ , and then,

x = \frac{2 R C}{(R + C + K)} \frac{1}{(1 + \sqrt{1 - \frac{4 R C}{(R + C + K)^{2}}})} .

By definition, the functional response and searching efficiency are:

F_{CP} (R, C) = \frac{k x}{C},

Ξ_{CP} (R, C) = \frac{k x}{R C} .

Hence, we obtain the functional response and searching efficiency in this chasing-pair scenario:

F_{CP} (R, C)_{(1)} = k \frac{(R + C + K)}{2 C} (1 - \sqrt{1 - \frac{4 R C}{(R + C + K)^{2}}}),

Ξ_{CP} (R, C)_{(1)} = k \frac{(R + C + K)}{2 R C} (1 - \sqrt{1 - \frac{4 R C}{(R + C + K)^{2}}}) .

Since $\frac{4 R C}{(R + C + K)^{2}} < 4 \frac{C}{R} ≪ 1$ , using first-order approximations in Equation S7a, Equation S7b, we obtain $\sqrt{1 - \frac{4 R C}{(R + C + K)^{2}}} \approx 1 - \frac{2 R C}{(R + C + K)^{2}}$ . Then the functional response and searching efficiency are:

F_{CP} (R, C)_{(2)} = k \frac{R}{R + C + K},

Ξ_{CP} (R, C)_{(2)} = \frac{k}{R + C + K} .

Evidently, there is no predator interference within the chasing-pair scenario, yet the functional response form is identical to the B-D model involving intraspecific interference (see Equation S2). Meanwhile, using first-order approximations in the denominator of Equation S5, we have $x \approx \frac{R C}{(R + C + K) - \frac{R C}{(R + C + K)}}$ . Hence,

F_{CP} (R, C)_{(3)} = k \frac{R}{(R + C + K) - \frac{R C}{(R + C + K)}},

Ξ_{CP} (R, C)_{(3)} = \frac{k}{(R + C + K) - \frac{R C}{(R + C + K)}} .

In the case that $R ≫ C$ , then $R ≫ C > x = R - R^{(F)}$ . By applying $R \approx R^{(F)}$ in Equation S3, we obtain $x \approx \frac{R C}{R + K}$ . Then,

F_{CP} (R, C)_{(4)} = k \frac{R}{R + K},

Ξ_{CP} (R, C)_{(4)} = \frac{k}{R + K} .

To compare these functional responses with that of the B-D model, we determine the parameters t_h and t_w in the B-D model by calculating their ensemble average values in a stochastic framework. Using the properties of waiting time distribution in the Poisson process, we obtain $⟨ t_{h} ⟩ = \frac{1}{k}$ and $⟨ t_{w} ⟩ = \frac{1}{d^{'}}$ (in the chasing-pair scenario, $a^{'} = 0$ ). By substituting these calculations into Equation S2, we have:

Ξ_{CP}^{B-D} (R, C) = \frac{a}{1 + R a / k} = \frac{k}{k / a + R},

F_{CP}^{B-D} (R, C) = \frac{k R}{k / a + R} .

In the special case with $d = 0$ and $R ≫ C$ , the B-D model is consistent with our mechanistic model: $Ξ_{B-D} (R, C) = Ξ_{CP} (R, C)_{(4)}$ . Outside the special region, however, the discrepancy can be considerably large (see Appendix 1—figure 2A–C for the comparison).

C Scenario involving chasing pairs and intraspecific interference

Here, we consider the scenario with additional involvement of intraspecific interference in the simple case of $S_{C} = 1$ and $S_{R} = 1$ :

\begin{array}{ll} C^{(F)} + R^{(F)} ⇌_{d}^{a} C^{(P)} \lor R^{(P)} \overset{k}{\to} C^{(F)} (+), \\ C^{(F)} + C^{(F)} ⇌_{d^{'}}^{a^{'}} C^{(P)} \lor C^{(P)} . \end{array}

Here, $C^{(P)} \lor C^{(P)}$ stands for the intraspecific predator interference pair, denoted as $y$ ; $a^{'}$ and $d^{'}$ represent the encounter rate and separation rate of the interference pair, respectively. Then, the total population of consumers and resources are $C \equiv C^{(F)} + x + 2 y$ and $R \equiv R^{(F)} + x$ . Hence the population dynamics of the consumers and resources can be described as follows:

{\begin{cases} \dot{x} = a C^{(F)} R^{(F)} - (k + d) x, \\ \dot{y} = a^{'} [C^{(F)}]^{2} - d^{'} y, \\ \dot{C} = w k x - D C, \\ \dot{R} = g (R, x, C) . \end{cases}

The consumption process and interference process are supposed to be in fast equilibrium (i.e., $\dot{x} = 0, \dot{y} = 0$ ), then we can solve for $x$ with:

x^{3} + ϕ_{2} x^{2} + ϕ_{1} x + ϕ_{0} = 0,

where $ϕ_{0} = - C R^{2}$ , $ϕ_{1} = 2 C R + K R + R^{2}$ , $ϕ_{2} = 2 β K^{2} - K - C - 2 R$ with $β = a^{'} / d^{'}$ . The discriminant of Equation S13 (denoted as $Λ$ ) is:

Λ = - 4 ψ^{3} - 27 φ^{2},

with $ψ = ϕ_{1} - (ϕ_{2})^{2} / 3$ and $φ = ϕ_{0} - ϕ_{1} ϕ_{2} / 3 + 2 (ϕ_{2})^{3} / 27$ . When $Λ < 0$ , there are one real solution $x_{(1)}$ and two complex solutions $x_{(2)}, x_{(3)}$ , which are:

x_{(1)} = θ_{1} + θ_{2} - ϕ_{2} / 3, x_{(2)} = ω θ_{1} + ω^{2} θ_{2} - ϕ_{2} / 3, x_{(3)} = ω^{2} θ_{1} + ω θ_{2} - ϕ_{2} / 3,

where $ω = - 1 / 2 + i \sqrt{3} / 2$ ( $i$ stands for the imaginary unit), $θ_{1} = (- φ / 2 + \sqrt{- Λ / 108})^{1 / 3}$ , and $θ_{2} = (- φ / 2 - \sqrt{- Λ / 108})^{1 / 3}$ . On the other hand, when $Λ > 0$ , there are three real solutions $x_{(1)}, x_{(2)}$ , and $x_{(3)}$ , which are:

x_{(1)} = ψ^{'} \cos φ^{'} - ϕ_{2} / 3, x_{(2)} = ψ^{'} \cos (φ^{'} + \frac{2 π}{3}) - ϕ_{2} / 3, x_{(3)} = ψ^{'} \cos (φ^{'} + \frac{4 π}{3}) - ϕ_{2} / 3,

where $ψ^{'} = (- 4 ψ / 3)^{1 / 2}$ and $φ^{'} = \arccos (- (- ψ / 3)^{- 3 / 2} φ / 2) / 3$ . Note that $x \in [0, min (R, C)]$ , then we obtain the exact feasible solution of $x$ (denoted as x_ext), and hence the functional response and searching efficiency are:

F_{intra} (R, C)_{(1)} = \frac{k x_{e x t}}{C},

Ξ_{intra} (R, C)_{(1)} = \frac{k x_{e x t}}{R C} .

In the case of $R ≫ C$ , then $R - R^{(F)} = x < C ≪ R$ , and thus $R^{(F)} \approx R$ . Still, the consumption process is supposed to be in fast equilibrium (i.e. $\dot{x} = 0, \dot{y} = 0$ ), and then we obtain:

x \approx \frac{R C}{\sqrt{[\frac{1}{2} (K + R)]^{2} + 2 C β K^{2}} + \frac{1}{2} (K + R)} .

Consequently,

F_{intra} (R, C)_{(2)} = k \frac{R}{\sqrt{[\frac{1}{2} (K + R)]^{2} + 2 C β K^{2}} + \frac{1}{2} (K + R)},

Ξ_{intra} (R, C)_{(2)} = k \frac{1}{\sqrt{[\frac{1}{2} (K + R)]^{2} + 2 C β K^{2}} + \frac{1}{2} (K + R)} .

When $β ≪ \frac{1}{8 C}$ or $8 β C / (1 + R / K)^{2} ≪ 1$ , using first-order approximations in the denominator of Equation S18, we have:

x \approx \frac{R C}{(K + R) + \frac{2 K}{(1 + R / K)} β C},

and then,

F_{intra} (R, C)_{(3)} = k \frac{R}{(K + R) + \frac{2 K}{(1 + R / K)} β C},

Ξ_{intra} (R, C)_{(3)} = k \frac{1}{(K + R) + \frac{2 K}{(1 + R / K)} β C} .

In the case that $8 β C / (1 + R / K)^{2} ≫ 1$ , using first-order approximations in Equation S18, we obtain:

x \approx \frac{R C}{K \sqrt{2 C β} + \frac{(K + R)^{2}}{8 K \sqrt{2 C β}} + \frac{1}{2} (K + R)},

and thus,

F_{intra} (R, C)_{(4)} = k \frac{R}{K \sqrt{2 C β} + \frac{(K + R)^{2}}{8 K \sqrt{2 C β}} + \frac{1}{2} (K + R)},

Ξ_{intra} (R, C)_{(4)} = k \frac{1}{K \sqrt{2 C β} + \frac{(K + R)^{2}}{8 K \sqrt{2 C β}} + \frac{1}{2} (K + R)} .

Meanwhile, the B-D model only fits to the cases with $d = 0$ . By calculating the average values of $t_{h}$ and $t_{w}$ in the stochastic framework, we have $⟨ t_{h} ⟩ = \frac{1}{k}, ⟨ t_{w} ⟩ = \frac{1}{d^{'}}$ . Thus, we obtain the searching efficiency and functional response in the B-D model:

Ξ_{intra}^{B-D} (R, C) = \frac{a}{1 + \frac{a}{k} R + \frac{a^{'}}{d^{'}} C} = \frac{a}{1 + R / K ∣_{d = 0} + β C},

F_{intra}^{B-D} (R, C) = \frac{a R}{1 + R / K ∣_{d = 0} + β C} .

Overall, the searching efficiency (and the functional response) of the B-D model is quite different from either the rigorous form $Ξ_{intra} (R, C)_{(1)}$ , the quasi rigorous form $Ξ_{intra} (R, C)_{(2)}$ , or the more simplified forms $Ξ_{intra} (R, C)_{(3)}$ and $Ξ_{intra} (R, C)_{(4)}$ (Appendix 1—figure 2D–F). Still, there is a region where the discrepancies can be small, namely $d \approx 0$ and $R ≫ C$ (Appendix 1—figure 2D–F). Intuitively, when $β ≪ \frac{1}{8 C}$ and $d = 0$ , then $Ξ_{intra} (R, C)_{(3)} = \frac{a}{1 + \frac{a}{k} R + \frac{2}{(1 + R / K)} β C}$ . Consequently, if $R / K = x / C^{(F)} < 1$ , then $\frac{2}{(1 + R / K)} \in [1, 2]$ . In this case, the difference between $Ξ_{intra}^{B-D} (R, C)$ and $Ξ_{intra} (R, C)_{(3)}$ is small.

In fact, the above analysis also applies to cases with more than one types of consumer species (i.e., for cases with $S_{C} > 1$ ).

D Scenario involving chasing pairs and interspecific interference

Next, we consider the scenario involving chasing pairs and interspecific interference in the case of $S_{C} = 2$ and $S_{R} = 1$ :

\begin{array}{ll} C_{i}^{(F)} + R^{(F)} ⇌_{d_{i}}^{a_{i}} C_{i}^{(P)} \lor R^{(P)} \overset{k_{i}}{\to} C_{i}^{(F)} (+), i = 1, 2; \\ C_{1}^{(F)} + C_{2}^{(F)} ⇌_{d_{12}^{'}}^{a_{12}^{'}} C_{1}^{(P)} \lor C_{2}^{(P)} . \end{array}

Here $C_{1}^{(P)} \lor C_{2}^{(P)}$ stands for the interspecific interference pair, denoted as $z$ ; $a_{12}^{'}$ and $d_{12}^{'}$ represent the encounter rate and separation rate of the interference pair, respectively. Then, the total population of consumers and resources are $C_{i} \equiv C_{i}^{(F)} + x_{i} + z$ and $R \equiv R^{(F)} + x_{1} + x_{2}$ . The population dynamics of the consumers and resources follows:

\begin{array}{ll} {\begin{cases} {\dot{x}}_{i} = a_{i} C_{i}^{(F)} R^{(F)} - (k_{i} + d_{i}) x_{i}, i = 1, 2; \\ \dot{z} = a_{12}^{'} C_{1}^{(F)} C_{2}^{(F)} - d_{12}^{'} z, \\ {\dot{C}}_{i} = w_{i} k_{i} x_{i} - D_{i} C_{i}, \\ \dot{R} = g (R, x_{1}, x_{2}, C_{1}, C_{2}) . \end{cases} \end{array}

where the functional form of $g (R, x_{1}, x_{2}, C_{1}, C_{2})$ is unspecific, while $D_{i}$ and $w_{i}$ represents the mortality rates of the two consumers species and biomass conversion ratios. Still, the consumption/interference process is supposed to be in fast equilibrium, that is ${\dot{x}}_{i} = 0, \dot{z} = 0$ . In the case that $R ≫ C_{1} + C_{2} > x_{1} + x_{2}$ , by applying $R^{(F)} \approx R$ , we obtain:

x_{1} \approx \frac{2 C_{1} (R / K_{2} + 1) R / K_{1}}{\sqrt{[γ (C_{2} - C_{1}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)]^{2} + 4 γ C_{1} (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)} + γ (C_{2} - C_{1}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)},

x_{2} \approx \frac{2 C_{2} (R / K_{1} + 1) R / K_{2}}{\sqrt{[γ (C_{1} - C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)]^{2} + 4 γ C_{2} (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)} + γ (C_{1} - C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)} .

Then, the searching efficiencies and functional responses are:

Ξ_{1}^{inter} (R, C_{1}, C_{2})_{(1)} = \frac{2 k_{2} (R / K_{1} + 1) R / K_{2}}{γ (C_{1} - C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1) + \sqrt{[γ (C_{1} - C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)]^{2} + 4 γ C_{2} (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)}},

Ξ_{2}^{inter} (R, C_{1}, C_{2})_{(1)} = \frac{2 k_{2} (R / K_{1} + 1) / K_{2}}{γ (C_{1} - C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1) + \sqrt{[γ (C_{1} - C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)]^{2} + 4 γ C_{2} (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)}},

F_{1}^{inter} (R, C_{1}, C_{2})_{(1)} = \frac{2 k_{1} (R / K_{2} + 1) R / K_{1}}{γ (C_{2} - C_{1}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1) + \sqrt{[γ (C_{2} - C_{1}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)]^{2} + 4 γ C_{1} (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)}},

F_{2}^{inter} (R, C_{1}, C_{2})_{(1)} = \frac{2 k_{2} (R / K_{1} + 1) R / K_{2}}{γ (C_{1} - C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1) + \sqrt{[γ (C_{1} - C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)]^{2} + 4 γ C_{2} (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)}} .

Since $\frac{4 γ^{2} C_{1} C_{2}}{[γ C_{1} + γ C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)]^{2}} < 1$ , by applying first-order approximation to the denominator of Equation S26b, we obtain:

x_{1} \approx \frac{C_{1} R}{(R + K_{1}) + \frac{γ K_{1} K_{2} C_{2}}{(R + K_{2})} - \frac{γ^{2} K_{1} K_{2} C_{1} C_{2}}{[γ (C_{1} + C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)] (R + K_{2})}},

x_{2} \approx \frac{C_{2} R}{(R + K_{2}) + \frac{γ K_{1} K_{2} C_{1}}{(R + K_{1})} - \frac{γ^{2} K_{1} K_{2} C_{1} C_{2}}{[γ (C_{1} + C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)] (R + K_{1})}},

and the searching efficiencies and functional responses are:

Ξ_{1}^{inter} (R, C_{1}, C_{2})_{(2)} = \frac{k_{1}}{(R + K_{1}) + \frac{γ K_{1} K_{2} C_{2}}{(R + K_{2})} - \frac{γ^{2} K_{1} K_{2} C_{1} C_{2}}{[γ (C_{1} + C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)] (R + K_{2})}},

Ξ_{2}^{inter} (R, C_{1}, C_{2})_{(2)} = \frac{k_{2}}{(R + K_{2}) + \frac{γ K_{1} K_{2} C_{1}}{(R + K_{1})} - \frac{γ^{2} K_{1} K_{2} C_{1} C_{2}}{[γ (C_{1} + C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)] (R + K_{1})}},

F_{1}^{inter} (R, C_{1}, C_{2})_{(2)} = \frac{k_{1} R}{(R + K_{1}) + \frac{γ K_{1} K_{2} C_{2}}{(R + K_{2})} - \frac{γ^{2} K_{1} K_{2} C_{1} C_{2}}{[γ (C_{1} + C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)] (R + K_{2})}},

F_{2}^{inter} (R, C_{1}, C_{2})_{(2)} = \frac{k_{2} R}{(R + K_{2}) + \frac{γ K_{1} K_{2} C_{1}}{(R + K_{1})} - \frac{γ^{2} K_{1} K_{2} C_{1} C_{2}}{[γ (C_{1} + C_{2}) + (\frac{R}{K_{1}} + 1) (\frac{R}{K_{2}} + 1)] (R + K_{1})}} .

Likewise, the B-D model only fits to cases with $d = 0$ . By calculating the average values in a stochastic framework, we obtain $⟨ t_{h}^{i} ⟩ = \frac{1}{k_{i}}$ , $⟨ t_{w}^{i} ⟩ = \frac{1}{d_{12}^{'}}$ ( $i = 1, 2$ ). Then, we obtain the searching efficiencies in the B-D model:

Ξ_{1}^{B-D (inter)} (R, C_{1}, C_{2}) = \frac{a_{1}}{1 + \frac{a_{1}}{k_{1}} R + \frac{a_{12}^{'}}{d_{12}^{'}} C_{2}} = \frac{a_{1}}{1 + R / K_{1} ∣_{d = 0} + γ C_{2}},

Ξ_{2}^{B-D (inter)} (R, C_{1}, C_{2}) = \frac{a_{2}}{1 + \frac{a_{2}}{k_{2}} R + \frac{a_{12}^{'}}{d_{12}^{'}} C_{1}} = \frac{a_{2}}{1 + R / K_{2} ∣_{d = 0} + γ C_{1}} .

Consequently, the functional responses in the B-D model are:

F_{1}^{B-D (inter)} (R, C_{1}, C_{2}) = \frac{a_{1} R}{1 + R / K_{1} ∣_{d = 0} + γ C_{2}},

F_{2}^{B-D (inter)} (R, C_{1}, C_{2}) = \frac{a_{2} R}{1 + R / K_{2} ∣_{d = 0} + γ C_{1}} .

Evidently, the searching efficiencies in the B-D model are overall different from either the quasi rigorous form $Ξ_{i} (R, C_{1}, C_{2})_{1}$ , or the simplified form $Ξ_{i} (R, C_{1}, C_{2})_{2}$ (Appendix 1—figure 2G–I). Still, the discrepancy can be small when $d \approx 0$ and $R ≫ C$ (Appendix 1—figure 2G–I). Intuitively, when $γ ≪ min (C_{1}^{- 1}, C_{2}^{- 1})$ , we have:

Ξ_{1}^{inter} (R, C_{1}, C_{2})_{(2)} \approx \frac{a_{1}}{(1 + \frac{a_{1}}{k_{1}} R) + \frac{γ C_{2}}{R / K_{2} + 1}},

Ξ_{2}^{inter} (R, C_{1}, C_{2})_{(2)} \approx \frac{a_{2}}{(1 + \frac{a_{2}}{k_{2}} R) + \frac{γ C_{1}}{R / K_{1} + 1}} .

Thus, if $R / K_{i} = x_{i} / C_{i}^{(F)} < 1$ ( $i = 1, 2$ ), then $\frac{1}{1 + R / K_{i}} \in [0.5, 1]$ . In this case, the difference between $Ξ_{i}^{B-D (inter)} (R, C_{1}, C_{2})$ and $Ξ_{i}^{inter} (R, C_{1}, C_{2})_{(2)}$ is small.

Appendix 4

Scenario involving chasing pairs and intraspecific interference

A Two consumers species competing for one resource species

We consider the scenario involving chasing pairs and intraspecific interference in the simple case of $S_{C} = 2$ and $S_{R} = 1$ :

\begin{array}{ll} C_{i}^{(F)} + R^{(F)} ⇌_{d_{i}}^{a_{i}} C_{i}^{(P)} \lor R^{(P)} \overset{k_{i}}{\to} C_{i}^{(F)} (+), \\ C_{i}^{(F)} + C_{i}^{(F)} ⇌_{d_{i}^{'}}^{a_{i}^{'}} C_{i}^{(P)} \lor C_{i}^{(P)}, i = 1, 2. \end{array}

Here, the variables and parameters are just extended from the case of $S_{C} = 1$ and $S_{R} = 1$ (see Appendix 3.C). The total number of consumers and resources are $C_{i} \equiv C_{i}^{(F)} + x_{i} + 2 y_{i}$ and $R \equiv R^{(F)} + \sum_{i = 1}^{2} x_{i}$ . Then, the population dynamics of the consumers and resources can be described as follows:

{\begin{cases} {\dot{x}}_{i} = a_{i} C_{i}^{(F)} R^{(F)} - (k_{i} + d_{i}) x_{i}, i = 1, 2; \\ {\dot{y}}_{i} = a_{i}^{'} [C_{i}^{(F)}]^{2} - d_{i}^{'} y_{i}, \\ {\dot{C}}_{i} = w_{i} k_{i} x_{i} - D_{i} C_{i}, \\ \dot{R} = g (R, x_{1}, x_{2}, C_{1}, C_{2}) . \end{cases}

The functional form of $g (R, x_{1}, x_{2}, C_{1}, C_{2})$ is unspecified. For simplicity, we limit our analysis to abiotic resources, while all results generically apply to biotic resources. Besides, we define $K_{i} \equiv (d_{i} + k_{i}) / a_{i}$ , $α_{i} \equiv D_{i} / (w_{i} k_{i})$ and $β_{i} \equiv a_{i}^{'} / d_{i}^{'}$ $(i = 1, 2)$ . At steady state, from ${\dot{x}}_{i} = 0, {\dot{y}}_{i} = 0$ , we have:

{\begin{cases} x_{i} = C_{i}^{(F)} R^{(F)} / K_{i}, i = 1, 2; \\ y_{i} = β_{i} [C_{i}^{(F)}]^{2} . \end{cases}

Note that $C_{i} \equiv C_{i}^{(F)} + x_{i} + 2 y_{i}$ , and $R \equiv R^{(F)} + \sum_{i = 1}^{2} x_{i}$ . Then,

R^{(F)} = R / (1 + C_{1}^{(F)} / K_{1} + C_{2}^{(F)} / K_{2}),

C_{i} = C_{i}^{(F)} + R^{(F)} C_{i}^{(F)} / K_{i} + 2 β_{i} [C_{i}^{(F)}]^{2}, i = 1, 2.

By substituting Equation S35a into Equation S35b, we have:

C_{2}^{(F)} = \frac{K_{2}}{K_{1}} [\frac{R C_{1}^{(F)}}{C_{1} - C_{1}^{(F)} - 2 β_{1} [C_{1}^{(F)}]^{2}} - K_{1} - C_{1}^{(F)}],

(C_{2} - C_{2}^{(F)} - 2 β_{2} [C_{2}^{(F)}]^{2}) (1 + C_{1}^{(F)} / K_{1} + C_{2}^{(F)} / K_{2}) = R C_{2}^{(F)} / K_{2} .

Then, we can present $C_{i}^{(F)}$ with $C_{1}$ , $C_{2}$ and $R$ ( $i = 1, 2$ ). By further combining with Equation S34, Equation S35a and Equation S36a, we express $R^{(F)}$ , $x_{i}$ and $y_{i}$ using $C_{1}$ , $C_{2}$ and $R$ . In particular, for $x_{i}$ , we have:

x_{i} = u_{i} (R, C_{1}, C_{2}), i = 1, 2.

If all species coexist, then the steady-state equations of $\dot{C_{i}} = 0$ and $\dot{R} = 0$ are:

{\begin{cases} Ω_{1} (R, C_{1}, C_{2}) - D_{1} = 0, \\ Ω_{2} (R, C_{1}, C_{2}) - D_{2} = 0, \\ G (R, C_{1}, C_{2}) = 0, \end{cases}

where $G (R, C_{1}, C_{2}) \equiv g (R, u_{1} (R, C_{1}, C_{2}), u_{2} (R, C_{1}, C_{2}), C_{1}, C_{2})$ , and $Ω_{i} (R, C_{1}, C_{2}) \equiv \frac{w_{i} k_{i}}{C_{i}} u_{i} (R, C_{1}, C_{2})$ . In practice, Equation S38 corresponds to three unparallel surfaces, which share a common point (Figure 1H and Appendix 1—figure 3G). Importantly, the fixed point can be stable, and hence two consumer species may coexist at constant population densities.

1 Stability analysis of the fixed-point solution

We use linear stability analysis to study the local stability of the fixed point. Specifically, for an arbitrary fixed point $E (x_{1}, x_{2}, y_{1}, y_{2}, C_{1}, C_{2}, R)$ , only when all the eigenvalues (defined as $λ_{i}, i = 1, \dots, 7$ ) of the Jacobian matrix at point $E$ own negative real parts would the point be locally stable.

To investigate whether there exists a non-zero measure parameter region for species coexistence, we set $D_{i}$ $(i = 1, 2)$ to be the only parameter that varies with species $C_{1}$ and $C_{2}$ , and then $Δ \equiv (D_{1} - D_{2}) / D_{2}$ reflects the completive difference between the two consumer species. As shown in Appendix 1—figure 4B, the region below the blue surface and above the red surface corresponds to stable coexistence. Thus, there exists a non-zero measure parameter region to promote species coexistence, which breaks CEP.

2 Analytical solutions of the species abundances at steady state

At steady state, since ${\dot{x}}_{i} = {\dot{y}}_{i} = {\dot{C}}_{i} = 0$ $(i = 1, 2)$ , then,

{\begin{cases} x_{i} = α_{i} C_{i}, \\ C_{i}^{(F)} = K_{i} α_{i} C_{i} / R^{(F)}, \\ y_{i} = β_{i} (K_{i} α_{i} C_{i})^{2} [R^{(F)}]^{- 2} . \end{cases}

Meanwhile, $C_{i} = C_{i}^{(F)} + x_{i} + 2 y_{i}$ , and $C_{i}, R > 0$ . Then, we have:

C_{i} = \frac{(1 - α_{i}) [R^{(F)}]^{2} - K_{i} α_{i} R^{(F)}}{2 β_{i} (K_{i} α_{i})^{2}} .

If the resource species owns a much larger population abundance than the consumers (i.e. $R ≫ C_{1} + C_{2}$ ), then $R ≫ x_{1} + x_{2}$ , and $R^{(F)} \approx R$ . Thus,

C_{i} = \frac{(1 - α_{i}) R^{2} - K_{i} α_{i} R}{2 β_{i} (K_{i} α_{i})^{2}} .

By further assuming that the population dynamics of the resources follow identical construction rule as the MacArthur’s consumer-resource model (MacArthur, 1970), we have:

g (R, x_{1}, x_{2}, C_{1}, C_{2}) = ζ (1 - R / κ) - (k_{1} x_{1} + k_{2} x_{2}),

Since $\dot{R} = 0$ , then

R = \frac{- o_{1} + \sqrt{o_{1}^{2} + 4 o_{2} ζ}}{2 o_{2}},

where $o_{1} \equiv \frac{ζ}{κ} - \frac{k_{1}}{2 β_{1} K_{1}} - \frac{k_{2}}{2 β_{2} K_{2}}$ and $o_{2} \equiv \frac{k_{1} (1 - α_{1})}{2 β_{1} α_{1} (K_{1})^{2}} + \frac{k_{2} (1 - α_{2})}{2 β_{2} α_{2} (K_{2})^{2}}$

Equations. S41, S43 are the analytical solutions of species abundances at steady state when $R ≫ C_{1} + C_{2}$ . As shown in Figure 1E, the analytical solutions agree well with the numerical results (the exact solutions). To conduct a systematic comparison for different model parameters, we assign $D_{i}$ to be the only parameter varying with species $C_{1}$ and $C_{2}$ ( $D_{1} > D_{2}$ ), and define $Δ \equiv (D_{1} - D_{2}) / D_{2}$ as the competitive difference between the two consumer species. The comparison between analytical solutions and numerical results is shown in Appendix 1—figure 3H. Clearly, they are close to each other, exhibiting very good consistency.

Furthermore, we test if the parameter region for species coexistence is predictable using the analytical solutions. Since $D_{i}$ is the only parameter that varies with the two-consumer species, the supremum of the competitive difference tolerated for species coexistence (defined as $\hat{Δ}$ ) corresponds to the steady-state solutions that satisfy $R, C_{2} > 0$ and $C_{1} = 0^{+}$ , where $0^{+}$ stands for the infinitesimal positive number. To calculate the analytical solutions at the upper surface of the coexistence region, where $Δ = \hat{Δ}$ and $C_{1} = 0^{+}$ , we further combine Equation S41 and then obtain (note that $R > 0$ ):

R = \frac{K_{1} α_{1}}{1 - α_{1}} .

Meanwhile, $α_{1} = α_{2} (Δ + 1)$ . Thus, for the upper surface of the coexistence region:

α_{1} = α_{2} (\hat{Δ} + 1) .

Combining Equations S43-S45, we have:

\hat{Δ} = \frac{1}{α_{2} (κ_{1} ϖ + 1)} - 1,

where $ϖ \equiv \frac{1}{2} (\frac{1}{κ} - \frac{k_{2}}{2 ζ β_{2} K_{2}}) + \frac{1}{2} \sqrt{(\frac{1}{κ} - \frac{k_{2}}{2 ζ β_{2} K_{2}})^{2} + 2 \frac{k_{2} (1 - α_{2})}{ζ β_{2} α_{2} (K_{2})^{2}}} .$ When $R ≫ C_{1} + C_{2}$ , the comparison of $\hat{Δ}$ obtained from analytical solutions with that from numerical results (the exact solutions) are shown in Appendix 1—figure 3I, which overall exhibits good consistency.

B $S_{C}$ consumers species competing for $S_{R}$ resources species

Here, we consider the scenario involving chasing pairs and intraspecific interference for the generic case with $S_{C}$ types of consumers and $S_{R}$ types of resources. Then, the population dynamics of the system can be described as follows:

{\begin{cases} {\dot{x}}_{i l} = a_{i l} C_{i}^{(F)} R_{l}^{(F)} - (k_{i l} + d_{i l}) x_{i l}, \\ {\dot{y}}_{i} = a_{i i}^{'} [C_{i}^{(F)}]^{2} - d_{i i}^{'} y_{i}, \\ {\dot{C}}_{i} = \sum_{l = 1}^{S_{R}} w_{i l} k_{i l} x_{i l} - D_{i} C_{i}, \\ {\dot{R}}_{l} = g_{l} ({R_{l}}, {x_{i}}, {C_{i}}), i = 1, \dots, S_{C}, l = 1, \dots, S_{R} . \end{cases}

Note that Equation S47 is identical with Equations 1-2, and we use the same variables and parameters as that in the main text. Then, the populations of the consumers and resources are $C_{i} = C_{i}^{(F)} + \sum_{l = 1}^{S_{R}} x_{i l} + 2 y_{i}$ and $R_{l} = R_{l}^{(F)} + \sum_{i = 1}^{S_{C}} x_{i l}$ . For convenience, we define $K_{i l} \equiv (d_{i l} + k_{i l}) / a_{i l}$ , $α_{i l} \equiv D_{i l} / (k_{i l} w_{i l})$ and $β_{i} \equiv a_{i i}^{'} / d_{i i}^{'}$ ( $i = 1, \dots, S_{C}$ , $l = 1, \dots, S_{R}$ ).

1 Analytical solutions of species abundances at steady state

At steady state, from ${\dot{x}}_{i l} = 0, {\dot{y}}_{i} = 0$ , and ${\dot{C}}_{i} = 0$ , we have:

{\begin{cases} x_{i l} = C_{i}^{(F)} R_{l}^{(F)} / K_{i l}, \\ y_{i} = β_{i} [C_{i}^{(F)}]^{2}, \\ C_{i} = \sum_{l = 1}^{S_{R}} x_{i l} / α_{i l} = \sum_{l = 1}^{S_{R}} C_{i}^{(F)} R_{l}^{(F)} / (K_{i l} α_{i l}) . \end{cases}

Meanwhile $C_{i} = C_{i}^{(F)} + \sum_{l = 1}^{S_{R}} x_{i l} + 2 y_{i}$ , and note that $C_{i} > 0$ , thus,

C_{i}^{(F)} = \frac{1}{2 β_{i}} [- 1 + \sum_{l = 1}^{S_{R}} (\frac{1}{α_{i l}} - 1) \frac{R_{l}^{(F)}}{K_{i l}}] .

Combined with Equation S49, and then,

C_{i} = \sum_{l = 1}^{S_{R}} \frac{R_{l}^{(F)}}{2 β_{i} α_{i l} K_{i l}} [- 1 + \sum_{l^{'} = 1}^{S_{R}} (\frac{1}{α_{i l^{'}}} - 1) \frac{R_{l^{'}}^{(F)}}{K_{i l^{'}}}] .

We further assume that the specific function of $g_{l} ({R_{l}}, {x_{i}}, {C_{i}})$ satisfies Equation 4, that is

g_{l} ({R_{l}}, {x_{i}}, {C_{i}}) = ζ_{l} (1 - R_{l} / κ_{l}) - \sum_{i = 1}^{S_{C}} k_{i l} x_{i l} .

By combining Equations S48, S49 and S51, we have:

ζ_{l} (1 - \frac{R_{l}}{κ_{l}}) = \sum_{i = 1}^{S_{C}} \frac{k_{i l}}{2 β_{i} K_{i l}} [- 1 + \sum_{l^{'} = 1}^{S_{R}} (\frac{1}{α_{i l^{'}}} - 1) \frac{R_{l^{'}}^{(F)}}{K_{i l^{'}}}] R_{l}^{(F)} .

If the population abundance of each resource species is much more than the total population of all consumers (i.e. $R_{l} ≫ \sum_{i = 1}^{S_{C}} C_{i} (l = 1, \dots, S_{R})$ ), then $R_{l} ≫ \sum_{i = 1}^{S_{C}} x_{i l}$ and $R_{l}^{(F)} \approx R_{l}$ . Thus,

(\frac{ζ_{l}}{κ_{l}} - \sum_{i = 1}^{S_{C}} \frac{k_{i l}}{2 β_{i} K_{i l}} + \sum_{l^{'} = 1}^{S_{R}} \sum_{i = 1}^{S_{C}} \frac{k_{i l}}{2 β_{i} K_{i l}} (\frac{1}{α_{i l^{'}}} - 1) \frac{R_{l^{'}}}{K_{i l^{'}}}) R_{l} = ζ_{l},

with $l = 1, \dots, S_{R}$ . Equation S53 is a set of second-order algebraic differential equations, which is clearly solvable.

In fact, when $S_{R} = 1, S_{C} \geq 1$ , and $R_{l} ≫ \sum_{i = 1}^{S_{C}} C_{i}$ $(l = 1)$ , we can explicitly present the analytical solution of the steady-state species abundances. To simplify the notations, we omit the ‘ $l$ ’ in the sub-/super-scripts since $S_{R} = 1$ . Then, we have:

{\begin{cases} R = \frac{- ι_{1} + \sqrt{ι_{1}^{2} + 4 ι_{2} ζ}}{2 ι_{2}}, \\ C_{i} = \frac{1}{2 β_{i} α_{i} K_{i}} [(\frac{1}{α_{i}} - 1) \frac{R}{K_{i}} - 1] R, i = 1, \dots, S_{C} . \end{cases}

Here $ι_{1} \equiv \frac{ζ}{κ} - \sum_{i = 1}^{S_{C}} \frac{k_{i}}{2 β_{i} K_{i}}$ and $ι_{2} \equiv \sum_{i = 1}^{S_{C}} \frac{k_{i} (1 - α_{i})}{2 β_{i} α_{i} (K_{i})^{2}}$

C Intuitive understanding: an underlying negative feedback loop

Intuitively, how can intraspecific predator interference promote biodiversity? Here, we solve this question by considering the case that $S_{C}$ types of consumers compete for one resource species. The population dynamics of the system are described in Equations S47 and S51 with $S_{R} = 1$ . To simplify the notations, we omit the ‘ $l$ ’ in the subscript since $S_{R} = 1$ . The consumption process and interference process are supposed to be in fast equilibrium (i.e. $\dot{x_{i}} = 0, \dot{y_{i}} = 0$ ). Then, we have a set of equations to solve for $x_{i}$ and $y_{i}$ given the population size of each species:

{\begin{cases} x_{i} = C_{i}^{(F)} R^{(F)} / K_{i}, \\ y_{i} = β_{i} [C_{i}^{(F)}]^{2}, \\ C_{i} = C_{i}^{(F)} + x_{i} + 2 y_{i}, \\ R = R^{(F)} + \sum_{i = 1}^{S_{C}} x_{i} . \end{cases}

In the first three sub-equations of Equation S55, by getting rids of $C_{i}^{(F)}$ , we have,

{\begin{cases} 2 β_{i} [\frac{x_{i}}{R^{(F)}}]^{2} + x_{i} + \frac{K_{i} x_{i}}{R^{(F)}} - C_{i} = 0, \\ y_{i} = \frac{1}{2} [C_{i} - x_{i} - \frac{x_{i}}{R^{(F)}}] . \end{cases}

Then, by regarding $R^{(F)}$ as a temporary parameter, we solve for $x_{i}$ and $y_{i}$ :

{\begin{cases} x_{i} = \frac{2 R^{(F)} C_{i}}{\sqrt{(R^{(F)} + K_{i})^{2} + 8 β_{i} K_{i}^{2} C_{i}} + R^{(F)} + K_{i}}, \\ y_{i} = \frac{1}{2} [C_{i} - (1 + \frac{K_{i}}{R^{(F)}}) x_{i}] . \end{cases}

If the total population size of the resources is much larger than that of consumers (i.e. $R ≫ \sum_{i = 1}^{S_{C}} C_{i}$ ), then $R ≫ \sum_{i = 1}^{S_{C}} x_{i}$ and $R^{(F)} \approx R$ , and thus we get the analytical expressions of $x_{i}$ and $y_{i}$ :

{\begin{cases} x_{i} \approx \frac{2 R C_{i}}{\sqrt{(R + K_{i})^{2} + 8 β_{i} K_{i}^{2} C_{i}} + R + K_{i}}, \\ y_{i} \approx \frac{C_{i}}{2} [1 - \frac{2 (R + K_{i})}{\sqrt{(R + K_{i})^{2} + 8 C_{i} β_{i} K_{i}^{2}} + R + K_{i}}] . \end{cases}

Note that the fraction of $C_{i}$ individuals engaged in chasing pairs is $x_{i} / C_{i}$ , while that for individuals trapped in intraspecific interference pairs is $y_{i} / C_{i}$ . With Equation S58, it is straightforward to obtain these fractions:

{\begin{cases} x_{i} / C_{i} \approx \frac{2 R}{\sqrt{(R + K_{i})^{2} + 8 β_{i} K_{i}^{2} C_{i}} + R + K_{i}}, \\ y_{i} / C_{i} \approx \frac{1}{2} (1 - \frac{2 (R + K_{i})}{\sqrt{(R + K_{i})^{2} + 8 β_{i} K_{i}^{2} C_{i}} + R + K_{i}}) . \end{cases}

where both $x_{i} / C_{i}$ and $y_{i} / C_{i}$ are bivariate functions of $R$ and $C_{i}$ . From Equation S59, it is clear that for a given population size of the resource species, $y_{i} / C_{i}$ is a monotonously increasing function of $C_{i}$ , while $x_{i} / C_{i}$ is a monotonously decreasing function of $C_{i}$ . In Appendix 1—figure 15A, B, we see that the analytical results are highly consistence with the exact numerical solutions. By definition, the functional response of $C_{i}$ species is $F \equiv k_{i} x_{i} / C_{i}$ , and thus,

F (R, C_{i}) \approx \frac{2 R}{\sqrt{(R + K_{i})^{2} + 8 β_{i} K_{i}^{2} C_{i}} + R + K_{i}},

Evidently, the function response of $C_{i}$ species is negatively correlated with the population size of itself, which effectively constitutes a self-inhibiting negative feedback loop (Appendix 1—figure 15C).

Then, we have a simple intuitive understanding of species coexistence through the mechanism of intraspecific interference. In an ecological community, consumer species that of higher/lower competitiveness tend to increase/decrease their population size in the competition process. Without intraspecific interference, the increasing/decreasing trend would continue until the system obeys CEP. In the scenario involving intraspecific interference, however, for species of higher competitiveness (e.g. $C_{i}$ ), with the increase of $C_{i}$ ’s population size, a larger portion of $C_{i}$ individuals are then engaged in intraspecific interference pair which are temporarily absent from hunting (Appendix 1—figure 15A, B). Consequently, the functional response of $C_{i}$ drops, which prevents further increase of $C_{i}$ ’s population size, results in an overall balance among the consumer species, and thus promotes species coexistence.

Appendix 5

Scenario involving chasing pairs and interspecific interference

Here, we consider the scenario involving chasing pairs and interspecific interference in the case of $S_{C} = 2$ and $S_{R} = 1$ , with all settings follow that depicted in Appendix 3. D. Then, $C_{i} \equiv C_{i}^{(F)} + x_{i} + z, R \equiv R^{(F)} + x_{1} + x_{2}$ , and the population dynamics follows (identical with Equation S25):

{\begin{cases} {\dot{x}}_{i} = a_{i} C_{i}^{(F)} R^{(F)} - (k_{i} + d_{i}) x_{i}, i = 1, 2; \\ \dot{z} = a_{12}^{'} C_{1}^{(F)} C_{2}^{(F)} - d_{12}^{'} z, \\ {\dot{C}}_{i} = w_{i} k_{i} x_{i} - D_{i} C_{i}, \\ \dot{R} = g (R, x_{1}, x_{2}, C_{1}, C_{2}) . \end{cases}

Here, the functional form of $g (R, x_{1}, x_{2}, C_{1}, C_{2})$ is unspecified. For convenience, we define $K_{i} \equiv (d_{i} + k_{i}) / a_{i}$ , $α_{i} \equiv D_{i} / (w_{i} k_{i})$ $(i = 1, 2)$ and $γ \equiv a_{12}^{'} / d_{12}^{'}$ . At steady state, from ${\dot{x}}_{i} = 0$ $(i = 1, 2)$ and $\dot{z} = 0$ , we have:

{\begin{cases} x_{i} = C_{i}^{(F)} R^{(F)} / K_{i}, i = 1, 2; \\ z = γ C_{1}^{(F)} C_{2}^{(F)} . \end{cases}

Note that $C_{i} \equiv C_{i}^{(F)} + x_{i} + z$ and $R \equiv R^{(F)} + x_{1} + x_{2}$ , then,

{\begin{cases} C_{1} = C_{1}^{(F)} + R^{(F)} C_{1}^{(F)} / K_{1} + γ C_{1}^{(F)} C_{2}^{(F)}, \\ C_{2} = C_{2}^{(F)} + R^{(F)} C_{2}^{(F)} / K_{2} + γ C_{1}^{(F)} C_{2}^{(F)}, \\ R = R^{(F)} (1 + C_{1}^{(F)} / K_{1} + C_{2}^{(F)} / K_{2}) . \end{cases}

Then, we can express $C_{1}^{(F)}, C_{2}^{(F)}$ and $R^{(F)}$ with $C_{1}, C_{2}$ and $R$ . Combined with Equation S62, $x_{i}$ and $z$ can also be expressed using $C_{1}, C_{2}$ and $R$ . In particular, for $x_{i}$ , we have:

x_{i} = u_{i}^{'} (R, C_{1}, C_{2}), i = 1, 2.

If all species coexist, by defining $Ω_{i}^{'} (R, C_{1}, C_{2}) \equiv \frac{w_{i} k_{i}}{C_{i}} u_{i}^{'} (R, C_{1}, C_{2})$ , then, the steady-state equations of $\dot{C_{i}} = 0$ $(i = 1, 2)$ and $\dot{R} = 0$ are:

{\begin{cases} Ω_{1}^{'} (R, C_{1}, C_{2}) - D_{1} = 0, \\ Ω_{2}^{'} (R, C_{1}, C_{2}) - D_{2} = 0, \\ G^{'} (R, C_{1}, C_{2}) = 0, \end{cases}

where $G^{'} (R, C_{1}, C_{2}) \equiv g (R, u_{1}^{'} (R, C_{1}, C_{2}), u_{2}^{'} (R, C_{1}, C_{2}), C_{1}, C_{2})$ .

Here, Equation S65 corresponds to three unparallel surfaces and share a common point (Figure 1G and Appendix 1—figure 3A). However, all the fixed points are unstable (Appendix 1—figure 3F), and hence the consumer species cannot stably coexist at steady state (Figure 1D).

A Analytical results of the fixed-point solution

We proceed to investigate the unstable fixed point where $R, C_{1}, C_{2} > 0$ . From ${\dot{x}}_{i} = 0$ , $\dot{z} = 0$ , ${\dot{C}}_{i} = 0$ , and note that $C_{i} \equiv C_{i}^{(F)} + x_{i} + z$ , we have:

{\begin{cases} C_{i} = K_{i} α_{i} C_{i} (R^{(F)})^{- 1} + α_{i} C_{i} + z, i = 1, 2; \\ z = γ K_{1} α_{1} K_{2} α_{2} (R^{(F)})^{- 2} C_{1} C_{2} . \end{cases}

Since $C_{i} > 0$ , then:

{\begin{cases} C_{1} = \frac{(1 - α_{2}) [R^{(F)}]^{2} - K_{2} α_{2} R^{(F)}}{γ K_{1} α_{1} K_{2} α_{2}}, \\ C_{2} = \frac{(1 - α_{1}) [R^{(F)}]^{2} - K_{1} α_{1} R^{(F)}}{γ K_{1} α_{1} K_{2} α_{2}} . \end{cases}

If $R ≫ C_{1} + C_{2}$ , then $R ≫ x_{1} + x_{2}$ and $R^{(F)} \approx R$ , we have:

{\begin{cases} C_{1} = \frac{(1 - α_{2}) R^{2} - K_{2} α_{2} R}{γ K_{1} α_{1} K_{2} α_{2}}, \\ C_{2} = \frac{(1 - α_{1}) R^{2} - K_{1} α_{1} R}{γ K_{1} α_{1} K_{2} α_{2}} . \end{cases}

Still, we assume that the population dynamics of the resource species follows Equation S42. At the fixed point, $\dot{R} = 0$ . We have:

ζ (1 - \frac{R}{κ}) = k_{1} α_{1} C_{1} + k_{2} α_{2} C_{2} .

Combined with Equation S68, we can solve for $R$ :

R = \frac{- ϱ_{1} + \sqrt{ϱ_{1}^{2} + 4 ϱ_{2} ζ}}{2 ϱ_{2}} .

where $ϱ_{1} \equiv \frac{ζ}{κ} - \frac{k_{1}}{γ K_{1}} - \frac{k_{2}}{γ K_{2}}$ and $ϱ_{2} \equiv \frac{k_{1} (1 - α_{2})}{γ K_{1} K_{2} α_{2}} + \frac{k_{2} (1 - α_{1})}{γ K_{1} K_{2} α_{1}}$ .

Equation S68, Equation S70 are the analytical solutions of the fixed point when $R ≫ C_{1} + C_{2}$ . As shown in Appendix 1—figure 3E, the analytical predictions agree well with the numerical results (the exact solutions).

Appendix 6

Scenario involving chasing pairs and both intra- and inter-specific interference

Here, we consider the scenario involving chasing pairs and both intra- and inter-specific interference in the simple case of $S_{C} = 2$ and $S_{R} = 1$ :

\begin{array}{ll} C_{i}^{(F)} + R^{(F)} ⇌_{d_{i}}^{a_{i}} C_{i}^{(P)} \lor R^{(P)} \overset{k_{i}}{\to} C_{i}^{(F)} (+), \\ C_{1}^{(F)} + C_{2}^{(F)} ⇌_{d_{12}^{'}}^{a_{12}^{'}} C_{1}^{(P)} \lor C_{2}^{(P)}, \\ C_{i}^{(F)} + C_{i}^{(F)} ⇌_{d_{i}^{'}}^{a_{i}^{'}} C_{i}^{(P)} \lor C_{i}^{(P)}, i = 1, 2. \end{array}

We adopt the same notations as that depicted in Appendix 4.A and Appendix 5. Then, $C_{i} \equiv C_{i}^{(F)} + x_{i} + 2 y_{i} + z$ and $R \equiv R^{(F)} + x_{1} + x_{2}$ , and the population dynamics of the system can be described as follows:

{\begin{cases} {\dot{x}}_{i} = a_{i} C_{i}^{(F)} R^{(F)} - (k_{i} + d_{i}) x_{i}, \\ \dot{z} = a_{12}^{'} C_{1}^{(F)} C_{2}^{(F)} - d_{12}^{'} z, \\ {\dot{y}}_{i} = a_{i}^{'} [C_{i}^{(F)}]^{2} - d_{i}^{'} y_{i}, \\ {\dot{C}}_{i} = w_{i} k_{i} x_{i} - D_{i} C_{i}, \\ \dot{R} = g (R, x_{1}, x_{2}, C_{1}, C_{2}), i = 1, 2. \end{cases}

Here, the functional form of $g (R, x_{1}, x_{2}, C_{1}, C_{2})$ follows Equation S42. For convenience, we define $K_{i} \equiv (d_{i} + k_{i}) / a_{i}, α_{i} \equiv D_{i} / (w_{i} k_{i}), β_{i} \equiv a_{i}^{'} / d_{i}^{'}$ , and $γ \equiv a_{12}^{'} / d_{12}^{'}, (i = 1, 2)$ . At steady state, from ${\dot{x}}_{i} = 0, {\dot{y}}_{i} = 0, \dot{z} = 0$ , and ${\dot{C}}_{i} = 0, (i = 1, 2)$ , we have:

{\begin{cases} x_{i} = α_{i} C_{i}, \\ C_{i}^{(F)} = K_{i} α_{i} C_{i} (R^{(F)})^{- 1}, \\ y_{i} = β_{i} (K_{i} α_{i} C_{i})^{2} [R^{(F)}]^{- 2}, \\ z = γ K_{1} α_{1} K_{2} α_{2} [R^{(F)}]^{- 2} C_{1} C_{2} . \end{cases}

Combined with $C_{i} \equiv C_{i}^{(F)} + x_{i} + 2 y_{i} + z$ , and since $C_{i} > 0$ $(i = 1, 2)$ , then,

{\begin{cases} (1 - α_{1}) (R^{(F)})^{2} - K_{1} α_{1} R^{(F)} = 2 β_{1} (K_{1} α_{1})^{2} C_{1} + γ K_{1} α_{1} K_{2} α_{2} C_{2}, \\ (1 - α_{2}) (R^{(F)})^{2} - K_{2} α_{2} R^{(F)} = 2 β_{2} (K_{2} α_{2})^{2} C_{2} + γ K_{1} α_{1} K_{2} α_{2} C_{1} . \end{cases}

A Analytical solutions of species abundances at steady state

If $R ≫ C_{1} + C_{2}$ , then $R ≫ x_{1} + x_{2}$ and thus $R^{(F)} \approx R$ . Combined with Equation S73, we obtain:

{\begin{cases} C_{1} = R \frac{(2 β_{2} K_{2} α_{2} (1 - α_{1}) - γ K_{1} α_{1} (1 - α_{2})) R + (γ - 2 β_{2}) K_{1} α_{1} K_{2} α_{2}}{K_{1}^{2} α_{1}^{2} K_{2} α_{2} (4 β_{1} β_{2} - γ^{2})}, \\ C_{2} = R \frac{(2 β_{1} K_{1} α_{1} (1 - α_{2}) - γ K_{2} α_{2} (1 - α_{1})) R + (γ - 2 β_{1}) K_{1} α_{1} K_{2} α_{2}}{K_{1} α_{1} K_{2}^{2} α_{2}^{2} (4 β_{1} β_{2} - γ^{2})} . \end{cases}

Using $\dot{R} = 0$ and $R > 0$ , we have:

R = \frac{- χ_{2} + \sqrt{(χ_{2})^{2} + 4 χ_{1} ζ}}{2 χ_{1}},

where $χ_{1} \equiv \frac{k_{2} γ (α_{1} - 1)}{K_{1} K_{2} α_{1} (4 β_{1} β_{2} - γ^{2})} + \frac{k_{1} γ (α_{2} - 1)}{K_{1} K_{2} α_{2} (4 β_{1} β_{2} - γ^{2})} - \frac{k_{1} 2 β_{2} (α_{1} - 1)}{K_{1}^{2} α_{1} (4 β_{1} β_{2} - γ^{2})} - \frac{k_{2} 2 β_{1} (α_{2} - 1)}{K_{2}^{2} α_{2} (4 β_{1} β_{2} - γ^{2})}$ , and $χ_{2} \equiv \frac{k_{1} (γ - 2 β_{2})}{K_{1} (4 β_{1} β_{2} - γ^{2})} + \frac{k_{2} (γ - 2 β_{1})}{K_{2} (4 β_{1} β_{2} - γ^{2})} + \frac{ζ}{κ}$ . Equations S74-S75 are the analytical solutions of the species abundances at steady state when $R ≫ C_{1} + C_{2}$ . As shown in Appendix 1—figure 5E, the analytical calculations agree well with the numerical results (the exact solutions).

B Stability analysis of the coexisting state

In the scenario involving chasing pairs and both intra- and inter-specific interference, the behavior of species coexistence is similar to that without interspecific interference. Evidently, the influence of interspecific interference would be negligible if $d_{12}^{'}$ is extremely large, and vice versa for intraspecific interference if both $d_{1}^{'}$ and $d_{2}^{'}$ are tremendous. In the deterministic framework, the two-consumer species may coexist at constant population densities (Appendix 1—figure 5B), and the fixed points are globally attracting (Appendix 1—figure 5C). Furthermore, there is a non-zero measure of parameter set where both consumer species can coexist at steady state with only one type of resources (Appendix 1—figure 5A). In the stochastic framework, just as the scenario involving chasing pairs and intraspecific interference, the coexistence state can be maintained along with stochasticity (Appendix 1—figure 5D).

Appendix 7

Dimensional analysis for the scenario involving chasing pairs and both intra- and inter-specific interference

The population dynamics of the system involving chasing pairs and both intra- and inter-specific interference are shown in Equations 1-4:

\begin{array}{ll} {\begin{cases} {\dot{x}}_{i l} = a_{i l} C_{i}^{(F)} R_{l}^{(F)} - (d_{i l} + k_{i l}) x_{i l}, \\ {\dot{y}}_{i} = a_{i}^{'} [C_{i}^{(F)}]^{2} - d_{i}^{'} y_{i}, \\ {\dot{z}}_{i j} = a_{i j}^{'} C_{i}^{(F)} C_{j}^{(F)} - d_{i j}^{'} z_{i j} \\ {\dot{C}}_{i} = \sum_{l = 1}^{S_{R}} w_{i l} k_{i l} x_{i l} - D_{i} C_{i}, \\ {\dot{R}}_{l} = ζ_{l} (1 - R_{l} / κ_{l}) - \sum_{i = 1}^{S_{C}} k_{i l} x_{i l}, \end{cases} \end{array}

with $l = 1, \dots, S_{R}$ ; $i, j = 1, \dots, S_{C}$ , and $i \neq j$ . Here, $C_{i} = C_{i}^{(F)} + \sum_{l} x_{i l} + 2 y_{i} + \sum_{i \neq j} z_{i j}$ and $R_{l} = R_{l}^{(F)} + \sum_{i} x_{i l}$ represent the population abundances of the consumers and resources in the system. In fact, there are already several dimensionless variables and parameter in Equation S76, namely $x_{i l}$ , $y_{i}$ , $z_{i j}$ , $C_{i}^{(F)}$ , $R_{l}^{(F)}$ , $C_{i}$ , $R_{l}$ , $w_{i l}$ , $κ_{l}$ . To make all terms dimensionless, we define $\tilde{t} = t / τ$ , where $τ = \tilde{D_{1}} / D_{1}$ and $\tilde{D_{1}}$ is a reducible dimensionless parameter which is freely to take any positive values. Besides, we define dimensionless parameters ${\tilde{a}}_{i l} = a_{i l} τ$ , ${\tilde{d}}_{i l} = d_{i l} τ$ , ${\tilde{k}}_{i l} = k_{i l} τ$ , ${\tilde{a}}_{i}^{'} = a_{i}^{'} τ$ , ${\tilde{d}}_{i}^{'} = d_{i}^{'} τ$ , ${\tilde{a}}_{i j}^{'} = a_{i j}^{'} τ$ , ${\tilde{d}}_{i j}^{'} = d_{i j}^{'} τ$ , ${\tilde{D}}_{i} = D_{i} τ$ and ${\tilde{ζ}}_{l} = ζ_{l} τ$ . By substituting all the dimensionless terms into Equation S76, we have:

{\begin{cases} {\dot{x}}_{i l} = {\tilde{a}}_{i l} C_{i}^{(F)} R_{l}^{(F)} - ({\tilde{d}}_{i l} + {\tilde{k}}_{i l}) x_{i l}, \\ {\dot{y}}_{i} = {\tilde{a}}_{i}^{'} [C_{i}^{(F)}]^{2} - {\tilde{d}}_{i}^{'} y_{i}, \\ {\dot{z}}_{i j} = {\tilde{a}}_{i j}^{'} C_{i}^{(F)} C_{j}^{(F)} - {\tilde{d}}_{i j}^{'} z_{i j} \\ {\dot{C}}_{i} = \sum_{l = 1}^{S_{R}} w_{i l} {\tilde{k}}_{i l} x_{i l} - {\tilde{D}}_{i} C_{i}, \\ {\dot{R}}_{l} = {\tilde{ζ}}_{l} (1 - R_{l} / κ_{l}) - \sum_{i = 1}^{S_{C}} {\tilde{k}}_{i l} x_{i l} . \end{cases}

For convenience, we omit the notation ‘ ˜ ’ and use dimensionless variables and parameters in the simulation studies unless otherwise specified.

Appendix 8

Approximations applied in the pairwise encounter model

For consumers within a paired state, either in a chasing pair or an interference pair, the consumer may die following the mortality rate. Thus, in the scenario involving chasing pairs and both intra- and inter-specific interference, the population dynamics of the system should be described as follows:

{\begin{cases} {\dot{x}}_{i l} = a_{i l} C_{i}^{(F)} R_{l}^{(F)} - (d_{i l} + k_{i l} + D_{i}) x_{i l}, \\ {\dot{y}}_{i} = a_{i}^{'} [C_{i}^{(F)}]^{2} - (d_{i}^{'} + D_{i}) y_{i}, \\ {\dot{z}}_{i j} = a_{i j}^{'} C_{i}^{(F)} C_{j}^{(F)} - (d_{i j}^{'} + D_{i} + D_{j}) z_{i j} \\ {\dot{C}}_{i} = \sum_{l = 1}^{S_{R}} w_{i l} k_{i l} x_{i l} - D_{i} C_{i}, i = 1, \dots, S_{C}, \\ {\dot{R}}_{l} = ζ_{l} (1 - R_{l} / κ_{l}) - \sum_{i = 1}^{S_{C}} k_{i l} x_{i l}, l = 1, \dots, S_{R} . \end{cases}

However, since predation or interference processes are generally much faster than birth and death processes, that is $D_{i} << k_{i l}, d_{i l}, d_{i}^{'}, d_{i j}^{'}$ , the influence of mortality rate in a paired state is negligible. Therefore, we have used the following approximations throughout our manuscript: $(k_{i l} + d_{i l} + D_{i}) \approx (k_{i l} + d_{i l})$ , $(d_{i}^{'} + D_{i}) \approx d_{i}^{'}$ , $(d_{i j}^{'} + D_{i} + D_{j}) \approx d_{i j}^{'}$ . Hence, the approximated population dynamics is described as follows:

{\begin{cases} {\dot{x}}_{i l} = a_{i l} C_{i}^{(F)} R_{l}^{(F)} - (d_{i l} + k_{i l}) x_{i l}, \\ {\dot{y}}_{i} = a_{i}^{'} [C_{i}^{(F)}]^{2} - d_{i}^{'} y_{i}, \\ {\dot{z}}_{i j} = a_{i j}^{'} C_{i}^{(F)} C_{j}^{(F)} - d_{i j}^{'} z_{i j} \\ {\dot{C}}_{i} = \sum_{l = 1}^{S_{R}} w_{i l} k_{i l} x_{i l} - D_{i} C_{i}, i = 1, \dots, S_{C}, \\ {\dot{R}}_{l} = ζ_{l} (1 - R_{l} / κ_{l}) - \sum_{i = 1}^{S_{C}} k_{i l} x_{i l}, l = 1, \dots, S_{R}, \end{cases}

which is identical to those shown in the main text.

Appendix 9

Simulation details of the main text figures

In Figure 1C and F: $a_{i} = 0.1$ , $d_{i} = 0.5$ , $w_{i} = 0.1$ , $k_{i} = 0.1$ , $(i = 1, 2)$ , $D_{1} = 0.002$ , $D_{2} = 0.001$ , $κ = 5$ , $ζ = 0.05$ . In Figure 1D and G: $a_{i} = 0.02$ , $a_{i j}^{'} = 0.021$ , $d_{i} = 0.5$ , $d_{i j}^{'} = 0.01$ , $w_{i} = 0.08$ , $k_{i} = 0.03$ , $i, j = 1, 2$ , $i \neq j$ , $D_{2} = 0.001$ , $D_{1} = 0.0011$ , $κ = 20$ , $ζ = 0.01$ . In Figure 1E and H: $a_{i} = 0.5$ , $a_{i}^{'} = 0.625$ , $d_{i} = 0.5$ , $d_{i}^{'} = 0.02$ , $w_{i} = 0.2$ , $k_{i} = 0.4$ , $(i = 1, 2)$ , $D_{1} = 0.0286$ , $D_{2} = 0.022$ , $κ = 10$ , $ζ = 0.5$ . Figure 1C and F were calculated or simulated from Equations 1, 4. Figure 1D and G were calculated or simulated from Equations 1, 3 and 4. Figure 1E and H were calculated or simulated from Equations 1, 2 and 4. The analytical solutions in Figure 1E were calculated from Equations S41 and S43.

In Figure 2A: $a_{i} = 0.02$ , $a_{i}^{'} = 0.025$ , $d_{i} = 0.7$ , $d_{i}^{'} = 0.7$ , $w_{i} = 0.4$ , $k_{i} = 0.05$ , $(i = 1, 2)$ , $D_{1} = 0.0160$ , $D_{2} = 0.0171$ , $κ = 2000$ , $ζ = 5.5$ . In Figure 2B–C: $L = 100$ , $r^{(C)} = 5$ , $r^{(I)} = 5$ , $v_{C_{i}} = 1$ , $v_{R} = 0.1$ , $a_{i} = 0.2010$ , $a_{i}^{'} = 0.2828$ , $d_{i} = 0.7$ , $d_{i}^{'} = 0.8$ , $w_{i} = 0.33$ , $k_{i} = 0.2$ , $(i = 1, 2)$ , $D_{1} = 0.0605$ , $D_{2} = 0.0600$ , $κ = 1000$ . In Figure 2D: $a_{i} = 0.3$ , $a_{i}^{'} = 0.33$ , $w_{i} = 0.018$ , $k_{i} = 4.8$ , $d_{i} = 5.5$ , $d_{i}^{'} = 5$ , $(i = 1, 2)$ , $D_{1} = 0.011$ , $D_{2} = 0.010$ , $κ = 10000$ , $ζ = 35$ . In Figure 2E: $a_{i} = 0.2$ , $a_{i}^{'} = 0.24$ , $d_{i} = 4.5$ , $d_{i}^{'} = 4$ , $w_{i} = 0.02$ , $k_{i} = 4.5$ , $(i = 1, 2)$ , $D_{1} = 0.0120$ , $D_{2} = 0.010$ , $κ = 10000$ , $ζ = 35$ . In Figure 2D and E: We set $τ = 0.4$ Day (see Appendix 7). This results in an expected lifespan of Drosophila serrata in the model settings of $τ / D_{2} = 40$ days and that of Drosophila pseudoobscura $τ / D_{1} = 36.4$ days, which roughly agrees with experimental data showing that the average lifespan of D. serrata is 34 days for males and 54 days for females (Narayan et al., 2022), and the average lifespan of D. pseudoobscura is around 40 days for females (Gowaty et al., 2010). The time averages ( $（ \bar{C_{i}} ）$ ) and standard deviations ( $δ C_{i}$ ) of the species' relative/absolute abundances for the experimental data or SSA results are as follows: $\bar{^{(R)} C_{D . s e r r a t a_A R_G r p 1}^{Exp (SSA)}} = 0.53 (0.55)$ , $δ^{(R)} C_{D . s e r r a t a_A R_G r p 1}^{E x p (S S A)} = 0.12 (0.09)$ , $\bar{^{(R)} C_{D . s e r r a t a_A R_G r p 2}^{Exp (SSA)}} = 0.59 (0.61)$ , $δ^{(R)} C_{D . s e r r a t a_A R_G r p 2}^{E x p (SSA)} = 0.10 (0.12)$ , $\bar{C_{T . c o n f u s u m {_24}^{\circ} C}^{E x p (SSA)}} = 29.1 (28.6)$ , $δ C_{T . c o n f u s u m {_24}^{\circ} C}^{E x p (SSA)} = 5.4 (5.2)$ , $\bar{C_{T . c a s t a n e u m a {_24}^{\circ} C}^{E x p (SSA)}} = 45.9 (54.5)$ , where the superscript ‘(R)’ represents relative abundances. A comparison of Shannon entropies in the time series between experimental data and SSA results is presented in Appendix 1—figure 7C and D. Figure 2A–E were simulated from Equations 1, 2, and 4. See Appendix 1—figure 7E and G for the long-term time series of all species in Figure 2D and E, respectively.

Model settings in Figure 3A–B and D (plankton): $a_{i l} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i l} = 0.5$ , $d_{i}^{'} = 0.2$ , $w_{i l} = 0.3$ , $k_{i l} = 0.2$ , $κ_{1} = 8 \times 10^{4}$ , $κ_{2} = 5 \times 10^{4}$ , $κ_{3} = 3 \times 10^{4}$ , $ζ_{1} = 280$ , $ζ_{2} = 200$ , $ζ_{3} = 150$ , $D_{i} = 0.03 \times N (1, 0.25)$ , $(i = 1, \dots, S_{C}$ , $l = 1, \dots, S_{R})$ , $S_{C} = 140$ and $S_{R} = 3$ . Model settings in Figure 3C (bird): $a_{i} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i} = 0.5$ , $d_{i}^{'} = 0.5$ , $w_{i} = 0.3$ , $k_{i} = 0.2$ , $κ = 10^{5}$ , $ζ = 110$ , $D_{i} = 0.02 \times N (1, 0.28)$ , $(i = 1, \dots, S_{C})$ , $S_{C} = 250$ and $S_{R} = 1$ . Model settings in Figure 3C (fish): $a_{i} = 0.1$ , $a_{i}^{'} = 0.14$ , $d_{i} = 0.5$ , $d_{i}^{'} = 0.5$ , $w_{i} = 0.2$ , $k_{i} = 0.1$ , $κ = 10^{6}$ , $ζ = 550$ , $D_{i} = 0.015 \times N (1, 0.32)$ , $(i = 1, \dots, 45)$ , $S_{C} = 45$ and $S_{R} = 1$ . Model settings in Figure 3C (butterfly): $a_{i} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i} = 0.5$ , $d_{i}^{'} = 0.3$ , $w_{i} = 0.3$ , $k_{i} = 0.2$ , $κ = 10^{5}$ , $ζ = 300$ , $D_{i} = 0.034 \times N (1, 0.35)$ , $(i = 1, \dots, S_{C})$ , $S_{C} = 150$ and $S_{R} = 1$ . Model settings in Figure 3D (bat): $a_{i} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i} = 0.5$ , $d_{i}^{'} = 0.5$ , $w_{i} = 0.2$ , $k_{i} = 0.1$ , $κ = 10^{6}$ , $ζ = 250$ , $D_{i} = 0.013 \times N (1, 0.34)$ , $(i = 1, \dots, S_{C})$ , $S_{C} = 40$ and $S_{R} = 1$ . Model settings in Figure 3D (lizard): $a_{i} = 0.1$ , $a_{i}^{'} = 0.125$ , $d_{i} = 0.5$ , $d_{i}^{'} = 0.5$ , $w_{i} = 0.2$ , $k_{i} = 0.1$ , $κ = 10^{6}$ , $ζ = 250$ , $D_{i} = 0.014 \times N (1, 0.34)$ , $(i = 1, \dots, S_{C})$ , $S_{C} = 55$ and $S_{R} = 1$ . In Figure 3A–D, the mortality rate $D_{i}$ is the only parameter that varies with the consumer species, which was randomly sampled from a Gaussian distribution $N (μ, σ)$ , where μ and σ are the mean and standard deviation of the distribution. The coefficient of variation of the mortality rates (i.e. $σ / μ$ ) was chosen to be around 0.3, or more precisely, the best-fit in the range of 0.15–0.43. This range was estimated from experimental results (Menon et al., 2003) using the two-sigma rule. These settings for the mortality rates also apply to those in Appendix 1—figures 8–14. Figure 3A–D were simulated from Equations 1, 2 and 4. See Appendix 1—figures 10K, C, D, H, I, J and 14C, Figure 3A and B for the time series of Figure 3C ( ${ODEs}_{bird}^{S_{R} = 1}$ ), Figure 3C ( ${ODEs}_{butterfly}^{S_{R} = 1}$ ), Figure 3C ( ${ODEs}_{fish}^{S_{R} = 1}$ ), Figure 3D ( ${ODEs}_{bat}^{S_{R} = 1}$ ), Figure 3D ( ${SSA}_{bat}^{S_{R} = 1}$ ), Figure 3D ( ${ODEs}_{lizard}^{S_{R} = 1}$ ), Figure 3D ( ${SSA}_{lizard}^{S_{R} = 1}$ ), Figure 3D ( ${ODEs}_{plankton}^{S_{R} = 3}$ ) and Figure 3D ( ${SSA}_{plankton}^{S_{R} = 3}$ ), respectively. The Shannon entropies of the experimental data and simulation results for each ecological community are: $H_{E x p (O D E s)}^{b i r d (1982)} = 5.67 (6.79)$ , $H_{E x p (O D E s)}^{b i r d (2018)} = 6.63 (6.79)$ , $H_{E x p (O D E s)}^{b u t t e r f l y} = 4.78 (4.12)$ , $H_{E x p (O D E s)}^{f i s h} = 3.78 (3.40)$ , $H_{E x p (O D E s, S S A)}^{l i z a r d} = 4.05 (3.57, 3.50)$ . Here the Shannon entropy $H = - \sum_{i = 1}^{S_{C}} P_{i} \log_{2} (P_{i})$ , where $P_{i}$ is the probability that a consumer individual belongs to species $C_{i}$ .

Data availability

All data and codes for this paper are available on GitHub (copy archived at SchordK, 2024).

References

1. Arditi R
2. Ginzburg LR
(1989) Coupling in predator-prey dynamics: Ratio-Dependence
Journal of Theoretical Biology 139:311–326.

https://doi.org/10.1016/S0022-5193(89)80211-5
- Google Scholar
1. Arthur RM
(1969) Species packing, and what competition minimizes
PNAS 64:1369–1371.

https://doi.org/10.1073/pnas.64.4.1369
- PubMed
- Google Scholar
1. Ayala FJ
(1969) Experimental invalidation of the principle of competitive exclusion
Nature 224:1076–1079.

https://doi.org/10.1038/2241076a0
- PubMed
- Google Scholar
1. Beddington JR
(1975) Mutual interference between parasites or predators and its effect on searching efficiency
The Journal of Animal Ecology 44:331.

https://doi.org/10.2307/3866
- Google Scholar
(2005) The evolution of male traits in social insects
Annual Review of Entomology 50:395–420.

https://doi.org/10.1146/annurev.ento.50.071803.130416
- PubMed
- Google Scholar
(2004) Intraspecific interference and consumer-resource dynamics
Discrete & Continuous Dynamical Systems - B 4:527–546.

https://doi.org/10.3934/dcdsb.2004.4.527
- Google Scholar
1. Chen S
2. Lee AY
3. Bowens NM
4. Huber R
5. Kravitz EA
(2002) Fighting fruit flies: a model system for the study of aggression
PNAS 99:5664–5668.

https://doi.org/10.1073/pnas.082102599
- PubMed
- Google Scholar
1. Chesson P
(1990) MacArthur’s consumer-resource model
Theoretical Population Biology 37:26–38.

https://doi.org/10.1016/0040-5809(90)90025-Q
- Google Scholar
1. Chesson P
(2000) Mechanisms of maintenance of species diversity
Annual Review of Ecology and Systematics 31:343–366.

https://doi.org/10.1146/annurev.ecolsys.31.1.343
- Google Scholar
(2005) Life after logging: post‐logging recovery of a neotropical bat community
Journal of Applied Ecology 42:409–420.

https://doi.org/10.1111/j.1365-2664.2005.01024.x
- Google Scholar
Book
1. Cody ML
2. Smallwood JA
(1996)
Long-Term Studies of Vertebrate Communities

Academic Press.
- Google Scholar
(2002) Chemical warfare between microbes promotes biodiversity
PNAS 99:786–790.

https://doi.org/10.1073/pnas.012399899
- PubMed
- Google Scholar
(2021) Collective behaviour can stabilize ecosystems
Nature Ecology & Evolution 5:1435–1440.

https://doi.org/10.1038/s41559-021-01517-w
- Google Scholar
1. Daniel R
(2005) The metagenomics of soil
Nature Reviews. Microbiology 3:470–478.

https://doi.org/10.1038/nrmicro1160
- PubMed
- Google Scholar
1. Dankert H
2. Wang L
3. Hoopfer ED
4. Anderson DJ
5. Perona P
(2009) Automated monitoring and analysis of social behavior in Drosophila
Nature Methods 6:297–303.

https://doi.org/10.1038/nmeth.1310
- PubMed
- Google Scholar
(1975) A model for tropic interaction
Ecology 56:881–892.

https://doi.org/10.2307/1936298
- Google Scholar
1. de Vargas C
2. Audic S
3. Henry N
4. Decelle J
5. Mahé F
6. Logares R
7. Lara E
8. Berney C
9. Le Bescot N
10. Probert I
11. Carmichael M
12. Poulain J
13. Romac S
14. Colin S
15. Aury JM
16. Bittner L
17. Chaffron S
18. Dunthorn M
19. Engelen S
20. Flegontova O
21. Guidi L
22. Horák A
23. Jaillon O
24. Lima-Mendez G
25. Lukeš J
26. Malviya S
27. Morard R
28. Mulot M
29. Scalco E
30. Siano R
31. Vincent F
32. Zingone A
33. Dimier C
34. Picheral M
35. Searson S
36. Kandels-Lewis S
37. Tara Oceans Coordinators
38. Acinas SG
39. Bork P
40. Bowler C
41. Gorsky G
42. Grimsley N
43. Hingamp P
44. Iudicone D
45. Not F
46. Ogata H
47. Pesant S
48. Raes J
49. Sieracki ME
50. Speich S
51. Stemmann L
52. Sunagawa S
53. Weissenbach J
54. Wincker P
55. Karsenti E
(2015) Ocean plankton. Eukaryotic plankton diversity in the sunlit ocean
Science 348:1261605.

https://doi.org/10.1126/science.1261605
- PubMed
- Google Scholar
1. De Vries P
(1997) Species diversity in vertical, horizontal, and temporal dimensions of a fruit-feeding butterfly community in an Ecuadorian rainforest
Biological Journal of the Linnean Society 62:343–364.

https://doi.org/10.1006/bijl.1997.0155
- Google Scholar
1. Fuhrman JA
2. Steele JA
3. Hewson I
4. Schwalbach MS
5. Brown MV
6. Green JL
7. Brown JH
(2008) A latitudinal diversity gradient in planktonic marine bacteria
PNAS 105:7774–7778.

https://doi.org/10.1073/pnas.0803070105
- PubMed
- Google Scholar
Book
1. Gause GF
(1934)
The Struggle for Existence

Baltimore: The Williams & Wilkins Company.
- Google Scholar
1. Gillespie DT
(2007) Stochastic simulation of chemical kinetics
Annual Review of Physical Chemistry 58:35–55.

https://doi.org/10.1146/annurev.physchem.58.032806.104637
- PubMed
- Google Scholar
1. Goldford JE
2. Lu N
3. Bajić D
4. Estrela S
5. Tikhonov M
6. Sanchez-Gorostiaga A
7. Segrè D
8. Mehta P
9. Sanchez A
(2018) Emergent simplicity in microbial community assembly
Science 361:469–474.

https://doi.org/10.1126/science.aat1168
- PubMed
- Google Scholar
(2010) Polyandry increases offspring viability and mother productivity but does not decrease mother survival in Drosophila pseudoobscura
PNAS 107:13771–13776.

https://doi.org/10.1073/pnas.1006174107
- PubMed
- Google Scholar
1. Goyal A
2. Maslov S
(2018) Diversity, stability, and reproducibility in stochastically assembled microbial ecosystems
Physical Review Letters 120:158102.

https://doi.org/10.1103/PhysRevLett.120.158102
- PubMed
- Google Scholar
(2017) Higher-order interactions stabilize dynamics in competitive network models
Nature 548:210–213.

https://doi.org/10.1038/nature23273
- Google Scholar
Book
1. Grimm V
(2013)
Individual-Based Modeling and Ecology

Princeton University Press.
- Google Scholar
1. Gupta D
2. Garlaschi S
3. Suweis S
4. Azaele S
5. Maritan A
(2021) Effective resource competition model for species coexistence
Physical Review Letters 127:208101.

https://doi.org/10.1103/PhysRevLett.127.208101
- PubMed
- Google Scholar
1. Hardin G
(1960) The competitive exclusion principle
Science 131:1292–1297.

https://doi.org/10.1126/science.131.3409.1292
- PubMed
- Google Scholar
(1986) Bird community dynamics in a temperate deciduous forest: Long‐term trends at hubbard brook
Ecological Monographs 56:201–220.

https://doi.org/10.2307/2937074
- Google Scholar
1. Hoorn C
2. Wesselingh FP
3. ter Steege H
4. Bermudez MA
5. Mora A
6. Sevink J
7. Sanmartín I
8. Sanchez-Meseguer A
9. Anderson CL
10. Figueiredo JP
11. Jaramillo C
12. Riff D
13. Negri FR
14. Hooghiemstra H
15. Lundberg J
16. Stadler T
17. Särkinen T
18. Antonelli A
(2010) Amazonia through time: Andean uplift, climate change, landscape evolution, and biodiversity
Science 330:927–931.

https://doi.org/10.1126/science.1194585
- PubMed
- Google Scholar
1. Hsu SB
2. Ruan S
3. Yang TH
(2013) On the dynamics of two-consumers-one-resource competing systems withBeddington-DeAngelis functional response
Discrete & Continuous Dynamical Systems - B 18:2331–2353.

https://doi.org/10.3934/dcdsb.2013.18.2331
- Google Scholar
1. Hu J
2. Amor DR
3. Barbier M
4. Bunin G
5. Gore J
(2022) Emergent phases of ecological diversity and dynamics mapped in microcosms
Science 378:85–89.

https://doi.org/10.1126/science.abm7841
- PubMed
- Google Scholar
Book
1. Hubbell SP
(2001)
The Unified Neutral Theory of Biodiversity and Biogeography

Princeton University Press.
- Google Scholar
1. Huisman G
2. De Boer RJ
(1997) A formal derivation of the “Beddington” functional response
Journal of Theoretical Biology 185:389–400.

https://doi.org/10.1006/jtbi.1996.0318
- Google Scholar
1. Huisman J
2. Weissing FJ
(1999) Biodiversity of plankton by species oscillations and chaos
Nature 402:407–410.

https://doi.org/10.1038/46540
- Google Scholar
1. Hutchinson GE
(1961) The paradox of the plankton
The American Naturalist 95:137–145.

https://doi.org/10.1086/282171
- Google Scholar
(2015) Counteraction of antibiotic production and degradation stabilizes microbial communities
Nature 521:516–519.

https://doi.org/10.1038/nature14485
- Google Scholar
1. Koch AL
(1974) Competitive coexistence of two predators utilizing the same prey under constant environmental conditions
Journal of Theoretical Biology 44:387–395.

https://doi.org/10.1016/0022-5193(74)90169-6
- PubMed
- Google Scholar
1. Levin SA
(1970) Community equilibria and stability, and an extension of the competitive exclusion principle
The American Naturalist 104:413–423.

https://doi.org/10.1086/282676
- Google Scholar
1. Levin SA
(1974) Dispersion and population interactions
The American Naturalist 108:207–228.

https://doi.org/10.1086/282900
- Google Scholar
1. Levins R
(1979) Coexistence in a variable environment
The American Naturalist 114:765–783.

https://doi.org/10.1086/283527
- Google Scholar
1. MacArthur R
2. Levins R
(1964) Competition, habitat selection, and character displacement in a patchy environment
PNAS 51:1207–1210.

https://doi.org/10.1073/pnas.51.6.1207
- PubMed
- Google Scholar
1. MacArthur R
(1970) Species packing and competitive equilibrium for many species
Theoretical Population Biology 1:1–11.

https://doi.org/10.1016/0040-5809(70)90039-0
- Google Scholar
1. Martínez AE
2. Ponciano JM
3. Gomez JP
4. Valqui T
5. Novoa J
6. Antezana M
7. Biscarra G
8. Camerlenghi E
9. Carnes BH
10. Huayanca Munarriz R
11. Parra E
12. Plummer IM
13. Fitzpatrick JW
14. Robinson SK
15. Socolar JB
16. Terborgh J
(2023) The structure and organisation of an Amazonian bird community remains little changed after nearly four decades in Manu National Park
Ecology Letters 26:335–346.

https://doi.org/10.1111/ele.14159
- PubMed
- Google Scholar
1. McGehee R
(1977) Some mathematical problems concerning the ecological principle of competitive exclusion
Journal of Differential Equations 23:30–52.

https://doi.org/10.1016/0022-0396(77)90135-8
- Google Scholar
(2003) Mortality rates of autochthonous and fecal bacteria in natural aquatic ecosystems
Water Research 37:4151–4158.

https://doi.org/10.1016/S0043-1354(03)00349-X
- PubMed
- Google Scholar
(2022) Genetic and social contributions to sex differences in lifespan in Drosophila serrata
Journal of Evolutionary Biology 35:657–663.

https://doi.org/10.1111/jeb.13992
- PubMed
- Google Scholar
1. Niehaus L
2. Boland I
3. Liu M
4. Chen K
5. Fu D
6. Henckel C
7. Chaung K
8. Miranda SE
9. Dyckman S
10. Crum M
11. Dedrick S
12. Shou W
13. Momeni B
(2019) Microbial coexistence through chemical-mediated interactions
Nature Communications 10:2052.

https://doi.org/10.1038/s41467-019-10062-x
- PubMed
- Google Scholar
1. Park T
(1954) Experimental studies of interspecies competition II. Temperature, humidity, and competition in two species of tribolium
Physiological Zoology 27:177–238.

https://doi.org/10.1086/physzool.27.3.30152164
- Google Scholar
1. Pennisi E
(2005) What determines species diversity?
Science 309:90.

https://doi.org/10.1126/science.309.5731.90
- Google Scholar
(2017) Metabolic trade-offs promote diversity in a model ecosystem
Physical Review Letters 118:028103.

https://doi.org/10.1103/PhysRevLett.118.028103
- PubMed
- Google Scholar
(2020) Strength of species interactions determines biodiversity and stability in microbial communities
Nature Ecology & Evolution 4:376–383.

https://doi.org/10.1038/s41559-020-1099-4
- Google Scholar
(1992) Interference and generation cycles
Theoretical Population Biology 42:235–253.

https://doi.org/10.1016/0040-5809(92)90014-K
- Google Scholar
Software
1. SchordK
(2024) Intraspecific-predator-interference-promotes-biodiversity-in-ecosystems, version swh:1:rev:854433188ffefabcca853aa82914abedcde75ebb
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:550c791a4ac2b73fc0d1bb9e031f9835d5f51016;origin=https://github.com/SchordK/Intraspecific-predator-interference-promotes-biodiversity-in-ecosystems;visit=swh:1:snp:db8a3cacc29de10ae0e396b74bbfa531665cb625;anchor=swh:1:rev:854433188ffefabcca853aa82914abedcde75ebb
(2018) Ubiquitous abundance distribution of non-dominant plankton across the global ocean
Nature Ecology & Evolution 2:1243–1249.

https://doi.org/10.1038/s41559-018-0587-2
- Google Scholar
1. Sunagawa S
2. Acinas SG
3. Bork P
4. Bowler C
5. Acinas SG
6. Babin M
7. Bork P
8. Boss E
9. Bowler C
10. Cochrane G
11. de Vargas C
12. Follows M
13. Gorsky G
14. Grimsley N
15. Guidi L
16. Hingamp P
17. Iudicone D
18. Jaillon O
19. Kandels S
20. Karp-Boss L
21. Karsenti E
22. Lescot M
23. Not F
24. Ogata H
25. Pesant S
26. Poulton N
27. Raes J
28. Sardet C
29. Sieracki M
30. Speich S
31. Stemmann L
32. Sullivan MB
33. Sunagawa S
34. Wincker P
35. Eveillard D
36. Gorsky G
37. Guidi L
38. Iudicone D
39. Karsenti E
40. Lombard F
41. Ogata H
42. Pesant S
43. Sullivan MB
44. Wincker P
45. de Vargas C
46. Tara Oceans Coordinators
(2020) Tara Oceans: towards global ocean ecosystems biology
Nature Reviews Microbiology 18:428–445.

https://doi.org/10.1038/s41579-020-0364-5
- Google Scholar
(1990) Structure and organization of an amazonian forest bird community
Ecological Monographs 60:213–238.

https://doi.org/10.2307/1943045
- Google Scholar
1. Thingstad TF
(2000) Elements of a theory for the mechanisms controlling abundance, diversity, and biogeochemical role of lytic bacterial viruses in aquatic systems
Limnology and Oceanography 45:1320–1328.

https://doi.org/10.4319/lo.2000.45.6.1320
- Google Scholar
1. Vetsigian K
(2017) Diverse modes of eco-evolutionary dynamics in communities of antibiotic-producing microorganisms
Nature Ecology & Evolution 1:0189.

https://doi.org/10.1038/s41559-017-0189
- Google Scholar
(2020) Ocean currents promote rare species diversity in protists
Science Advances 6:eaaz9037.

https://doi.org/10.1126/sciadv.aaz9037
- PubMed
- Google Scholar
1. Wang X
2. Xia K
3. Yang X
4. Tang C
(2019) Growth strategy of microbes on mixed carbon sources
Nature Communications 10:1279.

https://doi.org/10.1038/s41467-019-09261-3
- PubMed
- Google Scholar
1. Wang X
2. Liu Y-Y
(2020) Overcome competitive exclusion in ecosystems
iScience 23:101009.

https://doi.org/10.1016/j.isci.2020.101009
- PubMed
- Google Scholar
(2019) Spatial ecology of territorial populations
PNAS 116:17874–17879.

https://doi.org/10.1073/pnas.1911570116
- PubMed
- Google Scholar
1. Xue C
2. Goldenfeld N
(2017) Coevolution maintains diversity in the stochastic “kill the winner” model
Physical Review Letters 119:268101.

https://doi.org/10.1103/PhysRevLett.119.268101
- PubMed
- Google Scholar

Article and author information

Author details

Ju Kang

School of Physics, Sun Yat-sen University, Guangzhou, China

Contribution
Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Writing – original draft

Contributed equally with
Shijie Zhang

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-3862-4065
Shijie Zhang
1. School of Mathematics, Sun Yat-sen University, Guangzhou, China
2. Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, United States
Contribution
Software, Investigation, Visualization

Contributed equally with
Ju Kang

Competing interests
No competing interests declared
Yiyuan Niu

School of Physics, Sun Yat-sen University, Guangzhou, China

Contribution
Data curation, Software, Validation, Investigation, Visualization

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0009-0008-6078-9265
Fan Zhong

School of Physics, Sun Yat-sen University, Guangzhou, China

Contribution
Investigation

Competing interests
No competing interests declared
Xin Wang

School of Physics, Sun Yat-sen University, Guangzhou, China

Contribution
Conceptualization, Data curation, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Writing – original draft, Project administration, Writing – review and editing

For correspondence
wangxin36@mail.sysu.edu.cn

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-6479-395X

Funding

National Natural Science Foundation of China (12004443)

Xin Wang

Guangzhou Municipal Science and Technology Bureau (202102020284)

Xin Wang

Sun Yat-sen University (The Hundred Talents Program)

Xin Wang

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Roy Kishony, Eric D Kelsic and Yang-Yu Liu for helpful discussions. This work was supported by National Natural Science Foundation of China (Grant No.12004443), Guangzhou Municipal Innovation Fund (Grant No. 202102020284) and the Hundred Talents Program of Sun Yat-sen University.

Version history

Preprint posted: July 19, 2023
Sent for peer review: October 10, 2023
Reviewed Preprint version 1: January 11, 2024
Reviewed Preprint version 2: August 5, 2024
Version of Record published: October 30, 2024

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.93115. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

535

views
52

downloads
0

citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Article PDF

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Ju Kang
Shijie Zhang
Yiyuan Niu
Fan Zhong
Xin Wang

(2024)

Intraspecific predator interference promotes biodiversity in ecosystems

eLife 13:RP93115.

https://doi.org/10.7554/eLife.93115.3

Share this article

Cite this article

A model of pairwise encounters may naturally break CEP.

Intraspecific predator interference facilitates species coexistence regardless of stochasticity.

Intraspecific interference enables a wide range of consumer species to coexist with only one or a handful of resource species.

Illustrations of symbols in our generic model of pairwise encounters.

Illustrations of other symbols used in our manuscript.

Estimation of the encounter rates with the mean-field approximations.

Functional response in scenarios involving different types of pairwise encounter.

Numerical solutions in scenarios involving chasing pairs and different types of predator interference.

Intraspecific predator interference facilitates species coexistence regardless of stochasticity.

Outcomes of multiple consumers species competing for one resource species involving chasing pairs and intra- and inter-specific interference.

The influence of stochasticity on species coexistence.

A model of intraspecific predator interference explains two classical experiments that invalidate CEP.

A model of intraspecific interference semi-quantitatively illustrates the rank-abundance curve of a plankton community (SC≫SR).

A model of intraspecific interference illustrates the rank-abundance curves across different ecological communities (SC≫SR).

Time courses of the species abundances in the scenario involving chasing pairs and intraspecific interference.

A model of intraspecific interference illustrates the rank-abundance curves across different ecological communities (SC≫SR).

A model of intraspecific interference illustrates the rank-abundance curves across different plankton communities (SC≫SR).

Time courses of the species abundances in the scenario involving chasing pairs and intraspecific interference.

A model of intraspecific interference illustrates the rank-abundance curve of a bird community (SC≫SR).

Intraspecific interference results in an underlying negative feedback loop and thus promotes biodiversity.

Author details

Ju Kang

Contribution

Contributed equally with

Competing interests

Shijie Zhang

Contribution

Contributed equally with

Competing interests

Yiyuan Niu

Contribution

Competing interests

Fan Zhong

Contribution

Competing interests

Xin Wang

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Further reading

A model of intraspecific interference semi-quantitatively illustrates the rank-abundance curve of a plankton community ( $S_{C} ≫ S_{R}$ ).

A model of intraspecific interference illustrates the rank-abundance curves across different ecological communities ( $S_{C} ≫ S_{R}$ ).

A model of intraspecific interference illustrates the rank-abundance curves across different ecological communities ( $S_{C} ≫ S_{R}$ ).

A model of intraspecific interference illustrates the rank-abundance curves across different plankton communities ( $S_{C} ≫ S_{R}$ ).

A model of intraspecific interference illustrates the rank-abundance curve of a bird community ( $S_{C} ≫ S_{R}$ ).