Adaptation and migration of a population between patches

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.

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 x∈Rd. They can move between two or several patches, which are favorable habitable zones, with constant rates (that we denote by ν1 and ν2 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 Ri(x,Ii) that is different in each patch, for instance because of a change in the temperature, and depends on the integral parameter Ii, which corresponds to the the pressure exerted by the whole population within patch i 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 0. 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.

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

where BL(p) is a ball of radius L with center in p and →n(x) is the unit normal vector, at the point x∈∂BL(0), to the boundary of BL(0). The Neumann boundary condition is a way to express that mutants cannot be born in Rd∖BL(0).

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

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

(4)

Moreover there exist positive constants Im, IM, δ and C such that, for all x∈BL(0) and i,j=1,2,

Theorem 1.1

Assume 4–6. Then, as ε→0 along subsequences, both sequences (u1ε)ε and (u2ε)ε converge uniformly in BL(0) to a continuous function u∈C(BL(0)) and (I1ε,I2ε) goes to (I1,I2), with (u,I1,I2) such that u 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 H 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 (n1ε,n2ε) as ε→0.
We can interpret H(x,I1,I2) as the fitness
of the system in the limit of ε→0 (see [23] for the definition of fitness).

The difficulty here is to find appropriate regularity estimates on uiε, 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 uiε allows us to describe the limit of the densities niε as ε vanishes. We prove

Theorem 1.2

Assume 4–7. Consider a subsequence such that u1ε and u2ε converge uniformly to u∈C(BL(0)) and (I1ε,I2ε) goes to (I1,I2), as ε→0, with (u,I1,I2) solution of 10. Let niε, for i=1,2, converge weakly in the sense of measures to ni 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 ni, for i=1,2, is a set of distinct points: suppni⊂{x1,x2,⋯,xk}, 14 implies that

Condition 17 means that the vector (ρ1jρ2j) is the eigenvector corresponding to the largest eigenvalue of the matrix A at the point xj, which is 0. Thereby 17 implies once again that suppni⊂Γ.

We point out that since ni, for i=1,2, is such that the fitness H vanishes on the support of ni and is negative outside the support, we can interpret ni 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 Iε and some regularity properties on uε that allow us to pass to the limit as ε→0 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 (uiε)ε 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.

Lemma 2.1

In particular, along subsequences, (I1ε,I2ε)ε converges to (I1,I2), with Im≤I1,I2≤IM.

Remark 1

This is the only part, where we use Assumption 4. If (n1ε,n2ε) 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 ψ1≢ψ2.

Proof 1

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

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

Using now 5, 6 and the fact that I1ε>IM we deduce that, for ε≤ε0 small enough,

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

and thus

ν1ν2IM≤I2ε.

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

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 D=D(C,Im,IM,δ,ν1,ν2) such that for all y0∈BLε(0) such that B1(y0)⊂BLε(0) and for i,j=1,2,

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

Rewriting the latter in terms of n1ε and n2ε and replacing (y0,z) by (xε,z′ε) we obtain

(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 xε∈BL(0) and for i=1 (the case where the maximum is achieved for i=2 can be treated by similar arguments). From the Neumann boundary condition in 9 we know that xε is an interior point of BL(0). We define p=|∇u1ε|2 and notice that

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

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

From 5, 6 and 19 we find that (R1(x,I1ε))ε is uniformly bounded for ε≤ε0. We conclude that (u1ε)ε is uniformly Lipschitz for ε≤ε0.

To prove uniform bounds from below, we notice from 4 and 19 that, for i=1,2, there exists a point ¯¯¯xi∈BL(0) 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 (uiε)ε are bounded from below for ε≤ε0 and i=1,2.

(iii) We prove 21 for i=1 by contradiction. The proof for i=2 follows the same arguments. We assume that there exists a sequence (εk,xk) such that εk→0 as k→∞, xk∈BL(0) and u1εk(xk)>a. 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 εk small enough. Therefore 21 holds.

Proof 3

Convergence to the Hamilton-Jacobi equation:
From (ii) and (iii) in Theorem 2.2 we have that for i=1,2, the families (uiε)ε are uniformly bounded and Lipschitz. Therefore, from the Arzela-Ascoli Theorem we obtain that, along subsequences, (u1ε)ε and (u2ε)ε converge locally uniformly to some continuous functions ui∈C(BL(0);R), with i=1,2. Moreover, from (i) in Theorem 2.2 we deduce that u1=u2. Here we consider a subsequence of (I1ε,I2ε,u1ε,u2ε)ε that converges to (I1,I2,u,u).

Let H(x,I1ε,I2ε), be the largest eigenvalue of the matrix

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

and (χ1ε(x)χ2ε(x)) be the corresponding eigenvector. Since the non-diagonal terms in Aε are strictly positive, using the Perron-Frobinius Theorem, we know that such eigenvalue exist and that χ1ε and χ2ε are strictly positive. We write

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

We prove that u is a viscosity solution of

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

To this aim, suppose that u−φ has a maximum in x∈BL(0). Then, we consider a sequence xε∈BL(0), such that as ε→0, xε→x and

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

is attained at the point xε and for i=1 (The case with i=2 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

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

The constraint on the limit (maxx∈BL(0)u(x)=0): From 21 we obtain that u(x)≤0. To prove that 0≤maxx∈BL(0)u(x), we use the lower bounds on Iiε in 19. The proof of this property is classical and we refer to [2, 20] for a detailed proof.

Proof 4

Support of ni: From 19, we deduce that, along subsequences and for i=1,2, (niε)ε converges weakly to a measure ni. The fact that suppni⊂Ω, for i=1,2, 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 i=1,2, uiε are uniformly semi-convex. Recall that the smooth function v is semiconvex with constant C, 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 i=2 can be treated similarly. We differentiate twice the first equation in 9 with respect to η and obtain

From 25 we obtain that Δu1ε,ηη(xε)≥0, ∇u1ε,ηη(xε)=0 and u2ε,ηη(xε)−u1ε,ηη(xε)≥0. Using 7 It follows that

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

Since u1ε,ηη=∇u1ε,η⋅η, we have |u1ε,ηη|≤|∇u1ε,η|. We deduce that

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

and thus

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

This proves that uiε, for i=1,2 are semiconvex functions with constant −√D2. By passing to the limit in ε→0 we obtain that u is also semiconvex with the same constant.

A semiconvex function is differentiable at its maximum points. Therefore u is differentiable with ∇u=0 in the set Ω. From 10, we deduce, that for all x∈Ω, H(x,I1,I2)=0, and thus Ω⊂{x∈BL(0)|H(x,I1,I2)=0}. The fact that maxx∈BL(0)H(x,I1,I2)=0 is immediate from 10 and the facts that u is almost everywhere differentiable and H(x,I1,I2) is a continuous function.

Value of ni on the support: Let ξ∈C∞c(BL(0)), i.e. ξ is a smooth function with compact support in BL(0). We multiply 3 by ξ and integrate with respect to x in BL(0) to obtain, for {i,j}={1,2},

In 16–18 we give a description of (n1,n2), assuming that the support of ni, for i=1,2, is a set of distinct points, i.e. ni 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 K limiting factors (as nutrients or geographic parameters) at most K distinct species can generally survive. Here we have two limiting factors, represented by I1 and I2, 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 ni is included in the set of maximum points of H(x,I1,I2), Γ, with (I1,I2) limits of (I1ε,I2ε). If now H is such that, for fixed (I1,I2), the corresponding set Γ consists of isolated points, it follows that the supports of n1 and n2 consist also of isolated points. We give an example below where H 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<di,for i=1,2.

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

We first notice that in the case where there is no migration between patches (ν1=ν2=0), from the results in [20], we know that in patch i, the population concentrates in large time on the maximum points of Ri(⋅,Ii) with Ii the limit of Iiε. Since Ri is a quadratic function in x, it has a unique maximum and thus ni is a single Dirac mass on this maximum point. However, allowing migration by taking positive values for ν1 and ν2 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.