# 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 type 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 steady state. Here, we present a mechanistic model and show that the 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. The 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 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 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. The model matches observed rank-abundance curves in several natural ecosystems. However, a more detailed synthesis of relevant prior studies is needed to clarify the contribution of this manuscript in the context of existing knowledge.

# 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; 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 2,000-18,000 (Daniel, 2005). In the photic zone of the world ocean, there are roughly 150,000 eukaryotic plankton species (Vargas et al., 2015). Explaining the 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, the number of consumer species cannot exceed that of resources at steady state (MacArthur and Levins, 1964; Levin, 1970; McGehee and Armstrong, 1977). On the contrary, in the paradox of plankton, a limited type of resources support 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 (Lczárán et al., 2002), cross-feeding (Goyal and Maslov., 2018; Goldford et al., 2018; Niehaus et al., 2019), spatial circulation (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 the plankton (Pennisi, 2005; Sunagawa et al., 2020).

In this work, we present a mechanistic model of predator interference that extends the classical Beddington-DeAngelis (B-D) phenomenological model (Beddington, 1975; DeAngelis et al., 1975) with 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 (Devries et al., 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 encounter

Predator interference, i.e., the pairwise encounter among consumer individuals, is commonly described by the B-D model (Beddington, 1975; DeAngelis et al., 1975). From the mechanistic perspective, however, the functional response of the B-D model can be formally derived from the consumption between a consumer species and a resource species without involving any form of predator interference (Wang and Liu, 2020; Huisman and Boer, 1997) (see Appendix II B). To resolve this issue, we consider a mechanistic model of pairwise encounters (Fig. 1A). Specifically, ** S_{C}** consumer species competing for

**resource species . The consumers are biotic, while the resources can be either biotic or abiotic. For simplicity, we assume that all species are motile and each moves at a certain speed, namely, for consumer species**

*S*_{R}**(**

*C*_{i}**= 1, …,**

*i***) and for resource species**

*S*_{C}**, (**

*R*_{l}**= 1, …,**

*l***). For abiotic resources, they cannot propel themselves, yet may passively drift due to environmental factors. Each consumer is free to feast on one or multiple types of resources, while consumers do not directly interact with one another other than pairwise encounters.**

*S*_{R}Then, we proceed to explicitly consider the population structure of the consumers and resources: some are wandering around freely, taking Brownian motions; others are encountering with one another, forming ephemeral entangled pairs. Specifically, when a consumer individual ** C_{i}** and a resource

**, get close in space within a distance of (Fig. 1A), the consumer can chase the resource and form a chasing pair: (Fig. 1 B), where the superscript “(P)” represents pair. The resource can either escape with rate**

*R*_{l}**, or be caught and consumed by the consumer with rate**

*d*_{il}**. Meanwhile, when a**

*k*_{il}**individual encounters another consumer**

*C*_{i}**(**

*C*_{j}**= 1, …,**

*j***) of the same or different species within a distance of (Fig. 1A), they can stare at, fight against or play with each other and thus form an interference pair: (Fig. 1 B). This paired state is evanescent, and the two consumers separate from each other with rate**

*S*_{C}**. For simplicity, we assume that all and are identical, respectively, i.e., ∀**

*d’*_{il}*i*,

*j*,

*l*, and .

In a well-mixed system with the size of *L*^{2} (Appendix-fig. 1), the encounter rates among the species, ** a_{il}**, and

**(Fig. 1B), can be obtained using the mean-field approximation: and (see Materials and Methods, and Appendix-fig. 1 for details). Then, we proceed to analyze scenarios involving different types of pairwise encounters. For the scenario involving only chasing pair, the population dynamics can be described as follows:**

*a’*_{ij}where , *g _{l}* is an unspecified function, the superscript “(F)” represents the freely wandering population,

**denotes the mortality rate of**

*D*_{i}**, and**

*C*_{i}**is the mass conversion ratio (Wang and Liu, 2020) from resource**

*w*_{il}**to consumer**

*R*_{i}**. With the integration of intraspecific predator interference, we combine Eq. 1 and the following equation:**

*C*_{i}where ** a’_{i}** =

**,**

*a’*_{ii}**=**

*d’*_{i}**, and represents the intraspecific interference pair. For the scenario involving chasing pair and interspecific interference, we combine Eq. 1 with the following equation:**

*d’*_{ii}where stands for the interspecific interference pair. In the scenario where chasing pair and both intra- and inter-specific interference are all relevant, we combine Eqs. 1–3, and the populations of consumers and resources are given by and , respectively.

Generically, the consumption and interference processes are much quicker compared to the birth and death processes. Thus, in derivation of the functional response, ** ℱ(R_{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**

*k*_{ll}x_{il}/C_{i}**≈ 0 and (Appendix-fig. 2, see Appendix II for details).**

*d*_{il}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 (MacArthur, 1970; Chesson, 1990). Then,

In the absence of consumers, biotic resources exhibit logistic growth. Here, and represent the intrinsic growth rate and the carrying capacity of species ** R_{l}**. For abiotic resources, stands for the external resource supply rate of

**, and is the abundance of**

*R*_{l}**at 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 VI). For convenience, we retain the same notations below, with all parameters considered dimensionless unless otherwise specified.**

*R*_{l}## 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

**= 1. To simplify the notations, we omit the subscript/superscript “**

*S*_{R}**” since**

*l***= 1. For clarity, we assign each consumer species of unique competitiveness by setting that the mortality rate**

*S*_{R}**is the only parameter that varies with the consumer species.**

*D*_{i}First, we conduct the analysis in a deterministic framework with ordinary differential equations (ODEs). In the scenario involving only chasing pair, consumer species cannot coexist at 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 (** Ċ_{i}** = 0) yield

**(**

*f*_{i}

*R*^{(F)}) ≡

*R*^{(F)}/(

*R*^{(F)}+

**) =**

*K*_{i}**(**

*D*_{i}**= 1,2), with**

*i***= (**

*K*_{i}**+**

*d*_{i}**)/**

*k*_{i}**, which corresponds to two parallel surfaces in the (**

*a*_{i}

*C*_{1},

*C*_{2},

**) coordinates, making steady coexistence impossible (Wang and Liu, 2020) (Fig. 1C, F).**

*R*Meanwhile, in the scenario involving chasing pair and interspecific interference, if all species coexist, the steady-state equations correspond to three non-parallel surfaces Ω’_{i} (*R,C*_{1}, *C*_{2}) = ** D_{i}** (

**= 1,2),**

*i***’(**

*G*

*R, C*_{1},

*C*_{2}) = 0 (Fig. 1G and Appendix-fig. 3C, see Appendix IV for details), which can intersect at a common point (fixed point). However, this fixed point is unstable (Fig. 1G, Appendixfig. 3

*A*), and thus one of the consumer species is doom to extinct (Fig. 1 D).

Next, we turn to the scenario involving chasing pair and intraspecific interference. Likewise, steady coexistence requires that three non-parallel surfaces Ω_{i} (*R, C*_{1}, *C*_{2}) = ** D_{i}** (

**= 1, 2),**

*i***(**

*G*

*R, C*_{1},

*C*_{2}) = 0 (Fig. 1H and Appendix-fig. 3D, see Appendix III for details) cross at a common point. Indeed, this naturally happens, and encouragingly the fixed point can be stable. Therefore, two consumer species may stably coexist at steady state with only one type of resources, which obviously breaks CEP (Fig. 1E Appendix-fig. 4A). In fact, the coexisting state is globally attractive (Appendix-fig. 4A), and there exists a non-zero volume of parameter space where the two consumer species stably coexist at constant population densities (Appendix-fig. 4B), demonstrating that the violation of CEP does not depend on special parameter settings. We further consider the scenario involving chasing pair and both intra- and inter-specific interference (Appendix-fig. 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 (Appendix-fig. 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 pair and intraspecific interference, species may coexist indefinitely in the SSA simulations (Fig. 2A and Appendix-fig. 4C). In fact, the parameter region for species coexistence in this scenario is rather similar between the SSA and ODEs studies (Appendixfig. 6C-D). Similarly, in the scenario involving chasing pair and both inter- and intra-specific interference, all species may coexist indefinitely in company with stochasticity (Appendix-fig. 5D).

To further mimic a real ecosystem, we resort to individual-based modeling (IBM) (Grimm and Railsback, 2013; Vetsigian, 2017), an essentially stochastic simulation method. In the simple case of ** S_{C}** = 2 and

**= 1, we simulate the time evolution of a 2-D squared system in a size of**

*S*_{R}

*L*^{2}with periodic boundary conditions (see Materials and Methods for details). In the scenario involving chasing pair and intraspecific interference, two consumer species coexist for long with only one type of resources in the IBM simulations (Fig. 2B-C). Together with the SSA simulation studies, it is obvious that intraspecific interference promotes species coexistence along with stochasticity.

## 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 pair and intraspecific interference to simulate the two systems. Overall, our SSA results show good consistency with those of the experiments (Fig. 2D-E, see also Appendix-figs. 6C-D, 7). The fluctuations in experimental time series can be mainly accounted by stochasticity.

## A handful of resource species can support an unexpected 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}** consumers species compete for

**resource species (with**

*S*_{R}**) within the scenario involving chasing pair and intraspecific interference. The population dynamics is described by equations combining Eqs. 1, 2, 4. As with the cases above, each consumer species is assigned a unique competitiveness through a distinctive**

*S*_{C}> S_{R}**(**

*D*_{i}**= 1, —,**

*i***).**

*S*_{C}Strikingly, a plethora of consumer species may coexist at steady state with only one resource species (** S_{C} ≫ 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 are exemplified in Fig. 3 and Appendix-fig. 8–10, involving simulations with or without stochasticity. The number of consumer species in long-term coexistence can be up to hundreds or more (Fig. 3 and Appendixfig. 8). To mimic the real ecosystems, we further analyze the cases with more than one type of resources, such as systems with**

*S*_{R}**= 3 (**

*S*_{R}**). Just like the case of**

*S*_{C}≫ S_{R}**= 1 (**

*S*_{R}**), an extensive range of consumer species may coexist indefinitely regardless of stochasticity (Fig. 3 and Appendix-fig. 11–14).**

*S*_{C}≫ S_{R}## 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., ), the functional response

**≡**

*ℱ***and the steady-state population of each consumer and resource species can be obtained analytically (see Appendix III B-C for details). In fact, the functional response of a consumer species (e.g.,**

*k*_{i}x_{i}/C_{t}**) is negatively correlated with its own population size:**

*C*_{i}where ** β** ≡

**. The analytical steady-state solutions are highly consistent with the numerical results (Fig. 1E and Appendix-fig. 3E-F) and can even quantitatively predict the coexistence region of the parameter space (Appendix-fig. 3F).**

*a’*_{i}/d’_{i}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

**increases during competition, a larger portion of**

*C*_{i}**individuals are then engaged in intraspecific interference pairs which are temporarily absent from hunting (see Eq. S59 and Appendix-fig. 15A-B). Consequently, the fraction of**

*C*_{i}**individuals within chasing pairs decreases (see Eq. S59 and Appendix-fig. 15A-B) and thus form a self-inhibiting negative feedback loop through the functional response (see Eq. 5 and Appendix-fig. 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 III C for details).**

*C*_{i}## 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; Martinez et al., 2023; Clarke et al., 2005; Hubbell, 2001; Devries et al., 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 pair 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 round 0.3 (Menon et al., 2003) (see Appendix VII 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 (Fig. 3C-D, see also Appendix-figs. 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., Appendix-fig. 8C): those with experimental abundances less than 10 tend to 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; Lczaran et al., 2002; Goyal and Maslov., 2018; Goldford et al., 2018; 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; Devries et al., 1997).

In fact, predator interference has been introduced long ago by the B-D 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 pair without predator interference (Wang and Liu, 2020; Huisman and Boer , 1997) (see Eqs. S8 and S24). Therefore, it is questionable regarding the validity of applying the B-D model to break CEP. From mechanistic perspective, we resolve these issues and show that B-D model corresponds to a special case of our mechanistic model yet without the escape rate (Appendix-fig. 2, see Appendix II 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 evolution 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.

# Methods and Materials

## Derivation of the encounter rates with the mean-field approximation

In the model depicted in Fig. 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}**= 1, —,**

*i***) moves at speed with velocity , while a resource individual of species**

*S*_{C}**(**

*R*_{l}**= 1, …,**

*l***) moves at speed with velocity . Here and are two invariants, while the directions of and change constantly. The relative velocity between the two individuals is , with a relative speed of . Then, , where represents the angle between and . This system is homogeneous, thus, , where the overline stands for the temporal average. Then, we obtain the average relative speed between the**

*S*_{R}**and**

*C*_{i}**, individuals: . Likewise, the average relative speed between the**

*R*_{l}**and**

*C*_{i}**individuals is . Evidently, . Meanwhile, the concentrations of species**

*C*_{j}**and**

*C*_{i}**, in a squared system with a length of**

*R*_{l}**are and , while those of the freely wandering**

*L***and**

*C*_{i}**, are and .**

*R*_{l}Then, we use the mean-field approximation to calculate the encounter rates ** a_{il}** and

**in the well-mixed system. In particular, we estimate**

*a’*_{ij}**by tracking a randomly chosen consumer individual from species**

*a*_{il}**and counting its encounter frequency with the freely wandering individuals from resource species**

*C*_{i}**, (Appendix-fig. 1). At any moment, the consumer individual may form a chasing pair with a**

*R*_{l}**, individual within a radius of (Fig. 1A). Over a time interval of Δ**

*R*_{l}**, the number of encounters between the consumer individual and**

*t***individuals can be estimated by the encounter area and the concentration , which takes the value of (see Appendix-fig. 1). Combined with , for all freely wandering**

*R*_{l}**individuals, the number of their encounters with**

*C*_{i}

*R*^{(F)}during interval Δ

*t*is . Meanwhile, in the ODEs, this corresponds to . Comparing both terms above, for chasing pair, we have = . Likewise, for interference pair, we obtain In particular, .

## Stochastic simulations

To investigate the impact of stochasticity on species coexistence, we use stochastic simulation algorithm (SSA) (Gillespie, 2007) and individual-based modeling (IBM) (Vetsigian, 2017; Grimm and Railsback, 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 squared system in a length of ** L** with periodic boundary conditions. In the case of

**= 2 and**

*S*_{C}**= 1, consumers of species**

*S*_{R}**, (**

*C*_{i}*i*= 1, 2) move at speed , while the resources move at speed

**. The unit length is Δ**

*v*_{R}**= 1, while all individuals move probabilistically. Specifically, when Δ**

*l***is small so that ,**

*t***individuals jump a unit length with the probability . In practice, we simulate the temporal evolution of the model system following the procedures below.**

*C*_{i}### Initialization

We choose the initial position for each individual randomly from a uniform distribution in the squared space, which round to the nearest integer point in the ** x**-

**coordinates.**

*y*### Moving

We choose the destination of a movement randomly from four directions (** x**-positive,

**-negative,**

*x***-positive,**

*y***-negative) following a uniform distribution. The consumers and resources jump a unit length with probabilities and , respectively.**

*y*### Forming pairs

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

**and**

*C*_{i}**stand within a distance of**

*C*_{j}

*r*^{(I)}, they form an interference pair.

### Dissociation

We update the system with a small time step Δ** t** so that

**Δ**

*d*_{i}**,**

*t***Δ**

*k*_{i}**,**

*t***Δ**

*d’*_{ij}**≪ 1 (**

*t***= 1,2). In practice, a random number**

*i, j***is sampled from a uniform distribution between 0 and 1, i.e., U(0, 1). If**

*ζ***<**

*ζ***Δ**

*d*_{i}**or**

*t***<**

*ζ***Δ**

*d’*_{j}**, then, the chasing pair or interference pair dissociates into two separated individuals. One occupies the original position, while the other individual gets just out of the encounter radius in a uniformly distributed random angle. For a chasing pair, if**

*t***Δ**

*d*_{i}**< ζ < (**

*t***+**

*d*_{i}**)Δ**

*k*_{i}**, then, the consumer catch 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**

*t***> (**

*ζ***+**

*d*_{i}**)Δ**

*k*_{i}**or**

*t***>**

*ζ***Δ**

*d’*_{j}**, the chasing pair or interference pair maintain the current status.**

*t*### Birth and death

For each species, we use two separate counters with decimal precision to record the contributions of the birth and death processes, which both 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.

# Acknowledgements

We thank Roy Kishony, Eric D. Kelsic, Yang-Yu Liu and Fan Zhong 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.

# Data and materials availability

All study data are included in the article and/or appendices.

# Appendix I 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*= 1, i.e., two consumer species

_{R}*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:

Here *C _{i}* and

*R*represent the population abundances of consumers and resources, respectively, while the functional forms of

*f*and

_{i}(R)*g(R, C*

_{1},

*C*

_{2}) are unspecific.

*D*

_{i}stands for the mortality rate of the species

*C*. If all consumer species can coexist at steady state, then

_{i}*f*(

_{i}*R*)/

*D*= 1 (

_{i}*i*= 1, 2). In a 2-D representation, this requires that three lines

*y*=

*f*(

_{i}*R*)/

*D*(

_{i}*i*= 1, 2) 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*= 1 (

_{i}*i*= 1, 2) are parallel, and hence do not share a common point (see Ref. (Wang and Liu, 2020) for details).

# Appendix II 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*= 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 th and the wasting time

_{R}*t*. Both can be determined by specifying the scenario and using statistical physics modeling analysis. In fact, Beddington analyzed the searching efficiency Ξ

_{w}_{B-D}rather than the functional response

*ℱ*

_{B-D}, yet both can be reciprocally derived with Ξ

_{B-D}≡

*ℱ*

_{B-D}/

*R*. Here

*R*stands for the population abundance of the resources, and the specific form of Ξ

_{B-D}is (Beddington, 1975):

where *C*’ = *C* − 1, and *C* stands for the population abundance of the consumes. Generally, *C* ≫ 1, and thus *C*’ ≈ *C*.

## B Scenario involving only chasing pair

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*= 1). When an individual consumer is chasing a resource, they form a chasing pair:

_{R}where the superscript “(F)” stands for populations that are freely wandering, and “(+)” signifies gaining biomass (we count *C*^{(F)} (+) as *C*^{(F)}. *C*^{(P)} ∨ *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* ≡ C^{(F)} + *x* and *R* ≡ *R*^{(F)} + *x*. Then, the population dynamics of the system follows:

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., *ẋ* = 0). Then, we can solve for *x* with:

where , and then,

By definition, the functional response and search efficiency are:

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

Since , using first order approximations in Eq. S7, we obtain . Then the functional response and search efficiency are:

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 Eq. S2). Meanwhile, using first order approximations in the denominator of Eq. S5, we have . Hence,

In the case that *R* ≫ *C*, then *R* ≫ *C* > *x* = *R* − *R*^{(F)}. By applying *R* ≈ *R*^{(F)} in Eq. S3, we obtain . Then,

To compare these functional responses with that of the B-D model, we determine the parameters *t _{h}* and

*t*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 and (in the chasing-pair scenario, a’ = 0). By substituting these calculations into Eq. S2, we have

_{w}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-fig. 2A-B for the comparison).

## C Scenario involving chasing pair and intraspecific interference

Here we consider the scenario with additional involvement of intraspecific interference in the simple case of *S _{C}* = 1 and

*S*= 1

_{R}Here *C*^{(P)} ∨ *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* ≡ *C*_{(F)} + *x* + 2*y* and *R* ≡ *R*_{(F)} + *x*. Hence the population dynamics of the consumers and resources can be described as follows:

The consumption process and interference process are supposed to be in fast equilibrium (i.e., *ẋ* = 0, *ຏ* = 0), then we can solve for *x* with:

where *ϕ*_{0} = —*CR*^{2}, *ϕ*_{1} = 2*CR* + *KR* + *R*^{2}, *ϕ*^{2} = 2*βK*^{2} — *K* — *C* — 2*R*, with *β* = *a*’/*d*’. The discriminant of Eq. S13 (denoted as ∧) is

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

where (i stands for the imaginary unit), , and . On the other hand, when ∧ > 0, there are three real solutions *x*_{(1)}, *x*_{(2)}, and *x*(_{3}), which are

where *ψ*’ = (-4*ψ*/3)^{1/2} and *φ*’ = arccos(-(-*ψ*/3)^{-3/2}*φ*/2)/3. Note that *x* ∈ [0, min(*R, C*)], then we obtain the exact feasible solution of *x* (denoted as *x _{ext}*), and hence the functional response and search efficiency are

In the case of *R* ≫ *C*, then *R* − *R*^{(F)} = *x* < *C* *C* ≪ *R*, and thus *R*^{(F)} ≈ *R*. Still, the consumption process is supposed to be in fast equilibrium (i.e., *ẋ* = 0, *ຏ* = 0), and then we obtain

Consequently,

When or 8*βC*/(1 + *R/K*)^{2} ≫ 1, using first order approximations in the denominator of Eq. S18, we have

and then,

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

and thus,

Meanwhile, the B-D model only fits to the cases with *d* = 0. By calculating the average values of *t _{h}* and

*t*in the stochastic framework, we have , . Thus, we obtain the searching efficiency and functional response in the B-D model:

_{w}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-fig. 2C-D). Still, there is a region where the discrepancies can be small, namely d ≈ 0 and *R* ≫ *C* (Appendix-fig. 2C-D). Intuitively, when and d = 0, then Consequently, if , then . In this case, the difference between 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 pair and interspecific interference

