Camille Coron

Manon Costa
email: manon.costa@math.univ-toulouse.fr

Fabien Laroche
email: fabien.laroche@irstea.fr

Hélène Leman
email: helene.leman@inria.fr

Charline Smadi
email: charline.smadi@irstea.fr

Emergence of homogamy in a two-loci stochastic population model

Keywords: Birth and death processes with interactions, multitype branching processes, large population limits, mating preferences MSC 2000 subject classifications: 60J80, 60J27, 37N25, 92D25

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction and motivation

Assortative mating is a mating pattern in which individuals with similar phenotypes reproduce more frequently than expected under uniform random mating. Such a reproductive behaviour is widespread in natural populations and has an important role in the shape of their evolution (see for instance [START_REF] Mclain | Male choice, fighting ability, assortative mating and the intensity of sexual selection in the milkweed longhorn beetle, tetraopes tetraophthalmus (coleoptera, cerambycidae)[END_REF][START_REF] Herrero | Male and female synchrony and the regulation of mating in flowering plants[END_REF][START_REF] Savolainen | Sympatric speciation in palms on an oceanic island[END_REF] or the review [START_REF] Jiang | Assortative mating in animals[END_REF] on assortative mating in animals). In particular assortative mating is expected to be a driving force for speciation, which is the process by which several species arise from a single one [START_REF] Gregorius | A two-locus model of speciation[END_REF]. Here we ask the question of assortative mating emergence in a population: if at some point in time, one mutant starts mating preferentially with individuals of the same type, while the other individuals still choose their mate uniformly at random, can this mutant invade the population? A key feature to answer this question is how the assortative mating mutation affects the total reproduction rate of individuals. The existence of a preference for a given phenotype is often associated with a decay of reproductive success when mating with other phenotypes. As a consequence, if the proportion of preferred individuals is low in the population, the assortative mating mutation may be detrimental because it decreases the total reproduction success. Consequently, we expect that an assortative mating mutant will be able to invade only if its choosiness is compensated by an increased number of potential mates or if the advantage given by the preference is high enough.

In this work we aim at quantifying the conditions on the trade-off between advantage and cost for assortative mating and on the phenotype composition of the existing population needed for the mutation to invade the population. In order to reach this goal, we build a stochastic individual-based population model with varying size, which explores how the relationship between increase in the number of mates and mating bias towards individuals of the same type affects the long time number of individuals having mating preferences in the population.

The class of stochastic individual-based models with competition and varying population size we are extending have been introduced in the 90's in [START_REF] Bolker | Using moment equations to understand stochastically driven spatial pattern formation in ecological systems[END_REF][START_REF] Dieckmann | The dynamical theory of coevolution: a derivation from stochastic ecological processes[END_REF] and made rigorous in a probabilistic setting in the seminal paper of Fournier and Méléard [START_REF] Fournier | A microscopic probabilistic description of a locally regulated population and macroscopic approximations[END_REF]. Initially restricted to asexual populations, such models have evolved to incorporate the case of sexual reproduction, in both haploid [START_REF] Smadi | An eco-evolutionary approach of adaptation and recombination in a large population of varying size[END_REF][START_REF] Leman | A stochastic model for reproductive isolation under asymmetrical mating preferences[END_REF] and diploid [START_REF] Collet | A rigorous model study of the adaptive dynamics of mendelian diploids[END_REF][START_REF] Coron | Slow-fast stochastic diffusion dynamics and quasi-stationary distributions for diploid populations[END_REF][START_REF] Neukirch | Survival of a recessive allele in a mendelian diploid model[END_REF][START_REF] Smadi | Looking for the right mate in diploid species: How does genetic dominance affect the spatial differentiation of a sexual trait?[END_REF] populations. Taking into account varying population sizes and stochasticity is necessary if we aim at better understanding phenomena involving small populations like invasion of a mutant population [START_REF] Champagnat | A microscopic interpretation for adaptive dynamics trait substitution sequence models[END_REF] or population extinction time. Assuming that individuals initially have no preference and choose their mate uniformly at random, we suppose that a mutation arises in the population: individuals carrying the mutation (denoted P ) have a higher (resp. smaller) reproductive success when mating with individuals of the same (resp. different) phenotype than individuals without the mutation. We study under which conditions on the parameters (birth and death rates, competition, mutational effects, initial population state, ...) the mutation P has a positive probability to invade the population, and how to identify this probability. We also characterize the time needed for the mutation to get fixed in the population when it happens. Finally, we provide the invasion dynamics as well as the final population state, when the mutation gets fixed.

In order to obtain our results, we study the population process at two different scales. When one sub-population is of small size the stochasticity of its size has a major effect on the population long time behaviour, and we study its dynamics on N := {0, 1, 2, ...}. This is for example the case of the mutant population when it arises. When on the contrary all sub-populations sizes are large, we approximate the stochastic process by a mean field limit which is a dynamical system.

Model and main results

We consider a population of individuals that reproduce sexually and compete with each other for a common resource. Individuals are haploid and are characterized by their genotype at two loci located on different chromosomes. Locus 1 presents two alleles, denoted by A and a, and codes for phenotypes. Locus 2 presents two alleles denoted by P and p, and codes for assortative mating, which is defined relatively to the first locus (similar models were introduced in Biology, see for example [START_REF] Gregorius | A two-locus model of speciation[END_REF]). More precisely, we assume that all individuals try to reproduce at the same rate. To this aim, they choose a mate, uniformly at random among the other individuals of the population. Next, individuals carrying allele p reproduce indifferently with their chosen partner, while individuals carrying allele P reproduce with a higher probability with individuals carrying the same allele at locus 1. Note that reproduction is not completely symmetric: only the genotype of the individual initiating the reproduction determines the presence or not of assortative mating.

The genotype of each individual belongs to the set G := {AP, Ap, aP, ap} and the state of the population is characterized at each time t by a vector in N 4 giving the respective numbers of individuals carrying each of these four genotypes. The dynamics of this population is modeled by a multi-type birth-and-death process (N (t), t ≥ 0) := (N AP (t), N Ap (t), N aP (t), N ap (t), t ≥ 0) with values in N 4 , integrating competition, Mendelian reproduction and assortative mating. More precisely, when the population is in state n = (n AP , n Ap , n aP , n ap ) ∈ N 4 with size n = n AP + n Ap + n aP + n ap , then the rate at which the population looses an individual with genotype i ∈ G, is equal to

d i (n) = n i d + c K n .
(2.1)

The parameters d ∈ R + , c > 0 and K > 0 respectively model the natural and the competition death rates of individuals and a scaling parameter of the total population size. This parameter quantifies the environment's carrying capacity, which is a measure of the maximal population size that the environment can sustain for a long time. In the sequel we will be interested in the behaviour of the system for large but finite K.

When the population is in state n, the rate b i (n) at which an individual with genotype i ∈ G is born, is defined by 1 in Appendix B summarizes the different rates at which a pair of parents with given genotypes gives birth to an offspring with a given genotype. This explains how the birth rates (2.2) are obtained.

b AP (n) = b n AP + 1 n β 1 n AP n AP + n Ap 2 -β 2 n AP n aP + n ap 4 + n Ap n aP 4 + ∆ aP 2n b Ap (n) = b n Ap + 1 n β 1 n Ap n AP 2 -β 2 n Ap n aP 4 + n AP n ap 4 - ∆ aP 2n b aP (n) = b n aP + 1 n β 1 n aP n aP + n ap 2 -β 2 n aP n AP + n Ap 4 + n ap n AP 4 - ∆ aP 2n b ap (n) = b n ap + 1 n β 1 n ap n aP 2 -
Throughout the paper, we will make the following assumptions on the parameters:

(1) b > d (2) β 1 ≥ 0 (3) 0 ≤ β 2 ≤ 1
The first assumption ensures that a population of individuals mating uniformly at random is not doomed to a rapid extinction because of a natural death rate larger than the birth rate under uniform random mating. The second (resp. third) assumption means that choosy individuals have a higher (resp. smaller) probability to give birth when mating with an individual with the same (resp. different) trait (A or a).

We assume that at time 0, all individuals mate uniformly at random (no sexual preference, all individuals carry allele p), and that the population size is close to its long time equilibrium, (b -d)K/c (see on page 5 for details). A mutant (or a migrant) appears in the population, with genotype αP , where α ∈ A := {A, a}. The goal of our main theorem (Theorem 1) is to study a step in Darwinian evolution, that consists in the progressive invasion of the new allele P and loss of initial allele p in the population. The proof of this theorem relies on the study of three phases in the population dynamics trajectories (mutant survival or extinction, mean-field phase, and resident allele extinction) that are respectively defined and studied in Subsections 3.1, 3.2, 3.3. The statement of Theorem 1 requires the introduction of several quantities that we define now.

Our first goal is to determine conditions under which the mutant population has a positive probability to survive and invade the resident population. In order to answer this question, we will compare the mutant population with a branching process during the first times of the invasion. This comparison follows from the following observation that will be proved in Proposition 3.1: as long as the mutant population size is negligible with respect to the carrying capacity K, the dynamics of the resident population will not be affected by the presence of the mutants and will stay close to its initial state. In other words, the size and proportions of the resident population will remain almost constant and the dynamics of the mutant population will be close to the dynamics of the process N = ( NA , Na ), which is a bi-type branching process with the following transition rates:

( NA , Na ) → ( NA + 1, Na ) at rate βAA NA + βaA Na ( NA , Na ) → ( NA , Na + 1) at rate βAa NA + βaa Na ( NA , Na ) → ( NA -1, Na ) at rate b NA ( NA , Na ) → ( NA , Na -1) at rate b Na , (2.3) where for α ∈ A, ᾱ = A \ α, βαα := b 2 1 + (β 1 + 1)ρ α - β 2 2 ρ ᾱ , βαᾱ := b 2 1 - β 2 2 ρ ᾱ, (2.4)
and

ρ A := N Ap (0) N p (0) and ρ a := 1 -ρ A = N ap (0) N p (0) (2.5)
are the initial proportions in the resident population. The rates of this branching process have been obtained by considering the dynamics of (N AP , N aP ) described by (2.1) and (2.2)

when (N Ap , N ap ) = (Kρ A b-d c , K(1 -ρ A ) b-d c ), N = K b-d c
and the second order terms in N AP and N aP are neglected. We denote the extinction probabilities of the process N by

q α := P(∃t < ∞, N(t) = 0| N(0) = e α ), (2.6)
α ∈ A, e A = (1, 0) and e a = (0, 1), meaning that the process starts with only one individual of type A or a. Classical results of branching process theory (see [START_REF] Athreya | Branching processes[END_REF]) ensure that these extinction probabilities correspond to the smallest solution to the system of equations

u A (s A , s a ) := b(1 -s A ) + βAA (s 2 A -s A ) + βAa (s A s a -s A ) = 0 (2.7) u a (s A , s a ) := b(1 -s a ) + βaa (s 2 a -s a ) + βaA (s A s a -s a ) = 0.
Moreover, the process N is supercritical (and in this case q A and q a are not equal to one) if and only if the following matrix, which corresponds to the infinitesimal generator of the branching process, has a positive eigenvalue

J := βAA -b βAa βaA βaa -b (2.8)
that is to say if and only if

β 1 > β 2 or ρ A (1 -ρ A ) < β 1 (β 2 + 2) 2(β 1 + β 2 )(β 1 + 2)
(2.9)

(see the proof of Proposition 3.3). We denote by λ the maximal eigenvalue of (2.8), which is thus positive when (2.9) holds and which will be of interest to quantify the time before invasion. Notice that J can be written as b times a matrix only depending on (ρ A , β 1 , β 2 ).

