A model of pairwise encounters may naturally break CEP.

(A) Ageneric model of SC consumer species feeding on SR resource species with pairwise encounters. (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: Ṙ = 0 (orange surface), Ċ1 = 0 (blue surface), Ċ2 = 0 (green surface). 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 VII for simulation details.

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 two years with one type of food (flour) within a lab (Park, 1954).

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; Devries et al., 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 × 105 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: , , , ; , , , ; , . See Appendix VII for simulation details and the Shannon entropies.

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

To calculate ail in the chasing pair, we suppose that all individuals of species Rl stand still while a Ci individual moves at the speed of (the relative speed). Over a time interval of Δt, the length of zigzag trajectory of the Ci individual is approximately , while the encounter area (marked with dashed lines) is estimated to be . Then, we can estimate the encounter rate ail using the encounter area and concentrations of the species (see Materials and Methods for details). Similarly, we can estimate a’ij in the interference pair.

Functional response in scenarios involving different types of pairwise encounter.

(A-B) In the scenario involving only chasing pair, the red surface corresponds to the B-D model (calculated from Eq.S10), while the green surface represents the exact solutions to our mechanistic model (calculated from Eq. S7). The magenta (calculated from Eq. S8) and blue (calculated from Eq. S9) surfaces represent the approximate solutions to our model. (C-D) In the scenario involving chasing pair and intraspecific interference, the red surface corresponds to the B-D model (calculated from Eq.S24), while the green surface represents the exact solutions to our mechanistic model (calculated from Eq. S17). The blue surface (calculated from Eq. S19) and the magenta surface (calculated from Eq. S21) represent the quasi-rigorous and the approximate solutions to our model, respectively. (E-F) In the scenario involving chasing pair and interspecific interference, the red surface corresponds to the B-D model (calculated from Eq. S31), while the green surface represents the quasi-rigorous solutions to our mechanistic model (calculated from Eq. S27). The blue surface (calculated from Eq. S29) represents the approximate solutions to our model. In (A-B): k = 0.1, a = 0.25. In (A): d = 0. (see Appendix II B). In (C-D): a = 0.1, k = 0.1, d’ = 0.1, a’ = 0.12. In (C): d = 0 (see Appendix II C). In (E-F): a1 = a2 = 0.1, k1 = k2 = 0.1, d’12 = 0.1, a’12 = 0.12. In (E): di = 0 (i = 1,2) (see Appendix II D).

Fixed point solutions in scenarios involving chasing pair and different types of predator interference.

Here, SC = 2 and SR = 1. Di (i = 1, 2) is the only parameter varying with the consumer species (with D1 > D2), and Δ = (D1 − D2)/D2 represents the competitive difference between the two species. (A-C) Scenario involving chasing pair and interspecific interference. (D-F) Scenario involving chasing pair and intraspecific interference. (A, D) Positive solutions to the steady-state equations: Ṙ = 0 (orange surface), Ċ1 = 0 (blue surface), Ċ2 = 0 (green surface). The intersection point marked by black/red dots is an unstable/stable fixed point. (B, E) 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 (B) and (E) (marked with superscript “(A)”) were calculated from Eqs. S68, S70 and Eqs. S41, S43, respectively. (C) In this scenario, there is no parameter region for stable coexistence. The region below the red surface and above Δ = 0 represents unstable fixed points. (F) Comparisons between the numerical results and analytical solutions of the coexistence region. Here represents the maximum competitive difference tolerated for species coexistence. The red and cyan surfaces represent the analytical solutions (calculated from Eq. S46) and numerical results, respectively. The numerical results in (A-C) and (D-F) were calculated from Eqs. S42, S61 and Eqs. S33, S42, respectively. In (A): a1 = a2 = 0.05, d1 = d2 = 0.05, K0 = 20, a’12 = 0.06, d’12 = 0.01, k1 = k2 = 0.02, w1 = w2 = 0.08, D1 = 0.0011, D2 = 0.001, Ra = 0.01. In (B): a1 = a2 = 0.04, d1 = d2 = 0.2, K0 = 10, a’12 = 0.048, d’12 = 0.001, k1 = k2 = 0.1, w1 = w2 = 0.3, D2 = 0.0008, Ra = 0.2. In (C): a1 = a2 = 0.1, K0 = 100, k1 = k2 = 0.1, w1 = w2 = 0.1, D2 = 0.001, a’12 = 0.12, Ra = 0.05. In (D): a1 = a2 = 0.5, a’1 = a’2 = 0.625, d1 = d2 = 0.5, d’1 = d’2 = 0.5, k1 = k2 = 0.4, D2 = 0.02, w1 = w2 = 0.5, K0 = 10, D1 = 1.2D2, Ra = 0.1. In (E): a1 = a2 = 0.1, a’1 = a’2 = 0.12, k1 = k2 = 0.12, w1 = w2 = 0.3, D2 = 0.02, K0 = 100, Ra = 0.8, d1 = d2 = 0.5, d’1 = d’2 = 0.05. In (F): a1 = a2 = 0.5, k1 = k2 = 0.2, w1 = w2 = 0.2, D2 = 0.008, K0 = 60, Ra = 0.8, d1 = d2 = 0.8.

Intraspecific predator interference facilitates species coexistence regardless of stochasticity.

Here we consider the case of SC = 2, SR = 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) A 3D phase diagram in the ODEs studies. Di (i = 1, 2) is the only parameter that varies with the two consumer species, and Δ ≡ (D1 — D2)/D2 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). (C) Time courses of the species abundances simulated with ODEs or SSA. (D) Representative trajectories of species coexistence in the phase space simulated with SSA. The coexistence state is stable and globally attractive (see (C) for time courses, SSA results). (A-D) were calculated or simulated from Eqs. 1, 2, 4. In (A): ai = 0.1, a’i = 0.125, di = 0.1, d’i = 0.05, wi = 0.1, ki = 0.1 (i = 1, 2); D1 = 0.0035, D2 = 0.0038, K0 = 100, Ra = 0.3. In (B): ai = 0.1, di = 0.1, wi = 0.1, ki = 0.1 (i = 1, 2); K0 = 100, D2 = 0.001, Δ = (D1 — D2)/D2, K0 = 100, Ra = 0.1. In (C-D): ai = 0.02, a’i = 0.025, di = 0.7, d’i = 0.7, wi = 0.4, ki = 0.05 (i = 1, 2); D1 = 0.0160, D2 = 0.0171, K0 = 2000, Ra = 5.5.

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

