A birth-death model of ageing: from individual-based dynamics to evolutive differential inclusions
Sylvie Méléard Michael Rera Tristan Roget

High mobility InAs 2DEGs on GaSb substrates: a platform for mesoscopic quantum transport
    M. J. Manfra
    December 30, 2023
========================================================================================


Ageing’s sensitivity to natural selection has long been discussed because of its apparent negative effect on an individual's fitness. Thanks to the recently described (Smurf) 2-phase model of ageing (<cit.>) we propose a fresh angle for modeling the evolution of ageing. Indeed, by coupling a dramatic loss of fertility with a high-risk of impending death - amongst other multiple so-called hallmarks of ageing - the Smurf phenotype allowed us to consider ageing as a couple of sharp transitions.
The birth-death model (later called bd-model) we describe here is a simple life-history trait model where each asexual and haploid individual is described by its fertility period x_b and survival period x_d. We show that, thanks to the Lansing effect, the effect through which the “progeny of old parents do not live as long as those of young parents”, x_b  and x_d converge during evolution to configurations x_b-x_d≈ 0 in finite time. 

To do so, we built an individual-based stochastic model which describes the age and trait distribution dynamics of such a finite population. Then we rigorously derive the adaptive dynamics models, which describe the trait dynamics at the evolutionary time-scale. We extend the Trait Substitution Sequence with age structure to take into account the Lansing effect. Finally, we study the limiting behaviour of this jump process when mutations are small. We show that the limiting behaviour is described by a differential inclusion whose solutions x(t)=(x_b(t),x_d(t)) reach the diagonal { x_b=x_d} in finite time and then remain on it. This differential inclusion is a natural way to extend the canonical equation of adaptive dynamics in order to take into account the lack of regularity of the invasion fitness function on the diagonal { x_b=x_d}. 




§ INTRODUCTION

Ageing is commonly defined as an age-dependant increase of the probability to die after the maturation phase (<cit.>). It affects a broad range of organisms in various ways ranging from negligible senescence to fast post-reproductive death (reviewed in <cit.>). In the recent years, a new 2-phases model of ageing proposed by <cit.> described the ageing process not as being continuous but as made of at least 2 consecutive phases separated by a dramatic transition. This transition, dubbed “Smurf transition”, was first described in drosophila (<cit.>, <cit.>). In short, this transition occurs in every individuals prior to death and is marked by a series of associated phenotypes encompassing high-risk of impending death, increased intestinal permeability, loss of energy stores, reduced fertility (<cit.>). It was later showed to be evolutionarily conserved in Caenorhabditis elegans and Danio rerio (<cit.>, <cit.>). Such broad evolutionary conservation of a marker for physiological age raises the question of an active selection of the underlying mechanisms throughout evolution.
Since the beginning of ageing studies, the question of its ability to appear through evolution has been raised. In fact, since the Darwinian theory of evolution stipulates that species arise and develop thanks to the natural selection of small, inherited variations that increase an individual's ability to compete, survive, and reproduce (<cit.>), numerous researchers suggested that ageing - and more precisely senescence - could not be actively and directly selected thanks to evolution (<cit.>). One of the first to publicly address the question of the evolution of ageing was August Weismann who proposed in 1881 that the life expectancy was programmed by “the needs of the species” (<cit.>). Numerous theoretical works have been developed about ageing for the past 60 years in order to recenter the selection of an ageing process on the individuals more than on the population. 
Here we will focus our attention on the capability of a process such as ageing to be selected through evolution.

If fitness alone - as an individual's reproductive success or its average contribution to the gene pool of the next generation - were at play in the evolution process, the best adapted individuals would have infinite fertility as well as longevity. Nevertheless, this situation is never observed mainly because organisms adapted to constant variations of environmental conditions and physical limitations of resources availability. Thus, an active mechanism for the elimination of these fitness-excessive individuals would represent a selective advantage in an environment where scarcity is the rule. The Lansing effect is a good candidate for such a mechanism. It is the effect through which the “progeny of old parents do not live as long as those of young parents” first described in rotifers (<cit.>, <cit.>). More recently, it has been shown that older drosophila females and in some extent males tend to produce shorter lived offspring (<cit.>), zebra finch males give birth to offspring with shorter telomere lengths and reduced lifespans (<cit.>) and finally in humans, “Older father’s children have lower evolutionary fitness across four centuries and in four populations” (<cit.>).

In the present article, we decided to approach the problem of ageing selection and evolution by using an extremely simplified version of a living organism. It is an haploid and asexual organism carrying only two traits, x_b that defines the duration of its ability to reproduce and x_d that defines the duration of its ability to maintain its integrity - stay alive (see Figure <ref>). We will further discuss the properties of this simple model in the next part. Although quite simple, it allows the modeling of all types of observed ageing modes : negligible senescence, sudden post-reproductive death, or post-reproductive “menopause-like” survival as well as the smurf phase.

The main result of the present article is that a pro-senescence program can be selected through Darwinian mechanisms thanks to the Lansing effect. Indeed, our main mathematical result (see Theorem <ref>) shows that evolution drives the trait (x_b,x_d) towards configurations x_b = x_d. It means that the individuals can enjoy all their reproductive capacity, and then are quickly removed from the population. 
Moreover, this result shows that after reaching the configurations x_b=x_d, the traits x_b and x_d continue to increase with decreasing speed,  while maintaining x_b=x_d. This decrease in the speed of evolution is a consequence of the fitness gradients being decreasing functions of the traits (see Remark <ref>). It is related to the well-known fact that the strength of selection decreases with age, i.e that  a mutation having an effect on the reproduction or mortality rates at a given age will have all the more impact as this age is small (<cit.>, <cit.>, <cit.>). Indeed, in our model a perturbation of the trait x_b (resp. x_d) is equivalent to a perturbation of the birth rates at age x_b (resp. death rates at age x_d). 

We build an individual based stochastic model inspired by <cit.>. It describes an asexual and haploid population with a continuous age and a continuous life-history trait structure. In this model, the life-history trait of every individual is thus a pair of positive numbers (x_b, x_d)∈_+^2. An individual with trait (x_b,x_d) reproduces at rate one as long as it is younger than x_b and cannot die as long as it is younger than x_d (see Figure <ref>).

This model leads to three configurations: x_d<x_b, x_d=x_b and x_d>x_b (see Figure <ref>).

From one generation to the next, variation on the trait (x_b,x_d) is generated through genetic mutations. In addition, natural selection occurs through mortality due to competition for resources thanks to a logistic equation defining the maximum load of the medium. Finally, we model an epigenetic effect of senescence through the Lansing Effect. It introduces a source of phenotypic variation at a much faster time-scale than genetic mutations. In that aim, we assume that an individual that reproduces after age x_d transmits to his descendant a shorter life expectancy (see Section 2 for details). Therefore, only individuals with trait x_d < x_b are affected (see Figure <ref> (c)). That creates an adaptive trade-off which impacts the phenotypic evolution of the population. 

The purpose of the present article is to study the long-term evolution of the trait x=(x_b,x_d) and to determine whether it concentrates on x_b -x_d = 0. To do so, we are inspired by the theory of adaptive dynamics (<cit.>, <cit.>, <cit.>) which studies the phenotypic eco-evolution of large populations under the assumption that genetic mutations are rare and have small effects. A central tool in that theory is the concept of invasion fitness. The invasion fitness is a function 1-z(y,x) informally defined as the probability that an individual with trait y survives in a resident population with trait x at demographic  equilibrium. In section 4, we prove that the invasion fitness satisfies the simple relation 1-z(y,x)=(λ(y) -λ(x))∨ 0 where λ(x) is the Malthusian parameter, describing the adaptive value associated with the trait x (see Section 3.1 (<ref>) for the definition). This allows us to introduce the Trait Substitution Sequence process (TSS) which is a pure jump process describing the successive invasions of successfull mutants in monomorphic populations at the demographic equilibrium. The TSS has been heuristically introduced in <cit.>, <cit.> for a trait structured population. In <cit.>, it has been rigorously derived from an individual based model and generalised in <cit.>, to age-structured populations. Our case differs from <cit.> by mainly two aspects: the additional Lansing effect and the specific form of the mutations kernel which is not absolutely continuous with respect to the Lebesgue measure on ^2 (as assumed in <cit.>,  <cit.>). 
In the usual case, the TSS is approximated by the solution of the Canonical equation of adaptive dynamics when the size of mutations is small and on a longer time-scale (<cit.>, <cit.>, <cit.>). This limit theorem requires the Lipschitz regularity of the fitness gradient. In our model this assumption is not satisfied. Nonetheless, we prove that the limiting behaviour of the TSS when mutation are small is captured by a differential inclusion, using an approach  developed in <cit.>. A differential inclusion is an extension of ordinary differential equation to set-valued time-derivatives, which extends Cauchy-Lipschitz theory to non regular gradient cases. In our model, the gradient is smooth except on the diagonal { x_b=x_d}. We prove that the solutions are well-defined until they reach the diagonal which they do in finite time. Indeed, the drift towards the diagonal is due to on one side (x_b<x_d) to the fact that the individuals with larger x_b will reproduce more and thus tend to invade ; and on the other side (x_d<x_b) to the Lansing effect because old individuals produce short-lived offsprings. Hence, there is no advantage associated with an increasing x_b, only when associated with an increasing x_d that maximizes survival. Then the trait (x_b,x_d) evolves following the diagonal. The drift of the differential inclusion depends on the derivatives of the Malthusian parameter with respect to the trait variable. These derivatives are expressed as functions of the stable age distribution and reproductive value as in <cit.>, <cit.> (see Remark <ref>).



In Section 2, we present the individual-based  model. Thanks to simulations, we show what was suggested by observations: the trait distribution of the population stabilises on the diagonal x_b-x_d=0.

 In Section 3, we study the deterministic approximations of the stochastic dynamics under the assumption of large population and rare mutations. These approximations are non-linear systems of partial differential equations similar to the Gurtin-McCamy Equation (<cit.>). We study their long-time behaviour and give some results of convergence to the stationary states.

In Section 4,  we state and prove the main mathematical results of this article concerning the approximation by the TSS and the canonical inclusion of adaptive dynamics. (Theorem <ref> and Theorem <ref>).

 In Section 5, we give some comments on our results.


§ A STOCHASTIC MODEL FOR THE EVOLUTION OF LIFE-HISTORY TRAITS

At time t≥ 0, the population is described by a point measure on (_+^*)^2×_+

    Z_t^K(dx,da)=1/K∑_i=1^N_t^Kδ_(x^i(t),a^i(t))(dx,da)

