Adaptation and migration of a population between patches

# Adaptation and migration of a population between patches

Sepideh Mirrahimi CMAP, Ecole Polytechnique, CNRS, INRIA. Route de Saclay, 91128 Palaiseau Cedex, France. Email: mirrahimi@cmap.polytechnique.fr.
August 23, 2019
###### Abstract

A Hamilton-Jacobi formulation has been established previously for phenotypically structured population models where the solution concentrates as Dirac masses in the limit of small diffusion. Is it possible to extend this approach to spatial models? Are the limiting solutions still in the form of sums of Dirac masses? Does the presence of several habitats lead to polymorphic situations?

We study the stationary solutions of a structured population model, while the population is structured by continuous phenotypical traits and discrete positions in space. The growth term varies from one habitable zone to another, for instance because of a change in the temperature. The individuals can migrate from one zone to another with a constant rate. The mathematical modeling of this problem, considering mutations between phenotypical traits and competitive interaction of individuals within each zone via a single resource, leads to a system of coupled parabolic integro-differential equations. We study the asymptotic behavior of the stationary solutions to this model in the limit of small mutations. The limit, which is a sum of Dirac masses, can be described with the help of an effective Hamiltonian. The presence of migration can modify the dominant traits and lead to polymorphic situations.

Key-Words: Structured populations, phenotypical and spatial structure, Hamilton-Jacobi equation, viscosity solutions, Dirac concentrations, stationary solutions

AMS Class. No: 35B25, 47G20, 49L25, 92D15

## 1 Introduction

Non-local Lotka-Volterra equations arise in models of adaptive evolution of phenotypically structured populations. These equations have the property that the solutions concentrate generally, in the limit of small diffusion, on several isolated points, corresponding to distinct traits. Can we generalize these models by adding a spatial structure? How do the dominant traits evolve if we introduce a new habitat? To understand the interaction of ecological and evolutionary processes in population dynamics, spatial structure of the communities and adaptation of species to the environmental changes, it is crucial to dispose mathematical models that describe them jointly. We refer to [19] and the references therein for general literature on the subject. In this manuscript we consider a model where several distinct favorable habitable zones are possible. Population dynamics models structured by spatial patches have been studied using both deterministic and probabilistic methods (see for instance [30, 1]). Our model, in the case of two patches, is indeed very close to the one studied in [30] where the authors use adaptive dynamics theory (adaptive dynamics is a theory, based on dynamical systems and their stability, to study population dynamics [11]). Here we model similar phenomena, by adding a spatial structure to an earlier known integro-differential model describing the darwinian evolution. Integro-differential models have the advantage that the mutations can be considered directly in the model without assuming a separation of time scales of evolution and ecology. The present work provides a general description of the asymptotic stationary solutions, in the general case where two or several patches are possible.

We study the asymptotic behavior of solutions of a system of coupled elliptic integro-differential equations with small diffusion terms. These solutions are the stationary solutions to a parabolic system describing the dynamics of a population density. The individuals are characterized by phenotypical traits, that we denote by . They can move between two or several patches, which are favorable habitable zones, with constant rates (that we denote by and in the case of two patches). The mathematical modeling is based on the darwinian evolution and takes into account mutations and competition between the traits. There is a large literature for mathematical modeling and analysis on the subject of adaptive evolution, we refer the interested reader to [16, 15, 11, 12, 22, 10]. Here, we represent the birth and death term by a net growth term that is different in each patch, for instance because of a change in the temperature, and depends on the integral parameter , which corresponds to the the pressure exerted by the whole population within patch on the resource. To model the mutations, we use Laplace terms with a small rate that is introduced to consider only rare mutations. We study the asymptotic behavior of stationary solutions as the mutation rate goes to . The asymptotic solutions are generally concentrated on one or several Dirac masses. We describe the position and the weight of these Dirac masses using a Hamilton-Jacobi approach.

