Horizontal gene transfer: numerical comparison between stochastic and deterministic approaches
Vincent Calvez Susely Figueroa Iglesias Hélène Hivert Sylvie Méléard Anna Melnykova Samuel Nordmann







Horizontal gene transfer: numerical comparison between stochastic and deterministic approaches
    Vincent Calvez[CNRS & Institut Camille Jordan, Lyon.  ], 
Susely Figueroa Iglesias[Institut de Mathématiques de Toulouse; UMR 5219, Université de Toulouse, UPS, IMT, F-31062 Toulouse Cedex 9, France. e-mail: ], 
Hélène Hivert[Univ Lyon, École centrale de Lyon, CNRS UMR 5208, Institut Camille Jordan, F-69134 Écully, France. e-mail: ],  

Sylvie Méléard[CMAP, UMR 7641, Ecole polytechnique, France. ], Anna Melnykova[Université de Cergy-Pontoise, AGM CNRS UMR 8088 and Université de Grenoble Alpes, LJK CNRS UMR 5224, e-mail: ], Samuel Nordmann[CAMS, EHESS, PSL Université, Paris, e-mail: ]
    
==================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================



 Horizontal gene Transfer (HT) denotes the transmission of genetic material between two living organisms, while the vertical transmission refers to a DNA transfer from parents to their offspring. Consistent experimental evidence report that this phenomenon plays an essential role in the evolution of certain bacterias. In particular, HT is believed to be the main instrument of developing the antibiotic resistance. In this work, we consider several models which describe this phenomenon: a stochastic jump process (individual-based) and the deterministic nonlinear integrodifferential equation obtained as a limit for large populations. We also consider a Hamilton-Jacobi equation, obtained as a limit of the deterministic model under the assumption of small mutations. The goal of this paper is to compare these models with the help of numerical simulations. More specifically, our goal is to understand to which extent the Hamilton-Jacobi model reproduces the qualitative behavior of the stochastic model and the phenomenon of evolutionary rescue in particular.




Keywords: Horizontal gene transfer, stochastic individual-based models, integro-differential equations, Hamilton-Jacobi equation, evolution dynamics, resistance to antibiotics.




§ INTRODUCTION


Accurate mathematical description of the evolutionary mechanism is an open question in biology, medicine, and industry. In particular, transmission of pathogens, or antibiotic resistance of bacteria is directly linked to the ability of the bacteria population to mutate and exchange genetic material either vertically (from parents to offspring), or horizontally (from the interaction between non-parental individuals). 

Horizontal Gene Transfer was first described in bacteria when the antibiotic resistance was discovered. This resistance occurs when one bacterial cell becomes resistant to an antibiotic due to mutation, and then transfers resistance genes to other species of bacteria. 
However the Horizontal Transfer of biologic information is not restricted to genes, it also describes the transfer of plasmids and endosymbionts, see for example <cit.>. 
Some artificial applications of horizontal transfer include forms of genetic engineering (Gene Delivery) that result in an organism with its genes changed in some way, and, consequently, possessing new properties or functions (see for instance <cit.>). These applications are particularly useful for "Gene Therapy", which is an experimental procedure that may help treat or prevent genetic disorders and some types of cancer.

The primary goal of our work is to describe the mechanism of the transfer itself and explain how it affects the population dynamics.  Throughout the paper we abbreviate the Horizontal Transfer to HT.

Our study starts with finding a good model of a bacteria population. Several mathematical models for describing a population dynamics were proposed in literature. The first model we consider is a stochastic birth and death process (see, for reference, <cit.>), which describes the dynamics of reproduction, competition, and exchange of genetic material between individuals in a population. The phenotype of each individual is described by a numerical parameter, called trait. Numerical experiments show that the effect of a unilateral horizontal gene transfer may lead to a cyclic behavior of the population. Roughly speaking, while HT drives individuals towards a non-fit phenotype — and, consequently, to extinction, very few not affected by transfer fit individuals may eventually repopulate the environment, before being driven again to deleterious phenotypes. This phenomenon is called an evolutionary rescue of a small population. 

However, within a framework of stochastic jump processes, it is hard to define and study the observed cycling phenomena accurately. The second drawback of the stochastic system is that it is costly to compute, especially for a large time scale and population size. Thus, in the case of a large population, it is more practical to work with a deterministic PDE model, describing  the limiting behaviour of a stochastic system when the population size goes to infinity <cit.>.
In certain settings, the population dynamics involve concentration phenomena (i.e., the convergence of the population density to singular solutions, such as Dirac masses). In that case, the PDE formulation is not suitable. Thus, applying a limiting procedure for small mutations and time rescaling to the PDE model, we pass to a Hamilton-Jacobi type equation. 

The primary goal of our work is thus to conduct a numerical analysis of the population dynamics on a macroscopic individual-based model and to compare it with the deterministic system which is obtained as a limit for a large population. We are especially interested in determining to which extent the limiting Hamilton-Jacobi equation can grasp qualitative properties of the stochastic model.
This framework has already been successfully used to understand the concentration phenomena, and the location of the dominant trait (see for instance <cit.>). We aim to understand if the Hamilton-Jacobi approach is also well suited to describe the evolutionary rescue phenomena which crucially rely on an accurate description of the small populations.

On this step, the choice of an approximation scheme for simulating solutions of the PDE model is of tremendous importance. As we further explain in Section <ref>, classical explicit schemes do not preserve the asymptotic behavior of the solution if the time rescaling step goes to 0. From a numerical point of view, it involves operations with exponentially big values, which lead to non-negligible errors for explicit numerical schemes. We address this question by proposing an asymptotic preserving scheme for a Hamilton-Jacobi equation, adapting an approach proposed in <cit.>. More generally, the numerical approximation problem for solutions of Hamilton-Jacobi equations is treated in <cit.>.

This paper is structured as follows: in Section <ref> we introduce the model both in a stochastic and deterministic setting, and formally derive the limiting Hamilton-Jacobi equation. Then, we simulate a jump process, describing the bacteria population, and study its properties for different values of parameters. Numerical experiments are gathered in Section <ref>. We aim to numerically determine the critical HT rate, which leads to an almost sure extinction of the whole population. On the next step, we conduct the same analysis for a Hamilton-Jacobi equation with the help of an asymptotic preserving scheme and compare it with the stochastic model on an appropriate timescale, and explain why the classical scheme fails to work. We end our study with conclusions and discussion of yet unsolved numerical and theoretical questions.



§ MODEL



 §.§ Stochastic model

We consider a stochastic model describing the evolution of a population structured by phenotype. In a general case it is described at each time t by the point measure

    ν^K_t(dx) = 1/K∑_i=0^N^K_tδ_X_i(t)(dx),

where parameter K is a scaling parameter, referred to as the carrying capacity. It stands for the maximal number of individuals that the underlying environment is able to host (K can represent, for example, the amount of available resources). N^K_t = K∫ν^K_t(dx) is the size of the population at time t,  and X_i(t)∈R^n is the trait of i-th individual living at t, which summarizes the phenotype information. In this work we assume n=1, that is, the trait is given by a point on a real line. 

