Multi-species coexistence in Lotka-Volterra competitive systems with crowding effects

Classical Lotka-Volterra (LV) competition equation has shown that coexistence of competitive species is only possible when intraspecific competition is stronger than interspecific competition, i.e., the species inhibit their own growth more than the growth of the other species. Note that density effect is assumed to be linear in a classical LV equation. In contrast, in wild populations we can observed that mortality rate often increases when population density is very high, known as crowding effects. Under this perspective, the aggregation models of competitive species have been developed, adding the additional reduction in growth rates at high population densities. This study shows that the coexistence of a few species is promoted. However, an unsolved question is the coexistence of many competitive species often observed in natural communities. Here, we build an LV competition equation with a nonlinear crowding effect. Our results show that under a weak crowding effect, stable coexistence of many species becomes plausible, unlike the previous aggregation model. An analysis indicates that increased mortality rate under high density works as elevated intraspecific competition leading to the coexistence. This may be another mechanism for the coexistence of many competitive species leading high species diversity in nature.

usually treated as constant (i.e., linear), and crowding effect is not included. To consider the nonlinear density effects, aggregation models have been developed and studied extensively introducing 'mean crowding' [8][9][10]15,16 . These models show the coexistence of a few species, but not many species. The 'mean crowding' is a statistical feature of crowding affecting population growth rate (both birth rate and mortality). Here, we develop a simple LV type competition model with crowding effect on mortality only. We assume that the mortality rate of an individual increases with the density of a population. Moreover, using our modified LV competition model with nonlinear crowding effect, we show that multiple species are generally possible to coexist using LV system with crowding effect. This coexistence dynamics should be applicable to insect or animal species.

Results
We consider the LV competition system where all interspecific competition rates are unity, i.e., α ij = α ji = 1∀i. Setting dx i /dt = 0, we obtain the zero isoclines for the modified Lotka-Volterra (LV) competition equations with nonlinear crowding effect rate m i . These isoclines are straight lines in the classical LV competition equation when we sketch the graph of the population density 1 with respect to the population density 2 (Fig. 1, dotted lines). However, we observed that these isoclines turn out to be curved when we include a nonlinear crowding effect m i (Fig. 1, solid lines). All four cases of Lotka-Volterra model show convergent-stable coexistence by adding a crowding effect, where an inferior (superior) species always increases (decreases) in densities (Fig. 1). The equilibrium point is moving from E 1 (d i = 0.0) to E 2 (d i = 0.5) which is also caused by the inclusion of a crowding effect m i , irrespective of other constant parameters (Fig. 1).
Curved lines in the − x x 1 2 plane imply that coexistence only take place when the value of δ → 1.0 for all d i > 0 (Fig. 2). However, by increasing δ, the competitive exclusion reappears in all exclusion cases, where the equilibrium state in − x x 1 2 plane is strongly curved to return to the originated axis ( Fig. 2 except 2c). In contrast, an isocline is straight in − x x 1 2 plane, where the two species coexist without crowding effects, such that < α K K 2 1 12 and < α K K 2 2 12 (Fig. 2c). Furthermore, Fig. 3a shows that two competing species can coexist for any small positive value of δ for all d i > 0. Moreover, small positive real number δ will make coexistence possible but bigger value indicates the competitive exclusion principle again (Fig. 2 except 2c and Fig. 3). Figure 4 shows many-species Lotka-Volterra competition dynamics with crowding effects. Temporal dynamics of all species becomes convergent stable by adding crowding effects (Fig. 4a,b). The effect of δ in 5 or 10 species (Fig. 4c,d) is qualitatively same with the case of two species competition (Fig. 3b). Thus, many-species LV competition model with crowding effects leads to the stable coexistence of all species (Fig. 4).
We also build LV competition models with aggregation effects on mortality that are qualitatively equivalent to the aggregation model of Hartley and Shorrocks 8 . We compare them with the current crowding models using .
Colors indicate x 1 (red) and x 2 (blue). K i : carrying capacity of species i, and α ij : competition coefficient from species j to species i. Parameters value . Colors indicate x 1 (red) and x 2 (blue). Parameter value and software used: see Fig. 1. the same parameter conditions (Fig. 5). In the 2-and 5-species dynamics, all species survive and converge to a stable equilibrium in both the crowding model (Fig. 5a,c) and aggregation model (Fig. 5b,d). In the 10-species dynamics, all species survive and converge to stable equilibrium in the crowding models ( Fig. 5e), while only five species in the aggregation model. We also consider the effect of nonlinear competition terms 21 with or without crowding effect. We compare them with the current crowding models using the same parameter conditions (Fig. 6). In their original model, Taylor and Crizer consider a modified Lotka-Volterra model introducing the effect of nonlinear competition on growth rate. Here, we introduce the nonlinear competition on birth rate alone, since crowding effect is introduced only on mortality rate. In this manner, we can compare the effect of these changes separately. In the 2-species dynamics, both species survive and converge to a stable equilibrium in the crowding model with linear or nonlinear competition terms (Fig. 6a,c). In the current parameter conditions, the LV competition model with nonlinear competition terms lead to the extinction of one species (Fig. 6b). Note that by changing the conditions, this model lead to the coexistence of the two species 21 . In the 5-species dynamics, all five species survive and converge to a stable equilibrium (Fig. 6d,e and f). Interestingly, the model with nonlinear competition terms converges to the same density for all five species (Fig. 6e), while the crowding effect lead to different densities among all species (Fig. 6d,f). In the 10-species dynamics, seven species survive and converge to stable equilibrium in the crowding models (Fig. 6g). In contrast, all ten species in the LV competition model with nonlinear terms lead to the coexistence of all species with an equal density (Fig. 6h). By combining the nonlinear competition and the crowding effect, all ten species survive and converge to different densities (Fig. 6i).

