# Spatial eco-evolutionary dynamics along environmental gradients: multi-stability and cluster dynamics

Martín Andrade-Restrepo, Nicolas Champagnat, Régis Ferrière

► To cite this version:

Martín Andrade-Restrepo, Nicolas Champagnat, Régis Ferrière. Spatial eco-evolutionary dynamics along environmental gradients: multi-stability and cluster dynamics. Ecology Letters, 2019, 22 (5), pp.767-777. hal-01732325

HAL Id: hal-01732325

https://inria.hal.science/hal-01732325

Submitted on 14 Mar 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Spatial eco-evolutionary dynamics along environmental gradients:
multi-stability and cluster dynamics

Martín Andrade-Restrepo¹, Nicolas Champagnat²,³, Régis Ferrière⁴,⁵

¹ Institut Jacques Monod, CNRS UMR 7592, Université Paris Diderot, Paris Cité Sorbonne, F-750205, Paris,
France

² IECL, CNRS UMR 7502, Université de Lorraine, Vandœuvre-lès-Nancy, F-54506, France

³ Inria, TOSCA team, Villers-lès-Nancy, F-54600, France

⁴ Département of Ecology and Evolutionary Biology, University of Arizona

⁵ Eco-Evolution Mathématique, CNRS UMR 8197, Ecole Normale Supérieure, Paris Sciences & Lettres University

**Abstract**

How the interplay of local adaptation and dispersal determines species appearance, distribution and range dynamics is still incompletely understood. Here we combine individual-based simulations and mathematical analysis of large-population approximation models to advance the analysis of spatial spread and phenotypic diversification of a single-species population along a one-dimensional resource gradient. Local competition shapes selection on heritable variation in the individual ecological trait (niche position) and the evolutionary response feeds back on the local ecological state of the population (abundance). Key parameters of spatial spread and phenotypic diversification are the individual dispersal rate, the size of the spatial competition neighborhood, and the phenotype mutational variance. From a focal location the population spreads by forming clusters in space and/or trait, or by spreading along a continuous cline in both space and trait. The conditions for clustering are broader than previously known. The spacing of clusters is determined by the spatial scale of competition. When the space-trait domain is bounded, multi-stability occurs, whereby small initial differences can lead to alternative spatial and trait distributions. The transient dynamics involve adaptational lags which cause a slow-down in cluster formation and population range expansion.

Keywords:

Spatially structured population, Adaptive evolution, Adaptive diversification, Clustering,
Polymorphism, Individual-based process

**1. Introduction**

The concept of ‘eco-evolutionary feedback’ refers to the reciprocal influence that ecological and evolutionary change can exert on one another [1, 2, 3, 4]). Such ‘eco-evolutionary feedbacks’ can be described by using three ingredients: (1) heritable traits that affect some ecological properties of the system, (2) ecological modifications that are persistent and strong enough to alter selection on the traits, and (3) an actual adaptive response of these traits to the change in selection [1, 5]. When such a closed feedback loop operates between ecological and evolutionary processes the joint trajectories of ecological and evolutionary variables describe the system’s ‘eco-evolutionary dynamics’ [2, 6].

Email address: martinandraderestrepo@gmail.com (Martín Andrade-Restrepo¹)
Mathematical models have shown how eco-evolutionary feedbacks alter population dynamics [7, 8] and viability [9], species interactions such as predator-prey, host-pathogen and mutualism [10, 11, 12], community diversity and structure [13, 14], and ecosystem function [15]. However, the spatial scale of eco-evolutionary feedbacks, and the impact of these effects on spatial population dynamics and species spatial range and distribution remain poorly understood [16, 17, 18, 19, 20].

Seminal models of spatial eco-evolutionary dynamics [21] addressed the case of a single species in-
habiting a one-dimensional environmental gradient, assuming local interactions (competition) among
individuals and local dispersal; heritable variation in niche position fuels local adaptation, which
shapes spatial variation in population abundance, closing the feedback loop on trait evolution. Ex-
tending the earlier theory of character displacement and species coexistence along one-dimensional
niche axis and its non-spatial elaborations [22, 23, 24, 25, 26], these models made the distinctive
prediction of population clustering in geographic and phenotypic space. This was a shift from the
paradigm of local adaptation along environmental gradients resulting in gradual and smooth vari-
ation in the average value of a trait (cline-like phenotypic distributions [27, 28, 29, 30]). Spatial
eco-evolutionary feedbacks driving population clustering might provide a mechanism for evolution-
ary processes of diversification such as parapatric speciation. This is supported by empirical case
studies, such as speciation of rockfish in parapatry along a gradual change in oceanic depth [31].
Other relevant examples are reported in [32]. We refer to clustering as a collective phenomenon
emerging from local interactions in which individuals form high-density groups interspersed with
low density areas [33, 34, 35, 36, 37]. Clusters stand in contrast of smooth distributions across the
space-phenotype domain [38, 39, 40, 41].

The current theory of spatial eco-evolutionary dynamics leaves basic questions unanswered re-
garding cluster formation and dynamics. First, our understanding of the conditions required for
cluster formation and persistence is still incomplete. In particular, how does individual mobility
affect cluster formation, and how does the effect of individual mobility interact with the scale over
which individuals compete for resources? How does the bounded nature of both geographic and
phenotypic space influence the dynamics of clustering? And can we predict key characteristics of
the emerging population structure, such as inter-cluster distance, from individual-level parameters?
Here we answer these questions by taking a scaling approach whereby the construction and analysis
of a ‘microscopic’ individual-based, stochastic model of spatial eco-evolutionary dynamics is supple-
mented with the analysis of a ‘macroscopic’ deterministic approximation model [42].

Our analysis reveal two important properties of cluster dynamics. First, spatial eco-evolutionary
dynamics exhibits multi-stability. Thus, even small variations in initial conditions can lead to dif-
ferent attractors, manifesting as different cluster patterns. This shows that spatial eco-evolutionary
feedbacks can drive the same invading population to different spatial and trait distributions; also,
disturbances may cause a population to switch rapidly between alternate spatial distribution and
phenotypic composition. Second, clusters do not form at a constant rate. As the population spreads
from a focal location, the wait time of peripheral cluster formation increases, which results in an
invasion slow-down. This pattern arises from an increasing ‘adaptational lag’ in marginal clusters.
The transient dynamics of local adaptation in newly formed clusters is thus critical to predict the
long-term dynamics of invasion.

Our model is similar to [42] and differs from [43] in several respects. First, rather than limiting
individual mobility to dispersal at birth, we assume that individuals can change location throughout
their life. Intuitively, this favors mixing and sets a priori more restrictive conditions for cluster formation. Second, we use competition kernels that are box-shaped rather than Gaussian. This removes spurious effects of Gaussian competition kernels on population dynamics [22, 23, 24, 25, 44, 43] and allows for an unambiguous measurement of the competition range. Third, we exclude fitness frequency-dependence from our assumptions, by assuming that competition intensity does not depend on relative trait values. Fourth, we address both cases of unbounded and bounded trait space. Fifth, we do not include Allee effects, which would otherwise extinguish local populations of extremely low density. Not including Allee effects gives us the opportunity to compare and contrast the exact dynamics of the individual-based stochastic model with the approximated dynamics of the macroscopic deterministic model.

The paper unfolds as follows. In Section 2 we present the stochastic microscopic model and its deterministic macroscopic approximation. In Section 3.1 we perform a numerical analysis of the models' dynamics. A stability analysis is presented in Section 3.2. In Section 3.3 we use a Hamilton-Jacobi approach to analyse the asymptotic pattern and transient dynamics of population spread through cluster formation. Numerical schemes used in simulations and additional numerical and mathematical analyses are given in the Supplementary Material, where we also contrast the eco-evolutionary model with its purely ecological counterpart.

## 2. Models and Methods

### 2.1. Individual-based stochastic model of phenotypic evolution

We use an individual-based stochastic model of spatial eco-evolutionary dynamics which was first introduced in [42]. Given two smooth domains (open, connected) $\mathcal{X} \subseteq \mathbb{R}^d$ (spatial domain) and $\mathcal{U} \subseteq \mathbb{R}^k$ (phenotypic domain), we consider an asexual population where each individual is characterized by its physical location $x \in \bar{\mathcal{X}}$ and its phenotypic trait $u \in \mathcal{U}$. Individuals give birth at a rate which depends on how adapted they are to their local environment; the degree of (mal)adaptation is measured by the difference between the individual trait and the local optimum. Individuals die at a rate which increases with the intensity of local competition. Offspring inherit their parent's trait, unless a mutation occurs; the mutation probability is denoted by $\gamma$. Individuals mobility is modeled as spatial diffusion reflected at the boundary of $\mathcal{X}$. Unless specified otherwise, we take $\mathcal{X} = \mathcal{U} = (0, 1)$.

More precisely, the birth rate $B$ of an individual with trait $u$ located at position $x$ is given by

$$B(x, u) = \max \{b_0 - b_1(x - u)^2; 0\},$$

where $b_0 > 0$ and $b_1 > 0$. The individual birth rate is maximal on the line $x = u$, representing the environmental gradient. The width of the region with positive birth rate is $\sqrt{b_0/b_1}$. The death rate of an individual at $(x, u)$ in a population of $N_t$ individuals at positions $x_1, \dots, x_{N_t}$ and traits $u_1, \dots, u_{N_t}$ is given by

$$d_0 + d_1 \sum_{i=1}^{N_t} 1_{|x-x_i|<\delta},$$

where $d_0$ measures the natural death rate, and $d_1$ scales the mortality effect of competition. Applying the appropriate time-scaling, we assume $d_1 = 1$ without loss of generality. It is convenient to introduce the population counting process

$$f_t = \sum_{j \neq i}^{N_t} \delta_{(x_j, u_j)}$$
which is the sum of Dirac delta function at the points where individuals are located at time $t$. The death rate can then be written as

$$D(x, u, f_t) := d_0 + \int_{\mathcal{X} \times \mathcal{U}} 1_{|x-y|<\delta} f_t(dy, dw). \quad (1)$$

The distribution of mutational effects from an individual at position $(x, u)$ is Gaussian centered at $u$ (and independent of $x$) with variance $\sigma^2$, conditioned to remain within $\mathcal{U} = (0, 1)$. The spatial diffusion coefficient (mobility rate) is assumed constant with value $D_m$. Simulations of the individual-based stochastic model simulations are implemented as explained in Appendix D.1.

## 2.2. Large population approximation of the dynamics

Taking a large-population limit on the individual-based simulation model yields the following deterministic approximation:

$$\begin{align}
\frac{\partial n(x, u, t)}{\partial t} &= D_m \frac{\partial^2 n(x, u, t)}{\partial x^2} + n(x, u, t) \times \tag{2} \\
&\quad \left( (1 - \gamma) B(x, u) - d_0 - \int_{\mathcal{X}} \int_{\mathcal{U}} 1_{|x-y|<\delta} n(y, w, t) dw dy \right) \nonumber \\
&\quad + \gamma \int_{\mathcal{U}} n(x, w, t) B(x, w) \frac{1}{\sqrt{2\pi}\sigma} \exp \left( -\frac{(u-w)^2}{2\sigma^2} \right) dw, \nonumber \\
\left. \frac{\partial n(x, u, t)}{\partial x} \right|_{x=0} &= \left. \frac{\partial n(x, u, t)}{\partial x} \right|_{x=1} = 0, \nonumber \\
n(x, 0, t) &= n(x, 1, t) = 0, \quad \forall x \in \bar{\mathcal{X}}, \quad \forall t \in [0, \infty). \nonumber
\end{align}$$

Note that the Neumann boundary condition in physical space corresponds to reflection of spatial motion at the boundary of $\mathcal{X}$, and the Dirichlet boundary condition in phenotype space indicates that traits are excluded from the boundary of $\mathcal{U}$. The approximation of the individual-based model by Equation (2). can be formally justified in the limit of large population as follows: when assuming a fixed amount of total resources, a large system composed of the order of $N$ individuals may be sustained if the biomass of each individual scales as $\frac{1}{N}$; the intensity of competition must scale as $\frac{1}{N}$ as well. Using the martingale properties of the individual-level stochastic process, [42] proved that in the limit of large $N$ the renormalized population process converges to a macroscopic deterministic limit, in which the local population density is a weak solution to Equation (2). Numerical simulations of Eq. (2) are done by using an explicit finite difference scheme, as explained in Appendix D.2.

# 3. Results

## 3.1. Conditions for the emergence of clusters vs. cilne-like distribution

Here we examine the effect of parameters $\delta$ (spatial range of competition), $D_m$ (spatial motion), and $\sigma$ (mutation range) on the spatial eco-evolutionary dynamics. We also address the influence of the initial distribution. Fig. 1 shows the formation of clusters in the stochastic individual-based model with an initial population concentrated at a single space-phenotype position. An expanding wave of clusters propagates the population across the space-phenotype domain. Here, new clusters are established by individuals that colonize competition-free areas. Even though competition is defined with respect to physical location, the interplay between spatial competition and local adaptation results in a correlation between physical location and trait, which leads to fragmentation in
both dimensions.