The time-dependent model, in the case of two patches, is written as

 ⎧⎪ ⎪⎨⎪ ⎪⎩∂tn1ε−εΔn1ε=1εn1εR1(x,I1ε)+1εν2n2ε−1εν1n1ε,x∈Rd,∂tn2ε−εΔn2ε=1εn2εR2(x,I2ε)+1εν1n1ε−1εν2n2ε, (1)

with

 Iiε=∫ψi(x)niε(x)dx,for i=1,2. (2)

Such models, without the structure in space, have been derived from stochastic individual based models in the limit of large populations (see [7, 6]). This manuscript follows earlier works on parabolic Lotka-Volterra type equations to study concentration effects in models of phenotypically structured populations, that are based on a Hamilton-Jacobi formulation (see [12, 25, 2, 20]). The novelty of our work is that we add a spatial structure to the model by considering a finite number of favorable habitable zones. We thus have a system instead of a single equation. A Hamilton-Jacobi approach in the case of systems has also been introduced in [5] for an age structured model. See also [4] for a study of stationary solutions of the latter system. The Hamilton-Jacobi approach can also be used in problems other than adaptive evolution to prove concentration phenomena. See for instance [27, 26, 24] where related methods have been used to study the motion of motor proteins.

We are interested in the equilibria of 1 limited to a bounded domain, that are given by solutions of the following system

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩−ε2Δn1ε=n1εR1(x,I1ε)+ν2n2ε−ν1n1εin BL(0),−ε2Δn2ε=n2εR2(x,I2ε)+ν1n1ε−ν2n2εin BL(0),∇niε⋅→n=0in ∂BL(0) % and for i=1,2, (3)

where is a ball of radius with center in and is the unit normal vector, at the point , to the boundary of . The Neumann boundary condition is a way to express that mutants cannot be born in .

To formulate our results we introduce the assumptions we will be using throughout the paper. We assume that, there exist positive constants and such that

 ψ1=ψ2=ψ,am≤ψ(x)≤aM,∥ψ(x)∥W2,∞≤Aand∇ψ⋅→n=0 in ∂BL(0). (4)

Moreover there exist positive constants , , and such that, for all and ,

 (5)
 −C≤∂Ri∂I(x,I)≤−1C, (6)
 −D≤Riξξ(x,I),for x∈BL(0), I∈[Im,IM]% , ξ∈Rd, |ξ|=1 and i=1,2. (7)

We use the Hopf-Cole transformation

 niε=exp(uiεε),% for i=1,2, (8)

and replace the latter in the system satisfied by to obtain

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩−εΔu1ε=|∇u1ε|2+R1(x,I1ε)+ν2exp(u2ε−u1εε)−ν1,in BL(0),−εΔu2ε=|∇u2ε|2+R2(x,I2ε)+ν1exp(u1ε−u2εε)−ν2in BL(0),∇uiε⋅→n=0in ∂BL(0)and for i=1,2. (9)

We prove the following

###### Theorem 1.1

Assume 46. Then, as along subsequences, both sequences and converge uniformly in to a continuous function and goes to , with such that is a viscosity solution to the following equation

 ⎧⎪ ⎪⎨⎪ ⎪⎩−|∇u|2=H(x,I1,I2),in BL(0),maxx∈BL(0)u(x)=0, (10)

with

 H(x,I1,I2) the largest eigenvalue of the % matrixA=(R1(x,I1)−ν1ν2ν1R2(x,I2)−ν2). (11)

The function is indeed an effective Hamiltonian that contains information from the two patches and helps us in Theorem 1.2 to describe the support of the weak limits of as . We can interpret as the fitness of the system in the limit of (see [23] for the definition of fitness).

The difficulty here is to find appropriate regularity estimates on , that we obtain using the Harnack inequality [3] and the Bernstein method [9]. To prove convergence to the Hamilton-Jacobi equation, we are inspired from the method of perturbed test functions in homogenization [14].

The above information on the limit of allows us to describe the limit of the densities as vanishes. We prove