where N_t^K = K⟨ Z_t^K,1⟩ is the total population size, K is the order of magnitude of a population at equilibrium, x^i(t)=(x_b^i(t),x_d^i(t))∈ (_+^*)^2 is the trait of the individual i and a^i(t)∈_+ is its age. The dynamics is defined as a piecewise deterministic Markov process which jumps as follows:


  •An individual (x,a)∈ (ℝ_+^*)^2×_+ reproduces at rate 1_a≤ x_b. The trait of the descendant y=(y_b,y_d) is determined by the following two-steps mechanism (see Figure <ref> below for an illustration):


  -Step 1: If a≤ x_d, the offspring inherits the trait x=(x_b,x_d).

Lansing Effect: if a>x_d, we assume that the offspring carries the trait  (x_b,0).
Let us denote by x̃  the trait  defined as x̃=x if a≤ x_d and  x̃=(x_b,0) if a>x_d. We observe that x_b remains unchanged and that only individuals with configurations x_d<x_b are concerned by the second type (x_b,0) (see Figure <ref> (c)).

  -Step 2: Genetic mutations. A mutation appears instantaneously on each trait x_b and x̃_d independently with probability p_K∈]0,1[. If the trait x_b mutates, the trait y_b of the descendant is y_b=x_b +h_b where h_b∈ is chosen according to the probability measure k(x_b,h_b)dh_b; if the trait x̃_d mutates, the trait y_d of the newborn is y_d=x̃_d+h_d where h_d∈ is chosen according to k(x̃_̃d̃,h_d)dh_d, where the mutational kernel k is defined for all u∈ℝ_+ and v∈ by

    k(u,v)=1_[0∨(u-1),u+1](u + v)e^-v^2/σ^2/∫_1_[0∨(u-1),u+1](u + z)e^-z^2/σ^2dz

where σ^2>0. Note that for u>1, 

    k(u,v)=k(v)= 1_[-1,1]( v)e^-v^2/σ^2/∫_1_[-1,1](z)e^-z^2/σ^2dz.



  •An individual with trait (x,a) has a death rate 1_a>x_d+η N_t^K, with η>0, meaning that  each individual is subjected to the same competition pressure η from any individual in the population and whatever the value of its trait.

  •Between jumps, individuals age at speed one: an individual with age a at time t has an age a+s at time t+s.

Figure <ref> summarizes the trait dynamics described above. 



  
Numerical simulation. 
 The pictures in Figure <ref>  represent a simulation of the trait marginals dynamics  of the process Z_t^K(dx,da). We consider a monomorphic initial population with trait x=(1.2, 2.5) and N_0^K=10000. We consider a competition rate η=0.0005, a probability of mutation p_K=0.05 and a variance of mutations σ=0.05.

At time t=0, the population is monomorphic with trait (x_b,x_d)=(1.2,2.5). We observe that before the trait x_b reaches the value of 2.5, the trait x_d remains constant. When x_b reaches x_d≃ 2.5, the traits x_b and x_d continue to increase by maintaining x_b≃ x_d.


§ MONOMORPHIC AND BIMORPHIC DETERMINISTIC DYNAMICS

In this section, we study a deterministic approximation of the process Z^K when K goes to infinity. We also assume that p_K goes down to zero: almost no mutation occurs on a time interval [0,T]. Nonetheless, some phenotypic variation is created by Lansing Effect. Since our model is density-dependent, the deterministic approximation is a system of classical non-linear partial differential equations similar to the Gurtin MacCamy Equation  (<cit.>). In the monomorphic case, we show that the dynamics converges to the unique non-trivial equilibrium. In the bimorphic case, we show the convergence to a monomorphic equilibrium. We will consider that a monomorphic population with trait x is composed of two subpopulations with traits (x_b,x_d) and (x_b,0).

We define I=([ 1 0; 0 1 ]).



 §.§ Monomorphic dynamics

Let x=(x_b,x_d)∈ (_+^*)^2 be a phenotypic trait. We consider a monomorphic initial population such that the sequence (Z_0^K)_K weakly  converges as K→∞ to δ_x n_x(a)da which describes a monomorphic population with trait x. Then, as in <cit.>, we can prove that the sequence of processes (Z^K)_K converges in probability on any finite time interval to the weak solution (n_x(t,.),t≥ 0)=((n_x^1(t,.),n_x^2(t,.)),t≥ 0)∈ C(_+, L^1(_+)^2) of the following system of partial differential equations

    ∂_t n_x(t,a) +∂_a n_x(t,a) = -( D_x(a) +η‖ n_x(t,.)‖_1 I)n_x(t,a)
    
    n_x(t,0)=∫__+B_x(α)n_x(t,α)dα

with initial condition given by n_x(0,a)=n_x(a), where the densities n^1_x(t,.) and n^2_x(t,.) describe the population distributions with trait (x_b,x_d) and (x_b,0) respectively; 

    ‖ n_x(t,.)‖_1 = ∑_i∈{ 1,2}∫__+| n_x^i(t,α)| dα

is the total population size and

    B_x(a)=([ 1_a≤ x_b∧ x_d             0;  1_x_d<a≤ x_b      1_a≤ x_b ]),  D_x(a)=([ 1_a>x_d       0;       0       1 ])

are the birth and death interactions. We refer to <cit.> for the well-posedness theory of L^1(_+)^2 solutions of Equation (<ref>). 

Let us comment the different terms in Equation (<ref>). The transport term on the left-hand side of the first equation describes the aging of the individuals, the right-hand side describes the death of the individuals. The renewal condition in a=0 describes the births (i.e the individuals with age 0).


We introduce the set of viable traits

    𝒱={ x=(x_b,x_d)∈(_+^*)^2: x_b∧ x_d > 1}.
 
We show that for any trait x∈𝒱, there exists a unique non-trivial and globally stable stationary state of (<ref>). Indeed, the quantity x_b∧ x_d represents the mean number of descendants with trait (x_b,x_d) of an individual with trait (x_b,x_d) if there is no competition.

Assume that x∈𝒱. There exists a unique non-trivial stationary solution n_x∈ L^1(_+)^2 to Equation (<ref>). Moreover, any L^1(_+)^2 non negative solution n_x(t,.) of (<ref>) such that n^1_x(t,.)≠ 0 converges to n_x in L^1(_+)^2 as t→∞. 


In Proposition <ref> below, we will give an explicit expression for the stationary state n_x.

The proof of Proposition <ref> is based on the  study of the associated linear dynamics. We introduce the linear operator

    [ A:D(A)⊂ L^1(_+)^2                 →         L^1(_+)^2;                 u                 ↦     -u' - D_x(a)u ]

where D(A)={ u∈ L^1(_+)^2: u  abs. cont., u'∈ L^1(_+)^2, u(0)=∫__+B_x(α)u(α)dα}. It is well known that A is the infinitesimal generator of a strongly continuous semigroup of linear operators <cit.> which describes the solutions of the linear system of McKendrick Von-Foerster Equation

    ∂_t v_x(t,a) +∂_a v_x(t,a) = -D_x(a)v_x(t,a)
    
    v_x(t,0)=∫__+B_x(α)v_x(t,α)dα.

In <cit.>, a similar linear model is studied. The "entropy method" (<cit.>) allows the authors to prove the convergence of the normalised solutions to some stable distribution in some weighted L^1-space. We need stronger convergence in order to study the long-time behaviour of the masses of the solutions of (<ref>). Since the birth  matrix B_x is not irreducible (it is triangular) and the parameters B_x and D_x are not smooth, we cannot apply Theorems 4.9 and 4.11 in <cit.>. Nonetheless, we easily extend them to our reducible and non-smooth setting. 

We define the Malthusian parameter λ(x) associated with some trait x as the unique solution of the equation

    ∫_0^x_b∧ x_de^-λ(x) ada=1.

Proposition <ref> below justifies this definition and shows that λ(x) is the asymptotic growth rate of the dynamics defined by (<ref>). Its proof is presented in the Appendix.

Let us define for all z∈ℂ the 2× 2 matrix

    F(z) = ∫__+B_x(a)exp( -∫_0^a(D_x(α)+ z I)dα)da.

Note that it is well-defined since B_x has compact support.

Assume that x∈𝒱. Then the linear operator A admits a unique pair of simple principal eigenelements (λ(x),N_x)∈_+^*× D(A) where the stable age distribution N_x satisfies

    N_x^1(a)=e^-(λ(x) a + (a-x_d)∨ 0 ),   N_x^2(a)=[F(λ(x))]_21/1-[F(λ(x))]_22 e^-(1+λ(x))a.

Moreover, for any  non-negative solution v_x(t,a) of (<ref>) in L^1(_+)^2, there exists a positive constant c(v_x(0,.)) such that e^-λ(x) tv_x(t,.)→ c(v_x(0,.))N_x in L^1(_+)^2 as t→ +∞.

Let us now give a lemma which will be used to study the long-time behaviour of the
masses of the solutions of (<ref>) and whose proof is postponed to the Appendix.

Let (m_11,m_12,m_22)∈_+^*×_+×_-^*. Let 𝒟_11(t),𝒟_12(t),𝒟_22(t) be continuous functions from _+ to  which tend to zero as t→∞. Let us denote

    M=([ m_11    0; m_21 m_22 ]),  𝒟(t)=([ 𝒟_11(t)       0; 𝒟_21(t) 𝒟_22(t) ]).

Then any solution (z(t),t≥ 0) of the equation

    dz(t)/dt=(M +𝒟(t))z(t) - η‖ z(t)‖_1 z(t)

started at z(0)∈_+^*×_+ converges to a vector z which satisfies

    z_1= m_11/η1/1+m_21/m_11-m_22,  z_2 = m_21/m_11-m_22z_1.

 
We conclude this section by proving Proposition <ref>.

We prove the first assertion. Let x∈𝒱 and let λ(x) be the principal eigenvalue of A given by Proposition <ref>. Let n_x be the (unique) principal eigenvector of A which satisfies η ‖n_x‖_1 =λ(x). It is obvious that n_x is a non-trivial stationary state of (<ref>). Reciprocally, let n be a stationary state of (<ref>). Then we have necessarily  λ(x) =η ‖n‖_1 and that n is an eigenvector of A associated with the eigenvalue λ(x) that allows us to conclude. We now study the long-time behaviour of the solutions. Let us define 

    v_x(t,a)=exp(η∫_0^t‖ n_x(s,.)‖_1 ds)n_x(t,a).

It is straightforward to prove that v_x is a solution of the linear equation (<ref>). By Proposition <ref> we have e^-λ(x) tv_x(t,.)→ c(v_x(0,.))N_x in L^1(_+)^2 as t→∞. We deduce that for i∈{ 1,2} and denoting ρ_x^i(t)=‖ n_x^i(t,.)‖_1,

    n_x^i(t,.)/ρ_x^i(t)=e^-λ(x) tv_x^i(t,.)/∫__+e^-λ(x) tv_x^i(t,α)dα→N^i_x/∫__+N_x^i(α)dα

in L^1(_+) as t→ +∞. We now study the behaviour of the masses ρ_x(t). By taking the derivative under the integral, we obtain that for i∈{ 1,2}

    dρ_x^i(t)/dt   =∑_j=1^2∫__+([B_x(α)]_ij -[D_x(α)]_ij)n_x^j(t,α)dα - ηρ_x^i(t)‖ρ_x(t)‖_1
       = ∑_j=1^2ρ_x^j(t)∫__+([B_x(α)]_ij -[D_x(α)]_ij)n_x^j(t,α)/ρ_x^j(t)dα -ηρ_x^i(t)‖ρ_x(t)‖_1.

Hence we obtain by (<ref>) that

    dρ_x(t)/dt=(A +𝒜(t))ρ_x(t) - η‖ρ_x(t)‖_1ρ_x(t)

where A=(a_ij) and for (i,j)∈{ 1,2}^2,

    a_ij= ∫__+([B_x(α)]_ij -[D_x(α)]_ij)N^j_x(α)/∫__+N^j_xdα,

and 𝒜(t) is a continuous function decreasing to zero as t tends to infinity. Since a_11=λ(x) >0, a_21≥ 0, a_12=0 and a_22<0, Lemma <ref> allows us to conclude that ρ_x(t) converges to ρ_x, which is defined as the unique  solution of the equation Aρ_x - η‖ρ_x‖_1 ρ_x =0. We easily solve this system and we obtain that ρ_x satisfies 

    ρ^1_x=λ(x)/η1/1 + a_21/λ(x) -a_22,  ρ^2_x = a_21/λ(x) -a_22ρ^1_x.


We conclude this section by writing more explicit formulas for the stationary state n_x.

Let x∈𝒱 then we have

    n^1_x(a)=ρ^1_xN_x^1(a)/∫__+ N_x^1(α)dα,  n^2_x(a)=ρ^2_x N_x^2(a)/∫__+ N_x^2(α)dα

where

    ρ^1_x=λ(x)/η1/1 + a_21/λ(x) -a_22,  ρ^2_x = a_21/λ(x) -a_22ρ^1_x,

N_x is defined in Proposition <ref> and a_ij are defined in (<ref>).


It is a direct consequence of (<ref>) and (<ref>).


Equation (<ref>) describes the dynamics of a large monomorphic population with trait x. Proposition <ref> shows that the age distribution of the population stabilizes around the equilibrium n_x =(n_x^1,n_x^2). The equilibria n^1_x and n_x^2 describe the age equilibria of the population with trait (x_b,x_d) and (x_b,0) respectively. We observe that if x_b<x_d, then  a_21=0= ρ^2_x and  Proposition <ref> leads to the equilibrium  n_x=(n_x^1,0). 



 §.§ Bimorphic dynamics

Let x=(x_b,x_d)∈𝒱 and y=(y_b,y_d)∈𝒱. We consider a bimorphic initial population such that the sequence (Z_0^K)_K weakly converges to δ_x n_x(a)da +δ_y n_y(a)da as K tends to infinity. Using similar arguments as in <cit.>, we can prove that the sequence of processes (Z^K)_K converges in probability,  on any finite time interval, to the solution ((n_x(t,.),n_y(t,.)),t≥ 0)∈ L^1(_+)^4 of the following system of non-linear partial differential equations

    ∂_t n_x(t,a) +∂_a n_x(t,a) = -( D_x(a) +η (‖ n_x(t,.)‖_1 +‖ n_y(t,.)‖_1)I)n_x(t,a)
    
    n_x(t,0)=∫__+B_x(α)n_x(t,α)dα
    ∂_t n_y(t,a) +∂_a n_y(t,a) = -( D_y(a) +η (‖ n_x(t,.)‖_1 +‖ n_y(t,.)‖_1)I)n_y(t,a)
    
    n_y(t,0)=∫__+B_y(α)n_y(t,α)dα,

with initial condition given by (n_x(0,a),n_y(0,a))=(n_x(a),n_y(a)). Equations (<ref>) describe the dynamics of a dimorphic population with traits x and y interacting through competition. We prove the following proposition.

Let x,y∈𝒱 such that x_b∧ x_d< y_b∧ y_d.  Then any L^1(_+)^4-solution (n_x(t,.),n_y(t,.)) of (<ref>)  satisfying n_y^1(t,.)≠ 0 converges to (0,n_y) in L^1(_+)^4 as t→∞.


Let x,y∈𝒱 such that x_b∧ x_d < y_b∧ y_d. Then by (<ref>), we have λ(x)<λ(y). For u∈{ x,y} we define

    v_u(t,a)=exp(η∫_0^t(‖ n_x(s,.)‖_1 +‖ n_y(s,.)‖_1)ds)n_u(t,a).

The functions v_x and v_y are solutions of the linear systems (<ref>). We deduce from Proposition <ref> that e^-λ(u) tv_u(t,.)→ c(v_u(0,.))N_u in L^1(_+)^2 as t→  + ∞, for a positive constant c(v_u(0,.)). Hence, we obtain that e^-λ(x) t∫__+v_x^1(t,α)dα converges to  a positive limit. Since λ(y)>λ(x) it comes that e^-λ(y) t∫__+v_x^1(t,α)dα converges to 0 as t→+∞. We deduce that

    ρ_y^1(t)/ρ_x^1(t)=e^-λ(y) t∫__+v_y^1(t,α)dα/e^-λ(y) t∫__+v_x^1(t,α)dα→ +∞.

Since ρ_y^1(t) is bounded we deduce that ρ_x^1(t)→ 0 and similarly that ρ_x^2(t)→ 0. Then the population with trait x becomes extinct. Using similar arguments as in the previous proof we obtain that n_y(t,.)→n_y in L^1(_+)^2 which allows us to conclude.


Equation (<ref>) describes a competition dynamics between two infinite monomorphic populations with trait x and y. Proposition <ref> shows that if x_b∧ x_d < y_b∧ y_d, then the population with trait y invades and becomes fixed while  the population with trait x becomes extinct. That  gives us an invasion-implies-fixation criterion.



§ ADAPTIVE DYNAMICS ANALYSIS

In this section, we study the model introduced in Section 2 under the different scaling of the adaptive dynamics. We generalise the Trait Substitution Sequence with age structure (<cit.>) to take into account the Lansing Effect. Then we study the behaviour of the TSS on a large time-scale when mutations are small. We show that the limiting behaviour of the TSS is described by  a differential inclusion which generalises the canonical equation of adaptive dynamics (<cit.>, <cit.>, <cit.>) to non-regular fitness  functions. We first state some properties of the demographic parameters and  introduce the invasion fitness function.


 §.§ Malthusian parameter and invasion fitness

Let us introduce the following sets: U_1={ x∈𝒱:x_b<x_d}, U_2={ x∈𝒱:x_d<x_b} and ℋ={ x∈𝒱:x_b=x_d}.



  §.§.§ Malthusian parameter

We now give some properties of the Malthusian parameter λ(x) defined in (<ref>). 



  (i)For all x∈𝒱,  0≤λ(x)<1.

  (ii) The map x∈𝒱↦λ(x) is continuous. It is differentiable on U_1∪ U_2 and satisfies

    ∀ x∈ U_1,  ∇λ(x)= (e^-λ(x) x_b/G(x), 0)
       ∀ x∈ U_2,  ∇λ(x)= (0,e^-λ(x)x_d/G(x))




where G(x)=∫_0^x_b∧ x_da e^-λ(x) ada. 

  (iii)We have sup_x∈ U_1∪ U_2‖∇λ(x)‖ < +∞. Moreover, for all i∈{ 1,2}, the fitness gradient ∇λ is Lipschitz on U_i.



Note that the Malthusian parameter λ(x) is not differentiable on the diagonal ℋ, which can be easily obtained by computing left and right partial derivatives on ℋ.


For all h∈ we define (h)_1 = ([ h; 0 ]) and (h)_2=([ 0; h ]).


(i) Since x_b∧ x_d >1 we obtain from the definition (<ref>) that λ(x)> 0. Assume that λ(x)≥ 1. Then we obtain that 1=∫_0^x_b∧ x_de^-λ(x) ada ≤ 1 - e^-x_b∧ x_d which is a contradiction.


(ii) We prove the continuity. Let x∈𝒱 and let (x_n) be a sequence of 𝒱 such that x^n→ x. By (i), λ is bounded and we can extract a subsequence (still denoted x^n by simplicity)  such that λ(x^n)→λ^*. We deduce that 1 = ∫_0^x^n_b∧ x^n_de^-λ(x^n)ada→∫_0^x_b∧ x_de^-λ^*ada=1 which allows us to conclude. Differentiability properties are a direct consequence of the Implicit Function Theorem. For each i∈{ 1,2} we apply implicit function Theorem to the map

    (λ,x)∈× U_i↦∫_0^x_b∧ x_de^-λ ada - 1.

We deduce that λ is differentiable over U_1∪ U_2 and that

    ∀ i∈{ 1,2},  ∀ x ∈ U_i,  ∇λ(x)= (e^-λ(x) (x_b∧ x_d)/G(x))_i.

(iii) It is straightforward to check that for all x∈𝒱,

    e^-λ(x)(x_b∧ x_d)/G(x)≤1/∫_0^1a e^-ada

which allows us to obtain that sup_x∈ U_1∪ U_2‖∇λ(x)‖ < +∞. Moreover, the gradient ∇λ is obviously differentiable on U_i. Since G is bounded below by ∫_0^1ae^-ada, we deduce that ∇λ has bounded derivatives on U_i and that ∇λ is Lipschitz on U_i.


Formulae (<ref>) describe the sensitivity of the Malthusian parameter to small variations of the trait x as well as the strength of selection at ages x_b and x_d for a population with Lansing effect. The quantity G(x) can be interpreted as the mean generation time associated with the trait x. Moreover (<ref>) and Proposition <ref> yield

    ∀ i∈{ 1,2},  ∀ x ∈ U_i,  ∇λ(x)= (N_x^1(x_b∧ x_d)/G(x))_i

where N_x^1 is the stable age distribution. In <cit.>, Caswell obtains similar formulae for derivatives of the Malthusian parameter with respect to some little perturbations on the intensity of birth or death at some given age while our formulae are obtained considering a small perturbation on the duration of the reproduction phase (not on the intensity which remains constant equal to one).

The following proposition recalls a simple link between the Malthusian parameter and the stationary state of the monomorphic partial differential equation (<ref>).



  (i) For all x∈𝒱, we have λ(x) = η ‖n_x‖_1. 
  (ii) The map x∈𝒱↦n_x(0) is continuous and bounded.



(i)  has been proved at the beginning of the proof of Proposition <ref>.


(ii) By (i) we have 

    λ(x) =η( n_x^1(0)u_1(x) +n_x^2(0)u_2(x))

where

    u_1(x)   =∫_0^+∞exp(-∫_0^a1_α>x_ddα -λ(x) a)da ; 
    u_2(x)=∫_0^+∞exp(-(1+λ(x))a)da.

Moreover n_x(0) is a solution of  

    n_x(0)=F(λ(x)) n_x(0)

where F is defined in (<ref>). From (<ref>) and (<ref>) we obtain by simple computation that

    n_x^1(0)=λ(x)/η1/u_1(x) + u_2(x)[F(λ(x))]_21/1-[F(λ(x))]_22

which is a continuous function of x. Boundedness is obvious arguing that n_x^1(0)≤‖n_x‖_1.



  §.§.§ Invasion fitness

We extend the definition of the invasion fitness for age-structured populations introduced in <cit.> to take into account the Lansing effect. 

For all y∈_+^* and x∈𝒱, the invasion fitness 1-z(y,x) is defined as the survival probability of a bi-type age structured branching process with birth rates and death rates defined in (<ref>), respectively equal to

    B_y(a)   and   D_y(a) +η ‖n_x‖I.


The next proposition gives a precise and precious relation between the invasion fitness and the Malthusian parameter.

Let y∈_+^* and x∈𝒱, then the invasion fitness satisfies

    1-z(y,x)=(λ(y) - λ(x))∨ 0.



Let Z_t(da)=(Z^1_t(da), Z_t^2(da)) be an age- structured branching process with birth rates B_y(a) and death rates D_y(a) + η ‖n_x‖I= D_y(a) +λ(x) I. The process Z_t becomes extinct if and only if the process Z_t^1 becomes extinct. Indeed, if Z_0^1 = 0, the process Z_t^2 evolves as a sub-critical branching process. The process Z_t^1 is an age structured branching process with birth rates and death rates respectively

    [B_y(a)]_11 = 1_a≤ y_b∧ y_d   and [D_y(a)]_11+λ(x) = 1_a>y_d+λ(x).

We deduce that z(y,x) equals the smallest solution of the equation z = F(z) where

    F(z)   =∫__+e^(z-1)∫_0^a[B_y(α)]_11dα([D_y(a)]_11+λ(x)) e^-∫_0^a([D_y(α)]_11+λ(x))dαda
        =∫_0^y_b∧ y_dλ(x) e^(z-1-λ(x))ada + e^(z-1)(y_b∧ y_d)∫_y_b∧ y_d^y_dλ(x) e^-λ(x) ada
        +e^(z-1)(y_b∧ y_d)(1+λ(x))∫_y_d^+∞e^-(a-y_d)-λ(x) ada
       =1 - e^(z-1-λ(x))(y_b∧ y_d)+(z-1)∫_0^y_b∧ y_de^(z-1-λ(x))ada
       +e^(z-1)y_b∧ y_d(e^-λ(x) (y_b∧ y_d) -e^-λ(x) y_d)+e^(z-1)(y_b∧ y_d) +e^(z-1)(y_b∧ y_d) e^-λ(x) y_d
       = 1 +(z-1)∫_0^y_b∧ y_d e^(z-1-λ(x))a da.

We have obtained that the equation z=F(z) is equivalent to 

    z-1 = (z-1)∫_0^y_b∧ y_d e^(z-1-λ(x))ada.

If λ(y)>λ(x), Equation (<ref>) admits two solutions z=1 and z=λ(x) -λ(y) +1<1. If λ(y)<λ(x), Equation (<ref>) admits a  unique solution z=1. That allows us to conclude the proof.


It is interesting to note that in our model, the invasion fitness (which is a concept from adaptive dynamics theory) and the Malthusian parameter are connected thanks to the simple relation (<ref>).

The following proposition characterises the set of traits y which can invade some given trait x.

For all x∈𝒱, y∈_+^2, 

    1-z(y,x)>0 ⟺λ(y)>λ(x)⟺  y_b∧ y_d > x_b∧ x_d.



From the definition (<ref>) of the Malthusian parameter, we easily deduce the following equivalences:

    (λ(y)-λ(x))∨ 0>0    ⟺λ(y)-λ(x)>0
       ⟺∫_0^x_b∧ x_de^-λ(y)ada <1 
       ⟺ x_b∧ x_d<y_b∧ y_d,

which concludes the proof.



 §.§ Trait Substitution Sequence with age structure

We first generalise the definition of the  Trait Substitution Sequence (TSS) with age structure defined in <cit.> to take into account the Lansing Effect. 

We define the measure valued process
(T_t(dx,da),t≥ 0) by

    T_t(dx,da)=δ_X(t)(dx) n_X(t)^1(a)da +δ_(X_b(t),0)(dx) n_X(t)^2(a)da

where (X(t),t≥ 0)= ((X_b(t),X_d(t)),t≥ 0) is defined as the pure jump Markov process on 𝒱 with infinitesimal generator L defined for all measurable and bounded function φ:𝒱→ and x∈𝒱 by

    Lφ(x)=∫_(_+)^2(φ(x+h)-φ(x))(λ(x+h)-λ(x))n_x^1(0)μ(dh)

where

    μ(dh)=k(h_b)dh_b⊗δ_0(dh_d) +δ_0(dh_b)⊗ k(h_d)dh_d/2,

k being defined in (<ref>). 

The process X will be called the Trait Substitution Sequence. 

The process (T_t,t≥ 0) describes the evolution of the phenotypic structure of the population at the mutational time-scale. At each time, and because of the Lansing effect, the population is composed of two sub-populations: the first one corresponds to viable individuals with trait X(t)=(X_b(t), X_d(t))∈𝒱 whose age distribution is given by n_X(t)^1(a)da; the second one is composed of individuals generated through the Lansing effect, with trait (X_b(t), 0) and age distribution n_X(t)^2(a)da.


Figure <ref> describes the behaviour of the process (X(t),t≥ 0). 

Any trait x∈ℋ is an absorbing state for the process X. Indeed, by Proposition <ref>, for all φ:𝒱→ measurable and bounded, for all x∈ℋ, Lφ(x)=0. That means that when the trait (x_b,x_d) satisfies x_b=x_d, no mutation can invade. 

By definition of the measure μ(dh) the process evolves horizontally or vertically (which means that the two traits x_b and x_d do not mutate simultaneously).  By Proposition <ref>, we deduce easily that the process evolves from the left to the right on U_1 and from bottom to top on U_2. 

Since the jump rates are continuous and tend to zero on ℋ, the process slows down as it approaches ℋ. 


We now explain the heuristics, rigorously proved in <cit.>, which allow to obtain the TSS from the individual based model defined in Section 2.1. The main ideas  have been introduced in <cit.>, for a population without age-structure. They are based on the time-scale separation assumption on mutation probability p_K: as K→ +∞

    ∀ V>0,  exp(-KV)=o(p_K),   p_K=o(1/Klog(K)),

which allows to separate the effect of the natural selection and the appearance of new mutants. Let x∈𝒱 and consider a sequence (Z_0^K)_K  converging to δ_x n_x^1(0,a)da as K→ +∞.



1) Monomorphic approximation. For large K, the process Z_t^K stays close to the measure δ_x n_x^1(t,.) + δ_(x_b,0) n_x^2(t,.) where n_x(t,.)=(n_x^1(t,.),n_x^2(t,.)) satisfies the partial differential equation (<ref>). By Proposition <ref>, the dynamics n_x(t,.) converges to n_x= (n_x^1,n_x^2) as t tends to infinity and hence reaches a given neighbourhood of n_x in finite time. By using large deviation results (<cit.>), we obtain with probability tending to one as K tends to infinity that the process Z_t^K stays in this neighbourhood of δ_x n_x^1 + δ_(x_b,0)n_x^2 during a time e^CK for some C>0. The left-hand side in Assumption (<ref>) ensures that the next mutation appears before the process leaves this neighbourhood.




2) Appearance of a mutant.
We deduce that the monomorphic population with trait x creates a mutant with trait: 


  (i) y=(x_b +h_b,x_d) or y=(x_b ,x_d +h_d) at a rate approximately equal to
2K p_K(1-p_K)n^1_x(0);

  (ii)y=(x_b +h_b,x_d+h_d) at a rate approximately equal to Kp_K^2 n_x^1(0);

  (iii)y=(x_b, h_d) or y=(x_b+h_b, 0) at  at rate approximately equal to 2Kp_K(1-p_K)n_x^2(0);

  (iv)y=(x_b +h_b, h_d) at rate approximately equal to Kp_K^2 n_x^2(0).

where the variables h_b and h_d are chosen independently with distribution k. 



Since p_K^2=o(p_K(1-p_K)), the cases (ii) and (iv) cannot be observed on the mutation time-scale t/2K(1-p_K).



3) Effect of the natural selection.
In cases (i) and (iii), the mutant population dynamics is approximated by a bi-type age structured branching process with birth rates B_y and death rates D_y +η‖n_x‖_1 I. By Proposition <ref>, the mutant population survives with probability

    (λ(y)-λ(x))∨ 0.

