Access to the PDF text

Free Article !

Comptes Rendus Mathématique
Volume 355, n° 2
pages 155-160 (février 2017)
Doi : 10.1016/j.crma.2016.12.001
Received : 11 October 2016 ;  accepted : 7 December 2016
A Hamilton–Jacobi method to describe the evolutionary equilibria in heterogeneous environments and with non-vanishing effects of mutations
Méthode de Hamilton–Jacobi pour décrire des équilibres évolutifs dans les environnements hétérogènes avec des mutations non évanescentes

Sylvain Gandon a , Sepideh Mirrahimi b
a Centre d'écologie fonctionnelle et évolutive (CEFE), UMR CNRS 5175, 34293 Montpellier cedex 5, France 
b CNRS, Institut de mathématiques (UMR CNRS 5219), Université Paul-Sabatier, 118, route de Narbonne, 31062 Toulouse cedex, France 


In this note, we characterize the solution to a system of elliptic integro-differential equations describing a phenotypically structured population subject to mutation, selection, and migration. Generalizing an approach based on the Hamilton–Jacobi equations, we identify the dominant terms of the solution when the mutation term is small (but nonzero). This method was initially used, for different problems arisen from evolutionary biology, to identify the asymptotic solutions, while the mutations vanish, as a sum of Dirac masses. A key point is a uniqueness property related to the weak KAM theory. This method allows us to go further than the Gaussian approximation commonly used by biologists, and is an attempt to fill the gap between the theories of adaptive dynamics and quantitative genetics.

The full text of this article is available in PDF format.

Dans cette note, nous étudions un système d'équations intégro-différentielles elliptiques, décrivant une population structurée par trait phénotypique soumise à des mutations, à la sélection et à des migrations. Nous généralisons une approche basée sur des équations de Hamilton–Jacobi pour détérminer les termes dominants de la solution lorsque les effets des mutations sont petits (mais non nuls). Cette méthode était initialement utilisée, pour différents problèmes venant de la biologie évolutive, pour identifier les solutions asymptotiques, lorsque les effets des mutations tendent vers 0, sous forme de sommes de masses de Dirac. Un point-clé est une propriété d'unicité en rapport avec la théorie de KAM faible. Cette méthode nous permet d'aller au-delà des approximations gaussiennes habituellement utilisées par les biologistes, et contribue ainsi à relier les théories de la dynamique adaptative et de la génétique quantitative.

The full text of this article is available in PDF format.

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) integro-differential 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 non-vanishing 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+εlog⁡Niσε). 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 satisfiesmaxR⁡u(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 second-order 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 fourth-order 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]].

Model and motivation

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 two-habitat 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.

Theorem 3.1

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 zero-order terms  

The identification of the zero-order 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]].

Theorem 4.1

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∈R⁡u(z)=0. Moreover, we have the following condition on the zero-level set of u: suppn1⁎=suppn2⁎⊂{z|u(z)=0}⊂{z|W(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⁎={z|u(z)=0}={z|W(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

(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

How to compute the next-order 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 zero-order 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 first-order approximation of the integrals in terms of ε , it is enough to have a fourth-order approximation of  , a second-order approximation of  , and a zero-order 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 ANR-13-BS01-0004 and MODEVOL ANR-13-JS01-0009.


Contreras G. Action potential and weak KAM solutions Calc. Var. Partial Differ. Equ. 2001 ;  13 (4) : 427-458 [cross-ref]
Day T. Competition and the effect of spatial resource heterogeneity on evolutionary diversification Amer. Nat. 2000 ;  155 (6) : 790-803 [cross-ref]
Débarre F., Ronce O., Gandon S. Quantifying the effects of migration and mutation on adaptation and demography in spatially heterogeneous environments J. Evol. Biol. 2013 ;  26 : 1185-1202
Diekmann O., Jabin P.-E., Mischler S., Perthame B. The dynamics of adaptation: an illuminating example and a Hamilton–Jacobi approach Theor. Popul. Biol. 2005 ;  67 (4) : 257-271 [cross-ref]
Evans L.C., Souganidis P.E. A PDE approach to geometric optics for certain semilinear parabolic equations Indiana Univ. Math. J. 1989 ;  38 (1) : 141-172 [cross-ref]
C. Fabre, S. Méléard, E. Porcher, C. Teplitsky, A. Robert, Evolution of a structured population in a heterogeneous environment, preprint.
Fathi A. Weak Kam Theorem in Lagrangian Dynamics  Cambridge, UK: Cambridge University Press (2016). 
Freidlin M. Geometric optics approach to reaction–diffusion equations SIAM J. Appl. Math. 1986 ;  46 : 222-232 [cross-ref]
Hendry A., Day T., Taylor E.B. Population mixing and the adaptive divergence of quantitative traits in discrete populations: a theoretical framework for empirical tests Evolution 2001 ;  55 (3) : 459-466 [cross-ref]
Lions P.-L. Generalized Solutions of Hamilton–Jacobi Equations  Boston, MA, USA: Pitman Advanced Publishing Program (1982). 
Meszéna G., Czibula I., Geritz S. Adaptive dynamics in a 2-patch environment: a toy model for allopatric and parapatric speciation J. Biol. Syst. 1997 ;  5 (02) : 265-284
S. Mirrahimi, A Hamilton–Jacobi approach to characterize the evolutionary equilibria in heterogeneous environments, forthcoming.
Mirrahimi S. Concentration Phenomena in PDEs from BiologyPhD thesis.   Paris-6: Université Pierre-et-Marie-Curie (2011). 
Mirrahimi S. Migration and adaptation of a population between patches Discrete Contin. Dyn. Syst., Ser. B 2013 ;  18 (3) : 753-768
S. Mirrahimi, S. Gandon, The equilibrium between selection, mutation and migration in spatially heterogeneous environments, forthcoming.
Mirrahimi S., Roquejoffre J.-M. Uniqueness in a class of Hamilton–Jacobi equations with constraints C. R. Acad. Sci. Paris, Ser. I 2015 ;  353 : 489-494 [inter-ref]
Perthame B., Barles G. Dirac concentrations in Lotka–Volterra parabolic PDEs Indiana Univ. Math. J. 2008 ;  57 (7) : 3275-3301
Rice S.H. Evolutionary Theory: Mathematical and Conceptual Foundations  : Sinauer Associates, Inc. (2004). 
Ronce O., Kirkpatrick M. When sources become sinks: migration meltdown in heterogeneous habitats Evolution 2001 ;  55 (8) : 1520-1531 [cross-ref]
Yeaman S., Guillaume F. Predicting adaptation under migration load: the role of genetic skew Evolution 2009 ;  63 (11) : 2926-2938 [cross-ref]

© 2016  Académie des sciences@@#104156@@
EM-CONSULTE.COM is registrered at the CNIL, déclaration n° 1286925.
As per the Law relating to information storage and personal integrity, you have the right to oppose (art 26 of that law), access (art 34 of that law) and rectify (art 36 of that law) your personal data. You may thus request that your data, should it be inaccurate, incomplete, unclear, outdated, not be used or stored, be corrected, clarified, updated or deleted.
Personal information regarding our website's visitors, including their identity, is confidential.
The owners of this website hereby guarantee to respect the legal confidentiality conditions, applicable in France, and not to disclose this data to third parties.
Article Outline