Next, we consider the scenario involving chasing pair and interspecific interference in the case of *S _{C}* = 2 and

*S*= 1:

_{R}Here 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 and . The population dynamics of the consumers and resources follows:

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, i.e.,

*ẋ*= 0,

_{i}*ż*= 0. In the case that

*R*≫

*C*

_{1}+

*C*

_{2}>

*x*

_{1}+

*x*

_{2}, by applying

*R*

^{(F)}≈

*R*, we obtain

Then, the searching efficiencies and functional responses are:

Since , by applying first order approximation to the denominator of Eq. S26, we obtain:

and the searching efficiencies and functional responses are

Likewise, the B-D model only fits to cases with d = 0. By calculating the average values in a stochastic framework, we obtain , (*i* = 1, 2). Then, we obtain the searching efficiencies in the B-D model:

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

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-fig. 2E-F). Still, the discrepancy can be small when *d* ≈ 0 and *R* ≫ *C* (Appendix-fig. 2E-F). Intuitively, when , we have

Thus, if (*i* = 1, 2), then . In this case, the difference between and is small.

# Appendix III Scenario involving chasing pair and intraspecific interference

## A Two consumers species competing for one resource species

We consider the scenario involving chasing pair and intraspecific interference in the simple case of *S _{C}* = 2 and

