A rigorous model study of the adaptative dynamics of Mendelian diploids.

A rigorous model study of the adaptative dynamics of Mendelian diploids.

Pierre Collet, Sylvie Méléard, Johan A.J. Metz CPHT, Ecole Polytechnique, CNRS UMR 7644, route de Saclay, 91128 Palaiseau Cedex-France; e-mail: collet@cpht.polytechnique.frCMAP, Ecole Polytechnique, CNRS, route de Saclay, 91128 Palaiseau Cedex-France; e-mail: sylvie.meleard@polytechnique.edu.Institute of Biology & Department of Mathematics, Leiden University, & NCB Naturalis, Leiden, Netherlands & Ecology and Evolution Program, Institute of Applied Systems Analysis, Laxenburg. Austria; e-mail: j.a.j.metz@biology.leidenuniv.nl
July 5, 2019
Abstract

Adaptive dynamics so far has been put on a rigorous footing only for clonal inheritance. We extend this to sexually reproducing diploids, although admittedly still under the restriction of an unstructured population with Lotka-Volterra-like dynamics and single locus genetics (as in Kimura’s 1965 infinite allele model). We prove under the usual smoothness assumptions, starting from a stochastic birth and death process model, that, when advantageous mutations are rare and mutational steps are not too large, the population behaves on the mutational time scale (the ’long’ time scale of the literature on the genetical foundations of ESS theory) as a jump process moving between homozygous states (the trait substitution sequence of the adaptive dynamics literature). Essential technical ingredients are a rigorous estimate for the probability of invasion in a dynamic diploid population, a rigorous, geometric singular perturbation theory based, invasion implies substitution theorem, and the use of the Skorohod topology to arrive at a functional convergence result. In the small mutational steps limit this process in turn gives rise to a differential equation in allele or in phenotype space of a type referred to in the adaptive dynamics literature as ’canonical equation’.


MSC 2000 subject classification: 92D25, 60J80, 37N25, 92D15, 60J75

Key-words: individual-based mutation-selection model, invasion fitness for diploid populations, adaptive dynamics, canonical equation, polymorphic evolution sequence, competitive Lotka-Volterra system.

1 Introduction

Adaptive dynamics (AD) aims at providing an ecology-based framework for scaling up from the micro-evolutionary process of gene substitutions to meso-evolutionary time scales and phenomena (also called long term evolution in papers on the foundations of ESS theory, that is, meso-evolutionary statics, cf Eshel (1983, in press); Eshel, Feldman, and Bergman (1998); Eshel and Feldman (2001)). One of the more interesting phenomena that AD has brought to light is the possibility of an emergence of phenotypic diversification at so-called branching points, without the need for a geographical substrate Metz et al. (1996); Geritz et al. (1998); Doebeli and Dieckmann (2000). This ecological tendency may in the sexual case induce sympatric speciation Dieckmann and Doebeli (1999). However, a population subject to mutation limitation and initially without variation stays essentially uni-modal, closely centered around a type that evolves continuously, as long as it does not get in the neighborhood of a branching point. In this paper we focus on the latter aspect of evolutionary trajectories.

AD was first developed, in the wake of Hofbauer and Sigmund (1987); Marrow, Law and Cannings (1992); Metz, Nisbet and Geritz (1992), as a systematic framework at a physicist level of rigor by Diekmann and Law Dieckmann and Law (1996) and by Metz and Geritz and various coworkers Metz, Nisbet and Geritz (1992); Metz et al. (1996); Geritz et al. (1998). The first two authors started from a Lotka-Volterra style birth and death process while the intent of the latter authors was more general, so far culminating in Durinx, Metz and Meszéna (2008). The details for general physiologically structured populations were worked out at a physicist level of rigor in Durinx, Metz and Meszéna (2008) while the theory was put on a rigorous mathematical footing by Champagnat and Méléard and coworkers Champagnat, Ferrière and Méléard (2008); Champagnat (2006); Méléard and Tran (2009), and recently also from a different perspective by Peter Jagers and coworkers Klebaner et al. (2011). All these papers deal only with clonal models. In the meantime a number of papers have appeared that deal on a heuristic basis with special models with Mendelian genetics (e.g. Kisdi and Geritz (1999); Van Dooren (1999, 2000); Van Doorn and Dieckmann (2006); Proulx and Phillips (2006); Peischl and Bürger (2008)), while the general biological underpinning for the ADs of Mendelian populations is described in Metz (in press). In the present paper we outline a mathematically rigorous approach along the path set out in Champagnat, Ferrière and Méléard (2008); Champagnat (2006), with proofs for those results that differ in some essential manner between the clonal and Mendelian cases. It should be mentioned though that just as in the special models in Kisdi and Geritz (1999); Van Dooren (1999, 2000); Proulx and Phillips (2006); Peischl and Bürger (2008) and contrary to the treatment in Metz (in press) we deal still only with the single locus infinite allele case (cf Kimura Kimura (1965)), while deferring the infinite loci case to a future occasion.