The demography of the population is regulated, first of all, by its birth and death rates. An individual with a trait x gives birth to a new individual with rate b(x).  The trait y of the offspring is chosen from a probability distribution m(x-y)dy (by that we mean that  ∫_R m(x-y)dy = 1). We will refer to it as the mutation kernel. 
An individual with a trait x dies according to an intrinsic death rate d(x) plus an additional death rate CN_t^KK (independent of x) which stands for the competition between individuals. 

Finally, an individual with a trait x can induce a unilateral HT to an individual with trait y at rate h_K(x,y,ν), such that the pair (x,y) becomes (x,x).  In literature this kind of transfer is sometimes referred to as a conjugation.  For simplicity, we assume h_K(x,y,ν) to be in the particular form 

    h_K(x,y,ν) = h_K(x-y,N) =τ_0α(x-y)/N/K,

where N=K∫_Rν(dx) is the number of individuals, τ_0>0 is a constant and α is either a  Heaviside, or a smooth bounded function, such that for a small δ>0:

    α(z)=
    {   0       if z<-δ
       1       if z>+δ.,
      α'(0)=1/2δ,

where δ is the stiffness parameter. We introduce δ to have the advantage of working with a smooth function (which will be useful in the following parts), while mimicking the binary nature of the Heaviside function. 

For a population ν=1/K∑_i=1^Nδ_x_i and a generic measurable bounded function F, the generator of the process is then given by: 

    L^K F(ν) = ∑_i=1^N b(x_i)∫_R(F(ν+1/Kδ_y) - F(ν))m(x_i, dy)
    
    +∑_i=1^N (d(x_i) + CN/K)(F(ν-1/Kδ_x_i) - F(ν)) 
     
    +∑_i,j=1^N h_K(x_i, x_j, ν)( F(ν + 1/Kδ_x_i-1/Kδ_x_j) - F(ν) ).


It is standard to construct the measure-valued process ν^K as the solution of a stochastic differential equation driven by Poisson point measures and to derive moment and martingale properties (see for instance <cit.>). 




 §.§ The PDE model

It is proven (see in particular <cit.>) that for K→+∞ the stochastic process defined by a sequence of point measures given by (<ref>) converges  in probability to the unique solution of a non-linear integro-differential equation. This equation is given by:

    {[                                        ∂_t f(t,x)                                                 =   -(d(x)+Cρ_1(t))f(t,x)+ ∫_^n m(x-y)b(y)f(t,y)dy+;                                                                                                     f(t,x)∫_^nτ(x-y)f(t,y)/ρ_1(t)dy,    (t,x) ∈_+×^n,;                                            ρ_1(t)                                                 =                                     ∫_R f(t,x)dx,;                                            f(0,x)                                                 =                                         f^0(x)>0, ].

where f(t,x) is the macroscopic density of the population with trait x at time t and, accordingly to the previous section, b(x), d(x) and C are the birth, death and competition rate respectively, m is the mutation kernel, and 

    τ(y-x):=τ_0[α(x-y)-α(y-x)]
 
is the horizontal transfer flux.


Now our goal is to pass from micro- to a macroscopic scale with the help of a rescaling. 
On the one hand, we consider the case of small mutations: for a small parameter >0 we define

    m_(x-y)=1/^nm(x-y/).

With a change of variable z=x-y/ we can rewrite the mutation term at (t,x) as

    ∫_^n m_(x-y)b(y)f(t,y)dy= ∫_^n m(z)b(x+ z)f(t,x+ z)dz.

On the other hand, when  is small, the effect of mutations can only be observed in a larger time scale. Thus, we rescale time with t↦t/.

We end up with the following system, for >0, and (t,x)∈_+×^n:

    {[                                          ∂_t f_(t,x)                                                    = -(d(x)+Cρ_(t))f_(t,x)+ ∫_^n m(z)b(x+ z)f_(t,x+ z)dz+;                                                                                                                              f_(t,x)∫_^nτ(x-y)f_(t,y)/ρ_(t)dy,;                                                ρ_(t)                                                    =                                       ∫_R f_(t,x)dx,;                                              f_(0,x)                                                    =                                           f_^0(x)>0. ].




 §.§ The Hamilton-Jacobi limit

We now derive the limiting problem (<ref>) when →0. As we will see, the limiting problem allows us to give a rigorous mathematical framework and to perform useful formal calculations.

Equations in the form of (<ref>) often give rise to a concentration phenomenon, i.e the convergence of f_ towards a Dirac mass when →0 (see <cit.>). The usual way to deal with these asymptotics is to perform a Hopf-Cole transformation (or WKB ansatz), i.e to consider 

    u_(t,x):= ln(f_(t,x)).

This change of variable comes from the intuition that a Dirac mass is no more than a narrow Gaussian, and more precisely that f_ should behave like a Gaussian of variance  when →0.
Accordingly, we expect u_ to have a non singular limit when →0. Incidentally, this substitution also gives insights on the convenient scheme to use for numerical simulations,  as we will see in the following section.

Now our goal is to identify and derive the asymptotic properties of u_ when →0, which will be used for discussions in the sequel. The following computations are only formal, since rigorous proofs are often intricate in this context.
Substituting (<ref>) into (<ref>) we deduce that u_ satisfies 

    ∂_t u_ = -(d(x)+Cρ_(t))+∫_^n m(z)b(x+
     z)exp{u_(t,x+ z)-u_(t,x)/}dz
     +∫_^nτ(x-y)f_(t,y)/ρ_(t)dy.


Formally, at the limit →0, u_ converges to a continuous function u which satisfies the following Hamilton-Jacobi equation in the "viscosity" sense:

    ∂_t u = -(d(x)+Cρ(t))+b(x)∫_^n m(z)e^z·∇_x udz+τ(x-x(t)),

where 
ρ(t)≥0 is the weak limit of ρ_(t) and 

    x̅(t)=argmax u(t,·).
 
We formally assume here and in the following that the definition of x̅(t) is unambiguous, i.e that u reaches its maximum in a single point.
Note that the limiting function u is not expected to be C^1 for all time. We thus need to deal with a generalized notion of solutions, namely viscosity solution (see <cit.>).

This framework is convenient because most of the information is contained in the dynamics of x̅(t). See the next section for further analysis.



 §.§ Formal analysis on the Hamilton-Jacobi equation

Hamilton-Jacobi equations are particularly known in mathematical biology to be a good model to describe how a population concentrates around the dominant trait(s) when the mutations are small. However, here we are interested to use this model to describe a phenomenon of evolutionary rescue. In this subsection we attempt an analysis of the equation (<ref>). We point out that the calculations are only formal, since rigorous proofs are intricate and beyond the scope of this paper.



  §.§.§ Generality

From an integration of (<ref>) with respect to x and classical computations (under the assumptions of bounded functions for the birth, death and transfer rates), we deduce that our model satisfies a saturation property, i.e. ρ_(t) is bounded from above, uniformly in t≥0 and >0.
From this and ρ_(t)=∫_^ne^u_(t,x)/dx, we deduce that
∀ t>0, sup_x∈^n u(t,x)≤ 0 and the following constraint holds:

    sup_x∈^n u(t,x)=0  when  ρ(t)>0.

Note that our model allows the population to get extinct, thus we cannot expect ρ to be always positive. 
As a byproduct, we derive the concentration property, i.e the formal weak convergence of measures

    f_(t,x)⇀ρ(t)δ_x̅(t)(dx),  when → 0,

where δ_x̅(t)  denotes, as usually, the Dirac measure centered in x̅(t).
From (<ref>) it is possible to formally derive a formula for ρ. Indeed, either ρ(t)=0 or ρ(t)>0 and

    ∂_t u(t,x̅(t))=0,

which implies

    ρ(t)=b(x̅(t))-d(x̅(t))+τ(0)/C=b(x̅(t))-d(x̅(t))/C,

for τ defined in (<ref>).

Having above definitions in hand, we can now perform a formal analysis on the dynamics of x̅(t), defined below in (<ref>). Our aim is to show how the behaviour of the system can be analyzed within the framework of a Hamilton-Jacobi equation (<ref>). To fix ideas, we fix all constants but τ_0 and we assume (<ref>)-(<ref>) as follows:

    b(x)=b_r>0, 
    
    d(x)=d_r x^2,   d_r>0, 
    
    m(z)=1/√(2π)σe^-z^2/2σ^2,

and the transfer function h_K(x,y,ν) is defined in (<ref>).
Moreover we work under the following assumptions:



  (A1) u(t,·) reaches its maximum on a single point x̅(t),

  (A2) x̅(t) is a non-degenerate maximum, i.e ∇_x^2u(t,x)<0,

  (A3) x̅(t) is smooth with respect to t.



    u(t,·)  reaches its maximum on a single point x̅(t),
       x̅(t)  is a non-degenerate maximum, i.e ∇_x^2u(t,x)<0,
       x̅(t)  is smooth with respect to t.

Finally we assume that the initial condition f^0 is a given function of x which reads:

    f_^0(x) = 1/√()e^-x^2/2.









  §.§.§ Smooth dynamics x̅(t).


The following statement deals with the smooth dynamics of x̅(t), i.e in the regime where no jump occurs in the dynamics of x̅(t). 

Under assumptions (<ref>)-(<ref>), the function t↦x̅(t) is an increasing function which satisfies the following inequality ∀ t≥ 0:

    0≤x̅(t)≤τ_0/2dδ.

More precisely, x̅(t) satisfies the canonical equation 
    d/dtx̅(t)= [-∇^2_x u(t,x̅(t))]^-1·(∇_x r(x̅(t))+∇_xτ(0)),

where

    r(x):=b(x)-d(x),

and ∇^2_x u denotes the Hessian of u with respect to the x variable.

 Under the above assumptions we can derive the dynamics of x̅(t), referred to as the canonical equation in the literature (see for instance <cit.>). Indeed, starting from

    ∇_x u (t,x̅(t))=0,

a differentiation with respect to t gives (<ref>).
Equation (<ref>) has a unique singular point x_⋆, which satisfies r'(x_⋆)+τ'(0)=0, with τ defined in (<ref>) and r in (<ref>). 
We find

    x_⋆= τ_0/2d_rδ.
 Note that t↦x̅(t) is increasing when x̅(t)<x_⋆ and decreasing when x̅(t)>x_⋆. Besides, from the initial condition (<ref>), we have x̅(0)=0,
and consequently 0≤x̅(t)≤ x_⋆  ∀ t. 



  §.§.§ Evolutionary rescue.


In general, the canonical equation (<ref>) does not hold in every point of time. Indeed, a new maximum of u can arise in a finite time, which would cause a "jump" in the dynamics of x̅(t): this is what we call an evolutionary rescue. Formally, this is what happens (periodically in time) in the case of cycles, see Figure <ref>. We thus expect x̅(t) to possibly jump periodically, and to follow (<ref>) between two jumps. We now try to characterize the possible jumps.
For T>0, we denote 

    x̅(T^-):=lim_t→ T
     t<Tx̅(t),  x̅(T^+):=lim_t→ T
     t>Tx̅(t).



We assume that (<ref>)-(<ref>) hold until a time T>0, such that u(T,·) reaches its maximum on x̅(T^-) and on another point x̃. Then x̃=0 and x̅(t) will jump towards 0 at time T, i.e x̅(T^+)=0. 


From assumption (<ref>), we have  ∀ t∈[0,T] that u(t,·) is concave non-degenerate on [x̅(t)±θ], with θ>0. For simplicity, we further assume δ≤θ, where δ is defined in (<ref>). 

First, let us show that x̃=0. We define the fitness function of trait x in a population concentrated in x̅:

    F_x̅(x):= r(x)+τ(x-x̅),

where r and τ are respectively defined in (<ref>) and (<ref>).
Note that we have ∂_t u(t,x)=F_x̅(t)(x)-Cρ(t), for t<T. 
But x̃∉[x̅(t)±δ] and the choice of parameters (<ref>)-(<ref>)-(<ref>) implies x̃ must maximize F_x̅(T^-)(·), hence x̃=0.

The second step is to prove that there will be an actual jump towards 0, i.e x̅(T^+)=0. First, note that there exists a small η>0 such that ∀ t∈(T-η,T), u(t,x̅(t))=0 and u(t,0)<0. Let us fix t∈(T-η,T). 
We have F_x̅(t)(0)≥ F_x̅(t)(x̅(t)), and we claim that the inequality is strict. Indeed, since t↦ x(t) is increasing, F_x̅(t)(x̅(t)) is decreasing, whereas F_x̅(t)(0) is constant (as long as η is small enough such that x̅(T-η)>δ). We end up with 

    F_x̅(t)(0)> F_x̅(t)(x̅(t)).

The above inequality expresses the fact that 0 is fitter than x̅(t) in a population with trait x̅(t). In general, this does not allow to conclude that 0 will invade and become the new dominant trait (i.e., that the jump will occur) because it does not imply that 0 will remain fitter during all the process of invasion.
But the particular form of our problem, especially the fact that τ is an odd function, implies 

    F_0(0)> F_0(x̅(t)).

Indeed we have from the definition of F_x̅(x) that

    F_0(0)- F_0(x̅(t))=r(0)-r(x̅(t))+τ(x̅-x̅)-τ(0)=d_r x̅(t)^2>0.

Consequently that for all λ∈[0,1]

    λ F_0(0)+(1-λ)F_x̅(t)(0)>λ F_0(x̅(t))+(1-λ)F_x̅(t)(x̅(t)).

It shows that 0 remains the fittest trait during all the process of invasion, and therefore that 0 will actually invade, i.e that x̅(t) will actually jump towards 0 at time T^+. 





  §.§.§ Threshold for cycles

In the previous section, we described the possible evolutionary rescue, i.e the possible jumps in the dynamics of x̅(t) towards x=0. 
When a jump occurs, a new cycle begins: it leads to a periodical behavior of x̅(t), hence the cycling phenomenon.

We recall that a jump corresponds to a rescue of the population concentrated at x̅(t) by the small population with trait x=0. It is possible only if x̅(t)>δ and if 0 is fitter than x̅(t) during a sufficiently large interval of time (which is the time needed for the small population at x=0 to regrow). Note that 0 is fitter than x̅(t) if and only if

    F_x̅(t)(0)≥ F_x̅(t)(x̅(t))     iff   b_r-τ_0≥ b_r-d_rx̅(t)^2,
       iff  x̅(t)≥ x_resc:=√(τ_0/d_r).

But if no jump occurs, x̅(t) formally follows (<ref>), thus x̅(t)< x_⋆ and x̅(t) converges to x_⋆ when t→+∞ (with x_⋆ is defined in (<ref>)). 

Under assumptions (<ref>)-(<ref>),
 the evolutionary rescue phenomena occurs if and only if 

    τ_0>τ_cyc:= 4d_rδ^2.


 Note that the condition τ_0>τ_cyc is equivalent to x_resc < x_⋆ , which are defined respectively in (<ref>) and (<ref>).



  §.§.§ Threshold for extinction.

The population is said to be "extinct" at time t if ρ(t)=0. 

According to (<ref>), we define x_ext as to solve r(x_ext)=0, i.e

    x_ext:= √(%s/%s)b_rd_r,

that is, a population concentrated at trait x̅ is extinct iff x̅≥ x_ext.

The picture is simple in the case of stabilization without cycles, i.e when τ_0≤τ_cyc (see (<ref>)). In this case, we recall that x̅(t) formally follows (<ref>) for all t>0, thus x̅(t)< x_⋆ and x̅(t) converges to x_⋆ when t→+∞ (where x_⋆ is defined in (<ref>)). Thus, if x_⋆≤ x_ext, we have ρ(t)>0 for all t>0; on the contrary, if x_⋆> x_ext, there exists a time t_ext>0 for which ρ(t)=0 for all t≥ t_ext. It gives a sharp threshold for extinction of the population: indeed, the population eventually gets extinct if and only if x_⋆>x_ext, which naturally leads us to the following statement. 
 Under assumptions (<ref>)-(<ref>), if τ_0≤τ_cyc, then 
 the population eventually gets extinct if and only if 
 
    τ_0> τ_ext:=2√(b_rd_r)δ.



We point out that, surprisingly enough, τ_ext is an increasing function of the death rate d_r, meaning that under a higher death rate, the population can survive to a higher HT rate. The interpretation we propose is that if d_r is high, the population driven outward x=0 dies rapidly, thus the  population that remained closer to 0 undergoes a milder HT, which makes the overall population more resistant to a high HT rate. 

Let us now focus on the case where the cycling phenomenon occurs, i.e when τ_0>τ_cyc. In this case, x̅(t) will follow (<ref>) and will periodically jump to x=0. First, note that if x_⋆<x_ext, x̅(t) remains below x_ext for all t and the population does not get extinct:

    if  τ≤τ_ext,  then  ρ(t)>0,  ∀ t>0.

The most intricate case is when x_⋆>x_ext, which contains cases of extinction and non-extinction, depending on whether the jump of x̅(t) towards 0 happens before or after x̅(t) has passed beyond x_ext. In other words, extinction can be avoided if the evolutionary rescue happens before the dominant trait is led to extinction, i.e  if x̅(T^-)≤ x_ext, where T is the time where the jump of x̅(t) towards 0 occurs. However, we are not able to give a satisfactory formula or estimate on T. 

Besides, when the jump of x̅(t) occurs, it can happen that the trait x=0 is not fit enough to avoid extinction:  in this case the evolutionary rescue does not manage to sustain the population. It corresponds to the case x_resc>x_ext. We have the following threshold:
the evolutionary rescue is able to sustain the population iff r(0)+τ_0>0, which is equivalent to

    τ_0<τ_sus:=b_r.


If τ≥τ_sus, the population eventually gets extinct. If τ<τ_sus, the population is effectively rescued by the evolutionary rescue, even in the case where it passed through an episode of extinction during the previous cycle: in some cases the population is able to regrow after being extinct, which can be seen on Figure <ref>. We think this is an interesting feature that the Hamilton-Jacobi approach is able to grasp. Regarding the stochastic model, an episode of extinction on Hamilton-Jacobi corresponds to an interval of time where the population reaches extremely small values (of order e^-1/, with  the variance of the mutation kernel), and the probability that every individual dies is bigger than the survival of the population.  

Assume (<ref>)-(<ref>) and τ_0>τ_cyc.


  * if τ_0≤τ_ext, the population never gets extinct.

  * the evolutionary rescue effectively manages to sustain the population if and only if τ_0 <τ_sus:=b_r. 





  §.§.§ Characteristics of a Hamilton-Jacobi equation

Denoting

    -H(t,x,p):= -(d(x)+Cρ(t))+b(x)∫_m(z)e^pzdz+τ(x-x̅(t)),

from (<ref>) we have ∂_t u(t,x)+H(t,x,∇_x u(t,x))=0. Since H is convex in the p variable, we have the following representation formula (see <cit.>).

    u(t,x)=inf_γ∈ C^0(_+,)
    γ(t)=x[∫_0^t L(s,γ(s),γ̇(s))ds+u^0(γ(0))],

where L(t,x,v) is the Lagrangian of the equation, obtained through a Legendre transform (or a convex conjugate) of H.

Every γ which is admissible as a minimizer in (<ref>) is called a characteristic of the Hamilton-Jacobi equation (<ref>). 
Note that every characteristic γ formally satisfies the condition

    d/ds[∂_v L(s,γ,γ̇(s))]=∂_x L(s,γ(s),γ̇(s)).

 (<ref>) holds because γ is a critical point of the functional defined in (<ref>). 
Note that if we replace H by H̃(x,p)= -x^2/2+p^2/2+1,
the Legendre transform of H̃ can be computed explicitly:

    L̃(x,v)=x^2/2+v^2/2-1.

Then (<ref>) becomes

    γ̈(s)=γ(s).




§ NUMERICAL TESTS

In this section we perform several numerical tests for the presented models considering different values of parameters, replicating different scenarios: stabilization around an optimal value, cycles (occurring through the evolutionary rescue phenomena) and the extinction.We then compare the numerical results obtained for the stochastic and deterministic approaches, using in particular an asymptotic-preserving scheme which allows us to observe the population dynamics on the passage from the integro-differential equation (<ref>) to a limit (<ref>).
Throughout this section we define the birth, death rates and the mutation kernel to those given in (<ref>)-(<ref>) respectively, with the parameters fixed throughout all the experiments to b ≡ 1, d_r≡ 1, C≡ 0.5 respectively (unless otherwise stated). 



 §.§ Stochastic model



  §.§.§ The scheme


Our aim is to simulate the population dynamics over a fixed interval [0,T]. 

We begin by simulating an initial population of size N^0. We assume that the population is normally distributed around the mean trait x^0_mean with a standard deviation σ^0 so that the resulting vector X^0∈R^N^0. We know that in a time step Δ, an individual can die, give birth, or be a subject to HT. Each event happens according to a certain probability that we compute from the rates. More detailed description of the simulations is provided in Algorithm 1.

Note that in our setting it is possible that 1,2  or  3 events happen within the same time step. Keeping a discretization time step small helps us to keep a biological sense in our simulation: even if the event of horizontal transfer with an "already dead" individual is possible in our setting (if  T_d≤ T_HT≤Δ), this event is extremely rare. 



We simulate the population of initial size N_0 = 10000 up to time T=1000 with Δ = 0.01, with the parameters being defined at the beginning of the section, and α is a Heaviside function. Even if a Heaviside function is not the most easy to analyze when we pass to the deterministic limit of the system (see Subsections <ref> and <ref>), we use it for the stochastic simulation, since it is the most straightforward model for HT in biological context, and is much faster to compute than a smooth function.  We fix all constants but τ_0, which regulates the Horizontal Transfer, and study how it affects the dynamics.
Then we plot the density of the population at each moment of time (left side of each Figure): brighter colors on plot mean that there is a big amount of individuals with very similar traits. On the right top and right bottom we plot the normalized population size (ratio between the actual size and the carrying capacity of the system), and the mean trait.  
 
Depending on the parameters we may observe three types of behavior (see Figure <ref>). First possibility, for small values of τ_0, is the stabilization (Figure <ref>).
In this case the population rapidly reaches the equilibrium and concentrates around the optimal trait, which is close to 0.1

(with stochastic fluctuations). Note that in this case, the mean trait is shifted in comparison to the optimal trait without HT (which is x=0).

Second option, for intermediate values of τ_0, is the cycling behavior (Figure <ref>). Since the transfer rate is sufficiently large, the population is driven towards a deleterious trait, which is eventually less fit than the trait x=0. 
If the drift is not too strong, the very few individuals which were not affected by HT and remained fit (with x close to 0) manage to regrow and eventually repopulate the environment, which launches the cycle again. 

The last possibility, for large values of the horizontal transfer rate τ_0, is the extinction of the population (Figure <ref>). It occurs because too many individuals were affected by deleterious traits of their neighbors, so that they die faster than is needed for replicating the population. 



To understand better this phenomenon, we have to give a precise definition of what do we actually refer to, when we say "the critical value" of the transfer rate?  In stochastic setting the answer is not trivial, and that is where the individual-based model reaches its limit. What we observe experimentally is the following when we change the value of HT rate starting from zero, the cycles in the population dynamics become more clearly visible, the fluctuations of the mean trait and the population size become more ample, until at some point the probability of extinction overweights the probability of survival and, finally, at the value of τ_0, which we call "critical" we obtain an almost sure extinction.

But since we are working with a point process, giving a strict definition of a "critical value for an extinction" in terms of probability measures seems to be out of reach. Even in the experimental setting this notion is ambiguous: when the value of τ_0 is getting closer to a "critical" (numerically we observe an almost sure extinction at τ_0 = 0.49), in different repetitions of the same experiment we may observe 
different types of behavior: either cycles, or extinction, which occurs after several cycles. It is illustrated on Figure <ref>, where the computations, launched with exactly the same set of parameters, give very different results. Furthermore, it is not always clear how to differentiate between the stabilization and cycles, especially when the variance of the mutation kernel is large. To the best of our knowledge, there is no straightforward way to analytically measure the probability of each outcome under given initial conditions, which makes the model difficult to analyse.

This constraint of an individual-based model naturally leads us to studying a limiting system described in Subsection <ref>. 






  §.§.§ Lineages

With the help of the stochastic model we can keep track of the lineage of an individual i which lives at a final observed time T. More precisely, we are interested in a history of a phenotype which leads to a long-term survival of an individual. 


We illustrate some numerical experiments on Figure <ref>. The four simulations are done with the same parameters. In the background, every point with coordinates (t,x) represents an individual with trait x living at time t (as in Figure 1). The solid lines represent the lineages of the individuals that live at final time. Small fluctuations are the results of birth with mutation, while the large upwards jumps correspond to an occurrence of a HT.

First of all, we can see on the plot that all the lineages are gathered into one line up to t=400. It means that all individuals that live at final time t=700 emanate from one single ancestor of the initial population. This phenomenon is well known and referred to as coalescence in the literature (see for instance <cit.>, 
or <cit.>  for a mathematical description of a classical population genetics theory).

Besides, we see that the lineages remain centered around x=0 during almost all the observed time. It is explained by the fact that every lineage that goes to a high value of x (corresponding to deleterious phenotype) cannot recover (since the mutations are small), and eventually goes extinct. This illustrates that the population manage to sustain because of the very few individuals that were not affected by HT throughout the history.





 §.§ 
Numerical scheme for the PDE model



In this subsection,  a numerical scheme for (<ref>) is presented, and its properties are numerically investigated. 
For the discretization of (<ref>), we consider a bounded space of traits [X_min, X_max], discretized with N_x points. Denoting N_x the number of discretization points of the interval [X_min, X_max], we define 

    = X_min-X_max/N_x-1,

and 

    x_i=X_min+i,   0≤ i≤ N_x-1.

We consider the time interval [0,T_max], discretized with N_t points t_n=n, for 0≤ n≤ N_t-1, and where  is defined as

    =T_max/N_t-1.

The approximations of the solution f of (<ref>) at (t_n,x_i), and of its density ρ at t_n are denoted f^n_i and ρ^n respectively. 
We recall that the initial condition f^0 is a smooth function of x given in (<ref>) and the initial density ρ^0 is computed using a left-point quadrature rule for f^0 as follows: 

    ρ^0=∑_i=0^N_x-1 f^0(x_i).

The scheme is written with an explicit Euler scheme, in which the integrals are computed with a left-point quadrature rule. For n≥ 1 and 0≤ i≤ N_x-1, it reads

    εf^n+1_i-f^n_i/= 
    (d(x_i)+Cρ^n)f^n_i  
    + [ m*(bf) ]^n_i + f^n_i ∑_j=0^N_x-1τ(x_i-x_j) f^n_j/ρ^n.


In (<ref>), the convolution product [ m*(bf) ]^n_i is computed with a left-point quadrature rule, as well of the other integrals.
To do so, a grid in the z variable is defined as for the x variable. Let Z_min and Z_max, and the number N_z of discretization  points be given. The grid in z is defined as 

    ∀ 0≤ k≤ N_z-1, z_k=Z_min+k,

where =(Z_max-Z_min)/(N_z-1). When x_i+ε z_k ∈[X_min,X_max],  the value of f(t_n,x_i+ε z_k) is approximated by linear interpolation of the (f^n_i)_0≤ i ≤ N_x-1. When x_i+ε z_k< X_min, or x_i+ε z_k> X_max, it is computed with a linear extrapolation of the (f^n_i)_0≤ i≤ N_x-1, using the slope at the corresponding end of the X domain. Using the notation f^n(x_i+ε z_k) for the approximation of f(t_n,x_i+ε z_k), we then define


    [ m*(bf) ]^n_i= ∑_k=0^N_z-1 m(z_k) b(x_i+ε z_k) f^n(x_i+ε z_k).




  §.§.§ Case = 1: comparison with stochastic model


First thing that we are interested in is whether under identical parameters and initial conditions we may reproduce the same
behavior as in the stochastic model. 

Thus, we conduct several experiments, fixing parameter  to 1 (thus, we do not rescale time, nor mutation rate), leaving all the other parameters fixed to the same values as in the stochastic simulation case.



As we may see on Figure <ref>, simulations in overall correspond to those of the stochastic model. Indeed, when the HT rate τ_0 is small enough the population rapidly stabilizes around its equilibrium state (see  Figure <ref>), as in the stochastic simulations. Further similarity between two models is that in both cases the optimal trait is shifted a bit above 0. It is caused by the HT phenomenon.  

For larger values of τ_0, where we would expect to have distinguishable cycles, we observe indeed damped oscillations, see Figure <ref>. We stress out that for the stochastic model it is not the case, see Figure <ref>. The way we understand the damping in the oscillations is that the PDE model and the numerical algorithm that we use are not designed to have a precise grasp on the exponential small values of f, on which the cycling phenomenon relies. This limitation suggests to perform the change of variable (<ref>), and to write a numerical scheme which converges uniformly when → 0. This is what the next subsection is devoted to. 

On Figure <ref>, we observe that as τ_0 becomes larger the population gets extinct, and then, surprisingly enough, "reborns" after a period of extinction. This scenario can only be reproduced on density-based models, since in individual-based model any extinction is definitive. On Figure <ref> we observe a full extinction of the population without regrowth. We will give further insights on those two cases in the next subsection. 



 §.§ The scheme for the Hamilton-Jacobi equation



  §.§.§ Case → 0: description of the numerical scheme


As the rescaling parameter  goes to 0, the model given by (<ref>) gets closer to its limiting state (<ref>). However, numerical approximation of the (<ref>) for ≪ 1 is not a trivial task. 
 Indeed, for small ε, the solution f^ of (<ref>), is expected to concentrate around the dominant trait. To be able to catch its stiffness numerically, one has to refine the grid in x, to ensure enough precision in the computation of f. As a consequence, the computational cost of the numerical simulations increases when → 0, and reaching the asymptotic regime with this scheme is not possible. In this part, we present a numerical scheme for (<ref>) which enjoys stability properties in the limit ε→ 0.


To avoid the increase of computational cost when reaching the asymptotics, and to ensure the scheme approaches the limit Hamilton-Jacobi equation for small , a scheme for the solution u^ of (<ref>) which enjoys the Asymptotic Preserving (AP) property is proposed here. Such schemes have been introduced in <cit.>, their properties are often summarized by the following diagram: 

    [    P_         P_0;    90          90;  S_^h       S_0^h ]

It should be understood as follows: when the parameter >0 is fixed, the scheme S^h_ is consistent with the -dependent problem P_. When  goes to 0, the solution of P_ converges to the solution of the limit problem P_0. The AP scheme S^h_ is stable along the transition to the asymptotic regime. It means that, when  goes to 0 with fixed discretization parameters h, the scheme becomes a limit scheme S^h_0, which is consistent with the limit problem P_0.

As an AP scheme is required to enjoy stability properties when  is going to 0, one has to ensure that all the quantities that have to be computed enjoy this property. In the case we are considering, the main concerns are the computation of the integral containing the birth term, the computation of the integral containing the transfer term and the computation of ρ. If all of them are correctly defined, the scheme we propose reads 

    u^n+1_i-u^n_i/=-(d(x_i)+Cρ^n+1) 
    + B^n_i + T^n_i,

where B^n_i stands for an approximation of 

    ∫_ m(z) b(x_i+ z) e^(u_(t^n,x_i+ z)-u_(t^n,x_i))/dz,

and T^n_i is for 

    ∫_τ(x_i-y) f(t^n,y)/ρ(t^n)dy.

Here, we used the notations and discretization grids defined in the beginning of Section <ref>, and the dependencies in  are omitted to simplify the notations.
In what follows, we present how T^n_i, B^n_i and ρ^n+1 can be computed in a way that ensures they are consistent with their definition for fixed , that they can be computed with a constant computational cost with respect to , and that their asymptotic behavior when  goes to 0 is meeting the continuous one (<ref>).


  * Computation of T^n_i. The direct approximation of (<ref>) with a quadrature rule is consistent for ∼ 1. However, since f is expected to concentrate when → 0, it lacks precision in the asymptotic regime. Especially, the convergence of f/ρ to a Dirac is not ensured when the integral is approximated directly. Remarking that 

    f^(t^n,y)/ρ^(t^n)=e^u^(t^n,y)//∫_e^u^(t^n,z)/dz 
    = e^(u^(t^n,y)-max_x u^(t^n,x))//∫_e^(u^(t^n,z)-max_x u^(t^n,x))/dz ,

(<ref>) is computed with a left-point quadrature rule in the integrals of the previous expression. It reads


    T_i^n=∑_j=1^N_x-1τ(x_i-y_j) e^(u^n_j-max_l u^n_l)//∑_k=0^N_x-1e^(u^n_k-max_l u^n_l)/
    =∑_j=1^N_x-1τ(x_i-x_j) e^(u^n_j-max_l u^n_l)//∑_k=0^N_x-1e^(u^n_k-max_l u^n_l)/.

For fixed , (<ref>) is consistent with (<ref>). Since all the arguments of the exponentials are nonpositive, the limit of (<ref>) for small  can be read on that expression. Denoting j_0 the index such that 

    u^n_j_0=max_l u^n_l,

and supposing that there exists a unique such j_0, the limit of (<ref>) for small  is

    τ(x_i-x_j_0).

This is consistent with the last term in the limit Hamilton-Jacobi equation (<ref>).



  * Computation of B^n_i. Once again, the numerical approximation of (<ref>) is done with a 

quadrature
in the integral. Using the notations of  Section <ref>, a grid in z is defined. The functions m and b are respectively evaluated at z_k and x_i+ z_k, but the interpolation of u^n at x_i+ z_k has to be done with special care to make the scheme enjoy the expected asymptotic behavior. 

Using a left-point quadrature rule, (<ref>) is approximated by 

    ∑_k=0  |z_k| ≤ dx ^N_z-1
    m(z_k) b(x_i+ z_k) e^
    z_k∇^,small_n,i,k
    
    + 
    
    ∑_k=0  |z_k| > dx ^N_z-1
    m(z_k) b(x_i+ z_k) e^
    z_k∇^,large_n,i,k
    ,

where ∇_n,i,k^ stands for an approximation of 

    u^(t^n,x_i+ z_k)-u^(t^n,x_i)/ z_k.

In both cases, it is computed with a linear interpolation of the values u^n_i. Hence, ∇_n,i,k^, large is given by 

    ∇_n,i,k^, large=ũ^n_i,k-u^n_i/ z_k,

where ũ^n_i,k is computed as the linear interpolation of (u^n_i)_1≤ i≤ N_x at x_i+ z_k. If x_i+ z_k < X_min or x_i + z_k > X_max, the extrapolation is done linearly using the slope at the first or last point of the interval. Since z_k >, no stability issue is faced in this computation.  Still using a linear interpolation, when 0<  z_k ≤, it is worth noticing that 

    ũ^n_i,k-u^n_i/ z_k=u^n_i+1-u^n_i/,

and when 0> z_k ≥ -, 

    ũ^n_i,k-u^n_i/ z_k=u^n_i-u^n_i-1/.

as a consequence, 

we define:

    ∇^, small_n,i,k=
     {[ u^n_i+1-u^n_i/,   if   0< z_k ≤; u^n_i-u^n_i-1/, if   -≤ z_k < 0;              0,      if  z_k=0. ].

This definition of B^n_i is consistent with (<ref>). Moreover,
 when  goes to 0 with fixed numerical parameters, such as Z_min and Z_max, the expression ∇_n,i,k^, large is not used at all, and
 
    B^n_i → 0= B^n,0_i= ∑_k=0  z_k<0^N_z-1 m(z_k) b(x_i) e^ z_k u^n_i-u^n_i-1/ +  m(0)b(x_i)+
    ∑_k=0  z_k >0^N_z-1 m(z_k) b(x_i) e^z_k u^n_i+1-u^n_i/.




  * Computation of ρ^n+1. 
In (<ref>), ρ^n+1 is considered in an implicit way, to make the limit scheme be consistent with the limit equation (<ref>). Since 

    ρ(t)=∫_e^u(t,x)/dx,

for >0, we define

    ρ^n+1=∑_i=0^N_x-1e^u^n+1_i/.

A closed equation on ρ^n+1 can be deduced from (<ref>). Indeed, (<ref>) yields

    e^u^n+1_i/=e^-ρ^n+1/e^(u^n_i+[-d(x_i) +B^n_i+T^n_i ] )/,

and so 
 
    ρ^n+1= e^-ρ^n+1/∑_i=0^N_x-1e^
     A^n_i
     /,

 where A^n_i denotes u^n_i+( -d(x_i)+B^n_i+T^n_i ) to simplify the notations. Eventually, ρ^n+1 is the solution of h(y)=0, where
 
    h(y)=ye^ y/-e^A^n_i_0/∑_i=0^N_x-1e^(A^n_i-A^n_i_0)/,

where A^n_i_0=max_i A^n_i has been taken apart to get an uniform estimate with respect to  on the remaining sum. It is also a solution of the equivalent equation
g(y)=0, with
 
    g(y)=-ln(y)- y +ln() + A^n_i_0+ln( ∑_i=0^N_x-1e^(A^n_i-A^n_i_0)/).

 To find ρ^n+1, a Newton's method is applied on expression (<ref>) or on (<ref>). Both expressions are smooth convex functions of ρ, and are equivalent. Hence, the Newton's method converges whatever is used. Nevertheless, it must be chosen with care.   (<ref>) is to be chosen when ρ^n+1 is close to 0 (for large values it becomes less accurate), whereas (<ref>) is more adapted when ρ^n+1 is not small, since it is more prone to accumulate numerical errors when ρ^n+1→ 0. In the effective implementation of the method, either one formulation or the other is chosen, depending on the values reached during the iterations of the algorithm. Eventually, to ensure the stability of the numerical resolution of (<ref>) when → 0, the inverse of the derivatives of h and g are analytically computed and implemented as 
 
    1/h'(y)= /+e^- y/,     
     1/g'(y)=-y/+.

 Since y>0, these two expressions are uniformly bounded with respect to  when  is fixed. As a consequence, the cost of the numerical resolution of (<ref>) does not increase with .

When >0 is fixed, the scheme (<ref>) is consistent with  (<ref>), since only quadrature formula and interpolation methods have been used to write it. The way all the terms are computed, as well as the numerical resolution of the non-linear equation (<ref>), ensures the stability of the numerical computations in the small  regime. Hence, when → 0 with fixed discretization parameters, the scheme (<ref>)
becomes 

    u^n+1_i-u^n_i/=-( d(x_i)+Cρ^n+1) + 
    
    B^n,0_i
    + τ(x_i-x_j_0),

where j_0 is such that u^n_j_0=max_iu^n_i, and B^n,0_i has been defined in (<ref>). 



We do not give a strict proof of consistency of this scheme with respect to the limiting Hamilton-Jacobi equation (<ref>), since it is out of scope of the project. However, we draw the attention to few important points which need to be taken into account while working with the scheme. In particular, the behaviour of the quantity ρ(t) is not well understood in the case of an extinction. The problem is that intuitively ρ(t) must represent the density of the population — so that when it goes to zero, we expect an extinction. However, in a Hamilton-Jacobi case even when the ρ(t) reaches zero, the population can still regrow after some time. This can be explained by the fact that after two limiting procedures (passing first to the infinite system size, and then to the infinite time horizon), the "size" of the population can not be described straightforwardly. Accurate link between the quantities obtained as a result of stochastic and PDE simulation is also a question which requires further investigation when ρ(t)≪ 1. 



  §.§.§ Case → 0: the numerical results

In this subsection we simulate the dynamics of the population by considering a small value of  and discuss the obtained results in order to compare them with previous simulations. Note that, in order to compare both, the stochastic and the Hamilton-Jacobi behaviours, the first thing to do is to obtain the simulations for the stochastic model also in the case where the HT rate is a smooth function as we do for the Hamilton-Jacobi case. We recall that, in subsection <ref> simulations for stochastic model are done with a Heaviside function as HT rate since it is a more natural choice for simulation of a jump process.





On Figure <ref> we simulate the population dynamics for = 0.01. Upon rescaling time (for chosen  time scale T=10 corresponds, in fact, to T/=1000 in previous simulations) and the variance parameter, we see the same patterns, with few differences. 

On Figure <ref>, we observe a stabilization of the mean trait, as in Figure <ref>. Similarly, on Figure <ref>, we observe cycles, but on the contrary to the PDE model oscillations are not damped. Moreover, it is worth pointing out that the duration of a cycle here corresponds to what we observe in the corresponding stochastic plot (on Figure <ref>) multiplied by =0.01. 
On Figure <ref>, we also observe a cycling behavior, but the population goes periodically extinct (i.e the population reaches exponentially small value, of order e^1/), and then reborn. On the stochastic model, it corresponds to what is illustrated in Figure <ref>. It is not surprising that this behavior is difficult to observe on the stochastic model, since very small populations are likely to go extinct.

On Figure <ref>, we can see that the population goes completely extinct. The most interesting case to comment is probably the "partial" extinction seen on <ref>. Note that despite the fact that ρ remains at 0 for some time, the population regrows. The point is that, as it was already mentioned above, this numerical parameter has no 1:1 correspondence to the population size parameter N_t/K used in stochastic model. Also note that similar behaviour of stochastic and HJ model are reproduced under a bit different values of parameters. It is caused by the rescaled time and mutation kernel, so that the rigorous link between two models is still to be developed. 


Another interesting thing to comment is that on Figure <ref> we may notice that, from the dynamics of the mean trait and the density of the population, it is easy to estimate the periods of the system. Indeed, since the system is deterministic, we just have to compute the distances between local maxima on each curve. For the stochastic system this task is more difficult, especially for a small population, because it includes filtering problem of a noisy signal. To get more accurate results in stochastic model we have to increase the time scale and number of individuals, which is costly from computational point of view. However, if our goal is to study numerically the lineages which lead to the evolutionary rescue of the population, it is still more straightforward to use the individual-based model.

To finish with, let us give some flavor on the computational cost of the simulations for each type. In Table <ref> we give a short overview of the elapsed time for the same values of parameters, but for different schemes. As expected, individual-based model is the most expensive to compute. All the computations were performed in  library of Python on MacBook Pro (Intel Core i5 processor, 2,7GHz). 





 §.§ Comparison of the theoretical analysis of the Hamilton-Jacobi equation and the numerical simulations of the stochastic model




  §.§.§ Formal computations


In this section, we propose some formal computations on the stochastic model, based on the analysis of the Hamilton-Jacobi equation performed in the previous section. To fix ideas, we assume n=1 and (<ref>)-(<ref>)-(<ref>), and we fix all constants but τ_0, as in the previous section. However, we choose the function α as a Heaviside function (this is what has been used in the simulations), which is not a smooth function, and thus will lead to minor modifications compared to the previous section.

We make a strong formal assumption: taking K≫1, we assume that the population behaves like a normally distributed random variable all the time, i.e

    ν_t^K(dx)= ρ(t)1/√(2π)s(t)e^-| x-x̅(t)|^2/2s(t)^2dx,

for some standard deviation s(t) and for x̅(t) defined in (<ref>). We expect s(t) to be of the same order as σ, but giving a general estimate for s(t) in function of x̅(t) seems intricate. The normalized size of the population ρ(t):=N_t^K/K is approximately given by (see (<ref>))

    ρ(t)= 1/Cr(x̅(t)),

where r is defined in (<ref>).

We now formally compute the evolutionary singular state x_⋆. But as α is a Heaviside function (which formally corresponds to the case when δ→0 in  (<ref>)), our derivations must be slightly adapted. In particular,
τ(x-x̅(t)) in (<ref>) has to be replaced by

    ∫_τ(x-y)ν_t^K(dy)/ρ(t),

and accordingly, recalling that the weak derivative of a Heaviside is a Dirac mass at 0,
τ'(0) in (<ref>) has to be replaced by 
    ∫_τ'(x̅(t)-y)ν(dy)/ρ(t)=2τ_0/√(2π)s(t).

We find

    x_⋆= τ_0/√(2π)s_⋆ d_r,

where s_⋆ is an unknown corresponding to the standard deviation of the population at equilibrium concentrated at x=x_⋆. Note that it corresponds to (<ref>) with δ̃:= s_⋆√(π/2).

We now try to estimate s_⋆. Formally, s_⋆ should be such that u_⋆(x):=-(x-x_⋆)^2/2s_⋆^2  is a stationary solution of (<ref>). Differentiating twice, and applying at x=x_⋆ we find 

    0=b_rσ^2 (u_⋆”(x_⋆))^2-2d_r,

(with the reasonable assumption τ”(0)=0), which gives

    s_⋆=√(σ√(b_r/2d_r)).

Numerically, we find s_⋆= 0.12. We end up with the following formula:


    x_⋆= τ_0/√(2πσ)d_r√(2d_r/b_r).




  §.§.§ Stabilization


We run a numerical test on the stochastic model corresponding to stabilization, for τ_0=0.02, and the other parameters as in Figure (<ref>). In this case, x_⋆ correspond to the mean trait of the population for large time. From, (<ref>) we find x_⋆=0.067, and from (<ref>), we obtain ρ_⋆=1.99, which corresponds to what we can see on Figure (<ref>).



  §.§.§ Threshold for cycles

Since equation (<ref>) remains unchanged, we obtain the following threshold for cycles (corresponding to (<ref>)):

    τ_cyc= 2π d_rσ√(b_r/2d_r).

With our choice of parameters, we obtain τ_cyc=0.09. This threshold 
corresponds to the numerical simulations (however, characterizing precisely whether cycles occurs or not on the numerical simulations is not easy when τ_0 is close to the threshold).



  §.§.§ Threshold for extinction

Using (<ref>), we can also find a threshold for extinction:

    τ_ext:=√(2π b_r d_rσ)√(b_r/2d_r).
 For our choice of parameters, we obtain τ_ext=0.30. 

We now compare this formula with numerical experiments on the individual-based model. They are organized as follows: we fix the birth b_r or the death rate d_r, and save the first value of τ_0 under which the extinction occurs. Then, we increase the rate and save the next HT rate under which we have an extinction. Resulting curve for the birth rate is saved on Figure <ref> (for death rate: Figure <ref>). Non-concerned parameters remain fixed as in Subsection <ref>. 

The numerical results, in particular, justify at the first glance surprising fact that the extinction threshold depends on the birth and death rate in the same manner. It seems logical to assume that while the higher birth rate contributes to a bigger survival probability even with a relatively big horizontal transfer rate, higher death rate must have an opposite effect. However, in conditions of a very "harsh" environment individuals with non-fit traits die out before they manage to transfer their genetic information to the other individuals. As a consequence, value of the critical τ increases as the value of the birth (or death) rate constant increases.








§ CONCLUSIONS


First achievement of the paper consists in an accurate numerical study conducted on the stochastic model given by a point measure (<ref>). To the best of our knowledge, in-depth analysis of the influence of the HT rate on the evolutionary dynamics has not been yet attempted. Along with its accuracy, the stochastic model reveals its limitation: an accurate theoretical description of what actually happens in each observed scenario from a mathematical point of view seems to be out of reach. However, we show that this model can be used for tracing back the lineage of the survived individuals through several generations.

On the next step, in a numerical comparative study between the stochastic (individual based) and the PDE (density) model both models exhibit the same behavior for a given set of parameters, which justifies theoretical results from <cit.>. Minor differences (in particular, the presence of damping oscillations) can be explained by a choice of a numerical scheme. However, further analysis shows that the classical PDE model defined by (<ref>) leads to instabilities if we try to pass to an asymptotic setting under the small mutation assumption. Those instabilities are then resolved by a transformation of an initial model to a Hamilton-Jacobi type equation and using an asymptotic-preserving scheme. Further advantage of this approach is that the resulting equation (<ref>) makes an easier subject of a theoretical analysis. 

Finally, in a Hamilton-Jacobi setting we manage to numerically replicate the evolutionary rescue of a small population which we observe in the stochastic model. This phenomena is illustrated for stochastic, PDE and HJ simulation on Figure <ref>. On Figures <ref>-<ref>  we trace the moment of the regrowth for different models. Figure <ref> show the state of the population at certain moment of time: we see how the individuals are centered around a mean trait. For PDE and HJ model (red and green line respectively) we simply plot the density function, and on the first (blue) plot we approximate a histogram which describes ratio N_t/K sorted by traits in stochastic model.  Stochastic simulations show the evolutionary rescue in more distinct manner: we see how the very small number of non-mutated individuals rescues the whole population from extinction (transition from <ref> to <ref>). On the contrary, the transition on the PDE model is dumped, and the regrowth is not clearly visible. It is due to, again, numerical instability of the PDE scheme for small values of the density function. Finally, HJ explicitly shows how the cycle occurs: the regrow of the "fit" individuals we see in stochastic plot is reproduced by a change of the maximum point (see again <ref> to <ref>). 


We highlight again that in order to compare the models on a more applied level, we have to give a formal definition of a quantity represented by ρ in a Hamilton-Jacobian case. In this work we have made few steps toward the theoretical analysis of the limiting equation and an accurate description of each event (evolutionary rescue, extinction, etc) in terms of solutions of a PDE. Even though establishing a rigorous mathematical link between the behavior observed in the individual-based model and the Hamilton-Jacobi equation is out of scope of this project, obtained analytical results already give a flavor of how the analysis of the evolutionary dynamics can be simplified in the future.








§ ACKNOWLEDGEMENTS

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 639638). The project has received funding from  the Chair “Modélisation Mathématique et Biodiversité” of Veolia Environnement-Ecole Polytechnique-Museum National d’Histoire Naturelle-Fondation X. 

plain	