Let us detail the two different cases.


  • Case (i). With  probability 1/2, the trait of the mutant is  y=(x_b+h_b, x_d). From Proposition <ref>, we deduce that the mutant can survive if and only if

    λ(x_b+h_b,x_d)>λ(x)⟺ h_b>0 and x∈ U_1.

With probability 1/2, y=(x_b,x_d+h_d) and can survive if and only if

    λ(x_b,x_d +h_d)>λ(x)⟺ h_d>0 and x∈ U_2.


  • Case (iii) (Lansing effect). The mutant has the trait y=(x_b, h_d) or y=(x_b +h_b, 0). By (<ref>), we have y_b∧ y_d<1 which implies that λ(y)<0 and then (λ(y)-λ(x))∨ 0 =0. In this case,  the mutant population becomes extinct.

We deduce that the mutant can only survive (with positive probability) in case (i). The birth rate of such a mutant (on the time-scale t/2p_K(1-p_K)) is given by the intensity measure on 

    n_x^1(0)μ(dh)

that leads to the right hand side in (<ref>). The probability that such a mutant survives and reaches a size of  order K equals 

    (λ(y)-λ(x))∨ 0.
 
Moreover (<ref>) and (<ref>) imply that 

    ((λ(x+h)-λ(x))∨ 0)μ(dh)= (λ(x+h)-λ(x))1__+^2(h)μ(dh).