Our reference framework is a diploid population in which each individual’s ability to survive and reproduce depends only on a quantitative phenotypic trait determined by its genotype, represented by the types of two alleles on a single locus. Evolution of the trait distribution in the population results from three basic mechanisms: heredity, which transmits traits to new offsprings thus ensuring the extended existence of a trait distribution, mutation, generating novel variation in the trait values in the population, and selection acting on these trait values as a result of trait dependent differences in fertility and mortality. Selection is made frequency dependent by the competition of individuals for limited resources, in line with the general ecological spirit of AD. Our goal is to capture in a simple manner the interplay between these different mechanisms.

2 The Model

We consider a Mendelian population and a hereditary trait that is determined by the two alleles on but a single locus with many possible alleles (the infinite alleles model of Kimura Kimura (1965)). These alleles are characterized by an allelic trait . Each individual is thus characterized by its two allelic trait values , hereafter referred to as its genotype, with corresponding phenotype , with . In order to keep the technicalities to a minimum we shall below proceed on the assumption that . In the Discussion we give a heuristic description of how the extension to general and can be made. When we are dealing with a fully homozygous population we shall refer to its unique allele as and when we consider but two co-circulating alleles we refer to these as and .

We make the standard assumptions that and all other coefficient functions are smooth and that there are no parental effects, so that , which has as immediate consequence that if , , then and , i.e., the genotype to genotype map is locally additive, , and the same holds good for all quantities that smoothly depend on the phenotype.

Remark 2.1

The biological justification for the above assumptions is that the evolutionary changes that we consider are not so much changes in the coding regions of the gene under consideration as in its regulation. Protein coding regions are in general preceded by a large number of relatively short regions where all sorts of regulatory material can dock. Changes in these docking regions lead to changes in the production rate of the gene product. Genes are more or less active in different parts of the body, at different times during development and under different micro-environmental conditions. The allelic type should be seen as a vector of such expression levels. The genotype to phenotype map maps these expression levels to the phenotypic traits under consideration. It is also from this perspective that we should judge the assumption of smallness of mutational steps : the influence of any specific regulatory site among its many colleagues tends to be relatively minor.

The individual-based microscopic model from which we start is a stochastic birth and death process, with density-dependence through additional deaths from ecological competition, and Mendelian reproduction with mutation. We assume that the population’s size scales with a parameter tending to infinity while the effect of the interactions between individuals scales with . This allows taking limits in which we count individuals weighted with . As an interpretation think of individuals that live in an area of size such that the individual effects get diluted with area, e.g. since individuals compete for living space, with each individual taking away only a small fraction of the total space, the probability of finding a usable bit of space is proportional to the relative frequency with which such bits are around.

2.1 Model setup

The allelic trait space is assumed to be a closed and bounded interval of . Hence the phenotypic trait space is compact. For any , we introduce the following demographic parameters, which are all assumed to be smooth functions of the allelic traits and thus bounded. Moreover, these parameters are assumed to depend in principle on the allelic traits through the intermediacy of the phenotypic trait. Since the latter dependency is symmetric, we assume that all coefficient functions defined below are symmetric in the allelic traits.

: the per capita birth rate (fertility) of an individual with genotype .

: the background death rate of an individual with genotype .

: a parameter scaling the per capita impact on resource density and through that the population size.

: the competitive effect felt by an individual with genotype from an individual with genotype . The function is customarily referred to as competition kernel.

: the mutation probability per birth event (assumed to be independent of the genotype). The idea is that is made appropriately small when we let increase.