As a consequence, λ can be written λ = b λ(ρ A , β 1 , β 2 ). We will use this notation in Theorem 1 to make appear the dependence on the parameters, and use λ elsewhere for the sake of readability.

If the mutant population invades and its size reaches order K with K large, the population dynamics enters a second phase during which it is well approximated (see Proposition 2.1 for a rigorous statement) by a mean field process. More precisely, if we define the rescaled process

(Z K (t), t ≥ 0) := N AP (t) K , N Ap (t) K , N aP (t) K , N ap (t) K , t ≥ 0 ,
then it will be close to the solution of the dynamical system

żi = b i (z) -(d + cz)z i , i ∈ G, (2.10)
where z = z AP + z Ap + z aP + z ap is the total size of the population and the functions b i have been defined in Equation (2.2). This dynamical system has a unique solution, as the vector field is locally Lipschitz and that the solutions do not explode in finite time [START_REF] Chicone | Ordinary Differential Equations with Applications[END_REF]. If we denote by (z

(z 0 ) (t), t ≥ 0) = (z AP (t), z Ap (t), z aP (t), z ap (t), t ≥ 0)
this unique solution starting from z(0) = z 0 ∈ R 4 + , we have the following result, which derives from Theorem 2.1 p 456 in [START_REF] Ethier | Markov processes. Characterization and convergence[END_REF].

Lemma 2.1. Let T ∈ R * + .
Assume that the sequence (Z K (0), K ≥ 1) converges in probability to some deterministic vector z 0 = (z AP (0), z Ap (0), z aP (0), z ap (0)) ∈ R 4 + when K goes to infinity. Then lim

K→∞ sup s≤T ||Z K (s) -z (z 0 ) (s)|| ∞ = 0 in probability, where || • || ∞ denotes the L ∞ -Norm in R 4 .
Notice that when there are only individuals of type p in the population (no sexual preferences), the dynamical system (2.10) is

żAp = z Ap (b -d -c(z Ap + z ap )) żap = z ap (b -d -c(z Ap + z ap ))
This system admits an infinity of equilibria:

• (z Ap , z ap ) = (0, 0), which is unstable

• (z Ap , z ap ) = (ρ(b -d)/c, (1 -ρ)(b -d)/c
) for all ρ ∈ [0, 1], which are non hyperbolic. However, if we consider the equation giving the dynamics of the total population size z = z Ap + z ap , we get ż = z(b -d -cz). Its solution, with a positive initial condition, converges to its unique stable equilibrium, (b -d)/c. That is why we will assume that the initial population size, before the arrival of the mutant, is (b -d)K/c.

A fine study of the dynamics of the solutions to (2.10) with our particular initial conditions, that is to say few individuals mating assortatively at the beginning and a majority of A (or a) in both resident and mutant populations (see Section 3.2.2), will allow us to show that the dynamical system converges to an equilibrium where some of the variables z i , i ∈ G are equal to 0. When the population size of these i becomes too small (of order smaller than K before rescaling), the mean fields approximation stops being a good approximation, and we will again compare the dynamics of the small population sizes with these of branching processes (now subcritical). The birth and death rates of these branching processes will provide the time to extinction of these small populations (see Section 3.3).

Combining all these steps, we are able to describe the invasion/extinction dynamics of the mutant population, which is the subject of the main result of this paper, Theorem 1. Before stating it, we need to introduce some last notations: a set of interest for the rescaled process Z K , for any µ > 0

S µ := b(1 + β 1 ) -d c -µ, b(1 + β 1 ) -d c + µ × {0} × {0} × {0}, (2.11)
a stopping time describing the time at which Z K reaches this set,

T Sµ := inf{t ≥ 0, Z K (t) ∈ S µ }.
(2.12) as well as a stopping time which gives the first time when the rescaled P -mutant population size reaches any threshold (from below or above): for any ε ≥ 0,

T P ε := inf {t > 0, N P (t) = εK } , (2.13)
where x is the integer part of x.

Theorem 1. Assume that λ = 0,

Z K Ap (0), Z K ap (0) → K→∞ ρ A b -d c , (1 -ρ A ) b -d c
in probability with ρ A > 1/2 and that

(N αP (0), N ᾱP (0)) = (1, 0), α ∈ A.
Then there exists a Bernoulli random variable with parameter 1 -q α , B, such that for any

0 < µ < (b(1 + β 1 ) -d)/c: lim K→∞ T Sµ ∧ T P 0 ln K , 1 {T Sµ <T P 0 } = B × 1 b λ(ρ A , β 1 , β 2 ) + 2 bβ 1 , 1 , (2.14)
where the convergence holds in probability. Moreover,

1 {T P 0 <T Sµ } N(T P 0 ) K -(0, ρ A , 0, 1 -ρ A ) b -d c 1 -→ K→∞ 0 in probability, (2.15)
where (a, b)

1 = |a| + |b|.
Notice that if condition (2.9) does not hold, q α = 1, and the convergence in (2.14) is an almost sure convergence to (0, 0) meaning that the mutant population dies out in a time smaller than ln(K). In this case, the allelic proportions in the resident population do not vary. Condition (2.9) gives two possible sufficient conditions for the mutant population to invade with positive probability. The first one imposes that the trade-off between the advantage for homogamous reproduction (β 1 ) and the loss for heterogamous reproduction (β 2 ) has to be favourable enough. The second condition requires a low level of initial allelic diversity at locus 1 (alleles A and a). In particular, even if the advantage for homogamy is very low, very asymmetrical initial conditions (ρ A close to 0 or 1) will ensure the invasion of the mutation with positive probability. As expected, these conditions are the same as the conditions for the approximating branching process N defined on page 4 to be supercritical. In fact, as we will see later in the proof, the random variable B will be the indicator of survival of a version of N coupled with the mutant process.

Let us emphasize that our result ensures that when the mutant population invades, then the final population is monomorphic, and all individuals carry the allele a or A which was majority in the resident p-population. Only the mutant invasion probability depends on the allele carried by the first P individual.

We were not able to obtain an explicit formula in general for the extinction probability q α of the assortative mating mutation, solutions of (2.7). However, in the particular case when there are only A or a-individuals in the population before the arrival of the mutant, we can derive the invasion probability (see the proof in Section A.2) Proposition 2.1. Assume that there are only A individuals before the arrival of the mutant (ρ A = 1). In this case,

q A = 2 2 + β 1 and q a = 1 2 -β 2   6 -β 1 β 2 + 4β 1 -β 2 2 + β 1 - 6 -β 1 β 2 + 4β 1 -β 2 2 + β 1 2 -4(2 -β 2 )   .
Results obtained with the help of the software Mathematica show a complex dependency with respect to parameters. We performed numerical simulations of the extinction probabilities (q A , q a ) using Newton approximation scheme starting from (0, 0). We computed the values of q A as a function of ρ A for different values of β 1 and β 2 . Using the symmetry of our model, we have that q a (ρ A ) = q A (1 -ρ A ). We observe on Figure 1 that q A is a continuous function of ρ A but that it is not differentiable near criticality. On the left, β 2 is fixed to 0.7 and β 1 varies. On the right β 1 is fixed to 0.2 and β 1 varies. In both cases b = 1.

Remark 1. We assumed that the initial population state is close to the equilibrium state of the population when all individuals mate uniformly at random, because any neighbourhood of such an equilibrium is reached within a finite time as soon as the initial population size is of order K. We thus could relax this assumption and only assume that the p-population size is of order K and N Ap (0) > N ap (0). This would however require more complex notations.

The rest of the paper is devoted to the proof of Theorem 1. Notice that for the sake of readability, we will not indicate anymore the dependency of the rescaled process Z K on K and will instead write Z.

Proof of Theorem 1

3.1. Probability and time of the mutant invasion. The first step of the proof of Theorem 1 consists in studying the population dynamics when a mutant of type P appears in a well-established population of types ap and Ap. We would like to know under which conditions on the parameters the mutant population may invade the resident population and what is the probability that the invasion happens.

We will show in particular that when the mutant appears and as long as the mutant population size is negligible with respect to the carrying capacity K, its dynamics is close to the dynamics of the process N, which has been introduced in Section 1. Next, as long as the P -population size is small compared with K, that is, as long as its dynamics is close to the one of N, we prove that the p-population size and the proportion of A-individuals in the p-population will not vary considerably from their initial values.

In order to state rigorously these results, let us recall definition (2.13) and introduce two more stopping times. The first one gives the first time when the genetic proportions in the p-population deviate considerably from their starting values: for any ε > 0,

U ε := inf t ≥ 0, N Ap (t) N p (t) - N Ap (0) N p (0) > ε . (3.1)
The second one concerns the total p-population size: for any ε > 0,

R ε := inf t ≥ 0, N p (t) K - b -d c > ε . (3.2)
Note that these stopping times depend on the scaling parameter K. However, to avoid cumbersome notations, we drop the K dependency. We recall that q α is the extinction probability of this process and that λ is the principal eigenvalue of the matrix J defined in (2.8), which can be rewritten

J = b 2   ρ A β 1 -(1 -ρ A ) β 2 2 + 1 ρ A 1 -β 2 2 (1 -ρ A ) 1 -β 2 2 (1 -ρ A )β 1 -a 1 β 2 2 + 1   . (3.3)
The main result along the route of proving Theorem 1 can now be stated. It ensures that the probability that a mutant P generates a P -population whose size reaches the order K is close to 1-q α (which is the survival probability of the process N starting from an α-individual and has been defined in (2.6)), whereas its probability of extinction is close to q α . Moreover, the invasion or extinction of the mutant population occurs before the resident population size deviates substantially from its equilibrium, and the time of invasion is approximately log(K)/λ. Proposition 3.1. Let α be in A and assume that the initial condition satisfies

N P (0) = N αP (0) = 1, N p (0) = b -d c K and ρ A = N Ap (0) N p (0) ,
and moreover that λ = 0,

where λ is the principal eigenvalue of matrix (2.8). There exist a function η going to 0 at 0 and a positive constant A 0 such that for any ξ ∈ {1/2, 1}, lim sup

K→∞ P T P ε ξ < T P 0 ∧ R A 0 ε ∧ U ε 1/6 , T P ε ξ ln K - 1 λ ≤ η(ε) N P (0) = e α -(1 -q α ) = o ε (1),
and lim sup

K→∞ P T P 0 < T P ε ξ ∧ R A 0 ε ∧ U ε 1/6 N P (0) = e α -q α = o ε (1), (3.5)
where by convention, o ε (1) goes to 0 when ε goes to 0.

Remark 2. This proposition accounts for two opposite behaviours of the mutant process. Indeed Assumption (3.4) ensures that either the process N is supercritical (λ > 0 under condition (2.9)) and q α ∈ (0, 1), or the process N is subcritical and q α = 1.

The end of this section will be devoted to the proof of Proposition 3.1, which will be divided into three steps.

3.1.1. Control of the proportions in the resident population. We will first prove that the proportions in the resident population do not vary substantially before the mutant population goes extinct or invades. More precisely, we have the following lemma. Lemma 3.1. Suppose that the assumptions of Proposition 3.1 hold. For any A > 0, there exist ε 0 such that for any ξ ∈ {1/2, 1} and ε ≤ ε 0 , lim sup

K→∞ P U ε 1/6 < R Aε ∧ T P ε ξ ∧ T P 0 ≤ C(A, ξ)ε 1/12 ,
where C(A, ξ) is a positive constant.

Proof. The statement of Lemma 3.1 is a direct consequence of the following inequality:

lim sup K→∞ P   sup t≤U ε 1/8 ∧R Aε ∧T P ε ξ ∧T P 0 N Ap (t) N p (t) - N Ap (0) N p (0) > ε 1/6   ≤ Cε 1/12 . (3.6)
To prove (3.6), we decompose the process N Ap /N p as the sum of a square integrable martingale M p and of a finite variation process V p (see (3.10) and (3.11) for their expressions).

Using such a decomposition and introducing for the sake of readability the notation

τ ε := U ε 1/8 ∧ R Aε ∧ T P ε ξ ∧ T P 0 , (3.7)
we find that for ε small enough,

P sup t≤τε N Ap (t) N p (t) - N Ap (0) N p (0) > ε 1/6 ≤ P sup t≤τε |M p (t)| > ε 1/6 2 + P sup t≤τε |V p (t)| > ε 1/6 2 ≤ 2 ε 1/6 E [|M p (τ ε )|] + √ 2 ε 1/12 E sup t≤τε |V p (t)| ≤ 2 ε 1/6 E M 2 p (τ ε ) + E sup t≤τε |V p (t)| , (3.8)
where we applied Doob maximal, Markov, Cauchy-Schwarz and Jensen inequalities. Hence, it remains to bound the two last expectations of (3.8). In the vein of Fournier and Méléard [START_REF] Fournier | A microscopic probabilistic description of a locally regulated population and macroscopic approximations[END_REF] we represent the population process in terms of Poisson measures.

Let (Q

(ρ) αp (ds, dθ), α ∈ A, ρ ∈ {b, d}) be four independent Poisson random measures on R 2
+ with intensity dsdθ representing respectively the birth and death events of Ap and ap individuals. That is, for any α ∈ A, the p-population size processes can be written

N αp (t) = N αp (0) + t 0 R + 1 {θ≤bαp(N(s-))} Q (b) αp (ds, dθ) -1 {θ≤dαp(N(s-))} Q (d) αp (ds, dθ) (3.9)
where the quantities b αp and d αp have been defined in (2.1) and (2.2). Let us also denote by

Q( ) αp (ds, dθ) := Q ( )
αp (ds, dθ) -dsdθ the associated compensated measure, for any ∈ {b, d}, α ∈ A. From (3.9), we find, for t ≥ 0,

N Ap (t) N p (t) = N Ap (0) N p (0) + M p (t) + V p (t),
with M p and V p such that:

M p (t) = t 0 R + 1 {θ≤b Ap (N(s-))} N ap (s-) N p (s-)(N p (s-) + 1) Q(b) Ap (ds, dθ) (3.10) - t 0 R + 1 {θ≤d Ap (N(s-))} N ap (s-) N p (s-)(N p (s-) -1) Q(d)
Ap (ds, dθ)

- t 0 R + 1 {θ≤bap(N(s-))} N Ap (s-) N p (s-)(N p (s-) + 1) Q(b)
ap (ds, dθ)

+ t 0 R + 1 {θ≤dap(N(s-))} N Ap (s-) N p (s-)(N p (s-) -1) Q(d)
ap (ds, dθ), and

V p (t) = t 0 b Ap (N(s)) N ap (s) N p (s)(N p (s) + 1) -d Ap (N(s)) N ap (s) N p (s)(N p (s) -1) -b ap (N(s)) N Ap (s) N p (s)(N p (s) + 1) + d ap (N(s)) N Ap (s) N p (s)(N p (s) -1) ds = t 0 b Ap (N(s))N ap (s) -b ap (N(s))N Ap (s) ds N p (s)(N p (s) + 1) . (3.11)
Using Equation (2.2), we obtain the existence of a finite constant C such that

|b Ap (N )N ap -b ap (N )N Ap | ≤ C N P N 2 p N P + N p .
Hence,

sup t≤τε |V p (t)| ≤ C τε 0 N P (s) N P (s) + N p (s)
ds.

(3.12)

This will help us bounding the last term of inequality (3.8). On the other hand, to deal with the penultimate term in (3.8), we use the quadratic variation of the martingale M p which equals

M p τε = τε 0 b Ap (N(s)) N 2 ap (s) N 2 p (s)(N p (s) + 1) 2 ds + τε 0 d Ap (N(s)) N 2 ap (s) N 2 p (s)(N p (s) -1) 2 ds (3.13) + τε 0 b ap (N(s)) N 2 Ap (s) N 2 p (s)(N p (s) + 1) 2 ds + τε 0 d ap (N(s)) N 2 Ap (s) N 2 p (s)(N p (s) -1) 2 ds. = τε 0 b Ap (N(s))N 2 ap (s) + b ap (N(s))N 2 Ap (s) N 2 p (s)(N p (s) + 1) 2 + d Ap (N(s))N 2 ap (s) + d ap (N(s))N 2 Ap (s) N 2 p (s)(N p (s) -1) 2 ds
To handle the first term, let us remark that b Ap (N) and b ap (N) can be bounded from above by CN p for a positive constant C. Therefore

b Ap (N)N 2 ap + b ap (N)N 2 Ap N 2 p (N p + 1) 2 ≤ C N p .
For the second term we have

d Ap (N)N 2 ap + d ap (N)N 2 Ap N 2 p (N p -1) 2 ≤ (d + cN/K)N ap N Ap N p (N p -1) 2 . Since, for N p ≥ 2, N ap N Ap N p (N p -1) 2 ≤ 4 N p ,
we obtain that, if C and K are sufficiently large,

M p τε ≤ τε 0 4 N p (s) C + d + c K (N P (s) + N p (s)) ds ≤ C τε 0 1 N p (s) ds ≤ C τε 0 N P (s) N p (s) ds.
(3.14)

From (3.12) and (3.14), we get that there exists a finite C such that

E M 2 p (τ ε ) + E sup t≤τε |V p (t)| ≤ C K E τε 0 N P (s)ds . (3.15)
In view of (3.8), the problem is thus reduced to show the following property:

ε -1/6 C K E τε 0 N P (s)ds ≤ C ε 1/12 ,
for a finite C , or equivalently,

E τε 0 N P (s)ds ≤ C Kε 1/2 . (3.16)
To this aim, we will prove that there exist two real numbers γ 1 and γ 2 such that the function

f on N 4 defined by f (N) := γ 1 N AP + γ 2 N aP , (3.17)
satisfies that there exists ε sufficiently small such that for any t ≤ τ ε (recall equation (3.7)),

Lf (N(t)) ≥ N P (t). (3.18)
Here L is the infinitesimal generator of N. Indeed, if (3.18) holds, it will imply that

C K E τε 0 N P (s)ds ≤ C K E τε 0 Lf (N(s))ds = C K E [f (N(τ ε )) -f (N(0))] ≤ C K max{γ 1 , γ 2 }ε ξ K -min{γ 1 , γ 2 } , (3.19)
which is sufficient to obtain (3.16), whatever the signs of γ 1 and γ 2 .

The last step of the proof consists in proving the existence of γ 1 and γ 2 satisfying (3.17) and (3.18). Let us now apply the infinitesimal generator of N to the function f defined in (3.17):

Lf (N(t)) = γ 1 [b AP (N) -d AP (N)] + γ 2 [b aP (N) -d aP (N)] = N AP (t) γ 1 β (P ) AA (t) -δ(t) + γ 2 β (P ) Aa (t) + N aP (t) γ 1 β (P ) aA (t) + γ 2 β (P ) aa (t) -δ(t) , (3.20)
where δ(t) = d + cN (t)/K and for α ∈ A,

β (P ) αα (t) = b 1 + β 1 2 2N αP (t) + N αp (t) N P (t) + N p (t) - β 2 4 4N ᾱP (t) + N ᾱp (t) N P (t) + N p (t) - 1 2
N ᾱp (t) N P (t) + N p (t) (3.21)
and

β (P ) α ᾱ (t) = b 2 1 - β 2 2 N ᾱp (t) N P (t) + N p (t) . (3.22)
Then we see that to obtain (3.18) it is enough to choose Γ = (γ 1 , γ 2 ) such that ∀t ≤ τ ε

J (P ) (t)Γ T := β (P ) AA (t) -δ(t) β (P ) Aa (t) β (P ) aA (t) β (P ) aa (t) -δ(t) γ 1 γ 2 > 1 1
where the inequality is applied to each coordinate. Note that J (P ) (t) is not easy to study. We will thus approximate this matrix by a simpler one as soon as t ≤ τ ε . More precisely, we will prove that there exists a constant C such that for every t ≤ τ ε ,

|J (P ) (t) -J| ij ≤ Cε 1/8 , i, j ∈ {1, 2} (3.23)
where the matrix J given in (2.8) is the infinitesimal generator of the branching process N (defined in (2.3)) which approximates (N AP , N aP ) near the equilibrium of the resident population.

First, as t ≤ τ ε ≤ T P ε ξ K ∧ R Aε , we have b -d -c N P + N p K ≤ c b -d c - N p K + c N P K ≤ cAε + cε ξ . (3.24)
Secondly, using also that t ≤ U ε 1/8 and recalling that ρ A = N Ap (0)/N p (0), we find that for ε small enough,

2N AP + N Ap N P + N p -ρ A ≤ 2N AP N P + N p + ρ A N p N P + N p -1 + N p N P + N p N Ap N p -ρ A ≤ (2 + ρ A ) N P N P + N p + ε 1/8 ≤ 3 ε ξ c b -d -cAε + ε 1/8 ≤ C 1 ε 1/8 , (3.25)
where C 1 is a finite constant and ξ ∈ {1/2, 1}. Similarly, we prove that, if C 1 is sufficiently large,

