Horizontal gene transfer: numerical comparison between stochastic and deterministic approaches

TL;DR

This paper compares stochastic and deterministic models of horizontal gene transfer to evaluate how well the simplified Hamilton-Jacobi approach captures the complex dynamics observed in individual-based simulations, especially regarding evolutionary rescue.

Contribution

It provides a numerical comparison between stochastic jump processes, deterministic integrodifferential equations, and Hamilton-Jacobi models for horizontal gene transfer, highlighting their similarities and differences.

Findings

Hamilton-Jacobi model reproduces key qualitative behaviors of stochastic models

Deterministic models effectively approximate large population dynamics

Insights into the phenomenon of evolutionary rescue in bacterial populations

Abstract

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 integrod-ifferential 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.

Tables

Table 1. Table 1 : Elapsed time for simulation of population dynamics for different models (other parameters are fixed to values used throughout all the other simulations, τ = 0.5 𝜏 0.5 \tau=0.5 ).

	$\Delta =0.1,T=10$	$\Delta =0.01,T=10$
SM ($N=1000$)	3.883s	38.145s
SM ($N=10000$)	15.805s	153.255s
PDE ($\epsilon =1$)	0.186s	1.673s
HJ ($\epsilon ={10}^{-2}$)	0.191s	1.636s
HJ ($\epsilon ={10}^{-6}$)	0.195s	1.656s

Equations

$${\nu }_t^K(dx)=\frac{1}{K}\sum _{i=0}^{N_t^K}{\delta }_{X_i(t)}(dx),$$

$$h_K(x,y,\nu )=h_K(x-y,N)={\tau }_0\frac{\alpha (x-y)}{N/K},$$