###### Theorem 1.2

Assume 47. Consider a subsequence such that and converge uniformly to and goes to , as , with solution of 10. Let , for , converge weakly in the sense of measures to along this subsequence. We have

 suppni⊂Ω∩Γ,for i=1,2, (12)

with

 Ω={x∈BL(0)|u(x)=0},Γ={x∈BL(0)|H(x,I1,I2)=maxx∈BL(0)H(x,I1,I2)=0}. (13)

Moreover, we have

 (R1(x)−ν1)n1(x)+ν2n2(x)=0,(R2(x)−ν2)n2(x)+ν1n1(x)=0,in BL(0) (14)

in the sense of distributions. The above condition is coupled by

 ∫BL(0)ψi(x)ni(x)=Ii. (15)

Theorem 1.2 provides us with a set of algebraic constraints on the limit, which allows us to describe the latter. In particular, if the support of , for , is a set of distinct points: , 14 implies that

 ni=k∑j=1ρijδ(x−xj),for i=1,2, (16)

with

 ρ2j=ρ1j(ν1−R1(xj,I1)ν2)=ρ1j(ν1ν2−R2(xj,I2)). (17)

Furthermore, the weights satisfy the normalization condition

 k∑j=1ψi(xj)ρij=Ii, for i=1,2. (18)

Condition 17 means that the vector is the eigenvector corresponding to the largest eigenvalue of the matrix at the point , which is . Thereby 17 implies once again that .

We point out that since , for , is such that the fitness vanishes on the support of and is negative outside the support, we can interpret as evolutionary stable distribution of the model. In adaptive dynamics, evolutionary stable distribution (ESD) corresponds to a distribution that remains stable after introduction of small mutants (see [21, 13, 17] for a more detailed definition). See also [10, 28] for related works on stability and convergence to ESD for trait-structured integro-differential models.

The set of assumptions in Theorem 1.2 allows us to describe the asymptotics of the stationary solutions, in the limit of rare or small mutations. In Section 5 we provide some examples where based on this information we can describe the asymptotics. In particular, we notice that the introduction of a new environment can lead to dimorphic situations. We refer to [8] for a related work using the Hamilton-Jacobi approach, where polymorphic situations can also appear in a model with multiple resources.

The paper is organized as follows. In Section 2 we prove some bounds on and some regularity properties on that allow us to pass to the limit as and derive the Hamilton-Jacobi equation with constraint. Theorem 1.1 is proved in Section 3. Using the results obtained on the asymptotic behavior of we prove Theorem 1.2 in Section 4. In Section 5 we provide some examples where the information given by Theorem 1.1 and Theorem 1.2 allows us to describe the limit. The asymptotic behavior of the stationary solutions in a more general framework, where more than two habitable zones are considered, is given in Section 6. Finally in Section 7 we present some numerical simulations for the time-dependent problem and compare them with the behavior of stationary solutions.

## 2 Regularity results

###### Lemma 2.1

Under assumptions 46 we have, for chosen small enough,

 Im≤Iiε≤IM,for i=1,2. (19)

In particular, along subsequences, converges to , with .

###### Remark 1

This is the only part, where we use Assumption 4. If is a solution of 3 such that 19 is satisfied, then the results of Theorems 1.1 and 1.2 hold true without necessarily assuming 4. In particular, one can take .

###### Proof 1

We prove the result by contradiction. We suppose that (the case with , and the inequalities from below can be treated following similar arguments). We multiply the first equation in 3 by , integrate, and use 4 to obtain

 −ε2AamI1ε≤∫ψ(x)n1ε(x)R1(x,I1ε)dx+ν2I2ε−ν1I1ε.

Using now 5, 6 and the fact that we deduce that, for small enough,

 0≤(δ−ε2Aam)I1ε≤ν2I2ε−ν1I1ε,

and thus

 ν1ν2IM≤I2ε.