Simulations of the deterministic approximation model Eq. (2) give results that are remarkably consistent (Fig. 2) with the individual-based stochastic model (Fig. 1). The deterministic model clearly highlights the dichotomy between clustering patterns versus cline-like states which can be defined [43] as symmetrical distributions with respect to the line $x = u$, of the form $n(u, x, t) = \varphi(x-u)$. Comparing the stochastic and deterministic models shows that the domain boundaries and edge effects do not invalidate the cluster dynamics and population spread predicted by the deterministic model, even though deterministic diffusion creates non-zero density across the whole domain from the time of introduction. Boundary conditions may affect the long-term population state, but the transient dynamics of cluster formation and range expansion are robust to them.

Figures 3 and 4 illustrate the influence of parameter variation on cluster formation and dynamics. In all cases, snapshots shown for the longest simulation times represent the stationary state of the system. Three main effects are apparent. First (Fig. 3A, 4A), clustering may evolve in phenotypic space and not in geographic space, as a consequence of a larger mobility rate. Thus, individual mobility can spread the population out geographically, without preventing phenotypic differentiation in distinct trait clusters. In this case, the population distribution is geographically continuous, but distinct ranges of phenotypes evolve in different geographic areas.

Second (Fig. 3B, 4B), the interaction (competition) range $\delta$ is a critical determinant of population clustering. Populations with relatively short competition range evolve a cline-like pattern. Third, figures 3C, 4C show the effect of increasing the mutation range $\sigma$. Genetic variation fuels the process of local adaptation, with effects potentially conflicting with individual mobility, as alleles are exported to spatial neighborhoods where they may be poorly adapted. Larger mutational effects tend to blur phenotypic clustering without altering geographic clustering. This response is essentially opposite to the effect of increased mobility, with the difference that even very large mutation ranges may not completely offset phenotypic clustering. Considering very small mutation ranges, the stationary clustering pattern appears unaffected; only the time of cluster formation and population spread is changed, increasing as mutational variance decreases. In the limit of zero mutational variance, i.e. in the absence of genetic variation, range expansion is prevented all together, due to the inability of the species to maintain viable populations in geographic areas where optimal conditions are too different from its original niche.

Figures 5 and 6 further document the effect of large mobility rates. Large mobility rates drive fast expansion of the population across the spatial domain, in the form of two geographically broad but phenotypically narrow and diverging clusters. Each cluster is structured into a ‘hot spot’ of adaptation at the center, and ‘cold spots’ at the margins. Throughout the process of population spread, the adaptation hot spot of each cluster acts as a population source fueling the highly maladapted geographic margins. The spatial spread of each cluster is established early in the process of expansion from the site of introduction, and within each cluster, the adaptation hot spot moves along the environmental gradient from the site of introduction (early on) to asymptotic trait values that are close to, but distinct from the edges of the gradient. The long-term pattern is one of relatively uniform spatial distribution, with three (Fig. 5A, 6A) or only two (Fig. 5B, 6B) phenotypic clusters –a pattern known from models of sexually reproducing organisms, to promote parapatric speciation [21].

Finally, Figures 7 and 8 document the effect of the initial spread of the introduced population
on its spatial eco-evolutionary dynamics. In finite, relatively small populations (Fig. 7), clusters form and propagate as in the case of a narrow introduction range (Fig. 1), even for broad initial phenotypic and spatial coverage. In contrast, in large populations well approximated by the deterministic model (Fig. 8), a broad initial range can cause the population to invade as a traveling wave. Local competition and adaptation tend to create variation in density (Fig. 8A), but clusters do not separate until the population reaches the edges of the domain, where complete isolation of clusters initiate and propagates back toward the center of the gradient (Fig. 8B). Thus, with sufficiently broad initial conditions, the spatial eco-evolutionary dynamics drive a continuous range expansion, followed by cluster patterning once the environmental gradient has been fully invaded.

Simulations reported in Figs. 1, 2, 7 and 8 indicate the number of clusters evolving asymptotically may vary. This suggests that the spatial eco-evolutionary dynamics may be multi-stable, with slightly different initial conditions leading to alternate attractors. In the next section we investigate this phenomenon in more detail.

## 3.2. Multi-stability: Evidence for alternate stable clustering patterns

We use Turing's method for pattern formation analysis [45] to study the conditions under which the spatial eco-evolutionary dynamics converge to a clustering pattern or a cline-like pattern. To avoid artifactual boundary effects, we change the boundary conditions from Neumann and Dirichlet to periodic boundary conditions. In this case, our models are translation invariant both in space and trait.

The method consists in determining a cline-like invariant solution of (2), of the form $n(t, x, u) = \rho(x - u)$, and analyzing, either numerically or analytically, the stability of perturbations of the stationary solution $\rho(x - u)$. Note that, contrary to [43] where the state space is assumed unbounded, the fact that our domain is bounded imposes to consider periodic perturbations on $[0, 1]^2$, and so restricts the set of possible perturbations of the model. By the decomposition of periodic functions in Fourier series and because the model is translation invariant, it is enough to consider perturbations of the form:

$$n(0, x, u) = [1 + \varepsilon \cos(2\pi m_1 u + 2\pi m_2 x)]\rho(x - u) \quad (3)$$

for small $\varepsilon > 0$ and for any nonnegative integers $m_1, m_2$, usually called *frequencies*. If the perturbation does not grow for any couple of integers $(m_1, m_2)$, then the cline-like solution is stable. If, for a certain $(m_1, m_2)$, the perturbation grows, then one expects the attracting state of the system (if it exists) to be distributed among approximately $m_1 + m_2 + 1$ clusters, and their distance along the line $x = u$ can be expected to be close to multiples of $1/(m_1 + m_2)$, where $(m_1, m_2) \neq (0, 0)$ are the frequencies with higher growth rate of the perturbation (Fig. 9B and Fig. 9C). Hence, our analysis allows to characterize cases where the population stabilizes at cline-like or clustered equilibria and also provides an estimate of the number of clusters.

Note that the analysis of [43] was performed only for such perturbations with $m_1 = 0$ (Fig. 9C) but for any real value of $m_2$ due to the lack of boundary conditions. Hence, their stability analysis may miss perturbations acting both in space and trait directions, which could make the cline-like solution unstable in cases where they would predict it to be stable. Furthermore, we also extend there analysis by estimating as well the expected and the viable number of clusters.

Our analysis, detailed in Appendix C.1 shows the existence of the cline-like invariant solution in the torus. We also show that the Lyapunov exponent (i.e. the rate of exponential growth or decay of the perturbation) of a perturbation of the form (3) can be approximated by:
$$ \lambda_{m_1, m_2} = -4\pi^2(D_m m_2^2 + D_\gamma m_1^2) - \frac{b_0 - d_0 - \sqrt{b_1(D_m + D_\gamma)}}{2\pi\delta(m_1 + m_2)} \exp\left(-2\pi^2 m_1^2 \sqrt{\frac{D_m + D_\gamma}{b_1}}\right) \sin(2\pi\delta(m_1 + m_2)), \quad (4) $$

where $D_\gamma = \frac{b_0\gamma\sigma^2}{2}$. This equation should provide a decomposition of the dynamics in terms of the parameter space.

In the case of perturbation in the space direction only ($m_1 = 0$), the Lyapunov exponent takes the form:

$$ \lambda_{0,m_2} = -4\pi^2 D_m m_2^2 - \frac{b_0 - d_0 - \sqrt{b_1(D_m + D_\gamma)}}{2\pi\delta m_2} \sin(2\pi\delta m_2). \quad (5) $$

In Fig. 9A, we show the dependence of $\lambda_{m_1,m_2}$ on the values of $m_1$ and $m_2$ in $\mathbb{Z}$, for the same choice of parameters as in Fig. 1 and Fig.2. Instability of the cline-like solution can be predicted as a consequence of the existence of positive values of $\lambda_{m_1,m_2}$ for multiple pairs ($m_1, m_2$) of perturbations in the form of equation (3). Larger ranges of values of $m_1$ and $m_2$ are unnecessary since $\lambda_{m_1,m_2}$ decreases when $m_1$ and $m_2$ become too large due to the boundness of the exponential and the sine and the unboundness of the first term in equation (4). In this case, the highest values of $\lambda_{m_1,m_2}$ occur for $|m_1 + m_2| = 7$, thus suggesting an attracting distribution for the original system fragmented into 8 clusters at a distance 1/7 from each other, as observed in Fig. 1. In the torus only 7 clusters form since the ones at the boundary merge into a single one. Nevertheless, $\lambda_{m_1,m_2}$ can also be positive for $|m_1 + m_2| = 6$ and $|m_1 + m_2| = 8$ suggesting that a distribution fragmented into 7 or 9 clusters is also possible when varying characteristics not considered here (for instance, the initial coverage of the population). However, these distributions grow slower for a population close to the cline-like equilibrium. Furthermore, the maxima of $\lambda_{m_1,m_2}$ occur for non-zero values of $m_1$ at the points $(m_1, m_2) = (2, 5)$ and $(m_1, m_2) = (-2, -5)$, thus showing that the effect is strongest when periodic perturbations are made as in Fig. 9B instead of as in Fig. 9C. This is corroborated in Fig. 10A and Fig. 10B, where values of $D_m$ and $\delta$ exist for which a perturbation only in the spatial direction would not predict instability on the cline-like equilibrium. Indeed, between $D_m \approx 1 \times 10^{-4}$ and $D_m \approx 1.6 \times 10^{-4}$ and in the small range between $\delta \approx 0.076$ and $\delta \approx 0.08$, only perturbations with $(m_1, m_2) = (2, 5)$ and $(m_1, m_2) = (-2, -5)$ would predict instability. In Appendix B, we investigate this region and compare the analysis in [43] with ours for $D_m = 1.3 \times 10^{-4}$. In addition, we show how in absence of evolution clustering (although only spatial) would occur under smaller ranges of parameters.

From Fig. 10 we can predict that the cline-like equilibrium distribution becomes stable for $D_m > D_m^* \approx 1.6 \times 10^{-4}$ when all other parameters remain at their default values (Table S1). Therefore, the transition from the attractor observed in Fig. 2 to the one in Fig. 4A must occur close to this value. Contrary to $D_m$, an increase in $\delta$ above the threshold value $\delta^* \approx 0.076$ causes the cline-like equilibrium to become unstable. It is expected that the dynamics will then converge towards a clustering pattern. We thus argue that the transition observed in Fig. 4B occurs close to this threshold value $\delta^*$.

In Appendix D.3 we give the numerical scheme used for validating the results derived in this section with respect to the exponent $\lambda_{m_1,m_2}$. The numerical growth exponent $\lambda_{num}$ obtained with
this scheme is shown in Fig. 10. Both approaches yield similar values of mobility rate and competition range at which the cline-like equilibrium loses stability. The agreement between $\lambda_{m1,m2}$ and $\lambda_{num}$ is satisfactory only for values $\lambda_{m1,m2}$ which are not too negative. As explained in Appendix D.3, this is due to the fact that the computation of $\lambda_{num}$ is very sensitive to small errors in the approximation of the cline-like stationary solution of the deterministic approximation model. However, since the agreement is strong for values of $\lambda_{m1,m2}$ close to zero, our method accurately predicts the stability of the cline-like solution.

### 3.3. Cluster dynamics and asymptotics under low mobility and small mutations

Here we present a Hamilton-Jacobi approach (following the work in [46]) to analyze the dynamics and asymptotics of clustering the PDE model for low individual mobility and small mutations. The method is based on a concentration approximation of the population as a sum of Dirac delta functions, allowing to derive a Hamilton-Jacobi equation whose solution $\varphi(x, u, t)$ is non-positive and has its zeros exactly at the points where the population is concentrated.

We take $\mathcal{X} = \mathcal{U} = \mathbb{R}$ to make competition uniform along $\mathcal{X} \times \mathcal{U}$. Let $\varepsilon > 0$. In order to concentrate the population's distribution into Dirac delta functions (by letting $\varepsilon \to 0$), we will redefine the diffusion coefficient as $D_m\varepsilon^2$ and the standard deviation of the mutation transition measure as $\varepsilon\sigma$. In this case, the appropriate time scaling to observe a limit process in the limit of small mutations and slow motion is $\tau = \frac{t}{\varepsilon}$. Hereafter we will simplify notations and use $t$ for $\tau$. If $n^\varepsilon(x, u, t)$ represents the approximation of $n(x, u, t)$ according to this rescaling then one has that

$$ \varepsilon \frac{\partial}{\partial t} n^\varepsilon(x, u, t) = \varepsilon^2 D_m \Delta_x n^\varepsilon(x, u, t) + n^\varepsilon(x, u, t) \times \left( (1 - \gamma) B(x, u) - d_0 - \int_{\mathcal{X}} \int_{\mathcal{U}} 1_{|x-y|<\delta} \times n^\varepsilon(y, w, t) dw dy \right) + \gamma \int_{\mathcal{U}} B(x, w) n^\varepsilon(x, w, t) \Pi^\varepsilon(w \to u) dw. $$

Where

$$ \Pi^\varepsilon(w \to u) = \frac{1}{\varepsilon\sigma\sqrt{2\pi}} \exp\left(-\frac{(u-w)^2}{2\sigma^2\varepsilon^2}\right) $$

denotes the rescaled mutation-transition measure concentrating mutation events in a neighborhood of size $\varepsilon$ around $w$. Multiplying on both sides of the equation by $\frac{1}{n^\varepsilon(x, u, t)}$ yields