*S*= 1:

_{R}Here, the variables and parameters are just extended from the case of *S _{C}* = 1 and

*S*= 1 (see Appendix II. C). The total number of consumers and resources are and . Hence, the population dynamics of the consumers and resources can be described as follows:

_{R}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}* ≡ (

*d*+

_{i}*k*)/

_{i}*a*,

_{i}*α*≡

_{i}*D*/(

_{i}*w*), and

_{i}k_{i}*β*=

_{i}*a’*/

_{i}*d’*(

_{i}*i*= 1,2). At steady state, from

*ẋ*= 0,

_{i}*ຏ*= 0, we have

_{i}Note that , and . Then,

By substituting Eq. S35a into Eq. S35b, we have

Then, we can present with *C*_{1}, *C*_{2} and *R* (*i* = 1, 2). By further combining with Eqs. S34, S35a and S36a, we express *R*^{(F)}, *x _{i}*, and

*y*using

_{i}*C*

_{1},

*C*

_{2}and

*R*. In particular, for

*x*, we have:

_{i}If all species coexist, then the steady-state equations of *Ċ _{i}* = 0 (

*i*= 1, 2) and

*Ṙ*= 0 are:

where *G*(*R,C*_{1}, *C*_{2}) ≡ *g*(*R,u*_{1}(*R,C*_{1}, *C*_{2}),*u*_{2}(*R,C*_{1}, *C*_{2}),*C*_{1}, *C*_{2}), and . In practice, Eq. S38 corresponds to three unparallel surfaces, which share a common point (Fig. 1H and Appendix-fig. 3*D*). 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, …, 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 Δ ≡ (*D*_{1} − *D*_{2})/*D*_{2} reflects the completive difference between the two consumer species. As shown in Appendix-fig. 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 *ẋ _{i}* =

*ຏ*=

_{i}*Ċ*= 0 (

_{i}*i*=1, 2), then,

Meanwhile, , and *C _{i}*,

*R*> 0 (

*i*= 1, 2). Then, we have

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)} ≈ *R*. Thus,

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

Since *Ṙ* = 0, then

where and .

Eqs. S41, S43 are the analytical solutions of species abundances at steady state when *R* ≫ *C*_{1} + *C*_{2}. As shown in Fig. 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}* (

*i*= 1, 2) to be the only parameter varying with species

*C*

_{1}and

*C*

_{2}(

*D*

_{1}>

*D*

_{2}), and define Δ ≡ (

*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-fig. 3E. 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} (*i* = 1, 2) is the only parameter that varies with the two-consumer species, the supremum of the tolerated competitive difference for species coexistence (defined as ) 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 and *C*_{1} = 0^{+}, we further combine Eq. S41 and then obtain (note that *R* > 0)

Meanwhile,

Combining Eqs. S43–S45, we have

where . When *R* ≫ *C*_{1} + *C*_{2}, the comparison of obtained from analytical solutions with that from numerical results (the exact solutions) are shown in Appendix-fig. 3F, which overall exhibits good consistency.

## B *S*_{C} consumers species competing for *S*_{R} resources species

_{C}

_{R}

Here we consider the scenario involving chasing pair and intraspecific interference for the generic case with *S _{C}* types of consumers and

*S*types of resources. Then, the population dynamics of the system can be described as follows:

_{R}Note that Eq. S47 is identical with Eqs. 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 and . For convenience, we define and .

### 1 Analytical solutions of species abundances at steady state

At steady state, from , *ẏ _{i}* =0, and , we have,

Meanwhile , and note that *C _{i}* > 0, thus

Combined with Eq. S49, and then

We further assume that the specific function of satisfies Eq. 4, i.e.,

By combining Eqs. S48, S49 and S51, we have

If the population abundance of each resource species is much more than the total population of all consumers (i.e., , then and . Thus,

with l = 1, … , *S _{R}*. Eq. S53 is a set of second-order algebraic differential equations, which is clearly solvable.

In fact, when *S _{R}* = 1,

*S*≥ 1, and , 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

_{C}*S*= 1.Then, we have

_{R}Here and .

## 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 Eqs. S47 and S51 with

*S*= 1. To simplify the notations, we omit the “l” in the subscript since

_{R}*S*= 1. The consumption process and interference process are supposed to be in fast equilibrium (i.e.,

_{R}*ẋ*= 0,

_{i}*ẏ*= 0). Then, we have a set of equations to solve for

_{i}*x*and

_{i}*y*given the population size of each species:

_{i}In the first three sub-equations of Eq. S55, by getting rids of , we have,

Then, by regarding *R*^{(F)} as a temporary parameter, we solve for *x _{i}* and

*y*:

_{i}If the total population size of the resources is much larger than that of consumers (i.e., ), then and R^{(F)} ≈ R, and thus we get the analytical expressions of *x _{i}* and

*y*

_{i}Note that the fraction of *C _{i}* individuals engaged in chasing pairs is

*x*, while that for individuals trapped in intraspecific interference pairs is

_{i}/C_{i}*y*. With Eq. S58, it is straightforward to obtain these fractions:

_{i}/C_{i}where both *x _{i}/C_{i}* and

*y*are bivariate functions of

_{i}/C_{i}*R*and

*C*. From Eq. S59, it is clear that for a given population size of the resource species,

_{i}*y*is a monotonously increasing function of

_{i}/C_{i}*C*, while

_{i}*x*is a monotonously decreasing function of

_{i}/C_{i}*C*. In Appendix-fig. 15A-B, we see that the analytical results are highly consistence with the exact numerical solutions. By definition, the functional response of

_{i}*C*species is

_{i}*ℱ*≡

*k*, and thus,

_{i}x_{i}/C_{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-fig. 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*’s population size, a larger portion of

_{i}*C*individuals are then engaged in intraspecific interference pair which are temporarily absent from hunting (Appendix-fig. 15A-B). Consequently, the functional response of

_{i}*C*drops, which prevents further increase of

_{i}*C*‘s population size, results in an overall balance among the consumer species, and thus promotes species coexistence.

_{i}# Appendix IV Scenario involving chasing pair and interspecific interference

Here we consider the scenario involving chasing pair and interspecific interference in the case of *S _{C}* = 2 and

*S*= 1 , with all settings follow that depicted in Appendix II. D. Then, , and the population dynamics follows (identical with Eq. S25):

_{R}Here the functional form of *g*(*R, x*_{1}, *x*_{2}, *C*_{1}, *C*_{2}) is unspecified. For convenience, we define *K _{i}* ≡ (

*d*+

_{i}*k*)/

_{i}*a*,

_{i}*α*≡

_{i}*D*/(

_{i}*w*)(

_{i}k_{i}*i*= 1, 2), and