This allows to obtain the left hand side of (<ref>).

If the mutant population becomes extinct, the resident population remains close to its equilibrium n_x.

If the mutant population survives, then it reaches a size of order K with a probability that tends to one and the  population dynamics is approximated by the solution (n_x(t,.),n_y(t,.)) of the bimorphic system of partial differential equations (<ref>). In this case we have necessarily λ(y)>λ(x). Following Proposition <ref>, the deterministic dynamics (n_x(t,.),n_y(t,.))  reaches a neighbourhood of (0,n_y). By using branching processes approximations and  arguments introduced in <cit.>, we can deduce that the resident population with trait x becomes extinct. One can prove as in <cit.> that this competition phase has a duration of order log(K).

The right hand side of Assumption (<ref>) ensures that the three steps of invasion are completed before the next mutation occurs. 

•

The Markov property allows to reiterate the same reasoning for the next mutation occurence.  In summary, the following theorem holds. 

The following convergence holds in the sense of finite dimensional marginals:

    (Z_t/2Kp_K(1-p_K)^K, t≥ 0)⟶ (T_t, t≥ 0),   as K→∞,

where the process T is defined in Definition <ref>.



 §.§ A canonical inclusion for adaptive dynamics

In this section, we assume in that mutations are small. We study the behaviour of the process X defined in (<ref>) when the mutation scale equals ϵ>0 and the time is rescaled by 1/ϵ^2. To this aim, we define the rescaled trait substitution sequence process X^ϵ and study the limiting behaviour of the process X^ϵ as ϵ→ 0. In the usual cases (smooth fitness functions) the canonical equation introduced by Dieckmann-Law can be derived as limit of X^ϵ as ϵ→ 0 (<cit.>). As observed in Section 4.1, the fitness function λ(x) does not satisfy these regularity assumptions. To overpass this difficulty, we use the approach developed in <cit.> based on differential inclusions. We prove in Theorem <ref> that the set of limit points of the family X^ϵ is characterised as the set of solutions of a differential inclusion.

