Species competition: coexistence, exclusion and clustering

We present properties of Lotka–Volterra equations describing ecological competition among a large number of interacting species. First we extend previous stability conditions to the case of a non-homogeneous niche space, i.e. that of a carrying capacity depending on the species trait. Second, we discuss mechanisms leading to species clustering and obtain an analytical solution for a state with a lumped species distribution for a specific instance of the system. We also discuss how realistic ecological interactions may result in different types of competition coefficients.

characterizing the interaction among species i and j, more specifically the decrease in the growth rate of species i due to the presence of j. Competitive interactions are characterized by G ij ≥ 0, the situation to be considered here, whereas negative interactions may model situations of mutualism, predation or symbiosis.
In classical ecological niche theory, species are associated to points in an abstract niche space. Coordinates in this space represent relevant phenotypic characteristics, for example size of individuals in a species, or the size of preferred prey, such that intensity of competition is larger if species are closer in this space. For simplicity we assume this space to be one-dimensional (multi-dimensional generalizations are straightforward, as briefly mentioned later). If niche locations can be considered to be a continuum, we can write equation (1.1) as where ψ(u, t) is the population density at niche location u. The integral extends over the full niche space, which could be finite or infinite. For most purposes, equations (1.1) and (1.2) can be considered as equivalent, since the second is obtained from the first in the limit of many close interacting species, and equation (1.1) can be recovered from equation (1.2) for a discrete distribution of species, with G ij = G(u i , u j ), r i = r(u i ) and K i = K (u i ).
It is widely believed that equations (1.1) or (1.2) predict a competitive exclusion leading to a limiting similarity situation (Abrams 1983), in which a pair of species too close in niche space cannot coexist, and one of them would become extinct. However, it is known that the model allows for continuous coexistence of species in some situations (Roughgarden 1979) and refinements on the conditions for this coexistence have been developed with emphasis on the effect of the shape of the carrying capacity function K (u) Szabó & Meszéna 2006). In this context, a particularly surprising result was the finding by Scheffer & van Nes (2006) of a situation-for uniform carrying capacity-which was neither of full coexistence nor of full exclusion, but of clusters or lumps of tightly packed species which did not exclude each other, but were well separated from other clusters so that there was a type of limiting similarity leading to a minimum intercluster distance. Clustering of individuals or entities under competitive interactions of the LV type had already been observed in other contexts (Fuentes et al. 2003;, 2005Ramos et al. 2008), where the mechanism was the diffusive broadening of an otherwise zero-width species or entity. In contrast, the lumps in Scheffer & van Nes (2006) appeared even in the absence of diffusion in niche space, which is the situation also considered here.
The importance of the functional form of the interaction kernel G ij in equation (1.1) or G (u, v) in equation (1.2) was stressed by Pigolotti et al. (2007) for the case of uniform carrying capacity and interactions depending only on differences of niche positions, and found to be relevant in an evolutionary context by Leimar et al. (2008). For that case the positive-definite character of the Fourier transform of G(u, v) = G(u − v) is a condition implying the absence of limiting similarity. Species clustering was reported, but for interaction functions rather different from the Gaussian used in Scheffer & van Nes (2006). In fact, for the Gaussian interaction case most results are extremely sensitive to details such as the implementation of the boundary conditions or weak ecological second-order effects (Pigolotti et al. 2008). Thus, a clarification of the mechanisms leading to species clustering in LV models would be desirable. In addition, the results in Pigolotti et al. (2007Pigolotti et al. ( , 2008 were obtained under the unrealistic assumption of homogeneity in niche space, whereas the inhomogeneities in the carrying capacity are known to play relevant roles (Szabó & Meszéna 2006). For simplicity we restrict our description to the standard situation in which competition is stronger among species closer in niche space. The existence of studies of LV systems where non-local interactions beyond that type are considered is worth mentioning (Doebeli & Dieckmann 2000). That situation can also be described by the general formalism used here of an integral kernel function, and our general results therefore also apply to the situation with more general non-local interactions.
In this paper we analyse some mathematical properties of the LV model (1.1) or (1.2). In §2 we show that the positive-definiteness of the kernel G remains a determining condition for stable coexistence even for non-constant K (u). In §3, we discuss the mechanism producing lumped species distributions and explicitly provide an analytical expression for a particular interaction kernel. In the appendix we show that, in contrast with the earliest characterizations of the interaction kernel G (MacArthur & Levins 1967;Roughgarden 1979), both positive-and non-positive-definite kernels can arise from more detailed ecological models which consider the dynamics of the consumed resource. We use periodic boundary conditions in our numerical simulations. We expect the effects of this simplifying but unrealistic assumption to be unimportant at least when a nonconstant carrying capacity limits the presence of species to a limited region of niche space.

The stability of close coexistence
A simplifying assumption for the study of the LV model is that of homogeneity in niche space. In this case, the carrying capacity and growth rate are constants, K 0 and r 0 , and the interaction kernel depends only on differences of niche positions G(u, v) = G(|u − v|). Niche space could be infinite, but, in the case in which it is finite, homogeneity can only be achieved under periodic boundary conditions. Under these restrictions it is easy to see that a steady solution to equation (1.2) which is homogeneous and everywhere non-vanishing always exists: ψ 0 = K 0 /Ĝ 0 , whereĜ 0 ≡ duG(u). This solution represents coexistence of all possible species without a limit to their similarity. Its stability against small perturbations can be analysed by linearization of the equation resulting from substitution of ψ(u, t) = ψ 0 + δψ(u, t) in equation (1.2). The solution for the Fourier transform of the deviation from the homogeneous state, δψ q (t), is whereĜ q is the Fourier transform of G(u). Thus, the homogeneous solution ψ 0 is stable ifĜ q is positive ∀q, while a instability leading to pattern formation occurs whenĜ q may take negative values (Fuentes et al. 2004;Pigolotti et al. 2007). We note that many steady solutions to equation (1.2) exist besides ψ 0 (in particular, solutions of the form equation (1.3)). This is so because dynamics preserves ψ(u) = 0 at all places where there is no initial population. Notice also that ψ 0 is the only strictly positive solution. Among this multiplicity of solutions the ones that will be more relevant are those which are stable under perturbations or small immigration (Pigolotti et al. 2007). An interesting class of functions to be used as kernels and carrying capacities is the family {g p σ } given by which is parameterized by the value of p. The widely used Gaussian kernel is obtained for p = 2. When p < 2 the functions are more peaked around u = 0 (the case p = 1 is an exponential) and for p > 2 they become more box-like (g ∞ σ (u) is the flat box with value 1 in the interval [−σ , σ ] and zero outside). The width of the kernel σ gives the competition range in niche space. We have positivity of the Fourier transform if p ≤ 2. This implies that the homogeneous solution is stable under evolution with uniform K and kernel G of the form equation (2.2) if p ≤ 2. When p > 2, the homogeneous solution is unstable and the system approaches delta comb solutions of the type equation (1.3), with a spacing approximately 1.4σ (Pigolotti et al. 2007) that represents limiting similarity situations.
We now generalize the above stability analysis to the more realistic case in which there is no homogeneity in niche space. First, we consider the simpler case of a symmetric kernel G(u, v) = G (v, u), which in particular includes the previous case of kernels depending only on species distance: G(u, v) = G(|u − v|). Note that in this symmetric case one can write equation (1.2) in potential form, with the functional potential given by Stationary solutions of equation (1.2) are those for which the right-hand side of equation (2.3) equals 0. This has many possible solutions. We define the natural stationary solution, ψ N (u), as the one which is positive and non-vanishing for all u, so that δV that is, the one satisfying The solution ψ N (u) can be considered as the non-homogeneous generalization of ψ 0 introduced in the homogeneous case. In the particular case in which , the natural solution can be explicitly written in terms of Fourier transforms of the competition kernel and the carrying capacity, either in an infinite system or in a finite one with periodic boundary conditions, ( 2.7) This requires that these Fourier transforms and their inverses exist and lead to positive population densities. When this happens, a continuum species coexistence is obtained, and its existence is generally robust against small changes in G or K . Later we show that it is also an attractor of the dynamics when G q satisfy positivity requirements (p ≤ 2, for the family in equation (2.2), being p = 2 the marginal case). For a uniform carrying capacity, the natural solution equation (2.7) always exists and is uniform in phenotype space ψ N (u) = ψ 0 . But the natural solution may lose positivity or even cease to exist depending on the properties of G and K . For example, when both G(u) and K (u) are of the form equation (2.2) with p = 2, the inverse Fourier transform of equation (2.7) exists when the carrying capacity has a value of σ larger than the kernel G, but not in the opposite case. Figure 1 shows stationary solutions attained at long times by the dynamics in equation (1.2) illustrating the situations described above, starting from a smooth non-vanishing initial condition. In the first case we choose a kernel and carrying capacity functions (G(u) = g 1 σ (u), K (u) = sech(u/σ )), such that the natural solution exists and is positive everywhere. Thus it is stable, and it is the steady state attained at long times. In fact it can be analytically calculated, ψ N (u) = a −1 sech 3 (u/σ ). (2.8) In the second case the non-positiveness of the kernel used (with a carrying capacity of the type equation (2.2)) breaks down the initial configuration into lumps, which at long times approach zero-width delta functions with forbidden zones in between. In the third case, despiteĜ q being positive, a positive natural solution does not exist. Several outcomes are possible but, for the kernel and capacity used, the system approaches a single hump solution which vanishes in part of the niche space.
More in general, but still in the symmetric G case G(u, v) = G (v, u), writing the LV model in potential form (equation (2.3)) is of great use since one can show that, provided r(u) and K (u) are positive, dV /dt ≤ 0. This implies that V is a Lyapunov potential and dynamics proceeds towards its absolute minimum, or if ψ(u, t = 0) = 0 for some u, towards the minimum of V under such constrain. Notice that, since the potential is a quadratic form, ψ N is a global attractor (starting from non-vanishing initial conditions) when the competition kernel is a positive-definite quadratic form, which means that f (u) G(u, v) and shows that the stability result was global indeed. In a multi-dimensional niche space, the same analysis shows that the positive-definiteness of the quadratic form remains the condition for the global stability of the natural solution. In any case, the important consequence is that the stability of the natural solution depends uniquely on the competition kernel and not on the carrying capacity (provided the relation kernel-capacity is such that the natural solution exists and is positive). In particular, for competition kernels of the form equation (2.2), ψ N is always (if existing and positive) a globally stable solution of the dynamics for p ≤ 2, and unstable otherwise.
The crucial difference in the case of a non-symmetric competition kernel is that there is no obvious Lyapunov potential for the system. This implies that there are no global stability results available. However, local stability can be investigated. Let us consider a small perturbation of the positive natural solution ψ N (u) + δψ (u, t). To linear order, the perturbation evolves as (2.9)

We now consider the functional H (δψ) ≡ du (A(u)K (u)/r(u)ψ N (u))(δψ) 2 ,
where A(u) is a positive function so that H ≥ 0 and H (0) = 0. Let us compute its time derivative, (2.10) If for some choice of A(u) one has that A(u) G(u, v) is positive definite, then dH /dt < 0 and δψ = 0 will be approached. This shows that ψ N is linearly stable in such a case. We note that the case in which G (u, v) itself is positive definite trivially guarantees the positivity of A(u) G(u, v), with a constant A. Thus, even in this more general non-symmetric case, it is the character of the interaction kernel G, and not of the carrying capacity (provided it is such that the natural solution exists and is positive), which determines the stability of the natural solution. Scheffer & van Nes (2006) found transient but long-lived solutions of equation (1.1) consisting of periodically spaced lumps containing many close species. They used a Gaussian interaction kernel which turned out to introduce an excessive sensitivity of the results to the numerical implementation of the model and boundary conditions (Pigolotti et al. 2008). However, they found similar solutions as steady configurations when adding an extra predation term acting effectively only on species with high population. This can be thought of as an extra intraspecific competition since it decreases the growth of species with many individuals. Exploiting this idea, Pigolotti et al. (2007) checked the effect of using in equation (1.2) a kernel of the type equation (2.2) but with an enhanced interaction at u = 0, i.e. enhanced intraspecific competition. In particular, they used a constant carrying capacity K (u) = K 0 and a flat box kernel with an added delta function at the origin (figure 2),

Lumped species distributions
Lumped patterns were obtained numerically for a = 1. Because the dynamics of equation (1.2) usually involves very long transients, it is interesting to calculate analytically the steady lumped solution in the simple case of a kernel (3.1) and uniform carrying capacity K 0 (in the infinite line).
We begin with the steady-state condition valid at u such that ψ(u) = 0, that particularized to equation (3.1) and constant K reads This is transformed into a differential-difference equation after differentiation with respect to u, This steady equation has many solutions, including the natural one ψ 0 = K 0 /(a + 2σ ) which is non-vanishing everywhere, or delta combs such as equation (1.3).
We search for solutions of the type in figure 2, i.e. periodic arrays of lumps, of period d, each one having a symmetric hump shape We are assuming that the lumps do not overlap, so that d > 2L. We also note that if σ + 2L < d there is no interaction between different lumps, so that for u ∈ [−L, L] equation (3.4) reduces to f (u) = 0 and there is no lump solution. Moreover, analysis is much simplified if each of the lumps interacts only with its neighbours (σ + 2L < 2d). Thus we restrict to d < σ + 2L < 2d, for which equations (3.4) with (3.5) and u ∈ (−L, L) becomes The general solution of this linear equation is obtained as a sum of exponentials exp(λu), with aλ = sinh(λ(d − σ )).  with , (3.9) and the value of λ which is plotted in figure 3. Figure 2 shows the analytical solution equation (3.5) with equations (3.8) and (3.9). We have not studied the stability of this configuration. But the numerical results in Pigolotti et al. (2007) indicate that it is stable for some values of L and d.
We finally stress that the appearance of the lumped solution is not a consequence of the singularity of the delta function in the kernel. In fact, any kernel sufficiently peaked at the origin will favour the coexistence of close species. If the behaviour at larger distances of the kernel makes it not positive definite, then full coexistence will be unstable and the natural solution will split into disjoint lumps. An example of the final steady state in this situation is shown in figure 4, with a kernel G = g 4 0.2 + 0.8g 1 0.02 , which has the properties just described and contains no delta singularity.
where r i = α S iα Q α is the maximum growth rate and C ij = α S iα D αj . Thus, the result is an effective interaction among the predators which is of LV type. It is customary to write equation (A 5) in terms of the carrying capacity K i , defined as the equilibrium population N i attained in the absence of the other competitors, i.e. K i = r i /C ii . In terms of it, equation (A 5) becomes identical to equation (1.1), with Having a continuum R(x) of resources instead of a discrete set R α does not introduce important difficulties. Simply, one should replace sums by integrals, replacing the coefficients of equation (A 5) by One can also consider a continuum of species, labelled by their phenotypes u, so that equation (1.1) is replaced by equation (1.2) with K (u) = r(u)/C (u, u), G(u, v) = C (u, v)/C (u, u), and r(u) and C (u, v) given by obvious generalizations of equations (A 7) and (A 8).
It is clear that the presence in the kernel G ij of two different functions (compare with the most restrictive expression (A 1)) gives enough freedom to obtain a variety of kernel behaviours under different circumstances. A particularly popular choice is to assume that impact and sensitivity are proportional: S iα = D αi , with a constant efficiency . In the continuum resource case the functions can be written in terms of a single utilization function u i (x) as D i (x) = u i (x) and S i (x) = u i (x), leading to the classical expression (A 1). Slightly more general cases arise when the efficiency depends only on the resource, = α , or on the consumer, = i , or when dependence on the two types of species factorizes, = v i w α . In all these cases (if the efficiency is positive) one is led to a kernel G ij , which is positive definite. In more general cases, one can have a kernel leading to instability of the coexistence state.
We conclude with two instances of ecological interactions in equation (A 3) which allows the stability to be tuned. First, a homogeneous discrete and infinite niche space in which all resources have the same internal dynamics Q α = Q, β α = β, ∀α, as well as the consumers: d i = d, ∀i. The interactions are taken to be This models a situation in which the consumer k grows only by consuming its optimal resource R k , whereas it depletes also the neighbouring resources, R k+1 and R k−1 . We have r i = Qg, K i = β/a, C ij = (Qg/β)(aδ i,j + b(δ i,j−1 + δ i,j+1 )), and G ij = δ ij + (b/a)(δ i,j−1 + δ i,j+1 ) so that equation (1.1) is noẇ The natural solution, i.e. the one in which all species have positive non-zero population, is N i = β/(a + 2b), ∀i. Its linear stability can be studied by linearization, N l (t) = N l + δN l (t) and substitution of the ansatz δN l ≈ e λ q t e iql (here i = √ −1). We find λ q = −(Qg/β)(a + 2b cos q), q ∈ [−π, π]. λ q are the eigenvalues of −C ij , and stability of N i requires all these eigenvalues to be negative, i.e. C ij to be positive definite. When a > 2b, then λ q < 0 ∀q, and the natural coexistence solution is globally stable (see results in §2). It is unstable otherwise. In this example, there is no well-defined single utilization function and the positivity properties of the interaction kernel and thus the stability of the natural solution can be changed by varying the parameters.
As a second example, with non-constant carrying capacity, we consider a continuous distribution of resources and species on the line, and we take which implies that the consumers of phenotype u grow only from the resource at location x = u, but they deplete a wider range characterized by f . This leads (in the continuous formalism) to r(u) = sQg(u), C (u, v) = sQf (u − v)/β, K (u) = βg(u)/f (0) and G(u, v) = f (u − v)/f (0). In this example, by choosing the functions g and f , we can impose any desired combinations of carrying capacity and interaction kernel. Gaussianity or positive-definiteness are particular cases, no more natural in this generalization than alternative choices leading to non-positiveness, instability and thus exclusion zones between clumps of species.