Adaptive dynamics in logistic branching populations
aa r X i v : . [ m a t h . P R ] J un Adaptive dynamics in logistic branchingpopulations
Nicolas Champagnat ∗ , Amaury Lambert † October 25, 2018
MSC 2000 subject.
Primary 60J80; secondary 60J25, 60J70, 60J75,60J85, 60K35, 92D10, 92D15, 92D25, 92D40.
The recent biological theory of adaptive dynamics [6,9] proposes a descrip-tion of the long term evolution of an asexual population by putting emphasison the ecological interactions between individuals, in contrast with classi-cal population genetics models which focus on the genetic structure of thepopulation. The basic models are individual-based models in which thepopulation dynamics is precisely described and includes birth, death, com-petition and mutation. The basic idea of the theory of adaptive dynamicsis to try to get insights about the interplay between ecology and evolutionby studying the invasion of a single mutant type appearing in a simplified(monotype stable) resident population. The evolution of the population canthen be described as a sequence of mutant invasions in the population. Ifthe resident type goes extinct when the mutant type invades (we say thatthe mutant type fixates ), the evolution is described by the so-called ‘traitsubstitution sequence’ (TSS) [10]. This appraoch has revealed powerful topredict the qualitative behaviour of complicated evolutionary dynamics. Inparticular, it allows to determine the (local) direction of evolution in thespace of phenotypic traits (or simply traits ) from the individual ecological ∗ Project-team TOSCA, INRIA Sophia Antipolise-mail: [email protected] † Unit of Mathematical Evolutionary Biology, Laboratoire d’´Ecologie et ´Evolution, Uni-versit´e Pierre et Marie Curie-Paris 6 and ´Ecole Normale Sup´erieure, Parise-mail: [email protected] T n ) n and ( V n ) n , where T n is the n -th time where the population becomesmonomorphic (i.e. composed of only one type) and V n is the surviving typeat time T n . The sequence ( V n ) is the above-mentionned TSS. It is possibleto prove the convergence of an individual-based model to the TSS undertwo biologically motivated assumptions [10,1]. First, the assumption of raremutations guarantees that, in the timescale of mutations, the widths oftime intervals during which the population is polymorphic vanish, so thatthere is only one type surviving at any time t . To prevent the populationfrom rapidly becoming extinct in the new timescale, one also has to rescalepopulation sizes, thereby making the assumption of large populations .Subsequently, the TSS is a Markov jump process on the trait space whosesemigroup is shown [1] to depend on the invasion fitnesses (as defined in [9]) f ( x, y ), x, y ∈ X , where f ( x, y ) is defined as the expected growth rate ofa single individual of type y — the mutant — entering a monomorphicpopulation of type x ‘at equilibrium’ — the residents. Note that this fitnessis not given a priori , but derived from the microscopic model of individualinteractions. Because of the assumption of large population, the sign of thisfitness determines the possibility of invasion of a mutant type: if f ( x, y ) < y cannot invade a resident population of type x . Thus,evolution proceeds by successive invasions of (only) advantageous mutanttypes replacing the resident one.The ‘canonical equation of adaptive dynamics’ [3], which describes theevolution of a one-dimensional trait x as the solution of the following ODE,is obtained from the TSS in the limit of small mutations : dxdt = 12 σ ( x ) µ ( x )¯ n ( x ) ∂∂y f ( x, x ) , (1)where σ ( x ) stands for the (rescaled) variance of the mutation step law,2 n ( x ) for the equilibrium size of a pure x -type population, and f ( x, y ) for theinvasion fitness mentioned above. Note how only advantageous types getfixed (the trait follows the fitness gradient) and how the fitness landscape y f ( x, y ) depends on the current state x of the population.However, it is well-known that slightly deleterious types can be fixed bychance in finite populations. This phenomenon is known under the name of genetic drift . Depending on the strength of genetic drift, selection is said tobe strong (genetic drift has negligible effects) or weak . In the large popula-tion asymptotic from which the TSS of adaptive dynamics is derived, geneticdrift has negligible impact compared to the action of selection. Therefore,the fixation of slightly deleterious types cannot be observed. Our goal hereis to include genetic drift in the adaptive dynamics models by considering finite populations under weak selection . We continue using the bottom-upapproach of adaptive dynamics; that is, model (macroscopic) evolution from(microscopic) populations. In particular, we allow the population sizes tofluctuate randomly through time and we aim to reconstruct a fitness functionfrom the microscopic parameters.After the description of the model (Section 2), we derive a new TSS inthe limit of rare mutations (Section 3), from which a limit of small mu-tations gives what we call the ‘canonical diffusion of adaptive dynamics’(Section 4). The coefficients of this diffusion involve the first-order deriva-tives of the fixation probabilities, which are computed in Section 5 as a linearcombination of four fundamental components associated to fertility , defence , aggressiveness and isolation . New numerical results on the robustness of thepopulation with respect to these fundamental components are also given, aswell as some consequences on the canonical diffusion of adaptive dynamicsin large populations. We will restrict here to logistic interaction. More general models are con-sidered in [2].A monotype (binary) logistic branching process (LBP, see [7]) with dy-namical parameters ( b, c, d ) is a Markov chain in continuous time ( X t ; t ≥ q ij = bi if j = i + 1 ci ( i −
1) if j = i − − i ( b + c ( i − j = i ci ( i −
1) describes competition mortality due to randomencounters between individuals. Other terms correspond to independentbirth events with constant individual rates. This Markov chain is positive-recurrent and converges in distribution to a r.v. ξ , where ξ is a Poissonvariable of parameter θ := b/c conditioned on being nonzero P ( ξ = i ) = e − θ − e − θ θ i i ! i ≥ . (2)We consider in the sequel, a multitype asexual birth and death pro-cess with mutation, generalizing this LBP. At any time t , the population iscomposed of a finite number N ( t ) of individuals characterized by their phe-notypic traits x ( t ) , . . . , x N ( t ) ( t ) belonging to a given trait space X , assumedto be a closed subset of R k . The population state at time t is representedby the counting measure on X ν t = N ( t ) X i =1 δ x i ( t ) . The population dynamics is governed by the following parameters. • b ( x ) is the rate of birth from an individual of type x . The function b is assumed to be C b . • c ( x, y ) is the rate of death of an individual of type x due to the com-petition with another individual of type y . Therefore, the total deathrate of an individual of type x in a population ν may be written as R c ( x, y )( ν ( dy ) − δ x ( dy )). In this expression, the Dirac mass at x sub-stracted to the measure ν means that the individual does not competewith himself. The function c is assumed to be C b and bounded awayfrom 0 on X . • γµ ( x ) is the probability that a birth from an individual with trait x produces a mutant individual, where µ ( x ) ∈ [0 ,
1] and where γ ∈ (0 , γ → • M ( x, dh ) is the law of the trait difference h = y − x between a mutantindividual with trait y born from an individual with trait x . We assumethat M ( x, dh ) has 0 expectation (no mutation bias), i.e. R hM ( x, dh ) =0, and has a density on R k which is uniformly bounded in x ∈ X bysome function ¯ M ( h ) with finite third-order moment.4e will denote the dependence of ν t on the parameter γ with the notation ν γt . Observe that such a population cannot go extinct because the deathrate is 0 when there is only one individual in the population. Since we wantto apply a limit of rare mutations while keeping the population size finite,this is necessary to prevent the population to become extinct before anymutation occur.Let ξ ( x ) be a random variable whose law is the stationary distributionof a pure x -type population with no mutation ( µ ≡ θ is replaced by θ ( x ) := b ( x ) /c ( x, x ).The last notation needed concerns a population with initially only twotypes x and y and with no mutation. Then ν t = X t δ x + Y t δ y , where ( X t , Y t : t ≥
0) is a bivariate Markov chain. For this Markov chain, P ( T < ∞ ) = 1,where T is the first time where either X t or Y t reach 0. We call fixation (of the mutant y ) the event { X T = 0 } . The probability of fixation will bedenoted by u n,m ( x, y ) u n,m ( x, y ) := P ( X T = 0 | X = n, Y = m ) . In this section, we apply the limit of rare mutations ( γ →
0) to the process ν γ , in order to describe the evolution of the population as a TSS in finitepopulation. This limit requires to rescale time properly, as t/γ , to describethe evolution on the mutation timescale. Theorem 3.1
Fix x ∈ X . Assume that ν γ = N γ δ x where sup γ ∈ (0 , E (( N γ ) p ) < ∞ for some p > . Then, for any < t < . . . < t n , the n -tuple ( ν γt /γ , . . . , ν γt n /γ ) converges in law for the weak topology to ( N t δ S t , . . . , N t n δ S tn ) where (1) ( S t ; t ≥ is a Markov jump process on X with initial value S = x andwhose jumping rates q ( x, dh ) from x to x + h are given by q ( x, dh ) = β ( x ) χ ( x, x + h ) M ( x, dh ) , where β ( x ) = µ ( x ) b ( x ) E ( ξ ( x )) = µ ( x ) b ( x ) θ ( x ) / (1 − e − θ ( x ) ) and χ ( x, y ) = X n ≥ n P ( ξ ( x ) = n ) E ( ξ ( x )) u n, ( x, y ) = X n ≥ e − θ ( x ) θ ( x ) n − ( n − u n, ( x, y ) . (3)5 Conditional on ( S t , . . . , S t n ) = ( x , . . . , x n ) , the N t i are independentand respectively distributed as ξ ( x i ) . Therefore, in the limit of rare mutations, on the mutation timescale, thepopulation is always monomorphic and the dominant trait of the populationevolves as a jump process over the trait space, where a jump corresponds tothe appearance and fixation of a mutant type. Moreover, at any time, thepopulation size stationary (i.e. has the stationary distribution correspondingto the dominant trait of the population). The fixation rate of a mutant isgoverned by the function χ ( x, y ), which is therefore the random analogue ofthe traditional invasion fitness [9], defined as the probability of invasion ofa mutant type y is a resident population of type x at equilibrium. Observethat, as usual in adaptive dynamics, the fitness landscape depends on thecurrent state of the population. Moreover, in contrast with the classicalTSS [10], from a given monomorphic resident population, any mutant traithas a positive probability to invade (by genetic drift). Therefore, evolution ispossible in any direction of the trait space. However, a directional selectionstill exists, as will appear in the next section.We refer to [2] for the proof of this result. However, this convergence isnatural in view of the following interpretation of each parameter. β ( x ) canbe seen as the mean mutant production rate of a stationary x -type population(i.e. with size ξ ( x )), and χ ( x, y ) is the probability of fixation of a single y -type mutant entering a pure x -type population with size-biased stationarysize . The size bias comes from the fact that the mutant appears at a birthtime in the stationary population (since the birth rate is proportional tothe population size, the population size after a birth event in the stationarypopulation is given by the size-biaised stationary population size). Let us assume for simplicity that X = R k . Let σ ( x ) be the square rootmatrix of the covariance matrix of M ( x, · ). We also need to assume that thematrix σ ( x ) is a Lipschitz function of x .In order to obtain the equivalent of the canonical equation of adaptivedynamics in a finite population, we want to apply a limit of small mutationsteps (weak selection) to the TSS S . To this aim, we introduce a param-eter ǫ > M ( x, · ) by their image by theapplication h ǫh . Time also has to be rescaled in order to obtain anon-degenerate limit. The correct time scaling is 1 /ǫ , which leads to the6ollowing generator for the rescaled TSS ( S ǫt ; t ≥ A ǫ ϕ ( x ) = 1 ǫ Z R k ( ϕ ( x + ǫh ) − ϕ ( x )) β ( x ) χ ( x, x + ǫh ) M ( x, dh ) . Using the assumption that the mutation kernels M ( x, · ) have 0 expectation,it is elementary to compute the limit of this expression as ǫ → ϕ ). This limit, which takes the form of a diffusion generator,explains the following result (its full proof can be found in [2]). Theorem 4.1
If the family ( S ǫ ) ǫ> has bounded first-order moments andconverges in law as ǫ → to a random variable Z , then the process S ǫ withinitial state S ǫ converges in law for the Skorohod topology on D ( R + , R k ) tothe diffusion process ( Z t ; t ≥ with initial state Z unique solution to thestochastic differential equation dZ t = β ( Z t ) σ ( Z t ) · ∇ χ ( Z t , Z t ) dt + p β ( Z t ) χ ( Z t , Z t ) σ ( Z t ) · dB t (4) where ∇ χ denotes the gradient w.r.t. the second variable y of χ ( x, y ) and B is a standard k -dimensional Brownian motion. This result gives the equivalent of the canonical equation of adaptivedynamics (1) when the population is finite. It is no longer a deterministicODE, but a diffusion process, in which the genetic drift remains present (inthe form of a stochastic diffusion term), as a consequence of the populationfiniteness and of the limit of weak selection. The deterministic drift partof (4) is very similar to the standard canonical equation of adaptive dynam-ics (1), and involves in particular the gradient of the fitness function χ . Theprocess (4) provides a diffusion model describing the evolution of the dom-inant trait value in a population [8,5], grounded on a precise microscopicdensity-dependent modelling of the population dynamics. It also gives theprecise balance between directional selection and genetic drift as a functionof the individual’s dynamical parameters. The SDE (4) involves the fixation probability χ ( x, x ) and the fitness gradientwith respect to the second variable ∇ χ ( x, x ). In this section, we explainhow these quantities can be explicitly computed.We need to compute the derivatives of the fixation probabilities u n,m ( x, y )when y is close to x . Recall that the law of the two-types LBP without mu-tation ( X, Y ) used to define u n,m in the end of Section 2 is characterized by7he birth vector B and the competition matrix CB = (cid:18) b ( x ) b ( y ) (cid:19) , C = (cid:18) c ( x, x ) c ( x, y ) c ( y, x ) c ( y, y ) (cid:19) . We will say that the mutant is neutral if all individuals are exchangeable,i.e. when b ( y ) = b ( x ) and c ( x, y ) = c ( y, x ) = c ( y, y ) = c ( x, x ) (this holds inparticular when y = x ). As will appear below, using the notation ( b, c ) :=( b ( x ) , c ( x, x )), it is natural to focus on deviations from the neutral caseexpressed as B = b + (cid:18) λ (cid:19) , C = c − (cid:18) δ δ (cid:19) + (cid:18) α α (cid:19) − (cid:18) εε (cid:19) In words, deviations from the neutral case are a linear combination of four fundamental selection coefficients λ , δ , α , ε , that are chosen to be positivewhen they confer an advantage to the mutant. It is convenient to assessthese deviations to the neutral case in terms of1. fertility ( λ ): positive λ means increased mutant birth rate2. defence capacity ( δ ): positive δ means reduced competition sensitiv-ity of mutant individuals w.r.t. the total population size3. aggressiveness ( α ): positive α means raised competition pressureexerted from any mutant individual onto the rest of the population4. isolation ( ε ): positive ε means lighter cross-competition between thetwo different types, that would lead, if harsher, to a greater probabilityof exclusion of the less abundant oneUnder neutrality, an elementary martingale argument shows that the fixa-tion probability equals the initial mutant frequency p := m/ ( m + n ). Thisimplies in particular that χ ( x, x ) = e − θ ( x ) − θ ( x ) θ ( x ) . (5)The following theorem unveils the dependence of u upon λ , δ , α , ε , whenthey slightly deviate from 0, and explains why these four selection coefficientsprovide a natural basis to decompose the gradient of the fixation probability.8 heorem 5.1 As a function of the multidimensional selection coefficient s = ( λ, δ, α, ε ) , the probability u is differentiable, and in a neighbourhood of s = 0 (selective neutrality), u = p + v ′ .s + o ( s ) , (6) where the selection gradient v = ( v λ , v δ , v α , v ε ) can be expressed as v ιn,m = p (1 − p ) g ιn + m ι = ε,v εn,m = p (1 − p ) (1 − p ) g εn + m where the g ’s depend solely on the resident’s characteristics b, c , and on thetotal initial population size n + m . They are called the invasibility coefficients. The invasibility coefficients of a pure resident population are interestingto study, as they provide insights about how the fixation probability devi-ates from p as the selection coefficients of the mutant deviate from 0. Theyalso provide information about the robustness of the resident population,i.e. its resistance to mutant invasions. In particular, this allows to comparethe sensitivity of the invasion probability in a given monomorphic residentpopulation with respect to the four fundamental selection coefficients. Inthe simplest case where mutations in the parameter space are isotropic,the biggest invasibility coefficent gives the direction of the parameter spacewhere a mutant is more likely to invade. More generally, when there arecorrelations between mutations in the parameter space (either because ofthe phenotypic structure in the functions b ( · ) and c ( · , · ), or because the co-variance matrix of the mutation steps is non-diagonal), the likeliest directionof evolution in the trait space is given by the deterministic coefficient of thecanonical diffusion (4), in which the fitness gradient is given by ∇ χ ( x, x ) = a λ ( x ) ∇ b ( x ) − a δ ( x ) ∇ c ( x, x ) + a α ( x ) ∇ c ( x, x ) , (7)where, for ι = λ, δ, α , a ι ( x ) = e − θ ( x ) ∞ X n =1 ng ιn +1 ( x ) θ ( x ) n − ( n + 1) ( n − . It is possible to obtain explicit expressions for the invasibility coefficients g ι as series. We refer to [2] for the exact expressions. In particular, theseexpressions yield that a ι ( x ) = ˆ a ι ( θ ( x )) /c ( x, x ) for some function ˆ a ι . More-over, they allow one to compute numerically the invasibility coefficients, andtherefore the quantities ˆ a ι for ι = λ, δ, α as functions of the parameter θ ( x ).90.10.20.30.40.50.6 0 5 10 15 20 25 30 35 40ˆ a α ( θ )ˆ a λ ( θ )ˆ a δ ( θ )Figure 1: The functions ˆ a λ , ˆ a δ and ˆ a α as functions of θ .00.10.20.30.40.50.6 0 5 10 15 20 25 30 35 40 θ ˆ a α ( θ ) θ ˆ a λ ( θ )Figure 2: The functions θ θ ˆ a λ ( θ ) and θ θ ˆ a α ( θ ).10hese numerical results can be used to make simulations of the canonicaldiffusion of adaptive dynamics in various ecological examples. In particular,in contrast with the classical canonical equation of adaptive dynamics, thepresence of a genetic drift can induce the evolutionary dynamics to driftaway from evolutionary stable strategies, where the fitness gradient is zero.When the fitness gradient admits several zeros, it can visit several basins ofattraction on various timescales.This numerical study is a work in progress, that is quite delicate becausethe series involved in the computation of g ι are slowly converging, with firstterms that grow exponentially fast with θ . We shall give here our firstresults. Fig. 1 shows the functions ˆ a ι for ι = λ, δ, α . Several comments canbe made from this figure. First, for any θ >
0, ˆ a δ ( θ ) > ˆ a α ( θ ) > ˆ a λ ( θ ). Thismeans that, for equal mutation steps in the parameter space, a mutation isalways more advantageous in the direction δ than in the direction α , whichis itself more advantageous than in the direction λ . In other words, in agiven population, a better defence capacity is more beneficial than a betteraggressiveness, which is more beneficial than a better fertility.Moreover, as θ goes to infinity, these functions have different asymptoticbehaviors. ˆ a δ seems to converges to 1 /
2, whereas ˆ a λ ( θ ) and ˆ a α ( θ ) are bothequivalent to 1 / θ (see Fig. 2). However, this does not mean that mutationsare much more likely to fixate in the δ direction, because, when θ is large, b = θc is larger than c , and, in (7), a δ and a α are multiplied by ∇ c and ∇ c respectively, whereas a λ is multiplied by ∇ b . More formally, to computethe limit of the canonical diffusion when θ goes to infinity, one can dividethe competition kernel c ( · , · ) by a constant K in the microscopic model andthen let K go to infinity in the canonical diffusion. Denoting by χ K ( · , · ) thefitness function obtained this way and using the asymptotic behaviors givenabove, one gets thatlim K → + ∞ ∇ χ K ( x, x ) = 12 b ( x ) ( ∇ b ( x ) − θ ( x ) ∇ c ( x, x )) . Moreover, by (5), χ K ( x, x ) converges to 0 when K → ∞ . Now, as provedin [1], with our notation, the fitness function of the canonical equation ofadaptive dynamics (1) is given by f ( x, y ) = b ( y ) − c ( y, x ) θ ( x ). Therefore,as K → + ∞ , the canonical diffusion converges to a deterministic ODEwhich is precisely the canonical equation of adaptive dynamics. This givesa new justification of this equation, and this also allows one to study thefluctuations around the canonical diffusion when K is large. In particular,this diffusion equation with small diffusion term enters the framework ofFreidlin-Wentzell’s theory [4], which can be applied to predict the long time11ehaviour of the diffusion, and its chain of visit of basins of attractions when K is large. This kind of information is biologically very relevant, since itallows one to predict in which order all the evolutionary stable strategieswill be visited by the population and on which timescale. References [1] N. Champagnat,
A microscopic interpretation for adaptive dynamicstrait substitution sequence models , Stoch. Proc. Appl. 116 (2006), 1127–1160.[2] N. Champagnat and A. Lambert,
Evolution of discrete populations andthe canonical diffusion of adaptive dynamics , to appear in Ann. Appl.Prob.[3] U. Dieckmann and R. Law,
The dynamical theory of coevolution: Aderivation from stochastic ecological processes , J. Math. Biol. 34 (1996),579–612.[4] M.I. Freidlin and A.D. Wentzell, Random Perturbations of DynamicalSystems. Springer-Verlag, Berlin, 1984.[5] T.F. Hansen,
Stabilizing selection and the comparative analysis of adap-tation , Evolution 51 (1997), 1341–1351.[6] J. Hofbauer and R. Sigmund,
Adaptive dynamics and evolutionary sta-bility , Appl. Math. Letters 3 (1990), 75–79.[7] A. Lambert,
The branching process with logistic growth , Ann. Appl.Prob. 15 (2005), 1506–1535.[8] R. Lande,
Natural selection and random genetic drift in phenotypicevolution , Evolution 30 (1976), 314–334.[9] J.A.J. Metz, R.M. Nisbet and S.A.H. Geritz,
How should we define‘fitness’ for general ecological scenarios ? , Trends in Ecology and Evo-lution 7 (1992), 198–202.[10] J.A.J. Metz, S.A.H. Geritz, G. Mesz´ena, F.A.J. Jacobs and J.S. vanHeerwaarden,