*γ*=

*a*’

_{12}/

*d*’

_{12}. At steady state, from

*ẋ*= 0(

_{i}*i*= 1,2) and

*ż*= 0, we have

Note that and *R* ≡ *R*^{(F)} + *x*_{1} + *x*_{2}, then,

Then, we can express , and *R ^{(F)}* with

*C*

_{1},

*C*

_{2}and

*R*. Combined with Eq. S62,

*x*and

_{i}*z*can also be expressed using

*C*

_{1},

*C*

_{2}and

*R*. In particular, for

*x*, we have

_{i}If all species coexist, by defining , then, the steady-state equations of *Ċ _{i}* = 0 (

*i*= 1,2) and

*Ṙ*= 0 are:

_{i}where *G*’(*R, C*_{1},*C*_{2}) ≡ *g*(*R*, *u*’_{1}(*R, C*_{1},*C*_{2}), *u*’_{2}(*R, C*_{1},*C*_{2}),*C*_{1},*C*_{2}).

Here, Eq. S65 corresponds to three unparallel surfaces and share a common point (Fig. 1G and Appendix-fig 3A). However, all the fixed points are unstable (Appendix-fig. 3C), and hence the consumer species cannot stably coexist at steady state (Fig. 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 *ẋ _{i}* = 0 (

*i*= 1,2),

*ż*= 0,

*Ċ*= 0, and note that , we have

_{i}Since *C _{i}* > 0, then

If *R* ≫ *C*_{1} + *C*_{2}, then *R* ≫ *x*_{1} + *x*_{2} and *R*^{(F)} ≈ *R*, we have

Still, we assume that the population dynamics of the resource species follows Eq. S42. At the fixed point, *Ṙ* = 0. We have

Combined with Eq. S68, we can solve for *R*

where and .

Eqs. S68, S70 are the analytical solutions of the fixed point when *R* ≫ *C*_{1} + *C*_{2}. As shown in Appendix-fig. 3B, the analytical predictions agree well with the numerical results (the exact solutions).

# Appendix V Scenario involving chasing pair and both intra- and inter-specific interference

Here we consider the scenario involving chasing pair and both intra- and inter-specific interference in the simple case of *S _{C}* = 2 and

*S*= 1:

_{R}We adopt the same notations as that depicted in Appendix III. A and Appendix IV. Then, and *R* = *R*^{(F)} + *x*_{1} + *x*_{2}, and the population dynamics of the system can be described as follows:

Here, the functional form of *g*(*R*, *x*_{1}, *x*_{2}, *C*_{1}, *C*_{2}) follows Eq. S42. For convenience, we define *K _{i}* ≡ (

*d*+

_{i}*k*)/

_{i}*a*,

_{i}*α*≡

_{i}*D*/(

_{i}*w*),

_{i}k_{i}*β*≡

_{i}*a’*/

_{i}*d’*, and

_{i}*γ*≡

*a*’

_{12}/

*d*’

_{12}, (

*i*= 1,2). At steady state, from

*ẋ*= 0,

_{i}*ẏ*= 0,

_{i}*ż*= 0, and

*Ċ*= 0, (

_{i}*i*= 1, 2), we have

Combined with , and since *C _{i}* > 0(

*i*= 1, 2), then,

## 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)} ≈ *R*. Combined with Eq. S73, we obtain

Using *Ṙ* = 0 and *R* > 0, we have

where and . Eqs. S74–S75 are the analytical solutions of the species abundances at steady state when *R* ≫ *C*_{1} + *C*_{2}. As shown in Appendix-fig. 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 pair 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-fig. 5B), and the fixed points are globally attracting (Appendix-fig. 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-fig. 5A). In the stochastic framework, just as the scenario involving chasing pair and intraspecific interference, the coexistence state can be maintained along with stochasticity (Appendix-fig. 5D).