Discussion
Many species compete for precisely the same limited resources to survive 1,2 . Gause's exclusion principle show that multiple competing species cannot coexist in natural communities 1,4,18,19 . Only one species, the superior competitor, will survive and other competitors will eventually become extinct. We should note that frequency dependence does not promote the coexistence of multiple species 20 . Niche theory suggested that it will only become possible for the competing species to coexist if they have different niche 2,3 . Linear density effects show that coexistence becomes possible under very limited conditions. Hence, we search for mechanisms that will enable coexistence of competitive species 5,6 . As a universal and more biologically founded solution, we consider crowding effect, nonlinear density effects at high densities, in the LV competition systems. Classical LV competition model shows that it is only possible for two species to coexist together if intraspecific competition is stronger than interspecific competition 3 . However, our results have shown that inclusion of crowding effect to the classical Lotka-Volterra competition system guarantees the coexistence of two or more species. We have investigated that two species can coexist when we include crowding effect to the classical LV competition system ( Fig. 1 (solid lines), Fig. 2 as δ → 1.0, Fig. 3a as δ → 1.0, and Fig. 4a,b (solid lines)) compare to the results of the classical LV competition model which show the competitive exclusion principle (Fig. 1 (dotted lines) and Fig. 4a,b (dotted lines)).
Crowding effect has been recognized in many natural and experimental populations [7][8][9][10][11][12][13][14][15][16] . Note that if the density of population is increased a lot more than the carrying capacity, then crowding effect will kill all competing individuals. We here, introduced the intraspecific crowding effect into the Lotka-Volterra competition model. We have shown that a weak crowding effect make it possible to achieve a stable coexistence of multiple species. Our analysis implies that increased mortality under high density works as elevated intraspecific competition leading to the coexistence. This may be another ubiquitous mechanism for the coexistence of multiple species leading species diversity in nature.
We compare the current model with the aggregation model of coexistence 8 (Fig. 5). Unlike the original aggregation model of Hartley and Shorrocks in which the aggregation reduces growth rates, we only include the aggregation effect on mortality rate. These examples show that the aggregation model becomes difficult to maintain the coexistence of all or many species when the number of species is increased. In contrast, the coexistence of all or many species can be easily achieved in the current crowding model even when the number of species is increased (Fig. 5a,c and e). The reason why the coexistence becomes difficult when the number of species increased in the aggregation model seems to depend on the combination of the parameters, where the coexistence region is expressed as a polygon in the n-species parameter space which can disappear easily when the number of species is increased. The logic is same with linear programming in n-dimensional space. In contrast, in the crowding model, the mortality rate of any species increases when its density approaches its carrying capacity. Because of this, increasing mortality rate near carrying capacity will keep all the species at the densities much below their carrying capacity. Thus, the intraspecific competition becomes most severe at or near carrying capacity, resulting in the stable coexistence of all species.
We also compare the effects of the nonlinear competition terms with the current crowding effects (Fig. 6). Unlike the LV model with nonlinear competition effects on growth rates 21 , we only include the nonlinear competition effects on birth rate alone since crowding effect is included only on mortality rate, so that these effects can be easily distinguished. In the current parameter conditions, the 2-species model results in the exclusion of one species (Fig. 6b). However, when the number of species is increased, LV model with nonlinear competition terms will lead to the coexistence of all species with an equal density regardless of the parameter combinations (Fig. 6e,h). It is not sure whether the stability can be achieved easily when the number of species increased in the nonlinear competition model. However, the equilibrium density is identical for all coexisting species in the model with nonlinear competition terms. This means that the effects of other species-specific parameters are completely cancelled by the introduction of nonlinear competition effect. In contrast, the coexistence of all or many species can be easily achieved in the current crowding model even when the number of species is increased (Fig. 6a,d and g). By combining the nonlinear competition and the crowding effect, many species survive and converge to different densities (Fig. 6c,f and i). Thus, both mechanisms can promote the coexistence of many species differently.
The most important assumption in our model of crowding effect is the increase in mortality rate at high density. Here, d i represents the proportion of crowding mortality contribution and δ is the power of crowding effect. Therefore, when d i = 0, this model reduces to the classical LV competition model. This assumption should be one natural way to include the crowding effect. However, there may be many other natural ways to include crowding effect, e.g., the aggregation model of Hartley and Shorrocks 8 . We have to wait for empirical studies to verify which way is actually functioned in a natural ecosystem. Note that these functional mechanisms are not exclusive of each other. Thus, the valid mechanisms may be different depending on a natural ecosystem. We should also note that the nonlinearity in the functional responses should be an important factor driving population dynamics and resulting evolution. The nonlinearity in density effect may be biologically inherent and appears as crowding effect at high densities and as Allee effect at low densities. Our studies thus show that the real competitive communities have a much more complicated dynamical system than the classical LV competition system.

