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.

Lotka–Volterra competition and species distribution

Competitive interactions occur when entities in a system grow by consuming common finite resources. They are ubiquitous in many fields of science: examples include biological species competing for food (MacArthur & Levins 1967; Roughgarden 1979; Case 1981), mode competition in nonlinear optical systems (Benkert & Anderson 1991), or alternative technologies competing for a market (Pistorius & Utterback 1997). An early, simple, but powerful model for competitive interactions is the Lotka–Volterra (LV) set of competition equations (Volterra 1926; Lotka 1932), where m is the number of species, N is the population of species i, r is its maximum growth rate, K is its carrying capacity and G is the matrix 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≥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=G(u,u), r=r(u) and K=K(u).

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) (Meszéna ; 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 ; Hernández-García & López 2004, 2005; Ramos ), 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 in equation (1.1) or G(u,v) in equation (1.2) was stressed by Pigolotti 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 . 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 . Thus, a clarification of the mechanisms leading to species clustering in LV models would be desirable. In addition, the results in Pigolotti , 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 non-constant 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, K0 and r0, 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: , where . 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, , is where is the Fourier transform of G(u). Thus, the homogeneous solution ψ0 is stable if is positive ∀q, while a instability leading to pattern formation occurs when may take negative values (Fuentes ; Hernández-García & López 2004; López & Hernández-García 2004; Pigolotti ). 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 ).

An interesting class of functions to be used as kernels and carrying capacities is the family 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 ( 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 ) 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, ψ(u), as the one which is positive and non-vanishing for all u, so that that is, the one satisfying The solution ψ(u) can be considered as the non-homogeneous generalization of ψ0 introduced in the homogeneous case. In the particular case in which G(u,v)=G(u−v), 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, 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 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 ψ(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 (, 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,

Figure 1.

Long-time solutions of equation (1.2) for different kernels and carrying capacities. (a) , K=sech(u/σ), with σ=0.1. The natural steady solution (ψ= a−1sech3(u/σ)), which is positive and non-vanishing everywhere, is reached at long times. (b) , . Under this non-positive-definite competition kernel, the solution shown is still evolving and approaches a singular delta comb of the type (1.3) at long times. (c) , . A positive natural solution does not exist and the system approaches a single hump solution which vanishes in part of the niche space.

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 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, ψ is a global attractor (starting from non-vanishing initial conditions) when the competition kernel is a positive-definite quadratic form, which means that , ∀f (or , ∀{x} in the discrete case). This generalizes the previous stability condition on the Fourier transform to niche inhomogeneous cases, 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), ψ 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 ψ(u)+δψ(u,t). To linear order, the perturbation evolves as We now consider the functional , where A(u) is a positive function so that H≥0 and H(0)=0. Let us compute its time derivative,

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 ψ 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.

Lumped species distributions

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 ). 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 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)=K0 and a flat box kernel with an added delta function at the origin (figure 2), Lumped patterns were obtained numerically for a=1.

Figure 2.

The kernel in equation (3.1) (a), and the analytical steady solution given by equations (3.5) and (3.8–3.9) for a=K0=1, σ=0.8, L=0.3 and d=1 (b).

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 K0 (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=K0/(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 f(u) of width 2L (i.e. f(u)=0 if u∉[−L,L]), We are assuming that the lumps do not overlap, so that d>2L. We also note that if σ+2La the only solution is the constant one, but in the opposite case (the situation favoured by enhanced intraspecific competition a) the solution is a linear combination of three exponentials. Two of the constants of the combination are determined from f(L)=f(−L)=0. The third one, which gives the overall normalization, can be obtained by returning back to the original equation (3.3). Finally we get with 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 indicate that it is stable for some values of L and d.

Figure 3.

The solution λ (positive branch) of equation (3.7) giving the inverse width of species lumps. The width is finite for 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 , which has the properties just described and contains no delta singularity.

Figure 4.

The kernel (a), and the steady solution obtained numerically from it at long times with constant K0=1 (b).