The rescaled TSS process (X^ϵ(t),t≥ 0) is defined as a pure jump Markov process with infinitesimal generator L^ϵ defined for all measurable and bounded function φ:𝒱→ and x∈𝒱 by

    L^ϵφ(x)=1/ϵ^2∫_(_+)^2(φ(x+ϵ h)-φ(x))(λ(x+ϵ h)-λ(x))n_x^1(0)μ(dh)



The process X^ϵ shows a dynamics similar to the process X. The jump rates are of order 1/ϵ and the jump sizes are of order ϵ.

We first introduce the set-valued map F:𝒱→𝒫(_+^2) defined for any x∈𝒱 by

    ∀ i ∈{ 1,2},     ∀ x∈ U_i,   F(x)=( f(x,1) )_i 
       ∀ x∈ℋ,   F(x)={1/2([ f(x,u); f(x,u) ]), u∈[0,1]},

where for all (x,u)∈𝒱×[0,1],

    f(x,u)=(∫_0^uh^2k(h)dh +∫_u^1 h u k(h)dh)e^-λ(x) (x_b∧ x_d)/G(x)n_x^1(0)/2;

and G is defined in Proposition <ref>.



This set-valued map F somehow generalises the classical fitness gradient. It is represented by a picture in Figure <ref> (b). 

Let us explain the ideas leading to consider this function. Let us consider a compact subset K of U_i. Since the Malthusian parameter λ is differentiable on U_i, the following approximation holds: for all h∈_+, uniformly for x∈ K, we have

    λ(x+ϵ (h)_i)-λ(x)≈ϵ (h)_i.∇λ(x),

which leads to the definition of F on U_i. We analyse the case x∈ℋ for which the approximation (<ref>) is not true. Indeed, let x∈ℋ and let u∈[0,1] and let us consider a sequence x^ϵ=x-ϵ(u)_i, we have 

    λ(x^ϵ +ϵ (h)_i)-λ(x^ϵ)   =λ(x +ϵ(h-u)_i) -λ(x-ϵ(u)_i)
       =λ(x +ϵ(h-u)_i)  -λ(x) +λ(x)-λ(x-ϵ(u)_i)

Assume h<u, since λ is differentiable on U_i we obtain when ϵ tends to 0 that

    λ(x^ϵ +ϵ (h)_i)-λ(x^ϵ) ≈ϵ (h-u)_i.∇λ(x) + ϵ(u)_i.∇λ(x) ≈ϵ (h)_i.∇λ(x),

where ∇λ(x) is defined as the limit of ∇λ(y), y→ x, y∈ U_i. That leads to the first integral in (<ref>).
If h>u, we obtain that  λ(x +ϵ(h-u)_i)  -λ(x)=0 and 

    λ(x^ϵ +ϵ (h)_i)-λ(x^ϵ)≈ϵ(u)_i.∇λ(x)≤ϵ(h)_i.∇λ(x),

which leads to the second integral in (<ref>).
The inequality above means that when the process evolves near the diagonal ℋ the adaptation slows down. Let us now define what we do  mean by a solution of the differential inclusion driven by the set-valued map F.

A solution of the differential inclusion driven by the set-valued map F is an absolutely continuous function x :[0,T]→𝒱 which satisfies x(0)=c^0 and for almost all t∈[0,T],

    dx(t)/dt∈ F(x(t)).


It generalises the classical canonical equation for adaptive dynamics. For any T>0 and x^0∈𝒱, we denote by 𝒮_F(T,x^0) the set of solutions of the differential inclusion (<ref>). The following theorem characterises the limit of the process X^ϵ as the solution of the differential inclusion (<ref>).

Let x^0∈𝒱. Assume that X^ϵ(0)→ x^0 in probability as ϵ→ 0. For all T,δ>0,

    lim_ϵ→ 0ℙ( inf_x∈𝒮_F(T,x^0)sup_t∈[0,T]| X^ϵ(t) -x(t)| >δ) = 0.



Theorem <ref> justifies our complete study. Let (x(t),t∈[0,T]) be a solution of (<ref>). On each U_i, it satisfies

    dx(t)/dt=(f(x(t),1))_i,   x(t)∈ U_i.

The map x∈ U_i↦ (f(x,1))_i is Lipschitz and bounded below by a positive constant. Hence, uniqueness holds on U_i for (<ref>) and any solution  reaches in finite time the diagonal ℋ. On ℋ, the solution satisfies

    dx(t)/dt∈ F(x(t)),   x(t)∈ℋ.

Since for all x∈ℋ, F(x)⊂ℋ, we deduce that any solution stays in ℋ. 

Figure <ref> illustrates Theorem <ref>. We represent some trajectories of the process (X^ϵ(t),t≥ 0) started at X^ϵ(0)=(2, 1.5) for ϵ=0.001. 


The proof of Theorem <ref> is based on <cit.>  recalled in the Appendix, see Theorem <ref>. We start by writing the process X^ϵ as a time-changed Markov chain. We first re-write the generator L^ϵ as follows

    L^ϵφ(x) = ∫__+^2(φ(x+h)-φ(x))n_x^1(0)/ϵk^ϵ(x,dh)

where k^ϵ(x,dh)=

    λ(x+h)-λ(x)/ϵμ^ϵ(dh) +(∫_(_+)^2(1-λ(x+g)-λ(x)/ϵ)μ^ϵ(dg))δ_0(dh)

and μ^ϵ is the image measure of  μ by the map h↦ϵ h. The following lemma is proved in <cit.>.
[<cit.>]
Let τ = sup_x∈𝒱n_x^1(0). Let (Y^ϵ(k),k≥ 0) be a Markov chain with jump law

    k̃^ϵ(x,dh)=n_x^1(0)/τk^ϵ(x,dh) + (1-n_x^1(0)/τ)δ_0(dh)

and Λ^ϵ be a Poisson process with intensity τ/ϵ. Then the processes X^ϵ and Y^ϵ(Λ^ϵ) have the same law.

Then we are led to study the Markov chain (Y^ϵ(k),k≥ 0). We first define the drift of the Markov chain Y^ϵ by 

    g_ϵ(x)=𝔼[Y^ϵ(1)-Y^ϵ(0)| Y^ϵ(0)=x].

A simple calculation gives us 

    ∀ i ∈{ 1,2},     ∀ x∈ U_i,   g_ϵ(x)=ϵn_x^1(0)/2τ∫__+λ(x+ϵ(h)_i)-λ(x)/ϵk(h) (h)_i dh 
       ∀ x∈ℋ,   g_ϵ(x)=0.
 
Then, we write the Markov chain Y^ϵ as a stochastic approximation algorithm

    Y^ϵ(k+1)=Y^ϵ(k) + ϵ/τU^ϵ(k) + g^ϵ(Y^ϵ(k))

where U^ϵ is a martingale difference sequence. Assumptions of Theorem <ref>  are clearly satisfied. In order to apply it, we compute the following set-valued map

    ∀ x∈𝒱,   H(x)=conv{acc_ϵ→ 0τ g_ϵ(x^ϵ)/ϵ:x^ϵ→ x}

where conv(A) denotes the smallest convex set which contains A and acc_ϵ→ 0τ g_ϵ(x^ϵ)/ϵ is the set of accumulation points of the sequence τ g_ϵ(x^ϵ)/ϵ as ϵ tends to zero. 

The set-valued map H satisfies

    ∀ i ∈{ 1,2},     ∀ x∈ U_i,   H(x)={( f(x,1) )_i}
       ∀ x∈ℋ,   H(x)={α([ f(x,u);      0 ])+(1-α)([      0; f(x,v) ]):(u,v,α)∈[0,1]^3}

where 

    f(x,u)=e^-λ(x)(x_b∧ x_d)/G(x)n_x^1(0)/2(∫_0^uh^2 k(h)dh + ∫_u^1u h k(h)dh).



Let i∈{ 1,2}. Let K be a compact subset of U_i. Let δ>0. We fix ϵ_0>0 such that for all x∈ K, ϵ <ϵ_0 and h∈[0,1], we have x+ϵ(h)_i∈ U_i. The map λ is differentiable on U_i. Hence, for all (x,ϵ,h)∈ K×[0,ϵ_0]×[0,1], there exists θ∈[x,x+ϵ (h)_i] such that 

    λ(x+ϵ (h)_i)-λ(x) = ϵ (h)_i.∇λ(θ).