Now we multiply the equations in 3 by , integrate and add them and use 4 to obtain

 −ε2Aam(I1ε+I2ε)≤∫ψ(x)n1ε(x)R1(x,I1ε)dx+∫ψ(x)n2ε(x)R2(x,I2ε)dx.

From 5, 6 and the above bounds on and it follows that

 −ε2Aam(I1ε+I2ε)≤−δ(I1ε+I2ε),

which is not possible if is small enough. We conclude that .

###### Theorem 2.2

Assume 46. Then
(i) there exists a positive constant , such that for ,

 |uiε(x)−ujε(y)|≤Dε,% for all x,y∈BL(0), |x−y|≤ε and i,j∈{1,2}. (20)

(ii) For and all , the family is uniformly Lipschitz and uniformly bounded from below.
(iii) For all , there exists such that for all ,

 uiε(x)≤a,for x∈BL(0) and i=1,2. (21)
###### Proof 2

(i) We define

 ˜niε(y)=niε(εy),for i=1,2.

From 3 we have

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩−Δ˜n1ε=˜n1εR1(εx,I1ε)+ν2˜n2ε−ν1˜n1εin BLε(0),−Δ˜n2ε=˜n2εR2(εx,I2ε)+ν1˜n1ε−ν2˜n2εin BLε(0), (22)

Moreover, from 5, 6 and 19 we have, for ,

 δ−C(IM−Im)≤R(εx,Iε)≤−δ+C(IM−Im).

Therefore the coefficients of the linear elliptic system 22 are bounded uniformly in . It follows from the classical Harnack inequality ([3], Theorem 8.2) that there exists a constant such that for all such that and for ,

 supz∈B1(y0)˜niε(z)≤Dinfz∈B1(y0)˜njε(z).

Rewriting the latter in terms of and and replacing by we obtain

 supz′∈Bε(x)niε(z′)≤Dinfz′∈Bε(y0)njε(z′),

and thus from 8 we deduce 20.

(ii) To prove the Lipschitz bounds, we use the Bernstein method (see [9]). We assume that

 maxx∈BL(0)(|∇u1ε(x)|,|∇u2ε(x)|)=|∇u1ε(xε)|, (23)

that is the maximum is achieved at a point and for (the case where the maximum is achieved for can be treated by similar arguments). From the Neumann boundary condition in 9 we know that is an interior point of . We define and notice that

 Δp=2Tr(Hessu1ε)2+2∇(Δu1ε)⋅∇u1ε.

We now differentiate the first equation in 9 with respect to and multiply it by to obtain

 −ε∇(Δu1ε)⋅∇u1ε=∇p⋅∇u1ε+∇R1⋅∇u1ε+ν2(∇u2ε−∇u1εε)⋅∇u1εexp(u2ε−u1εε).

From 23 we have

 (∇u2ε(xε)−∇u1ε(xε))⋅∇u1ε(xε)≤0,

and thus

 −ε2Δp(xε)+εTr(Hessu1ε(xε))2≤∇p(xε)⋅∇u1ε(xε)+∇R1(xε)⋅∇u1ε(xε).

Moreover from 23 we have and . It follows that

 ε(Δu1ε(xε))2≤εdTr(Hessu1ε(xε))2≤d∇R1(xε)⋅∇u1ε(xε).

Using again 9 we obtain

 (|∇u1ε|2+R1(xε,I1ε)+ν2exp(u2ε−u1εε)−ν1)2≤εd∇R1(xε,I1ε)⋅∇u1ε(xε).

From 5, 6 and 19 we find that is uniformly bounded for . We conclude that is uniformly Lipschitz for .

To prove uniform bounds from below, we notice from 4 and 19 that, for , there exists a point such that

 εln(ImaM|BL(0)|)≤uiε(¯¯¯xi).

From the latter and the Lipschitz bounds we obtain that

 −2LC1+εln(ImaM|BL(0)|)≤uiε,in BL(0) and for i=1,2.

It follows that the families are bounded from below for and .