$$\alpha (z)=\{\begin{array}{rlrl}& 0& & \text{if }z<-\delta \\ & 1& & \text{if }z>+\delta \end{array},{\alpha }^{\prime }(0)=\frac{1}{2\delta },$$

$$\begin{gathered} L^{K}F(\nu )=\sum _{i=1}^{N}b(x_i){\int }_{\text{mathds}R}\left(F\left(\nu +\frac{1}{K}{\delta }_y\right)-F(\nu )\right)m(x_i,dy) \\ +\sum _{i=1}^N\left(d(x_i)+C\frac{N}{K}\right)\left(F\left(\nu -\frac{1}{K}{\delta }_{x_i}\right)-F(\nu )\right) \\ +\sum _{i,j=1}^N{h}_K(x_i,x_j,\nu )\left(F\left(\nu +\frac{1}{K}{\delta }_{x_i}-\frac{1}{K}{\delta }_{x_j}\right)-F(\nu )\right). \end{gathered}$$

$$\left\{\begin{array}[]{rcl}\partial_{t}f(t,x)&=&-(d(x)+C\rho_{1}(t))f(t,x)+\displaystyle\int_{\mathbb{R}^{n}}m(x-y)b(y)f(t,y)dy+\\ &&\qquad f(t,x)\int_{\mathbb{R}^{n}}\tau(x-y)\frac{f(t,y)}{\rho_{1}(t)}dy,\quad(t,x)\in\mathbb{R}_{+}\times\mathbb{R}^{n},\\ \rho_{1}(t)&=&\displaystyle\int_{\mathds{R}}f(t,x)dx,\\ f(0,x)&=&f^{0}(x)>0,\end{array}\right.$$

$$\tau (y-x):={\tau }_0\left[\alpha (x-y)-\alpha (y-x)\right]$$

$$m_{\epsilon }(x-y)=\frac{1}{{\epsilon }^n}m\left(\frac{x-y}{\epsilon }\right).$$

$${\int }_{{\mathbb{R}}^n}{m}_{\epsilon }(x-y)b(y)f(t,y)dy={\int }_{{\mathbb{R}}^n}m(z)b(x+\epsilon z)f(t,x+\epsilon z)dz.$$

$$\left\{\begin{array}[]{rcl}\varepsilon\partial_{t}f_{\varepsilon}(t,x)&=&-(d(x)+C\rho_{\varepsilon}(t))f_{\varepsilon}(t,x)+\int_{\mathbb{R}^{n}}m(z)b(x+\varepsilon z)f_{\varepsilon}(t,x+\varepsilon z)dz+\\ &&\qquad f_{\varepsilon}(t,x)\int_{\mathbb{R}^{n}}\tau(x-y)\frac{f_{\varepsilon}(t,y)}{\rho_{\varepsilon}(t)}dy,\\ \rho_{\varepsilon}(t)&=&\displaystyle\int_{\mathds{R}}f_{\varepsilon}(t,x)dx,\\ f_{\varepsilon}(0,x)&=&f_{\varepsilon}^{0}(x)>0.\end{array}\right.$$

$$u_{\epsilon }(t,x):=\epsilon \ln (f_{\epsilon }(t,x)).$$

$$\begin{gathered} {\partial }_t{u}_{\epsilon }=-(d(x)+C{\rho }_{\epsilon }(t))+{\int }_{{\mathbb{R}}^n}m(z)b(x+\epsilon z)\exp \left\{\frac{u_{\epsilon }(t,x+\epsilon z)-u_{\epsilon }(t,x)}{\epsilon }\right\}dz \\ +{\int }_{{\mathbb{R}}^n}\tau (x-y)\frac{f_{\epsilon }(t,y)}{{\rho }_{\epsilon }(t)}dy. \end{gathered}$$

$${\partial }_{t}u=-(d(x)+C\rho (t))+b(x){\int }_{{\mathbb{R}}^n}m(z)e^{z\cdot {\nabla }_{x}u}dz+\tau (x-\overline{x}(t)),$$

$$\bar{x}(t)=\mathrm{argmax}u(t,\cdot ).$$

$$\underset{x\in {\mathbb{R}}^n}{\sup }u(t,x)=0\text{when}\rho (t)>0.$$

$$f_{\epsilon }(t,x)\rightharpoonup \rho (t){\delta }_{\bar{x}(t)}(dx),\text{when }\epsilon \to 0,$$

$${\partial }_{t}u(t,\bar{x}(t))=0,$$

$$\rho (t)=\frac{b(\bar{x}(t))-d(\bar{x}(t))+\tau (0)}{C}=\frac{b(\bar{x}(t))-d(\bar{x}(t))}{C},$$

$$b(x)=b_r>0,$$

$$d(x)=d_r{x}^2,d_r>0,$$

$$m(z)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{z^2}{2{\sigma }^2}},$$

$$u(t,\cdot )\text{ reaches its maximum on a single point }\bar{x}(t),$$

$$\bar{x}(t)\text{ is a non-degenerate maximum, i.e }{\nabla }_x^{2}u(t,x)<0,$$

$$\bar{x}(t)\text{ is smooth with respect to }t.$$

$$f_{\epsilon }^0(x)=\frac{1}{\sqrt{\epsilon }}{e}^{-\frac{x^2}{2\epsilon }}.$$

$$0\le \bar{x}(t)\le \frac{{\tau }_0}{2d\delta }.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\bar{x}(t)={\left[-{\nabla }_x^{2}u(t,\bar{x}(t))\right]}^{-1}\cdot \left({\nabla }_{x}r(\bar{x}(t))+{\nabla }_x\tau (0)\right),$$

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

$${\nabla }_{x}u(t,\bar{x}(t))=0,$$

$$x_{\star }=\frac{{\tau }_0}{2d_r\delta }.$$

$$\bar{x}(T^{-}):=\underset{\begin{array}{c}t\to T\\ t<T\end{array}}{\lim }\bar{x}(t),\bar{x}(T^{+}):=\underset{\begin{array}{c}t\to T\\ t>T\end{array}}{\lim }\bar{x}(t).$$

$$F_{\bar{x}}(x):=r(x)+\tau (x-\bar{x}),$$

$$F_{\bar{x}(t)}(0)>F_{\bar{x}(t)}(\bar{x}(t)).$$

Full text

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 M Henry et al. (2013), Lili et al. (2007). 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 Kamimura et al. (2011)). 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, Billiard et al. (2015), Fournier and Méléard (2004)), 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 Billiard et al. (2018, 2016), Ferrière and Tran (2009). 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 Lorz et al. (2011), Perthame and Barles (2008)). 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 2, classical explicit schemes do not preserve the asymptotic behavior of the solution if the time rescaling step goes to [math]. 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 Crandall and Lions (1984). More generally, the numerical approximation problem for solutions of Hamilton-Jacobi equations is treated in Achdou et al. (2013).

This paper is structured as follows: in Section 1 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 2. 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.

1 Model

1.1 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