: a parameter scaling the mutation amplitude.

: the mutation law of a mutant allelic trait from an individual with allelic trait , with a probability measure with support . As a result the support of is of size .

Notational convention: When only two alleles and co-circulate, we will use the shorthand:


To keep things simple we take our model organisms to be hermaphrodites which in their female role give birth at rate and in their male role have probabilities proportional to to act as the father for such a birth.

We consider, at any time , a finite number of individuals, each of them with genotype in . Let us denote by the genotypes of these individuals. The state of the population at time , rescaled by , is described by the finite point measure on

(2.1)

where is the Dirac measure at .

Let denote the integral of the measurable function with respect to the measure and the support of the latter. Then and for any , the positive number is called the density at time of genotype .

Let denote the set of finite nonnegative measures on , equipped with the weak topology, and define

An individual with genotype in the population reproduces with an individual with genotype at a rate .

With probability reproduction follows the Mendelian rules, with a newborn getting a genotype with coordinates that are sampled at random from each parent.

At reproduction mutations occur with probability and then change one of the two allelic traits of the newborn from to with drawn from .

Each individual dies at rate

The competitive effect of individual on an individual is described by an increase of of the latter’s death rate. The parameter scales the strength of competition: the larger , the less individuals interact. This decreased interaction goes hand in hand with a larger population size, in such a way that densities stay well-behaved. Appendix A summarizes the long tradition of and supposed rationale for the representation of competitive interactions by competition kernels.

For measurable functions and , symmetric, let us define the function on by .

For a genotype and a point measure , we define the Mendelian reproduction operator

(2.2)

and for a measure on parametrized by , we define the Mendelian reproduction-cum-mutation operator

(2.3)
(2.4)

The process is a -valued Markov process with infinitesimal generator defined for any bounded measurable functions from to

and by

(2.5)

The first term describes the deaths, the second term describes the births without mutation and the third term describes the births with mutations. (We neglect the occurrence of multiple mutations in one zygote, as those unpleasantly looking terms will become negligible anyway when goes to zero.) The density-dependent non-linearity of the death term models the competition between individuals and makes selection frequency dependent.

Let us denote by (A) the following three assumptions

(A1)

The functions , , and are smooth functions and thus bounded since is compact.Therefore there exist such that

(A2)

for any , and there exists such that .

(A3)

For any , there exists a function , , such that for any and .

For fixed , under (A1) and (A3) and assuming that , the existence and uniqueness in law of a process on with infinitesimal generator can be adapted from the one in Fournier-Méléard Fournier and Méléard (2004) or Champagnat, Ferrière and Méléard (2008). The process can be constructed as solution of a stochastic differential equation driven by point Poisson measures describing each jump event. Assumption (A2) prevents the population from exploding or going extinct too fast.

3 The short term large population and rare mutations limit: how selection changes allele frequencies

In this section we study the large population and rare mutations approximation of the process described above, when tends to infinity and tends to zero. The limit becomes deterministic and continuous and the mutation events disappear.

The proof of the following theorem can be adapted from Fournier and Méléard (2004).

Theorem 3.1

When tends to infinity and if converges in law to a deterministic measure , then the process converges in law to the deterministic continuous measure-valued function solving

Below we have a closer look at the specific cases of genetically mono- and dimorphic initial conditions.

3.1 Monomorphic populations

Let us first study the dynamics of a fully homozygote population with genotype corresponding to a unique allele and genotype . Assume that the initial condition is , with converging to a deterministic number when goes to infinity.

In that case the population process is where is a logistic birth and death process with birth rate and death rate . The process converges in law when tends to infinity to the solution of the logistic equation

(3.1)

with initial condition . This equation has a unique stable equilibrium equal to the carrying capacity:

(3.2)

3.2 Genetic dimorphisms

Let us now assume that there are two alleles and in the population (and no mutation). Then the initial population has the three genotypes , and . We use to denote the respective numbers of individuals with genotype , and at time , and to indicate the typical state of the population. Let

be the relative frequency of in the gametes. Then the population dynamics is a birth and death process with three types and birth rates and death rates defined as follows.

(3.3)
(3.4)

To see this, it suffices to consider the generator (2.5) with ; for instance, .

Proposition 3.2

Assume that the initial condition converges to a deterministic vector when goes to infinity. Then the normalized process converges in law when tends to infinity to the solution of