4N aP + N ap N -(1 -ρ A ) ≤ C 1 ε 1/8 and N ap N P + N p -(1 -ρ A ) ≤ C 1 ε 1/8 . ( 3
β (P ) AA -δ -( βAP -b) ≤ b -d -c N P + N p K + bβ 1 2 2N AP + N Ap N P + N p -ρ A + β 2 4 4N aP + N ap N P + N p -(1 -ρ A ) + 1 2 N ap N P + N p -(1 -ρ A ) ≤ C 2 ε 1/8 , and
β (P ) Aa -βAa ≤ b(2 -β 2 ) 4 N ap N P + N p -(1 -ρ A ) ≤ C 2 ε 1/8 .
The last terms in (3.20) can be bounded using similar computations, which yields (3.23).

Let us finally choose (γ 1 , γ 2 ). Recall the definition of J in (2.8). In particular, the matrix J + bId has positive coefficients and we can apply Perron-Frobenius Theorem: J + bId possesses a positive right eigenvector Γ = (γ 1 , γ2 ) associated to the positive principal eigenvalue λ + b and thus (J + bId) ΓT = (λ + b) ΓT , and

J ΓT = λ ΓT .
Since both coordinates of Γ are positive and λ = 0 , we can define Γ = 2 Γ(λ(γ 1 ∧ γ2 )) -1 = (γ 1 , γ 2 ). It is solution to

JΓ T = λΓ T where λγ i ≥ 2, ∀i ∈ {1, 2}.
Combining with (3.23), we deduce

γ 1 β (P ) AA (t) -δ + γ 2 β (P ) Aa (t) -λγ 1 = γ 1 β (P ) AA (t) -δ -( βAA -b) + γ 2 β (P ) Aa (t) -βAa ≤ (|γ 1 | + |γ 2 |) C 2 ε 1/8 .
Finally, as λγ 1 ≥ 2, if ε is sufficiently small, for any t ≤ τ ε ,

γ 1 β (P ) AA (t) -δ(t) + γ 2 β (P )
Aa (t) ≥ 1, which leads to (3.18) and ends the proof of Lemma 3.1 using similar computations for the second term.

3.1.2. Control of the resident population size. Lemma 3.1 ensures that the proportions of types A and a in the p-population stay almost constant during the time interval under consideration. We now prove that it is also the case for the total p-population size. Lemma 3.2. Under the assumptions of Proposition 3.1 there exist two finite constants A 0 and ε 0 such that for any ξ ∈ {1/2, 1} and ε ≤ ε 0 , lim sup

K→∞ P R A 0 ε < U ε 1/6 ∧ T P ε ξ ∧ T P 0 = 0.
Proof. Recall that Z = N/K. As long as t ≤ T P 0 ∧ T P ε ξ , we couple the process Z p , which describes the total p-population size dynamics, with two birth and death processes, Z 1 p and Z 2 p such that Z 1 p (t) ≤ Z p (t) ≤ Z 2 p (t), a.s. ∀t ≤ T P 0 ∧ T P ε ξ . To this aim, we use bounds on the birth and death rates of Z p . Once again, everything depends on K, but for the sake of readability, we drop the K dependency. Since β 2 ≤ 1, processes Z 1 p and Z 2 p may be chosen with the following birth and death rates

Z 1 p : i K → i+1 K at rate K b i K -bε ξ i K → i-1 K at rate K i K d + cε ξ + c i K and Z 2 p : i K → i+1 K at rate K b i K + bε ξ (β 1 +1) 2 i K → i-1 K at rate K i K d + c i K .
We will first prove that processes Z 1 p and Z 2 p stay close to the value ζ := (b -d)/c for at least an exponential (in K) time with a probability close to one when K is large. To this aim, we will study the following stopping times The differential equation (3.28) admits two positive equilibria:

R i η := inf t ≥ 0, Z i p ∈ [ζ -η, ζ + η] , (3.27
ζ 1,± (ε) := b -d -cε ξ ± (b -d -cε ξ ) 2 -4bε ξ c 2c .
A direct analysis of the sign of z(b -d -cε ξ -cz) -bε ξ shows that for any fixed ε > 0, any solution with initial condition on ]ζ 1,-(ε), +∞[, converges to the stable equilibrium ζ 1,+ (ε) when t goes to infinity. Moreover, using

( √ a + √ b)( √ a - √ b) = a -b (3.29) yields, if ε is sufficiently small (such that ζ 1,+ (ε) ≥ (b -d)/2c), ζ 1,-(ε) = bε ξ cζ 1,+ (ε) ≤ 2b b -d ε ξ ≤ 2ε ξ 0 .
Then, we can find two constants A 0 and ε 0 such that, for any ε ≤ ε 0 ,

|ζ 1,+ (ε) -ζ| ≤ (A 0 -1)ε ξ and 2ε ξ 0 ∈ [ζ -A 0 ε ξ , ζ + A 0 ε ξ ].
Moreover, using a reasoning similar to the one in the proof of Theorem 3(c) in [START_REF] Champagnat | A microscopic interpretation for adaptive dynamics trait substitution sequence models[END_REF] (see also Proposition 4.1 in [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]), we construct a family (over K) of Markov jump processes Z 1 p whose transition rates are positive, bounded, Lipschitz and uniformly bounded away from 0, and for which the following estimate holds (Chapter 5 of Freidlin and Wentzell [START_REF] Freidlin | Random perturbations of dynamical systems[END_REF]): there exists V > 0 such that,

P(R 1 A 0 ε > e KV ) = P( R1 A 0 ε > e KV ) → K→+∞ 1, (3.30)
where R1 η is defined similarly as R 1 η but for the process Z 1 . Using a similar reasoning for process Z 2 p and if ε and V are small enough and A 0 is large enough, we have that

P(R 2 A 0 ε > e KV ) → K→+∞ 1. (3.31)
Finally, note firstly that

R A 0 ε ≥ R 1 A 0 ε ∧ R 2 A 0 ε on the set {R A 0 ε ≤ T P 0 ∧ T P ε ξ }.
In addition with (3.30) and (3.31), we deduce that

P(R A 0 ε ≤ e KV , R A 0 ε ≤ T P 0 ∧ T P ε ξ ) → K→+∞ 0.
Secondly, for ε small enough,

P(R A 0 ε ≤ T P 0 ∧ T P ε ξ ∧ U ε 1/6 ) ≤ P(R A 0 ε ≤ e KV , R A 0 ε ≤ T P 0 ∧ T P ε ξ ∧ U ε 1/6 ) + P(R A 0 ε ∧ T P 0 ∧ T P ε ξ ∧ U ε 1/6 ≥ e KV ) ≤ P(R A 0 ε ≤ e KV , R A 0 ε ≤ T P 0 ∧ T P ε ξ ) + e -KV E R A 0 ε ∧ T P ε ξ ∧ T P 0 ∧ U ε 1/8
. Thirdly, Equation (3.19) implies that for a finite C and ε small enough,

E R A 0 ε ∧ T P ε ξ ∧ T P 0 ∧ U ε 1/8 ≤ E U ε 1/8 ∧R A 0 ε ∧T P ε ξ ∧T P 0 0 N P (s)ds ≤ C(ε ξ K + 1).
The three last equations imply the statement of Lemma 3.2, which ends its proof.

3.1.3. Proof of Proposition 3.1. Lemmas 3.1 and 3.2 give us a control on the p-population size and the proportions of A and a individuals in this population. It will allow us to approximate the mutant population size by a bitype branching process at the beginning of the invasion process. We will assume along the proof that N P (0) = N αP (0) = 1, with α ∈ A, but we drop the conditioning notation for the sake of readability. Combining Lemmas 3.1 and 3.2, we obtain that the p-population size and the genotypic proportions in the p-population stay almost constant as long as the P -mutant population size is small. More precisely, if (3.4) is satisfied, there exist two constants A 0 and ε 0 such that for any ξ ∈ {1/2, 1} and ε ≤ ε 0 , lim inf

K→∞ P T P ε ξ ∧ T P 0 < R A 0 ε ∧ U ε 1/6 ≥ 1 -C(A 0 , ξ)ε 1/12 , (3.32)
where C(A 0 , ξ) is a positive constant. Hence, in what follows, we study the process on the event

Σ ε = T P ε ξ ∧ T P 0 < R A 0 ε ∧ U ε 1/6
, which has a probability close to 1. On this event, the death or invasion of the mutant population will occur before that the p-population deviates substantially from its initial composition. Thus, we can study the mutant population dynamics by approximating the resident population dynamics with a constant dynamics. More precisely, we couple the process (N AP , N aP ) on Σ ε with two multitype birth and death processes N (ε,-) and N (ε,+) with values in N 2 such that almost surely, for any t

≤ T P ε ξ ∧ T P 0 ∧ R A 0 ε ∧ U ε 1/6 and α ∈ A, N (ε,-) α (t) ≤ Nα (t) ≤ N (ε,+) α (t), N (ε,-) α (t) ≤ N αP (t) ≤ N (ε,+) α (t).
(3.33)

For * ∈ {+, -}, the process N (ε, * ) may be chosen with the rates:

N (ε, * ) → N (ε, * ) + e α at rate β (ε, * ) Aα N (ε, * ) A + β (ε, * ) aα N (ε, * ) a N (ε, * ) → N (ε, * ) -e α at rate δ (ε, * ) α N (ε, * ) α
where ) are increasing (resp. decreasing) when ε increases. Let us denote for * ∈ {-, +} and α ∈ A by q (ε, * ) α the extinction probability of the process N (ε, * ) with initial state e α . As the extinction probability of a supercritical branching process is continuous (see Appendix A.3) with respect to the birth and death rates of this process, increases with the death rate and decreases with the birth rate, we find for α ∈ A that

β (ε,+) αα = b 1 + β 1 2 ρ α + ε 1/6 + 2ε ξ ζ -A 0 ε - β 2 4 + 1 2 ρ ᾱ -ε 1/6 ζ -A 0 ε ζ + A 0 ε + ε ξ β (ε,+) α ᾱ = b 2 1 - β 2 2 (ρ ᾱ + ε 1/6 ), δ (ε,+) α = b -cA 0 ε, β (ε,-) αα = b 1 + β 1 2 (ρ α -ε 1/6 )(ζ -A 0 ε) ζ + A 0 ε + ε ξ - β 2 4 + 1 2 (ρ ᾱ + ε 1/6 ) - β 2 ε ξ ζ -A 0 ε β (ε,-) α ᾱ = b 2 1 - β 2 2 (ρ ᾱ -ε 1/6 )(ζ -A 0 ε) ζ + A 0 ε + ε ξ δ (ε,-) α = b + c A 0 ε + ε ξ
0 ≤ q (ε,-) α -q (ε,+) α → ε→0 0, (3.34)
and

q (ε,+) α ≤ q α ≤ q (ε,-) α ,
where we recall that q α has been defined by (2.6) for the process N. In other words, for * ∈ {-, +},

q (ε, * ) α -q α = o ε (1). (3.35)
Since the coupling is only valid on Σ ε , we still need to prove that the probabilities of extinction and invasion of the actual process N are also given by q α and 1 -q α respectively. To this aim, let us introduce the following stopping times, for * ∈ {-, +},

∀x ∈ R + , T (ε, * ) x := inf{t > 0, N (ε, * ) (t) = Kx }. (3.36)
Recall that, on Σ ε , the coupling (3.33) is satisfied and thus

P T (ε,-) ε ξ < T (ε,-) 0 , Σ ε ≤ P T P ε ξ < T P 0 , Σ ε ≤ P T (ε,+) ε ξ < T (ε,+) 0 , Σ ε . (3.37)
Indeed, if a process reaches the size ε ξ K before extinction, it is also the case for a larger process. However, Σ ε is independent from N (ε,-) and N (ε,+) , thus with (3.32), lim sup

K→∞ P T (ε, * ) ε ξ < T (ε, * ) 0 , Σ ε = lim sup K→∞ P T (ε, * ) ε ξ < T (ε, * ) 0 P (Σ ε ) ≥ (1 -q (ε, * ) (α) )(1 -C(A 0 , ξ)ε 1/12 ).
Then letting K go to infinity in (3.37), we find

(1 -q (ε,-) (α) )(1 -C(A 0 , ξ)ε 1/12 ) ≤ lim sup K→+∞ P T P ε ξ < T P 0 , Σ ε ≤ (1 -q (ε,+) (α) )(1 -C(A 0 , ξ)ε 1/12 ).
Finally, adding (3.32) and (3.35) we get lim sup It remains to prove that in the case of invasion (which happens with probability 1 -q α ), the time before reaching size Kε ξ is of order log K/λ, where we recall that λ is the maximal eigenvalue of the matrix J defined in (2.8), and that in the case of invasion, λ is positive.

K→∞ P(T P ε ξ < T P 0 ∧ R A 0 ε ∧ U ε 1/6 ) -(1 -q α ) ≤ lim sup K→∞ P(T P ε ξ < T P 0 , Σ ε ) -(1 -q α ) + lim sup K→∞ |P(Σ ε c )| = o ε (1). ( 3
We denote by λ (ε, * ) the maximal eigenvalue of the mean matrix for the process N (ε, * ) . This eigenvalue is positive for ε small enough, and converges to λ when ε converges to 0. In other words there exists a nonnegative function η going 0 at 0 such that, for any ε small enough,

λ (ε, * ) λ -1 ≤ η(ε) 2 . (3.39)
Thus, let us fix ε small enough such that the previous inequality holds. Then from the coupling (3.33), which is true on Σ ε ,

P T (ε,-) ε ξ ≤ T (ε,-) 0 ∧ ln K λ (1 + η(ε)), Σ ε ≤ P T P ε ξ ≤ T P 0 ∧ ln K λ (1 + η(ε)), Σ ε . (3.40)
Once again, with independence between Σ ε and N (ε, * ) , using classical results on bitype branching processes (see [START_REF] Athreya | Branching processes[END_REF]), and (3.39), yields that for ε small enough (at least such that η(ε) < 1), lim sup

K→∞ P T (ε,-) ε ξ ≤ T (ε,-) 0 ∧ ln K λ (1 + η(ε)), Σ ε = lim sup K→∞ P T (ε,-) ε ξ ≤ ln K λ (1 + η(ε)) P (Σ ε ) ≥ lim sup K→∞ P T (ε,+) ε ξ ≤ ln K λ (ε,-) 1 - η(ε) 2 (1 + η(ε)) P (Σ ε ) ≥ lim sup K→∞ P T (ε,+) ε ξ ≤ ln K λ (ε,-) 1 + η(ε) -η 2 (ε) 2 P (Σ ε ) ≥ 1 -q (ε,-) (α) 1 -C(A 0 , ξ)ε 1/12 .
In addition with (3.40), we deduce that lim sup

K→∞ P T P ε ξ < T P 0 ∧ R A 0 ε ∧ U ε 1/6 , T P ε ξ ≤ ln K λ (1 + η(ε)) ≥ 1 -q α -3κ,
as soon as ε is small enough. We can prove in a similar way, using the upper bound Z (ε,+) of the coupling, that lim inf

K→∞ P T P ε ξ < T P 0 ∧ R A 0 ε ∧ U ε 1/6 , T P ε ξ ≤ ln K λ (1 -η(ε)) ≤ 1 -q α + 3κ.
Putting all pieces together, we conclude the proof of Proposition 3.1.

Mean-field phase.

Once the mutant population size has reached an order K, the meanfield approximation (2.10) becomes a good approximation for the population dynamics (cf Lemma 2.1). An important question however is the initial condition of the dynamical system used as an approximation. Indeed, depending on the initial state, the system (2.10) may converge to various equilibria (see Appendix A.1 for a study of these equilibria). The initial state to be considered for the dynamical system and the convergence to a stable equilibrium are the subjects of Sections 3.2.1 and 3.2.2, respectively.

3.2.1.

Mutant A/a proportions. We have seen that when (2.9) is satisfied, then the mutant population dynamics is close to that of the supercritical bitype branching process N defined in (2.3). For such a process we are able to control the long time proportion of the different types of individuals. More precisely, Kesten-Stigum theorem (see [START_REF] Georgii | Supercritical multitype branching processes: the ancestral types of typical individuals[END_REF] for instance) ensures the following property, if λ is positive:

( NA (t), Na (t)) NA (t) + Na (t) → t→∞ (π A , π a ) almost surely
on the event of survival of N, where π is the positive left eigenvalue of J associated to λ such that π A + π a = 1.

The next proposition states that with a probability close to one for large K, if the mutant population reaches the size εK, we may choose a time when the proportion of type A individuals in the P -population belongs to [π A -δ, π A + δ], with δ > 0 small. Proposition 3.2. Let C > 2 be such that

C max{q A , q a } C -1 1-1/C < 1.
Assume that ρ A > 1/2 and that (2.9) holds. Let δ > 0 such that π A -δ > 1/2. Then under the same assumptions as Proposition 3.1,

lim inf K→∞ P ∃t ∈ T P ε , T P √ ε , εK C ≤ N P (t) ≤ √ εK, π A -δ < N AP (t) N P (t) < π A + δ T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 ≥ 1 -o ε (1). Proof. If π A -δ < N AP (T P ε ) N P (T P ε ) < π A + δ
there is nothing to show. Thus we assume that

N AP (T P ε ) N P (T P ε ) ≤ π A -δ.
The symmetric case, when N AP (T P ε )/N P (T P ε ) ≥ π A + δ can be treated with similar arguments. Then, we introduce the event

Σε := {T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6
} on which all calculus will be done. Our first aim is to prove that the time interval [T P ε , T P √ ε ] is large when ε is small and that the mutant population size is not too small on this interval. Precisely, we introduce, for any ε > 0, the stopping time On the other hand, we obtain from Lemma A.2 that

T (ε) ε/C := inf{t ≥ T P ε , N P (t) ≤ εK/
lim K→∞ P T P √ ε ≤ T P ε + ln ln 1/ε| Σε ≤ √ ε (ln 1/ε) b(1+β 1 ) , (3.42)
since the process of the total size of P -individuals is always stochastically bounded from above by a Yule process with birth rate b(1 + β 1 ). Notice that Lemmas A.1 and A.2 can be applied here because, as we assumed (2.9), the mutant P invades with a positive probability and the approximating process N is supercritical. With this in mind, we are now interested in the dynamics of the fraction of A-individuals in the P -population. Our aim is to find a suitable lower bound to N AP (t)/N P (t) to prove that this fraction cannot stay below π A -δ on the interval [T P ε , T P √ ε ] with a probability close to 1. The fraction is a semi-martingale and can be decomposed as

N AP (t) N P (t) = N AP (T P ε ) N P (T P ε ) + M P (t) + V P (t),
for any t ≥ T P ε , with M P a martingale and V P a finite variation process. Let us start with the martingale part, M P . Its predictable quadratic variation can be obtained as in (3.13) with P replacing p and by integrating between T P ε and t instead of 0 and T P ε . It gives the bound

M P t ≤ C 0 (t -T P ε ) sup T P ε ≤s≤t 1 N P (s) -1 ,
where C 0 is a finite constant. Hence

M P (T P ε +ln ln 1/ε)∧T (ε) ε/C ≤ C 0 ln ln 1/ε εK/C -1 (3.43)
and lim sup K→∞ P sup

T P ε ≤t≤(T P ε +ln ln 1/ε) |M P (t)| ≥ ε Σε ≤ lim sup K→∞   P   sup T P ε ≤t≤(T P ε +ln ln 1/ε)∧T (ε) ε/C |M P (t)| ≥ ε Σε   + P T (ε) ε/C < T P ε + ln ln 1/ε Σε   ≤ lim sup K→∞ 1 ε 2 E M P T P ε +ln ln 1/ε∧T (ε) ε/C Σε + √ ε (ln 1/ε) b(1+β 1 ) = √ ε (ln 1/ε) b(1+β 1 ) ,
(3.44) using Doob's martingale inequality to obtain the third line, and (3.41), (3.42) and (3.43) for the last one. In particular, the martingale is larger than -ε with a probability close to one.

It remains to deal with the finite variation process V P . Itô's formula with jumps gives the following formulation of V P :

V P (t) = t T P ε P (s) N AP (s) N P (s) N P (s) N P (s) + 1 ds, with P (s) [X] := β (P ) AA (N(s))X+β (P ) aA (N(s))(1-X) (1-X)-β (P ) aa (N(s))(1-X)+β (P ) Aa (N(s))X X, and
β (P )
AA and β

(P )
Aa are defined by (3.21) and (3.22). Notice that, when ε is small, the polynomial function P (s) is close, on the interval [0, 1], to the polynomial function

P [X] := βAA X + βaA (1 -X) (1 -X) -βaa (1 -X) + βAa X X
where the functions βi , i ∈ G are defined in (2.4). Since P [0] > 0, P [1] ≤ 0 the equation ẋ = P [x] has a unique positive equilibrium in (0, 1]. Since (π A , 1 -π A ) is a left eigenvector of matrix (2.8), direct computation ensures that π A is a root of P and thus corresponds to this equilibrium. Moreover, since ρ A > 1/2 we obtain that P [1/2] > 0, and we deduce that π A > 1/2 (therefore, δ is well defined) and that there exists a positive θ such that for any x < π A -δ, P [x] > θ. Using the continuity of polynomial functions with respect to their coefficients, we deduce the following property conditioning on Σε and for ε small enough:

∀s ∈ T P ε , T P √ ε , ∀ x ∈ (0, π A -δ) , P (s) [x] ≥ θ 2 > 0. (3.45) Let us introduce τ (ε) A := inf t ≥ T P ε , N AP (t) N P (t) ≥ π A -δ .
From (3.44) and (3.45), we thus obtain that, conditioning on Σε , for any t ∈ [T P ε , (

T P ε + ln ln 1/ε) ∧ τ (ε) A ] π A -δ ≥ N AP (t) N P (t) ≥ θ 4 ln ln 1/ε ∧ (τ (ε) A -T P ε ) -ε, (3.46)
with a probability higher than 1 -√ ε (ln 1/ε) b(1+β 1 ) . Since θ 4 ln ln 1/ε -ε converges to +∞ with ε, τ (ε)

A is smaller than T P ε + ln ln 1/ε and so it is smaller than T P √ ε with a probability close to one (conditioning on Σε ), as soon as ε is sufficiently small, according to (3.42).

Finally, notice that each step of the process N AP (t)/N P (t) is smaller than (εK/C + 1) -1 , hence it is smaller than δ as soon as K is sufficiently large. Thus, after time τ (ε) A , the process will belong to the interval [π A -δ, π A + δ], at least for some times, if K is sufficiently large. This ends the proof of Proposition 3.2.

3.2.2.

Convergence of the dynamical system. In this section, we will study the behaviour of the dynamical system (2.10) after the 'stochastic' phase.

The following proposition states that the equilibrium without mutant is unstable under condition (2.9), and Proposition 3.4 states the convergence of the solution to (2.10) under suitable conditions. Proposition 3.3. Assume that (2.9) holds. Then for every ρ A ∈ [0, 1],

• the equilibrium

(0, ρ A (b -d)/c, 0, (1 -ρ A )(b -d)/c) is unstable
• the branching process N whose transition rates are given in (2.4) is supercritical On the opposite, if (2.9) does not hold, the largest eigenvalue of the Jacobian matrix for (2.10) is 0. In any case, the equilibrium (0,

ρ A (b -d)/c, 0, (1 -ρ A )(b -d)/c) is non-hyperbolic.
Proof. We compute the Jacobian matrix of system (2.10) at the equilibrium (0, ρ A (bd)/c, 0, (1 -ρ A )(b -d)/c), and obtain when reordering lines and columns (z

AP , z aP , z Ap , z ap )       J = b 2   ρ A β 1 -(1 -ρ A ) β 2 2 + 1 -ρ A β 2 2 -1 -(1 -ρ A ) β 2 2 -1 (1 -ρ A )β 1 -ρ A β 2 2 + 1   * 0 J = -(b -d) ρ A ρ a ρ A ρ a      
Therefore the eigenvalues of this matrix are the eigenvalues of the two sub-matrices J and J.

The eigenvalues of J are 0 and -(b -d) < 0.

Let us notice that the matrix J admits a positive eigenvalue if and only if either T r(J) > 0 or ∆(J) < 0 where

   T r(J) = b 2 β 1 -β 2 2 + 1 ∆(J) = b 2 4 ρ A (1 -ρ A )(β 1 + β 2 )(β 1 + 2) -β 1 2 (β 2 + 2)
We thus obtain that the equilibrium under consideration is unstable if one of the following conditions is satisfied:

β 1 > β 2 2 + 1 or ρ A (1 -ρ A ) < β 1 (β 2 + 2) 2(β 1 + β 2 )(β 1 + 2)
.

But from a functional study, we can prove that the function

β 1 → β 1 (β 2 +2)/(2(β 1 +β 2 )(β 1 + 2)) is larger than 1/4 for any β 1 ∈]β 2 , β 2 /2 + 1]
. This concludes the proof for the stability of the equilibrium point

(0, ρ A (b -d)/c, 0, (1 -ρ A )(b -d)/c).
Concerning the bitype branching process N, recall that J is also the mean matrix associated to it. As a consequence, N is supercritical if and only if the maximal eigenvalue of J is positive, and the conditions are the same.

Proposition 3.4. Let us consider an initial condition z 0 such that z Ap (0) > z ap (0) and z AP (0) > z aP (0). Let us furthermore assume that one of the following conditions is satisfied:

β 1 > β 2 or z A (0)z a (0) z(0) 2 < β 1 (β 2 + 2) 2(β 1 + β 2 )(β 1 + 2) . (3.47)
Then the solution z (z 0 ) of the system (2.10) converges as t → ∞ toward

(((1 + β 1 )b -d)/c, 0, 0, 0).
Proof. To prove the convergence we will consider the diversity at locus A/a using the quantity D := z A z a z 2 and prove that this quantity converges to 0. The differential equation followed by D is:

Ḋ = b z 2 (z AP z a + z A z aP ) D(β 1 + 2β 2 ) - β 2 2 -Dβ 1 (z AP z A + n aP z a ) = b z 2 (z AP z a + z A z aP ) β 2 2 (4D -1) -Dβ 1 (z AP -z aP ) 2 -Dβ 1 (z AP -z aP )(z Ap -z ap ) ≤ -D bβ 1 z 2 (z AP -z aP )(z Ap -z ap ), (3.48
) since D is always less than 1/4. Let us introduce the function

Π(t) := (z AP (t) -z aP (t))(z Ap (t) -z ap (t)).
Under the assumption of Proposition 3.4, Π(0) > 0. We want to prove that Π(t) > 0 for all t > 0. We start by computing the derivative of this quantity :

dΠ dt =2(b -d -cz)Π + b(z Ap -z ap ) z β 1 (z AP (z AP + z Ap 2 ) -z aP (z aP + z ap 2 )) + (z aP z Ap -z AP z ap ) + b(z AP -z aP ) z β 1 (z Ap z AP 2 -z ap z aP 2 ) -(z aP z Ap -z AP z ap ) . (3.49)
By reorganizing the terms, we find

dΠ dt = 2b -2d + bβ 1 3z AP + 2z aP + z ap 2z -b z aP + z ap z -2cz Π + bβ 1 2z [z AP (z Ap -z ap ) 2 + z ap (z AP -z aP ) 2 ] + b z [z aP (z Ap -z ap ) 2 + z ap (z AP -z aP ) 2 ] ≥ (b -2d -2cz) Π.
(3.50) We thus need information on the dynamics of z to conclude. From (2.10), we obtain

dz(t) dt = z(b -d -cz) + b β 1 z (z AP z A + z aP z a ) -b β 2 z (z AP z a + z aP z A ) ≤ z(t) (b(1 + β 1 ) -d -cz(t)) .
In other words, for any t ≥ 0,

z(t) ≤ z(0) ∨ b(1 + β 1 ) -d c .
Combining with (3.50) we deduce that as long as Π(t) > 0,

dΠ(t) dt ≥ b -2d -2c z(0) ∨ b(1 + β 1 ) -d c Π(t),
and thus

Π(t) ≥ Π(0)e -Ct > 0, ∀t ≥ 0.
Combining this result with (3.48), we deduce that D is a positive and decreasing quantity and converges to a nonnegative value where its derivative Ḋ vanishes. We deduce from the fact that all three terms of the second line of (3.48) are negative that D(z AP -z aP ) 2 = 0.

From Proposition A.1, the possible limits are the points (0, 0, 0, 0),

χ AP := (1 + β 1 )b -d c , 0, 0, 0 , b(1 + (β 1 -β 2 )/2) -d 2c , 0, b(1 + (β 1 -β 2 )/2) -d 2c , 0 ,
and the line

0, π b -d c , 0, (1 -π) b -d c with π ∈ [1/2, 1].
The proof of Proposition A.1 (i) ensures that no positive trajectory converges to the null point. Moreover, we proved that D is decreasing. As it starts from

D(0) = z A (0)(z(0) - z A (0))/z(0) 2 < 1/4
, the set of possible limits is thus restricted to the point χ AP or the line

0, π b -d c , 0, (1 -π) b -d c , π ∈ [z A (0)/z(0), 1]. (3.51)
As D is decreasing, the trajectory cannot oscillate close to the line of (3.51). Hence, if it approaches the line in large time, it should converge to a point of this line. Assume that it converges to (0, π(b -d)/c, 0, (1 -π)(b -d)/c). Note that, from Assumption (3.47),

β 1 > β 2 or π(1 -π) ≤ D(0) < β 1 (β 2 + 2) 2(β 1 + β 2 )(β 1 + 2) meaning that the equilibrium (0, π(b -d)/c, 0, (1 -π)(b -d)/c
) is unstable. Hence, from Perron-Frobenius Theorem and Proposition 3.3, there exists (γ 1 , γ 2 ) left positive principal eigenvector of the matrix J (see the proof of Proposition 3.3) which is positive and associated to a positive eigenvalue λ. Using similar computations, we obtain that in the neighbourhood of (0,

π(b -d)/c, 0, (1 -π)(b -d)/c) γ 1 żAP + γ 2 żaP ≥ λ 2 (γ 1 z AP + γ 2 z aP ).
Thus, as soon as z AP and z aP are not equal to 0, the quantity γ 1 z AP + γ 2 z aP will grow exponentially fast when the trajectory is close to (0, π(b -d)/c, 0, (1 -π)(b -d)/c), and therefore it cannot converge to this state.

3.3. Extinction. After the deterministic phase, the process is close to the state ((b(β 1 + 1)d)/c, 0, 0, 0). In this subsection, we are interested in estimating the time before the extinction of all but AP -individuals in the population. We also need to check that the AP -population size stays close to its equilibrium during this extinction time. We recall here the definition of the set S ε and the stopping time T Sε in (2.11) and (2.12), respectively:

S ε := b(1 + β 1 ) -d c -ε, b(1 + β 1 ) -d c + ε × {0} × {0} × {0}, T Sε := inf{t ≥ 0, Z K (t) ∈ S ε }.
Proposition 3.5. There exist two positive constants ε 0 and C 0 such that for any ε ≤ ε 0 , if there exists η ∈]0, 1/2[ that satisfies

Z AP (0) - b(1 + β 1 ) -d c ≤ ε and ηε/2 ≤ Z Ap (0) + Z ap (0) + Z aP (0) ≤ ε/2, then for all C > 2/(bβ 1 ) + C 0 ε, P (T Sε ≤ C log(K)) → K→+∞ 1, for all 0 ≤ C < 2/(bβ 1 ) -C 0 ε, P (T Sε ≤ C log(K)) → K→+∞ 0.
Proof. This proof is very similar to the proof of Proposition 4.1 in [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]. We thus only detail parts of the proof that are significantly different. Following these ideas, we prove that as long as the sum Z Ap + Z ap + Z aP is small (lower than ε), the process Z AP stays close to (b(1 + β 1 ) -d)/c.

Then, we can bound the death and birth rates of Z Ap , Z ap and Z aP under the previous approximation and compare the dynamics of these three processes with the ones of

N Ap (t) K , N ap (t) K , N aP (t) K , t ≥ 0 ,
where (N Ap , N ap , N aP ) ∈ N 3 is a three-types branching process with types Ap, ap and aP and such that The goal is thus to estimate the extinction time of such a sub-critical three type branching process. According to [START_REF] Athreya | Branching processes[END_REF] p. 202 and Theorem 3.1 in [START_REF] Heinzmann | Extinction times in multitype markov branching processes[END_REF],

•
P (N Ap (t), N ap (t), N aP (t)) = (0, 0, 0) = (1-c 1 e rt ) Z Ap (0)K (1-c 2 e rt ) Zap(0)K (1-c 3 e rt ) Z aP (0)K , (3.52)
where c 1 , c 2 and c 3 are three positive constants and r is the largest eigenvalue of

b 2   -β 1 1 -β 2 2 0 0 -2β 1 -1 -β 2 2 0 0 1 -β 2 2 -2β 1 -2β 2   ,
which is r = -bβ 1 /2. From (3.52), we deduce that the extinction time is of order (2/bβ 1 ) log K when K tends to +∞ by arguing as in step 2 in the proof of Proposition 4.1 in [START_REF] Coron | A stochastic model for speciation by mating preferences[END_REF]. This concludes the proof.

3.4.

Proof of Theorem 1. The proof strongly relies on the coupling (3.33). More precisely, we consider a trajectory of N (defined in (2.3)) coupled with the mutant process. The random variable B of Theorem 1 is then defined as

B := 1 { T0 =∞} ,
which equals 1 if the process N survives and 0 otherwise. In particular, B is indeed a Bernoulli random variable with parameter 1 -q α where α ∈ A is the genotype of the first mutant individual. Let the function η be defined as in Proposition 3.1. The convergence in probability claimed in (2.14) is equivalent to

lim inf K→∞ P T Sµ ∧ T P 0 ln K , 1 {T Sµ <T P 0 } - 1 λ + 2 bβ 1 , 1 B 1 ≤ η(ε) ≥ 1 + o ε (1). (3.53)
As ε is as small as we need, we can assume without loss of generality that η(ε) < 1. In the sequel, we divide the probability into two terms according to the values of B using that {B = 1} = { T0 = +∞}, we obtain

P T Sµ ∧ T P 0 ln K , 1 {T Sµ <T P 0 } - 1 λ + 2 bβ 1 , 1 B 1 ≤ η(ε) = P T Sµ ln K - 1 λ + 2 bβ 1 ≤ η(ε), T Sµ < T P 0 , T0 = +∞ + P T P 0 ln K ≤ η(ε), T P 0 < T Sµ , T0 < +∞ =: F(K, ε) + G(K, ε). (3.54)
Let us first consider G(K, ε), which is simpler to deal with and which represents the case of extinction of P -individuals. We introduce A 0 , C > 2, δ > 0 and z(0) as in Propositions 3.1, 3.2 and 3.4. First of all, notice that

G(K, ε) ≥ P T P 0 ln K ≤ η(ε), T P 0 < T Sµ , T0 < +∞, T P 0 < T P ε ∧ R A 0 ε ∧ U ε 1/6 .
Then if ε is chosen small enough such that ε < ((b(1 + β 1 ) -d)/c -µ) ∧ z p (0)/A 0 and considering our initial conditions, we have

T P ε ∧ R A 0 ε ∧ U ε 1/6 < T Sµ a.s., hence
T P 0 < T Sµ , T P 0 < T P ε ∧ R A 0 ε ∧ U ε 1/6 = T P 0 < T P ε ∧ R A 0 ε ∧ U ε 1/6 a.s.
and

G(K, ε) ≥ P T P 0 ln K ≤ η(ε), T0 < +∞, T P 0 < T P ε ∧ R A 0 ε ∧ U ε 1/6 . ( 3
P T0 < ∞ T P 0 < T P ε ∧ R A 0 ε ∧ U ε 1/6 = o ε (1) (3.56)
and

lim sup K→∞ P T0 < ∞ T (ε,+) 0 < T (ε,+) ∞ = o ε (1).
In addition with (3.55), we thus deduce

lim inf K→∞ G(K, ε) ≥ lim inf K→∞ P λT P 0 ln K ≤ η(ε), T (ε,+) 0 < T (ε,+) ∞ + o ε (1) ≥ lim inf K→∞ P λT (ε,+) 0 ln K ≤ η(ε), T (ε,+) 0 < T (ε,+) ∞ + o ε (1) ≥ lim inf K→∞ P T (ε,+) 0 < T (ε,+) ∞ + o ε (1) = q α + o ε (1), (3.57)
where the second inequality comes from coupling (3.33). This allows the case of extinction to be processed.

Let us now deal with F(K, ε), which represents the case of survival and invasion of Pindividuals. Firstly, reasoning as for (3.56) but with ξ = 1/2, we can get

lim sup K→∞ P T0 = ∞ T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 = o ε (1). Hence lim inf K→∞ F(K, ε) = lim inf K→∞ P T Sµ ln K - 1 λ + 2 bβ 1 ≤ η(ε), T Sµ < T P 0 , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 + o ε (1). (3.58)
We introduce two sets for any ε > 0, µ > 0,

K 1 ε = [π A -δ, π A + δ] × ε C , √ ε × ρ A -ε 1/6 , ρ A + ε 1/6 × b -d c -A 0 ε, b -d c + A 0 ε , K 2 µ = b(1 + β 1 ) -d c - µ 2 , b(1 + β 1 ) -d c + µ 2 × 0, µ 2 3
.

as well as the stopping times

T 1 ε = inf t ≥ 0, N AP (t) N P (t) , N P (t) K , N Ap (t) N p (t) , N p (t) K ∈ K 1 ε , T 2 µ = inf t ≥ T 1 ε , Z K (t) ∈ K 2 µ .
Our aim is essentially to prove that the only path to S µ is through K 1 ε and K 2 µ , as presented in the introduction of the paper. Then, using the Markov property and the previous propositions, we want to estimate T Sµ , by dividing [0, T Sµ ] into three parts:

[0, T 1 ε ], [T 1 ε , T 2 µ ] and [T 2 µ , T Sµ ]. From (3.58), lim inf K→∞ F(K, ε) ≥ lim inf K→∞ P T Sµ ln K - 1 λ + 2 bβ 1 ≤ η(ε), T Sµ < T P 0 , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 , T 2 µ < T Sµ + o ε (1) ≥ lim inf K→∞ P T 1 ε ln K - 1 λ ≤ η(ε) 3 , T 2 µ -T 1 ε ln K ≤ η(ε) 3 , T Sµ -T 2 µ ln K - 2 bβ 1 ≤ η(ε) 3 , T Sµ < T P 0 , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 , T 2 µ < T P 0 ∧ T Sµ + o ε (1)
Then, since for ε sufficiently small, R A 0 ε ∧ U ε 1/6 ≤ T Sµ a.s. and using the Markov property at times T 2 µ and T 1 ε we obtain lim inf

K→∞ F(K, ε) ≥ lim inf K→∞ P T 1 ε ln K - 1 λ ≤ η(ε) 3 , T 1 ε < T P 0 , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 × inf z(0)∈K 1 ε P T 2 µ ln K ≤ η(ε) 3 , T 2 µ < T P 0 Z(0) = z(0) × inf z(0)∈K 2 µ P T Sµ ln K - 2 bβ 1 ≤ η(ε) 3 , T Sµ < T P 0 Z(0) = z(0) + o ε (1).
(3.59)

To complete the proof it remains to show that r.h.s of (3.59) is close to 1 -q α when K goes to ∞ and ε is small. Let us start with the first term. Our aim is to prove that lim inf

K→∞ P T 1 ε ln K - 1 λ ≤ η(ε) 3 , T 1 ε < T P 0 , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 = 1 -q α + o ε (1). (3.60)
To this aim, let us notice that the following series of inequalities holds:

P λT 1 ε ln K -1 ≤ η(ε), T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 ≥ P λT 1 ε ln K - λT P √ ε ln K ≤ η(ε) 2 , λT P √ ε ln K -1 ≤ η(ε) 2 , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 ≥ P λT P √ ε ln K - λT P ε ln K ≤ η(ε) 2 , T P ε ≤ T 1 ε ≤ T P √ ε , λT P √ ε ln K -1 ≤ η(ε) 2 , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6
Now, if A, B, C and D are events, we have

P(A ∩ B ∩ C ∩ D) = P(C ∩ D) -P (A ∩ B) C ∩ C ∩ D) ≥ P(C ∩ D) -P A C ∩ D -P B C ∩ D .
Applying this to the previous series of inequalities yields

P λT 1 ε ln K -1 ≤ η(ε), T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 ≥ P λT √ ε ln K -1 ≤ η(ε) 2 , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 -P λT P √ ε ln K - λT P ε ln K ≥ η(ε) 2 , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 -P T 1 ε / ∈ T P ε , T P √ ε , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 (3.61)
Proposition 3.1 implies that the first term in the right hand side of (3.61) satisfies lim inf

K→∞ P T P √ ε ln K - 1 λ ≤ η(ε) 2 , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 ≥ 1 -q α -o ε (1).
From Lemma A.2, we deduce that the second term of the right hand side of (3.61) satisfies lim inf

K→∞ P λT P √ ε ln K - λT P ε ln K ≥ η(ε) 2 , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 = o ε (1).
Finally Proposition 3.2 implies that the last term of the right hand side of (3.61) satisfies:

lim inf K→∞ P T 1 ε / ∈ T P ε , T P √ ε , T P √ ε < T P 0 ∧ R A 0 ε ∧ U ε 1/6 = o ε (1).
This leads to (3.60).

Then we deal with the second term of (3.59) by using Proposition 3.4. Using the continuity of flows of the dynamical system (2.10) with respect to the initial condition and the convergence given by Proposition 3.4, we get that there exist ε 0 , δ 0 > 0 such that for all ε ≤ ε 0 , δ ≤ δ 0 , there exists a t µ,δ,ε > 0 such that for all t ≥ t µ,δ,ε

z (z 0 ) (t) - b(1 + β 1 ) -d c , 0, 0, 0 ∞ ≤ µ 4 ,
for every initial condition z 0 = (z 0 AP , z 0 aP , z 0 Ap , z 0 ap ) such that (z AP /z P , z P , z Ap /z p , z p ) belongs to K 1 ε . Now using Lemma 2.1, we get that for any µ > 0, and ε < ε 0 ,

lim K→∞ P T 2 µ -T 1 ε ≤ t ν,δ,ε |Z(0) ∈ K 1 ε = 1.
In other words, the second term of (3.59) is close to 1 when K converges to ∞ and ε is small.

Finally, we deal with the third term of (3.59). Applying Proposition 3.5 we obtain that there exists µ 0 (defined by ε 0 in Proposition 3.5) such that for µ < µ 0 and ε small enough,

lim K→∞ P T Sµ ln K - 2 bβ 1 ≤ η(ε) 3 Z(0) ∈ K 2 µ = 1. (3.62)
By combining (3.60), the convergence of the second term of (3.59) to 1 and (3.62), we get

lim inf K→∞ F(K, ε) ≥ 1 -q α + o ε (1).
In addition with (3.57) and (3.54), we deduce (3.53). Finally (2.15) derives from (3.5) which ends the proof of Theorem 1.

Thus, we obtain that ∀t ≤ t ε ,

żA -ża ≥ (z A -z a )(b -d -cε) -bβ 1 ε.
Let φ be the unique solution to the linear differential equation

φ = φ(b -d -cε) -bβ 1 ε. Then φ(t) = (φ(0) -ε(b -d -cε) -1 )e (b-d-cε)t + bβ 1 ε(b -d -cε) -1 .
Using classical results on differential inequalities we deduce that if

(z A -z a )(0) = φ(0) then for all t ≤ t ε , z(t) ≥ (z A -z a )(t) ≥ φ(t).
Since for ε small enough φ(t) → ∞ as t → ∞, we deduce that t ε is finite. In other words, (0, 0, 0, 0) is unstable.

żAP = 0 = z AP (b -d -cz) + bβ 1 z z AP - bβ 2 z z aP , (A.2)
and

żaP = 0 = z aP (b -d -cz) + bβ 1 z z aP - bβ 2 z z AP . (A.3)
If z AP = 0, then z aP = ((1 + β 1 )b -d)/c, and similarly when exchanging A and a. For these equilibria, the eigenvalues of the Jacobian matrix are:

- bβ 1 2 , -b(β 1 + β 2 ), - b 4 (2 + 4β 1 + β 2 ), -b(1 + β 1 ) + d .
Since b > d, these eigenvalues are negative and these equilibria are therefore stable. If z AP > 0 and z aP > 0 then by dividing (A.2) by z AP and (A.3) by z aP and making the difference between both expressions, we get:

b z (β 1 + β 2 ) (z AP -z aP ) = 0. Then z AP = z aP = b(1 + (β 1 -β 2 )/2) -d 2c from (A.2).
The eigenvalues of the Jacobian matrix in this equilibrium are:

b 2 (β 1 + β 2 ), b 4 (β 2 -β 1 ), - b 4 (2 + β 1 -2β 2 ), - b 2 (2 + β 1 -β 2 ) + d .
The first eigenvalue is positive, therefore this equilibrium is unstable. We finally show that there is no other equilibrium with at least a null coordinate. To this aim, we first consider the case where z a = 0. Then

żA = (b -d -cz)z A + bβ 1 z z AP z A = 0,
and z A = 0 (which is the trivial equilibrium (i)) or b -d -cz + bβ 1 z AP /z = 0. In the second case, from the equation satisfied by z AP , we deduce that bβ 1 z AP z Ap /(2z) = 0. Hence either z AP = 0 or z aP = 0 (which corresponds to Equilibrium (ii)). Similar equilibria are retrieved by assuming z A = 0.

Finally consider the case where z ap = 0. Then from the equation satisfied by z ap given in (A.1),

żap = 0 = b 2z z aP z Ap (1 - β 2 2
).

Therefore z aP = 0 (then z a = 0 which corresponds to the case that has just been considered) or z Ap = 0 (which corresponds to Equilibrium (iii)). Similar arguments can be made assuming z Ap = 0 or z aP = 0 or z AP = 0.

A.2. Proof of Proposition 2.1. In the particular case where ρ A = 1, the transition rates of the bitype branching process N are equal to

βAA = b 2 (2 + β 1 ), βaa = βaA = b 2 1 - β 2 2
and βAa = 0, and the system (2.7) giving the extinction probabilities of the branching process N takes the simpler form:

u A (s A , s a ) = b(1 -s A ) + b 2 (2 + β 1 )(s 2 A -s A ) u a (s A , s a ) = b(1 -s a ) + b 2 1 - β 2 2 (s 2 a -s a ) + b 2 1 - β 2 2 (s A s a -s a ).
Recall that the extinction probabilities we are looking for are the smallest solution to u α (s A , s a ) = 0, α ∈ A. We easily obtain from the first, linear, equation an expression of q A . Then replacing it with its expression in the second equation gives that q a is the root of a second order polynomial function. This gives the result.

A.3. Probabilistic technical results.

Lemma A.1. Let ( NA , Na ) be a two type supercritical birth and death process. We recall that, for α ∈ A, q α is the extinction probability of the process when the initial individual is of type α q α = P ∃t > 0, NA (t) + Na (t) = 0 ( NA (0), Na (0)) = e α < 1,

where we denote by (e A , e a ) the canonical basis of R 2 . Let C > 2 satisfying

C max{q A , q a } C -1 1-1/C < 1. (A.4) Then lim k→∞ P S k/C < ∞| NA (0) + Na (0) = k = 0,
where the stopping time S l is defined for any l ∈ N by S l := inf{t ≥ 0, NA (t) + Na (t) = l}.

Proof. Let us first remark that it is possible to choose such a constant C > 2 since the map x → x(max(q A , q a )/(x -1)) 1-1/x is continuous and goes to max(q A , q a ) < 1 as x → ∞.

There are initially k individuals and we want to lower bound the probability that the population size reaches k/C . If this happens at a finite time, then it means that, at least, k -k/C individuals alive at time 0 have a finite line of descent. But we know that each individual has a finite line of descent with a probability smaller than q := max(q A , q a ). Then using the branching property, the probability that exactly i initial individual out of k have a finite line of descent is smaller than k i q i (1 -q) k-i . Hence

P S k/C < ∞| NA (0) + Na (0) = k ≤ k i=k-k/C k i q i (1 -q) k-i .
Since i → q i and i → k i are decreasing functions as soon as i ≥ k/2 , we deduce that for C > 2 and k large

P S k/C < ∞| NA (0) + Na (0) = k ≤ k k -k/C q k-k/C k i=k-k/C (1 -q) k-i ≤ 1 q k! ( k/C )!(k -k/C )! q k-k/C . Moreover, using Stirling's formula we get k! ( k/C )!(k -k/C )! q k-k/C ∼ k→∞ √ kk k q k-k/C 2π k/C (k -k/C ) k/C k/C (k -k/C ) k-k/C ∼ k→∞ C 2πk(1 -1/C) kq k -k/C k k -k/C k/C (1 -λ A ∧ λ a ) k/C → k→∞ 0,
under assumption (A.4). This ends the proof. Proof. Using the Markov property of the process, we find P τ √ εK < τ εK + ln ln 1/ε = P τ √ εK < ln ln 1/ε X(0) = εK = P X(τ √ εK )e -bτ √ εK > √ εK (ln 1/ε) -b X(0) = εK . Now using Markov Inequality and the fact that conditioning on {X(0) = εK }, X(t)e -bt is a martingale with expectation εK , we obtain

P τ √ εK < τ εK + ln ln 1/ε ≤ εK √ εK (ln 1/ε) -b ∼ √ ε(ln 1/ε) b , for K → ∞.
This concludes the proof.

Lemma A.3. Let us consider a family of two-type branching processes ( N ε A , N ε a , ε ∈ R) whose transition rates are given by

( N ε A , N ε a ) → ( N ε A + 1, N ε a ) at rate b ε AA N ε A + b ε aA N ε a , ( N ε A , N ε a ) → ( N ε A , N ε a + 1) at rate b ε Aa N ε A + b ε aa N ε a , ( N ε A , N ε a ) → ( N ε A -1, N ε a ) at rate d ε N ε A , ( N ε A , N ε a ) → ( N ε A , N ε a -1)
at rate d ε N ε a , and let us denote by q ε = (q ε A , q ε a ) the extinction probabilities of the process N ε with initial state an individual of type A or a.

(i) Let us assume that the functions ε → b ε i > 0 for i ∈ A 2 (resp ε → d ε > 0) are of class C k for k ≥ 0 in ε = 0 and that the process ( N 0 A , N 0 a ) is supercritical. Then the application ε → q ε is of class C k in ε = 0. (ii) Let us assume furthermore that the functions ε → b ε i > 0 for i ∈ A 2 (resp ε → d ε > 0) are non decreasing (resp. non increasing), and consider ε 1 ≤ ε 2 then the extinction probabilities q i = (q ε i A , q ε i a ) (i ∈ {1, 2}) of the two branching processes ( N ε i A , N ε i a ) satisfy q 1 ≤ q 2 , where the inequality applies to both coordinates.

Proof. (i) The proof relies on Theorem 6.2 of [START_REF] Alili | On the regularity of the extinction probability of a branching process in varying and random environments[END_REF] that considers multi-type discrete time branching processes. The process ( N ε A , N ε a , ε ∈ R) is a continuous time multi-type linear birth-and-death process in which for all α 1 , α 2 ∈ A 0 , individuals with genotype α 1 die at rate d ε and produce an offspring with genotype α 2 at rate b ε α 1 ,α 2 . For all α 1 , α 2 ∈ A 0 , the random variable N α 1 ,α 2 giving the number of offsprings of type α 2 of a given individual of type α 1 satisfies

P(N α 1 ,α 2 = k) = b ε α 1 α 2 b ε α 1 A + b ε α 1 a + d ε k d ε b ε α 1 A + b ε α 1 a
+ d ε which is assumed to be C k in at 0. Let us consider the discrete time stochastic process with values in N 2 giving the number of individuals of each type at each generation, whose extinction probability is equal to q ε . Then for any (s A , s a ) ∈ [0, 1) 2 , To apply Theorem 6.2 of [START_REF] Alili | On the regularity of the extinction probability of a branching process in varying and random environments[END_REF], we therefore need to find s = (s A , s a ) ∈ [0, 1) 2 such that E(s

E(s

N α,A A s
Nα,a a

) < s α for all α ∈ A.

This condition is sufficient to check Assumption 6.1 of [START_REF] Alili | On the regularity of the extinction probability of a branching process in varying and random environments[END_REF] (in which α is now denoted (s 1 , s 2 )) because we use here the particular framework of constant environment. Therefore, by taking We have φ(1, 1) = (0, 0) and we seek (s 1 , s 2 ) ∈ [0, 1) 2 such that φ(s 1 , s 2 ) > (0, 0) where the inequality applies to both coordinates. This is possible if the jabobian matrix of the application φ in and N 1 a ≥ N 2 a . Therefore, for every t ≥ 0, P( N 1 (t) = 0) ≤ P( N 2 (t) = 0), which gives the result, by letting t → ∞.

Appendix B. Table for birth rates

In this table we present the birth rates and the possible offspring of every couples in the population. When a P individual is involved, we differentiate whether it is the choosing parent (1st parent) or the chosen one (2nd parent). Let us briefly recall how the table is constructed. For the possible offspring, we assume Mendelian reproduction meaning that for each gene independently an allele is chosen at random among the two alleles of the parent. As an example, in a mating AP × Ap the offspring will necessary receive allele A and then choose with equal probability between p and P , and we note in the third column 1/2AP and 1/2Ap. For the birth rate of the same couple, since the choosing parent carries allele P , mating occurs with a preference at rate b(1 + β 1 ) since both parent carry allele A. Table 1. This table gives the rates at which two parents with given genotypes give birth to an offspring with given genotype, for all possible values of these genotypes. By convention, the first parent is assumed to be responsible for homogamy, when carrying allele P .

  ) for η > 0 and i ∈ {∅, 1, 2} (by convention Z ∅ p = Z p ). Let us first consider the process Z 1 p . When K is large and according to [10, Chapter 11, Theorem 2.1 p. 456], the dynamics of Z 1 p is close to the dynamics of the unique solution to d dt z = z(b -d -cε ξ -cz) -bε ξ . (3.28)

  and ζ = (b -d)/c. Note that for (α, α ) ∈ A 2 and * ∈ {-, +}, the applications ε → β (ε, * ) αα and ε → δ (ε, * ) α are continuous and converge respectively as ε → 0 to βα,α and b which are the birth and death rates of the process N introduced in (2.3). Moreover β

  .38) Equation (3.5) is derived similarly.

  C}, where C is the constant introduced in Proposition 3.2, and we want to prove that the stopping time T P √ ε is larger than T P ε + ln ln(1/ε) and smaller than T (ε) ε/C . On the one hand, we obtain from coupling (3.33), satisfied on Σε , and Lemma A.1, that lim

  any Ap-individual gives birth to a Ap-individual at rate b(2 + β 1 )/2, • any aP -individual gives birth to a aP -individual at rate b(1 -β 2 ), • any individual gives birth to a ap-individual at rate b(2 -β 2 )/4, • any individual dies at rate b(1 + β 1 ).

  (ii): Let us assume that z AP = z aP = 0. Then, (A.1) can be reduced to żAp = (b -d -cz)z Ap żap = (b -d -cz)z ap . Therefore, the set of points (0, z Ap , 0, (b -d)/c -z Ap ) with z Ap ∈ [0, (b -d)/c] corresponds to the set of non null equilibria such that z AP = z aP = 0. (iii): Let us assume that z Ap = z ap = 0. Then from (A.1),