# Appendix VI Dimensional analysis for the scenario involving chasing pair and both intra- and inter-specific interference

The population dynamics of the system involving chasing pair and both intra- and inter-specific interference are shown in Eqs. 1–4:

with *l* = 1, … , *S _{R}*;

*i,j*= 1, … ,

*S*, and

_{C}*i*≠ =

*j*. Here and represent the population abundances of the consumers and resources in the system. In fact, there are already several dimensionless variables and parameter in Eq. S76, namely

*x*,

_{il}*y*,

_{i}*z*, , ,

_{ij}*C*,

_{i}*R*,

_{l}*w*, . To make all terms dimensionless, we define , where and is a reducible dimensionless parameter which is freely to take any positive values. Besides, we define dimensionless parameters , , , , , , , and . By substituting all the dimensionless terms into Eq. S76, we have

_{il}For convenience, we omit the notation “a” and use dimensionless variables and parameters in the simulation studies unless otherwise specified.

# Appendix VII Simulation details of the main text figures

In Fig. 1C, F: *a _{i}* = 0.1,

*d*= 0.5,

_{i}*w*= 0.1,

_{i}*k*= 0.1 (

_{i}*i*= 1, 2);

*D*

_{1}= 0.002,

*D*

_{2}= 0.001,

*K*