Methods
Mathematical models. We consider the modified Lotka-Volterra (LV) competition equations with crowding effect rate m i for species i. In our model, we only consider competition between two species. In addition, carrying capacity K i is set to be equal to 1, i.e., K i = 1. The modified LV competition model is shown on the following equations:  where x i represents the population density of species i where i = 1, 2. In this model, parameter b i represents the birth rate of species i while α ij represents the effect of species j on i where i,j = 1, 2 and i ≠ j. The crowding effect rate m i is given by where parameter m i0 represents the initial mortality factor of species i. Parameter d i represents the density-dependent factor of species i. In addition, the sum of the initial mortality and density-dependent factor of species i must be greater than 0 but less than or equal to 1, i.e.,0 < m i0 + d i ≤ 1. Note that, if the initial mortality factor m i0 is zero then nonlinear crowding effect rate m i will imply that the intraspecific competition is perfectly density-dependent. In addition, if d i = 0 ∀i then m i = m i 0 ∀i which will imply that the modified LV competition model is the same with the classical LV competition model.
Following equations 4 and 11 in the paper of Hartley and Shorrocks 8 , we arrived with the Lotka-Volterra competition model adding the effect of a few more individuals, shown on the following equations: where ε is any positive real number and is the effect of a few more individuals for all species i. Note that, we do not include the crowding effect on birth rates unlike the aggregation model of Hartley and Shorrocks 8 .
In addition, we also used the modified LV competition model of Taylor and Crizer 21 with the inclusion of nonlinear crowding effect for two species. In their model, they add nonlinear competition terms to prevent the population of species 2 to have a smaller effect on the population of species 1 when the population density of species 1 is very small compare to the population density of species 2 and vice versa. Taylor and Crizer's competition model with nonlinear crowding effect is shown on the following equations: Numerical simulations. In order to determine the impact of the inclusion of nonlinear crowding effects to the classical Lotka-Volterra equation we simulate the modified LV competition equations using Anaconda package of the software Python 3.6 22 . Initially, we determine its effect if there are two competing species in a community and later extend it up to 10 competing species. We also determine the effect when we use small and large values of δ. Moreover, we identify the right combination of d i and δ that will allow competing species to coexist. Without losing essential qualitative dynamics, we considered the following parameter ranges in our numerical simulations: • 0 ≤ Initial of x i ≤ 1 for all i; • 0 < α ij ≤ 1 for all i, j; where K i is the carrying capacity of species i and α ij represents the effect of species j on i where i, j = 1, 2 and i ≠ j.