Lemma A. 2 .

 2 Let us consider a one dimensional pure birth process X with birth rate b. Denote for k > 0 by τ k the hitting time of k by the process X. Then there exists a finite C such that lim sup K→∞ P τ √ εK < τ εK + ln ln 1/ε ≤ C √ ε(ln 1/ε) b .

  + b αa + d -b αA s 1 -b αa s 2 .

)

  < sα < s α for all α ∈ A 2 which is exactly Assumption 6.1 of[START_REF] Alili | On the regularity of the extinction probability of a branching process in varying and random environments[END_REF]. For any(s 1 , s 2 ) ∈ [0, 1) 2 , let φ(s 1 , s 2 ) = (s 1 (b 11 + b 12 + d -s 1 b 11 -s 2 b 12 ) -d, s 2 (b 21 + b 22 + d -s 1 b 21 -s 2 b 22 ) -d).

( 1 , 1 ) 1 A ≥ N 2 A

 1112 which is equal to d -b 11 -b 12 -b 21 d -b 22 has a negative eigenvalue and this condition is equivalent to the supercriticality of the process ( N 0 A , N 0 a ). (ii) The proof relies on a coupling argument. Let us construct the two processes using the same Poisson point measures. Then we have that almost surely N

1 + β 1 ) n AP n AP n aP aP aP b( 1 + β 1 ) n aP n aP n aP AP 1 2 aP b( 1 -β 2 ) n aP n AP n 1 2 1 -β 2 ) n aP n AP n 1 2

 11111122122 AP b(1 -β 2 ) n aP n AP n AP b(1 -β 2 ) n aP n AP n

  Values of q A as a function of ρ A for different values of β 1 and β 2 .

		1.00								1.00
		0.95
										0.98
		0.90
	q_A	0.85		beta_1=0.2000					q_A	0.96	beta_2=0.0000
		0.80		beta_1=0.3333 beta_1=0.4667						0.94	beta_2=0.1667 beta_2=0.3333
		0.75		beta_1=0.6000							beta_2=0.5000
		0.70		beta_1=0.7333 beta_1=0.8667						0.92	beta_2=0.6667 beta_2=0.8333
				beta_1=1.0000							beta_2=1.0000
		0.65	0.0	0.2	0.4	rho_A	0.6	0.8	1.0	0.0	0.2	0.4	rho_A	0.6	0.8	1.0
			Figure 1.