$$ \frac{\varepsilon}{n^\varepsilon(x, u, t)} \frac{\partial}{\partial t} n^\varepsilon(x, u, t) = (1 - \gamma) B(x, u) - d_0 - \int_{\mathcal{X}} \int_{\mathcal{U}} 1_{|x-y|<\delta} \times n^\varepsilon(y, w, t) dw dy + \frac{\varepsilon^2 D_m \Delta_x n^\varepsilon(x, u, t)}{n^\varepsilon(x, u, t)} + \frac{\gamma}{n^\varepsilon(x, u, t)} \times \int_{\mathcal{U}} B(x, w) n^\varepsilon(x, w, t) \frac{1}{\varepsilon\sigma\sqrt{2\pi}} \exp\left(-\frac{(u-w)^2}{2\sigma^2\varepsilon^2}\right) dw. \quad (6) $$

Consider now the substitution $\varphi^\varepsilon(x, u, t) := \varepsilon \log(n^\varepsilon(x, u, t))$, or, expressed differently, $n^\varepsilon(x, u, t) = \exp\left(\frac{\varphi^\varepsilon(x, u, t)}{\varepsilon}\right)$. Intuitively, what we aim for is that, when setting $\varepsilon \to 0$, $\varphi^\varepsilon(x, u, t)$ will converge to a function $\varphi(x, u, t)$ which will be negative at the points where $n^\varepsilon(x, u, t)$ converges to 0, and that will never be positive as this would make $n^\varepsilon(x, u, t)$ blow up. Thus, assuming that $\int_{\mathcal{X}} \int_{\mathcal{U}} n^\varepsilon(y, w, t) dw dy$ is strictly positive and bounded from above (both conditions independent of $\varepsilon$), we will have that, as $\varepsilon \to 0$, $n^\varepsilon(x, u, t)$ will concentrate in Dirac masses at the points where $\varphi^\varepsilon(x, u, t) \to \varphi(x, u, t) = 0$.
As anticipated, making the substitution in (6) and introducing the change of variable $z = \frac{u-w}{\epsilon}$ we obtain

$$ 
\begin{aligned}
\frac{\partial}{\partial t} \varphi^{\epsilon}(x, u, t) ={}& (1 - \gamma) B(x, u) - d_0 - \int_{\mathcal{X}} \int_{\mathcal{U}} 1_{|x-y|<\delta} \times \exp\left(\frac{\varphi^{\epsilon}(y, w, t)}{\epsilon}\right) dw dy \\
& + \epsilon D_m \Delta_x \varphi^{\epsilon}(x, u, t) + D_m \left|\frac{\partial}{\partial x} \varphi^{\epsilon}(x, u, t)\right|^2 \\
& - \gamma \int_{\frac{u}{\epsilon}}^{\frac{u-1}{\epsilon}} B(x, u + z\epsilon) \exp\left(-\frac{(\varphi^{\epsilon}(x, u, t) - \varphi^{\epsilon}(x, u + z\epsilon, t))}{\epsilon}\right) \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{z^2}{2\sigma^2}\right) dz.
\end{aligned}
$$

From this equation, assuming $\epsilon$ small suggests that $\varphi^{\epsilon}(x, u, t)$ converges to a function $\varphi(x, u, t)$ satisfying the Hamilton–Jacobi equation (of first order with a Hamiltonian non-linearity in $(\frac{\partial \varphi(x, u, t)}{\partial x}, \frac{\partial \varphi(x, u, t)}{\partial u}))$:

$$ 
\begin{aligned}
\frac{\partial}{\partial t} \varphi(x, u, t) ={}& B(x, u) - d_0 - \int_{\mathcal{X}} \int_{\mathcal{U}} 1_{|x-y|<\delta} \times \mu_t(dy, dw) \\
& + D_m \left|\frac{\partial}{\partial x} \varphi(x, u, t)\right|^2 + \gamma B(x, u) H\left(\frac{\partial \varphi(x, u, t)}{\partial u}\right),
\end{aligned}
\quad (7) $$

where

$$ H\left(\frac{\partial \varphi(x, u, t)}{\partial u}\right) := \int_{-\infty}^{\infty} \left( \exp\left(-\frac{\partial \varphi(x, u, t)}{\partial u}z\right) - 1 \right) \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{z^2}{2\sigma^2}\right) dz. $$

Since the mutation kernel is a Gaussian we can simplify the function $H(\frac{\partial \varphi(x, u, t)}{\partial u})$ to be:

$$ H\left(\frac{\partial \varphi(x, u, t)}{\partial u}\right) = \exp\left(\frac{\sigma^2}{2}\left(\frac{\partial \varphi(x, u, t)}{\partial u}\right)^2\right) - 1. \quad (8) $$

The measure $\mu_t(dy, dw)$ is the weak limit of $n^{\epsilon}(x, u, t) = \exp(\frac{\varphi^{\epsilon}(y, w, t)}{\epsilon})dwdy$ when $\epsilon$ tends to zero. Since $\varphi^{\epsilon}(x, u, t) = \epsilon \log(n^{\epsilon}(x, u, t))$, $\varphi^{\epsilon}(x, u, t) \to \varphi(x, u, t)$ implies that $\mu_t(dy, dw)$ should be concentrated at the zeros of $\varphi$ at a time $t$, set which we denote by $\Omega_t$. In general, the proof of $\varphi^{\epsilon}(x, u, t) \to \varphi(x, u, t)$ and $\exp(\frac{\varphi^{\epsilon}(y, w, t)}{\epsilon})dwdy \to \mu_t(dy, dw)$ is difficult and has only been obtained for a few models [47, 48, 49, 50, 51, 52, 53]; here we will assume that convergence holds and focus on the limiting equation. Thus, we have that $\varphi(x, u, t)$ must satisfy the following conditions:

$$ (i) \quad \varphi(x, u) \le 0 \quad \forall (x, u) \in \mathcal{X} \times \mathcal{U}, $$

$$ (ii) \quad \sup_{(x,u) \in \mathcal{X} \times \mathcal{U}} \varphi(x, u, t) = 0, $$

$$ (iii) \quad \mathrm{Supp}(\mu) \subseteq \Omega_t := \{(x, u) : \varphi(x, u, t) = 0\}, $$

$$ (iv) \quad \forall (x, u) \in \Omega_t, \quad (1 - \gamma)B(x, u) - d_0 - $$

$$ \int_{\mathcal{X}} \int_{\mathcal{U}} 1_{|x-y|<\delta} \mu_t(dy, dw) \le 0, $$

$$ (v) \quad \forall (x, u) \in \mathrm{Supp}(\mu_t), \quad \{(1 - \gamma)B(x, u) - d_0 - \int_{\mathcal{X}} \int_{\mathcal{U}} 1_{|x-y|<\delta} \mu_t(dy, dw)\} = 0, $$

where (i) is needed to prevent the solution $u_{\epsilon}$ to explode as $\epsilon \to 0$, (ii) follows from the fact that $\int_{\mathcal{X} \times \mathcal{U}} u_{\epsilon}(t, x, u)dx du$ cannot vanish when $\epsilon \to 0$, and (iii) is due to the fact that $\mu_t$ cannot give mass
at points where $\varphi$ is negative. Conditions (iv), (v) are consequences of the elements of $\Omega_t$ being local extrema of $\varphi$ and indicates that $\mu_t$ is a quasi-equilibrium for the dynamics without mutation for all $t \ge 0$ (see [52]). In particular, when $\Omega_t = \{(x_i, u_i); i \in I_t\}$ with $I_t$ finite or countable, condition (iii) implies that $\mu_t(dy, dw)$ has the form

$$ \mu_t(dy, dw) = \sum_{i \in I_t} \alpha_i \delta_{(x_i, u_i)} dw dy \qquad (9) $$

for some appropriate weights $(\alpha_i)_{i \in I_t}$. Thus, the distance between two successive zeros of $\varphi$ at time $t$ can be interpreted as the distance between clusters.

In Appendix C.2 we prove that in a periodic, invariant solution $\varphi(x, u)$ of (7) which is not identically 0 along the line $x = u$, the zeros of $\varphi$ along the line $x = u$ are necessary separated by a distance $\delta/n$ for some $n \in \mathbb{N}^*$. We also show that this invariant solution is stable only if the Dirac masses are at a distance $\delta$. This suggest that, in cases where clustering occurs (as discussed in 3.2), provided individual mobility is low and mutations are small, the clusters are spaced by a distance close to $\delta$. This is consistent with the simulations of Fig. 1 and Fig. 2.

Equation (7) also provides information on the transient dynamics of clusters – on their formation times, location, motion, and shape. As discussed in [54], one has to be cautious about the conclusions regarding time scales since the speed of evolution in the Hamilton-Jacobi model is very sensitive to the initial condition away from zero. However, the location of clusters in the Hamilton-Jacobi equation brings interesting biological insights, as illustrated by the simulation of Fig. 11. In this simulation, the dynamics of $\varphi$ are studied numerically by integrating equation (7). The algorithm built for this purpose is non trivial and is described in Appendix D.4. We consider the parameter values in Table S1, which, as observed, yield clustering. For the initial population density we take $n^\varepsilon(x, u, 0) = \exp(-\frac{2(x-0.5)^2}{\varepsilon} - \frac{2(u-0.5)^2}{\varepsilon})$, i.e. an initial population close to the Dirac mass at $(0.5, 0.5)$. This choice for a gaussian distribution in the original model translates into an initial condition where $\varphi$ is quadratic and non-positive with one zero at $(0.5, 0.5)$.

Figure 11A shows snapshots of the simulation at the times where new pairs of clusters emerge (when new zeros of $\varphi$ appear). The irregular shape of the function $\varphi$ is due to the fact that the interaction kernel is not continuous. Otherwise, one would expect $\varphi$ to be smoother, except at local minima of $\varphi$, where singularities are expected in Hamilton-Jacobi equations. Fig. 11 confirms that, when the population is initially concentrated at the mid-point of the 2-dimensional domain, and the dynamics converge to clustering patterns, the expected distance between the clusters (in physical space) will be $\delta$. In addition, the rate at which zeros appear (which is equivalent to the rate at which new clusters form) is not constant: the time between two consecutive clustering events increases. Even though time has been re-scaled, properties such as the ratio of clustering times are left unaffected.

We also observe that spatial clustering and phenotypical clustering need not be simultaneous. In fact, spatial clustering occurs faster. This is due to the fact that the local speed of diffusion $D_m = 5 \times 10^{-5}$ is larger than the local maximal speed of trait motion $b_0\gamma\sigma^2/2 = 10^{-5}$, obtained by a local expansion of $\exp(-\frac{\partial\varphi}{\partial u}z)$ in the mutation term in (7) (valid since $\sigma^2$ is small). This is shown in Fig. 11B, where clusters first occur away from the diagonal, and then slowly shift towards the line $x = u$ as the local populations adapt. Because our rescaling of ?motion? was the same both in the spatial and phenotypic directions, clusters in the original system are expected to behave similarly.
# 4. Discussion

Previous studies of spatial eco-evolutionary dynamics focused either on stochastic individual-based simulations or on phenomenological deterministic population models. In this study, we combine numerical simulations of an individual-based model with the mathematical analysis of a deterministic model which was derived rigorously by appropriate rescaling of the individual-level processes [42]. Our analysis confirms the previous finding [43] that asymptotically the system evolves either to a clustered pattern or a cline-like distribution, and that clustering is not an artifact of boundary conditions. Our results shed new light on the conditions under which clustering versus cline-like patterns emerge; on the transient eco-evolutionary dynamics of population structure and expansion; and on the expected long-term distribution of clusters in geographic and trait space.

## 4.1. Model assumptions

Like previous models of spatial eco-evolutionary dynamics along an environmental gradient, our model assumes a continuous one-dimensional habitat and one trait (niche position) evolving with a single, constant optimum at any location. The smooth continuous variation of the optimum represents the environmental gradient. The trait affects the birth rate directly. Competition increases the individual death rate. In contrast with earlier models, we do not describe competition with a Gaussian kernel; rather, competition occurs uniformly within a specific range interval of constant width (competition range), and the intensity of competition is assumed independent of the competitors' traits. In other words, there is no frequency dependent selection built in the model. Negative frequency-dependence can generate disruptive selection on a trait like niche position [44, 55, 56, 57]; in our model, negative frequency-dependence emerges from competition that operates solely as a function of the geographic distance between individuals.

## 4.2. Conditions for the evolution of clustered versus cline-like distributions

Quantitative genetics model in which there are no eco-evolutionary feedbacks predict the evolution of cline-like distributions along environmental gradients. By generating local negative frequency-dependence, spatial eco-evolutionary feedbacks set the stage for the divergence of sub-populations in trait and geographic space. Albeit questioned on the basis that clustering might be an artifact of boundary conditions [58], this expectation was backed up by Leimar et al. (2008) [43] who performed a stability analysis of a deterministic model similar to ours. Our work extends their approach and reveals even broader parameter ranges for the evolution of clustering. Two parameters have a major influence: the individual mobility rate and the competition range. The spatial eco-evolutionary dynamics change from clustering to cline-like as the mobility increases through a threshold, which is detected in the individual-based stochastic simulations and predicted accurately by the deterministic macroscopic model. At even higher mobility rates, clustering re-occurs, but only in phenotype, whereas the spatial distribution is continuous. In this case, geographic areas of high local adaptation alternate with areas of low adaptation. This is the situation in which individual-based models that include sexual reproduction predict parapatric speciation [21].

