Asymptotic analysis of selection-mutation models in the presence of multiple fitness peaks

TL;DR

This paper analyzes the long-term behavior of phenotype-structured models with multiple fitness peaks, revealing how populations concentrate around peaks and how phenotypic variations influence this process.

Contribution

It provides a rigorous analysis of the asymptotic limits of models with and without phenotypic variation, including explicit formulas and criteria for population distribution convergence.

Findings

Long-time solutions concentrate at fitness peaks as Dirac masses.

Explicit formulas for weights of population at each peak.

Conditions under which phenotypic variation affects the number of peaks.

Abstract

We study the long-time behaviour of phenotype-structured models describing the evolutionary dynamics of asexual species whose phenotypic fitness landscape is characterised by multiple peaks. First we consider the case where phenotypic variations do not occur, and then we include the effect of heritable phenotypic changes. In the former case the model is formulated in terms of an integrodifferential equation for the phenotype distribution of the individuals of the species, whereas in the latter case the evolution of the phenotype distribution is governed by a non-local parabolic equation whereby a linear diffusion operator captures the presence of phenotypic variations. We prove that the long-time limit of the solution to the integrodifferential equation is unique and given by a measure consisting of a weighted sum of Dirac masses centred at the peaks of the phenotypic fitness landscape. We also derive an explicit formula to compute the weights in front of the Dirac masses. Moreover, we demonstrate that the long-time solution of the non-local parabolic equation exhibits a qualitatively similar behaviour in the asymptotic regime where the diffusion coefficient modelling the rate of phenotypic variations tends to zero. However, we show that the limit measure of the non-local parabolic equation may consist of less Dirac masses, and we provide a sufficient criterion to identify the positions of their centres. Finally, we provide a detailed characterisation of the speed of convergence of the integral of the solution (i.e. the population size) to its long-time limit for both models. Taken together, our results support a more in-depth theoretical understanding of the conditions leading to the emergence of stable phenotypic polymorphism in asexual species.

Equations

$$\frac{\partial n}{\partial t}=R(x,\rho (t))n,n\equiv n(t,x),(t,x)\in (0,\infty )\times \overline{\Omega }.$$

$$\rho (t):={\int }_{\Omega }n(t,x)\mathrm{d}x,$$

$$\frac{\partial n}{\partial t}=R(x,\rho (t))n+\beta \Delta n,n\equiv n(t,x),(t,x)\in (0,\infty )\times \Omega ,$$

$$R(x,\rho (t)):=r(x)-\rho (t).$$