[TABLE]

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\int\nu^{K}_{t}(dx)$ is the size of the population at time $t$, and $X_{i}(t)\in\mathds{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 $\int_{\mathds{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 $C\dfrac{N_{t}^{K}}{K}$ (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,\nu)$, 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,\nu)$ to be in the particular form

[TABLE]

where $N=K\int_{\mathds{R}}\nu(dx)$ is the number of individuals, $\tau_{0}>0$ is a constant and $\alpha$ is either a Heaviside, or a smooth bounded function, such that for a small $\delta>0$:

[TABLE]

where $\delta$ is the stiffness parameter. We introduce $\delta$ 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 $\nu=\frac{1}{K}\sum_{i=1}^{N}\delta_{x_{i}}$ and a generic measurable bounded function $F$, the generator of the process is then given by:

[TABLE]

It is standard to construct the measure-valued process $\nu^{K}$ as the solution of a stochastic differential equation driven by Poisson point measures and to derive moment and martingale properties (see for instance Fournier and Méléard (2004)).

1.2 The PDE model

It is proven (see in particular Billiard et al. (2018), Champagnat et al. (2008)) that for $K\to+\infty$ the stochastic process defined by a sequence of point measures given by (1) converges in probability to the unique solution of a non-linear integro-differential equation. This equation is given by:

[TABLE]

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

[TABLE]

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 $\varepsilon>0$ we define

[TABLE]

With a change of variable $z=\frac{x-y}{\varepsilon}$ we can rewrite the mutation term at $(t,x)$ as

[TABLE]

On the other hand, when $\varepsilon$ is small, the effect of mutations can only be observed in a larger time scale. Thus, we rescale time with $t\mapsto\frac{t}{\varepsilon}$.

We end up with the following system, for $\varepsilon>0$, and $(t,x)\in\mathbb{R}_{+}\times\mathbb{R}^{n}$:

[TABLE]

1.3 The Hamilton-Jacobi limit

We now derive the limiting problem (5) when $\varepsilon\to 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 (5) often give rise to a concentration phenomenon, i.e the convergence of $f_{\varepsilon}$ towards a Dirac mass when $\varepsilon\to 0$ (see Perthame and Barles (2008), Diekmann et al. (2005)). The usual way to deal with these asymptotics is to perform a Hopf-Cole transformation (or WKB ansatz), i.e to consider

[TABLE]

This change of variable comes from the intuition that a Dirac mass is no more than a narrow Gaussian, and more precisely that $f_{\varepsilon}$ should behave like a Gaussian of variance $\varepsilon$ when $\varepsilon\to 0$. Accordingly, we expect $u_{\varepsilon}$ to have a non singular limit when $\varepsilon\to 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_{\varepsilon}$ when $\varepsilon\to 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 (6) into (5) we deduce that $u_{\varepsilon}$ satisfies

[TABLE]

Formally, at the limit $\varepsilon\to 0$, $u_{\varepsilon}$ converges to a continuous function $u$ which satisfies the following Hamilton-Jacobi equation in the ”viscosity” sense:

[TABLE]

where $\rho(t)\geq 0$ is the weak limit of $\rho_{\varepsilon}(t)$ and

[TABLE]

We formally assume here and in the following that the definition of $\bar{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 Barles (1994)).

This framework is convenient because most of the information is contained in the dynamics of $\bar{x}(t)$. See the next section for further analysis.

1.4 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 (8). We point out that the calculations are only formal, since rigorous proofs are intricate and beyond the scope of this paper.

1.4.1 Generality

From an integration of (5) 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. $\rho_{\varepsilon}(t)$ is bounded from above, uniformly in $t\geq 0$ and $\varepsilon>0$. From this and $\rho_{\varepsilon}(t)=\int_{\mathbb{R}^{n}}e^{\frac{u_{\varepsilon}(t,x)}{\varepsilon}}dx$, we deduce that $\forall t>0$, $\sup\limits_{x\in\mathbb{R}^{n}}u(t,x)\leq 0$ and the following constraint holds:

[TABLE]

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

[TABLE]

where $\delta_{\bar{x}(t)}$ denotes, as usually, the Dirac measure centered in $\bar{x}(t)$. From (10) it is possible to formally derive a formula for $\rho$. Indeed, either $\rho(t)=0$ or $\rho(t)>0$ and

[TABLE]

which implies

[TABLE]

for $\tau$ defined in (4).

Having above definitions in hand, we can now perform a formal analysis on the dynamics of $\bar{x}(t)$, defined below in (15). Our aim is to show how the behaviour of the system can be analyzed within the framework of a Hamilton-Jacobi equation (8). To fix ideas, we fix all constants but $\tau_{0}$ and we assume (12)-(14) as follows:

[TABLE]

and the transfer function $h_{K}(x,y,\nu)$ is defined in (2). Moreover we work under the following assumptions:

[TABLE]

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

[TABLE]

1.4.2 Smooth dynamics $\bar{x}(t)$.

The following statement deals with the smooth dynamics of $\bar{x}(t)$, i.e in the regime where no jump occurs in the dynamics of $\bar{x}(t)$.

Statement 1**.**

Under assumptions (12)-(15), the function $t\mapsto\bar{x}(t)$ is an increasing function which satisfies the following inequality $\forall t\geq 0$:

[TABLE]

More precisely, $\bar{x}(t)$ satisfies the canonical equation

[TABLE]

where

[TABLE]

and $\nabla^{2}_{x}u$ denotes the Hessian of $u$ with respect to the $x$ variable.

Proof.

Under the above assumptions we can derive the dynamics of $\bar{x}(t)$, referred to as the canonical equation in the literature (see for instance Mirrahimi and Roquejoffre (2016)). Indeed, starting from

[TABLE]

a differentiation with respect to $t$ gives (17). Equation (17) has a unique singular point $x_{\star}$, which satisfies $r^{\prime}(x_{\star})+\tau^{\prime}(0)=0$, with $\tau$ defined in (2) and $r$ in (18). We find

[TABLE]

Note that $t\mapsto\bar{x}(t)$ is increasing when $\bar{x}(t)<x_{\star}$ and decreasing when $\bar{x}(t)>x_{\star}$. Besides, from the initial condition (16), we have $\bar{x}(0)=0$, and consequently $0\leq\bar{x}(t)\leq x_{\star}\quad\forall t$. ∎

1.4.3 Evolutionary rescue.

In general, the canonical equation (17) 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 $\bar{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 5(b). We thus expect $\bar{x}(t)$ to possibly jump periodically, and to follow (17) between two jumps. We now try to characterize the possible jumps. For $T>0$, we denote

[TABLE]

Statement 2**.**

We assume that (12)-(15) hold until a time $T>0$, such that $u(T,\cdot)$ reaches its maximum on $\bar{x}(T^{-})$ and on another point $\tilde{x}$. Then $\tilde{x}=0$ and $\bar{x}(t)$ will jump towards [math] at time $T$, i.e $\bar{x}(T^{+})=0$.

Proof.

From assumption (15), we have $\forall t\in[0,T]$ that $u(t,\cdot)$ is concave non-degenerate on $[\bar{x}(t)\pm\theta],$ with $\theta>0$. For simplicity, we further assume $\delta\leq\theta$, where $\delta$ is defined in (3).

First, let us show that $\tilde{x}=0$. We define the fitness function of trait $x$ in a population concentrated in $\bar{x}$:

[TABLE]

where $r$ and $\tau$ are respectively defined in (18) and (4). Note that we have $\partial_{t}u(t,x)=F_{\bar{x}(t)}(x)-C\rho(t)$, for $t<T$. But $\tilde{x}\not\in[\bar{x}(t)\pm\delta]$ and the choice of parameters (12)-(13)-(3) implies $\tilde{x}$ must maximize $F_{\bar{x}(T^{-})}(\cdot)$, hence $\tilde{x}=0$.

The second step is to prove that there will be an actual jump towards [math], i.e $\bar{x}(T^{+})=0$. First, note that there exists a small $\eta>0$ such that $\forall t\in(T-\eta,T)$, $u(t,\bar{x}(t))=0$ and $u(t,0)<0$. Let us fix $t\in(T-\eta,T)$. We have $F_{\bar{x}(t)}(0)\geq F_{\bar{x}(t)}(\bar{x}(t))$, and we claim that the inequality is strict. Indeed, since $t\mapsto x(t)$ is increasing, $F_{\bar{x}(t)}(\bar{x}(t))$ is decreasing, whereas $F_{\bar{x}(t)}(0)$ is constant (as long as $\eta$ is small enough such that $\bar{x}(T-\eta)>\delta$). We end up with

[TABLE]

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

[TABLE]

Indeed we have from the definition of $F_{\bar{x}}(x)$ that

[TABLE]

Consequently that for all $\lambda\in[0,1]$

[TABLE]

It shows that [math] remains the fittest trait during all the process of invasion, and therefore that [math] will actually invade, i.e that $\bar{x}(t)$ will actually jump towards [math] at time $T^{+}$. ∎

1.4.4 Threshold for cycles

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

We recall that a jump corresponds to a rescue of the population concentrated at $\bar{x}(t)$ by the small population with trait $x=0$. It is possible only if $\bar{x}(t)>\delta$ and if [math] is fitter than $\bar{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 [math] is fitter than $\bar{x}(t)$ if and only if

[TABLE]

But if no jump occurs, $\bar{x}(t)$ formally follows (17), thus $\bar{x}(t)<x_{\star}$ and $\bar{x}(t)$ converges to $x_{\star}$ when $t\to+\infty$ (with $x_{\star}$ is defined in (19)).

Statement 3**.**

Under assumptions (12)-(15), the evolutionary rescue phenomena occurs if and only if

[TABLE]

Note that the condition $\tau_{0}>\tau_{cyc}$ is equivalent to $x_{resc}<x_{\star}$ , which are defined respectively in (19) and (20).

1.4.5 Threshold for extinction.

The population is said to be ”extinct” at time $t$ if $\rho(t)=0$. According to (11), we define $x_{ext}$ as to solve $r(x_{ext})=0$, i.e

[TABLE]

that is, a population concentrated at trait $\bar{x}$ is extinct iff $\bar{x}\geq x_{ext}$.

The picture is simple in the case of stabilization without cycles, i.e when $\tau_{0}\leq\tau_{cyc}$ (see (22)). In this case, we recall that $\bar{x}(t)$ formally follows (17) for all $t>0$, thus $\bar{x}(t)<x_{\star}$ and $\bar{x}(t)$ converges to $x_{\star}$ when $t\to+\infty$ (where $x_{\star}$ is defined in (19)). Thus, if $x_{\star}\leq x_{ext}$, we have $\rho(t)>0$ for all $t>0$; on the contrary, if $x_{\star}>x_{ext}$, there exists a time $t_{ext}>0$ for which $\rho(t)=0$ for all $t\geq t_{ext}$. It gives a sharp threshold for extinction of the population: indeed, the population eventually gets extinct if and only if $x_{\star}>x_{ext}$, which naturally leads us to the following statement.

Statement 4**.**

Under assumptions (12)-(15), if $\tau_{0}\leq\tau_{cyc}$, then the population eventually gets extinct if and only if

[TABLE]

We point out that, surprisingly enough, $\tau_{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 [math] 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 $\tau_{0}>\tau_{cyc}$. In this case, $\bar{x}(t)$ will follow (17) and will periodically jump to $x=0$. First, note that if $x_{\star}<x_{ext}$, $\bar{x}(t)$ remains below $x_{ext}$ for all $t$ and the population does not get extinct:

[TABLE]

The most intricate case is when $x_{\star}>x_{ext}$, which contains cases of extinction and non-extinction, depending on whether the jump of $\bar{x}(t)$ towards [math] happens before or after $\bar{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 $\bar{x}(T^{-})\leq x_{ext}$, where $T$ is the time where the jump of $\bar{x}(t)$ towards [math] occurs. However, we are not able to give a satisfactory formula or estimate on $T$.

Besides, when the jump of $\bar{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)+\tau_{0}>0$, which is equivalent to

[TABLE]

If $\tau\geq\tau_{sus}$, the population eventually gets extinct. If $\tau<\tau_{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 5(c). 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^{-\frac{1}{\varepsilon}}$, with $\varepsilon$ the variance of the mutation kernel), and the probability that every individual dies is bigger than the survival of the population.

Statement 5**.**

Assume (12)-(15) and $\tau_{0}>\tau_{cyc}$.

• •

if $\tau_{0}\leq\tau_{ext}$, the population never gets extinct.

• •

the evolutionary rescue effectively manages to sustain the population if and only if $\tau_{0}<\tau_{sus}:=b_{r}$.

1.4.6 Characteristics of a Hamilton-Jacobi equation

Denoting

[TABLE]

from (8) we have $\partial_{t}u(t,x)+H(t,x,\nabla_{x}u(t,x))=0$. Since $H$ is convex in the $p$ variable, we have the following representation formula (see Lions (1982)).

[TABLE]

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

Every $\gamma$ which is admissible as a minimizer in (26) is called a characteristic of the Hamilton-Jacobi equation (8). Note that every characteristic $\gamma$ formally satisfies the condition

[TABLE]

(27) holds because $\gamma$ is a critical point of the functional defined in (26). Note that if we replace $H$ by $\tilde{H}(x,p)=-\frac{x^{2}}{2}+\frac{p^{2}}{2}+1,$ the Legendre transform of $\tilde{H}$ can be computed explicitly:

[TABLE]

Then (27) becomes

[TABLE]

2 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 (5) to a limit (7). Throughout this section we define the birth, death rates and the mutation kernel to those given in (12)-(14) respectively, with the parameters fixed throughout all the experiments to $b\equiv 1$, $d_{r}\equiv 1$, $C\equiv 0.5$ respectively (unless otherwise stated).

2.1 Stochastic model

2.1.1 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 $\sigma^{0}$ so that the resulting vector $X^{0}\in\mathds{R}^{N^{0}}$. We know that in a time step $\Delta$, 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\text{ 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}\leq T_{HT}\leq\Delta$), this event is extremely rare.

We simulate the population of initial size $N_{0}=10000$ up to time $T=1000$ with $\Delta=0.01$, with the parameters being defined at the beginning of the section, and $\alpha$ 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 1.2 and 1.3), 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 $\tau_{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 1). First possibility, for small values of $\tau_{0}$, is the stabilization (Figure 1(a)). 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 $\tau_{0}$, is the cycling behavior (Figure 1(b)). 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 [math]) manage to regrow and eventually repopulate the environment, which launches the cycle again.

The last possibility, for large values of the horizontal transfer rate $\tau_{0}$, is the extinction of the population (Figure 1(c)). 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 $\tau_{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 $\tau_{0}$ is getting closer to a ”critical” (numerically we observe an almost sure extinction at $\tau_{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 2, 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 1.2.

2.1.2 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 3. 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 Kingman (1982), or Arenas and Posada (2014, 2010) 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.

2.2

Numerical scheme for the PDE model

In this subsection, a numerical scheme for (5) is presented, and its properties are numerically investigated. For the discretization of (5), we consider a bounded space of traits $[X_{\mathrm{min}},X_{\mathrm{max}}]$, discretized with $N_{x}$ points. Denoting $N_{x}$ the number of discretization points of the interval $[X_{\mathrm{min}},X_{\mathrm{max}}]$, we define

[TABLE]

and

[TABLE]

We consider the time interval $[0,T_{\mathrm{max}}]$, discretized with $N_{t}$ points $t_{n}=n{\Delta t}$, for $0\leq n\leq N_{t}-1$, and where ${\Delta t}$ is defined as

[TABLE]

The approximations of the solution $f$ of (5) at $(t_{n},x_{i})$, and of its density $\rho$ at $t_{n}$ are denoted $f^{n}_{i}$ and $\rho^{n}$ respectively. We recall that the initial condition $f^{0}$ is a smooth function of $x$ given in (16) and the initial density $\rho^{0}$ is computed using a left-point quadrature rule for $f^{0}$ as follows:

[TABLE]

The scheme is written with an explicit Euler scheme, in which the integrals are computed with a left-point quadrature rule. For $n\geq 1$ and $0\leq i\leq N_{x}-1$, it reads

[TABLE]

In (29), the convolution product $\left[m*(bf)\right]^{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_{\mathrm{min}}$ and $Z_{\mathrm{max}}$, and the number $N_{z}$ of discretization points be given. The grid in $z$ is defined as

[TABLE]

where ${\Delta z}=\left(Z_{\mathrm{max}}-Z_{\mathrm{min}}\right)/\left(N_{z}-1\right)$. When $x_{i}+\varepsilon z_{k}\in\left[X_{\mathrm{min}},X_{\mathrm{max}}\right]$, the value of $f(t_{n},x_{i}+\varepsilon z_{k})$ is approximated by linear interpolation of the $(f^{n}_{i})_{0\leq i\leq N_{x}-1}$. When $x_{i}+\varepsilon z_{k}<X_{\mathrm{min}}$, or $x_{i}+\varepsilon z_{k}>X_{\mathrm{max}}$, it is computed with a linear extrapolation of the $(f^{n}_{i})_{0\leq i\leq N_{x}-1}$, using the slope at the corresponding end of the $X$ domain. Using the notation $f^{n}(x_{i}+\varepsilon z_{k})$ for the approximation of $f(t_{n},x_{i}+\varepsilon z_{k})$, we then define

[TABLE]

2.2.1 Case $\varepsilon=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 $\varepsilon$ 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 4, simulations in overall correspond to those of the stochastic model. Indeed, when the HT rate $\tau_{0}$ is small enough the population rapidly stabilizes around its equilibrium state (see Figure 4(a)), as in the stochastic simulations. Further similarity between two models is that in both cases the optimal trait is shifted a bit above [math]. It is caused by the HT phenomenon.

For larger values of $\tau_{0}$, where we would expect to have distinguishable cycles, we observe indeed damped oscillations, see Figure 4(b). We stress out that for the stochastic model it is not the case, see Figure 1(b). 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 (6), and to write a numerical scheme which converges uniformly when $\varepsilon\to 0$. This is what the next subsection is devoted to.

On Figure 4(c), we observe that as $\tau_{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 4(d) we observe a full extinction of the population without regrowth. We will give further insights on those two cases in the next subsection.

2.3 The scheme for the Hamilton-Jacobi equation

2.3.1 Case $\varepsilon\to 0$: description of the numerical scheme

As the rescaling parameter $\varepsilon$ goes to [math], the model given by (7) gets closer to its limiting state (8). However, numerical approximation of the (5) for $\varepsilon\ll 1$ is not a trivial task. Indeed, for small $\varepsilon$, the solution $f^{\varepsilon}$ of (5), 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 $\varepsilon\to 0$, and reaching the asymptotic regime with this scheme is not possible. In this part, we present a numerical scheme for (5) which enjoys stability properties in the limit $\varepsilon\to 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 $\varepsilon$, a scheme for the solution $u^{\varepsilon}$ of (7) which enjoys the Asymptotic Preserving (AP) property is proposed here. Such schemes have been introduced in Klar (1998, 1999), Jin (1999), their properties are often summarized by the following diagram:

[TABLE]

It should be understood as follows: when the parameter $\varepsilon>0$ is fixed, the scheme $S^{h}_{\varepsilon}$ is consistent with the $\varepsilon$-dependent problem $P_{\varepsilon}$. When $\varepsilon$ goes to [math], the solution of $P_{\varepsilon}$ converges to the solution of the limit problem $P_{0}$. The AP scheme $S^{h}_{\varepsilon}$ is stable along the transition to the asymptotic regime. It means that, when $\varepsilon$ goes to [math] 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 $\varepsilon$ is going to [math], 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 $\rho$. If all of them are correctly defined, the scheme we propose reads

[TABLE]

where $B^{n}_{i}$ stands for an approximation of

[TABLE]

and $T^{n}_{i}$ is for

[TABLE]

Here, we used the notations and discretization grids defined in the beginning of Section 2.2, and the dependencies in $\varepsilon$ are omitted to simplify the notations. In what follows, we present how $T^{n}_{i}$, $B^{n}_{i}$ and $\rho^{n+1}$ can be computed in a way that ensures they are consistent with their definition for fixed $\varepsilon$, that they can be computed with a constant computational cost with respect to $\varepsilon$, and that their asymptotic behavior when $\varepsilon$ goes to [math] is meeting the continuous one (8).

• •

Computation of $T^{n}_{i}$. The direct approximation of (32) with a quadrature rule is consistent for $\varepsilon\sim 1$. However, since $f$ is expected to concentrate when $\varepsilon\to 0$, it lacks precision in the asymptotic regime. Especially, the convergence of $f/\rho$ to a Dirac is not ensured when the integral is approximated directly. Remarking that

[TABLE]

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

[TABLE]

For fixed $\varepsilon$, (33) is consistent with (32). Since all the arguments of the exponentials are nonpositive, the limit of (33) for small $\varepsilon$ can be read on that expression. Denoting $j_{0}$ the index such that

[TABLE]

and supposing that there exists a unique such $j_{0}$, the limit of (33) for small $\varepsilon$ is

[TABLE]

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

• •

Computation of $B^{n}_{i}$. Once again, the numerical approximation of (31) is done with a quadrature in the integral. Using the notations of Section 2.2, a grid in $z$ is defined. The functions $m$ and $b$ are respectively evaluated at $z_{k}$ and $x_{i}+\varepsilon z_{k}$, but the interpolation of $u^{n}$ at $x_{i}+\varepsilon z_{k}$ has to be done with special care to make the scheme enjoy the expected asymptotic behavior. Using a left-point quadrature rule, (31) is approximated by

[TABLE]

where $\nabla_{n,i,k}^{\varepsilon}$ stands for an approximation of

[TABLE]

In both cases, it is computed with a linear interpolation of the values $u^{n}_{i}$. Hence, $\nabla_{n,i,k}^{\varepsilon,large}$ is given by

[TABLE]

where $\tilde{u}^{n}_{i,k}$ is computed as the linear interpolation of $(u^{n}_{i})_{1\leq i\leq N_{x}}$ at $x_{i}+\varepsilon z_{k}$. If $x_{i}+\varepsilon z_{k}<X_{\mathrm{min}}$ or $x_{i}+\varepsilon z_{k}>X_{\mathrm{max}}$, the extrapolation is done linearly using the slope at the first or last point of the interval. Since $\varepsilon z_{k}>{\Delta x}$, no stability issue is faced in this computation. Still using a linear interpolation, when $0<\varepsilon z_{k}\leq{\Delta x}$, it is worth noticing that

[TABLE]

and when $0>\varepsilon z_{k}\geq-{\Delta x}$,

[TABLE]

as a consequence, we define:

[TABLE]

This definition of $B^{n}_{i}$ is consistent with (31). Moreover, when $\varepsilon$ goes to [math] with fixed numerical parameters, such as $Z_{\mathrm{min}}$ and $Z_{\mathrm{max}}$, the expression $\nabla_{n,i,k}^{\varepsilon,large}$ is not used at all, and

[TABLE]

• •

Computation of $\rho^{n+1}$. In (30), $\rho^{n+1}$ is considered in an implicit way, to make the limit scheme be consistent with the limit equation (8). Since

[TABLE]

for $\varepsilon>0$, we define

[TABLE]

A closed equation on $\rho^{n+1}$ can be deduced from (30). Indeed, (30) yields

[TABLE]

and so

[TABLE]

where $A^{n}_{i}$ denotes $u^{n}_{i}+{\Delta t}\left(-d(x_{i})+B^{n}_{i}+T^{n}_{i}\right)$ to simplify the notations. Eventually, $\rho^{n+1}$ is the solution of $h(y)=0$, where

[TABLE]

where $A^{n}_{i_{0}}=\max\limits_{i}A^{n}_{i}$ has been taken apart to get an uniform estimate with respect to $\varepsilon$ on the remaining sum. It is also a solution of the equivalent equation $g(y)=0$, with

[TABLE]

To find $\rho^{n+1}$, a Newton’s method is applied on expression (37) or on (38). Both expressions are smooth convex functions of $\rho$, and are equivalent. Hence, the Newton’s method converges whatever is used. Nevertheless, it must be chosen with care. (37) is to be chosen when $\rho^{n+1}$ is close to [math] (for large values it becomes less accurate), whereas (38) is more adapted when $\rho^{n+1}$ is not small, since it is more prone to accumulate numerical errors when $\rho^{n+1}\to 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 (36) when $\varepsilon\to 0$, the inverse of the derivatives of $h$ and $g$ are analytically computed and implemented as

[TABLE]

Since $y>0$, these two expressions are uniformly bounded with respect to $\varepsilon$ when ${\Delta t}$ is fixed. As a consequence, the cost of the numerical resolution of (36) does not increase with $\varepsilon$.

When $\varepsilon>0$ is fixed, the scheme (30) is consistent with (7), 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 (36), ensures the stability of the numerical computations in the small $\varepsilon$ regime. Hence, when $\varepsilon\to 0$ with fixed discretization parameters, the scheme (30) becomes

[TABLE]

where $j_{0}$ is such that $u^{n}_{j_{0}}=\max\limits_{i}{u^{n}_{i}}$, and $B^{n,0}_{i}$ has been defined in (35).

We do not give a strict proof of consistency of this scheme with respect to the limiting Hamilton-Jacobi equation (8), 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 $\rho(t)$ is not well understood in the case of an extinction. The problem is that intuitively $\rho(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 $\rho(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 $\rho(t)\ll 1$.

2.3.2 Case $\varepsilon\to 0$: the numerical results

In this subsection we simulate the dynamics of the population by considering a small value of $\varepsilon$ 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 2.1 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 5 we simulate the population dynamics for $\varepsilon=0.01$. Upon rescaling time (for chosen $\varepsilon$ time scale $T=10$ corresponds, in fact, to $\frac{T}{\varepsilon}=1000$ in previous simulations) and the variance parameter, we see the same patterns, with few differences.

On Figure 5(a), we observe a stabilization of the mean trait, as in Figure 1(a). Similarly, on Figure 5(b), 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 1(b)) multiplied by $\varepsilon=0.01$. On Figure 5(c), we also observe a cycling behavior, but the population goes periodically extinct (i.e the population reaches exponentially small value, of order $e^{1/\varepsilon}$), and then reborn. On the stochastic model, it corresponds to what is illustrated in Figure 2. 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 5(d), we can see that the population goes completely extinct. The most interesting case to comment is probably the ”partial” extinction seen on 5(c). Note that despite the fact that $\rho$ remains at [math] 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 $\frac{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 5(b) 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 1 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 numpy library of Python on MacBook Pro (Intel Core i5 processor, 2,7GHz).

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

2.4.1 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 (3)-(12)-(13), and we fix all constants but $\tau_{0}$, as in the previous section. However, we choose the function $\alpha$ 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\gg 1$, we assume that the population behaves like a normally distributed random variable all the time, i.e

[TABLE]

for some standard deviation $s(t)$ and for $\bar{x}(t)$ defined in (9). We expect $s(t)$ to be of the same order as $\sigma$, but giving a general estimate for $s(t)$ in function of $\bar{x}(t)$ seems intricate. The normalized size of the population $\rho(t):=\frac{N_{t}^{K}}{K}$ is approximately given by (see (11))

[TABLE]

where $r$ is defined in (18).

We now formally compute the evolutionary singular state $x_{\star}$. But as $\alpha$ is a Heaviside function (which formally corresponds to the case when $\delta\to 0$ in  (19)), our derivations must be slightly adapted. In particular, $\tau(x-\bar{x}(t))$ in (8) has to be replaced by

[TABLE]

and accordingly, recalling that the weak derivative of a Heaviside is a Dirac mass at [math], $\tau^{\prime}(0)$ in (19) has to be replaced by

[TABLE]

We find

[TABLE]

where $s_{\star}$ is an unknown corresponding to the standard deviation of the population at equilibrium concentrated at $x=x_{\star}$. Note that it corresponds to (19) with $\tilde{\delta}:=s_{\star}\sqrt{\pi/2}$.

We now try to estimate $s_{\star}$. Formally, $s_{\star}$ should be such that $u_{\star}(x):=\frac{-(x-x_{\star})^{2}}{2s_{\star}^{2}}$ is a stationary solution of (8). Differentiating twice, and applying at $x=x_{\star}$ we find

[TABLE]

(with the reasonable assumption $\tau^{\prime\prime}(0)=0$), which gives

[TABLE]

Numerically, we find $s_{\star}=0.12.$ We end up with the following formula:

[TABLE]

2.4.2 Stabilization

We run a numerical test on the stochastic model corresponding to stabilization, for $\tau_{0}=0.02$, and the other parameters as in Figure (1(a)). In this case, $x_{\star}$ correspond to the mean trait of the population for large time. From, (42) we find $x_{\star}=0.067$, and from (40), we obtain $\rho_{\star}=1.99$, which corresponds to what we can see on Figure (1(a)).

2.4.3 Threshold for cycles

Since equation (20) remains unchanged, we obtain the following threshold for cycles (corresponding to (22)):

[TABLE]

With our choice of parameters, we obtain $\tau_{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 $\tau_{0}$ is close to the threshold).

2.4.4 Threshold for extinction

Using (23), we can also find a threshold for extinction:

[TABLE]

For our choice of parameters, we obtain $\tau_{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 $\tau_{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 6(a) (for death rate: Figure 6(b)). Non-concerned parameters remain fixed as in Subsection 2.1.

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 $\tau$ 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 (1). 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 Billiard et al. (2018, 2015). 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 (1.2) 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 (7) 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 7. On Figures 7(a)-7(c) we trace the moment of the regrowth for different models. Figure 7(a) 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 $\frac{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 7(b) to 7(c)). 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 7(b) to 7(c)).

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 $\rho$ 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.