## 4.3. Cluster dynamics and multi-stability

Even though our model does not include an Allee effect, the dynamics of the deterministic macroscopic approximation appear to be remarkably consistent with the individual-based simulations. Taken together, deterministic and stochastic simulations show two alternate transient regimes for evolution of clustered distributions: *cluster projection*, whereby new clusters are borne out at the edge of the current distribution, versus *front fragmentation*, whereby a continuous front expands while gaps progressively form and deepen in the core of the distribution. As with the evolution of
clusters in general, the fragmentation mode specifically is not induced by border effects; indeed the
process of gap formation can begin before the traveling front reaches the edges of the domain. One
key factor determining the transient regime is the spread of the initial focal population. Narrow pop-
ulation foci spread by *cluster projection*, whereas wide population foci spread by *front fragmentation*.
In the long run, such different initial conditions can lead to different attractors, namely clustered
distributions that differ in the number of clusters. Multi-stability is predicted by adaptive dynamics
models for well-mixed populations [11, 59, 60]; our results show that alternate stable states may also
occur as a consequence of eco-evolutionary feedbacks generated by spatially localized competition.

4.4. *Adaptational lags and invasion slow-down.*

Assuming narrow population focus and small mobility rates, we used a Hamilton-Jacobi ap-
proach [46] to gain quantitative insights in the cluster dynamics. The Hamilton-Jacobi approach
has been used both in various ecological settings and to study phenotypic diversification in well-
mixed population models [47, 48, 49, 50, 51, 52, 53]; this work provides an extension to the case of
spatially-extended populations. The analysis of the spatial Hamilton-Jacobi equation shows that the
competition range is a good approximation of inter-cluster distance. Also, cluster projection tends
to slow down as the population spread. This is due to the seeding of new clusters by individuals
that enjoy reduced competition ahead of the population edge, but suffer severe maladaptation. The
time needed for a new cluster to reach a locally adapted state appears to increase with the number
of previously formed clusters, resulting in longer adaptational lags and the invasion slow-down. This
increase in adaptational lag is due to the seeding of new clusters from previously formed leading
clusters that were themselves incompletely adapted to their local environment.

4.5. *Conclusions*

Spatially localized competition along an environmental gradient drives local adaptation and leads to the emergence of clustered versus cline-like population distributions. Clusters can form by projection (from the previously formed leading clusters) or fragmentation (behind an expanding population front). Alternate clustered distributions can be reached by a population depending on its initial distribution. When the population range expands by cluster projection, the invasion process slows down as a consequence of the increasing adaptational lag experienced by newly formed clusters. Spatial eco-evolutionary feedbacks are thus important to predict invasion dynamics and range structure in stable environments, as well as range shifts (possibly to different space-trait attractors) in response to environmental change.

**Data accessibility**

The supporting information, methods and results of this article have been uploaded as part of
the supplementary material.

Authors’ contributions

M.A.R, N.C and R.F designed and supervised the study. M.A.R performed the numerical sim-
ulations. M.A.R and N.C performed the mathematical analysis and wrote initial drafts of the
manuscript. M.A.R, N.C and R.F wrote the paper.

Competing interests

We declare we have no competing interests.
Acknowledgments

Martin Andrade Restrepo would like to thank Khashayar Pakdaman, the CNRS, and the Ad-
ministrative Department of Science, Technology and Innovation of Colombia *Colciencias*.

Nicolas Champagnat benefited from the support of the Chair “Modélisation Mathématique et Biodiversité” of Veolia Environnement - École Polytechnique - Muséum National d’Histoire Naturelle - Fondation X and from the support of the French National Research Agency for the project ANR-14-CE25-0013, “Propagation phenomena and nonlocal equations” (ANR NONLOCAL).

The authors would like to thank Feng Gao as well for the helpful discussions and his comments on an initial version of the manuscript.

References

[1] Ferrière R, Dieckmann U, Couvet D. Evolutionary conservation biology. vol. 4. Cambridge University Press; 2004.

[2] Fussmann G, Loreau M, Abrams P. Eco-evolutionary dynamics of communities and ecosystems. Functional Ecology. 2007;21(3):465-477.

[3] Schoener TW. The newest synthesis: understanding the interplay of evolutionary and ecological dynamics. science. 2011;331(6016):426-429.

[4] Reznick DN. A critical look at reciprocity in ecology and evolution: introduction to the symposium. The American Naturalist. 2013;181(S1):S1-S8.

[5] Kokko H, López-Sepulcre A. The ecogenetic link between demography and evolution: can we bridge the gap between theory and data? Ecology Letters. 2007;10(9):773-782.

[6] Smallegange IM, Coulson T. Towards a general, population-level understanding of eco-evolutionary change. Trends in ecology & evolution. 2013;28(3):143-148.

[7] Ferriere R, Gatto M. Chaotic population dynamics can result from natural selection. In: Proc. R. Soc. Lond. B. vol. 251. The Royal Society; 1993. p. 33-38.

[8] Dercole F, Ferrière R, Rinaldi S. Ecological bistability and evolutionary reversals under asymmetrical competition. Evolution. 2002;56(6):1081-1090.

[9] Ferriere R, Legendre S. Eco-evolutionary feedbacks, adaptive dynamics and evolutionary rescue theory. Phil Trans R Soc B. 2013;368(1610):20120081.

[10] Ferriere R, Bronstein JL, Rinaldi S, Law R, Gauduchon M. Cheating and the evolutionary stability of mutualisms. Proceedings of the Royal Society of London B: Biological Sciences. 2002;269(1493):773-780.

[11] Dercole F, Ferriere R, Rinaldi S. Chaotic Red Queen coevolution in three-species food chains. Proceedings of the Royal Society of London B: Biological Sciences. 2010;p. rspb20100209.

[12] Carval D, Ferriere R. A unified model for the coevolution of resistance, tolerance, and virulence. Evolution. 2010;64(10):2988-3009.

[13] Stegen JC, Enquist BJ, Ferriere R. Advancing the metabolic theory of biodiversity. Ecology letters. 2009;12(10):1001-1015.
[14] Stegen JC, Enquist BJ, Ferriere R. Eco-evolutionary community dynamics: covariation between diversity and invasibility across temperature gradients. The American Naturalist. 2012;180(4):E110-E126.

[15] Loeuille N, Loreau M, Ferrière R. Consequences of plant-herbivore coevolution on the dynamics and functioning of ecosystems. Journal of theoretical biology. 2002;217(3):369-381.

[16] Hanski I. Eco-evolutionary dynamics in a changing world. Annals of the New York Academy of Sciences. 2012;1249(1):1-17.

[17] De Meester L, Pantel J. Eco-evolutionary dynamics in freshwater systems. Journal of Limnology. 2014;73(s1).

[18] Bocedi G, Palmer SC, Pe'er G, Heikkinen RK, Matsinos YG, Watts K, et al. RangeShifter: a platform for modelling spatial eco-evolutionary dynamics and species' responses to environmental changes. Methods in Ecology and Evolution. 2014;5(4):388-396.

[19] Tack AJ, Laine AL. Spatial eco-evolutionary feedback in plant-pathogen interactions. European journal of plant pathology. 2014;138(3):667-677.

[20] You L, Brown JS, Thuijsman F, Cunningham JJ, Gatenby RA, Zhang J, et al. Spatial vs. non-spatial eco-evolutionary dynamics in a tumor growth model. Journal of theoretical biology. 2017;435:78-97.

[21] Doebeli M, Dieckmann U. Speciation along environmental gradients. Nature. 2003;421(6920):259-264.

[22] Roughgarden J. Evolution of niche width. American Naturalist. 1972;p. 683-718.

[23] MacArthur RH. Geographical ecology: patterns in the distribution of species. Princeton University Press; 1972.

[24] May RM. Stability and complexity in model ecosystems. vol. 6. Princeton University Press; 1973.

[25] Roughgarden J. Theory of population genetics and evolutionary ecology: an introduction. 1979;.

[26] Scheffer M, van Nes EH. Self-organized similarity, the evolutionary emergence of groups of similar species. Proceedings of the National Academy of Sciences. 2006;103(16):6230-6235.

[27] Slatkin M. Spatial patterns in the distributions of polygenic characters. Journal of Theoretical Biology. 1978;70(2):213-228.

[28] Kirkpatrick M, Barton NH. Evolution of a species' range. The American Naturalist. 1997;150(1):1-23.

[29] Barton NH. Clines in polygenic traits. Genetic research. 1999;74(03):223-236.

[30] Case TJ, Taper ML, et al. Interspecific competition, environmental gradients, gene flow, and the coevolution of species' borders. The American Naturalist. 2000;155(5):583-605.

[31] Ingram T. Speciation along a depth gradient in a marine adaptive radiation. Proceedings of the Royal Society of London B: Biological Sciences. 2011;278(1705):613-618.
[32] Keller I, Seehausen O. Thermal adaptation and ecological speciation. Molecular ecology. 2012;21(4):782-799.

[33] Hassell M, May R. Aggregation of predators and insect parasites and its effect on stability. The Journal of Animal Ecology. 1974;p. 567-594.

[34] Hassell MP, Pacala SW, Tingley G. Heterogeneity and the dynamics of host-parasitoid interactions [and discussion]. Philosophical Transactions of the Royal Society of London B: Biological Sciences. 1990;330(1257):203-220.

[35] Niwa HS. Self-organizing dynamic model of fish schooling. Journal of theoretical Biology. 1994;171(2):123-136.

[36] Flierl G, Grünbaum D, Levins S, Olson D. From individuals to aggregations: the interplay between behavior and physics. Journal of Theoretical biology. 1999;196(4):397-454.

[37] Young W, Roberts A, Stuhne G. Reproductive pair correlations and the clustering of organisms. Nature. 2001;412(6844):328-331.

[38] Haldane J. The theory of a cline. Journal of genetics. 1948;48(3):277-284.

[39] Fisher RA. Gene frequencies in a cline determined by selection and diffusion. Biometrics. 1950;6(4):353-361.

[40] Bazykin A. Hypothetical mechanism of speciation. Evolution. 1969;23(4):685-687.

[41] Endler JA. Geographic variation, speciation, and clines. 10. Princeton University Press; 1977.

[42] Champagnat N, Méléard S. Invasion and adaptive evolution for individual-based spatially structured populations. Journal of Mathematical Biology. 2007;55(2):147-188.

[43] Leimar O, Doebeli M, Dieckmann U. Evolution of phenotypic clusters through competition and local adaptation along an environmental gradient. Evolution. 2008;62(4):807-822.

[44] Metz JA, Geritz SA, Meszéna G, Jacobs FJ, Van Heerwaarden JS. Adaptive dynamics: a geometrical study of the consequences of nearly faithful reproduction. 1995;

[45] Turing AM. The chemical basis of morphogenesis. Philosophical Transactions of the Royal Society of London B: Biological Sciences. 1952;237(641):37-72.

[46] Diekmann O, Jabin PE, Mischler S, Perthame B. The dynamics of adaptation: an illuminating example and a Hamilton-Jacobi approach. Theoretical population biology. 2005;67(4):257-271.

[47] Barles G, Perthame B. Concentrations and constrained Hamilton-Jacobi equations arising in adaptive dynamics. In: Recent developments in nonlinear partial differential equations. vol. 439 of Contemp. Math. Providence, RI: Amer. Math. Soc.; 2007. p. 57-68.

[48] Perthame B, Barles G. Dirac concentrations in Lotka-Volterra parabolic PDEs. Indiana Univ Math J. 2008;57(7):3275-3301. Available from: http://dx.doi.org/10.1512/iumj.2008.57.3398.

[49] Carrillo JA, Cuadrado S, Perthame B. Adaptive dynamics via Hamilton-Jacobi approach and entropy methods for a juvenile-adult model. Mathematical Biosciences. 2007;205(1):137-161.
[50] Barles G, Mirrahimi S, Perthame B, et al. Concentration in Lotka-Volterra parabolic or integral equations: a general convergence result. Methods and Applications of Analysis. 2009;16(3):321-340.

[51] Lorz A, Mirrahimi S, Perthame B. Dirac mass dynamics in multidimensional nonlocal parabolic equations. Communications in Partial Differential Equations. 2011;36(6):1071-1098.

[52] Champagnat N, Jabin PE. The evolutionary limit for models of populations interacting competitively via several resources. Journal of Differential Equations. 2011;251(1):176-195.

[53] Mirrahimi S, Perthame B, Wakano JY. Evolution of species trait through resource competition. Journal of mathematical biology. 2012;64(7):1189-1223.

[54] Perthame B, Gauduchon M. Survival thresholds and mortality rates in adaptive dynamics: conciliating deterministic and stochastic simulations. Math Med Biol. 2010;27(3):195-210.

[55] Geritz SA, Metz JA, Kisdi É, Meszéna G. Dynamics of adaptation and evolutionary branching. Physical Review Letters. 1997;78(10):2024.

[56] Geritz SA, Mesze G, Metz JA, et al. Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree. Evolutionary ecology. 1998;12(1):35-57.

[57] Dieckmann U, Doebeli M. On the origin of species by sympatric speciation. Nature. 1999;400(6742):354-357.

[58] Polechová J, Barton NH. Speciation through competition: a critical review. Evolution. 2005;59(6):1194-1210.

