Self-localized states in species competition

We study the conditions under which species interaction, as described by continuous versions of the competitive Lotka-Volterra model (namely the nonlocal Kolmogorov-Fisher model, and its differential approximation), can support the existence of localized states, i.e. patches of species with enhanced population surrounded in niche space by species at smaller densities. These states would arise from species interaction, and not by any preferred niche location or better fitness. In contrast to previous works we include only quadratic nonlinearities, so that the localized patches appear on a background of homogeneously distributed species coexistence, instead than on top of the no-species empty state. For the differential model we find and describe in detail the stable localized states. For the full nonlocal model, however competitive interactions alone do not allow the conditions for the observation of self-localized states, and we show how the inclusion of additional facilitative interactions lead to the appearance of them.


I. INTRODUCTION
The interactions among the biological entities integrating an ecosystem give rise to surprising emergent behaviors. Competition is one of the most important and ubiquitous of these interactions: if there is an increase in the population of one species, due to the consumption of common resources or shared predators [1], there is a decrease in the growth rate of the others. Because of this interaction it is usually argued that a given ecosystem can host only a limited number of species that should be sufficiently separated from each other in the so-called niche space. This is the d-dimensional space whose coordinates (x, y, . . . ) quantify the traits of the species relevant for the utilization of the resources distributed as a function of these coordinates. The competitive exclusion principle [2] is a formulation of this situation, in which species can not coexist too close in niche space (limiting similarity). Despite this reasoning, however, it should be said that even the most traditional mathematical model of competitive species, the Lotka-Volterra (LV) model [3], is known to allow solutions characterized by a continuous distribution of species [4] under some circumstances (see reviews in [5,6]). More remarkable in this context are recent results on the LV model (or closely related ones) showing the existence of solutions that do not represent purely continuous coexistence, nor are typical of a limiting similarity situation [5,[7][8][9]: clusters of species around particular niche positions well separated from each other and filling out the niche space. For the competitive LV model, these lumped distributions appear due to pattern forming instabilities triggered by the shape of the interaction function [8][9][10][11]. Besides these many-cluster * Email address: damia@ifisc.uib-csic.es; telephone: +34 971259837; fax: +34 971173248. species configurations, an ecologically relevant question is under which conditions solitary clusters of species may appear. These would arise from an evolutionary or random drift towards a particular niche position, or simply from an advantageous initial condition. In this Paper, we address this question in the context of pattern formation in continuous versions of the LV model, both in an integral formulation as in its differential approximation. Our focus is on competitive interactions, but we will be forced to consider also some facilitative (i.e. mutualistic or symbiotic) situations. Through the paper we will keep in mind the situation of species competition in niche space, but we stress that the concepts and type of models used here are equally valid to describe organisms randomly moving in physical space and nonlocally competing for resources with other individuals in their spatial neighborhood [12,13], or rather in evolutionary situations [14,15].
A pattern-forming instability or bifurcation is a source of great complexity, and many different scenarios may arise from it. One of the simplest cases is the formation of a periodic structure, that in two or higher dimensions can have different geometries depending on the nature of the nonlinearities. In some cases the bifurcation can be subcritical so that periodic patterns can coexist with homogeneous distributions. In this case localized solutions consisting of one or more isolated lumps on top of a homogeneous distribution might exist, being supported by the nonlinearity and the spatial coupling, as shown in general amplitude equations [16,17]. If this mechanism turns out to be present in the context of biological competition, then a stable localized lump could be formed in a given stable ecosystem supporting a continuous coexistence of species. Such lumps can be formed at any position in niche space triggered by particular perturbations or initial conditions. This means that species with no special advantage with respect to their competitors might prevail at some point due to a particular initial condition. These high values of the population of certain species would be supported by the nonlinear dynamics and the spatial interaction, and not by a better fitness to the ecosystem. The LV competition model in niche space turns out to be a nonlocal model, i.e., species interact with others not located closely in the niche axis. Population dynamics has revealed many different interesting phenomena due to nonlocal competition [8,[18][19][20][21][22], such as periodic patterns, discrete clusters, defects and fronts in space, etc. Self-localization has been broadly studied in physical systems [23][24][25] but much less in the context of population dynamics [21,22,26]. Previous works have already found self-localization of biological entities by inclusion of the Allee effect, i.e. a tendency to extinction when population numbers are too small, in nonlocal competition models [21,22,26] with cubic nonlinearity. The Allee effect naturally induces bistability between the empty or extinct state and the natural occupation determined by the carrying capacity. This bistability allows the existence of self-localized patches of densities close to the carrying capacity surrounded by empty space. In this case the bistability involves two different spatially homogeneous states [27]. In this Paper, in contrast, we address the situation involving coexistence of a spatially homogeneous state and a spatially periodic pattern [16,17]. Also, we consider always positive linear growth rates, so that the Allee effect is absent and a small population will always grow, and we use only quadratic nonlinearities. In consequence, we are looking for localized structures on top of a non-zero homogeneous density, instead of the localization on top of an unpopulated background as described previously [21,22,26]. Thus, we are considering the possibility that the interaction enhances the density locally, but without driving to extinction the rest of the system.
The structure of this article is as follows: In Section II our nonlocal model for species competition is introduced. In section III we approximate this model by a partial differential equation (PDE) which reproduces the basic original results, and allows us to use methods for the analysis of self-localized solutions in PDEs. In Section IV we present the results of our analysis for this differential model showing under which conditions localized solutions can be found. Then, in Section V we discuss the features the nonlocal interaction kernel must have to observe localized states in the full nonlocal model. Finally, in Section VI we give some concluding remarks. The Appendix briefly summarizes details of the numerical methods.

II. THE NONLOCAL KOLMOGOROV-FISHER MODEL FOR SPECIES COMPETITION
The classical Lotka-Volterra model of N species in competition, each utilizing a common distributed re-source x is given by [4,7,8] where the dot denotes temporal derivative, n i is the population of species i, r is the growth rate (that we assume the same for all species), and x i is the position of the species i in the niche axis (for simplicity we work in one dimension). G(x) is the interaction kernel, which unless explicitly said will take positive values to model competitive interaction. We also assume G to depend only on the modulus of the relative difference |x i − x j |, meaning that resources are homogeneously distributed in niche space and interactions are isotropic there. G sets the scale of the carrying capacity, which is then also the same for all niche positions. More complex situations are reviewed in [5].
If niche locations are considered to form a continuum (the infinite real line), we can write the former equation as:Ψ where Ψ(x) is now the population density (always positive), andG is an integral operator describing the competition term:G A further step in the modeling is considering diffusion in niche space, that may account, for instance, for mutations in an evolutionary context, or random phenotypic changes [11,28]: where D is the diffusion coefficient. Note that Eq. (4) is a type of nonlocal Kolmogorov-Fisher-like equation [12,18,[29][30][31][32]. This type of equation may also describe organisms randomly moving in physical space and nonlocally competing with other individuals for resources [12,13]. In that case D is a true diffusion coefficient modeling random dispersion in space. It has been shown that arbitrarily small structural perturbations of this model away from having a constant r may destroy the continuous all-positive solution for zero diffusion [6,33]. The presence of diffusion in our case ensures, however, that the homogeneous solution only deforms continuously under small perturbations from a constant r. In this case the existence and dynamical properties of self-localized states are not drastically altered, as shown, for instance, in a nonlinear optical system [25].

III. TRUNCATION OF THE NONLOCAL OPERATOR
In order to analyze the existence of localized states in Eq. (4), we first approximate the nonlocal operator (3) by a simpler differential operator. This will allow us to apply standard techniques for PDEs to find localized states. To do so we Taylor expand the function G in the nonlocal operator to obtain a series of derivatives of Ψ: Because of the assumed isotropy of G, G = G(|z|), all terms G n with odd values of n are zero. The k-Fourier component of the convolution integral operator can be written as: where the hat indicates Fourier transform. From Eq. (5), one can also find (for isotropic systems) that: We illustrate the above manipulations with a relevant class of competition kernels widely discussed in [8,9,11]: p describes how steep the edge of the kernel is and σ is the range of the competition. Note that p = 2 is the Gaussian kernel, and p = 1 the exponential one. Species consuming very different resources, i.e. with a large distance between them in niche space (|x−s| >> σ) interact very weakly, while species which are close (|x − s| < σ) compete significantly. Finally, a accounts for the strength of the competition and sets the scale of the carrying capacity. In Figure 1 we plot examples of the typical kernel (10) for two different values of parameters, showing the important role of parameter p. The Fourier transform of the function (10) is positive and tends to zero monotonously for k → ∞ when p < 2 (see dash-dotted line), however for p > 2 negative components appear in the Fourier transform (see solid line) [34]. This leads to a modulational instability of the homogeneous solution [8], as detailed later. Figure 1 also shows, for p = 6, how the Taylor decomposition approaches the full convolution kernel. The line marked by squares shows function (9) for a series of only three terms (G 0 , G 2 , and G 4 ). The line marked by crosses shows the approximation with two more terms (G 6 , G 8 ), and circles show the approximation by terms up to G 12 . The major differences between Fourier transform of the full nonlocal operatorĜ(k) and the approximation (9) occur at high values of k. We keep however only three terms in the series and perform the analysis as an intermediate step towards understanding of the original model. In this approximation the opera-torG becomes the Swift-Hohenberg operator or shifted diffusion, and Eq. (4) reduces to: At difference with the original Swift-Hohenberg equation, however, in this model the spatial operator appears in nonlinear terms.
Within this approximation we interpret coefficients G 0 , G 2 , and G 4 , characterizing the kernel, as independent parameters. This allows us to analyze the solutions of this model more accurately and extract later the features a kernel must have to access a given region of this parameter space.
Truncating the nonlocal operator to obtain a local model can be a rough approximation, however, it describes appropriately stationary distributions Ψ(x) provided G n k nΨ (k) tends to zero fast enough as k tends to infinity. Actually, we find that a set of parameters G 0 , G 2 , G 4 close to the ones obtained from (6) gives a good qualitative agreement between the stationary solutions of the two models, as presented in Figure 2. The patterns are calculated by solving Eqs. (4) and (11), starting from slightly (randomly) perturbed unstable homogeneous solutions as initial conditions. In both cases we observe the formation of "lumps", separated by less populated regions. The similarity of the results justifies the consideration of (11) in the following sections. One can note also that function (10) decays very fast for |x − s| → ∞. This means that the more narrow is the kernel the more local is the system, and the validity of the truncation of the Taylor series is better. A quantitative evaluation of the effect of the truncation at a certain order can be obtained for each k by comparing Eq. (7) with Eq. (8).
Equations (4) and (11) have two homogeneous solutions: Ψ = 0 and Ψ = Ψ 0 = r/G 0 . The zero solution corresponds to the situation in which niche space is not occupied, and it is always unstable for positive growth rates r. Any small number of individuals will be able to reproduce and the population will grow approaching the steady homogeneous state Ψ 0 . In the absence of diffusion (D = 0), this solution is modulationally unstable for kernels whose Fourier transform contains negative components [8,13,35], and lumped distributions over niche space arise instead. For kernels given by Eq. (10) negative Fourier components appear when p > 2 [34].
Adding diffusion D > 0, the stability condition is changed and a threshold value appears. In the case of Eq. (11), Ψ 0 is stable for G 2 < G th 2 , with For G 2 > G th 2 a pattern with periodicity determined by the critical wavenumber appears. For D = 0, condition (12) is equivalent to the condition of appearance of negative components in the Fourier transform of the kernel. Since localized solutions are usually found in parameter regions where a periodic pattern coexists with the homogeneous solution [16,17], we look for the conditions in which the pattern-forming bifurcation is subcritical. The way to do it (a weakly nonlinear analysis) is described in [36,37]. Introducing formally a small parameter ε, we write the solution and control parameter (we choose here G 2 ) as: where Ψ 0 is the homogeneous steady state. Substituting (14) and (15) into the stationary version of (11) and collecting terms at different orders of ε we obtain: wherẽ is the Jacobian of Eq. (11) evaluated at G 2 = G th 2 , and At first order Ψ 1 = A cos(k c x), which is the periodic solution bifurcating at G 2 = G th 2 . At second order, the solvability condition leads to G 21 = 0 and Ψ 2 = B + C cos(2k c x), with B = DG 0 √ G 0 A 2 /2r 2 √ G 4 , and C = DG 0 √ G 0 A 2 /18r 2 √ G 4 . Finally, the solvability condition at third order 2π/kc 0 f 3 Ψ 1 dx = 0 (21) leads to the following equation for the stationary amplitude A of the critical mode found at the first order: where The transition from a super to a sub-critical bifurcation occurs when the coefficient κ changes sign (κ = 0). This happens for For the full nonlocal operator, this condition is equivalent to setting the coefficient κ of Eq. (28) in Ref. [13] to zero.

IV. SELF-LOCALIZED SOLUTIONS
In the following we study the existence of localized solutions in Eq. (11). This equation has only two independent parameters, so that by rescaling t, x, and Ψ we can consider G 0 = 1, G 4 = 1, and r = 1 without loss of generality, and take G 2 and D as control parameters. The condition for instability of the homogeneous solution (12) becomes then: and the pattern appears subcritically if To illustrate the change from a supercritical to a subcritical bifurcation we plot the bifurcation diagram of the stationary pattern solution of (11) arising at G th 2 for two values of D, one below and one above the critical value D s (See Fig. 3).
The codimension-2 point indicated by D = D s and G 2 = G th 2 is called in the spatial dynamics parlance a Degenerate Hamiltonian-Hopf bifurcation, and it is known to be the origin of the existence of localized states in pattern forming systems [16]. In the following we focus on the existence of self-localized states consisting on a number of stable lumps of the pattern solution on top of the homogeneous solution Ψ 0 . For this we move well into the  parameter region where the pattern forming bifurcation is subcritical by increasing D.
Using a shooting method (see Appendix) where the spatial coordinate x is used as a dynamical variable in the stationary version of Eq. (11), we have found a selflocalized solution. Taking it as an initial guess we have computed its branch by continuation techniques using a Newton method. Figure 4 shows the bifurcation diagram of localized states with an even number of peaks for D = 1.7. Another analogous curve for localized states with odd number of peaks (not shown) also exists. The curve shows a characteristic snaking structure. The branch follows a series of saddle-node bifurcations that transform unstable solutions into stable localized states, adding each time a peak at each side of the structure as one moves up. Typical localized states, indicated by bold dots in the inset of Figure 4, are presented in Figure 5. Solutions (b) and (d) in figure 5 are stable, while the rest are unstable. As can be seen from the inset in Figure  4, the region of existence of stable self-localized states is very narrow due to the proximity to the codimension-2 point. This region becomes larger as one moves away from this point in the direction of increasing the subcriticality of the pattern, i.e. increasing D. However we can not go much further into that region, because the minimum of the population distribution approaches the trivial zero solution too much, and our simulations diverge. In this case further analysis is not possible and saturating terms should probably be included in Eq. (4) in order to observe stable localized states in wider parameter regions. We are unable to determine if the difficulties are only of numerical origin or if there is some more fundamental change of behavior or bifurcation when increasing subcriticality, perhaps associated to some spatial analog of the paradox of enrichment [38]. Further investigation is needed to clarify this point.
The localized solutions found (Fig. 5) consist on a few lumps of species, of a very high population density, which locally deplete close niche positions but do not make them completely empty. Further apart the effect of the lumps becomes unimportant and the density in the rest of the system consists on the stable homogeneous coexistence of species given by the homogeneous solution Ψ = 1. The spacing between the lumps forming the localized patch is of the order of the periodicity of the extended pattern (Fig. 2). These regions can be then considered as portions of the periodic pattern embedded inside the stable homogeneous solution. To illustrate the stability of the self-localized states (b) and (d) in Fig. 5, we show in Figure 6 the switching dynamics of localized states starting from suitable initial conditions. The stability of the states has also been checked with respect to small additive noise.

V. LOCALIZED STATES IN THE FULL NONLOCAL MODEL
Once the precise conditions for the observation of subcritical patterns and localized structures have been determined for Eq. (11), we can discuss the implications for the kernel in the full nonlocal model (4). The main assumption is that we can apply to (4) the results of the previous section by using the values G 0 , G 2 , etc. arising from the expansion of the nonlocal kernel.
The necessary conditions for the existence of selflocalized states were (for G 0 = G 4 = r = 1) the subcriticality criterion (26), D > D s = 27/23, and the instability condition (25), G 2 > G th 2 = D + 2. In addition, the subcritical region increases with increasing D, but as commented above our numerical results were unable to probe large values of D without divergences. Figure ( But it happens that this range of values of G 2 is far from what is achievable with an interaction kernel of competitive nature exclusively (i.e. one taking only positive values). To see this we note that a positively-defined and normalized (G 0 = 1) kernel can be interpreted as a probability density, so that G 2 and G 4 are its moments (see Eq. (6)): G 2 = x 2 /2, G 4 = x 4 /24. If G 4 = 1, then x 4 = 24. Applying the moment monotonicity inequality: |x| r 1/r ≤ |x| s 1/s , where 0 < r ≤ s, and using r = 2 and s = 4, we have (2G 2 ) 1/2 ≤ 24 1/4 , or G 2 ≤ √ 6 ≈ 2.449. This limiting value is well below the ones needed to observe self-localized solutions of Eq. (11) without encountering divergences. We can not completely discard in a rigorous manner the possibility of stable localized structures to exist for the nonlocal model at sufficiently large values of D, nor the presence or other localized solution branches not captured within the differential approximation. But the fact is that we have been unable to find numerically self-localized solutions of (4) when using a purely competitive (i.e. positive) kernel G(x).
A natural way to achieve the larger values of G 2 needed is to allow the kernel to take negative values close to x = 0 or for large values of x. This means the presence of facilitative interactions (mutualism, symbiosis, ...) together with the competitive ones. We note that such combination of positive and negative interactions at different length scales was already proposed from biological reasoning in an early paper [29], an it is an important ingredient in the modeling of vegetation patterns [39].
Here we consider an integral kernel G I of the form: were a 1 can take negative values modelling cooperative or facilitative interactions. To find localized states in the full nonlocal model we perform then a continuation of the localized states from the differential to the integral kernel. The main difference between these two cases consists in the behavior of G(k) for k → ∞: in the differential case, G(k) → ∞, while for the integral case G(k) → 0, as illustrated in Fig. (1). Since the stability range of localized states is so small it is a challenge to find the parameters of the nonlocal kernel that support stable localized states. We show now, however, that at least for the lowest unstable localized state marked by a dashed red line in Figure (4) our continuation strategy is able to find them. To do so, we consider a linear combination of the integral kernel G I and the differential approximation in the truncated model (11) G D in the form: where γ is a parameter characterizing how differential or how integral the resulting kernel G(k) is. So, for γ = 1, the kernel is purely differential, while for γ = 0, the kernel is purely integral.
We choose the parameters in such a way that the Fourier transform of G I is very close to the one of G D used in Fig. 4, except for high values of k [see Figure  7 (b)]. This implies that in real space the kernel G(x) takes negative values close to x = 0 [ Figure 7(a)] . Now, using a continuation method we follow the self-localized solution from the differential case γ = 1 to the integral case γ = 0. The corresponding modification of the kernel (28) and of the profile of the localized separatrix solution is shown in Fig. 7(c). In such a way we demonstrate the existence of localized states in the original model (4) for kernels fulfilling appropriate conditions.

VI. CONCLUSIONS
By studying a differential truncation of the nonlocal Kolmogorov-Fisher model we have rigorously calculated the conditions by which periodic patterns are subcritical and we have demonstrated numerically the existence of the stable localized states in the differential approximation. These are patches of finite extent containing a number of lumps of species and arise on top of the homogeneous distribution. In contrast to other works, our results show that the necessary ingredient to observe stable self-localized states, namely the presence of subcritical patterns, is already present in systems with spatial coupling in the quadratic nonlinearity only, rather than nonlinearities of different orders being necessary. In consequence, the localized patches appear on a background of homogeneously distributed species coexistence, instead than on top of the no-species empty state. Extending the results obtained for the truncation to the full nonlocal model we find, however, that competitive interactions alone can not lead to the conditions for the observation of localized states, and facilitative interactions at x = 0 or with distant locations in niche space, modeled by negative values of the kernel, are needed to observe this phenomenon.
From a biological point of view, the self-localization indicates that species with no particular advantage may predominate to competitors. A patch of species can be formed at any position of niche space by a particular initial condition or temporary perturbation. One should note that inhomogeneities in r could increase or decrease the stability of the considered states. Although the results have been obtained in one dimensional space, and there are important differences with higher dimensional cases, we expect that the conditions for the observation of localized states will be qualitatively similar.
We acknowledge financial support from FEDER and MINECO (Spain) through grant FIS2012-30634 IN-TENSE@COSYP, and from Comunitat Autónoma de les Illes Balears. DG acknowledges support from CSIC (Spain) through grant number 201050I016.

Appendix: Numerical methods
To find the self-localized solution of (11) we write first the steady state condition ∂ ∂t = 0: introducing auxiliary quantities a, b, c equation (29) is transformed to the system of ordinary differential equations: Interpreting now x as a dynamical variable we solve the system (30) with initial conditions a = 0, b = 0, c = 0, and Ψ = Ψ 0 = 1 plus small perturbations. For parameters close to the subcritical bifurcation, trajectories showing localized pulses as the ones shown in Figure  8 are easily found. Taking one of the chirped pulses of this figure as a initial guess we can compute the branch shown in Fig. 4 using a Newton method and continuation techniques [40].