(iii) We prove 21 for by contradiction. The proof for follows the same arguments. We assume that there exists a sequence such that as , and . Using the uniform Lipschitz bounds obtained in (ii) we have

 n1εk(x)>exp(a2εk),in [xk−a2C1,xk+a2C1]∩BL(0).

This is in contradiction with the bound from above in 19, for small enough. Therefore 21 holds.

## 3 Convergence to the Hamilton-Jacobi equation

In this section we prove Theorem 1.1.

###### Proof 3

Convergence to the Hamilton-Jacobi equation: From (ii) and (iii) in Theorem 2.2 we have that for , the families are uniformly bounded and Lipschitz. Therefore, from the Arzela-Ascoli Theorem we obtain that, along subsequences, and converge locally uniformly to some continuous functions , with . Moreover, from (i) in Theorem 2.2 we deduce that . Here we consider a subsequence of that converges to .

Let , be the largest eigenvalue of the matrix

 Aε=(R1(x,I1ε)−ν1ν2ν1R2(x,I2ε)−ν2),

and be the corresponding eigenvector. Since the non-diagonal terms in are strictly positive, using the Perron-Frobinius Theorem, we know that such eigenvalue exist and that and are strictly positive. We write

 ϕiε(x)=lnχiε(x),for i=1,2.

We prove that is a viscosity solution of

 −|∇u|2=H(x,I1,I2),in BL(0).

To this aim, suppose that has a maximum in . Then, we consider a sequence , such that as , and

 u1ε(xε)−φ(xε)−εϕ1ε(xε)=maxx∈BL(0)i=1,2uiε(x)−φ(x)−εϕiε(x)

is attained at the point and for (The case with can be treated similarly). In this case, we have in particular that

 u2ε(xε)−u1ε(xε)≤ε(ϕ2ε(xε)−ϕ1ε(xε)).

Using the latter and the viscosity criterion for the first equation in 9 we obtain that

 −ε(Δφ(xε)+εΔϕ1ε(xε))−|∇φ(xε)+ε∇ϕ1ε(xε)|2−R1(x,I1ε)−ν2exp(ϕ2ε(xε)−ϕ1ε(xε))+ν1≤0. (24)

We notice that, by definition of and , we have

 −R1(x,I1ε)−ν2exp(ϕ2ε(xε)−ϕ1ε(xε))+ν1=−H(x,I1ε,I2ε).

From the latter and by letting in 24 we deduce that

 −|∇φ(x)|2≤H(x,I1,I2),

and thus is a subsolution of 10 in the viscosity sense. The supersolution criterion can be proved in a similar way.

The constraint on the limit (): From 21 we obtain that . To prove that , we use the lower bounds on in 19. The proof of this property is classical and we refer to [2, 20] for a detailed proof.

## 4 Asymptotic behavior of stationary solutions

In this section we prove Theorem 1.2.

###### Proof 4

Support of : From 19, we deduce that, along subsequences and for , converges weakly to a measure . The fact that is a consequence of the Hopf-Cole transformation 8. To prove 12 it is enough to prove . To this aim following the idea in [25] we first prove that, for , are uniformly semi-convex. Recall that the smooth function is semiconvex with constant , if we have

 vξξ≥−C,for all |ξ|=1.

Let

 min{uiε,ξξ(x)|x∈BL(0),i=1,2,ξ∈Rd,|ξ|=1}=u1ε,ηη(xε). (25)

The case where the minimum is achieved for can be treated similarly. We differentiate twice the first equation in 9 with respect to and obtain

 −εΔu1ε,ηη=2∇u1ε⋅∇u1ε,ηη+2|∇u1ε,η|2+R1ηη+ν2⎛⎝(u2ε,η−u1ε,ηε)2+u2ε,ηη−u1ε,ηηε⎞⎠exp(u2ε−u1εε).

From 25 we obtain that , and . Using 7 It follows that

 |∇u1ε,η(xε)|2≤D2.