[59] Geritz SA, van der Meijden E, Metz JA. Evolutionary dynamics of seed size and seedling competitive ability. Theoretical population biology. 1999;55(3):324-343.

[60] Kisdi É, Geritz SA. Adaptive dynamics in allele space: evolution of genetic polymorphism by small mutations in a heterogeneous environment. Evolution. 1999;53(4):993-1008.
<figure><img src="image_2.png" /><figcaption>**Figure 1:** *Spatial eco-evolutionary dynamics in the individual-based stochastic model. Initially, a population of N = 3000 is concentrated at the point (0.5, 0.5). Parameters set to the default clustering values (Table S1).*</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2:** *Spatial eco-evolutionary dynamics in the deterministic large population limit. At T = 0 the distribution is a Gaussian with standard deviation σ₀ = 0.1 (in both dimensions) centered at (0.5, 0.5). Parameters set to the default clustering values (Table S1).*</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3:** *Effect of mobility, competition, and mutation in the individual-based stochastic model.* A. Increased mobility rate: $D_m = 5 \times 10^{-4}$. B. Reduced interaction range: $\delta = 0.05$. C. Increased mutation variance $\sigma = 0.01$. Initially, a population of $N = 3000$ is concentrated at the point (0.5, 0.5). Other parameters set to the default clustering values in Table S1.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4:** *Effect of mobility, competition, and mutation in the deterministic large population-size limit.*

A. Increased diffusion rate: $D_m = 5 \times 10^{-4}$. B. Reduced interaction range: $\delta = 0.05$. C. Increased mutation variance $\sigma = 0.01$. At $T = 0$ the distribution has the form of a Gaussian with standard deviation $\sigma_0 = 0.1$ (in both dimensions) centered at (0.5, 0.5). Other parameters set to the default clustering values in Table S1.</figcaption></figure>
A

<figure><img src="image_6.png" /></figure>

B

<figure><img src="image_7.png" /><figcaption>**Figure 5:** Effect of mobility in the individual-based stochastic model. A. $D_m = 5 \times 10^{-3}$. B. $D_m = 2 \times 10^{-2}$. Initially, a population of $N = 3000$ is concentrated at the point $(0.5, 0.5)$. Other parameters set to the default clustering values in Table S1.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6:** Effect of mobility in the deterministic large population limit. A. $D_m = 5 \times 10^{-3}$. B. $D_m = 2 \times 10^{-2}$. At $T = 0$ the distribution is a Gaussian with standard deviation $\sigma_0 = 0.1$ (in both dimensions) centered at (0.5, 0.5). Other parameters set to the default clustering values (Table S1).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7:** Effect of initial conditions on cluster formation and dynamics in the individual-based stochastic model. Initially, a population of $N = 3000$ is uniformly distributed in the range [0.25, 0.75] (panel A) or in the range [0.1, 0.9] (panel B) in both dimensions. Parameters set to the default clustering values (Table S1).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8:** Cluster formation in the deterministic large population limit. At $T = 0$ the distribution has the form of a Gaussian with standard deviations $\sigma_0 = 0.04$ (panel A) and $\sigma_0 = 0.1$ (panel B) in both dimensions centered at (0.5, 0.5). Parameters set to the default clustering values (Table S1).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9:** A. The exponent $\lambda_{m_1,m_2}$ in equation (4) as a function of the integer frequencies $m_1$ and $m_2$. The black dots denote the maxima: $(m_1, m_2) = (2, 5)$ and $(m_1, m_2) = (-2, -5)$. Parameter values set to the default clustering parameter values in Table S1. B. Dashed line: maxima of a perturbation in the form of equation (3) of $\rho(x-u)$ on the torus $[0,1]^2$ with $m_1 = 2$ and $m_2 = 5$. Black line: $line x = u$. Black points denote the peaks of $n(0, x, u)$. C. Dashed line: maxima of a perturbation in the form of equation (3) of $\rho(x-u)$ on the torus $[0,1]^2$ with $m_1 = 0$ and $m_2 = 7$. Black line: $line x = u$. Black points denote the peaks of $n(0, x, u)$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10:** Bifurcations between clustering and cline-like stable solutions. The growth exponents $\lambda_{m_1,m_2}$ (continuous line) and $\lambda_{num}$ (dashed line) are shown as functions of $D_m$ (panel A) and $\delta$ (panel B) for perturbations with $(m_1, m_2) = (0,7)$ (red) and $(m_1, m_2) = (2,5)$ (blue). In the computation of $\lambda_{num}$, $\varepsilon = 0.001$. Other parameter values set to the default clustering parameter values in Table S1.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11:** Simulations of the Hamilton-Jacobi equation. A. Snapshots of two different profiles of $\varphi(x, u, T_i)$, $i = 1, 2, 3$ at the times where $\varphi$ develops a new pair of zeros (the times where new clusters emerge). Top: $\varphi(x, u, T_i)$ along the line $x = u$. Bottom: Contour plot of $\varphi(x, u, T_i)$ and position of its local maxima in yellow. B. Time-evolution of the angle (in degrees) between each pair of clusters and the midpoint of the domain. Parameters set to the default clustering values (Table S1).</figcaption></figure>
# Appendix

## Appendix A. Model parameters and values

The parameters of the model and their description, together with the default clustering values are presented in Table S1.

<table><caption>**Table S1: Model parameters, description and default clustering values.**</caption><thead><tr><td>Parameter</td><td>Description</td><td>Default value</td></tr></thead><tbody><tr><td>b<sub>0</sub></td><td>Maximal birth rate</td><td>b<sub>0</sub> = 2</td></tr><tr><td>b<sub>1</sub></td><td>Quadratic coefficient in the rate of decay in the birth rate</td><td>b<sub>1</sub> = 20</td></tr><tr><td>d<sub>0</sub></td><td>Natural death rate</td><td>d<sub>0</sub> = 1</td></tr><tr><td>D<sub>m</sub></td><td>Spatial diffusion coefficient</td><td>D<sub>m</sub> = 5 × 10<sup>-5</sup></td></tr><tr><td>δ</td><td>Spatial competition range</td><td>δ = 0.1</td></tr><tr><td>σ</td><td>Standard deviation of mutation transition measure</td><td>σ = 0.01</td></tr><tr><td>γ</td><td>Mutation probability</td><td>γ = 0.1</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>**Figure S1:** Cluster formation in the simulations of the individual-based stochastic model in the absence of evolution ($\gamma = 0$). Initially, a population of $N = 3000$ is concentrated at the point (0.5, 0.5). Parameter values set to the default clustering parameter values in Table S1.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure S2:** Cluster formation in the simulations of the deterministic large population-size limit in the absence of evolution ($\gamma = 0$). At $T = 0$ the distribution has the form of a Gaussian with standard deviation $\sigma_0 = 0.1$ (in the spatial direction) centered at (0.5, 0.5). Parameter values set to the default clustering parameter values in Table S1.</figcaption></figure>

# Appendix B. Dynamics of the system in the absence of evolution

## Appendix B.1. Clustering in the individual-based stochastic model and in the deterministic large-population limit in the absence of evolution

We analyse the dynamics of the individual-based stochastic model and its deterministic large population-size limit under the absence of evolution to study the dependence of clustering and adaptive diversification on phenotypic mutations when competition is only spatial. As observed in Figs. S1 and S2, clustering occurs—in both cases—in the form of a multimodal distribution under the default clustering parameter values. Complete separation of clusters does not occur however due to the adaptation constraints.

## Appendix B.2. Effect of parameters in the individual-based stochastic model and in the deterministic large-population limit in the absence of evolution

As observed in Figs. S3 and S4 for both the stochastic IB model and its deterministic limit and similarly as in the original model, both an increase in the diffusion rate and a reduction in the
A

<figure><img src="image_16.png" /></figure>

B

<figure><img src="image_17.png" /><figcaption>**Figure S3:** Parameter dependence in the simulations of the individual-based stochastic model in absence of evolution ($\gamma = 0$). A. Increased diffusion rate: $D_m = 5 \times 10^{-4}$. B. Reduced interaction range: $\delta = 0.05$. Initially, a population of $N = 3000$ is concentrated at the point (0.5, 0.5). Other parameter values set to the default clustering parameter values in Table S1.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure S4:** Parameter dependence in the simulations of the deterministic large population-size limit in absence of evolution ($\gamma = 0$). A. Increased diffusion rate: $D_m = 5 \times 10^{-4}$. B. Reduced interaction range: $\delta = 0.05$. At $T = 0$ the distribution has the form of a Gaussian with standard deviation $\sigma_0 = 0.1$ (in the spatial dimension) centered at $(0.5, 0.5)$. Other parameter values set to the default clustering parameter values in Table S1</figcaption></figure>

interaction range hinder clustering and multimodality in absence of evolution. In both cases there is a transition to unimodality (with respect to Figs. S1 and S2) with a larger range in the case of faster spatial motion.

## Appendix B.3. Dependence of clustering on phenotypic evolution in the individual-based stochastic model and in the deterministic large-population limit

To further emphasize the role of phenotypic evolution in the appearance of clustering and in the stability analysis of the cline-like equilibrium solution, we show an example where phenotypic mutations are essential for the appearance and evolution of clustering.

For a value of $D_m$ in the interval emphasized in Fig. 10A, where perturbations only at the spatial component would predict stability of the cline-like equilibrium while perturbations in both dimensions would instead predict instability, clustering occurs in presence but not in absence of phenotypic evolution (Figs. S5 and S6). In the case of the individual-based stochastic model, although unimodality is not as stable as in Fig. S3, the distribution does not show a clear multimodality as in Fig. S1.

In the deterministic large population-size limit (Fig. S6), the difference is more visible. While clustering occurs in presence of evolution (although the convergence time is much larger than in Fig. 2 and the dynamics are dependent on the boundary conditions); clustering does not occur in absence
A

<figure><img src="image_19.png" /></figure>

T=5

<figure><img src="image_20.png" /></figure>

T=90

<figure><img src="image_21.png" /></figure>

T=150

<figure><img src="image_22.png" /></figure>

T=450

B

<figure><img src="image_23.png" /></figure>

T=5

<figure><img src="image_24.png" /></figure>

T=90

<figure><img src="image_25.png" /></figure>

T=150

<figure><img src="image_26.png" /></figure>

T=450

<caption>**Figure S5:** Clustering in presence of evolution (panel A) and no clustering in the absence of evolution (panel B with $\gamma = 0$) in the individual-based stochastic model. In both cases $D_m = 1.35 \times 10^{-4}$. Other parameter values set to the default clustering parameter values in Table S1.</caption>

of phenotypic mutations and the unimodality of the distribution is stable.

These results show that multimodality in the population's distribution and clustering is facilitated and occurs under a wider range of spatially dependent parameters when in presence of phenotypic evolution even under the assumption of only spatial competition. Indeed, phenotypic evolution becomes essential in some cases as a consequence of the interplay and the evolving correlation between selection (which is spatial and phenotype dependent) and competition.
<figure><img src="image_27.png" /><figcaption>**Figure S6:** Clustering in presence of evolution (panel A) and no clustering in the absence of evolution (panel B with $\gamma = 0$) in the deterministic large population-size limit. In both cases $D_m = 1.35 \times 10^{-4}$. Other parameter values set to the default clustering parameter values in Table S1.</figcaption></figure>
# Appendix C. Analytical Complements

## Appendix C.1. Analysis of Pattern Formation

In the infinite population-size limit, the distribution of the population in the model behaves according to

$$ 
\begin{aligned}
\frac{\partial n(x, u, t)}{\partial t} &= n(x, u, t) \times \\
& \quad \left( (1 - \gamma)B(x, u) - d_0 - \int_X \int_U \mathbb{1}_{|x-y|<\delta} \times n(y, w, t) dw dy \right) + \\
& \quad D_m \frac{\partial^2 n(x, u, t)}{\partial x^2} + \gamma \int_U n(x, w, t) B(x, w) \frac{1}{\sqrt{2\pi}\sigma} \exp\left(-\frac{(u-w)^2}{2\sigma^2}\right) dw, \\
\left. \frac{\partial n(x, u, t)}{\partial x} \right|_{x=0} &= \left. \frac{\partial n(x, u, t)}{\partial x} \right|_{x=1} = 0, \\
n(x, 0, t) &= n(x, 1, t) = 0, \quad \forall x \in \mathcal{X}, \quad \forall t \in [0, \infty),
\end{aligned}
$$

where

$$ B(x, u) = \max\{b_0 - b_1(x - u)^2, 0\}. $$

Our first goal is to construct an approximate model for which explicit computations are possible. If we consider $\sigma^2$ to be small, and $(x, u)$ to be near the line $x = u$ where the birth rate is $b_0$, by making a second order expansion of the function

$$ M(u) = \gamma \left( \int_{\mathcal{U}} n(x, w, t) B(x, w) \frac{1}{\sqrt{2\pi}\sigma} \exp\left(-\frac{(u-w)^2}{2\sigma^2}\right) dw - B(x, u) n(x, u, t) \right) $$

while assuming that mutations occur at a constant rate (independent of $(x, u)$), we can approximate the function $M(u)$ with a diffusion equation with rate $D_\gamma = \frac{b_0\gamma\sigma^2}{2}$. Thus we can replace equation (2) by the simpler reaction-diffusion equation:

$$ 
\begin{aligned}
\frac{\partial n(x, u, t)}{\partial t} &= n(x, u, t) \times \\
& \quad \left( B(x, u) - d_0 - \int_X \int_U \mathbb{1}_{|x-y|<\delta} \times n(y, w, t) dw dy \right) + \\
& \quad D_m \frac{\partial^2 n(x, u, t)}{\partial x^2} + D_\gamma \frac{\partial^2 n(x, u, t)}{\partial u^2}.
\end{aligned}
\quad (S1) $$

Since we are interested in studying time-constant cline-like solutions, we will modify the boundary conditions in order to make competition uniform along the line $x = u$. Let $\bar{\mathcal{X}} = \mathcal{U} = T^1$ and remove the Neumann and Dirichlet boundary conditions. Here we denote by $T^1$ the one-dimensional torus of length 1. From this point on, all operations on $x$ or $u$ will be defined on $T^1$. Hence, we will impose the condition $\delta < 0.5$ and redefine $B(x, u)$ as

$$ B(x, u) = \max\{b_0 - b_1 \times d(x, u)^2, 0\} $$

where $d(x, u) := \min\{|x - u|, 1 - |x - u|\}$ is the distance function in $T^1$.

Suppose then that $n(x, u, t)$ is a cline-like distribution, that is of the form $n(\beta, t)$, $\beta := d(x, u) \in [0, \frac{1}{2}]$. We wish to write equation (S1) in terms of the new parameter $\beta$. We proceed by analysing
the equation term by term. First, $B(x, u)$ can be substituted by $B(\beta)$ by simply redefining (with a slight abuse of notation) $B: [0, \frac{1}{2}] \to [0, b_0]$ as $B(\beta) = \max\{b_0 - b_1(\beta)^2, 0\}$. Consider now the term in equation (S1) representing death by competition. Let $D_c(x, n) := \int_X \int_U 1_{d(x-y)<\delta} \times n(y, w, t) dw dy$. We want to show that if $n(x, u, t)$ is of the form $n(\beta, t)$, $D_c(x, n)$ is independent of $x$. Let $x_1 \le x_2$, $x_1, x_2 \in [0, 1)$. Because of the new boundary conditions and the symmetry of $n$ with respect to $x$ and $u$, one has that $n(x_2, u, t) = n(x_1 + \alpha, u, t) = n(x_1, u - \alpha, t)$ for $\alpha = x_2 - x_1$. Thus $\int_U n(x_2, u, t) du = \int_U n(x_1, w, t) dw$. Hence, $\int_U n(y, w, t) dw$ is constant for all $y \in [0, 1)$. Denoting this constant value by $C$ we can deduce that

$$ 
\begin{aligned}
D_c(x, n) &= \int_X \int_U 1_{|x-y|<\delta} \times n(y, w, t) dw dy \\
&= \left| \int_{x-\delta}^{x+\delta} 1_{|x-y|<\delta} \times C dy \right| = 2\delta C,
\end{aligned}
$$

which has no dependence on $x$. Finally, for the diffusion terms in equation (S1), one obtains that $\frac{\partial^2 n(x, u, t)}{\partial x^2} = \frac{\partial^2 n(x, u, t)}{\partial u^2} = \frac{\partial^2 n(x, u, t)}{\partial \beta^2}$ for a cline-like solution $n(x, u, t) = n(\beta, t)$.

Hence, for a cline-like distribution $n(\beta, t)$, observing that $\int_{T_1} n(\beta) d\beta = 2 \int_0^{1/2} n(\beta) d\beta$, for all $\beta \in [0, 1/2]$,

$$ \frac{\partial n(\beta, t)}{\partial t} = n(\beta, t) \left( B(\beta) - d_0 - 4\delta \int_0^{1/2} n(a, t) da \right) + (D_m + D_\gamma) \frac{\partial^2 n(\beta, t)}{\partial \beta^2}. $$

Since we are interested in proving the existence of a cline-like stationary solution, we wish to study the solvability of the second order differential equation given by

$$ n(\beta) \left( B(\beta) - d_0 - 4\delta \int_0^{1/2} n(a) da \right) + (D_m + D_\gamma) n''(\beta) = 0 \quad (S2) $$

with appropriate boundary conditions: due to the form of the birth rate function and the symmetry of the torus, $n(\beta)$ must be minimum at $\beta = \frac{1}{2}$ and a maximum at $\beta = 0$, thus we have the boundary conditions $n'(0) = n'(\frac{1}{2}) = 0$.

Let $\rho = d_0 + 4\delta \int_0^{1/2} n(a) da$ (an unknown quantity) and $M = D_m + D_\gamma$. We can rewrite equation (S2) as

$$ Mn''(\beta) + n(\beta)B(\beta) = +\rho n(\beta), \quad \forall \beta \in [0, 1/2]. \quad (S3) $$

With the boundary conditions $n'(0) = n'(1/2) = 0$, this is a Sturm-Liouville (S-L) regular problem [61, 62, 63], which means that it is solvable for countably many values of $\rho$, $\rho_0 < \rho_1 < \dots$ and that the linear space of solutions associated to $\rho_n$ has dimension 1 and is generated by a function $\psi_n$ having exactly $n$ zeroes on $[0, 1/2]$. In particular, we can assume without loss of generality that $\psi_0 > 0$, and the only biologically relevant solutions are of the form $a\psi_0$ for some $a > 0$. In order to recover a solution of (S2), $a$ needs to be chosen such that $\rho_0 = d_0 + 4\delta a \int_0^{1/2} \psi_0(\beta) d\beta$. This is only possible if $\rho_0 > d_0$, hence for sufficiently small death rate. Note also that $\rho_0$ depends in a non-trivial way on the birth rate $B$ and hence on $b_0$ and $b_1$.

In order to make more precise the corresponding assumptions on $b_0$, $b_1$ and $d_0$, let us first integrate (S3) over $[0, 1/2]$. We obtain $\int_0^{1/2} \psi_0(\beta) B(\beta) d\beta = \rho_0 \int_0^{1/2} \psi_0(\beta) d\beta$, from which we deduce that $\rho_0 \in (0, b_0)$. Note also that, if the positive part in the definition of $B$ is removed, the principal eigenvalue problem (S3) on $[0, +\infty)$ with boundary conditions $n'(0) = 0$ and $n(+\infty) = 0$ becomes
explicitly solvable, with eigenvalue $\bar{\rho}_0 = b_0 - \sqrt{b_1 M}$ and eigenfunction

$$ \bar{\psi}_0(\beta) = \bar{\psi}_0(0) \exp \left( - \frac{\beta^2}{2\sqrt{M/b_1}} \right). $$

This gives a particular solution of (S3) with $\rho = \bar{\rho}_0$ and $n'(0) = 0$

$$ n(\beta) = \begin{cases} n(0) \exp \left( -\frac{\beta^2}{2\sqrt{M/b_1}} \right) & \text{if } \beta \le \sqrt{b_0/b_1} \\ n(0) \exp \left( -\frac{b_0}{2\sqrt{Mb_1}} \right) \cos \left( \sqrt{\bar{\rho}_0/M} \beta + \bar{\varphi}_0 \right) & \text{if } \sqrt{b_0/b_1} \le \beta \le 1/2, \end{cases} $$

where $\bar{\varphi}_0$ is a solution to

$$ \sin \left( \sqrt{\frac{\bar{\rho}_0 b_0}{Mb_1}} + \bar{\varphi}_0 \right) = \left( 1 - \frac{\sqrt{Mb_1}}{b_0} \right)^{-1/2}. $$

Hence, assuming $\frac{b_0}{\sqrt{Mb_1}} \gg 1$, the previous function satisfies $n'(1/2) \approx 0$ and hence is a good approximation to a solution of the Sturm-Liouville problem (S3). This function is not positive but takes very small negative values. Therefore, $\rho_0 \approx \bar{\rho}_0$ and $\psi_0 \approx \bar{\psi}_0$. This leads to the equation

$$ \begin{aligned} b_0 - \sqrt{b_1 M} &\approx \rho_0 \approx d_0 + 4\delta \int_0^{1/2} \bar{\psi}_0(\beta) d\beta \approx d_0 + 4\delta \int_0^{\infty} \bar{\psi}_0(\beta) d\beta \\ &= d_0 + 2\delta \bar{\psi}_0(0) \sqrt{2\pi \sqrt{M/b_1}}. \end{aligned} $$

Therefore, under the assumptions $d_0 < b_0 - \sqrt{Mb_1}$ and $b_0 \gg \sqrt{Mb_1}$, we obtain the following approximation of the cline-like stationary solution of (S1):

$$ n_0(\beta) \approx \frac{b_0 - d_0 - \sqrt{b_1 M}}{2\delta\sqrt{2\pi\sqrt{M/b_1}}} \exp\left(-\frac{\beta^2}{2\sqrt{M/b_1}}\right). $$

We will write this in the form:

$$ n_0(\beta) = C \frac{1}{\sqrt{2\pi\sigma_0^2}} \exp\left(-\frac{\beta^2}{2\sigma_0^2}\right), \qquad (S4) $$

with $C = \frac{b_0 - d_0 - \sqrt{b_1 M}}{2\delta}$ and $\sigma_0^2 = \sqrt{M/b_1}$. Notice that $C$ is the same as above when $M \ll b_1$.

Now, consider a perturbation of $n_0(x, u)$ of the form

$$ \hat{n}_0(x, u) = (1 + \varepsilon\mu(x, u))n_0(x, u), \qquad (S5) $$

where $\varepsilon \ll 1$ and $\mu(x, u) = \mu(m_1 u + m_2 x)$. Here, $\mu(x, u)$ is a periodic function on the torus $(T^1)^2$. This requires $\mu(m_1(u+1) + m_2x) = \mu(m_1u + m_2x) = \mu(m_1u + m_2(x+1))$ for all $x, u \in T^1$, i.e. $m_1, m_2 \in \mathbb{Z}$. Fourier analysis tells us that any periodic function on $(T^1)^2$ is a linear combination of such functions, so it is enough to study each perturbation with $m_1, m_2 \in \mathbb{Z}$. By symmetry, we can restrict to $m_1, m_2 \in \mathbb{N}$. We will study the dynamics of $\nu(x, u, t) = \hat{n}(x, u, t) - n_0(x, u)$ when $\hat{n}(x, u, 0) = \hat{n}_0(x, u)$. Furthermore, by local linearization of the PDE it is enough to look at solutions of the form $\nu(x, u, t) = \exp(\lambda t)\nu(x, u, 0) = \exp(\lambda t)(\varepsilon\mu(x, u)n_0(x, u))$. Turing's stability
method consists in computing the value of $\lambda$ and looking for values that satisfy $\lambda > 0$. Since, $\frac{\partial\nu(x,u,t)}{\partial t} = \frac{\partial\hat{n}(x,u,t)}{\partial t}$, we get:

$$ \frac{\partial\nu(x, u, t)}{\partial t} = \nu(x, u, t) \times \left( B(x, u) - d_0 - \int_X \int_U \mathbb{1}_{d(x,y)<\delta} \times n_0(y, w) dw dy \right) + D_m \frac{\partial^2\nu(x, u, t)}{\partial x^2} + D_\gamma \frac{\partial^2\nu(x, u, t)}{\partial u^2} - \hat{n}(x, u, t) \int_X \int_U \mathbb{1}_{d(x,y)<\delta} \times \nu(y, w, t) dw dy. \quad (S6) $$

Hence, simplifying by $\varepsilon \exp(\lambda t)$ and considering $\varepsilon \ll 1$ we have

$$ \begin{aligned} \lambda\mu(x, u)n_0(x, u) ={}& \mu(x, u)n_0(x, u)(B(x, u) - d_0 - 2\delta C) + \\ & D_m \left( \frac{\partial^2\mu(x, u)}{\partial x^2}n_0(x, u) + 2\frac{\partial\mu(x, u)}{\partial x}\frac{\partial n_0(x, u)}{\partial x} + \mu(x, u)\frac{\partial^2 n_0(x, u)}{\partial x^2} \right) + \\ & D_\gamma \left( \frac{\partial^2\mu(x, u)}{\partial u^2}n_0(x, u) + 2\frac{\partial\mu(x, u)}{\partial u}\frac{\partial n_0(x, u)}{\partial u} + \mu(x, u)\frac{\partial^2 n_0(x, u)}{\partial u^2} \right) - \\ & \hat{n}(x, u, t) \int_X \mathbb{1}_{d(x,y)<\delta} \int_U \mu(y, w)n_0(y, w) dw dy. \end{aligned} $$

Since $\varepsilon \ll 1$, we substitute the last term in the previous equation with:

$$ n_0(x, u) \int_{x-\delta}^{x+\delta} \int_{U} \mu(y, w)n_0(y, w) dw dy $$

Setting $u = x$ ($\beta = 0$) and abusing notation by denoting $n_0(x, u)$ as $n_0(\beta)$ we obtain

$$ \lambda\mu((m_1 + m_2)x) = (D_m m_2^2 + D_\gamma m_1^2)\mu''((m_1 + m_2)x) - \int_{x-\delta}^{x+\delta} \int_U \mu(y, w)n_0(y, w) dw dy. \quad (S7) $$

Suppose $\mu(m_1 u + m_2 x) = \cos(2\pi(m_1 u + m_2 x))$. Using (S4), we can compute