Here, Di (i = 1, 2) is the only parameter that varies with the consumer species (with D1 > D2), and Δ = (D1 − D2)/D2 measures the competitive difference between the two species. (A-C, E) ODEs results. (D) ODEs and SSA results. (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 tra jectories 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 Eqs. S74–S75. The numerical results in (A-E) were calculated or simulated from Eqs. 1–4. In (A-C): ai = a2 = 0.1, a’1 = a’2 = 0.12, ki = k2 = 0.1, wi = w2 = 0.1, D2 = 0.004, K0 = 100, a’12 = 0.12, Ra = 0.8, d’12 = 0.5. In (D): a1 = a2 = 0.1, a’1 = a’2 = 0.14, k1 = k2 = 0.12, w1 = w2 = 0.15, D1 = 0.0125, D2 = 0.012, K0 = 300, Ra = 5.5, d1 = d2 = 0.3, d’1 = d’2 = 0.5, a’12 = 0.14, d’12 = 5. In (E): a1 = a2 = 0.05, a’1 = a’2 = 0.06, k1 = k2 = 0.1, w1 = w2 = 0.2, D2 = 0.008, K0 = 100, Ra = 0.8, d1 = d2 = 0.5, d’1 = d’2 = 0.15, a’12 = 0.06, d’12 = 0.0005.

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 Ref. (Koch, 1974) and (Huisman and Weissing, 1999) for the parameters. (C-D) Phase diagrams in the scenario involving chasing pair and intraspecific interference. Here, SC = 2 and SR = 1. Di (i = 1, 2) is the only parameter varying with the consumer species (with D1 > D2), and Δ = (D1 — D2)/D2 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 Eqs. 1, 2, 4. In (C-D): ai = 0.1, a’i = 0.125, di = 0.5, wi = 0.1, ki = 0.1 (i = 1, 2); K0 = 100, Ra = 5, D2 = 0.0014.

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

(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. In all cases, two consumer species coexist for long with only one type of resources. (A) In Ayala’s experiment (Ayala, 1969), two Drosophila species (consumers), D. serrata and D. pseudoobscura, compete for the same type of abiotic resources within a laboratory bottle. (B) In Park’s experiment (Park, 1954), two Tribolium species, T. confusum and T. castaneum, compete for the same food (flour). (C-F) Time courses of the species abundances in the scenario involving chasing pair and intraspecific interference. The time series in (C-F) correspond to the long-term version of that shown in Fig. 2D, Appendix-fig. 7A, Fig. 2E, Appendix-fig. 7B, respectively. (A-F) were simulated from Eqs. 1, 2, 4. In (A): ai = 0.3, a’i = 0.33, wi = 0.018, ki = 4.8, d’i = 5, di = 5.5 (i = 1, 2); D2 = 0.010, Ra = 35, K0 = 10000, D1 = 0.0132. In (B): wi = 0.02, ki = 4.5, d’i = 4, di = 4.5 (i = 1, 2); D2 = 0.010, Ra = 35, K0 = 10000, ai = 0.3, a’i = 0.36 (i = 1, 2); D1 = 0.0122. In(A-B): τ = 0.4Day (see Appendix VI).

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

(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 × 105 in the time series (see (A) and (B)), respectively. The Shannon entropies of the experimental data and simulation results for the plankton community are: = 2.85(2.18, 2.00). In the model settings, SC = 200 and SR = 1. Di (i = 1, …, SC) 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 Eqs. 1,2,4. In (A-C): ai = 0.1, a’i = 0.125, d’i = 0.5, di = 0.2, wi = 0.2, ki = 0.1, K0 = 105 , Di = N(1, 0.38) × 0.008 (i = 1, …, 200), Ra = 150.

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

The solid dots represent the experimental data (marked with “Exp”) reported in existing studies (Hubbell, 2001), (Holmes et al., 1986), (Cody et al., 1996), while the hollow dots and those with “+” center are the ODEs and SSA results constructed from timestamp t = 1.0 × 105 in the time series (see Appendix-fig. 10A-G), respectively. In the model settings, SR = 1, SC = 20 (in (A)), 35 (in (B)) or 45 (in (C)). Di (i = 1, …, SC) 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: = 2.98(3.06, 2.98), = 4.04(4.02, 4.02), = 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: = 0.89, = 0.88; = 0.47, = 0.75; = 0.88, = 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 Eqs. 1, 2, 4. In (A): ai = 0.1, a’i = 0.125, di = 0.5, d’i = 0.6, wi = 0.22, ki = 0.1, Di = 0.016 × N(1, 0.35) (i = 1, …, 20); Ra = 350, K0 = 106. In (B): ai = 0.1, a’i = 0.125, di = 0.5, d’i = 0.6, wi = 0.22, ki = 0.1, Di = 0.012 × N(1, 0.35) (i = 1, …, 35); Ra = 350, K0 = 106. In (C): ai = 0.1, a’i = 0.14, di = 0.5, d’i = 0.5, wi = 0.2, ki = 0.1, Di = 0.015 × N(1, 0.32) (i = 1, …, 45); Ra = 550, K0 = 106.

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

The time series in (A, E), (B, F), (C, G), (D, H), (I, J) and (K) correspond to that shown in Appendixfig. 9A-C, Fig. 3D (bat), Fig. 3D (lizard) and Fig. 3C (butterfly), respectively.

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

The solid dots represent the experimental data (marked with “Exp”) reported in existing studies (Hubbell, 2001), (Holmes et al., 1986), (Cody et al., 1996), (Clarke et al., 2005), while the hollow dots and those with “+” center are the ODEs and SSA results constructed from timestamp t = 1.0 × 105 in the time series (see Appendix-fig. 13), respectively. In the model settings, SR = 3, SC = 20 (in (A)), 35 (in (B)), 40 (in (C)), 45 (in (D)) or 50 (in (E)). Di (i = 1, …, SC) 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: = 2.98(2.98, 3.26), = 4.04(4.34, 4.35), = 3.00(3.00, 3.00), = 3.78(3.28, 3.64), = 4.05(3.94, 3.94). In the K-S test, the p values that the simulation results and the corresponding experimental data come from identical distributions are: = 0.59, = 0.43, = 0.47, = 0.33, = 0.42, = 0.27, = 0.22, = 0.06, = 0.56, = 0.36. With a significance threshold of 0.05, none of the p values suggest there exists a statistically significant difference. The numerical results in (A-E) were simulated from Eqs. 1, 2, 4. In (A-E): ail = 0.1, a’i = 0.125, dil = 0.5. In (a): d’i = 0.3, wil = 0.2, kil = 0.12, , , , Di = 0.021 × N(1, 0.28) (i = 1, …, 20, l = 1, 2, 3); , , . In (B): d’i = 0.6, wil = 0.2, kil = 0.12, = 8 × 104, , , Di = 0.017 × N(1, 0.3) (i = 1, …, 35, l = 1,2, 3); , = 160, . In (C): d’i = 0.4, wil = 0.3, kil = 0.12, , , , Di = 0.023 × N(1, 0.34) (i = 1, …, 40, l = 1, 2, 3); , , . In (d): d’i = 0.3, wil = 0.3, kil = 0.12, , = 5 × 104, , Di = 0.027 × N(1, 0.32) (i = 1, …, 45, l = 1, 2, 3); = 80, , . In (E): d’i = 0.3, wil = 0.3, kil = 0.2, , , , Di = 0.034 × N(1, 0.34) (i = 1, …, 50, l = 1, 2, 3); , , .

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

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 × 105 in the time series (see Appendix-fig. 13), respectively. The plankton community data were obtained separately from the Norwegian Sea (NS) and Pacific Station (PS). In the model settings, SR = 1 (in (B-C)), 3 (in (A)); SC = 50 (in (A, C)), 150 (in (B)). Di (i = 1, …, SC) 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: for SR = 3; = 4.68(6.53, 6.16), for SR. In the K-S test, the p values that the simulation results and the corresponding experimental data come from identical distributions are: , for SR = 3; ; , , for SR = 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 Eqs. 1, 2, 4. In (A): ail = 0.1, a’i = 0.125, dil = 0.5, d’i = 0.2, wil = 0.3, kil = 0.2, = 8 × 104, , , = 280, , , Di = 0.035 × N(1, 0.25) (i = 1, …,50, l = 1, 2, 3). In (b): ai = 0.1, a’i = 0.125, di = 0.3, d’i = 0.3, wi = 0.3, ki = 0.2, Di = 0.025 × N(1, 0.25) (i = 1, …, 150, l = 1, 2, 3); Ra = 350, K0 = 104. In (c): ai = 0.1, a’i = 0.125, di = 0.3, d’i = 0.3, wi = 0.3, ki = 0.2, Ra = 350, K0 = 104, Di = 0.03 × N(1, 0.27) (i = 1, …, SC).

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

The time series in (A, E), (B, F), (C, G), (D, H), (I, J), (K, L), (M, N) and (O, P) correspond to that shown in Appendix-figs. 11A-E and 12A-C, respectively.

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

(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 × 105 in the time series (see (B)). (B) Time courses of the species abundances simulated with ODEs. In the model settings, SR = 3 and SC = 250. Di (i = 1, …, SC) 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 . In the K-S test, the p values that the simulation results and the corresponding experimental data come from identical distributions are: , . 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 Fig. 3A (bird), and the simulation parameters are the same as Fig. 3A (bird). The numerical results in (A-C) were simulated from Eqs. 1, 2, 4. In (A-B): ail = 0.1, a’i = 0.125, dil = 0.5, d’i = 0.6, wil = 0.3, kil = 0.2, Di = 0.032 × N(1, 0.17) (i = 1, …,250, l = 1, 2, 3); , , , , , .

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, SC = 40 and SR = 1. xi represents , and yi represents . Xi/Ci and yi/Ci stand for the fractions of consumer individuals within a chasing pair and an interference pair, respectively. The numerical results were calculated from Eq. S55, while the analytical solutions (marked with superscript “(A)”) were calculated from Eq. 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): ai = 0.0015, a’i = 0.0021, di = 0.1, d’i = 0.05, ki = 5. In (B): R = 2000.