(3.5)

where

(3.6)

with

and similar expressions for the other terms.

Due to its special functional form, the vector field has some particular properties. We summarize some of them in the following Propositions.

Proposition 3.3

The vector field (3.6) has two fixed points and (denoted below by and ) where

The () Jacobian matrix has the eigenvalues (negative by assumption ), , and

An analogous result holds for .

This result follows from a direct computation left to the reader.


As we will see later on, the eigenvalue will play a key role in the dynamics of trait substitutions. It describes the initial growth rate of the number of individuals in a resident population of individuals and is called the invasion fitness of an mutant in an resident population. It is a function of the allelic traits and .


Notation: When we wish to emphasize the dependence on the two allelic traits , we use the notation

(3.7)

Note that the function is not symmetric in and and that moreover

(3.8)

In Appendices B and C the long term behavior of the flow generated by the vector field (3.6) is analyzed in more detail. The main conclusions are:

Proposition 3.4

First consider the case when the mutant and resident traits are precisely equal. Then the total population density goes to a unique equilibrium and the relative frequencies of the genotypes go to the Hardy-Weinberg proportions , i.e., there exists a globally attracting one-dimensional manifold filled with neutrally stable equilibria parametrized by , with as stable manifolds the populations with the same .

For the mutant and resident sufficiently close, this attracting manifold transforms into an invariant manifold connecting the pure resident and pure mutant equilibria. When the pure resident equilibrium attracts only in the line without any mutant alleles and its local unstable manifold is contained in the aforementioned invariant manifold (Theorem C.1). When moreover the traits are sufficiently far from an evolutionarily singular point (defined by ) the movement on the invariant manifold is from the pure resident to the pure mutant equilibrium, and any movement starting close enough to the invariant manifold will end up in the pure mutant equilibrium (Theorem C.2).

4 The long term large population and rare mutations limit: trait substitution sequences

In this section we generalize the clonal theory of adaptive dynamics to the diploid case. We again make the combined large population and rare mutation assumptions, except that we now change the time scale to stay focused on the effect of the mutations. Recall that the mutation probability for an individual with genotype is . Thus the time scale of the mutations in the population is . We study the long time behavior of the population process in this time scale and prove that it converges to a pure jump process stepping from one homozygote type to another. This process will be a generalization of the simple Trait Substitution Sequences (TSS) that for the haploid case were heuristically derived in Dieckmann and Law (1996), and Metz et al. (1996) where they were called ’Adaptive Dynamics’, and rigorously underpinned in Champagnat (2006), Champagnat and Méléard (2011).

Let us define the set of measures with single homozygote support.

We will denote by the subset of where vanishes. We make the following hypothesis.

Hypothesis 4.1

For any we have

This hypothesis implies that the zeros of are isolated (see Dieudonné (1969)), and since is closed and compact, is finite.

Definition 4.2

The points such that are called evolutionary singular strategies (ess).

Note that because of (3.8),

Let us now define the TSS process which will appear in our asymptotic.

Definition 4.3

For any , we define the pure jump process with values in , as follows: its initial condition is and the process jumps from to with rate

(4.1)
Remark 4.4

Under our assumptions, the jump process is well defined on . Note moreover that the jump from to only happens if the invasion fitness .

We can now state our main theorem.

Theorem 4.5

Assume (A). Assume that with converging in law to uniformly bounded in and such that . (That is, the initial population is monomorphic for a type that is not an ess). Assume finally that

(4.2)

For introduce the stopping time

(4.3)

where is the distance on the allelic trait space.

Extend with the cemetery point .

Then there exists such that for all , the process converges (in the sense of finite dimensional distributions on equipped with the topology of the total variation norm) to the -valued Markov pure jump process with

where

The process is defined as follows: and jumps

with and infinitesimal rate (4.1).

Remark 4.6

Close to singular strategies the convergence to the TSS slows down. To arrive at a convergence proof it is therefore necessary to excise those close neighborhoods. This is done by means of the stopping times and : we only consider the process for as long as it stays sufficiently far away from any singular strategies. Assumptions (A) imply that the thus stopped TSS is well defined on . Since its jump measure is absolutely continuous with respect to the Lebesgue measure, it follows that converges almost surely to when tends to (for any fixed ).