Since , we have . We deduce that

 |u1ε,ηη(xε)|2≤D2,

and thus

 min{uiε,ξξ(x)|x∈BL(0),i=1,2,ξ∈Rd,|ξ|=1}≥−√D2.

This proves that , for are semiconvex functions with constant . By passing to the limit in we obtain that is also semiconvex with the same constant.

A semiconvex function is differentiable at its maximum points. Therefore is differentiable with in the set . From 10, we deduce, that for all , , and thus . The fact that is immediate from 10 and the facts that is almost everywhere differentiable and is a continuous function.

Value of on the support: Let , i.e. is a smooth function with compact support in . We multiply 3 by and integrate with respect to in to obtain, for ,

 −ε2∫BL(0)niε(x)Δξ(x)dx=∫BL(0)ξ(x)niε(x)Ri(x,Iiε)dx−νi∫BL(0)ξ(x)niε(x)dx+νj∫BL(0)ξ(x)njε(x)dx.

Since weakly and , for , as , we obtain that, for ,

 ∫BL(0)ξ(x)ni(x)Ri(x,Ii)dx−νi∫BL(0)ξ(x)ni(x)dx+νj∫BL(0)ξ(x)nj(x)dx=0,

and thus 14. Finally, 15 follows from 2.

## 5 Examples of application

In 1618 we give a description of , assuming that the support of , for , is a set of distinct points, i.e. is a sum of Dirac masses and does not have a continuous distribution. This is what we expect naturally in the models based on darwinian evolution. More precisely, from Volterra-Gause s competitive exclusion principle (see [18, 29]) it is known in theoretical biology that in a model with limiting factors (as nutrients or geographic parameters) at most distinct species can generally survive. Here we have two limiting factors, represented by and , that correspond to the environmental pressures in the two patches. We thus expect to observe only monomorphic or dimorphic situations. This is also the case in the numerical simulations represented in Section 7.

From 12 we know that the support of is included in the set of maximum points of , , with limits of . If now is such that, for fixed , the corresponding set consists of isolated points, it follows that the supports of and consist also of isolated points. We give an example below where has clearly this property.

###### Example 5.1

(monomorphism towards dimorphism) Consider a case with the following values for the parameters of the system

 R1(x,I)=a1x2+b1x+c1−d1I,R2(x,I)=a2x2+b2x+c2−d2I, (26)

with

 ai,bi,ci,di∈R,ai<0

Then the supports of and consist at most of two single points.

We first notice that in the case where there is no migration between patches (), from the results in [20], we know that in patch , the population concentrates in large time on the maximum points of with the limit of . Since is a quadratic function in , it has a unique maximum and thus is a single Dirac mass on this maximum point. However, allowing migration by taking positive values for and the population can become dimorphic. In Section 7 we give a numerical example where a dimorphic situation appears (see Figure 2). This is in accordance with the competitive exclusion principle since we have introduced a new limiting factor, which is the choice of habitable zones.

Next, we prove the result:

###### Proof 5 (Proof of Example 5.1. )

From 12 we have that the stationary solutions concentrate asymptotically on the maximum points of defined as below

 H(x,I1,I2)=12F+12√F2−4G,

with

 F(x,I1,I2):=R1(x,I1)−ν1+R2(x,I2)−ν2,G(x,I1,I2):=(R1(x,I1)−ν1)(R2(x,I2)−ν2)−ν1ν2, (27)

Since , we deduce that

 minx∈BL(0)G(x,I1,I2)=0, (28)

and

 Γ={x∈BL(0)|H(x,I1,I2)=0}={x∈BL(0)|G(x,I1,I2)=0}. (29)

For fixed , is a polynomial of order . Therefore it has at most two maximum points. It follows that consists of one or two distinct points.

###### Example 5.2

(An asymmetric case) We assume that the parameters are such that the support of , for , consists of isolated points, and we have

 Ri(x,I)=Ri(x)−