_{0}= 5,

*R*= 0.05. In Fig. 1D, G:

_{a}*a*= 0.02,

_{i}*a’*= 0.021,

_{ij}*d*= 0.5,

_{i}*d’*= 0.01,

_{ij}*w*= 0.08,

_{i}*k*= 0.03,

_{i}*i, j*= 1,2,

*i*≠

*j*,

*D*

_{2}= 0.001,

*D*

_{1}= 0.0011,

*K*

_{0}= 20,

*R*= 0.01. In Fig. 1E, H:

_{a}*a*= 0.5,

_{i}*a’*= 0.525,

_{i}*d*= 0.5,

_{i}*d’*= 0.5,

_{i}*w*= 0.2,

_{i}*k*= 0.4 (

_{i}*i*= 1, 2),

*D*

_{1}= 0.022,

*D*

_{2}= 0.020,

*K*

_{0}= 10,

*R*= 0.1. Fig. 1C, F were calculated or simulated from Eqs. 1, 4. Fig. 1D, G were calculated or simulated from Eqs. 1, 3, 4. Fig. 1E, H were calculated or simulated from Eqs. 1, 2, 4. The analytical solutions in Fig. 1E were calculated from Eqs. S41 and S43.

_{a}In Fig. 2A: *a _{i}* = 0.02,

*a’*= 0.025,

_{i}*d*= 0.7,

_{i}*d’*= 0.7,

_{i}*w*= 0.4,

_{i}*k*= 0.05 (

_{i}*i*= 1, 2);

*D*

_{1}= 0.0160,

*D*

_{2}= 0.0171,

*K*

_{0}= 2000,

*R*= 5.5. In Fig. 2B-C:

_{a}*L*= 100,

*r*

^{(C)}= 5,

*r*

^{(I)}= 5, = 1,

*v*= 0.1,

_{R}*a*= 0.2010,

_{i}*a’*= 0.2828,

_{i}*d’*= 0.8,

_{i}*d*= 0.7,

_{i}*w*= 0.33,

_{i}*k*= 0.2 (

_{i}*i*= 1, 2);

*D*

_{1}= 0.0605,

*D*

_{2}= 0.0600,

*K*

_{0}= 1000,

*R*= 100. In Fig. 2D:

_{a}*a*= 0.3,

_{i}*a’*= 0.33,

_{i}*w*= 0.018,

_{i}*k*= 4.8,

_{i}*d’*= 5,

_{i}*d*= 5.5 (

_{i}*i*= 1, 2);

*D*

_{2}= 0.010,

*R*= 35,

_{a}*K*

_{0}= 10000,

*D*

_{1}= 0.011. In Fig. 2E:

*w*= 0.02,

_{i}*k*= 4.5,

_{i}*d’*= 4,

_{i}*d*= 4.5 (

_{i}*i*= 1, 2);

*D*

_{2}= 0.010,

*R*= 35,

_{a}*K*

_{0}= 10000,

*a*= 0.2,

_{i}*a*= 0.24 (

_{i}*i*= 1, 2);

*D*