Let x∈ K, we have

    ‖τ g^ϵ(x)/ϵ - (f(x,1))_i‖    =n^1_x(0)/2|∫__+[∇λ(θ)]_i h^2 k(h)dh - [∇λ(x)]_i∫__+h^2 k(h)dh|
       ≤sup_x∈𝒱n_x^1(0)/2∫__+|[∇λ(θ)]_i - [∇λ(x)]_i| h^2 k(h)dh.

By Proposition <ref>, the map ∇λ is Lipschitz on U_i with some Lipschitz constant C. We deduce that

    ‖τ g^ϵ(x)/ϵ - (f(x,1))_i‖≤ϵsup_x∈𝒱n_x^1(0)C/2∫__+h^3k(h)dh.

We obtain that τ g^ϵ/ϵ converges uniformly on all compact subsets of U_i and we conclude that for all x∈ U_i, H(x)=(f(x,1))_i.


Let x∈ℋ. We first show that 

    {acc_ϵ→ 0τ g_ϵ(x^ϵ)/ϵ:x^ϵ→ x} = ⋃_i∈{ 1,2}{ (f(x,u))_i:u∈[0,1]}.

If the fitness gradient was smooth (i.e Lipschitz), the sets in (<ref>) would be reduced to one element and we would be in the cases studied in <cit.>, <cit.> or <cit.>. We prove the inclusion from right to left. Let i∈{ 1, 2}, let u∈[0,1], we define the sequence x^ϵ= x -ϵ(u)_i. We have

    τ g^ϵ(x^ϵ)/ϵ    =(∫__+λ(x -ϵ(u)_i +ϵ (h)_i) -λ(x-ϵ(u)_i)/ϵk(h)(h)_i dh)n_x^ϵ^1(0)/2
       =(∫_0^uλ(x -ϵ(u)_i +ϵ (h)_i) -λ(x-ϵ(u)_i)/ϵk(h)(h)_i dh .
       . + ∫_u^1λ(x -ϵ(u)_i +ϵ (h)_i) -λ(x-ϵ(u)_i)/ϵk(h)(h)_i dh)n_x^ϵ^1(0)/2.

For all h∈[0,u], we have [x-ϵ(u)_i,x - ϵ(u-h)_i]⊂ U_i. So we can find θ∈[x-ϵ(u)_i,x - ϵ(u-h)_i] such that λ(x -ϵ(u-h)_i) -λ(x-ϵ(u)_i) = ϵ (h)_i.∇λ(θ). By Proposition <ref> (ii), we deduce that 

    ∫_0^uλ(x -ϵ(u)_i +ϵ (h)_i) -λ(x-ϵ(u)_i)/ϵk(h)(h)_i dh →e^-λ(x) (x_b∧ x_d)/G(x)∫_0^u(h^2)_ik(h)dh

as ϵ tends to zero. For all h∈[u,1] we have λ(x +ϵ(h-u)_i)=λ(x). We deduce similarly that 

    ∫_u^1λ(x -ϵ(u)_i +ϵ (h)_i) -λ(x-ϵ(u)_i)/ϵk(h)(h)_i dh →e^-λ(x) (x_b∧ x_d)/G(x)∫_u^1 (h)_i uk(h)dh
 
as ϵ tends to zero. We conclude the proof of the first inclusion arguing that n_x^ϵ^1(0)→n_x^1(0) as ϵ tends to zero. 
We prove the inclusion from left to right. Let x∈𝒱. If x∈ℋ then we have 

    g^ϵ(x)=0.
 
If x∈ U_i, for some i∈{ 1,2}, then we have

    τ g^ϵ(x)/ϵ=n_x^1(0)/2∫_0^1λ(x+ϵ (h)_i)-λ(x)/ϵk(h)(h)_i dh

Moreover for all ϵ,h there exists θ∈[x,x+ϵ(h)_1] such that λ(x+ϵ (h)_i)-λ(x)≤∂λ(θ)/∂ x_iϵ h. We deduce that 

    [τ g^ϵ(x)/ϵ]_i≤n_x^1(0)/2∫_0^1∂λ(θ)/∂ x_i h^2 k(h) dh

From (<ref>) and (<ref>), we deduce easily the second inclusion in (<ref>). We conclude that

    conv{acc_ϵ→ 0τ g_ϵ(x^ϵ)/ϵ:x^ϵ→ x} = conv{∪_i∈{ 1,2}{ (f(x,u))_i:u∈[0,1]}}=H(x).
 

In Lemma <ref>, we prove that differential inclusions associated with H and F have identical solutions. Before, we give a technical lemma.

Let x∈ℋ. Let (u,v,α)∈[0,1]^3 and let 

    m=α([ f(x,u);      0 ]) +(1-α)([      0; f(x,v) ])∈ H(x).

Assume that m∉ F(x). Then we have α  f(x,u)≠ (1-α) f(x,v).


We prove the lemma by contradiction. Assume that α  f(x,u)= (1-α) f(x,v). Then we obtain 

    m=1/2([ 2f(x,u)f(x,v)/f(x,u)+f(x,v); 2f(x,u)f(x,v)/f(x,u)+f(x,v) ]).

We assume without loss of generality that f(x,u)≤ f(x,v). Then we obtain 0≤2f(x,u)f(x,v)/f(x,u)+f(x,v)≤ f(x,v). Finally we remark that the map s↦ f(x,s) is a bijection from [0,1] to [0,f(x,1)]. Hence there exists w∈[0,1] such that f(x,w)=2f(x,u)f(x,v)/f(x,u)+f(x,v) that allows us to obtain the contradiction.


Any solution of 

    dx(t)/dt∈ F(x(t)),   t∈[0,T]
    
    x(0)=x^0.

is a solution of 

    dx(t)/dt∈ H(x(t)),   t∈[0,T]
    
    x(0)=x^0.

and conversely.


For all x∈𝒱, we have F(x)⊂ H(x). We deduce that if (x(t),t∈[0,T]) is a solution of (<ref>) then it is a solution of (<ref>). Assume conversely that there exists a solution (x(t),t∈[0,T]) of (<ref>) which is not a solution of (<ref>). We deduce that there exists t_0∈[0, T] such that x is differentiable at t_0, dx(t_0)/dt∈ H(x(t_0)) and dx(t_0)/dt∉ F(x(t_0)). Then we have  x(t_0)∈ℋ and dx(t_0)/dt∉ F(x(t_0)).

