Positive density dependence acting on mortality can help maintain species-rich communities

Conspecific negative density dependence is ubiquitous and has long been recognized as an important factor favoring the coexistence of competing species at local scale. By contrast, a positive density-dependent growth rate is thought to favor species exclusion by inhibiting the growth of less competitive species. Yet, such conspecific positive density dependence often reduces extrinsic mortality (e.g. reduced predation), which favors species exclusion in the first place. Here, using a combination of analytical derivations and numerical simulations, I show that this form of positive density dependence can favor the existence of equilibrium points characterized by species coexistence. Those equilibria are not globally stable, but allow the maintenance of species-rich communities in multispecies simulations. Therefore, conspecific positive density dependence does not necessarily favor species exclusion. On the contrary, some forms of conspecific positive density dependence may even help maintain species richness in natural communities. These results should stimulate further investigations into the precise mechanisms underlying density dependence.


Differential equations
The systems of differential equations describing the changes in density of the two species (N 1 and N 2 ) are: The two species have the same carrying capacity K, the same growth rate r, the same mortality rate m and the same density-dependence factor S.
Following Williams and Banyikwa (1981), these equations are equivalent to: which is the usual form of Lotka-Volterra model. This form of the systems make clear that the term of conspecific positive density dependence changes the effective carrying capacities, K effi = K 1 − m r (1+S Ni) (as discussed by Holt 1985, in the case of density-independent mortality).

Scaling
All variables are rescaled to growth rate and carrying capacity, as follows: n 1 = N 1 /K, rescaled density of species 1; n 2 = N 2 /K, rescaled density of species 2; d = m/r, rescaled mortality; s = S K, rescaled density-dependence factor; t = τ r, rescaled time unit. The rescaled variables are dimensionless. From Equation 1 above, the dynamical system becomes:        d n 1 d t = n 1 1 − n 1 − d 1 + s n 1 d n 2 d t = n 2 1 − n 2 − α 1 n 1 − d 1 + s n 2 (3) This is the system of equations analyzed (Equations 1a and 1b in the main text). Equations 2 and 3 in the main text are scaled in a similar manner.
Supplementary file 1B: Analytical derivation of the two-species model with asymmetric competition for resources I analyze the two-species model with positive density-dependent mortality and asymmetric competition for resources. Using basic linear algebra, I find equilibria and determine their stability. I also determine the threshold value of s (s + ) above which the coexistence equilibrium can be recovered, and I prove that the least competitive species is always less abundant than the most competitive species at coexistence equilibrium. Finally, I show that the least competive species is at higher density at coexistence equilibrium if s → ∞ than if s = 0.
Here, I assume that 0 < d < 1 and 0 < α 1 < 1. Note however that the system is characterized by more unstable equilibria for d > 1; even without competitors, species may not invade if they are too rare due to high mortality at low density.

Calculation:
Given that 1 − n * 1 − d 1+s n * 1 = 0 (equilibrium condition), then d (1+s n * 1 ) 2 = (1−n * 1 ) 2 d and therefore: Function f is a concave polynomial function, which maximum is negative: Therefore, function f is negative and: 1 − 2 n * 1 − d (1+s n * 1 ) 2 < 0 4 Conditions of existence of the coexistence equilibrium n o 1 The conditions of existence of the coexistence equilibrium n o 1 are: Calculation: n * 2 ≥ 0 if the product of the roots of the polynomial is negative or the sum of the roots is positive (n * 2 is the highest root at coexistence equilibrium n o 1). Using Vieta's formulas, this is the case if (see preliminary calculations) , then: With X = 1 − α 1 n * 1 , this is equivalent to: This polynomial function has only one positive root which is equivalent to: with Given that s d < s, s which is a positive polynomial function. Therefore: A(s) > 0.
5 Conditions of existence of the coexistence equilibrium n o 2 The conditions of existence of the coexistence equilibrium n o 2 are: Calculation: However, n * 2 ≥ 0 if the product of the roots of the polynomial is positive and the sum of the roots is positive (n * 2 is the lowest root at coexistence equilibrium n o 2). Using Vieta's formulas, this is the case if with 6 Local stability analysis Jacobian matrix Complete Extinction: unstable equilibrium First eigenvalue: α1+d−1 and α 1 + d − 1 > 0 (see preliminary calculations) Therefore, for s ≤ 1/d, the equilibrium 'only species 1' is unstable if the coexistence equilibrium exists. For s > 1/d, the equilibrium 'only species 1' may be stable even if the coexistence equilibrium exists; this occur when the coexistence equilibrium n o 2 exists.
7 From s = 0 to s = 1/d: loss of the coexistence equilibrium n o 1 if high α 1 For s ≤ 1/d, there is no coexistence equilibrium n o 1 if: There is a coexistence equilibrium n o 1 if s = 0. Nonetheless, if s ≤ 1/d, increasing s can lead to the loss of the coexistence equilibrium n o 1. This occurs for strong asymetric competition. More precisely, increasing s to a value of 1/d can lead to extinction of the least competitive species if: Calculation: For s = 1/d, extinction of the least competitive species occurs if two conditions are fulfilled. The first condition on α 1 is: and the second condition on α 1 must fulfill s > α1 (α1−1) 1−α1−d at s = 1/d: Given that 1 − α 1 − d < 0, extinction occurs at s = 1/d if: and therefore the second condition on α 1 is: We can easily show the second condition on α 1 necessarily fulfill the first condition: Therefore, for s = 1/d, extinction of the least competitive species occurs if: 8 From s = 1/d to +∞: loss and recovery of the coexistence equilibrium n o 1 For s ≥ 1/d, there is a coexistence equilibrium n o 1 if: If function α 1 is decreasing, the coexistence equilibrium n o 1 may vanish with increased s (the condition necessary for coexistence gets more restricted). On the contrary, if function α 1 is increasing, the coexistence equilibrium n o 1 may be recovered with increased s.
More precisely, for s > (threshold value represented by a dashed red line in Figure 4 in the main text), the coexistence equilibrium n o 1 can be recovered with increased s.
Similarly, the sign of A (s) for s > s + is the same than the sign of lim s→∞ A (s): We can conclude that for s > s + (represented by a dashed red line in Figure 4 in the main text), A (s) > 0, i.e., the coexistence equilibrium n o 1 can be recovered with increased s.
Therefore: n * 2 < n * 1 10 At coexistence equilibrium n o 1, species 2 reaches a higher density if s → ∞ For s = 0: and for s → ∞: Supplementary file 1C: Numerical method for analyzing the coexistence equilibrium in two-species models Identifying the coexistence equilibrium I run numerical simulations in Python (Release 2.7.15, 2018) on the interval of time [0,5,000]. For each combination of parameters, I run 121 simulations with different initial species abundances (n 0 1 , n 0 2 ), with n 0 1 and n 0 2 having values belonging to the discrete set { 0, 0.1, 0.2, ... 1.0 }. At the end of each simulation, the equilibrium reached is identified as a 'coexistence equilibrium' if species abundances are higher than 10 −3 . Two coexistence equilibrium points are considered as equivalent if the Euclidean distance between the points is lower than 10 −3 .
For all combinations of parameters tested, only one coexistence equilibrium is identified at most.