_{1}= 0.0120. In Fig. 2D-E:

*τ*= 0.4Day (see Appendix VI). Fig. 2A-E were simulated from Eqs. 1, 2, 4. See Appendix-fig. 7C, E for the long-term time series of all species in Fig. 2D-E, respectively.

Model settings in Fig. 3A-B, D (plankton): *a _{il}* = 0.1,

*a’*= 0.125,

_{i}*d*= 0.5,

_{il}*d’*= 0.2,

_{i}*w*= 0.3,

_{il}*k*= 0.2, = 8 × 10

_{il}^{4}, = 5 × 10

^{4}, = 3 × 10

^{4}, = 280, = 200, = 150,

*D*= 0.03 ×

_{i}*N*(1,0.25) (

*i*= 1, …,

*S*,

_{C}*l*= 1, …,

*S*),

_{R}*S*= 140 and

_{C}*S*= 3. Model settings in Fig. 3C (bird):

_{R}*a*= 0.1,

_{i}*a’*= 0.125,

_{i}*d*= 0.5,

_{i}*d’*= 0.5,

_{i}*w*= 0.3,

_{i}*k*= 0.2,

_{i}*D*= 0.02 ×

_{i}*N*(1, 0.28) (

*i*= 1, …,

*S*);

_{C}*R*= 110,

_{a}*K*

_{0}= 10

^{5},

*S*= 250 and

_{C}*S*= 1. Model settings in Fig. 3C (fish):

_{R}*a*= 0.1,

_{i}*a’*= 0.14,

_{i}*d*= 0.5,

_{i}*d’*= 0.5,

_{i}*w*= 0.2,

_{i}*k*= 0.1,

_{i}*D*= 0.015 ×

_{i}*N*(1, 0.32) (

*i*= 1, …, 45);

*R*= 550,

_{a}*K*

_{0}= 10

^{6},

*S*= 45 and

_{C}*S*= 1. Model settings in Fig. 3C (butterfly):

_{R}*a*= 0.1,

_{i}*a’*= 0.125,

_{i}*d*= 0.5,

_{i}*d’*= 0.3,

_{i}*w*= 0.3,

_{i}*k*= 0.2,

_{i}*D*= 0.034 ×

_{i}*N*(1, 0.35) (

*i*= 1, …,

*S*);

_{C}*R*= 300,

_{a}*K*

_{0}= 10

^{5},

*S*= 150 and

_{C}*S*= 1. Model settings in Fig. 3D (bat):

_{R}*a*= 0.1,

_{i}*a’*= 0.125,

_{i}*d*= 0.5,

_{i}*d’*= 0.5,

_{i}*w*= 0.2,

_{i}*k*= 0.1,

_{i}*D*= 0.013 ×

_{i}*N*(1, 0.34) (

*i*= 1, …,

*S*);

_{C}*R*= 250,

_{a}*K*

_{0}= 10

^{6},

*S*= 40 and

_{C}*S*= 1. Model settings in Fig. 3D (lizard):

_{R}*a*= 0.1,

_{i}*a’*= 0.125,

_{i}*d*= 0.5,

_{i}*d’*= 0.5,

_{i}*w*= 0.2,

_{i}*k*= 0.1,

_{i}*D*= 0.014 ×

_{i}*N*(1, 0.34) (

*i*= 1, …,

*S*);

_{C}*R*= 250,

_{a}*K*

_{0}= 10

^{6},

*S*= 55 and

_{C}*S*= 1. In Fig. 3A-D, the mortality rate

_{R}*D*(

_{i}*i*= 1, …,

*S*) is the only parameter that varies with the consumer species, which was randomly sampled from a Gaussian distribution

_{C}*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-figs. 8–14. Fig. 3A-D were simulated from Eqs. 1, 2, 4. See Appendix-figs. 14C, 10K, C, D, H, I, J, Fig. 3A, Fig. 3B for the time series of Fig. 3C ), 3C (), 3C (), 3D (), 3D (), 3D (), 3D (), 3D () and 3D (), respectively. The Shannon entropies of the experimental data and simulation results for each ecological community are: = 5.67(6.79), = 6.63(6.79), = 4.78(4.12), = 3.78(3.40); = 3.00(2.95, 2.84), = 4.05(3.57, 3.50); = 4.68(6.43, 6.48). Here the Shannon entropy , where

*P*is the probability that a consumer individual belongs to species

_{i}*C*.

_{i}# References