Acknowledgments

The authors thank the CNRS for its financial support through its competitive funding programs on interdisciplinary research. This work was partially funded by the Chair "Modélisation Mathématique et Biodiversité" of VEOLIA-Ecole Polytechnique-MNHN-F.X. H.L. acknowledges support from CONACyT-MEXICO and the foundation Sofía Kovalevskaia of SMM.

Appendix A. Technical results

A.1. Equilibria of dynamical system (2.10). In this section, we study the existence and the stability of some equilibria of the dynamical system (2.10). For the sake of readability, we explicitly rewrite the dynamical system below

where z = z AP + z Ap + z aP + z ap is the total size of the population.

Proposition A.1. The dynamical system (A.1) admits the following equilibria, with at least a null coordinate: (i): The state (0, 0, 0, 0), which is unstable.

(ii): Any state where only allele p remains at locus 2,

The stability of these equilibria has been studied in Proposition 3.3. (iii): The three following states for which only allele P remains at locus 2

The first two equilibria are stable, whereas the last one is unstable.

Proof.

(i): The state (0, 0, 0, 0) is an equilibrium, from Equation (A.1). To prove that it is unstable, let us consider ε > 0 and assume that the initial condition z 0 satisfies z(0) = ||z 0 || 1 ≤ ε. We denote by t ε = inf{t ≥ 0, z(t) > ε} which would be infinite if z lies in the basin of attraction of (0, 0, 0, 0). From (A.1), we find żA -ża = (b -d -cz) (z A -z a ) + bβ 1 z (z AP z A -z aP z a ).