We now deduce the contradiction. By Lemma <ref>, we obtain that
dx_b(t_0)/dt≠dx_d(t_0)/dt. 
Without loss of generality, we may assume that 
dx_b(t_0)/dt < dx_d(t_0)/dt.
Since x_b(t_0)=x_d(t_0), there exists an interval ]t_0,t_1[ such that for all s∈]t_0,t_1[, x(s)∈ U_1. Assume that for all s∈]t_0,t_1[ such that x is differentiable at s, we have 

    dx_d(s)/dt=0.

The solution x of the differential inclusion (<ref>)  is absolutely continuous. Hence for all s∈]t_0, t_1[,  x_d(s)=x_d(t_0)=x_b(t_0). Since  x(s)∈ U_1, we obtain that x_b(s)<x_b(t_0). It is absurd since x_b is non-decreasing and hence it contradicts (<ref>). So, there exists s_0∈]t_0,t_1[ such that x(s_0)∈ U_1, x is differentiable at s_0  and satisfies
dx_d(s_0)/dt>0.
However, x is a solution of (<ref>) that leads to the final contradiction.


We now give the proof of Theorem <ref>. It is a direct consequence of <cit.> recalled in Appendix A.3.

Let T>0 be fixed. By Lemma <ref> and Theorem <ref>, we deduce that for all δ>0 we have

    lim_ϵ→ 0ℙ( inf_x∈𝒮_F(T,x^0)sup_t∈[0,T]| Y^ϵ(⌊ tτ/ϵ⌋) -x(t)| >δ) = 0.

We conclude by using similar arguments as in the proof of <cit.>. Since Λ^ϵ is a Poisson process with parameter τ/ϵ, we obtain that for all δ>0,

    ℙ(sup_t≤ T|Λ^ϵ(t) -tτ/ϵ|≥τδ/ϵ)≤Tϵ/τδ.

We have ℙ(inf_x∈𝒮_F(T,x^0)sup_t≤ T| X^ϵ(t) - x(t)| > δ)

    =ℙ(inf_x∈𝒮_F(T,x^0)sup_t≤ T| Y^ϵ(Λ^ϵ(t)) - x(t)| > δ)
       ≤ℙ(inf_x∈𝒮_F(T,x^0)sup_t≤ T{| Y^ϵ(Λ^ϵ(t))-x(ϵΛ^ϵ(t)/τ)| + | x(ϵΛ^ϵ(t)/τ) -x(t)|}>δ).

Let x∈𝒮_F(T,x^0) be a solution of the differential inclusion (<ref>). For almost all t∈[0,T], we have dx(t)/dt∈ F(x(t)). Since sup_x∈𝒱sup{ F(x)} < +∞, we deduce that there exists C_T>0 such that for all y∈𝒮_F(T,x^0), for all s,t∈[0,T], | y(t) - y(s)|≤ C_T | t - s|. We deduce that (<ref>) is less than

    ℙ(inf_x∈𝒮_F(T,x^0){sup_t≤ T| Y^ϵ(Λ^ϵ(t))-x(ϵΛ^ϵ(t)/τ)|} +  C_T sup_t≤ T|ϵΛ^ϵ(t)/τ - t| >δ)
       ≤ℙ(inf_x∈𝒮_F(T,x^0)sup_t≤ T| Y^ϵ(Λ^ϵ(t))-x(ϵΛ^ϵ(t)/τ)|>δ) +  ℙ( C_T sup_t≤ T|ϵΛ^ϵ(t)/τ - t| >δ)

and we conclude by using (<ref>) and (<ref>).




§ DISCUSSION



 §.§ General discussion and comments

In the present article, we study the evolution of a population with a trait structure describing a simple class of life-histories. We build a stochastic individual-based model in a framework that is continuous in time, age and trait. The trait is a pair of parameters (x_b,x_d) characterising the age at end-of-reproduction x_b and the age at transition to a non-zero mortality risk x_d.  The model sees two origins of phenotypic variation. First, the genetic mutations that are supposed to be rare and modify the traits symmetrically - equal probability to increase or decrease the value of the parameter. Second, we model the Lansing effect which can be considered as an epigenetic mutation affecting the progeny of an “old” individual. It is acting on a much faster time-scale - one generation - than genetic mutations and has only a negative effect on the life expectancy of the progeny. We highlight that we have chosen to model the Lansing effect by an extremely strong effect on the descendant. Indeed, it acts on each generation and degrades dramatically the life-expectancy of the progeny. Nevertheless, a more general and realistic study of epigenetic modifications on the genetic evolution has to be investigate. Some aspects of this question have been studied in <cit.>. The reasoning is based on the fact that epigenetic modifications are more frequent than genetic mutations (<cit.>). It would be interesting to extend adaptive dynamics tools in order to take into account epigenetics modifications (involving intermediary time scales).

We study the long term evolution of the trait distribution using adaptive dynamics theory. We extend the age-structured TSS to take into account the Lansing effect. The main mathematical result of the present work concerns the behaviour of the TSS when mutations are small. We show in Theorem <ref> that the behaviour of the rescaled TSS in the limit ϵ→ 0 is characterised by a differential inclusion whose solutions are not unique on the diagonal ℋ={ x_b=x_d}. This differential inclusion allows to generalize the Canonical Equation of Adaptive Dynamics <cit.>,<cit.> in order to consider the non-smooth fitness gradient. The proof is based on <cit.>. Thanks to this approach, we show that the evolution of our model, whatever its initial configuration, leads to the apparition and maintenance of configurations (x_b,x_d) satisfying x_b-x_d=0.

To our knowledge, differential inclusions have never been used before in the adaptive dynamics theory. In <cit.>, <cit.>,  <cit.>, the fitness gradient is assumed to be a Lipschitz function, which ensures the uniqueness of the solutions. Our approach seems useful for generalizing the canonical equation to situations where the fitness gradient is neither Lipschitz nor continuous.  

The drift associated with the differential inclusion depends  on the fitness gradient that satisfies (see Proposition <ref>) 

    ∀ i∈{ 1,2},  ∀ x∈ U_i,  ∇λ(x)=(e^-λ(x)(x_b∧ x_d)/G(x))_i

where λ(x) is the Malthusian parameter and G(x) the mean generation time associated to the trait x. Hence the fitness gradient ∇λ(x) describes the speed of evolution of the trait x. It can be related to the seminal work of Hamilton (<cit.>) on the moulding of senescence. In that article, Hamilton states that senescence is unavoidable because the strength of selection decreases with age. To show it, he defines the strength of selection at some given age a_0 as the sensitivity of the Malthusian parameter with respect to some little perturbation on the birth or death intensities at age a_0. He concludes arguing that these quantities decrease to zero as a_0 tends to infinity. Similarly, formula (<ref>) describes the sensitivity of the Malthusian parameter with respect to some perturbation on the duration of the reproduction or the survival phase at ages x_b and x_d. Then, they can be interpreted as the strength of selection at ages x_b and x_d and describe the speed of evolution of the traits x_b and x_d in the canonical inclusion (<ref>).



To conclude, the present article studies a  case of  birth and death model with a strong Lansing effect and constant competition applied to asexual and haploid individuals in order to validate mathematically the convergence of x_b and x_d observed in the numerical simulations. Further characterisation of this model is in progress, in order to better understand the influence of its different parameters on the evolution of (x_b, x_d).

Our initial motivation for developing the bd-model was an attempt to understand whether a phenomenon leading to a dramatic decrease of an individual’s fitness could be selected through evolution with simple and no explicitely constraining trade-offs. Indeed, in the past years, Rera and collaborators have identified and characterised a dramatic transition  preceding death in drosophila (<cit.>, <cit.>) as well as other organisms (<cit.>, <cit.>). We show here that, under uniform competition - i.e. environmental limitation equally affecting all genotypes - a mechanism coupling the end of reproductive capabilities and organismal homeostasis can and will be selected thanks to evolution. Thus, at least these two characteristics of senescent organisms can positively be selected through evolution. Concerning the biological interpretation of our model, our thesis is that individuals with a senescence mechanism associated to the Lansing effect tend to produce more genetic variants than those without senescence. Hence, these individuals could show a higher evolvability. This question is being investigated in a work in progress.

Here we show the positive selection of a property that limits reproduction. It is reminiscent of the article by <cit.> 
who proposed a new selective mechanism for post-reproductive life span. The latter relies on the hypothesis that the post-reproductive lifespan can be selected as an insurance against indeterminacy; a longer life expectancy reducing the risk of dying by chance before the cessation of reproductive activity. In the present article, the maintenance of individuals showing Lansing effect is the counterpart of individuals with post-reproductive survival. As discussed in <cit.>, one of August Weismann’s concepts that  persisted without changes throughout his life is a conviction that “life is endowed with a fixed duration, not because it is contrary to its nature to be unlimited, but because unlimited existence of individuals would be a luxury without any corresponding advantage” (<cit.>). This is what we showed in the present work.




 §.§ Generalizations of the model



  §.§.§ Lower Lansing Effect

In this article, we model the Lansing effect by a very strong effect since we assume that the descendants of old individuals with trait (x_b,x_d) inherits of the trait (x_b,0). It could be interesting to study a more general case like (x_b,α x_d) for some α∈[0,1]. 


  §.§.§ Mutation kernel


It would be interesting to consider two different mutation kernels k_b(x,h) and k_d(x,h) for the traits x_b and x_d.  This change should not modify  the behaviour of the process. Indeed, on the sides off the diagonal (x_b ≠ x_d), only the speed of evolution of the traits will be modified. Hence, the trait (x_b,x_d) will converge to the diagonal and will evolve on it at some speed depending on the variances of the kernels k_b and k_d. This case is being studied in an ongoing work mentionned in Section <ref>.
 
 
 


§ APPENDIX



 §.§ Proof of Proposition <ref>

The proof is based on classical arguments of spectral theory for strongly
continuous semi-groups. Let us denote by (P_t)_t≥ 0 the semi-group on L^1(_+)^2 associated with the infinitesimal generator (A,D(A)). Let us denote by σ(A) and σ_e(A) the spectrum and the essential spectrum of the operator A respectively. Let us denote by α[P_t] the measure of non-compactness of P_t <cit.>, and define w_1(A):=lim_t→∞t^-1log(α[P_t]). We show that there exists ω>0 such that

    max(w_1(A),sup_z∈(σ(A)∖σ_e(A))∖{λ(x)}ℛe(z))<ω<λ(x).

By using arguments similar to <cit.>, we obtain that for all t large enough

    α[P_t]≤ e^-t.

and that

    w_1(A):=lim_t→∞t^-1log(α[P_t])≤ -1.

By <cit.>, we obtain that σ_e(A)⊂{ z∈ℂ:ℛe(z)≤ -1}. Let z ∈σ(A)∖σ_e(A). Then there exists a non-zero u∈ L^1(_+,ℂ)^2 such that

    -u'(a) -D_x(a)u = z u(a)
    
    u(0)=∫__+B_x(α)u(α)dα.

By solving the first equation in (<ref>) and by injecting the result in the second equation, we obtain that u(0) satisfies

    u(0)= F(z)u(0)

where 

    F(z) = ∫__+B_x(a)exp( -∫_0^a(D_x(α)+ z I)dα)da.

Equation (<ref>) admits a non-trivial solution u(0) if and only if [F(z)-I] =0. Since the matrix F(z) is triangular we have

    [F(z)-I]= ([F(z)]_11 - 1)([F(z)]_22-1)

where [F(z)]_11=∫_0^x_b∧ x_de^-zada and [F(z)]_22=∫_0^x_be^-(1+z)ada.
We deduce that the Malthusian parameter λ(x) is the largest real solution of [F(z)-I]=0. By using similar analytical arguments as in the proof of <cit.>, we deduce that there exists only finitely many z∈ℂ such that ℛe(z)>0 and [F(z)-I]=0 which allows us to conclude for (<ref>). We prove that λ(x) is a simple eigenvalue of A by showing that λ(x) is a simple zero of the equation [F(z)-I]=0. Indeed we have 

    d[F(λ(x))-I]/dλ=d[F(λ(x))]_11/dλ([F(λ(x))]_22-1)>0.

Let N_x be a principal eigenvector associated with the eigenvalue λ(x). By (<ref>), we have

    N_x^1(a)=N_x^1(0)e^-((a-x_d)∨ 0) -λ(x) a,   N_x^2(a)=N_x^2(0) e^-(1+λ(x))a

and Equation (<ref>) gives that 

    N_x^2(0)= [F(λ(x))]_21/1-[F(λ(x))]_22N_x^1(0).

We conclude for the convergence by using  arguments similar to proof of <cit.>. 


 §.§ Proof of Lemma <ref>


Equation (<ref>) has the form du/dt=f(t,u) with f(t,u)→ g(u) as t→∞. So (<ref>) is called an asymptotically autonomous differential equation (<cit.>, <cit.>) with the limit equation

    dy(t)/dt= My(t) - η‖ y(t)‖_1 y(t).

We first show that any solution y(t) of (<ref>) started at y(0)∈_+^*×_+ converges to a stationary state. In <cit.>, the proof is given when M is irreducible. We give a slightly different proof. Let z be defined in (<ref>). It is straightforward to prove that z is the eigenvector of M associated with the simple eigenvalue m_11, which satisfy the condition ‖z‖_1=m_11/η. Since z(0)∈_+^*×_+ there exists a positive constant c(z(0)) such that e^-m_11te^Mtz(0)→ c(z(0))z as t→∞. We now write

    y_1(t)/y_2(t)=[e^-m_11te^Mty(0)]_1/[e^-m_11te^Mty(0)]_2→z_1/z_2.

We deduce that the ω-limit set of any solution of (<ref>) is a subset of Δ = { z∈_+^2:z_1=z_1/z_2z_2}. We conclude by proving that any solution starting from Δ converges to z. Let us consider such a solution (always denoted by y(t)). We have 

    dy_1(t)/dt   =y_1(t)(m_11 - η‖ y(t)‖_1),
    dy_2(t)/dt   =y_2(t)(m_21y_1(t)/y_2(t)+m_22- η‖ y(t)‖_1).

Since the ω-limit set is an invariant subset, we deduce that 

    dy_2(t)/dt= y_2(t)(m_11 - η‖ y(t)‖_1),

that ‖ y(t)‖_1→ m_11/η  and y(t)→z as t→∞. In order to conclude about the solutions of (<ref>) we use <cit.> arguing that z is an asymptotically stable equilibrium of (<ref>) and that for any y(0)∈_+^*, the ω-limit set of any solution y(t) of (<ref>) started at y(0) is not a subset of { 0}×_+. The first claim is easily proved by showing that the Jacobian matrix has negative eigenvalues. For the second claim, let us assume it is not satisfied. Then y_1(t)→ 0 as t→∞. Let us show the contradiction. Let ϵ be sufficiently small and t_0 such that for any t≥ t_0

    0<m_11-ϵ   ≤ m_11+𝒟_11(t)≤ m_11 +ϵ
    
    m_22-ϵ   ≤ m_22 +𝒟_22(t)≤ m_22 +ϵ <0
    
    0<m_12 -ϵ   ≤ m_12+𝒟_12(t)≤ m_12+ϵ.

We introduce 

    P^ϵ(y_1,y_2)   =y_1(m_11-ϵ -η(y_1+y_2))
    
    Q^ϵ(y_1,y_2)   =y_2(m_22+ϵ -η(y_1+y_2))+(m_12+ϵ)y_1

and 

    𝒜^ϵ   ={ y∈_+^2:y_1+y_2 ≤ m_11 -ϵ}
    ℬ^ϵ   ={ y∈_+^2:y_2≥1/2η(m_22+ϵ -η y_1 +√((m_22+ϵ -η y_1)^2 + 4η (m_12 +ϵ)y_1))}.

We deduce that there exists t_1 such that for any t≥ t_1, dy_2(t)/dt<0 on ℬ^ϵ and y_1(t)<ϵ. We deduce that there exists t_2 such that for all t≥ t_2, y_1(t)+y_2(t)≤m_11-ϵ/η. So for all t≥ t_2, dy_1(t)/dt≥ 0 which is absurd. 



 §.§ Differential inclusions

In this appendix, we recall the results of <cit.> which concern the approximation of Markov chains by differential inclusions.

Let ϵ>0 be a scale parameter. Let (Y^ϵ(k),k∈) be a Markov chain with values in ^d. The drift of the Markov chain Y^ϵ is defined by 

    g^ϵ(x)=𝔼[Y^ϵ(k+1)-Y^ϵ(k)| Y^ϵ(k)=x].

Let (γ^ϵ)_ϵ>0 be such that lim_ϵ→ 0γ^ϵ = 0 and let us denote

    f^ϵ(x)=g^ϵ(x)/γ^ϵ.

One can write the evolution of the Markov chain as a stochastic approximation algorithm with constant step size γ^ϵ

    Y^ϵ(k+1)=Y^ϵ(k)+γ^ϵ(f^ϵ(Y^ϵ(k))+U^ϵ(k+1))

where U^ϵ is a martingale difference sequence with respect to the filtration associated with the process Y^ϵ. 

Let us define 

    F(x)=conv({acc_ϵ→ 0f^ϵ(x^ϵ) for all x^ϵ such that lim_ϵ→ 0x^ϵ =x})

where conv(A) denotes the convex hull of the set A and acc_ϵ→ 0f^ϵ(x^ϵ) denotes the set of accumulation points of the sequence f^ϵ(x^ϵ) as ϵ→ 0.
Let us denote by 𝒮_F(T,x^0) the set of solutions (x(t), t∈[0,T]) of the differential inclusion

    dx(t)/dt∈ F(x(t)),   t∈[0,T]
    
    x(0)=x^0.

Let us recall the definition of  a solution of (<ref>).

A map x:[0,T]↦^d is a solution of (<ref>) if there exists a map φ:[0,T]↦^d such that:


  (i)For all t∈[0,T], x(t)=x^0 +∫_0^tφ(s)ds,

  (ii)For almost every t∈[0,T], φ(t)∈ F(x(t)).


In particular (i) is equivalent to saying that x is absolutely continuous. (i) and (ii) imply that x is differentiable at almost every t∈[0,T] with dx(t)/dt∈ F(x(t)).



We define the continuous process Y^ϵ(t) as the piecewise interpolation of Y^ϵ whose time has been accelerated by 1/γ^ϵ: for all k∈ℕ, Y^ϵ(kγ^ϵ)=Y^ϵ(k) and Y^ϵ is linear on [kγ^ϵ, (k+1)γ^ϵ]. We have the following theorem proved in <cit.>.

Assume that:


  •There exists a constant c>0 such that for all y∈^d,  ‖ f^ϵ(y)‖≤ c(1 +‖ y‖).

  •U^ϵ is a martingale difference sequence which is uniformly integrable.

If Y^ϵ(0) tends to x^0 in probability as ϵ tends to zero, then

    inf_y∈𝒮_F(T,x^0)sup_t∈[0,T]‖Y^ϵ(t) - x(t)‖⟶ 0

in probability as ϵ tends to zero.




 §.§ Numerical simulation

We give the Python script for the numerical simulations of the individual based model described in Section 2. The algorithm is based on a classical acceptation/reject method (<cit.>, <cit.>)






#parameters

#intensity of competition
competition = 0.0005  

#probability of mutation
p_mut = 0.05                    

#variance of mutations
var_mut = 0.05					

#number of jumps
number_of_jumps = 1000000           

#initial population size
population_size = 10000		


#traits in the initial population

trait = numpy.ones((population_size, 2))
for k in range(len(trait[:, 0])):
    trait[k, :] = [1.2, 1.6]          


#initial matrix population: [x_b,x_d, living/dead, birth date, death date, id, id parent, parent senescent]


population_0 = numpy.zeros((len(trait[:, 0]), 8), order='C', dtype=numpy.float32)
population_0[:, 0:2] = trait         
population_0[:, 2] = 1                      
population_0[:, 3: 5] = 0  
for l in range(len(population_0[:, 0])):
    population_0[l, 5] = l+1             


#birth rates

def b(x, a):
    if a <= x[0]:
        u = 1       
    else:
        u = 0.0
    return u

#death rates

def d(x, a):
    if a > x[1]:
        u = 1         
    else:
        u = 0.0
    return u


# maximal jump intensity by individual  (for the acceptation/reject method)       

intmax = 2.0 + competition 



#Lansing effect

def lansing_effect(x, u):
    if x[1] > 0:                                                               
        if u < x[1]:                                                           
            r = x                                                              
        else:                                                                  
            r = [x[0], 0]
    else:                                                                      
        r = x
    return r


#mutation kernel

def dm(x):
    g = random.gauss(0, var_mut)
    while x + g < 0:
        g = random.gauss(0, var_mut)
    return x + g


#functions for acceptation/reject method

def acceptation_rejet_clone(t, a, n):
    return b(t, a)*(1-(p_mut*p_mut+2*p_mut*(1-p_mut)))/(intmax*n)

def acceptation_rejet_mutant_1(t, a, n):
    return acceptation_rejet_clone(t, a, n)+(b(t, a)*2*p_mut*(1-p_mut))/(intmax*n)

def acceptation_rejet_mutant_2(t, a, n):
    return acceptation_rejet_mutant_1(t, a, n) + (b(t, a)*p_mut*p_mut)/(intmax*n)

def acceptation_rejet_mort(t, a, n):
    return acceptation_rejet_mutant_2(t, a, n) + (d(t, a) + (n-1)*intensite_competition)/(intmax*n)



#transition of the process: it is based on a classical acceptation/reject method

def transition(p, time):

#initialisation for acceptation/reject method
	
#living individuals	
	viv = p[(p[:, 2] == 1), :]     
    
#population size    
    n = len(viv[:, 0])                                                   
    
#jump time    
    jump_time = random.expovariate(1)/(intmax*n*n)        
    
#uniform law on (0,1)    
    u = random.uniform(0, 1)    
    
#random sampling of one individual   
    ind = random.randint(0, n-1)   
    
    
    w = time - viv[ind, 3] + jump_time #age of this individual
	
#acceptation/reject method        
    while u > acceptation_rejet_mort(viv[ind, :2], w, n):          
        jump_time += (random.expovariate(1)/(intmax*n*n))                
        u = random.uniform(0, 1)                                             
        ind = random.randint(0, n-1)                                         
        w = time - viv[ind, 3] + jump_time
	
#when accepted
    s = viv[ind, :]                                                          
    a1 = acceptation_rejet_clone(s[:2], w, n)                                
    a2 = acceptation_rejet_mutant_1(s[:2], w, n)                             
    a3 = acceptation_rejet_mutant_2(s[:2], w, n)                             

    a = 0                                                                    

    if w > s[1]:                                                             
        a = 1
	
#Lansing effect
    s[:2] = lansing_effect(s[:2], w)

#if clonal birth
    if u <= a1:                                                                                                                                             
        c1 = [1, time + jump_time, 0.0, p[-1, 5]+1, s[5], a]                                                                                      
        p = numpy.vstack((p, numpy.append(numpy.array(s[:2]), numpy.array(c1))))                                                      
	
#if birth with mutation (1 trait)    
    elif a1 < u <= a2:                                                                                                                                      
        x1 = int(numpy.random.randint(1, 3, 1))                                                                                                                  
        z1 = [[dm(s[0]), s[1]], [s[0], dm(s[1])]]                                                       
        c2 = [1, time + jump_time, 0.0, p[-1, 5]+1, s[5], a]                    
        p = numpy.vstack((p, numpy.append(numpy.array(z1[x1-1]), numpy.array(c2))))                                                        
    
#if birth with mutation (2 traits)     
    elif a2 < u <= a3:                                                                                               
        c3 = [1, time + jump_time, 0.0, p[-1, 5]+1, s[5], a]
        p = numpy.vstack((p, numpy.append(numpy.array([dm(s[0]), dm(s[1])]), numpy.array(c3))))
    
#if death    
    else:                                                                                              
        p[(p[:, 5] == s[5]), 2] = 0                                                                      
        p[(p[:, 5] == s[5]), 4] = time + jump_time            

    time += jump_time

    return p, time


#for simulating one trajectory of the process of number_of_jumps jumps

def trajectoire(p):
    population_size = sum(p[:, 2])						time = 0                                       
    for cs in range(0, number_of_jumps):
        if population_size >= 1:
            x = transition(p, time)                
            p = x[0]                                
            time = x[1]                            
            population_size = sum(p[:, 2])
            print cs
        else:
            break
    p[(p[:, 4] == 0), 4] = time                   
    return p, time                               










  
acknowledgements. We acknowledge partial support by the Chaire Modélisation Mathématique et Biodiversité of Veolia Environment - École Polytechnique - Museum National d'Histoire Naturelle - FX. We acknowledge partial support by CNRS and the ATIP/Avenir-Aviesan kick-starting group leaders program. This work was also supported by a public grant as part of the investissement d'avenir project, reference ANR-11-LABX-0056-LMH, LabEx LMH. Finally, we acknowledge the reviewers for their interesting and constructive suggestions.
  




plain 
 