- Coupling in predator-prey dynamics: ratio-dependence
*Journal of Theoretical Biology***139**:311–326https://doi.org/10.1016/S0022-5193(89)80211-5 - Experimental invalidation of the principle of competitive exclusion
*Nature***224**:1076–1079https://doi.org/10.1038/2241076a0 - Mutual interference between parasites or predators and its effect on searching efficiency
*Journal of Animal Ecology***44**:331–340https://doi.org/10.2307/3866 - The evolution of male traits in social insects
*Annual Review of Entomology***50**:395–420https://doi.org/10.1146/annurev.ento.50.071803.130416 - Fighting fruit flies: a model system for the study of aggression
*PNAS***99**:5664–5668https://doi.org/10.1073/pnas.082102599 - Macarthur’s consumer-resource model
*Journal of Theoretical Biology***37**:26–38https://doi.org/10.1016/0040-5809(90)90025-Q - Mechanisms of maintenance of species diversity
*Annual Review of Ecology and Systematics***31**:343–336 - Life after logging: post-logging recovery of a neotropical bat community
*Journal of Applied Ecology***42**:409–420https://doi.org/10.1111/j.1365-2664.2005.01024.x - Long-Term Studies of Vertebrate Communities
- Collective behaviour can stabilize ecosystems
*Nature Ecology & Evolution***5**:1435–1440https://doi.org/10.1038/s41559-021-01517-w - The metagenomics of soil
*Nature Reviews Microbiology***3**:470–478https://doi.org/10.1038/nrmicro1160 - Automated monitoring and analysis of social behavior in drosophila
*Nature Methods***6**:297–303https://doi.org/10.1038/nmeth.1310 - A model for tropic interaction
*Ecology***56**:881–892https://doi.org/10.2307/1936298 - Species diversity in vertical, horizontal, and temporal dimensions of a fruitfeeding butterfly community in an ecuadorian rainforest
*Biological Journal of the Linnean Society***62**:343–364https://doi.org/10.1006/bijl.1997.0155 - A latitudinal diversity gradient in planktonic marine bacteria
*PNAS***105**:7774–7778https://doi.org/10.1073/pnas.0803070105 - The Struggle for Existence
- Stochastic simulation of chemical kinetics
*Annual Review of Physical Chemistry***58**:35–55https://doi.org/10.1146/annurev.physchem.58.032806.104637 - Emergent simplicity in microbial community assembly
*Science***361**:469–474https://doi.org/10.1126/science.aat1168 - Diversity, stability, and reproducibility in stochastically assembled microbial ecosystems
*Physical Review Letters***120**https://doi.org/10.1103/PhysRevLett.120.158102 - Higher-order interactions stabilize dynamics in competitive network models
*Nature***548**:210–213https://doi.org/10.1038/nature23273 - Individual-based Modeling and Ecology
- Effective resource competition model for species coexistence
*Physical Review Letters***127**https://doi.org/10.1103/PhysRevLett.127.208101 - The competitive exclusion principle
*Science***131**:1292–1297https://doi.org/10.1126/science.131.3409.1292 - Bird Community Dynamics in a Temperate Deciduous Forest: Long-Term Trends at Hubbard Brook
*Ecol. Monogr***56**:201–220 - Amazonia through time: Andean uplift, climate change, landscape evolution, and biodiversity
*Science***330**:927–931https://doi.org/10.1126/science.1194585 - Emergent phases of ecological diversity and dynamics mapped in microcosms
*Science***378**:85–89https://doi.org/10.1126/science.abm7841 - The Unifled Neutral Theory of Biodiversity and Biogeography
- A formal derivation of the “beddington” functional responset
*Journal of Theoretical Biology***185**:389–400https://doi.org/10.1006/jtbi.1996.0318 - Biodiversity of plankton by species oscillations and chaos
*Nature***402**:407–410https://doi.org/10.1038/46540 - The paradox of the plankton
*The American Naturalist***95**:137–145https://doi.org/10.1086/282171 - Counteraction of antibiotic production and degradation stabilizes microbial communities
*Nature***521**:516–519https://doi.org/10.1038/nature14485 - Competitive coexistence of two predators utilizing the same prey under constant environmental conditions
*Journal of Theoretical Biology***44**:387–395https://doi.org/10.1016/0022-5193(74)90169-6 - Chemical warfare between microbes promotes biodiversity
*PNAS***99**:786–790https://doi.org/10.1073/pnas.012399899 - Community equilibria and stability, and an extension of the competitive exclusion principle
*The American Naturalist***104**:413–423 - Dispersion and population interactions
*The American Naturalist***108**:207–228https://doi.org/10.1086/282900 - Coexistence in a variable environment
*The American Naturalist***114**:765–783https://doi.org/10.1086/283527 - Species packing and competitive equilibrium for many species
*Journal of Theoretical Biology***1**:1–11https://doi.org/10.1016/0040-5809(70)90039-0 - Competition, habitat selection, and character displacement in a patchy environment
*PNAS***51**:1207–1210https://doi.org/10.1073/pnas.51.6.1207 - Ocean currents promote rare species diversity in protists
*Science Advances***6**https://doi.org/10.1126/sciadv.aaz9037 - The structure and organisation of an Amazonian bird community remains little changed after nearly four decades in Manu National Park
*Ecology Letters***26**:335–346https://doi.org/10.1111/ele.14159 - Some mathematical problems concerning the ecological principle of competitive exclusion
*Journal of Differential Equations***23**:30–52https://doi.org/10.1016/0022-0396(77)90135-8 - Mortality rates of autochthonous and fecal bacteria in natural aquatic ecosystems
*Water Research***37**:4151–4158https://doi.org/10.1016/S0043-1354(03)00349-X - Microbial coexistence through chemical-mediated interactions
*Nature Communications***10**:2052–474https://doi.org/10.1038/s41467-019-10062-x - Experimental studies of interspecies competition ii. temperature, humidity, and competition in two species of tribolium
*Physiological Zoology***27**:177–238https://doi.org/10.1086/physzool.27.3.30152164 - What determines species diversity?
*Science***309**:90–90https://doi.org/10.1126/science.309.5731.9 - Metabolic trade-offs promote diversity in a model ecosystem
*Physical Review Letters***118**https://doi.org/10.1103/PhysRevLett.118.028103 - Strength of species interactions determines biodiversity and stability in microbial communities
*Nature Ecology & Evolution***4**:376–383https://doi.org/10.1038/s41559-020-1099-4 - Ubiquitous abundance distribution of non-dominant plankton across the global ocean
*Nature Ecology & Evolution***2**:1243–1249https://doi.org/10.1038/s41559-018-0587-2 - Tara oceans: towards global ocean ecosystems biology
*Nature Reviews Microbiology***18**:428–445https://doi.org/10.1038/s41579-020-0364-5 - Structure and organization of an amazonian forest bird community
*Ecological Monographs***60**:213–238https://doi.org/10.2307/1943045 - Elements ofa theory for the mechanisms controlling abundance, diversity, and biogeochemical role of lytic bacterial viruses in aquatic systems
*Limnology and Oceanography***45**:1320–1328https://doi.org/10.4319/lo.2000.45.6.1320 - Eukaryotic plankton diversity in the sunlit ocean
*Science***348**https://doi.org/10.1126/science.1261605 - Diverse modes of eco-evolutionary dynamics in communities of antibiotic producing microorganisms
*Nature Ecology & Evolution***1**https://doi.org/10.1038/s41559-017-0189 - Overcome competitive exclusion in ecosystems
*iScience***23**https://doi.org/10.1016/j.isci.2020.101009 - Growth strategy of microbes on mixed carbon sources
*Nature Communications***10**https://doi.org/10.1038/s41467-019-09261-3 - Spatial ecology of territorial populations
*PNAS***116**:17874–17879https://doi.org/10.1073/pnas.1911570116 - Coevolution maintains diversity in the stochastic “kill the winner” model
*Physical Review Letters***119**https://doi.org/10.1103/PhysRevLett.119.268101

# Article and author information

## Version history

- Preprint posted:
- Sent for peer review:
- Reviewed Preprint version 1:
- Reviewed Preprint version 2:

## Copyright

© 2024, Kang et al.

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

- views
- 313
- downloads
- 16
- citations
- 0

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