$$ \int_{x-\delta}^{x+\delta} \int_{-\infty}^{\infty} \cos(2\pi(m_1 w + m_2 y)) n_0(y, w) dw dy = \frac{C}{\pi(m_1 + m_2)} \exp(-2\pi^2 \sigma_0^2 m_1^2) \sin(2\pi(m_1 + m_2)\delta) \cos(2\pi(m_1 + m_2)x). $$

Hence, $\mu(m_1 u + m_2 x) = \cos(2\pi(m_1 u + m_2 x))$ is a solution of equation (S7) if and only if

$$ \lambda = -4\pi^2(D_m m_2^2 + D_\gamma m_1^2) - \frac{C}{\pi(m_1 + m_2)} \exp(-2\pi^2 \sigma_0^2 m_1^2) \sin(2\pi\delta(m_1 + m_2)). $$

Furthermore, since the Fourier decomposition of periodic functions is unique we argue that every periodic solution of equation (S7) is a linear combination functions of the form $\mu(m_1 u + m_2 x) = \cos(2\pi(m_1 u + m_2 x))$ with $m_1, m_2 \in \mathbb{Z}$ giving identical values of $\lambda$.

Using the expressions of $C$ and $\sigma_0$, we finally obtain

$$ \lambda = -4\pi^2(D_m m_2^2 + D_\gamma m_1^2) - \frac{b_0 - d_0 - \sqrt{b_1(D_m + D_\gamma)}}{2\pi\delta(m_1 + m_2)} \exp\left(-2\pi^2 m_1^2 \sqrt{\frac{D_m + D_\gamma}{b_1}}\right) \sin(2\pi\delta(m_1 + m_2)). $$
Appendix C.2. The Hamilton–Jacobi approach for density dependent models with slow motion

Suppose we have an invariant solution $\hat{\varphi}(x, u)$ to equation (7). Since selection is strongest along the line $x = u$, we will assume that the points for which $\hat{\varphi}(x, u) = 0$ are located along this line i.e that $u_i = x_i, i \in \mathcal{I}$. We will try to determine the characteristics of this solution (the number of clusters and their spacing). To do this, let us first remove the boundary conditions ($\mathcal{X} = \mathcal{U} = \mathbb{R}$) and suppose, without loss of generality, that $\hat{\varphi}(0, 0) = 0$ (i.e that the population is concentrated in a Dirac delta at $(0, 0)$). Since we are interested in solutions exhibiting periodic clustering patterns we make the assumption of constant positive spacing between clusters (supported by simulations). Hence, $\Omega_t$ is countable and there exists $T > 0$ minimal for which $\hat{\varphi}(x, u) = \hat{\varphi}(x + T, u + T)$. Combining both conditions we have, $\Omega_t = \{(nT, nT), n \in \mathbb{Z}\}$. Notice that the value of $T$ is the distance between the Dirac delta peaks and hence the distance between clusters. Since we are interested in the behaviour along the line $x = u$, when imposing this condition, equation (7) for $\hat{\varphi}(x, x)$, which we denote by $\hat{\varphi}(x)$, becomes

$$0 = b_0 - d_0 - \sum_{n \in \mathbb{Z}} \alpha_n 1_{|x-nT|<\delta} + D_m \left| \frac{\partial \hat{\varphi}}{\partial x}(x, x) \right|^2 + \gamma b_0 H \left( \frac{\partial \hat{\varphi}}{\partial u}(x, x) \right). \quad (S8)$$

Since the points in $\Omega_t$ are local maxima, we have that, for all $n \in \mathbb{Z}$, $\frac{\partial \hat{\varphi}(x, u)}{\partial u}|_{(nT, nT)} = \frac{\partial \hat{\varphi}(x, u)}{\partial x}|_{(nT, nT)} = 0$. This last property requires some care, since actually, due to the singularity of the competition kernel, the function $\hat{\varphi}$ may not admit derivatives at its local maxima. However, we can approximate the competition kernel $1_{[-\delta, \delta]}$ by a sequence of smooth kernels $(K_k)_{k \ge 1}$ and argue that the corresponding solution $\varphi_k$ converges to $\hat{\varphi}$, and $\frac{\partial \varphi_k(x, u)}{\partial u}|_{(nT, nT)} = \frac{\partial \varphi_k(x, u)}{\partial x}|_{(nT, nT)} = 0$ for all $k \ge 1$. To make this rigorous, one would need to use stability properties of the Hamilton-Jacobi problem with constraints $(i)$ to $(v)$. However, this is not known in general since for such equations, even the problem of uniqueness is difficult [64, 65]. If we leave aside this difficulty, we obtain that, for $(x, u) \in \Omega_t$,

$$b_0 - d_0 - \sum_{n \in \mathbb{Z}} \alpha_n 1_{|x-nT|<\delta} = 0. \quad (S9)$$

Suppose now that $T > \delta$. Consider the point $(0, 0)$. Because of the condition stated above and the fact that all other $(x_i, u_i) \in \Omega_t$ satisfy $|x - x_i| > \delta$, one has that

$$b_0 - d_0 - \alpha_0 = 0.$$

Now, let $0 < \varepsilon < T - \delta$. By evaluating equation (S8) at the point $(\varepsilon, \varepsilon)$, one has that

$$\begin{aligned} 0 &= b_0 - d_0 - \alpha_0 + D_m \left| \frac{\partial \hat{\varphi}}{\partial x}(\varepsilon, \varepsilon) \right|^2 + \gamma b_0 H \left( \frac{\partial \hat{\varphi}}{\partial u}(\varepsilon, \varepsilon) \right) \\ &= D_m \left| \frac{\partial \hat{\varphi}}{\partial x}(\varepsilon, \varepsilon) \right|^2 + \gamma b_0 H \left( \frac{\partial \hat{\varphi}}{\partial u}(\varepsilon, \varepsilon) \right). \end{aligned} \quad (S10)$$

We observe that, by Jensen's inequality

$$H(p) = \mathbb{E}(\exp(pG) - 1) \ge \exp(p\mathbb{E}(G)) - 1 = 0$$

The inequality is strict unless $p = 0$. Therefore, equation (S10) implies $\frac{\partial \hat{\varphi}}{\partial u}(\varepsilon, \varepsilon) = \frac{\partial \hat{\varphi}}{\partial x}(\varepsilon, \varepsilon) = 0$. Since $\varepsilon$ was arbitrary, we have that for every $\kappa \in [0, T - \delta)$, $\frac{\partial}{\partial u} \hat{\varphi}(\kappa, \kappa) = \frac{\partial}{\partial x} \hat{\varphi}(\kappa, \kappa) = 0$. Hence, $\hat{\varphi}(\kappa, \kappa) = 0$ $\forall \kappa \in [0, T - \delta)$ which is a contradiction since $\Omega_t$ was taken to be countable.
Suppose now that $T < \delta$, and that $\frac{\delta}{T} \notin \mathbb{N}$. Then, again for $(0, 0)$ one has

$$b_0 - d_0 - \sum_{n \in \mathbb{Z}, \ |n|T < \delta} \alpha_n 1_{|x-nT| < \delta} = b_0 - d_0 - \sum_{n=-\lfloor \frac{\delta}{T} \rfloor}^{\lfloor \frac{\delta}{T} \rfloor} \alpha_n = 0.$$

Letting $0 < \varepsilon < \delta - \lfloor \frac{\delta}{T} \rfloor \times T$ and using the same analysis as in the previous case one can arrive to the same contradictions. This proves that $T \le \delta$ and that $\frac{\delta}{T} \in \mathbb{N}$.

We denote a solution with period $T = \frac{\delta}{m}$, $m \in \mathbb{N}^*$ along the axis $x = u$ by $\hat{\varphi}_m$. We now prove that solutions $\hat{\varphi}_m$ with $m > 1$ are unstable. We prove this result only for $\hat{\varphi}_2$, the proof is analogous for all $m > 2$. We suppose without loss of generality that the zeros of $\hat{\varphi}_2$ are located at the points $(\frac{n\delta}{2}, \frac{n\delta}{2})$, $n \in \mathbb{Z}$, associated each with a weight $\alpha_n$. The stability condition (v) implies that, for all $k \in \mathbb{Z}$,

$$b_0 - d_0 = \alpha_{k-1} + \alpha_k + \alpha_{k+1}.$$

In particular, this implies that the sequence $(\alpha_k)_{k \in \mathbb{Z}}$ is 3-periodic. Let $\hat{\varphi}_2^*$ be a small perturbation of $\hat{\varphi}_2$ where the only modification with respect to $\hat{\varphi}_2$ is that it is modified on a small ball of radius $\eta > 0$ centered at $(0, 0)$ such that the point $(0, 0)$ is still a local maximum but with a value slightly below zero. Hence the measure $\mu(dx, du)$ associated to $\hat{\varphi}_2^*$ satisfying conditions (i)–(v) has no Dirac mass close to $(0, 0)$. We denote the new weights at the zeros of $\hat{\varphi}_2^*$ located at the points $(n\frac{\delta}{2}, n\frac{\delta}{2})$, $n \in \mathbb{Z}$, $n \ne 1$, by $\alpha_n^*$. As above, these weights satisfy

$$b_0 - d_0 = \alpha_{k-1}^* + \alpha_k^* + \alpha_{k+1}^*, \quad (S11)$$

for all $k \in \mathbb{Z}$ with $k \ne 0$, with the convention that $\alpha_0^* = 0$. Again, the sequence $(\alpha_k^*)_{k \in \mathbb{Z}}$ is 3-periodic and hence $\alpha_{3k}^* = 0$ for all $k \in \mathbb{Z}$ and $\alpha_1^* + \alpha_2^* = b_0 - d_0$. In particular, (S11) is also true for $k = 0$. Since $\alpha_0 > 0$, the fact that $\alpha_1^* + \alpha_2^* = b_0 - d_0$ implies that we either have $\alpha_1^* > \alpha_1$ or $\alpha_2^* > \alpha_2$ (maybe both). Let us assume (without loss of generality) that $\alpha_1^* > \alpha_1$.

We can now consider the solution $\hat{\varphi}_2^*(x, u, t)$ of (7) with initial condition $\hat{\varphi}_2^*$ (with a slight abuse of notation). For all $x \in \mathbb{R}$ such that $|x| > \eta$ and $x \in (k\delta/2, (k+1)\delta/2)$ for some $k \in \mathbb{Z}$,

$$\begin{align*}
\frac{\partial \hat{\varphi}_2^*}{\partial t}(x, x, 0) &= b_0 - d_0 - \alpha_{k-1}^* - \alpha_k^* - \alpha_{k+1}^* - \alpha_{k+2}^* \\
&\quad + D_m \left| \frac{\partial \hat{\varphi}_2^*}{\partial x}(x, x, 0) \right|^2 + \gamma b_0 H \left( \frac{\partial \hat{\varphi}_2^*}{\partial x}(x, x, 0) \right) \\
&= b_0 - d_0 - \alpha_{k-1}^* - \alpha_k^* - \alpha_{k+1}^* - \alpha_{k+2}^* \\
&\quad + D_m \left| \frac{\partial \hat{\varphi}_2}{\partial x}(x, x) \right|^2 + \gamma b_0 H \left( \frac{\partial \hat{\varphi}_2}{\partial x}(x, x) \right) \\
&= \alpha_{k-1} + \alpha_k + \alpha_{k+1} + \alpha_{k+2} - \alpha_{k-1}^* - \alpha_k^* - \alpha_{k+1}^* - \alpha_{k+2}^* \\
&= \alpha_{k-1} - \alpha_{k-1}^*,
\end{align*}$$

where we used in the third equality that $\hat{\varphi}_2$ is a stationary solution of (7). Hence we have proved that $\frac{\partial \hat{\varphi}_2^*}{\partial t}(x, x, 0) < 0$ for all $|x| > \eta$ such that $x \in (k\delta/2, (k+1)\delta/2)$ for some $k \in 3\mathbb{Z} + 2$. In particular, the perturbation does not converge back to $\hat{\varphi}_2$, which is hence unstable.

To conclude we need to check that the solution $\hat{\varphi}_1$ with period $T = \delta$ is stable. This is a non-trivial problem, and we will only check a weak form of stability, assuming only specific perturbations of $\hat{\varphi}_1$. We assume without loss of generality that the zeroes of $\hat{\varphi}_1$ are located at the points $(n\delta, n\delta)$.
Suppose a perturbation is made in such a way that the new function $\hat{\varphi}_1^*$ is different from the original solution only in a small ball of radius $\eta > 0$ centered at $(0,0)$ (this point is arbitrary) such that $\hat{\varphi}_1^*(0,0) = -\epsilon < 0$ is a local maximum. Hence the measure $\mu$ associated to $\hat{\varphi}_1^*$ loses a Dirac mass at $(0,0)$. In this case, $\alpha_k^* = \alpha_k = b_0 - d_0$ for all $k \neq 1$, so that $\frac{\partial \hat{\varphi}_1^*}{\partial t}(x, u, 0) = 0$ for all $(x, u)$ at a distance larger than $\delta$ from $(0,0)$ and

$$ \left. \frac{\partial \hat{\varphi}_1^*(x, u, t)}{\partial t} \right|_{t=x=u=0} = b_0 - d_0 > 0, $$