Identifying the nature of the coexistence equilibrium
The coexistence equilibrium (n * 1 , n * 2 ) identified numerically is locally stable. Otherwise, it would not have been reached at the end of the simulation.
To identify if the coexistence equilibrium equilibrium is a global attractor, I run a simulation on the interval of time [0, 500] with initial conditions (n * 1 , 10 −7 ), i.e., the least competitive species is very rare. At the end of the simulation, the coexistence equilibrium is considered as a global attractor if n 2 > 10 −7 ; the least competitive species is invading even if it is rare initially. Otherwise, the coexistence equilibrium is considered as a local attractor ; the least competitive species does not invade and the coexistence equilibrium is not reached. If the least competitive species can invade when it is rare, so does the most competitive species. Therefore, the criterion of invasion of the least competitive species is sufficient to identify the coexistence equilibrium as a global attractor.

Validity of the numerical method
To test the validity of this numerical method, I applied it to the two-species model with asymmetric competition for resources and conspecific positive density dependence acting only on mortality. I compared the output with analytical derivations. As shown in Figure 4-Figure supplement 2, the numerical method predicts very well the condition of existence of the coexistence equilibrium and its nature (local or global attractor). Therefore, we can confidently use this numerical method to analyze other two-species models that we cannot analyze analytically.

Supplementary file 1D:
Simulations of the multi-species models

Sampling initial communities
For each species i, α i,i and α i,i are set to one. Parameters α j,i and α j,i < 1 (∀j = i) therefore represents the strength of competition for resources and interspecific reproductive interference, relative to intraspecific interactions. Four different scenarios are considered successively.
(1) Only asymmetric competition. -All parameters linked to interspecific reproductive interference (α i,j with i = j) and to differences in basal mortality rates (δ i ) are set to zero. For each pair of species (i, j), α i,j and α j,i are drawn randomly from a uniform distribution within the range [0, 0.1].
(2) Only asymmetric reproductive interference. -I assume there is symmetric competition for resources (α i,j = α if i = j). All parameters linked to differences in basal mortality rates (δ i ) are set to zero. For each pair of species (i, j), α i,j and α j,i are drawn randomly from a uniform distribution within the range [0, 0.1].
(3) Only differences in basal mortality. -I assume there is symmetric competition for resources (α i,j = α if i = j). All parameters linked to interspecific reproductive interference (α i,j with i = j) are set to zero. For each species i, δ i is drawn randomly from a uniform distribution within the range [0, 0.5].
(4) Random communities with all forms of asymmetry. -Species asymmetries act synergistically to promote species exclusion. Therefore, for clarity purpose, the maximum value for each parameter is reduced by half. For each pair of species (i, j), α i,j , α j,i , α i,j and α j,i are drawn randomly from a uniform distribution within the range [0, 0.05]. Likewise, δ i is drawn randomly from a uniform distribution within the range [0, 0.25].
In supplementary analyses, parameters are drawn randomly from uniform distributions within other ranges; the results are qualitatively similar to those obtained with the arbitrary ranges chosen in the main analysis (

Numerical simulations
For each scenario and for each combination of parameters, 500 species pools of 100 species each were constructed by drawing species parameter values from appropriate probability distributions. Numerical simulations are run in Python (Release 2.7.15, 2018) on the interval of time [0, 5000]. All species were equally abundant initially, and simulations were run long enough for initial transients to dissipate. Species were declared extinct if their density fell below 10 −3 . At the end of each simulation, the number of species remaining was recorded.

Community assembly simulations
Following Gross (2008), I repeated simulations in Figure 7 of the main text with the modification that species were introduced one at a time instead of all at once. At every t = 200 time units, one of the species in the species pool was selected at random and introduced with an abundance of 10 −2 or 10 −5 depending on the simulation. After 100 such introductions, the simulation ran to time t = 5, 000. Species were declared extinct if their density fell below 10 −3 . Figure 7-Figure supplement 9, this second set of simulations in which species were introduced one at a time produced different results. Species richness is lower because less competitive species (in the broad sense) may not invade the population. In particular, conspecific positive density dependence does not help produce species-rich communities when species are introduced at low density.