We now roughly describe the successive steps of the mutation, invasion and substitution dynamics making up the jump events of the limit process, following the biological heuristics of Dieckmann and Law (1996); Metz et al. (1996); Metz (in press). The details of the proof are described in Appendix D, based on the technical Appendices B and C.

The time scale separation that underlies the limit in Theorem 4.5 both simplifies the processes of invasion and of the substitution of a new successful mutant on the population dynamical time scale and compresses it to a point event on the evolutionary time scale. The two main simplifications of the processes of mutant invasion and substitution are the stabilization of the resident population before the occurrence of a mutation, simplifying the invasion dynamics, and the restriction of the substitution dynamics to a competition between two alleles. In the jumps on the evolutionary time scale these steps occur in opposite order. First comes the attempt at invasion by a mutant, then, if successful, followed by its substitution, that is, the stabilization to a new monomorphic resident population. After this comes again a waiting time till the next jump.

To capture the stabilization of the resident population, we prove, on the assumption that the starting population is monomorphic with genotype , that for arbitrary fixed for large the population density with high probability stays in the -neighborhood of until the next allelic mutant appears. To this aim, we use large deviation results for the exit problem from a domain (Freidlin and Wentzel (1984)) already proved in Champagnat (2006) to deduce that with high probability the time needed for the population density to leave the -neighborhood of is bigger than for some . Therefore, until this exit time, the rate of mutation from in the population is close to and thus, the first mutation appears before this exit time if one assumes that

Hence, on the time scale the population level mutation rate from parents is close to

To analyze the fate of these mutants , we divide the population dynamics of the mutant alleles into the three phases shown in Fig. 4.1, in a similar way as was done in Champagnat (2006).

Figure 4.1: Simulation of the three phases of mutant invasion.

In the first phase (between time 0 and in Fig. 4.1), the number of mutant individuals of genotype or is small, and the resident population with genotype stays close to its equilibrium density . Therefore, the dynamics of the mutant individuals with genotypes and is close to a bi-type birth and death process with birth rates and and death rates and for a state . If the fitness is positive (i.e. the branching process is super-critical), the probability that the mutant population with genotype or reaches at some time is close to the probability that the branching process reaches , which is itself close to its survival probability when is large.

Assuming the mutant population with genotype or reaches , a second phase starts. When , the population densities are close to the solution of the dynamical system (3.5) with the same initial condition, on any time interval . The study of this dynamical system (see Appendices B and C) implies that, if the mutation step is sufficiently small, then any solution to the dynamical system starting in some neighborhood of converges to the new equilibrium as time goes to infinity. Therefore, with high probability the population densities reach the -neighborhood of at some time . Applying the results in Theorems C.1 and C.2 for the deterministic system to the approximated stochastic process, is justified by observing that the definition of the stopping times and implies that the allelic trait stays at all times away from the set .

Finally, in the last phase, we use the same idea as in the first phase: since is a strongly locally stable equilibrium, we can approximate the densities of the traits and by a bi-type sub-critical branching process. Therefore, they reach in finite time and the process comes back to where we started our argument (a monomorphic population), until the next mutation.

In Champagnat and Méléard (2011) it is proved that the duration of these three phases is of order . Therefore, under the assumption

the next mutation occurs after these three phases with high probability. Then the time scale Assumption (4.2) allows us to conclude, taking the limits tending to infinity and then to . Then we repeat the argument using the Markov property.


Note that the convergence cannot hold for the usual Skorohod topology and the space equipped with the corresponding weak topology. Indeed, it can be checked that the total mass of the limit process is not continuous, which would be in contradiction with the -tightness of the sequence , which would hold in case of convergence in law for the Skorohod topology (since the jump amplitudes are equal to and thus tend to as tends to infinity).

However, certain functionals of the process converge in a stronger sense. Let us for example consider the average over the population of the phenotypic trait . This can be easily extended to more general symmetric functions of the allele.

Theorem 4.7

Assume that is strictly monotone. Define

where .

Under the assumptions of Theorem 4.5, the process

converges in law in the sense of the Skorohod topology to the process where .

The Skorohod topology is a weaker topology than the usual topology, allowing processes with jumps tending to to converge to processes with jumps (see Skorohod (1956)). For a càd-làg function <