$$\{\begin{array}{l}\frac{\partial n}{\partial t}=\left(r(x)-\rho (t)\right)n,n\equiv n(t,x),(t,x)\in (0,\infty )\times \overline{\Omega },\\ \\ n(0,x)=n^0(x),n^0\in C(\overline{\Omega }),n^0⩾0,n^0\not\equiv 0\end{array}$$

$$\{\begin{array}{l}\frac{\partial n_{\beta }}{\partial t}-\beta \Delta n_{\beta }=\left(r(x)-{\rho }_{\beta }(t)\right)n_{\beta },n_{\beta }\equiv n_{\beta }(t,x),(t,x)\in (0,\infty )\times \Omega ,\\ \\ \nabla n_{\beta }(t,x)\cdot \nu (x)=0,(t,x)\in (0,\infty )\times \partial \Omega ,\\ \\ n_{\beta }(0,x)=n^0(x),n^0\in C(\Omega ),n^0⩾0,n^0\not\equiv 0,\end{array}$$

$$r:\overline{\Omega }\to \mathbb{R},r\in C(\overline{\Omega }),\underset{x\in \Omega }{\max }r(x)=r_M>0$$

$$\arg \max (r)=\{{\bar{x}}_1,\dots ,{\bar{x}}_N\}\subset \Omega \text{with}N⩾2.$$

$$\mathrm{supp}(n^0)\cap \arg \max (r)\ne \varnothing$$

$$r\in C^2(\overline{\Omega }),\det (H_i)<0\text{for}i=1,\dots ,N,$$

$${\int }_{B({\bar{x}}_i,\epsilon )}f(x)e^{r(x)t}\mathrm{d}x\sim (2\pi {)}^{d/2}\frac{f({\bar{x}}_i)}{\sqrt{|\det (H_i)|}}\frac{e^{r_{M}t}}{t^{\frac{d}{2}}}\text{as }t\to \infty ,$$

$${\int }_{B({\bar{x}}_i,\epsilon )}f(x)e^{r(x)t}\mathrm{d}x=o\left(\frac{e^{r_{M}t}}{t^{\frac{d}{2}}}\right)\text{as }t\to \infty .$$

$${\int }_{B({\bar{x}}_i,\epsilon )}f(x)e^{r(x)t}\mathrm{d}x=\left(A_i+\frac{B_i}{t}+o\left(\frac{1}{t}\right)\right)\frac{e^{r_{M}t}}{t^{\frac{d}{2}}}\text{as }t\to \infty ,$$

$$L_{\beta }=\beta \Delta +r$$

$${\lambda }_{\beta }=\underset{\varphi \in H^1(\Omega )\setminus \{0\}}{\inf }\mathcal{R}(L_{\beta },\varphi ),$$

$$\mathcal{R}(L_{\beta },\varphi ):=\frac{\beta {\int }_{\Omega }{\left|\nabla \varphi (x)\right|}^2\mathrm{d}x-{\int }_{\Omega }r(x){\varphi }^2(x)\mathrm{d}x}{{\int }_{\Omega }{\varphi }^2(x)\mathrm{d}x}.$$

$${\gamma }_{\beta }:={\lambda }_{2,\beta }-{\lambda }_{\beta }<0.$$

$$S:=\{x_1=0\},$$

$$x=(x_1,x_2,\dots ,x_d)\in \Omega ⟹(-x_1,x_2,\dots ,x_d)\in \Omega .$$

$$\hat{f}(x)=f(-x_1,x_2,\dots ,x_d)\text{for }x\in \Omega .$$

$$\hat{r}(x)=r(x)\text{for all }x\in \Omega ,$$

$$\beta \Delta {\psi }_{\beta }+r{\psi }_{\beta }=-{\lambda }_{\beta }{\psi }_{\beta }\text{in }\Omega ,$$

$$\beta \Delta {\hat{\psi }}_{\beta }+\hat{r}{\hat{\psi }}_{\beta }=-{\lambda }_{\beta }{\hat{\psi }}_{\beta }\text{in }\Omega .$$

$$\beta \Delta {\hat{\psi }}_{\beta }+r{\hat{\psi }}_{\beta }=-{\lambda }_{\beta }{\hat{\psi }}_{\beta }\text{in }\Omega ,$$

$${\int }_{\Omega }\phi (x)\mu (x)\mathrm{d}x={\int }_{\Omega }\phi \mathrm{d}\mu .$$

$$\int_{\Omega}\varphi(x)\,\mu_{k}(x)\,{\rm d}x\;\xrightarrow[k\rightarrow\infty]{\;}\int_{\Omega}\varphi(x)\,\mu(x)\,{\rm d}x\quad\forall\,\varphi\in C\big{(}\overline{\Omega}\big{)}.$$

$$\int_{\Omega}\varphi(x)\,\mu_{k}(x)\,{\rm d}x\xrightarrow[k\rightarrow\infty]{\,}0\quad\forall\,\varphi\in C\big{(}\overline{\Omega}\big{)}\;\text{ s.t. }\operatorname{supp}(\varphi)\cap\omega=\emptyset,$$

$$\int_{\Omega}\varphi(x)\,\mu(x)\,{\rm d}x\;=\;\int_{\Omega}\hat{\varphi}(x)\,\mu(x)\,{\rm d}x\quad\forall\,\varphi\in C\big{(}\overline{\Omega}\big{)}\quad\text{s.t.}\quad\hat{\varphi}(x)=\varphi(x)\;\;\forall\,x\in\Omega$$

$$n(t,x)=n^0(x)e^{r(x)t-{\int }_0^t\rho (s)\mathrm{d}s}$$

$$\rho (t)\underset{t\to \infty }{\to }{r}_M$$

Full text

1 Introduction

Phenotype-structured models formulated as integrodifferential equations (IDEs) or non-local partial differential equations (PDEs) have been increasingly used as a theoretical framework to study evolutionary dynamics in a variety of asexual populations [2, 4, 7, 8, 9, 10, 11, 14, 13, 12, 16, 18, 20, 24, 30, 32, 39, 35, 36, 37, 40, 41, 42, 44, 45, 46, 47, 50]. In these models, the phenotypic state of each individual is represented by a continuous real variable $x$, and the phenotypic distribution of the individuals within the population at a given time $t\geqslant 0$ is described by a function $n(t,x)~{}\geqslant~{}0$. In many scenarios of biological and ecological interest one can assume $x\in\Omega$, where $\Omega$ is a smooth bounded domain of $\mathbb{R}^{d}$, $d\geqslant 1$.

We focus here on the case where, in the absence of phenotypic changes, the evolution of the population density function $n(t,x)$ is governed by an IDE of the form

[TABLE]

The function $R(x,\rho(t))$ represents the net per capita growth rate of individuals in the phenotypic state $x$, under the environmental conditions determined by the population size

[TABLE]

and models the effect of natural selection. In fact, depending on the sign of $R(x,\rho(t))$ the number density of individuals in the phenotypic state $x$ will either grow or decay at time $t$. The function $R(x,\rho(t))$ can thus be seen as the phenotypic fitness landscape of the population [26, 49].

In the ecological and biological scenarios whereby the effect of heritable, spontaneous phenotypic changes need to be taken into account, a linear diffusion operator can be included in the IDE (1.1). This leads to a non-local parabolic PDE of the form

[TABLE]

where the diffusion coefficient $\beta>0$ models the rate of phenotypic change. Since phenotypic changes preserve the total number of individuals within the population, zero Neumann is the most natural choice of boundary conditions for the non-local PDE (1.3).

The way in which the fitness function $R(x,\rho(t))$ is defined depends largely on the underlying application problem, and we consider here the prototypical definition

[TABLE]

In definition (1.4), the function $r(x)$ is the net per capita growth rate of the individuals in the phenotypic state $x$ (i.e. the difference between the rate of proliferation through asexual reproduction and the rate of death under natural selection). Hence the maximum points of this function correspond to the fitness peaks (i.e. the peaks of the phenotypic fitness landscape of the population). Moreover, the saturating term $-\rho(t)$ models the limitations on population growth imposed by carrying capacity constraints (e.g. limited availability of space and resources).

In the framework of the IDE (1.1), or the non-local PDE (1.3), a mathematical depiction of phenotypic adaptation can be obtained by studying the long-time behaviour of the population density function $n(t,x)$. In this regard, whilst the case of one single fitness peak has been broadly studied [1, 5, 15, 17, 22, 43, 48], there is a paucity of literature concerning the case where multiple fitness peaks are present, with the exception of the asymptotic results presented in [3, 11, 19, 21, 23, 31, 39, 38, 47, 51]. Based on these few previous results, we expect the solution to the IDE (1.1) complemented with (1.4) to converge to a limit measure given by a sum of weighted Dirac masses centred at the maximum points of the function $r(x)$ as $t\to\infty$, and we envisage the long-time solution of the non-local PDE (1.3) to exhibit a qualitatively similar behaviour in the asymptotic regime where $\beta\to 0$. This represents a mathematical formalisation of the biological notion that the phenotypic variants corresponding to the fitness peaks (i.e. the fittest phenotypic variants) are ultimately selected. However, the following key questions still remain open: is there a unique weighted sum of Dirac masses defining the limit of the solution to the IDE (1.1) when $t\to\infty$? If so, what are the only admissible values of the weights in front of the different Dirac masses? Are there any differences between the long-time solution to the IDE (1.1) and the asymptotic limit when $\beta\to 0$ of the long-time solution to the non-local PDE (1.3)? If so, how does the presence of the diffusion term change the values of the weights associated with the Dirac masses? How does the speed of convergence of the population size $\rho(t)$ to its long-time limit differs between the IDE (1.1) and the non-local PDE (1.3)?

In this paper, we address these questions focussing on the case where the fitness function is defined via (1.4). In summary, using Laplace’s method we derive an explicit formula to compute the weights in front of the Dirac masses that constitute the long-time limit of the solution to the IDE (1.1), and show that the weights are uniquely determined by the initial condition $n(0,x)$ and the Hessian $H(r)$ at the maximum points of the function $r(x)$ (vid. Theorem 2). Moreover, exploiting the properties of the principal eigenpair of the elliptic differential operator $\beta\Delta+r$, we prove that when $\beta\to 0$ the long-time limit of the non-local PDE (1.3) converges to a limit measure consisting of a weighted sum of Dirac masses with non-negative weights centred at the maximum points of the function $r(x)$ (vid. Proposition 1). We also derive sufficient conditions for the limit measure to be unique and demonstrate that, ceteris paribus, there can be a significant difference between this limit measure and the limit measure of the IDE (1.1). In particular, we show that, unless ad hoc symmetry assumptions are made (vid. Proposition 2), the limit measure for the non-local PDE (1.3) may consist of less Dirac masses than the limit measure for the IDE (1.1) (i.e. a smaller number of weights will be strictly positive), and we provide a sufficient criterion to identify the maximum points of the function $r(x)$ corresponding to the Dirac masses with positive weights (vid. Proposition 3). This criterion relies on a suitable multidimensional characterisation of the concavity of the function $r(x)$ at the maximum points which is borrowed from semiclassical analysis. Finally, we carry out a detailed characterisation of the speed of convergence of the population size $\rho(t)$ to its long-time limit both for the IDE (1.1) and for the non-local PDE (1.3). Taken together, our results support a more in-depth theoretical understanding of the conditions leading to the emergence of stable polymorphism in asexual populations.

The remainder of the paper is organised as follows. In Section 2 we introduce our main assumptions and a few technical preliminaries. In Section 3 we carry out a qualitative and quantitative characterisation of the solution to the IDE (1.1) complemented with (1.4) when $t\to\infty$, while in Section 4 we study the asymptotic properties of the solution to the non-local PDE (1.3) complemented with (1.4) by letting first $t\to\infty$ and then $\beta\to 0$. In both sections, we present a sample of numerical solutions that confirm the analytical results obtained. Section 5 concludes the paper by providing a brief overview of possible research perspectives.

2 Main assumptions, notation and preliminaries

In this paper, we will consider the following Cauchy problem for the IDE (1.1)

[TABLE]

and the following initial-boundary value problem for the non-local parabolic PDE (1.3)

[TABLE]

where $\nu(x)$ is the outward normal to the boundary $\partial\Omega$ at the point $x\in\partial\Omega$.

Main assumptions on the function $r(x)$.

In order to prevent $n(t,\cdot)$ from vanishing as $t\to\infty$, we will assume

[TABLE]

and, being interested in the case where $r(x)$ has multiple maximum points, we will also assume

[TABLE]

Notice that we assume all points $\bar{x}_{i}$ to belong to the interior of $\Omega$ in order to simplify the presentation. However, part of our results can be extended to the case where the set $\arg\max(r)$ contains some boundary points and, when appropriate, we will comment on how the proofs presented here could be adapted to such a case. Where necessary, we will make the additional assumptions

[TABLE]

and

[TABLE]

where $H_{i}$ is the Hessian $H(r)$ evaluated at the point $\bar{x}_{i}$. Assumptions (2.10) ensure that each maximum point $\bar{x}_{i}$ is nondegenerate.

Laplace’s method.

We recall some useful results on the asymptotic expansion of integrals involving exponentials, usually referred to as Laplace’s method. We refer the reader to [54] for a general presentation of such an asymptotic method and for the related proofs. Let $r$ satisfy assumptions (2.7), (2.8) and (2.10), and assume $f\in C(\Omega)$. For $\bar{x}_{i}\in\arg\max(r)$ and $\varepsilon$ small enough so that $\bar{x}_{i}$ is the only maximum point of $r(x)$ in the ball $B(\bar{x}_{i},\varepsilon)$, Laplace’s method ensures that if $\bar{x}_{i}\in\operatorname{supp}\left(f\right)$ then

[TABLE]

whereas if $\bar{x}_{i}\notin\operatorname{supp}\left(f\right)$ then

[TABLE]

Moreover, if $f(x)$ and $r(x)$ are, respectively, of class $C^{1}$ and $C^{3}$ on the ball $B(\bar{x}_{i},\varepsilon)$, higher order terms of the asymptotic expansion can be computed. In this case, Laplace’s method ensures that

[TABLE]

where $A_{i}$ and $B_{i}$ are real constants, the values of which are not relevant for our purposes.

Preliminaries about the operator $\beta\Delta+r$.

We will consider the elliptic differential operator

[TABLE]

acting on functions defined on $\Omega$, and we will denote by $(\lambda_{\beta},\psi_{\beta})$ the principal eigenpair of $L_{\beta}$ (i.e. $L_{\beta}\psi_{\beta}=-\lambda_{\beta}\psi_{\beta}$) with zero Neumann boundary condition. A useful result is established by the following lemma, which follows from the Krein-Rutman theorem [25, 33].

Lemma 1**.**

The principal eigenvalue $\lambda_{\beta}$ is simple and there is a unique normalised positive eigenfunction $\psi_{\beta}$ associated with $\lambda_{\beta}$. The eigenfunction $\psi_{\beta}$ is smooth and the principal eigenvalue $\lambda_{\beta}$ is given by

[TABLE]

where $\mathcal{R}$ denotes the Rayleigh quotient, i.e.

[TABLE]

The infimum is attained only when $\phi$ is a multiple of $\psi_{\beta}$.

We will also make use of the second eigenvalue $\lambda_{2,\beta}<\lambda_{\beta}$ of the operator $L_{\beta}$. For a function $u\in L^{2}(\Omega)$, we will denote by $\alpha_{2}(u)$ the $L^{2}(\Omega)$-projection of $u$ onto the finite-dimensional eigenspace associated to $\lambda_{2,\beta}$, and we will denote the opposite to the spectral gap of the operator $L_{\beta}$ by

[TABLE]

If the set $\Omega$ is symmetric with respect to some hyperplane $S$, which without loss of generality we will define as

[TABLE]

then

[TABLE]

Under assumption (2.18) and given a function $f:\overline{\Omega}\to\mathbb{R}$, we will use the notation

[TABLE]

With this notation, we will say $f$ to be symmetric with respect to the hyperplane $S$ if $\hat{f}(x)=f(x)$ for all $x\in\Omega$. A useful result is established by the following lemma.

Lemma 2**.**

If the set $\Omega$ satisfies assumption (2.18) and the function $r(x)$ is symmetric with respect to the hyperplane $S$, i.e. if

[TABLE]

then the principal eigenfunction $\psi_{\beta}$ is symmetric with respect to the hyperplane $S$.

Proof.

Since $L_{\beta}\psi_{\beta}=-\lambda_{\beta}\,\psi_{\beta}$, that is,

[TABLE]

if assumption (2.18) is satisfied then

[TABLE]

Under the additional assumption (2.19) the latter elliptic equation implies that

[TABLE]

that is, $L_{\beta}\hat{\psi}_{\beta}=-\lambda_{\beta}\,\hat{\psi}_{\beta}$. Hence $\hat{\psi}_{\beta}$ is a normalised positive eigenfunction associated with the principal eigenvalue $\lambda_{\beta}$, and the uniqueness of the principal eigenfunction of the operator $L_{\beta}$ ensures that $\hat{\psi}_{\beta}=\psi_{\beta}$. ∎

Function space framework

Throughout the paper, we will consider the space of Radon measures \mathcal{M}^{1}\big{(}\overline{\Omega}\big{)} as the dual of the space of continuous functions C\big{(}\overline{\Omega}\big{)}. With a slight abuse of notation, the integral of a function \varphi\in C\big{(}\overline{\Omega}\big{)} against a measure \mu\in\mathcal{M}^{1}\big{(}\overline{\Omega}\big{)} will be denoted in the same way as the integral of the product of the function $\varphi$ with an $L^{1}$-function, i.e.

[TABLE]

Sequences of functions $f_{k}$ in $L^{1}(\Omega)$ will be regarded as elements of the bigger space \mathcal{M}^{1}\big{(}\overline{\Omega}\big{)}. Given a sequence $\mu_{k}$ in \mathcal{M}^{1}\big{(}\overline{\Omega}\big{)}, we will write $\displaystyle{\mu_{k}\xrightharpoonup[k\rightarrow\infty]{}\mu}$ to indicate the weak-$*$ convergence of $\mu_{k}$ to $\mu$, namely that

[TABLE]

Moreover, we will say that the sequence $\mu_{k}$ concentrates on a set $\omega\subset\Omega$ if

[TABLE]

and we will use the result given by the following lemma, which is a well-known fact in measure theory.

Lemma 3**.**

If a sequence $\mu_{k}$ in \mathcal{M}^{1}\big{(}\overline{\Omega}\big{)} concentrates on a finite set $\{\bar{x}_{1},\ldots,\bar{x}_{N}\}$ and $\displaystyle{\mu_{k}\xrightharpoonup[k\rightarrow\infty]{}\mu}$ then the limit measure $\mu$ must be a linear combination of Dirac masses centred at the points $\bar{x}_{1},\ldots,\bar{x}_{N}$.

Finally, we will say that a measure \mu\in\mathcal{M}^{1}\big{(}\overline{\Omega}\big{)} is symmetric with respect to the hyperplane $S$ if

[TABLE]

and we will use the result given by the following lemma, the proof of which is straightforward.

Lemma 4**.**

Let $\mu_{k}$ be a sequence of symmetric measures in \mathcal{M}^{1}\big{(}\overline{\Omega}\big{)}. If $\displaystyle{\mu_{k}\xrightharpoonup[k\rightarrow\infty]{}\mu}$ then the limit measure $\mu$ is symmetric as well.

3 Long-time behaviour of the Cauchy problem (2.5)

In this section, we study the asymptotic behaviour of the solutions to the Cauchy problem (2.5) when $t\to\infty$ (Section 3.1), and we provide a sample of numerical solutions that confirm the analytical results obtained (Section 3.2).

3.1 Asymptotic analysis

Under assumptions (2.7) there exists a unique non-negative solution $n\in C([0,+\infty);L^{1}(\Omega))$ of the Cauchy problem (2.5) [21, 47]. Moreover, solving the Cauchy problem (2.5) yields the semi-explicit formula

[TABLE]

and, if assumption (2.8) is satisfied as well, the solution $n(t,x)$ is known to concentrate on the set $\arg\max(r)$ when $t\to\infty$, as established by the following theorem.

Theorem 1**.**

Under assumptions (2.7)-(2.9), the solution to the Cauchy problem (2.5) is such that

[TABLE]

and, up to extraction of subsequences,

[TABLE]

In the case where

[TABLE]

Theorem 1 can be proved through a few simple calculations building upon the method presented in [48], as shown by the proof provided in the Appendix for the sake of completeness. Alternatively, when assumption (3.23) is satisfied, one can prove Theorem 1 using the method of proof presented in [47], which relies on the observation that the semi-explicit solution (3.20) can be made explicit as

[TABLE]

On the other hand, in the case where assumption (3.23) is not satisfied, the proof is much more intricate and requires the use of a Lyapunov functional of the form

[TABLE]

with $n^{\infty}$ being any measure concentrated on the set $\arg\max(r)$ and having total mass $r_{M}$. We refer the interested reader to [31, 51] for a proof of Theorem 1 in such a more general case.

We prove here that the coefficients $a_{1},\ldots,a_{N}$ that define the limit measure (3.22) are uniquely determined by the initial condition and by the Hessian of $r(x)$ at the points $\bar{x}_{1},\ldots,\bar{x}_{N}$, which entails the convergence of the whole trajectory $n(t,\cdot)$ to a unique limit point as $t\to\infty$. These results are summarised by the following theorem, which also provides a characterisation of the rate of convergence of the total mass $\rho(t)$ to the long-term limit $r_{M}$.

Theorem 2**.**

If assumptions (2.7)-(2.10) are satisfied, then the solution to the Cauchy problem (2.5) is such that

[TABLE]

with

[TABLE]

where $A>0$ is a normalising constant such that $\displaystyle{\sum_{i=1}^{N}a_{i}=1}$.

Moreover, if the functions $n^{0}(x)$ and $r(x)$ are, respectively, of class $C^{1}$ and $C^{3}$ in a neighbourhood of each maximum point $\bar{x}_{1},\ldots,\bar{x}_{N}$ then

[TABLE]

Proof.

Throughout the proof we will use the following notation

[TABLE]

The results established by Theorem 1 ensure that we can extract a subsequence $n(t_{k},\cdot)$ such that

[TABLE]

Moreover, considering $\varepsilon>0$ small enough so that $\bar{x}_{i}$ is the only maximum point of the function $r(x)$ in the ball $B(\bar{x}_{i},\varepsilon)$ for every $i=1,\ldots,N$ and integrating over $B(\bar{x}_{i},\varepsilon)$ both sides of the expression (3.20) for $n(t,x)$ with $t=t_{k}$ we find that

[TABLE]

The long-time behaviour of the integral on the right-hand side of the above equation can be characterised using Laplace’s method. In particular, the asymptotic relation (2.11) ensures that if $\bar{x}_{i}\in\operatorname{supp}(n^{0})$ then

[TABLE]

where $J$ is defined according to (3.27). On the other hand, the asymptotic relation (2.12) ensures that if $\bar{x}_{i}\notin\operatorname{supp}(n^{0})$ then

[TABLE]

Taken together, the integral identity (3.28) and the asymptotic relations (3.29) and (3.30) allow us to conclude that there exist some real constants $K>0$ and $A>0$ such that

[TABLE]

and

[TABLE]

We remark that $K$ cannot be [math] because otherwise $\rho(t)$ would converge to [math], which cannot be since $\rho(t)\to r_{M}$ [cf. the asymptotic result (3.21)]. Hence, the coefficients $a_{i}$ that define the limit measure are uniquely determined and the limit measure is unique. This ensures that the whole trajectory $n(t,\cdot)$ converges to the limit point given by (3.24) and (3.25) as $t\to\infty$.

To prove claim (3.26) we proceed as follows. We note that the asymptotic result (3.31) now holds true for $J(t)$ and not for a mere subsequence $J(t_{k})$. This implies that

[TABLE]

To conclude we only need to show that $\rho(t)$ has an asymptotic expansion of the form $a+\frac{b}{t}+o\left(\frac{1}{t}\right)$ as $t\to\infty$. The coefficients $a$ and $b$ are then necessarily $r_{M}$ and $-\frac{d}{2}$, owing to (3.32).

From now on, $C_{1}$ and $C_{2}$ will denote some generic real constants which might vary from line to line.

We choose again $\varepsilon>0$ small enough so that $\bar{x}_{i}$ is the only maximum point of the function $r(x)$ in the ball $B(\bar{x}_{i},\varepsilon)$ for every $i=1,\ldots,N$. Integrating over $B(\bar{x}_{i},\varepsilon)$ both sides of the expression (3.20) for $n(t,x)$ we find that

[TABLE]

The second term on the right-hand side of the above equation decays exponentially to [math] as $t\to\infty$, since $\rho(t)\to r_{M}$. Moreover, choosing $\varepsilon$ small enough so that – under the additional assumption that the functions $n^{0}(x)$ and $r(x)$ are, respectively, of class $C^{1}$ and $C^{3}$ in a neighbourhood of each maximum point $\bar{x}_{1},\ldots,\bar{x}_{N}$ – we have $n^{0}\in C^{1}\left(B(\bar{x}_{i},\varepsilon)\right)$ and $r\in C^{3}\left(B(\bar{x}_{i},\varepsilon)\right)$ for every $i=1,\ldots,N$, we can use the asymptotic expansion (2.13) and in so doing obtain the following asymptotic expression for the first integral on the right-hand side of the latter equation

[TABLE]

Taken together, these results yield

[TABLE]

with $J$ defined according to (3.27).

We now prove that

[TABLE]

from which we can infer that $\rho(t)$ satisfies an asymptotic expansion of the same form, thus concluding the proof. In order to prove (3.34), we notice that a sufficient condition for this to hold is that the function $u(t):=e^{\int_{0}^{t}\rho(s)\,ds}$ satisfies the estimate

[TABLE]

In order to prove (3.35), we differentiate $u$ to obtain

[TABLE]

where the last equality has been established above. Since

[TABLE]

for any $\alpha>0$, we can conclude that estimate (3.36) still holds true after integration and, therefore, estimate (3.35) is satisfied. ∎

Remark 1**.**

In the case where $d=1$, expression (3.25) reads as

[TABLE]

Remark 2**.**

The results of Theorem 2 can be extended to the case where some maximum points of the function $r(x)$ belong to the boundary $\partial\Omega$, a case that might be relevant for applications. In particular, letting $\partial\Omega$ be sufficiently smooth and using Laplace’s method one can prove that if $\bar{x}_{j}\in\arg\max(r)\cap\partial\Omega$ is a stationary point of $r(x)$ (i.e. $\nabla r(x)=0$), then

[TABLE]

This implies that, all other things being equal, the weight in front of a Dirac mass centred at a boundary point will be half that of a Dirac mass centred at an interior point. On the other hand, if $\bar{x}_{j}\in\arg\max(r)\cap\partial\Omega$ is a nonstationary point of the restriction of $r(x)$ to the boundary, and there is at least one maximum point of $r(x)$ that belongs to the interior of $\Omega$, then the weight $a_{j}$ in front of $\delta_{\bar{x}_{j}}(x)$ will be zero (i.e. the mass in a neighbourhood of $\bar{x}_{j}$ will vanish as $t\to\infty$).

Remark 3**.**

The asymptotic relation (3.26) shows that the integral $\rho(t)$ of the solution to the IDE (1.1) complemented with (1.4) converges to $r_{M}$ more slowly than the solution of the related logistic ordinary differential equation

[TABLE]

which converges exponentially to $r_{M}$ as $t\to\infty$. Moreover, whilst

[TABLE]

which means that $N(t)$ will approach the asymptotic value $r_{M}$ from below if $N(0)<r_{M}$ and from above if $N(0)>r_{M}$, the asymptotic relation (3.26) indicates that $\rho(t)$ will always approach $r_{M}$ from below independently from the value of $\rho(0)$. This also implies that if $\rho(0)>r_{M}$ then the derivative of $\rho(t)$ will change sign at least once on $\mathbb{R}^{+}$.

It is also worth comparing the speed of convergence given by the asymptotic relation (3.26) with the one that can be obtained based on the best estimates currently available in the literature, namely

[TABLE]

Furthermore, this cannot vanish as $\displaystyle{\mathcal{O}\left(\frac{1}{t^{\alpha}}\right)}$ for any $\alpha>1$ [31, 51].

3.2 Numerical solutions

To confirm the asymptotic results established by Theorem 2, we solve numerically the Cauchy problem (2.5). In particular, we approximate the IDE (1.1) complemented with (1.4) using the forward Euler method with step size $0.01$. We select a uniform discretisation of the interval $\Omega:=[-1,2]$ consisting of $1000$ points as the computational domain of the independent variable $x$, and we consider $t\in[0,200]$. All numerical computations are performed in Matlab.

We choose the initial condition

[TABLE]

and we use the following definition

[TABLE]

which satisfies the assumptions of Theorem 2. As shown by the plot in Figure 1, the function $r(x)$ defined according to (3.38) has two maximum points, that is, $\bar{x}_{1}\in[-1,0]$ and $\bar{x}_{2}\in[0,2]$.

We compute numerically the integrals

[TABLE]

and the coefficients $a_{1}$ and $a_{2}$ given by (3.25). In the one-dimensional setting considered here, the expressions given by (3.25) read as

[TABLE]

The results obtained are summarised in Figure 2 and Figure 3. As we would expect based on Theorem 1, the numerical results displayed in Figure 2 show that $n(t,x)$ becomes concentrated as a sum of two Dirac masses centred at the points $\bar{x}_{1}$ and $\bar{x}_{2}$ (left panel), while $\rho(t)$ converges to $r_{M}$ (right panel).

Furthermore, the curves displayed in the left panel of Figure 3 show that, in agreement with the results established by Theorem 2, the integrals $\rho_{1}(t)$ and $\rho_{2}(t)$ defined according to (3.39) converge, respectively, to the values $a_{1}\,r_{M}$ and $a_{2}\,r_{M}$, with $a_{1}$ and $a_{2}$ given by (3.40), while the curves in the right panel of Figure 3 show that $\displaystyle{\left(r_{M}-\rho(t)\right)t\to\frac{1}{2}}$ as $t\to\infty$, i.e. the asymptotic relation (3.26) is verified.

4 Long-time behaviour of the initial-boundary value problem (2.6)

In this section, we study the asymptotic behaviour of the solutions to the initial-boundary value problem (2.6) as $t\to\infty$ (Section 4.1), and we provide a sample of numerical solutions that confirm the analytical results obtained (Section 4.2).

4.1 Asymptotic analysis

Under assumptions (2.7) there exists a unique non-negative classical solution $n_{\beta}\in C([0,+\infty),L^{1}(\Omega))\cap C^{1}((0,+\infty),C^{2,\alpha}(\Omega))$ of the initial-boundary value problem (2.6) [19]. Moreover, the behaviour of $n_{\beta}(t,x)$ in the asymptotic regime $t\to\infty$ is known to be governed by the principal eigenpair $(\lambda_{\beta},\psi_{\beta})$ of the elliptic differential operator $L_{\beta}$ (i.e. $L_{\beta}\psi_{\beta}=-\lambda_{\beta}\psi_{\beta}$) defined according to (2.14) with zero Neumann boundary condition. The following theorem builds on the method of proof presented in [19, 33] and extends previous results by characterising the speed of convergence of $\rho_{\beta}$ towards its limit as $t\to\infty$.

Theorem 3**.**

Under assumptions (2.7), the solution of the initial-boundary value problem (2.6) is such that if $\lambda_{\beta}\geqslant 0$ then

[TABLE]

whereas if $\lambda_{\beta}<0$ then

[TABLE]

Furthermore, when $\lambda_{\beta}<0$, if

[TABLE]

where $\alpha_{2}(n^{0})$ is the the $L^{2}(\Omega)$-projection of $n^{0}$ onto the finite-dimensional eigenspace associated to $\lambda_{2,\beta}$, then there exist some real constants $K_{1}\neq 0$ and $K_{2}~{}\neq~{}0$ such that

[TABLE]

with $\gamma_{\beta}$ being the opposite to the spectral gap of the operator $L_{\beta}$, which is defined via (2.16).

Proof.

Let $u_{\beta}(t,x):=w(t)\,n_{\beta}(t,x)$ where $w(t)$ solves the Cauchy problem

[TABLE]

We have

[TABLE]

Using the fact that $w(t)$ solves the Cauchy problem (4.45) and $n_{\beta}(t,x)$ is the solution of the initial-boundary value problem (2.6), we obtain the following initial-boundary value problem for $u_{\beta}(t,x)$

[TABLE]

This is a standard parabolic problem, the solution of which is such that

[TABLE]

Since

[TABLE]

the above convergence result yields

[TABLE]

Hence the ordinary differential equation for $w(t)$ can be rewritten as

[TABLE]

with the function $f(t):=\rho_{\beta}(t)w(t)-K$ being such that $f(t)\to 0$ as $t\to\infty$. Using (4.47) we conclude that

[TABLE]

These asymptotic results along with the asymptotic result (4.46) ensure that

[TABLE]

Moreover, recalling that $u_{\beta}(t,x)=w(t)\,n_{\beta}(t,x)$ we find that if $\lambda_{\beta}<0$ then

[TABLE]

We now turn our attention to estimate (4.44). We recall that $\rho_{\beta}(t)=\frac{\int_{\Omega}u_{\beta}(t,x)\,{\rm d}x}{w(t)}$, and we estimate the numerator and the denominator separately.

For the numerator, we expand $u_{\beta}$ further in the orthonormal basis associated to the operator $L_{\beta}$ and integrate to find that there exists some constant $C$ such that

[TABLE]

If assumption (4.43) is satisfied then $C\neq 0$. The latter estimate for $\int_{\Omega}u_{\beta}(t,x)\,{\rm d}x$ gives also a more detailed characterisation of the behaviour of the function $f(t)$ in (4.47) when $t\to\infty$, that is,

[TABLE]

For the denominator, we note that solving (4.47) subject to the initial condition $w(0)=1$ gives

[TABLE]

This along with estimate (4.48) allows us to conclude that if $\gamma_{\beta}>\lambda_{\beta}$ then

[TABLE]

whereas if $\gamma_{\beta}\leqslant\lambda_{\beta}$ then

[TABLE]

for some constant $\tilde{C}\neq 0$.

Combining these estimates for $\int_{\Omega}u_{\beta}(t,x)\,{\rm d}x$ and $w(t)$ yields (4.44), and concludes the proof of Theorem 3. ∎

Remark 4**.**

Comparing the results established by Theorem 3 [cf. asymptotic relation (4.44)] with the results established by Theorem 2 [cf. asymptotic relation (3.26)], one can see that there is a clear difference between the speed of convergence of $\rho_{\beta}(t)$ to its long-time limit $-\lambda_{\beta}$ and the speed of convergence of $\rho(t)$ to its long time-limit $r_{M}$.

Remark 5**.**

Some results on the behaviour of $\gamma_{\beta}$ in the asymptotic regime $\beta\to 0$ are available in the existing literature. In particular, it is known that typically $\gamma_{\beta}\sim e^{-C/\sqrt{\beta}}$ if $\beta$ tends to [math] [53]. Hence, when $\beta$ is small, the exponent in (4.44) will be exponentially small.

The results established by Theorem 3 show that analysing the long-time behaviour of the solution to the initial-boundary value problem (2.6) comes down to studying the properties of the principal eigenpair $(\lambda_{\beta},\psi_{\beta})$. In particular, a general characterisation of $\lambda_{\beta}$ and $\psi_{\beta}(x)$ can be obtained when $\beta\to 0$. To illustrate this, we define $\beta:=\varepsilon$, where $\varepsilon$ is a small positive parameter, and study the behaviour of $\lambda_{\varepsilon}$ and $\psi_{\varepsilon}(x)$ in the asymptotic regime $\varepsilon\to 0$. Proposition 1 shows that the limit of $-\lambda_{\varepsilon}\psi_{\varepsilon}$ when $\varepsilon\to 0$ is given by a measure which has total mass equal to $r_{M}$ and consists of a weighted sum of Dirac masses centred at the points $\bar{x}_{1},\ldots,\bar{x}_{N}$. This kind of concentration result is standard in semiclassical analysis (see for instance [29]) but we give here a short proof that applies to our case for the sake of self-containedness.

Proposition 1**.**

If assumptions (2.7) and (2.8) are satisfied, then

[TABLE]

and, up to extraction of subsequences,

[TABLE]

Proof.

We divide the proof of Proposition 1 into two steps. We first prove claim (4.49) and then claim (4.50).

Step 1: proof of (4.49).

Since $\lambda_{\varepsilon}$ is the principal eigenvalue of the differential elliptic operator $L_{\varepsilon}$, we have [cf. equation (2.15)]

[TABLE]

We start by noting that

[TABLE]

and, therefore, $\lambda_{\varepsilon}\geqslant-r_{M}$ for all $\varepsilon>0$. Thus it suffices to show that $\lim\,\lambda_{\varepsilon}\leqslant-r_{M}$ as $\varepsilon\to 0$ in order to prove (4.49). To do this we construct a sequence of positive normalised $H^{1}$-functions $\phi_{\varepsilon}$ such that $\mathcal{R}(L_{\varepsilon},\phi_{\varepsilon})$ converges to $-r_{M}$ as $\varepsilon\to 0$. We introduce the function

[TABLE]

where $|\cdot|$ denotes the Euclidean norm on $\mathbb{R}^{d}$ and $C$ is a normalising constant such that $G$ has integral $1$. Recalling the classical result

[TABLE]

we choose $\displaystyle{\phi_{\varepsilon}^{2}:x\longmapsto\frac{1}{\varepsilon^{\frac{d}{4}}}G\left(\frac{x-\bar{x}_{i}}{\varepsilon^{\frac{1}{4}}}\right)}$ so that $\phi_{\varepsilon}^{2}\xrightharpoonup[\varepsilon\rightarrow 0]{}\delta_{\bar{x}_{i}}$. Using the fact that

[TABLE]

and

[TABLE]

we obtain

[TABLE]

This concludes the proof of (4.49).

Step 2: proof of (4.50).

The pair $(\lambda_{\varepsilon},\psi_{\varepsilon})$ satisfies the eigenvalue problem

[TABLE]

and integrating over $\Omega$ we find

[TABLE]

Hence,

[TABLE]

Finally, for any \varphi\in C\big{(}\overline{\Omega}\big{)} with $\operatorname{supp}(\varphi)\cap\arg\max(r)=\emptyset$ we have

[TABLE]

The latter integral inequality along with the asymptotic result (4.53) yields

[TABLE]

This concludes the proof of (4.50). ∎

Remark 6**.**

The results of Proposition 1 imply that if $\arg\max(r)=\left\{\bar{x}_{1}\right\}$ then the limit of $-\lambda_{\varepsilon}\psi_{\varepsilon}$ when $\varepsilon\to 0$ is given by the measure $r_{M}\,\delta_{\bar{x}_{1}}(x)$.

Remark 7**.**

The proof of Proposition 1 can be adapted to the case where the function $r(x)$ attains its maximum only on the boundary $\partial\Omega$. We expect that this can be done, as in Laplace’s method, by adjusting the normalising constant $C$ in (4.52), depending on the nature of the maximum at the boundary (stationary or not). However, we consider here only the simpler case corresponding to Proposition 1, which suffices for our purposes.

Remark 8**.**

Since $r_{M}>0$, expression (4.51) for the Rayleigh quotient is such that if $\varepsilon$ is small enough then $-\lambda_{\varepsilon}>0$. Hence, based on the result of Theorem 3, we have that $\rho_{\varepsilon}(t)$ will not vanish as $t\to\infty$ when $\varepsilon$ is sufficiently small.

In the framework of the results established by Proposition 1, to fully characterise the long-time limit of $n_{\varepsilon}(t,x)$ when $\varepsilon\to 0$ it is necessary to assess whether there exists a unique set of admissible coefficients $a_{1},\ldots,a_{N}$ (i.e. if the limit measure is unique); if so, one needs to identify the values of the coefficients that define the only admissible limit measure.

A case where we expect the limit measure to be unique is when the set $\Omega$ and the function $r(x)$ are symmetric with respect to the hyperplane $S$ defined according to (2.17). In this case, a complete characterisation of the limit measure is given by the following proposition.

Proposition 2**.**

Under assumptions (2.7) and (2.8), letting $N=2$ and making the additional symmetry assumptions (2.18) and (2.19), we have

[TABLE]

Proof.

Under the symmetry assumptions (2.18) and (2.19) the points $\bar{x}_{1}\in\Omega$ and $\bar{x}_{2}\in\Omega$ are symmetric with respect to the hyperplane $S$, i.e.

[TABLE]

Moreover, in the case where $N=2$, the result established by Proposition 1 implies that

[TABLE]

Finally, Lemma 2 ensures that $\psi_{\varepsilon}$ is symmetric with respect to the hyperplane $S$, and Lemma 4 in turn ensures that the weak limit point of $\psi_{\varepsilon}$ for $\varepsilon\to 0$ is symmetric with respect to the hyperplane $S$ as well. Hence, $a=\frac{1}{2}$ and there is a unique limit point given by

[TABLE]

Since the sequence $\psi_{\varepsilon}$ is bounded in $L^{1}(\Omega)$, the Banach-Alaoglu Theorem ensures that it is relatively (weakly-$*$) compact in $\mathcal{M}^{1}(\overline{\Omega})$. This along with the uniqueness of the limit point gives the convergence of the whole sequence $\psi_{\varepsilon}$, i.e.

[TABLE]

which concludes the proof of Proposition 2. ∎

Furthermore, an almost exhaustive characterisation of the limit measure in the absence of particular symmetries is provided by the following proposition, whereby the function $\zeta:\arg\max\left(r\right)\to\mathbb{R}^{+}$, which is defined as

[TABLE]

where $(\lambda_{j}^{i})_{1\leqslant j\leqslant d}$ are the eigenvalues of $H_{i}$ (each counted with its multiplicity), is used to characterise the concavity of the function $r(x)$ at the maximum points, as it was done in previous papers on semiclassical analysis [27, 28].

Proposition 3**.**

Under assumptions (2.7) and (2.8),

[TABLE]

In particular, if $\arg\min\left(\zeta\right)=\left\{\bar{x}_{m}\right\}$ for some $1\leqslant m\leqslant N$, then

[TABLE]

Proof.

We note that studying the asymptotic behaviour of $\psi_{\varepsilon}$ when $\varepsilon\to 0$ is equivalent to studying the asymptotic behaviour of the principal eigenfunction $\varphi_{\varepsilon}$ of the differential elliptic operator $\varepsilon\Delta-V$ with $V:=-r$. The result of Proposition 1 ensures that the support of the weak limit of $\varphi_{\varepsilon}$ as $\varepsilon\to 0$ will be a (possibly improper) subset of the set $\arg\min\left(V\right)=\left(\bar{x}_{1},\ldots,\bar{x}_{N}\right)$. Investigating at which points of this discrete set the weak limit point of the sequence $\varphi_{\varepsilon}$ will actually be concentrated is a fundamental question in semiclassical analysis. Such a question arises in the study of the dynamics of a particle confined within a region of space surrounding a minimum point of the potential $V$ (i.e. a potential well) in the asymptotic regime of small noise (i.e. when $\varepsilon\to 0$) [52, 27, 28, 29]. Recasting the problem in this way, we can use the asymptotic results presented in [27, 28] which ensure that, under assumptions (2.7) and (2.8), when $\varepsilon$ tends to [math], the principal eigenfunction $\varphi_{\varepsilon}$ concentrates on the set $\displaystyle{\arg\min(\zeta)}$, with $\zeta$ defined via (4.55). Hence, under the additional assumption that the set $\displaystyle{\arg\min(\zeta)}$ coincides with the singleton $\left\{\bar{x}_{m}\right\}$ for some $1\leqslant m\leqslant N$, we find that $\varphi_{\varepsilon}$ concentrates at the point $\bar{x}_{m}$ as $\varepsilon\to 0$, whence (4.57). ∎

Remark 9**.**

In the case where $d=1$, the assumption $\arg\min\left(\zeta\right)=\left\{\bar{x}_{m}\right\}$ reads as

[TABLE]

Remark 10**.**

Comparing the results of Proposition 3 with the results established by Theorem 2 one can see that there is a stark difference between the long-time behaviour of the solution to the initial-boundary value problem (2.6) for $\beta\to 0$ and the long-time behaviour of the solution to the Cauchy problem (2.5). In fact, in the case where the set $\arg\min\left(\zeta\right)$ is reduced to a singleton and $n^{0}(x)>0$ for all $x\in\arg\max(r)$, Theorem 2 shows that the long-time limit of $n(t,x)$ will be given by a sum of multiple Dirac masses with different positive weights, whereas Proposition 3 shows that the long-time limit of $n_{\beta}(t,x)$ for $\beta\to 0$ will consist of one single Dirac mass. From the point of view of evolutionary dynamics, this implies that, all else being equal, individuals in the phenotypic states $x\in\arg\max(r)$ will coexist in the absence of phenotypic changes, whereas only individuals in the phenotypic state $\bar{x}_{m}$ will ultimately survive when heritable phenotypic changes occur. This provides a mathematical formalisation of the idea that, while being historically assumed to play a neutral role in evolutionary outcome, heritable phenotypic changes can shape the equilibrium phenotypic distribution of asexual populations with multi-peaked fitness landscapes, even if there is no bias in the generation of novel phenotypic variants.

Remark 11**.**

We remark that when the set $\displaystyle{\arg\min\left(\zeta\right)}$ is not a singleton, it is still possible to go further in reducing the support of the limit point of the sequence $\psi_{\varepsilon}$ as $\varepsilon\to 0$. However, the conditions determining which of the coefficients $a_{i}$ will be different from zero become rather convoluted, as shown by the results of semiclassical analysis presented in [29]. Therefore, we consider here only the simpler case corresponding to Proposition 3, which suffices for our purposes.

4.2 Numerical solutions

To confirm the asymptotic results established by Propositions 1-3, we solve numerically the initial-boundary value problem (2.6). Numerical solutions are constructed by approximating the diffusion term via a second-order central difference scheme [34] and then using the forward Euler method with step size $0.01$ to approximate the resulting system of ordinary differential equations. We select a discretisation of the interval $\Omega:=[-1,2]$ consisting of $1000$ points as the computational domain of the independent variable $x$ and let $t\in[0,t_{f}]$ with $t_{f}$ being either $200$ or $800$. All numerical computations are performed in Matlab.

We define $\beta=10^{-6}$, choose the initial condition (3.37), and use either the following definition

[TABLE]

or definition (3.38). Definition (4.59) satisfies the assumptions of Proposition 2 with $S~{}:=~{}\{0.5\}$ (cf. the plot in Figure 4), whereas definition (3.38) satisfies the assumptions of Proposition 3 and, as previously noted, it has two maximum points $\bar{x}_{1}\in[-1,0]$ and $\bar{x}_{2}\in[0,2]$ (cf. the plot in Figure 1).

We compute numerically the following integrals

[TABLE]

The results obtained are summarised in Figure 5 and Figure 6. As we would expect based on Proposition 1 and Proposition 2, the numerical results displayed in Figure 5 show that when $r(x)$ is defined according to (4.59) the solution $n_{\beta}(t,x)$ becomes concentrated as a sum of two Dirac masses centred at the points $\bar{x}_{1}$ and $\bar{x}_{2}$, the integral $\rho_{\beta}(t)$ converges to $r_{M}$, and the integrals $\rho_{1}(t)$ and $\rho_{2}(t)$ defined via (4.60) both converge to $\displaystyle{\frac{r_{M}}{2}}$.

On the other hand, the numerical results displayed in Figure 6 show that, in agreement with the results of Proposition 1 and Proposition 3, when $r(x)$ is defined according to (3.38) the integral $\rho_{\beta}(t)$ converges to $r_{M}$ while the solution $n_{\beta}(t,x)$ becomes concentrated as one single Dirac mass centred at the point $\bar{x}_{2}$ (left panel), which is the maximum point of the function $r(x)$ that satisfies condition (4.58) – i.e. $|r^{\prime\prime}(\bar{x}_{2})|<|r^{\prime\prime}(\bar{x}_{1})|$. As a consequence, the integral $\rho_{1}(t)$ converges to zero, whereas the integral $\rho_{2}(t)$ converges to $r_{M}$.

5 Research perspectives

There are several possible generalisations of the prototypical selection model (1.1) and selection-mutation model (1.3) for which suitable developments of the methods used here would be relevant.

More general saturating non-local terms.

A natural way to extend our study would be to replace the saturating term $\rho(t)$ with a more general non-local term of the form $\int_{\Omega}K(x,y)\,n(t,y)\,{\rm d}y$, where the kernel $K(x,y)$ models the effect of competitive interactions between individuals in the phenotypic state $x$ and other individuals in a generic phenotypic state $y$. While the long-time behaviour of the IDE model with such a more general saturating non-local term was extensively studied in [31], where the convergence of the solution to a weighted sum of Dirac masses was also investigated depending on the properties of the kernel $K$, the existing literature still lacks a precise characterisation of the long-time behaviour of the solutions for the corresponding PDE model, with the exception of the particular case when $K(x,y)\equiv k(y)$ or similar cases [19]. While we expect that extending our results to this particular case would be relatively easy, the case of a generic kernel $K(x,y)$ is an open problem that requires a different approach compared to the one undertaken here.

Integral kernel modelling phenotypic changes.

Our results could be extended to the case where the linear diffusion operator in the non-local PDE (1.3) is replaced by an integral term of the form $\int_{\Omega}\left(M(x,y)\,n(t,y)-M(y,x)\,n(t,x)\right)\,{\rm d}y$, where the kernel $M(x,y)$ models the transition of individuals from a generic phenotypic state $y$ to the phenotypic state $x$. In [6] it was shown that, when phenotypic changes are modelled through such an integral kernel, the solution of selection-mutation models like the one considered here will typically converge to a measure as $t\to\infty$, and a criterion was derived to determine whether the limit measure would be singular or absolutely continuous. Suitable developments of our methods would make it possible to investigate the dependence of such a criterion on the weight of mutations compared to selection, which would be captured by a scaling parameter analogous to our parameter $\varepsilon$.

Systems of equations.

It would also be interesting to extend our results to the case of systems of IDEs of the form of (1.1) and systems of non-local PDEs of the form of (1.3). In this regard, the results presented in [51] for a specific system of IDEs could prove useful, since they establish the convergence of the solution to a measure as $t\to\infty$ and provide a characterisation of its support. As for systems of corresponding non-local PDEs, the convergence of the components of the solution to the principal eigenfunctions of the related elliptic differential operators when $t\to\infty$ has been proved for a two-by-two competitive system [33]. Apart from these particular cases, the long-time behaviour of the solutions of these systems of IDEs and non-local PDEs is still an open problem, which requires a different approach compared to the one undertaken here.

Acknowledgments

The authors would like to thank Sepideh Mirrahimi for the interesting discussions, and gratefully acknowledge Bernard Helffer for a fruitful exchange of emails regarding his semiclassical analysis results obtained in collaboration with Johannes Sjöstrand [27, 28]. C.P. acknowledges support from the Swedish Foundation of Strategic Research grant AM13-004. T.L. acknowledges support of the project PICS-CNRS no. 07688.

Appendix. Proof of Theorem 1 under assumption (3.23)

We start by noting that $\rho(t)$ is uniformly bounded in $L^{\infty}([0,\infty))$. In fact, integrating both sides of the IDE for $n(t,x)$ over $\Omega$ and using Grönwall’s lemma yields

[TABLE]

and

[TABLE]

Then we prove that $\rho\in BV([0,+\infty))$. In order to do this, we define

[TABLE]

so that differentiating we obtain

[TABLE]

Multiplying both sides of the latter differential equation by \displaystyle{-\big{(}{\rm sgn}(q)\big{)}_{-}} and estimating the right-hand side of the resulting differential equation from above we find that

[TABLE]

Moreover, for any $T>0$ we have

[TABLE]

Using estimates (5.63) and (5.64) we obtain

[TABLE]

and letting $T\to\infty$ we find

[TABLE]

This estimate along with the fact that $\rho\in L^{\infty}([0,\infty))$ ensures that $\rho\in BV([0,+\infty))$.

Since $\rho\in BV([0,+\infty))$ we conclude that $\rho(t)$ admits a limit $\rho^{\infty}$ as $t\to\infty$. The fact that $\rho^{\infty}=r_{M}$ can be proved via contradiction. Suppose that $\rho^{\infty}<r_{M}$ and consider $\varepsilon>0$ such that $r(x)>r_{M}-\varepsilon$ for all $x\in B(x_{i},\varepsilon)$, where $B(x_{i},\varepsilon)$ is the ball of centre $\bar{x}_{i}\in\arg\max(r)$ and radius $\varepsilon$. Since $\rho(t)\to\rho^{\infty}$ as $t\to\infty$, if $\rho^{\infty}<r_{M}$ then for $\varepsilon$ small enough there exists $\tau_{\varepsilon}>0$ such that $\rho(t)<r_{M}-2\varepsilon$ for all $t\geqslant\tau_{\varepsilon}$. Solving the IDE (1.1) complemented with (1.4) for $t\geqslant\tau_{\varepsilon}$ gives

[TABLE]

Integrating both sides of (5.65) over $\Omega$ and estimating from below we find

[TABLE]

which implies that $\rho(t)\to\infty$ as $t\to\infty$. Thus we arrive at a contradiction. Now suppose that $\rho^{\infty}>r_{M}$. If so, there exist $\varepsilon>0$ sufficiently small and $\tau_{\varepsilon}>0$ sufficiently large so that $\rho(t)>r_{M}+\varepsilon$ for all $t\geqslant\tau_{\varepsilon}$. Solving the IDE (1.1) complemented with (1.4) for $t\geqslant\tau_{\varepsilon}$ gives (5.65). Moreover, integrating both sides of (5.65) over $\Omega$ and estimating from above yields

[TABLE]

which implies that $\rho(t)\to 0$ as $t\to\infty$. Thus we arrive again at a contradiction. In so doing we have proved that $\rho^{\infty}=r_{M}$.

Since the sequence $(n(t,\cdot))_{t>0}$ is bounded in $L^{1}(\Omega)$, the Banach-Alaoglu Theorem ensures that it is relatively (weakly-$*$) compact in $\mathcal{M}^{1}(\overline{\Omega})$. Thus we can extract a subsequence $n(t_{k},\cdot)\in\mathcal{M}^{1}(\overline{\Omega})$ such that

[TABLE]

where the measure $n^{\infty}$ is non-negative and its total mass is $r_{M}$. Finally, since solving the Cauchy problem (2.5) yields

[TABLE]

and $\rho(t)\to r_{M}$ as $t\to\infty$, we have

[TABLE]

which implies that

[TABLE]

This concludes the proof of Theorem 1 in the case where assumption (3.23) is satisfied.