During the last decade, an approach based on Hamilton–Jacobi equations with constraints has been developed to describe the asymptotic evolutionary dynamics of phenotypically structured populations, in the limit of vanishing mutations. Mathematical modeling of such phenomena leads to parabolic (or elliptic for the steady case) integrodifferential equations, whose solutions tend, as the diffusion term vanishes, toward a sum of Dirac masses, corresponding to dominant traits. These asymptotic solutions can be described using the Hamilton–Jacobi approach. There is a large literature on this method. We refer to [[4], [13], [17]] for the establishment of the basis of this approach for problems from evolutionary biology. Note that related tools were already used in the case of local equations (for instance KPP type equations) to describe the propagation phenomena (see, for instance, [[5], [8]]).
In almost all the previous works, the Hamilton–Jacobi approach has been used to describe the limit of the solution, corresponding to the population's phenotypical distribution, as the mutations steps vanish. However, from the biological point of view, it is sometimes more relevant to consider nonvanishing mutation steps. A recent work [[16]] has pointed out that such tools can also be used, for a simple model with homogeneous environment, to characterize the solution, while mutation steps are small, but nonzero. In this note, we show how such results can be obtained in a more complex situation with a heterogeneous environment.
Our purpose in this note is to study the solutions to the following system, for ,
(1){−ε2nε,1″(z)=nε,1(z)R1(z,Nε,1)+m2nε,2(z)−m1nε,1(z),−ε2nε,2″(z)=nε,2(z)R2(z,Nε,2)+m1nε,1(z)−m2nε,2(z),Nε,i=∫Rnε,i(z)dz,for i=1,2, with(2)Ri(z,Ni)=ri−gi(z−θi)2−κiNi,with θ1=−θ and θ2=θ. This system represents the equilibrium of a phenotypically structured population under mutation, selection and migration between two habitats. For more details on the modeling and the biological motivations, see Section 2.
Note that the asymptotic behavior, as and along subsequences, of the solutions to this system, under the assumption , for , and for bounded domains, was already studied in [[14]]. In the present work, we go further than the asymptotic limit along subsequences, and we obtain uniqueness of the limit and identify the dominant terms of the solution when ε is small but nonzero.
The main elements of the method: To describe the solutions , we use a WKB ansatz
nε,i(z)=12πεexp(uε,i(z)ε). Note that a first approximation, which is commonly used in the theory of ‘quantitative genetics’ (a theory in evolutionary biology that investigates the evolution of continuously varying traits, see [[18]], chapter 7), is a Gaussian distribution of the form:nε,i(z)=Ni2πεσexp(−(z−z⁎)2εσ2)=12πεexp(−12σ2(z−z⁎)2+εlogNiσε). Here, we try to go further than this a priori Gaussian assumption and to approximate directly . To this end, we write an expansion for in terms of ε :(3)uε,i=ui+εvi+ε2wi+O(ε3). We prove that is the unique viscosity solution to a Hamilton–Jacobi equation with constraint. The uniqueness of the viscosity solution to such Hamilton–Jacobi equation with constraint is related to the uniqueness of the Evolutionary Stable Strategy (ESS), see Section 3 for a definition and for the result on the uniqueness of the ESS, and to the weak KAM theory [[7]]. In section 4, we compute explicitly u , which indeed satisfiesmaxRu(z)=0, with the maximum points attained at one or two points corresponding to the ESS points of the problem. We then notice that, while , is exponentially small. Therefore, only the values of and at the points close to the zero level set of u matter, i.e. the ESS points. In section 5, we provide the main elements to compute formally and hence its secondorder Taylor expansion around the ESS points and the value of at those points. Then, we show, in section 6, that these approximations together with a fourthorder Taylor expansion of u around the ESS points are enough to approximate the moments of the population's distribution with an error in the order of .
The mathematical details of our results will be provided in [[12]]. The biological applications will be detailed in [[15]].
The solution to ((1)) corresponds to the steady solution to the following system, for ,
(4){∂tni(t,z)−ε2∂2∂z2ni(t,z)=ni(t,z)Ri(z,Ni(t))+mjnj(t,z)−mini(t,z),i=1,2,j=2,1,Ni(t)=∫Rni(t,z)dz,for i=1,2. This system represents the dynamics of a population that is structured by a phenotypical trait z , and lives in two habitats. We denote by the density of the phenotypical distribution in habitat i , and by the total population's size in habitat i . The growth rate is given by ((2)), where represents the maximum intrinsic growth rate, is the strength of the selection, is the optimal trait in habitat i , and represents the intensity of the competition. The constants are the migration rates between the habitats. In this note we assume that there is positive migration rate in both directions, i.e.(5)mi>0,i=1, 2. However, the source and sink case, where for instance , can also be analyzed using similar tools. We refer to [[15]] for the analysis of this case. We additionally assume that(6)max(r1−m1,r2−m2)>0. This guarantees that the population does not get extinct.
Such phenomena have already been studied by several approaches. A first class of results are based on the adaptive dynamics approach, where one considers that the mutations are very rare, so that the population has time to attain its equilibrium between two mutations and hence the population's distribution has discrete support (one or two points in a twohabitat model) [[2], [6], [11]]. A second class of results is based on an approach known as ‘quantitative genetics’, which allows more frequent mutations and does not separate the evolutionary and the ecological time scales. A main assumption in this class of works is that one considers that the population's distribution is a Gaussian [[9], [19]] or, to take into account the possibility of dimorphic populations, a sum of one or two Gaussian distributions [[3], [20]].
In our work, as in the quantitative genetics framework, we also consider continuous phenotypical distributions. However, we do not assume any a priori Gaussian assumption. We compute directly the population's distribution and in this way we correct the previous approximations. To this end, we also provide some results in the framework of adaptive dynamics and, in particular, we generalize previous results on the identification of the ESS to the case of nonsymmetric habitats. Furthermore, our work makes a connection between the two approaches of adaptive dynamics and quantitative genetics.


The adaptive dynamics framework 
In this section, we introduce some notions from the theory of adaptive dynamics that we will be using in the next sections [[11]]. We also provide our main result in this framework.
Effective fitness. The effective fitness is the largest eigenvalue of the following matrix:
(7)A(z;N1,N2)=(R1(z;N1)−m1m2m1R2(z;N2)−m2). This indeed corresponds to the effective growth rate associated with trait z in the whole metapopulation when the total population sizes are given by .
Demographic equilibrium. Consider a set of points . The demographic equilibrium corresponding to this set is given by , with the total population sizes , such that
ni(z)=∑j=1mαi,jδ(z−zj),Ni=∑j=1mαi,j,W(zj,N1,N2)=0, and such that is the right eigenvector associated with the largest eigenvalue of .
Evolutionary stable strategy. A set of points is called an evolutionary stable strategy (ESS) if
W(z,N1⁎,N2⁎)=0,for z∈A and,W(z,N1⁎,N2⁎)≤0,for z∉A, where are the total population sizes corresponding to the demographic equilibrium associated with the set .
Since, there are only two habitats for which we expect that at most two distinct traits coexist at the evolutionary stable equilibrium. We prove indeed the following.
Assume ((5))–((6)). There exists a unique set of points which is an evolutionary stable strategy. Such set has at most two elements.
We call an evolutionary stable strategy which has one (respectively two) element(s), a monomorphic (respectively dimorphic) ESS. We can indeed give a criterion to have monomorphic or dimorphic ESS, and we can identify the dimorphic ESS in the general case (see [[12]] for more details).


How to compute the zeroorder terms 
The identification of the zeroorder terms is based on the following result. Note that the part (ii) of the theorem below is a variant of Theorem 1.1 in [[14]].
Assume ((5))–((6)).
(i) As , converges to , the demographic equilibrium of the unique ESS of the model. Moreover, as , converges to , the total population size in patch i corresponding to this demographic equilibrium.
(ii) As , both sequences , for , converge along subsequences and locally uniformly in to a continuous function , such that u is a viscosity solution to the following equation
(8){−u′2=W(z,N1⁎,N2⁎),in R,maxz∈Ru(z)=0. Moreover, we have the following condition on the zerolevel set of u: suppn1⁎=suppn2⁎⊂{zu(z)=0}⊂{zW(z,N1⁎,N2⁎)=0}. (iii) There exists constants , for , which can be determined explicitly from , , , , , and θ, such that, under the condition (9)r2≠λ1r1+ν1,r1≠λ2r2+ν2, we have (10)suppn1⁎=suppn2⁎={zu(z)=0}={zW(z,N1⁎,N2⁎)=0}. The solution to ((8))–((10)) is unique, and hence the whole sequence converges locally uniformly in to u.
Note that a Hamilton–Jacobi equation of type ((8)) in general might admit several viscosity solutions. Here, the uniqueness is obtained thanks to ((10)) and a property from the weak KAM theory, which is the fact that viscosity solutions are completely determined by one value taken on each static class of the Aubry set ([[10]], Chapter 5 and [[1]]). In what follows, we assume that ((9)) and hence ((10)) always hold. We then give an explicit formula for u considering two cases.
(i) Monomorphic ESS. We consider the case where there exists a unique monomorphic ESS and where the corresponding demographic equilibrium is given by . Then u is given by
(11)u(z)=−∫z⁎z−W(x;N1⁎,N2⁎)dx.
(ii) Dimorphic ESS. We next consider the case where there exists a unique dimorphic ESS with the demographic equilibrium: , and . Then u is given by
u(z)=max(−∫za⁎z−W(x;N1⁎,N2⁎)dx,−∫zb⁎z−W(x;N1⁎,N2⁎)dx).


How to compute the nextorder terms 
In this section, we give the main elements to compute formally and the value of at the ESS point, with and the correctors introduced by ((3)), in the case of a monomorphic population. For the details of the computations for both monomorphic and dimorphic populations, we refer the interested reader to [[12]].
We consider the case of monomorphic population where the demographic equilibrium corresponding to the monomorphic ESS is given by . One can compute, using ((11)), a Taylor expansion of order 4 around the ESS point :
u(z)=−A2(z−z⁎)2+B(z−z⁎)3+C(z−z⁎)4+O(z−z⁎)5. To provide an approximation of the moments of the population's distribution, we have to compute constants , and such thatvi(z)=vi(z⁎)+Di(z−z⁎)+Ei(z−z⁎)2+O(z−z⁎)3,wi(z⁎)=Fi. A first element of the computations is obtained by replacing the functions u , and by the above approximations to compute . This leads tovi(z⁎)=log(Ni⁎A),Nε,i=Ni⁎+εKi⁎+O(ε2),with Ki⁎=Ni⁎(3CA2+EiA+Fi). Note also that writing ((1)) in terms of we obtain(12){−εuε,1″(z)=uε,1′2+R1(z,Nε,1)+m2exp(uε,2−uε,1ε)−m1,−εuε,2″(z)=uε,2′2+R2(z,Nε,2)+m1exp(uε,1−uε,2ε)−m2. A second element is obtained by keeping the zeroorder terms in the first line of ((12)) and using ((8)) to obtain(13)v2(z)−v1(z)=log(1m2(W(z,N1⁎,N2⁎)−R1(z,N1⁎)+m1)). The last element is derived from keeping the terms of order ε in ((12)), which leads to(14)−u″=2u′vi′−κiKi⁎+mjexp(vj−vi)(wj−wi),for {i,j}={1,2}. The functions and the coefficients , and can be computed by combining the above elements.


Approximation of the moments 
The above approximations of u , and around the ESS points allow us to estimate the moments of the population's distribution. In the monomorphic case these approximations are given below:
{Nε,i=∫nε,i(z)dz=Ni⁎(1+ε(Fi+EiA+3CA2))+O(ε2),με,i=1Nε,i∫znε,i(z)dz=z⁎+ε(3BA2+DiA)+O(ε2),σε,i2=1Nε,i∫(z−με,i)2nε,i(z)dz=εA+O(ε2),sε,i=1σε,i3Nε,i∫(z−με,i)3nε,i(z)dz=6BA32ε+O(ε32). One can obtain similar approximations in the case of dimorphic ESS. To compute the above integrals, replacing the approximation ((3)) in the integrals, a natural change of variable is to take . Therefore, each term can be considered as of order in the integration. This is why, to obtain a firstorder approximation of the integrals in terms of ε , it is enough to have a fourthorder approximation of , a secondorder approximation of , and a zeroorder approximation of , in terms of z around .
S. Mirrahimi has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 639638), and from the French ANR projects KIBORD ANR13BS010004 and MODEVOL ANR13JS010009.