since there are no neighboring maxima located at a distance smaller than $\delta$. Hence the solution has an initial tendency to approach $\hat{\varphi}_1$. However, the dynamics becomes more complicated after time 0 since the local maxima initially at $(n\delta, n\delta)$ for $n \neq 0$ might move. Of course, this is not sufficient to prove stability, but this strongly suggests that stability should hold true.

# Appendix D. Numerical Complements

## Appendix D.1. Simulation of the individual-based stochastic model

The algorithm used to simulate the individual-based model is based on [42]. Nevertheless, We use a slightly different acceptance/rejection procedure. We proceed as in this reference by constructing recursively sequences $N_k, T_k, X_k, U_k$ respectively of integers, positive real numbers, $N_k$-dimensional vectors in $\bar{\mathcal{X}}$ and $N_k$-dimensional vectors in $\mathcal{U}$, which represent respectively the number of individual in the $k$-th step of the algorithm, the end time of the $k$-th step and the vector of positions and traits of the $N_k$ living individuals at the end of the $k$-th step. We set the potential event rate to $C_{tot} := N_{k-1}(b_0 + d_0 + N_{k-1}/K)$ ($K$ being the carrying capacity) and then decide, based on a parameter $\theta_k$ (a uniformly distributed random variable on $[0, C_{tot}]$), which event can potentially occur at time $T_k = T_{k-1} + E_{k-1}$, where $E_{k-1}$ is an exponential random variable with parameter $C_{tot}$. If $\theta_k < N_{k-1}b_0$, we randomly select an individual $I_k = i$, which will give birth (to a clone or a mutant) with probability $\frac{B(X_{T_k}^i, U_{T_{k-1}}^i)}{b_0}$. If, instead, $N_{k-1}b_0 \le \theta_k \le N_{k-1}b_0 + N_{k-1}d_0$, a randomly selected individual $I_k = i$ dies no matter its position and trait (since the rate of natural death is the same for all individuals). Finally, if $\theta_k \ge N_{k-1}b_0 + N_{k-1}d_0$, two individuals $I_k = i, J_k = j$ are selected randomly. If the distance among them at time $T_k$ is smaller than $\delta$, then individual $i$ dies. Otherwise nothing happens.

## Appendix D.2. Numerical integration of the model's infinite population size limit

To simulate the dynamics of the infinite population-size limit of the model we compute $n(x, u, t)$ according to equation (2) using an explicit scheme where $n(x, u, 0)$ has the form of a Gaussian centered at $(0.5, 0.5)$ with standard deviation $\sigma_0$. For the numerical Turing's stability analysis we modify our boundary conditions so that $\bar{\mathcal{X}} \times \mathcal{U} = T^1 \times T^1$ and use an explicit scheme before and after introducing the perturbation. In both cases we considered a step size $\Delta_x = \Delta_u = 0.025$ and a time-step size $\Delta_t = 0.001$.

## Appendix D.3. Numerical validation of section 3.2

We validate the analytical results from Appendix C.1 with the following numerical scheme. We set periodic boundary conditions and start from the cline-like equilibrium's approximation computed in Appendix C.1. Since this is only an approximation of the cline-like equilibrium, we simulate the dynamics of this system for a time $T_c$ to correct for imperfections. Due to the periodic boundary conditions and the constant value of the distribution along the diagonal lines, this will not cause an infringement of the cline-like condition (a distribution only dependent on $\beta = |x - u|$). We set the
time $T_c$ to be the minimum time that satisfies that the maximum difference between two consecutive time steps at all points is less than a threshold value we set to be equal to the time step size $\Delta_t$. Based on the simulations, we observe that normally $T_c \sim 1$.

The obtained population distribution becomes the approximation of the cline-like equilibrium used to compute numerically the Lyapunov exponents of perturbations. Since the cline-like equilibrium may be unstable this might introduce errors. More precisely, since the approximate cline-like solution is actually a small perturbation of the exact cline-like solution, it has a Lyapunov exponent $\bar{\lambda}$ which might be positive. In particular, the Lyapunov exponent of a perturbation of the approximate equilibrium is not close to $\lambda_{m_1,m_2}$, but should rather be close to $\max\{\lambda_{m_1,m_2}, \bar{\lambda}\}$, hence introducing a threshold effect in the numerical computations. This threshold effect is indeed observed for small values of $\lambda_{m_1,m_2}$ in Fig. 10(B). Note also that we cannot use the clever method of computation of the cline-like equilibrium as done in [43] since, because of our boundary conditions, we are unable to reduce the spatial window of numerical resolution of the PDE so that possible perturbations have very strong frequencies and do not destabilize the equilibrium. Nevertheless, the approximation was good enough so that the dynamics always converged towards a distribution that only depends on values of $|x - u|$.

After we retrieve an approximation of the cline-like equilibrium distribution $n_0(x, u) = n_0(|x-u|)$, we introduce a perturbation of the form:

$$ \hat{n}_0(x, u) = (1 + \varepsilon \cos(2(m_1 u + m_2 x)))n_0(x, u), \qquad (D.1) $$

$\varepsilon \ll 1$, for several integer values of $m_1$ and $m_2$. Next, we simulate the dynamics of $\hat{n}(x, u, t)$, $\hat{n}(x, u, 0) = \hat{n}_0(x, u)$ and measure the growth or decay of the perturbation by calculating the exponent:

$$ \lambda_{num} = \max_{i \in I} \left\{ \left\langle \frac{\log(|n_0(x_i, u_i) - \hat{n}(x_i, u_i, n\delta t)| / |n_0(x_i, u_i) - \hat{n}_0(x_i, u_i)|)}{n\delta t} \right\rangle_{n\delta t \le T_p} \right\}. $$

where $\delta t$ is the time step size and $I$ is an enumeration of the points along the line $x = u$ of the grid on which we simulate the dynamics of the system. Again —although $\lambda_{num}$ is not necessarily equal to $\lambda_{m_1,m_2}$— if $\lambda_{num}$ is positive, the cline-like equilibrium is unstable for perturbations of the form in equation (D.1). If this is the case for at least one pair $(m_1, m_2)$ it is expected that the attracting state of the system will be distributed among isolated clusters.

In Fig. 10 we show the computation of $\lambda_{num}$ when varying $D_m$ and $\delta$ for perturbations with $(m_1, m_2) = (0, 7)$ and $(m_1, m_2) = (2, 5)$, this last pair resulted in the highest registered value of $\lambda_{num}$ in the considered ranges of both parameters. As observed, results in figure Fig. 10(A) show an almost perfect consistency with their analytical counterpart when approximating the transition value $D_m^*$ of $D_m$ for both pairs $(m_1, m_2)$ considered. Thus we corroborate that when increasing the dispersal rate $D_m$, the transition from the clustering attractor to the cline-like one occurs at a threshold value $D_m^* \approx 1.6 \times 10^{-4}$. Again this value would have been miscalculated if only perturbations along the spatial component had been considered.

The numerical simulations for determining $\lambda_{num}$ as function of $\delta$ are harder due to practical issues. More precisely, the spatial discretization of space as cells of width $\Delta_x$ needs to be chosen so that $\delta/\Delta_x$ is an integer. Given $\Delta_x$, for values of $\delta$ such that $\delta/\Delta_x$ is an integer, the numerical scheme is exactly the same as for $\delta' = \Delta_x[\delta/\Delta_x]$ and so the computed values of $\lambda$ are the same. Hence, to have a precise estimate of $\lambda$ as a function of $\delta$, we need to take very small $\Delta_x$. However, this is impractical in terms of numerical cost since the stability condition of explicit finite difference scheme requires to
take time-steps $\Delta_t$ even smaller. Therefore, our numerical validation for the dependence with respect to $\delta$ has been validated for $\Delta_x = 0.025$, and hence for the values $\delta = 0.025, 0.05, 0.075, 0.1, 0.125$, shown in Fig. 10(B). Of course, these five values are not sufficient to determine precisely the location of the transition from cline-like to clustering patterns, predicted to be close to $\delta^* \approx 0.076$ from the analytical approximation (5). Nevertheless, results are consistent, indicating that the transition occurs at a value larger but close to $\delta = 0.075$ and that $\lambda_{num}$ is smaller when $(m_1, m_2) = (0, 7)$ than when $(m_1, m_2) = (2, 5)$.

## Appendix D.4. Numerical integration of the Hamilton-Jacobi equation

We take as initial condition $\varphi(x, u, 0) = -(x - 0.5)^2 - (u - 0.5)^2$ and we compute the dynamics of $\varphi$ according to the Hamilton-Jacobi equation (7) using an explicit scheme. In order to always satisfy condition (i) we need to use a scheme with variable time-step where the time step is set to $\Delta t$ unless the condition of negativity is violated for at least one point. In this case, we take the minimum time step for which the nonpositivity condition is maintained and used it instead. We then obtain a function $\varphi(x, u, t + \Delta t')$ admitting at least one more zero that $\varphi(x, u, t)$. In the next time step, some of the zeroes might disappear if $\partial_t \varphi(x, u, t + \Delta t') < 0$ at some of these points. To compute the Hamiltonian function $H(\frac{\partial \varphi(x, u, t)}{\partial u})$, we use the expression in Eq. (8).

To compute the values of the $\alpha_i$'s in the density function $\mu(dy, dw, t)$ (introduced in equation (9)) which satisfy conditions (iv), (v), we use the following scheme. If there is just one pair $(x, u)$ satisfying $\varphi(x, u, t) = 0$, i.e $|\Omega_t| = 1$, we set its respective coefficient $\alpha$ to the unique value satisfying condition (v).

Consider now the case when $|\Omega_t| \ge 2$. Let $(x_i, u_i) \in \Omega_t$, $1 \le i \le m := |\Omega_t|$ denote the elements of $\Omega_t$. In order to find the appropriate values of the $\alpha_i$'s we first divide the set $\Omega_t$ into smaller sets $A_k$ constructed inductively by using the following method. Initially we construct $A_1$ as the set containing the first element $(x_1, u_1)$ of $\Omega_t$. Then, for every $1 < i \le m$, the pair $(x_i, u_i)$ is added to a set $A_k$ if there exists a pair $(x_j, u_j) \in A_k$ such that $|x_i - x_j| < \delta$. If there exists more than one $A_k$ satisfying this condition, then $A_k$ will be redefined as the union of such sets. One can picture the sets $A_k$ as sets containing the pairs $(x_i, u_i)$ with correlated values of $\alpha_i$'s in the sense that elements of different sets do not compete among them. Likewise, in the end, no elements from two different sets will be at a distance less than $\delta$. Once this is done, in each subset $A_k$ we verify whether conditions (iv), (v) can be satisfied when setting one of the values of the $\alpha_i$'s to the value imposed by condition (v) and setting the values of the other $\alpha_j$'s to zero. If this is the case for at least one element $(x_i, u_i)$ in each $A_k$, we take such values of the $\alpha$'s to be coefficients of $\mu(dy, dw, t)$. Otherwise, if this conditions are not satisfied we compute numerically the values of the $\alpha_i$'s as the solutions of the system appearing from imposing condition (v) to every element in $|\Omega_t|$;

$$b_0 - b_1(x_1 - u_1)^2 - d_0 - \sum_{i \le m} 1_{|x_1 - x_i| < \delta \alpha_i} = 0,$$

$$b_0 - b_1(x_2 - u_2)^2 - d_0 - \sum_{i \le m} 1_{|x_2 - x_i| < \delta \alpha_i} = 0,$$

$\vdots$

$$b_0 - b_1(x_m - u_m)^2 - d_0 - \sum_{i \le m} 1_{|x_m - x_i| < \delta \alpha_i} = 0.$$

In principle we should look for other solutions for which the subset $S_k$ of elements $(x_i, u_i)_k$ associated to non-zero $\alpha$'s in each set $A_k$ has size 2, then 3 and so on. In practice, however, we only
look at cases where $|S_k| = 1$ or $|\cup_k S_k| = m$ and these ones were enough at short time scales like the ones we use in our simulations.

## References Supplementary Material

[61] Everitt WN. A catalogue of Sturm-Liouville differential equations. In: Sturm-Liouville Theory. Springer; 2005. p. 271–331.

[62] Levitan B, Sargsjan I. Sturm-Liouville and Dirac operators, volume 59 of Mathematics and its Applications (Soviet Series). Kluwer Academic Publishers Group, Dordrecht; 1991.

[63] Teschl G. Ordinary differential equations and dynamical systems. vol. 140. American Mathematical Society Providence, RI; 2012.

[64] Mirrahimi S, Roquejoffre JM. A class of Hamilton-Jacobi equations with constraint: uniqueness and constructive approach. J Differential Equations. 2016;260(5):4717–4738. Available from: http://dx.doi.org/10.1016/j.jde.2015.11.027.

[65] Mirrahimi S, Roquejoffre JM. Uniqueness in a class of Hamilton-Jacobi equations with constraints. C R Math Acad Sci Paris. 2015;353(6):489–494. Available from: http://dx.doi.org/10.1016/j.crma.2015.03.005.

[66] Champagnat N, Méléard S. Invasion and adaptive evolution for individual-based spatially structured populations. Journal of Mathematical Biology. 2007;55(2):147–188.

[67] Leimar O, Doebeli M, Dieckmann U. Evolution of phenotypic clusters through competition and local adaptation along an environmental gradient. Evolution. 2008;62(4):807–822.