The recovery of a recessive allele in a Mendelian diploid model

TL;DR

This paper demonstrates how a recessive allele can re-emerge in a diploid population through a formal model, highlighting the genetic robustness enabled by sexual reproduction and mutation dynamics.

Contribution

It provides a formal description of the mechanism allowing recessive alleles to reappear in populations, extending previous work on allele survival in sexual diploid models.

Findings

Recessive alleles can re-invade populations after mutations.

Sexual reproduction enhances genetic robustness and allele persistence.

The model formalizes the re-emergence process of recessive alleles.

Abstract

We study the large population limit of a stochastic individual-based model which describes the time evolution of a diploid hermaphroditic population reproducing according to Mendelian rules. In [Neukirch, Bovier, 2016] it is proved that sexual reproduction allows unfit alleles to survive in individuals with mixed genotype much longer than they would in populations reproducing asexually. In the present paper we prove that this indeed opens the possibility that individuals with a pure genotype can reinvade in the population after the appearance of further mutations. We thus expose a formal description of a mechanism by which a recessive allele can re-emerge in a population. This can be seen as a statement of genetic robustness exhibited by diploid populations performing sexual reproduction.

Tables

$f_{u_1{u}_2}\in {ℝ}_{+}$	be the per capita birth rate (fertility) of an individual with genotype
	$u_1{u}_2$,
$D_{u_1{u}_2}\in {ℝ}_{+}$	be the per capita natural death rate of an individual with genotype $u_1{u}_2$,
$K\in \mathbb{N}$	be the carrying capacity, a parameter which scales the population size,
$\frac{c_{u_1{u}_2,v_1{v}_2}}{K}\in {ℝ}_{+}$	be the competition effect felt by an individual with genotype $u_1{u}_2$ from an individual of genotype $v_1{v}_2$,
$R_{u_1{u}_2}\left(v_1{v}_2\right)\in \left\{0,1\right\}$	be the reproductive compatibility of the genotype $v_1{v}_2$ with $u_1{u}_2$,
${\mu }_K\in {ℝ}_{+}$	be the mutation probability per birth event. Here it is independent of
	the genotype,
$m(u,dh)$	be the mutation law of a mutant allelic trait $u+h\in 𝒰$, born from an individual with allelic trait $u$.

$T_1$	$:=T_{{\epsilon }_0}^{aa+aA+aB+AB+BB}$,
$T_{=}$	$:=T^{aA=aB}$,
$T_2$	$:=T^{aA=\delta AA}\wedge T^{aB=\delta AB}\wedge T^{aa=aA\wedge aB}$,
$T_3$	$:=T_{{\overline{n}}_a-{\epsilon }^{\gamma /2}}^{aa}$.

Equations

$${\nu }_t=\frac{1}{K}\sum _{i=1}^{N_t}{\delta }_{u_1^i(t)u_2^i(t)}.$$

$$\displaystyle\mathcal{M}^{K}=\left\{\frac{1}{K}\sum_{i=1}^{n}\delta_{u_{1}^{i}u_{2}^{i}}\;\Big{|}\;n\geq 0,u_{1}^{1}u_{2}^{1},...,u_{1}^{n}u_{2}^{n}\in\mathcal{U}^{2}\right\},$$

$$0\le f_{u_1{u}_2}\le \bar{f},0\le D_{u_1{u}_2}\le \bar{D}\text{and}0\le c_{u_1{u}_2,v_1{v}_2}\le \bar{c}.\text{vspace}0pt$$

$$D_{aa}=D,D_{aA}=D-\Delta ,$$

$$D_{aa}=D,D_{aA}=D-\Delta ,D_{AA}=D-2\Delta ,$$

$$D_{aa}=D,D_{aA}=D-\Delta ,D_{AA}=D-\Delta ,$$

$$\ln K\ll \frac{1}{{\mu }_{K}K}\ll K^{1/2-\alpha }\text{as }K\to \infty ,\text{for some }\alpha >0,$$

$${\bar{n}}_u=\frac{f_{uu}-D_{uu}}{c_{uu,uu}}.$$

$$S_{uv,uu}=f_{uv}-D_{uv}-c_{uv,uu}{\bar{n}}_u,$$

$$f_i\equiv f.$$

$$D_{aa}$$

$$D_{AA}\equiv D_{aA}$$

$$D_{aB}\equiv D_{AB}\equiv D_{BB}$$

$$\left(c_{i,j}\right)_{\{i,j\}\in\mathcal{G}\times\mathcal{G}}=:\begin{tabular}[]{c|c|c|c}&$aa$&$aA\quad AA$&$aB\quad AB\quad BB$\\ \hline\cr$aa$&$c$&$c\quad c$&$0\quad 0\quad 0$\\ \hline\cr$aA$&$c$&$c\quad c$&$c\quad c\quad c-\eta$\\ $AA$&$c$&$c\quad c$&$c\quad c\quad c$\\ \hline\cr$aB$&$0$&$c\quad c$&$c\quad c\quad c$\\ $AB$&$0$&$c\quad c$&$c\quad c\quad c$\\ $BB$&$0$&$c-\eta\quad c$&$c\quad c\quad c$\\ \hline\cr\end{tabular}$$

$$\left(R_{i}(j)\right)_{\{i,j\}\in\mathcal{G}\times\mathcal{G}}\equiv\begin{tabular}[]{c|c|c|c}&$aa$&$aA\quad AA$&$aB\quad AB\quad BB$\\ \hline\cr$aa$&$1$&$1\quad 1$&$0\quad 0\quad 0$\\ \hline\cr$aA$&$1$&$1\quad 1$&$1\quad 1\quad 1$\\ $AA$&$1$&$1\quad 1$&$1\quad 1\quad 1$\\ \hline\cr$aB$&$0$&$1\quad 1$&$1\quad 1\quad 1$\\ $AB$&$0$&$1\quad 1$&$1\quad 1\quad 1$\\ $BB$&$0$&$1\quad 1$&$1\quad 1\quad 1$\\ \hline\cr\end{tabular}$$

$$S_{AB,AA}$$

$$S_{aa,BB}$$

$$c{\bar{n}}_B$$

$${\Sigma }_3:=n_{aa}+n_{aA}+n_{AA},$$

$${\Sigma }_6:=n_{aa}+n_{aA}+n_{AA}+n_{aB}+n_{AB}+n_{BB},$$

$${\Sigma }_5:=n_{aA}+n_{AA}+n_{aB}+n_{AB}+n_{BB}.$$

$$b_{aa}=$$

$$b_{aA}=$$

$$+f\frac{\left(\frac{1}{2}{n}_{aA}+n_{AA}\right)\left(n_{aa}+n_{aA}+\frac{1}{2}{n}_{aB}\right)+\frac{1}{4}{n}_{aA}{n}_{AB}}{{\Sigma }_6},$$

$$b_{AA}=$$

$$b_{aB}=$$

$$b_{AB}=$$

$$b_{BB}=$$

$$d_{aa}$$

$$d_{aA}$$

$$d_{AA}$$

$$d_{aB}$$

$$d_{AB}$$

$$d_{BB}$$

$$\dot{n}(t)$$

$$i.e.{\dot{n}}_i(t)$$

$$\underset{K\to \infty }{\lim }\underset{t\in [0,T]}{\sup }|n_i^K(t)-{\tilde{n}}_i(t,x^0))|=0,\text{a.s.},$$

$${\bar{n}}_A\ge n_{AA}(0)$$

Full text

We acknowledge financial support from the German Research Foundation (DFG) through the Hausdorff Center for Mathematics, the Cluster of Excellence ImmunoSensation, and the Priority Programme SPP1590 Probabilistic Structures in Evolution. L.C. has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) through the Exploratory Project CanDyPop, as well as by the Swiss National Science Foundation through the grant No. P300P2_161031, and by the Chair ”Modélisation Mathématique et Biodiversité” of Veolia Environnement - École Polytechnique - Muséum National d’Histoire Naturelle - Fondation X.

We would like to thank Pierre Collet and Vincent Beffara for their help on the theory of dynamical systems and fruitful discussions.

1. Introduction

In population genetics, the study of Mendelian diploid models of fixed population size began more than a century ago (see e.g. [30, 16, 29, 18, 19, 12, 26, 15, 3]), while their counterparts of variable population size models were studied in the context of adaptive dynamics from 1999 onwards [22]. The approach of adaptive dynamics is to introduce competition kernels to regulate the population size instead of maintaining it constant, see [21, 23, 24].

Stochastic individual-based versions of these models appeared in the 1990s, see [13, 5, 4, 17, 6, 7]. They assume single events of reproduction, mutation, natural death, and death by competition happen at random times to each individual in the population. An important and interesting feature of these models is that different limiting processes on different time-scales appear as the carrying capacity tends to infinity while mutation rates and mutation step-size tend to zero (see [4, 13, 25, 7, 1]). One of the major results in this context is the convergence of a properly rescaled process to the so called Trait Substitution Sequence (TSS) process, which describes the evolution of a monomorphic population as a jump process between monomorphic equilibria. More generally, Champagnat and Méléard [7] obtained the convergence to a Polymorphic Evolution Sequence (PES), where jumps occur between equilibria that may include populations that have multiple co-existing phenotypes. The appearance of co-existing phenotypes is, however, exceptional and happens only at so-called evolutionary singularities. From a biological point of view, this is somewhat unsatisfactory, as it apparently fails to explain the biodiversity seen in real biological systems.

Most of the models considered in this context assume haploid populations with a-sexual reproduction. Exceptions are the paper [8] by Collet, Méléard and Metz from 2013 and a series of papers by Coron and co-authors [11, 9, 10] following it. In [8], the Trait Substitution Sequence is derived in a Mendelian diploid model under the assumption that the fitter mutant allele and the resident allele are co-dominant.

The main reason why both in haploid models and in the model considered in [8] the evolution along monomorphic populations is typical is that the time scales for the fixation of a new trait and the extinction of the resident trait are the same (both of order $\ln K$) (unless some very special fine-tuning of parameters occurs that allows for co-existence). This precludes (at least in the rare mutation scenarios considered) that an initially less fit trait survives long enough for several new mutations to occur, creating a situation where this trait may become fit again and recover.

In a follow-up paper to [8], two of the present authors [27], it was shown that, if instead one assumes that the resident allele is recessive, the time to extinction of this allele is dramatically increased. This will be discussed in detail in Section 1.2 and paves the way for the appearance of a richer limiting process.

The general framework in [8] and [27] is the following. Each individual is characterised by a reproduction and death rate which depend on a phenotypic trait determined by its genotype, which here is determined by two alleles (e.g. $A$ and $a$) on one single locus. The evolution of the trait distribution of the three genotypes $aa,aA$ and $AA$ is studied under the action of (1) heredity, which transmits traits to new offsprings according to Mendelian rules, (2) mutation, which produces variations in the trait values in the population onto which selection is acting, and (3) competition for resources between individuals.

The paper [27] proves that sexual reproduction allows unfit alleles to survive in individuals with mixed genotype much longer than they would in populations reproducing asexually. This opens the possibility that while this allele is still alive in the population, the appearance of new mutants alters the fitness landscape in such a way that is favourable for this allele and allows it to reinvade in the population, leading to a new equilibrium with co-existing phenotypes. The goal of this paper is to rigorously prove that such a scenario indeed occurs under fairly natural assumptions.

Recently, Billiard and Smadi [2] considered related questions for haploid individuals (performing clonal reproduction). The authors show that a deleterious allele can reinvade after a new mutation, but the range of parameters allowing this behaviour is though very small.

1.1. The stochastic model

The individual-based microscopic Mendelian diploid model is a non-linear birth-and-death process. We consider a model for a population of a finite number of hermaphroditic individuals which reproduce sexually. Each individual $i$ is characterised by two alleles, $u_{1}^{i}u_{2}^{i}$, taken from some allele space $\mathcal{U}\subset{\mathbb{R}}$. These two alleles define the genotype of the individual $i$. We suppress parental effects, which means that we identify individuals with genotype $u_{1}u_{2}$ and $u_{2}u_{1}$. Each individual has a Mendelian reproduction rate with possible mutations and a natural death rate. Moreover, there is an additional death rate due to ecological competition with the other individuals in the population. Let

Scaling the competition function $c$ down by a factor $1/K$ amounts to scaling the population size to order $K$. We are interested in asymptotic results when $K$ is large. We assume rare mutations, i.e. $\mu_{K}\ll 1$. If a mutation occurs at a birth event, only one allele changes from $u$ to $u+h$ where $h$ is a random variable with law $m(u,dh)$.

At any time $t$, there is a finite number, $N_{t}$, of individuals, each with genotype in $\mathcal{U}^{2}$. We denote by $u_{1}^{1}(t)u_{2}^{1}(t),...,u_{1}^{N_{t}}(t)u_{2}^{N_{t}}(t)$ the genotypes of the population at time $t$. The population, $\nu_{t}$, at time $t$ is represented by the rescaled sum of Dirac measures on $\mathcal{U}^{2}$,

[TABLE]

Formally, $\nu_{t}$ takes values in the set of re-scaled point measures

[TABLE]

on $\mathcal{U}^{2}$, equipped with the vague topology.

For a formal construction of the corresponding measure valued Markov process, see [27]. We just insist on the reproduction rate of an individual of genotype $u_{1}u_{2}$ with an individual of genotype $v_{1}v_{2}$, which takes the form $f_{u_{1}u_{2}}\frac{f_{v_{1}v_{2}}R_{u_{1}u_{2}}(v_{1}v_{2})}{K\langle\nu R_{u_{1}u_{2}},f\rangle}$, where $\nu R_{u_{1}u_{2}}$ is the population restricted to the pool of potential partners of an individual of genotype $u_{1}u_{2}$. At the individual level, this form of reproduction rate means that each individual of genotype $u_{1}u_{2}$ reproduces at a rate which depends on its genotype through its fertility $f_{u_{1}u_{2}}$; more precisely, it bears an exponentially distributed clock with parameter $f_{u_{1}u_{2}}$, and when this rings, it chooses a partner at random in the pool determined by the reproductive compatibility $R$, with a weight proportional to the fertility of the partner. The reproductive compatibility can for example be thought as a way of coding if two individuals are or not of the same species. Within a species, the reproduction rate of a pair of (compatible) individuals is given by the product of their fertilities.

We make the following Assumptions (A):

• (A1)

The functions $f,D$ and $c$ are measurable and bounded, which means that there exist constants $\bar{f},\bar{D},\bar{c}<\infty$ such that for all $u_{1}u_{2},v_{1}v_{2}\in\mathcal{U}^{2}$,

[TABLE]

• (A2)

There exists a constant $\underline{c}>0$ such that for all $u_{1}u_{2}\in\mathcal{U}^{2}$, $c_{u_{1}u_{2},u_{1}u_{2}}\geq\underline{c}$ and $f_{u_{1}u_{2}}-D_{u_{1}u_{2}}>0$,

• (A3)

There exists a function, $\bar{m}:{\mathbb{R}}\rightarrow{\mathbb{R}}_{+}$, such that $\int\bar{m}(h)dh<\infty$ and $m(u,h)\leq\bar{m}(h)$ for any $u\in\mathcal{U}$ and $h\in{\mathbb{R}}$.

For fixed $K$, under the Assumptions (A1)+(A3) and assuming that ${\mathbb{E}}(\langle\nu_{0},\mathbbm{1}\rangle)<\infty$, Fournier and Méléard [17] have shown existence and uniqueness in law of the process. For $K\rightarrow\infty$, under mild restrictive assumptions, they prove the convergence of the process in the space $\mathbb{D}({\mathbb{R}}_{+},\mathcal{M}^{K})$ of càdlàg functions from ${\mathbb{R}}^{+}$ to $\mathcal{M}_{K}$, to a deterministic process, which is the solution to a non-linear integro-differential equation. Assumption (A2) ensures that the population does not tend to infinity in finite time or becomes extinct too fast.

1.2. Previous works

Consider the process starting with a monomorphic $aa$ population, with one additional mutant individual of genotype $aA$. Assume that the phenotype difference between the mutant and the resident population is small. The phenotype difference is assumed to be a slightly smaller death rate compared to the resident population, namely

[TABLE]

for some small enough $\Delta>0$. The mutation probability for an individual with genotype $u_{1}u_{2}$ is given by $\mu_{K}$. Hence, the time until the next mutation in the whole population is of order $\frac{1}{K\mu_{K}}$. Now assume that the demographic parameters introduced in Section 1.1 depend continuously on the phenotype. In particular, they are the same for individuals bearing the same phenotype.

In [8] it is proved that if the two alleles $a$ and $A$ are co-dominant and if the allele $A$ is slightly fitter than the allele $a$, namely

[TABLE]

then in the limit of large population and rare mutations ($\ln K\ll\frac{1}{\mu_{K}K}\ll e^{VK}$ for some $V>0$), the suitably time-rescaled process converges to the TSS model of adaptive dynamics, essentially as shown in [4] in the haploid case. In particular, the genotypes containing the unfit allele $a$ decay exponentially fast after the invasion of $AA$ (see Figure 1).

If in place of co-dominance we assume, as in [27], that the fittest phenotype $A$ is dominant, namely

[TABLE]

then this has a dramatic effect on the evolution of the population and, in particular, leads to a much prolonged survival of the unfit phenotype $aa$. Indeed, it was known for some time (see e.g. [26]) that in this case the unique stable fixed point $(0,0,\bar{n}_{AA})$ corresponding to a monomorphic $AA$ population is degenerate, i.e. its Jacobian matrix has zero-eigenvalue. This implies that in the deterministic system, the $aa$ and $aA$ populations decay in time only polynomially fast to zero, namely like $1/t^{2}$ and $1/t$, respectively. This is in contrast to the exponential decay in the co-dominant scenario (see Figure 1). In [27] it was shown that the deterministic system remains a good approximation of the stochastic system as long as the size of the $aA$ population remains much larger than $K^{1/2}$ and therefore that the $a$ allele survives for a time of order at least $K^{1/2-{\alpha}}$, for any ${\alpha}>0$111The article [27] only state that survival occurs up to time $K^{1/4-{\alpha}}$. However, taking into account that it is really only the survival of the $aA$ population that needs to be ensured, one can easily improve this to $K^{1/2-{\alpha}}$.. Note that this statement is a non trivial fact, since it is not a consequence of the law of large numbers, because the time window diverges as $K$ grows. In summary, the unfit recessive $a$ allele survives in the population much longer due to the slow decay of the $aA$ population.

It is argued in [27] that if we choose the mutation time scale in such a way that there remain enough $a$ alleles in the population when a new mutation occurs, i.e.

[TABLE]

and if the new mutant can coexist with the unfit $aa$ individuals, then the $aa$ population can potentially recover. This is the starting point of the present paper.

1.3. Goal of the paper

The goal of this paper is to show that under reasonable hypotheses, the prolonged survival of the $a$ allele after the invasion of the $A$ allele can indeed lead to a recovery of the $aa$-type. To do this, we assume that there will occur a new mutant allele (on the same gene), $B$, that on the one hand has a higher fitness than the $AA$-phenotype but that (for simplicity) has no competition with the $aa$-type. The possible genotypes after this mutation are $aa,aA,AA,aB,AB$, and $BB$, so that even for the deterministic system we have now to deal with a $6$-dimensional dynamical system whose analysis if far from simple.

Under the assumption of dominance of the fittest phenotype, and mutation rate satisfying (1.7), we consider the model described in Section 1.1 starting at the time of the second mutation, that is (with probability converging to 1 as $K\to\infty$) the $AA$ population being close to its equilibrium and the $aA$ population having decreased to a size of order ${\varepsilon}_{K}=O(K\mu_{K})$, while the $aa$ population is of the order of the square of the $aA$ population. We assume that there just occurred a mutation to a fitter (and more dominant) allele $B$: we thus start with a quantity $\frac{1}{K}$ of genotype $AB$. We will start with a population where $AA$ is close to its equilibrium, the populations of $aA$ and $aa$ are already small (of order ${\varepsilon}_{K}$ and ${\varepsilon}_{K}^{2}$), and by mutation a single individual of genotype $AB$ appears.

By using well known techniques [4, 7, 8], we know that the $AB$ population behaves as a super-critical branching process and reaches a level ${\varepsilon}_{0}>0$ with positive probability in a time of order $\ln K$, without perturbing the 3-system $(aa,aA,AA)$.

We see in numerical solutions to the deterministic system that a reduced fertility together with a reduced competition between $a$ and $B$ phenotypes constitutes a sufficient condition for the recovery of the $aa$ population. For simplicity and in order to prove rigorous results, we suppose that there can be no reproduction between individuals of phenotypes $a$ and $B$, nor competition between them, and we reduce the number of remaining parameters as much as possible (see Section 2). We study the deterministic system which corresponds to the large population limit of the stochastic counterpart, and we show that (for an initial quantity ${\varepsilon}$ of $aA$, ${\varepsilon}^{2}$ of $aa$ and ${\varepsilon}^{3}$ of $AB$) the system converges to a fixed point denoted by $p_{aB}$ consisting of the two coexisting populations $aa$ and $BB$. If no further assumptions are made, we show that the number of individuals bearing an $a$ allele decreases to level ${\varepsilon}^{1+{\Delta}/(1-{\Delta})}$ (where $\Delta$ is defined in (1.4)) before $aa$ grows and stabilises at order 1.

Remark*.*

Note that our study considers an initial quantity of $aA$ individuals $n_{aA}(0)={\varepsilon}>0$, which should be thought as ${\varepsilon}={\varepsilon}_{K}=O(K\mu_{K})$ for the following discussion about the associated stochastic system. The quantity ${\varepsilon}_{K}$ satisfies $(\ln K)^{-1}\ll{\varepsilon}_{K}\ll O(K^{-\frac{1}{2}+\alpha})$ if the mutation time satisfies (1.7).

Let us discuss the (conjectural) implications of our results on the stochastic system:

If $\Delta<\frac{\alpha}{1-2\alpha}$, our control on the $a$ allele is in principle sufficient in order for the stochastic system to exhibit the recovery of $aa$ with positive probability in the large population limit. Indeed, if the mutation time is of order $K^{\frac{1}{2}-\alpha}$, then the initial amount of $aA$ and $aa$ genotypes is close to the typical fluctuations of those populations, that is respectively of orders $O(K^{-\frac{1}{2}+\alpha})$ and $O(K^{-1+2\alpha})$. Following the heuristics of [27] (although the six-dimensional stochastic process is surely much more tedious to study), the deterministic system should constitute a good approximation of the process if the typical fluctuations of populations containing an $a$ allele do not bring them to extinction. If $\Delta<\frac{\alpha}{1-2\alpha}$ this ensures that the population containing an $a$ allele is not falling below order $K^{-1/2}$ at any time.

In order to go deeper and control the speed of recovery of the $aa$ population, we look for a parameter regime which ensures that the $aa$ population always grows after the invasion of $B$. Ensuring this lower bound on $aa$ is not trivial at all, and the solution we found is to introduce an additional parameter $\eta$, which lowers the competition between the $aA$ and $BB$ populations, compared to the one between $AA$ and $BB$. Note that the competition does not depend only on the phenotype, and can be interpreted as a refinement of a phenotypic competition for resources: the strength (or ability to get resources) of an individual not only depends on its phenotype but also on the dominance of its genotype. A biological interpretation for this kind of competition could be that it is coded in the alleles which food an individual with a given genotype prefers. We show that for $\eta$ larger than some positive value (of order $\Delta$), the $aa$ population always grows after the invasion of $B$. The time of convergence to the coexistence fixed point is thus lowered, see Figure 5. Moreover, we point out the existence of a bifurcation: for $\eta$ larger than some threshold, the co-existence fixed point $p_{aB}$ becomes unstable and the system converges to another fixed point where all populations coexist.

Our contribution is a rigorous description of a mechanism by which a recessive allele can re-emerge in a population. This can be seen as a statement of genetic robustness exhibited by diploid populations performing sexual reproduction.

The structure of the paper is the following. In Section 2 we describe our assumptions on the parameters of the model, and compute the large population limit; in Section 3 we present our results on the evolution of the deterministic system towards the co-existence fixed point $p_{aB}$, and we give a heuristic of the proof. Section 4 is dedicated to the proof of these results. The closing Section 5 contain a heuristic considerations and numerical simulations of the model with relaxed assumption on the parameters.

Notation.

We write $x={\Theta}(y)$ whenever $x=O(y)$ and $y=O(x)$ as ${\varepsilon}\to 0.$

2. Model setup

Let $\mathcal{G}=\{aa,aA,AA,aB,AB,BB\}$ be the genotype space. Let $n_{i}(t)$ be the number of individuals with genotype $i\in\mathcal{G}$ in the population at time $t$ and set $n^{K}_{i}(t)\equiv\frac{1}{K}n_{i}(t)$.

Definition 2.1**.**

The equilibrium size of a monomorphic $uu$ population, $u\in\{a,A,B\}$, is the fixed point of a 1-dimensional Lotka-Volterra equation and is given by

[TABLE]

Definition 2.2**.**

For $u,v\in\{a,A,B\}$, we call

[TABLE]

the invasion fitness of a mutant $uv$ in a resident $uu$ population.

We take the phenotypic viewpoint and assume that the $B$ allele is the most dominant one. That means the ascending order of dominance (in the Mendelian sense) is given by $a<A<B$, i.e.

1. (1)

phenotype $a$ consists of the genotype $aa$, 2. (2)

phenotype $A$ consists of the genotypes $aA,AA$, 3. (3)

phenotype $B$ consists of the genotypes $aB,AB,BB$.

For simplicity, we assume that the fertilities are the same for all genotypes, and that natural death rates are the same within the three different phenotypes. Moreover, we assume that there can be no reproduction between $a$ and $B$ phenotypes.

To sumarize, we make the following Assumptions (B) on the rates:

• (B1)

Fertilities. For all $i\in\mathcal{G}$, and some $f>0$

[TABLE]

• (B2)

Natural death rates. The difference in fitness of the three phenotypes is realised by choosing a slightly higher natural death-rate of the $a$-phenotype and a slightly lower death-rate for the $B$-phenotype. For some $0<{\Delta}<D$,

[TABLE]

• (B3)

Competition rates. We require that phenotypes $a$ and $B$ do not compete with each other. Moreover, we introduce a parameter $\eta\geq 0$ which lowers the competition between $BB$ and $aA$. For some $0\leq\eta<c$,

[TABLE]

Remark*.*

If $\eta>0$, the competition does not depend only on the phenotype, and can be interpreted as a refinement of a phenotypic competition for resources: the strength (or ability to get resources) of an individual not only depends on its phenotype but also on its genotype. Genetically, it makes sense to assume that (positive) competition rates are decreasing in the ”genetic distance” between two individuals. This is the case for the above competition matrix if we assume that the mutations $a\to A\to B$ can only occur successively. Indeed, the Hamming distance between genotypes is then the following graph distance, where $d(aA,BB)>d(AA,BB)$.

• (B4)

Reproductive compatibility. We require that phenotypes $a$ and $B$ do not reproduce with each other.

[TABLE]

Observe that, under Assumptions (B),

[TABLE]

Therefore, the mutant $AB$ has a positive invasion fitness in the population $AA$, as well as $aa$ in the $BB$ population (due to the absence of competition between them).

2.1. Birth rates

Since we assume that there is no recombination between phenotypes $a$ and $B$. Thus,

1. (1)

the pool of possible partners for the phenotype $a$ consists of phenotypes $a$ and $A$; the total population of this pool is denoted by

[TABLE] 2. (2)

the pool of possible partners for the phenotype $A$ consists of the three phenotypes $a$, $A$, and $B$; the total population of this pool is denoted by

[TABLE] 3. (3)

the pool of possible partners for the phenotype $B$ consists of phenotypes $A$ and $B$; the total population of this pool is denoted by

[TABLE]

Computing the reproduction rates with the Mendelian rules as described in [27] leads to the following (time-dependent) birth-rates $b_{i}=b_{i}(n(t))$:

[TABLE]

2.2. Death rates

The death rates are the sum of the natural and competition death rates:

[TABLE]

2.3. Large population limit

By [17] or [8], for large populations, the behaviour of the stochastic process is close to the solution of a deterministic equation.

Proposition 2.3** (Theorem 2.1 in [14]).**

**

*Let $T>0$ and $C\subset{\mathbb{R}}_{+}^{6}$ be a compact set. Assume that the initial condition $n^{K}(0)=\frac{1}{K}(n_{aa}(0),n_{aA}(0),n_{AA}(0),n_{aB}(0),n_{AB}(0),n_{BB}(0))$ converges almost surely to a deterministic vector $x^{0}=(x^{0}_{1},x^{0}_{2},x^{0}_{3},x^{0}_{4},x^{0}_{5},x^{0}_{6})\in C$, as $K\to\infty$.

Let $\tilde{n}(t,x^{0})$ denote the solution to*

[TABLE]

with initial condition $x^{0}$, where $(b_{i})_{i\in{\mathcal{G}}}$ and $(d_{i})_{i\in{\mathcal{G}}}$ are given in (2.13)-(2.18) and (2.19)-(2.24). Then, for all $T>0$,

[TABLE]

for all $i\in{\mathcal{G}}$.

2.4. Initial condition

Fix ${\varepsilon}>0$ sufficiently small. For the results below, we will consider the dynamical system (2.25) starting with the initial condition:

[TABLE]

which corresponds to the long time behaviour of the dynamical system considered in [27] plus a quantity ${\varepsilon}^{3}$ of the new mutant $AB$.

Remark*.*

In all the figures below, the choice of parameters is the following:

[TABLE]

and the parameter $\eta$ is specified on each picture.

3. Results

We are working with a 6-dimensional dynamical system, and computing all the fixed points analytically is impossible for a general choice of the parameters. We can, however, compute those which are relevant for our study. We will call $p_{A}$ (resp. $p_{B}$) the fixed points corresponding to the monomorphic $AA$ (resp. $BB$) population at equilibrium, and $p_{aB}$ the fixed point corresponding to the coexisting $aa$ and $BB$ populations. Setting the relevant populations to 0 and solving $\dot{n}(t)=0$, we get:

[TABLE]

where $\bar{n}_{a}=\frac{f-D-\Delta}{c},\bar{n}_{A}=\frac{f-D}{c},$ and $\bar{n}_{B}=\frac{f-D+\Delta}{c}.$ Note that the $BB$ equilibrium population is the same in $p_{B}$ and $p_{aB}$. This is due to the non-interaction between phenotypes $a$ and $B$.

Our general result is that starting with initial conditions (2.28)-(2.33), that is close to $p_{A}$ (with small coordinates in directions $aa$, $aA$ and $AB$), and under minimal assumptions on the parameters, the system gets very close to $p_{B}$ before finally converging to $p_{aB}$, see Figure 3.

Theorem 3.1**.**

Consider the dynamical system (2.25) started with initial conditions (2.28)-(2.33). Suppose the following Assumptions (C) on the parameters hold:

1. (C1)

$\Delta$* sufficiently small,* 2. (C2)

$f$* sufficiently large,* 3. (C3)

$0\leq\eta<c/2$.

*Then the system converges to the fixed point $p_{aB}$. More precisely, for any fixed $\delta>0$, as ${\varepsilon}\to 0$, it reaches a $\delta$-neighbourhood of $p_{aB}$ in a time of order ${\Theta}({\varepsilon}^{-1/(1+\eta\bar{n}_{B}-{\Delta})})$.

Moreover:*

1. (1)

for $\eta=0$, the amount of allele $a$ in the population decays to ${\Theta}({\varepsilon}^{1+{\Delta}/(1+{\Delta})})$ before reaching ${\Theta}(1)$, 2. (2)

for $\eta>\frac{4{\Delta}}{\bar{n}_{B}}$, the amount of $a$ allele in the population is bounded below by ${\Theta}({\varepsilon})$ for all $t>0$.

Remark*.*

For $\eta$ large, we prove that the fixed point $p_{aB}$ is unstable. We observe numerically that the system is attracted to a fixed point where all the 6 populations coexist, but we do not prove this.

Let us now briefly discuss the linear stability of the relevant fixed points and give a heuristics of the proof of Theorem 3.1.

3.1. Linear stability analysis

The Jacobian matrix $J_{F}:=(\partial F_{i}/\partial n_{j})_{ij}$ of the map $F$ defined in (2.25) can be explicitly computed at $p_{A}$ and $p_{aB}$ and the situation is as follows:

• •

The eigenvalues of $J_{F}(p_{A})$ are $0,\Delta>0$ and $-(f-D),-(f+\Delta),-(f-\Delta)$ (double) which are all strictly negative under Assumptions (C). The fixed point $p_{A}$ is thus unstable.

• •

The eigenvalues of $J_{F}(p_{aB})$ are [math] (double), and $-(2f-D),-(f-D+\Delta),-(f-D-\Delta),-((f-D)(5f-4D)+f\Delta)/(4(f-D)+\eta\bar{n}_{B})$ which are strictly negative under Assumptions (C). The linear analysis thus does not imply the stability of $p_{aB}$ but the Phase 4 of the proof does (see Section 4.5) .

It turns out that $J_{F}(p_{B})$ is singular but as the invasion fitness of $aa$ is positive, i.e. $S_{aa,BB}>0$ (see (2.7)), this implies that a small perturbation in the first coordinate will be amplified, and thus implies the instability of the fixed point $p_{B}$.

3.2. Heuristics of the proof

Recall we start the dynamical system (2.25) with initial conditions (2.28)-(2.33). A numerical solution of the system is provided on Figure 4.

Remark*.*

Assumption C1 of Theorem 3.1 is needed throughout the proof in order to be able to use the results of [27] which rely on the Center Manifold Theorem (a line of fixed points becomes an invariant line under small enough perturbation).

1. Phase 1.

Time period: until $n_{AB}={\varepsilon}_{0}$.

The mutant population, consisting of all individuals of phenotype $B$, first grows up to ${\varepsilon}_{0}$ exponentially fast with rate $\Delta$ without perturbing the behaviour of the 3-system $(aa,aA,AA)$. The rate of growth corresponds to the invasion fitness of $AB$ in the resident population $AA$, see (2.7). Following [27], $AA$ stays close to $\bar{n}_{A}$, while $aA$ and $aa$ continue to decay like $1/t$ and $1/t^{2}$ respectively. The duration $T_{1}$ of this phase is such that ${\Theta}({\varepsilon}^{3})e^{t{\Delta}}={\Theta}(1)\Leftrightarrow T_{1}={\Theta}(|\log{\varepsilon}|)$. 2. Phase 2.

Time period: until $n_{aA}={\Theta}(n_{AA})$.

The evolution is a perturbation of an effective 3-system $(AA,AB,BB)$ which behaves exactly the same as in [27], since the parameters satisfy the same hypotheses (slightly lower death rate for phenotype $B$ than for phenotype $A$, and constant competition parameters). A comparison result (following Theorem 4.6 below) shows that this 3-system is almost unperturbed until $n_{aA}={\Theta}(n_{AA})$. If that happens in a time $T_{2}$ diverging with ${\varepsilon}$ (which we ensure throughout the calculation), we thus know that $BB$ approaches $\bar{n}_{B}$, while $n_{AB}\propto 1/t$ and $n_{AA}\propto 1/t^{2}$.

The important fact in this phase is that the amount of allele $a$ in the population decays for $\eta$ small while it increases for large enough $\eta$. Indeed, let us derive some bounds on $\Sigma_{aA,aB}=n_{aA}+n_{aB}$. The population $\Sigma_{aA,aB}$ reproduces by taking the dominant allele in a population of order ${\Theta}(1)$ and the allele $a$ in itself. Thus its birth rate satisfies $b_{\Sigma_{aA,aB}}\approx f\Sigma_{aA,aB}$. We can compute its death rate exactly and use that $n_{BB}\approx\Sigma_{5}\approx\bar{n}_{B}$:

[TABLE]

The last equality comes from the fact that $aA$ newborns have mainly their $a$ allele coming from ${\Sigma}_{aA,aB}$ and their $A$ allele coming from $AB$. Using the $1/t$ decay of $AB$ we get:

[TABLE]

As ${\Sigma}_{aA,aB}(T_{1})={\Theta}({\varepsilon})$ we deduce that ${\Sigma}_{aA,aB}(t)={{\Theta}({\varepsilon})}({\Theta}(1)+{\Theta}(1)t)^{{\Theta}(\eta\bar{n}_{B}-{\Delta})}$, and thus $n_{aA}={\Theta}(n_{AB}\cdot{\Sigma}_{aA,aB})={{{\Theta}({\varepsilon})}({\Theta}(1)+{\Theta}(1)t)^{{\Theta}(\eta\bar{n}_{B}-{\Delta})}}/{({\Theta}(1)+{\Theta}(1)t)}$. By solving $n_{aA}={\Theta}(n_{AA})={\Theta}(n_{AB}^{2})$ we get the order of magnitude of $T_{2}={\Theta}\left({\varepsilon}^{-1/(1+\eta\bar{n}_{B}-{\Delta})}\right)$. Note that for $\eta=0$, ${\Sigma}_{aA,aB}(T_{2})={\Theta}({\varepsilon}^{1+{{\Delta}}/{(1-{\Delta})}})$. Moreover, (3.5) implies that for $\eta>\Delta/\bar{n}_{B}$, we have $\dot{\Sigma}_{aA,aB}>0$, which proves points 1 and 2 of Theorem 3.1. 3. Phase 3.

Time period: until ${aa}$ reaches equilibrium.

The fact that $n_{aA}={\Theta}(n_{AA})$ has a crucial effect on the birth rate of $aa$ (see (2.13)) since the term $(n_{aa}+\frac{1}{2}n_{aA})/({n_{aa}+n_{aA}+n_{AA}})$ becomes of order ${\Theta}(1)$. As long as $AA$ stays smaller than ${\Theta}({\varepsilon})$, we get a lower bound on $n_{aa}$ which grows exponentially fast since $f$ is chosen large enough (Assumption C2):

[TABLE]

As $aa$ grows, it makes ${\Sigma}_{aA,aB}$ grow, and thus $AA$ and $AB$ as well. We have to show that this could not prevent $aa$ from reaching equilibrium. We do not give a detailed argument here, but essentially, the presence of the macroscopic $BB$ population prevents all the non-$aa$ populations to grow too much. Note that if $\eta$ is too large, then $aA$ could get a positive fitness and grow to a macroscopic level. That is why we have to impose Assumption C3, which will become clearer heuristically in the next phase. We recall that $aa$ does not compete with $BB$ and thus it grows exponentially fast with rate ${f-(D+\Delta)}$ until an ${\varepsilon}_{0}$-neighbourhood of the fixed point where $aa$ and $BB$ coexist. The rate of growth corresponds to the invasion fitness of $aa$ in the resident population $BB$, see (2.7). Note that, due to Assumption C2, this rate is much larger than the invasion rate of $BB$ into $AA$. That is why the fourth phase looks very steep on Figure 4, see the stretched version on Figure 6. This phase lasts a time $T_{3}={\Theta}(|\log{\varepsilon}|)$. 4. Phase 4.

The Jacobian matrix of the field (2.25) at the fixed point $p_{aB}$ has two zero, and 4 negative eigenvalues. $p_{aB}$ is thus a non-hyperbolic equilibrium point of the system and linearisation fails to determine its stability properties. Instead, we use the result of center manifold theory ([20, 28]) that asserts that the qualitative behaviour of the dynamical system in a neighbourhood of the non-hyperbolic critical point $p_{aB}$ is determined by its behaviour on the center manifold near $p_{aB}$. Using the Center Manifold Theorem, we show that asymptotically as $f\to\infty$, the field is attractive for $\eta<c\cdot r_{max}$ where $r_{max}\simeq 0.593644$ is the maximum of the rational function (4.306). Thus $p_{aB}$ is a stable fixed point which is approached with speed $\frac{1}{t}$ as long as $\eta<c\cdot r_{max}$. For higher values of $\eta$, numerical solutions show that the system converges to a fixed point where the 6 populations co-exist, but we do not prove this.

4. Proof

Definition 4.1**.**

Let $x,y,z\in\{aa,aA,AA,aB,AB,BB\}$ and $h\in{\mathbb{R}}$. We define

[TABLE]

Moreover, let

[TABLE]

The value ${\varepsilon}_{0}$ is the small order 1 level in the Phase 1, see the proof heuristics (Section 3.2). We consider $\Delta$ fixed and sufficiently small, and will first send ${\varepsilon}\to 0$ and then ${\varepsilon}_{0}\to 0$. A summary of the proof structure can be found in the Appendix A together with the implications between intermediary results.

4.1. Preliminaries

We first state an elementary fact that will be useful throughout the following analysis.

Lemma 4.2**.**

*Let $T>0$ and $f:[0,T]\to{\mathbb{R}}$ be a differentiable function such that $f(0)\leq 0$. If, for all $t\in[0,T]$, $\left\{f(t)=0\Rightarrow\dot{f}(t)<0\right\}$ then $f(t)<0$, for all $t\in(0,T]$. *

In the sequel, a part of our majorations and comparison results will rely on this lemma, which ensures that the trajectories of two dynamical systems do not cross as long as the derivative of the upper trajectory is larger than the derivative of the lower one, where these trajectories meet. We will often encounter the situation where $c,T>0$ and $n,g:[0,T]\to{\mathbb{R}}$ are two differentiable functions such that $cn(0)\leq g(0)$.

Then, Lemma 4.2 implies that if, for all $t\in[0,T]$, $\left\{cn(t)=g(t)\Rightarrow c\dot{n}(t)<\dot{g}(t)\right\}$ then $cn(t)\leq g(t)$, for all $t\in[0,T]$.

Proposition 4.3**.**

We have the following ordering relations:

1. (1)

If $n_{AA}(0)\leq\bar{n}_{A}$ then $n_{AA}(t)\leq\bar{n}_{A}$, for all $t>0$. 2. (2)

If $n_{BB}(0)<\bar{n}_{B}$ then $n_{BB}(t)\leq\bar{n}_{B}$, for all $t>0$. 3. (3)

If $n_{aB}(0)<n_{AB}(0)$ then $n_{aB}(t)\leq n_{AB}(t)$, for all $t>0$.

Proof.

1. (1)

We use Lemma 4.2 and show that $\dot{n}_{AA}<0$ at $n_{AA}=\bar{n}_{A}$. For this we construct a majorising process on $n_{AA}$ by comparing the birth and the death rates. Since $\Sigma_{6}\geq\Sigma_{5}$ we get:

[TABLE]

At $n_{AA}=\bar{n}_{A}$ we have

[TABLE]

since ${\textstyle{f\over 4\Sigma_{5}}}n_{aA}$ and ${\textstyle{f\over 4\Sigma_{5}}}n_{AB}$ are smaller than ${\textstyle{f\over 4}}$. 2. (2)

We proceed similarly using Lemma 4.2. We show that $\dot{n}_{BB}<0$ at $n_{BB}=\bar{n}_{B}$ by constructing a majorising process on $n_{BB}$:

[TABLE]

At $n_{BB}=\bar{n}_{B}$ we have

[TABLE]

since ${\textstyle{f\over 4\Sigma_{5}}}n_{aB}$ and ${\textstyle{f\over 4\Sigma_{5}}}n_{AB}$ are smaller than ${\textstyle{f\over 4}}$ and $\eta\leq c$. 3. (3)

Intuitively, this inequality comes from the fact that phenotype $a$ individuals cannot reproduce with phenotype $B$. Indeed, if we consider the couples that could give rise to an $AB$ (resp. $aB$) individual, they are of the form $(Ag_{1},Bg_{2})$ (resp. $(ag_{1},Bg_{2})$), with $g_{1},g_{2}\in\{a,A,B\}$ and the combination $(AA,Bg_{2})$ is possible, whereas $(aa,Bg_{2})$ is impossible. Here is the rigorous derivation of the result: We compare the birth- and the death-rates of $n_{AB}$ and $n_{aB}$, introduced in (2.17),(2.23), and (2.16),(2.22) respectively:

[TABLE]

We see that the death-rates of the two populations are the same, whereas the birth-rates differ only in a factor which comes from the reproduction of the other populations. If we take a closer look to these factors $I_{aB},I_{AB}$ under the assumption that $n_{aB}=n_{AB}$ we see that

[TABLE]

Thus $I_{AB}>I_{aB}$. Hence, $\dot{n}_{AB}>\dot{n}_{aB}$ and $n_{AB}(t)$ stays above $n_{aB}(t)$ for all $t>0$.

∎

4.2. Phase 1: Perturbation of the 3-system $(aa,aA,AA)$ until $AB$ reaches ${\Theta}(1)$

**

We start with initial conditions given by (2.28)-(2.33). We will show that the mutant population, consisting of all individuals of phenotype $B$, grows up to some ${\varepsilon}_{0}>{\varepsilon}$ without perturbing the behaviour of the 3-system $(aa,aA,AA)$ in this time. Let

[TABLE]

Proposition 4.4**.**

With the initial conditions (2.28)-(2.33),

[TABLE]

Proof.

Until $T_{1}$ the perturbation of the dynamics of the 3-system $(aa,aA,AA)$ is at most of order ${\varepsilon}_{0}$, by continuity of the solutions of an ODE with respect to initial conditions. We thus have $\bar{n}_{A}-{\Theta}({\varepsilon}_{0})\leq n_{AA}(t)\leq\bar{n}_{A}+{\Theta}({\varepsilon}_{0})$, as well as $n_{aa},n_{aA}\leq{\Theta}({\varepsilon}_{0})$. With this rough bounds we will find finer bounds. Heuristically, the newborns of genotype $aA$ are still in majority produced by recombination of $AA$ and $aA$, because the mutant population is not large enough to contribute. The newborns of genotype $aB$ are in majority produced by reproduction of the $aA$ population with the $B$ population. Finally, the newborns of genotype $aa$ are in majority produced by recombination of $aA$ and $aA$, because the only mutant population which could perturb their dynamics is $aB$, which is of smaller order.

1. (1)

We show that $n_{BB}\leq n_{AB}^{2}$. Indeed, we use Lemma 4.2, note that $0=n_{BB}(0)\leq n_{AB}(0)^{2}$, and prove that when $n_{BB}=n_{AB}^{2}$ we have

[TABLE]

Indeed, using Proposition 4.3 we can bound (4.21) by

[TABLE]

Thus $T_{1}=T_{{\varepsilon}_{0}}^{aa+aA+aB+AB}$. 2. (2)

We show that $n_{aB}\leq n_{aA}n_{AB}$. Indeed, we use Lemma 4.2, note that $0=n_{aB}(0)\leq n_{aA}(0)n_{AB}(0)$, and prove that when $n_{aB}=n_{aA}n_{AB}$ we have

[TABLE]

Indeed, using (1) we can bound (4.23) by

[TABLE]

Thus $T_{1}=T_{{\varepsilon}_{0}}^{aa+aA+AB}$. 3. (3)

We show that $n_{aa}\leq n_{aA}^{2}$. Indeed, we use Lemma 4.2, note that $n_{aa}(0)\leq n_{aA}(0)^{2}$, and prove that when $n_{aa}=n_{aA}^{2}$ we have

[TABLE]

Indeed, using (1) and (2), we can bound (4.26) by

[TABLE]

Thus $T_{1}=T_{{\varepsilon}_{0}}^{aA+AB}$. 4. (4)

Now, we show that $T_{1}=T_{{\varepsilon}_{0}}^{AB}$. First, we estimate the maximal time $n_{AB}$ would need to reach the level ${\varepsilon}_{0}$. For this upper bound on the time $T_{{\varepsilon}_{0}}^{AB}$, we have to construct a minorising process for $n_{AB}$. Indeed, let us compare the birth and death rates using the results (1)-(3):

[TABLE]

Hence, we get for the minorising process

[TABLE]

and the time $T_{{\varepsilon}_{0}}^{AB}$ is at most of order ${\Theta}\left(\log(({{\varepsilon}_{0}}/{{\varepsilon}^{3}})^{\frac{1}{\Delta-{\Theta}({\varepsilon}_{0})}}\right)$.

In a second step we show that $n_{aA}(t)\leq{\varepsilon}e^{{\Theta}({\varepsilon}_{0})t}$ and that the minimal time $n_{aA}$ would need to reach ${\varepsilon}_{0}$ would be bigger than ${\Theta}\left(\log(({{\varepsilon}_{0}}/{{\varepsilon}^{3}})^{\frac{1}{\Delta-{\Theta}({\varepsilon}_{0})}}\right)$. For this we construct a majorising process on $n_{aA}$. The birth rate and the death rate can be bounded by using (1)-(3):

[TABLE]

and we get that

[TABLE]

Thus the minimal time $n_{aA}$ would need to reach ${\varepsilon}_{0}$ is ${\Theta}\left(\log({\varepsilon}_{0}/{\varepsilon})^{{\textstyle{1\over{\Theta}({\varepsilon}_{0})}}}\right)$ which is bigger than the maximal time $n_{AB}$ would need to reach the level ${\varepsilon}_{0}$. Thus $T_{1}=T_{{\varepsilon}_{0}}^{AB}$. 5. (5)

It is left to show that $T_{1}\geq{\Theta}\left(\log(({{\varepsilon}_{0}}/{{\varepsilon}^{3}})^{\frac{1}{\Delta-{\Theta}({\varepsilon}_{0})}}\right)$. For this lower bound on the time $T_{1}$, we have to construct a majorising process for $n_{AB}$:

[TABLE]

Hence, we get for the majorising process

[TABLE]

and the time $T_{1}$ is at least of order ${\Theta}\left(\log(({{\varepsilon}_{0}}/{{\varepsilon}^{3}})^{\frac{1}{\Delta-{\Theta}({\varepsilon}_{0})}}\right)$.

∎

Corollary 4.5**.**

With the initial conditions (2.28)-(2.33), for all $t\in[0,T_{1}]$, the following statements hold:

1. (1)

$n_{AB}$* grows exponentially with rate ${\Delta}$. It reaches the level ${\varepsilon}_{0}$ in a time at most of order ${\Theta}\left(\log\left(({{\varepsilon}_{0}}/{{\varepsilon}^{3}})^{\frac{1}{\Delta-{\Theta}({\varepsilon}_{0})}}\right)\right)$.* 2. (2)

$n_{BB}\leq n_{AB}^{2}$, 3. (3)

$n_{aB}\leq{\Theta}({\varepsilon}^{1-{\Theta}({\varepsilon}_{0})}{\varepsilon}_{0})$, $n_{aA}\leq{\Theta}({\varepsilon}^{1-{\Theta}({\varepsilon}_{0})})$, $n_{aa}\leq{\Theta}({\varepsilon}^{2-{\Theta}({\varepsilon}_{0})})$ and

$\bar{n}_{A}-{\Theta}({\varepsilon}_{0})\leq n_{AA}\leq\bar{n}_{A}+{\Theta}({\varepsilon}_{0})$.

Proof.

1. (1)

is a direct consequence of Proposition 4.4, see proof of points (4) and (5). 2. (2)

is a direct consequence of Proposition 4.4, see proof of points (1). 3. (3)

follows from the proof of Proposition 4.4 together with $\bar{n}_{A}-{\Theta}({\varepsilon})\leq\Sigma_{5}\leq\bar{n}_{A}+2{\Delta}{\varepsilon}_{0}$, which we prove now.

Using the results of Proposition 4.4 we construct a minorising and a majorising processes on $\Sigma_{5}$:

[TABLE]

We use Lemma 4.2: at time 0 the bounds are fulfilled; moreover, at the upper bound we have $\dot{\Sigma}_{5}<0$ and at the lower bound $\dot{\Sigma}_{5}>0$, which ensures the claimed bounds.

∎

Note that Corollary 4.5 implies that

[TABLE]

4.3. Phase 2: Perturbation of the 3-system $(AA,AB,BB)$ until $n_{aA}={\Theta}(n_{AA})$

**

The initial conditions at the beginning of the second phase are:

[TABLE]

Let $\delta>0$ (to be chosen sufficiently small in the sequel). Let

[TABLE]

As the process stays uniformly bounded in time, observe that until $T_{2}$ we automatically have

[TABLE]

We will show that for $t\in[T_{1},T_{2}]$ the system behaves as a main 3-system $(AA,AB,BB)$ plus perturbations of order $\delta$. The 3-system $(AA,AB,BB)$ behaves exactly the same as in [27] since the parameters satisfy the same hypotheses (slightly lower death rate for phenotype $B$ than for phenotype $A$ individuals, and constant competition parameters).

Moreover, the crucial role of the parameter $\eta$ is that the population containing an allele $a$ only continues to grow in this phase when $\eta$ is large enough. This is due to the smaller competition that $aA$ feels from $BB$, the $aA$ population is thus higher and induces the growth of $aB$.

We start by considering how the growth of $aa,aA$- and $aB$ populations can perturb the 3-system $(AA,AB,BB)$.

Lemma 4.6**.**

Let $n_{.}^{up}(t)$ be the population of the unperturbed 3- system $(AA,AB,BB)$, that is the solution to (2.25) with $n_{aa}=n_{aA}=n_{aB}=0$. The 3-system $(AA,AB,BB)$ satisfies

[TABLE]

Proof.

We consider the rates of $AA,AB$ and $BB$ under the perturbation of $aa,aA$ and $aB$:

[TABLE]

Thus,

[TABLE]

For the $AB$ population we get:

[TABLE]

And finally for the $AA$ population:

[TABLE]

∎

As solutions of a dynamical system are continuous with respect to its parameters (in particular with respect to $\delta$), the latter lemma shows that until $T_{2}$, the 3-system $(AA,AB,BB)$ is at most perturbed by ${\Theta}(\delta)$. We will show that $T_{2}$ diverges with ${\varepsilon}$. Thus, for small enough ${\delta}$, $AB$ will have time to reach the small fixed value $\sqrt{{\varepsilon}_{0}}>0$ in this phase, and we can use the asymptotic decay $1/t$ of the $AB$ and $1/t^{2}$ of the $AA$ populations, which is proved in [27]. We now start to analyse the growth of the small $aa$-, $aA$- and $aB$ populations. The sum-process $\Sigma_{5}$ plays a crucial role for the behaviour of the system in this phase and we need finer bounds on it.

Remark*.*

In the sequel we often use that $\Sigma_{5}\leq\Sigma_{6}\leq\Sigma_{5}+{\Theta}({\delta})$, which follows from (4.55) and Lemma 4.6.

Proposition 4.7**.**

The sum-process $\Sigma_{5}=n_{aA}+n_{AA}+n_{aB}+n_{AB}+n_{BB}$ satisfies for all $t\in[T_{1},T_{2}]$:

[TABLE]

Proof.

We estimate a minorising process and a majorising process on $\Sigma_{5}$:

[TABLE]

We get

[TABLE]

We start with the proof of the upper bound. We use Lemma 4.2 and show that when $\Sigma_{5}$ reaches the upper-bound, it decays faster than the latter. Using (4.79) we compute $\dot{\Sigma}_{5}$ at the bound. Note that if $\Sigma_{5}=\bar{n}_{B}-{\textstyle{{\Delta}\over c\bar{n}_{B}}}n_{AA}+{\textstyle{{\Delta}^{2}\over c\bar{n}_{B}}}n_{AA}$, then $\Sigma_{5}^{2}=\bar{n}_{B}^{2}-{\textstyle{2{\Delta}\over c}}n_{AA}+{\textstyle{{\Delta}^{2}\over c^{2}\bar{n}_{B}^{2}}}n_{AA}^{2}+{\textstyle{2{\Delta}^{2}\over c}}n_{AA}+{\Theta}({\Delta}^{4})n_{AA}^{2}$, thus

[TABLE]

It is left to show that $\dot{\Sigma}_{5}\leq-{\textstyle{{\Delta}\over c\bar{n}_{B}}}\dot{n}_{AA}+{\textstyle{{\Delta}^{2}\over c\bar{n}_{B}}}\dot{n}_{AA}$. Since we already know (cf. Lemma 4.6) that $(AA,AB,BB)$ behaves like a 3-system with ${\Theta}(\delta)$ perturbations, then $\dot{n}_{AA}\leq 0$, this finishes the proof of the upper bound.

Now we check the lower bound. If $\Sigma_{5}=\bar{n}_{B}-\frac{{\Delta}}{c\bar{n}_{B}}n_{AA}-{\textstyle{{\Delta}^{2}\over c\bar{n}_{B}}}n_{AA}$ then $\Sigma_{5}^{2}=\bar{n}_{B}^{2}-{\textstyle{2{\Delta}\over c}}n_{AA}-{\textstyle{{\Delta}^{2}\over c^{2}\bar{n}_{B}^{2}}}n_{AA}-{\textstyle{2{\Delta}^{2}\over c}}n_{AA}$. Using (4.80), the derivative of $\Sigma_{5}$ at the lower bound is thus bounded by

[TABLE]

For the second inequality we use that $n_{AA}\leq\bar{n}_{A}$ (Proposition 4.3). By (4.74) and Lemma 4.2, it is enough to show that at the lower bound $\dot{\Sigma}_{5}\geq-{\textstyle{{\Delta}\over c\bar{n}_{B}}}\dot{n}_{AA}$. Since $AA$ is decreasing we have to calculate a majorising process on $AA$:

[TABLE]

Hence we have to show that the slope of $\Sigma_{5}$ at the lower bound is bigger than the one of the lower bound, namely we show that ${\Delta}^{2}n_{AA}\left(1-{\textstyle{1+{\Theta}({\Delta})\over c\bar{n}_{B}}}\right)-{\Theta}({\delta})\geq{\textstyle{{\Delta}f\over c\bar{n}_{B}\bar{n}_{A}}}\left(n_{AA}n_{BB}-{\textstyle{1\over 4}}n_{AB}^{2}\right)-{\Theta}({\delta}{\Delta})$, in the case $n_{AA}n_{BB}>{\textstyle{1\over 4}}n_{AB}^{2}$. This is equivalent to show that $\chi:=n_{AA}n_{BB}-{\textstyle{1\over 4}}n_{AB}^{2}\leq{\textstyle{{\Delta}\bar{n}_{A}\over f}}\left(c\bar{n}_{B}-1\right)n_{AA}$. For this we use once again Lemma 4.2 and estimate the derivative of $\chi$ from above with the help of minorising processes on $AA$ and $BB$ and a majorising process on $AB$. For bounding the death rates we use Proposition 4.3:

[TABLE]

The derivative is given by:

[TABLE]

At the upper bound we get:

[TABLE]

It is left to show that at the upper bound $\dot{\chi}\leq{\textstyle{{\Delta}\bar{n}_{A}\over f}}(c\bar{n}_{B}-1)\dot{n}_{AA}$. Using the minorising process $\dot{n}_{AA}\geq-{\Delta}n_{AA}-{\textstyle{f\over\bar{n}_{A}}}\chi-{\Theta}({\delta})$ (see (4.88)), (4.3) and Proposition 4.3 we show that

[TABLE]

This finishes the proof of the lower bound. ∎

Lemma 4.8**.**

For $t\in[T_{1},T_{2}]$ and for ${\Delta}$ sufficiently small,

[TABLE]

Proof.

Using Proposition 4.7, we have the following bound on the process:

[TABLE]

∎

Lemma 4.9**.**

For all $t\in[T_{1},T_{2}]$ the $aa$ population is bounded by

[TABLE]

Observe that this implies $T_{2}=T^{aA=\delta AA}\wedge T^{aB=\delta AB}$.

Proof.

First observe that the inequality is satisfied at $t=T_{1}$. We start with the upper bound and show that $n_{aa}$ would decrease at this bound. For this we estimate a majorising process on $aa$:

[TABLE]

We calculate the slope of this process at the upper bound using Proposition 4.7 and (4.55):

[TABLE]

By Lemma 4.2, to ensure that (4.102) stays an upper bound it is enough to show that the right-hand side satisfies

[TABLE]

This is a consequence of Lemma 4.8.

For the lower bound we proceed similarly. Since ${\textstyle{1\over\Sigma}}_{6}={\textstyle{1\over\Sigma}}_{5}-{\textstyle{n_{aa}\over\Sigma_{5}\Sigma_{6}}}$ and with the knowledge of the upper bound, we estimate a minorising process on $aa$:

[TABLE]

With Proposition 4.7 we see that at the lower bound the process increases:

[TABLE]

By Lemma 4.2, it is left to show that $\dot{n}_{aa}\geq{\textstyle{f\over 8\bar{n}_{B}(f+{\Delta}+{\Theta}({\delta}))}}\dot{\Sigma}_{aA,aB}\Sigma_{aA,aB}$ at the lower bound. Thus we have to calculate a majorising process on $\Sigma_{aA,aB}$:

[TABLE]

Besides we get at the lower bound using Proposition 4.7

[TABLE]

This finishes the proof of the lower bound. ∎

Let

[TABLE]

Proposition 4.10**.**

For all $t\in[T_{1},T_{=}]$,

[TABLE]

Proof.

In this time interval the newborns of genotype $aA$ are in majority produced by reproductions of a population of order one, namely $AB$ or $AA$, with the population $aA$. Since $n_{aA}$ feels competition from a macroscopic population ($AA$, $AB$ or $BB$) the $aA$ population stays of order ${\Theta}({\varepsilon})$. We make this more rigorous. To show this we consider a majorising process on $aA$ and use Proposition 4.7, and Lemma 4.9:

[TABLE]

By Proposition 4.6 and Proposition 3.4 in [8] there exists a time $t_{0}={\Theta}(1)$ such that the expression in the first bracket becomes bigger than the expression in the second bracket. Thus $n_{aA}$ decreases after $t_{0}$ and since $aA$ does not exceed ${\Theta}({\varepsilon})$ until $t_{0}$ it will stay smaller or equal to ${\Theta}({\varepsilon})$ until $T_{=}$. ∎

We show that as soon as $aB$ crosses $aA$ the $BB$ population is already bigger than or equal to the $AA$ population. First we construct a process that provide an upper bound on $aB$:

Lemma 4.11**.**

For all $t\in[T_{1},T_{2}]$ the $aB$ population is upper bounded by

[TABLE]

Proof.

First observe that the bound is fulfilled at $t=T_{1}$. Similarly to the proof of Lemma 4.9 we estimate a majorising process on $aB$ given by:

[TABLE]

By Lemma 4.2, we have to show that as soon as $aB$ reaches the upper bound it decreases faster than the bound, thus we calculate the slope of the majorising process at this value:

[TABLE]

We have to show that, at the upper bound, $\dot{n}_{aB}<C(t)\dot{n}_{aA}+\dot{C}(t)n_{aA}$. Since the 3-system $(AA,AB,BB)$ converges towards $(0,0,\bar{n}_{B})$, we know from Proposition 3.4 in [8] that $C(t)$ is a monotone increasing function and hence $\dot{C}(t)\geq 0$. Thus if we can show that $\dot{n}_{aB}<C(t)\dot{n}_{aA}$ we are done. For this we have to calculate the slope of the minorising process on $aA$ when $aB$ reaches the upper bound. This process is given by:

[TABLE]

The slope at the upper bound is:

[TABLE]

Since $C(t)>0$ this finishes the proof. ∎

Lemma 4.12**.**

We have $T_{=}<T_{2}$. Moreover,

[TABLE]

Proof.

We first show that $T_{=}<T_{2}$. Using Proposition 4.7 we construct two processes that provide an upper bound and a lower bound on $n_{aB}$. For the birth rates of these processes we use that until time $T_{=}$, $n_{aB}\leq n_{aA}={\Theta}({\varepsilon})$ (Proposition 4.10) and Lemma 4.9:

[TABLE]

We first show that $T_{=}<\infty$. We know that the 3-system $(AA,AB,BB)$ converges to $(0,0,\bar{n}_{B})$ and that $n_{aB}\leq n_{aA}={\Theta}({\varepsilon})$ (Proposition 4.10), for $t\leq T_{=}$. If we assume that $n_{aB}<n_{aA}$, then (4.132) implies that at some time $t_{0}$, where $n_{AB}+2n_{BB}$ is already macroscopic, we have

[TABLE]

Thus the time $n_{aB}$ needs to reach $n_{aA}={\Theta}({\varepsilon})$ is of order ${\Theta}(1)$. This time is shorter than $T_{aA=\delta AA}$. Indeed, suppose the contrary, then by Proposition 4.10 $n_{aA}$ does not exceed ${\Theta}({\varepsilon})$ before $T_{2}$, and thus $T^{aA=\delta AA}\geq T^{AA}_{{\Theta}({\varepsilon}/\delta)}={\Theta}\left((\delta/{\varepsilon})^{2}\right)$ which diverges with ${\varepsilon}$. A similar reasoning shows that $T_{=}<T^{aB=\delta AB}$. Hence $T_{=}<T_{2}$.

It is left to show that $n_{AA}(T_{=})\leq n_{BB}(T_{=})+{\Theta}({\Delta})$. From Lemma 4.11 we deduce that at $T_{=}$,

[TABLE]

∎

Lemma 4.13**.**

For all $t\in[T_{1},T_{2}]$ the $AB$ population is bounded by

1. (1)

$n_{AB}\geq 2\sqrt{\bar{n}_{B}n_{AA}}-2n_{AA}\left(1+{\textstyle{{\Delta}\over c\bar{n}_{B}}}\right)$, 2. (2)

$n_{AB}\leq 2\sqrt{\bar{n}_{B}n_{AA}\left(1+{\textstyle{{\Delta}\over f}}\right)}-2n_{AA}$.

Proof.

The proof works like the one of Lemma 4.9. First observe that the bound holds at $t=T_{1}$ which can be proven by using Corollary 4.5 and Proposition 4.7. Then we calculate a minorising process on $AB$:

[TABLE]

We use Proposition 4.7 and show that this minorising process would increase quicker than the lower bound if $AB$ reaches it:

[TABLE]

From Lemma 4.2 it is left to show that at the lower bound,

[TABLE]

For this we calculate a majorising process on $AA$ using Proposition 4.7 and (4.55):

[TABLE]

If we now insert the lower bound and use Proposition 4.7 we get,

[TABLE]

Thus (4.139) is fulfilled.

First, observe that the upper bound is fulfilled at $t=T_{1}$. We then have to estimate a majorising process on $AB$. For the birth rate we use (4.55) and for the death rate we use Proposition 4.7.

[TABLE]

As before we calculate the slope of this majorising process if it would reach the upper bound:

[TABLE]

By Lemma 4.2 we have to show that

[TABLE]

For this we calculate the slope of a minorising process on $AA$:

[TABLE]

At the upper bound $AA$ would start to increases:

[TABLE]

Thus we get

[TABLE]

This finishes the proof of (2). ∎

The following Proposition is a statement for the 3-system $(AA,AB,BB)$ but it remains valid until time $T_{2}$ in the 6-system $(aa,aA,AA,aB,AB,BB)$, for ${\delta}<{\Delta}$.

Proposition 4.14**.**

The maximal value $n_{AB}^{max}$ of $n_{AB}$ in $[T_{1},T_{2}]$ is bounded by

[TABLE]

Moreover, let $T_{AB}^{max}$ be the time when $n_{AB}$ takes its maximum, then $n_{AA}$ and $n_{BB}$ are bounded by

[TABLE]

Proof.

From Lemma 4.13 (1) we get that

[TABLE]

We look for the value of $AA$ where the expression on the right hand side takes on its minimum, thus we have to derivate $n_{AA}$ and set it to zero:

[TABLE]

If we insert this in $n_{AB}$ we get the lower bound:

[TABLE]

For the upper bound on $n_{AB}$ we proceed similarly. From Lemma 4.13 (2) we get

[TABLE]

Setting the derivation of the rhs to zero gives:

[TABLE]

Finally we get

[TABLE]

∎

Remark*.*

Note that by Proposition 4.7, (4.166) and (4.55) $n_{AA}=n_{BB}\pm{\Theta}({\Delta})={\textstyle{\bar{n}_{B}\over 4}}\pm{\Theta}({\Delta})$ as soon as $n_{AB}$ reaches its maximal value.

Proposition 4.15**.**

For all $t\in[T_{1},T_{2}]$,

[TABLE]

Proof.

For $t\leq T_{=}$ this follows from Proposition 4.10. For $t>T_{=}$ we show this by constructing a majorising process of $n_{aA}(t)$. For the birth rate we use (4.55) and Lemma 4.9:

[TABLE]

By Lemma 4.2, it is left to show that $\dot{n}_{aA}\leq\dot{n}_{aB}$ whenever $n_{aA}=n_{aB}$. At this upper bound we have $\dot{n}_{aA}\leq{n_{aB}}(\frac{f}{\bar{n}_{A}}(n_{AA}-n_{BB})+\eta n_{BB}+{\Theta}({\delta}))$. We now calculate a minorising process on $n_{aB}$. The birth rate can be estimate by using Proposition 4.7, (4.55) and that ${\textstyle{1/2n_{AB}+n_{BB}\over\Sigma_{5}}}\geq{\textstyle{1\over 2}}-{\Theta}({\Delta})$, for $t>T_{=}$:

[TABLE]

Thus $\dot{n}_{aB}\geq{\textstyle{f\over\Sigma_{5}}}n_{aB}(n_{BB}-n_{AA}+n_{AB})$ whenever $n_{aA}=n_{aB}$, and hence $\dot{n}_{aB}-\dot{n}_{aA}\geq{\textstyle{f\over\bar{n}_{A}}}n_{aB}(2n_{BB}-2n_{AA}+\eta n_{BB}-{\Theta}(\Delta))>0$ by Proposition 4.14. This finishes the proof. ∎

Now we show that the time $T^{aA=\delta AA}$ is finite and prove that it is smaller than or equal to $T^{aB=\delta BB}$. To estimate the order of magnitude of the time $T_{2}$ we need bounds on $n_{aA}$ which depends on $\Sigma_{aA,aB}$.

Lemma 4.16**.**

For all $t\in[T_{1},T_{2}]$ the $aA$ population is bounded by

[TABLE]

Proof.

We start with the upper bound. First observe that it holds at $t=T_{1}$. By Lemma 4.2 it is enough to show that if $n_{aA}$ would reach the upper bound it would decrease faster than the bound. Using Lemma 4.9 and (4.55) the birth rate of a majorising process on $aA$ is given by

[TABLE]

We calculate the slope of the majorising process at the upper bound using Proposition 4.7:

[TABLE]

We have to show that at the upper bound,

[TABLE]

To do this we calculate minorising processes on $n_{AB}$ and $n_{AA}$. For the birth rates we use Lemma 4.9 and (4.55):

[TABLE]

From (4.184) and (4.185) we get that

[TABLE]

By Lemma 4.8, we know that $\dot{\Sigma}_{aA,aB}\geq-{\Delta}\Sigma_{aA,aB}$. Thus, using (4.186) the right-hand side minus the left-hand side of (4.179) is bounded from below by

[TABLE]

This finishes the proof of (1).

For the lower bound we proceed similarly (using Lemma 4.2). This time we show that if $n_{aA}$ would reach the lower bound it would start to increase faster than the bound. Thus we need a minorising process on $n_{aA}$. Using Lemma 4.9 and (4.55) the birth rate of such a process is given by

[TABLE]

We calculate the slope of the minorising process at the lower bound:

[TABLE]

Thus the minorising process on $n_{aA}$ would increase when the $aA$ population would reach the lower bound. To ensure this lower bound we have to show

[TABLE]

For this we consider a majorising process on $\Sigma_{aA,aB}$. The birth rate is the same as in (4.112) and the death rate can be lower bounded by using Proposition 4.7 and Lemma 4.9:

[TABLE]

Using that $\eta<c$, the slope of this process if $n_{aA}$ reaches the lower bound is estimated by

[TABLE]

Moreover we need majorising processes on $AA$ and $AB$. Using (4.55) the birth rates of these processes can be bounded by:

[TABLE]

From (4.200) and (4.201) we get that

[TABLE]

Thus for (4.192) it is enough to show that

[TABLE]

using (4.195) and (4.191) the lhs minus the rhs of (4.203) is given by

[TABLE]

For the first inequality we use the rough estimations ${\textstyle{f-D\over f+{\Delta}}},{\textstyle{f\over f+{\Delta}}}<1$ and ${\textstyle{2n_{AA}+n_{AB}\over\bar{n}_{B}}}<2$. This concludes the proof. ∎

Proposition 4.17**.**

For all $t\in[T_{1},T_{2}]$ the process $\Sigma_{aA,aB}$ is bounded by

• (1)

$\dot{\Sigma}_{aA,aB}\leq n_{aA}\left(\eta n_{BB}-{\Delta}\frac{n_{AB}+{\Theta}({\Delta}n_{AA})}{n_{AB}+2n_{AA}}\right).$**

• (2)

$\dot{\Sigma}_{aA,aB}\geq n_{aA}(\eta n_{BB}-{\Delta}-{\Theta}({\delta})).$**

Proof.

We construct a majorising process on $\Sigma_{aA,aB}$. For the birth rate we use Lemma 4.9:

[TABLE]

To bound ${n}_{aB}$ ,we use Lemma 4.11:

[TABLE]

This time we construct a minorising process on $\Sigma_{aA,aB}$ by using Proposition 4.7 and Lemma 4.9:

[TABLE]

∎

From this Proposition we can deduce

Corollary 4.18**.**

There exists a $t^{*}\in[T_{1},T_{2}]$, such that for all $t\in[t^{*},T_{2}]$ and $\eta>\frac{4{\Delta}}{\bar{n}_{B}}=:\eta^{\star}$,

[TABLE]

Proof.

A fine calculation will show that the competition $c-\eta$ felt by an $aA$ individual from a $BB$ individual allows the sum $\Sigma_{aA,aB}$ to grow when $\eta$ is large enough, whereas it decreases when $\eta=0$. Note that we consider here the sum ${\Sigma}_{aA,aB}$ because the influence of $\eta$ cannot be seen in the rates of the $aB$ population alone. Heuristically, the growth of the $aB$ population happens due to the indirect influence (source of $a$ allele) of the less decaying $aA$ population. We prove that the minorising process of $\Sigma_{aA,aB}$ estimated in Proposition 4.17 starts to increase:

[TABLE]

As soon as $n_{BB}>{\Delta}/\eta$, the sum-process $\Sigma_{aA,aB}$ increases. With the knowledge of the behaviour of the 3-system $(AA,AB,BB)$, (see [8]), Lemma 4.12 and Proposition 4.14 we know that, for $t\geq T_{=}$, we have $n_{BB}\geq{\textstyle{\bar{n}_{B}\over 4}}-{\Theta}({\Delta})$. Hence, if we choose $\eta>{\textstyle{4{\Delta}\over\bar{n}_{B}}}$ the sum-process $\Sigma_{aA,aB}$ increases. ∎

Now we calculate the time $T^{aA=\delta AA}\land T^{aB=\delta AB}$ and we will see that $T^{aA=\delta AA}\land T^{aB=\delta AB}=T^{aA=\delta AA}$.

Theorem 4.19**.**

The time $T_{2}={\Theta}({\varepsilon}^{-1/(1+\eta\bar{n}_{B}-{\Delta})})$.

Proof.

From Proposition 4.17 (2) we have a lower bound on $\dot{\Sigma}_{aA,aB}$, and with Lemma 4.16 we can further bound $n_{aA}$ (either from above or from below depending on the sign of the prefactor):

[TABLE]

where in the last line, the estimation on $n_{BB}$ comes from Proposition 4.14, and the estimation on $n_{AB}$ comes from [27] where it is proved that $n_{AB}$ starts to decrease like $1/t$ after a time of order ${\Theta}(1)$. As ${\Sigma}_{aA,aB}(T_{1})={\Theta}({\varepsilon})$, the solution of the ODE that gives a lower bound is:

[TABLE]

By using Proposition 4.17 (1), we get the same kind of solution as an upper bound on ${\Sigma}_{aA,aB}$ (note on the last step we can upper bound $n_{BB}$ by $\bar{n}_{B}$):

[TABLE]

Using (4.216) and the lower bound in Lemma 4.16 (together with the trivial estimation $n_{AA}\geq 0$) we get a minorising process on $aA$:

[TABLE]

The corresponding majorising process has an $\bar{n}_{B}$ instead of $\bar{n}_{B}/4$. By solving $n_{aA}=\delta n_{AA}={\Theta}(n_{AB}^{2})$ we get the order of magnitude of $T_{aA=\delta AA}$:

[TABLE]

Note that $1+\eta\bar{n}_{B}-{\Delta}>0$ for ${\Delta}$ small enough, and thus $T_{aA=\delta AA}$ diverges with ${\varepsilon}$ and the order calculations above are justified.

It is left to ensure that $aB$ does not exceed $\delta n_{AB}$ in this time. It follows from Lemma 4.16 that during the time interval $[T_{1},T_{2}]$, we have ${\Sigma}_{aA,aB}={\Theta}(n_{aB})$. Thus, solving $n_{aB}=\delta n_{AB}$ amounts to solving ${\Theta}({\Sigma}_{aA,aB})={\Theta}(1)/({\Theta}(1)+{\Theta}(1)t)$ which gives the same order of magnitude as for $T_{aA=\delta AA}$. Thus the two times are of the same order.

Note that for $\eta=0$, ${\Sigma}_{aA,aB}(T_{2})={\Theta}\left({\varepsilon}^{1+{{\Delta}}/{(1-{\Delta})}}\right)$. This proves point 1 of Theorem 3.1.

∎

Proposition 4.20**.**

$T_{2}=T^{aA=\delta AA}$.

Proof.

This follows from Theorem 4.19 and Lemma 4.9. ∎

Proposition 4.21**.**

At time $t=T_{2}$ and if $f$ is taken sufficiently large (Assumption C2), $n_{aa}$ starts to grow out of itself: there exists some positive constant $c_{T_{2}}>0$ such that

[TABLE]

Proof.

We have $n_{AA}(T_{2})={\Theta}\left({\varepsilon}^{2/(1+\eta\bar{n}_{B}-\Delta)}\right)$. Thus, at the end of the second phase,

[TABLE]

the right-hand side is positive for $f$ large enough. ∎

4.4. Phase 3: Exponential growth of $aa$ until co-equilibrium with $BB$

**

Since $aa$ is growing now also out of itself it will influence the sum-process $\Sigma_{5}=n_{aA}+n_{AA}+n_{aB}+n_{AB}+n_{BB}$ and we need new lower bounds on $\Sigma_{5}$ in the following steps, the proof of this works similar to the one of Proposition 4.7 by taking into account all contributing populations. Let us compute the ODE to which $\Sigma_{5}$ is the solution:

Proposition 4.22**.**

The sum-process $\Sigma_{5}$ is the solution to

[TABLE]

Proof.

We calculate the birth- and the death-rate of $\Sigma_{5}$ under consideration of the $aa$ population:

[TABLE]

which gives the result. ∎

We introduce some notation for the order of magnitude of $n_{AA}(T_{2})$. We write $n_{AA}(T_{2})={\Theta}({\varepsilon}^{\gamma})$ with

[TABLE]

Thus the initial condition of the third Phase can be written as:

[TABLE]

Let

[TABLE]

We define the stopping time

[TABLE]

where:

[TABLE]

Note that for the following arguments, ${\gamma/10}$ can be replaced by any positive power smaller than $\gamma/2$.

We have to ensure that the $aa$ population grows until a neighbourhood of its equilibrium. From Theorem 4.19 and Proposition 4.20, at the end of the second Phase we have that $n_{aA},n_{AA},n_{aB}\leq n_{AB}<{\Delta}$. We start to bound $\Sigma_{5}$:

Lemma 4.23**.**

For all $t\in[T_{2},T_{+}]$, the sum-process $\Sigma_{5}$ is bounded from above and below by

[TABLE]

Proof.

We use Lemma 4.2 and construct minorising and majorising processes on $\Sigma_{5}$ and $n_{AB}$. First observe that the bounds are satisfied at $t=T_{2}$ by Proposition 4.7. From Proposition 4.22 we can deduce that

[TABLE]

At the lower bound we get that the sum-process would increase:

[TABLE]

and at the upper bound it would decrease:

[TABLE]

It is left to show that

[TABLE]

For this we construct a minorising process on $AB$ using that $n_{aA},n_{AA},n_{aB}\leq n_{AB}$, for the birth rate:

[TABLE]

At the lower bound we have $\dot{n}_{AB}\geq-{\textstyle{5f+12{\Delta}\over 4\bar{n}_{B}}}n_{AB}^{2}+{\Theta}\left(n_{AB}^{3}\right)$.

The lhs of (4.245) $\leq\frac{(5f+12{\Delta})(f+3{\Delta})}{4c\bar{n}_{B}^{2}}n_{AB}^{2}+{\Theta}\left(n_{AB}^{3}\right)<{\Delta}n_{AB}-{\Theta}\left(n_{AB}^{2}\right)$. At the upper bound we get $\dot{n}_{AB}\geq-{\textstyle{25f+16{\Delta}\over 4\bar{n}_{B}}}n_{AB}^{2}-{\Theta}\left(n_{AB}^{3}\right)$ and that the rhs of (4.245) is larger or equal to $-{\textstyle{4(f+{\Delta})(21f+4{\Delta})\over 4c\bar{n}_{B}^{2}}}n_{AB}^{2}-{\Theta}\left(n_{AB}^{3}\right)>-3{\Delta}n_{AB}-{\Theta}\left(n_{AB}^{2}\right)$. ∎

The following lemma ensures that the $aA$- and the $AA$ populations stay smaller than ${\Theta}\left(n_{AB}^{2}\right)$ for all $t\in[T_{2},T_{+}]$:

Lemma 4.24**.**

1. (1)

For all $t\in[T_{2},T_{+}]$, the $AA$ population is bounded from above by

[TABLE] 2. (2)

For all $t\in[T_{2},T_{+}]$ and $\eta<{\textstyle{c\over 2}}\left(1-{\textstyle{{\Theta}({\Delta})\over f-D+{\Delta}}}\right)$, the $aA$ population is bounded from above by

[TABLE]

In particular this implies $T_{+}=T^{AB}_{{\varepsilon}^{\gamma/10}}$.

Proof.

The proof uses again Lemma 4.2. First observe that by Lemma 4.13 (1) and since $n_{aA}(T_{2})={\delta}n_{AA}(T_{2})$ the bounds are satisfied at $t=T_{2}$.

1. (1)

We have to show that

[TABLE]

We construct a majorising process on $n_{AA}$:

[TABLE]

At the upper bound we thus have: $\dot{n}_{AA}\leq-{\textstyle{f\over\bar{n}_{B}}}n_{AB}^{2}+{\Theta}\left({\Delta}n_{AB}^{2}\right)$.

We now construct a minorising process on $n_{AB}$. As in (4.246) until $T_{+}$ we have a minoration of the birth rate:

[TABLE]

Thus the rhs of (4.251) is larger or equal to $-{\textstyle{25f+{\Theta}({\Delta})\over\bar{n}_{B}^{2}}}n_{AB}^{3}$. Since $n_{AB}<{\varepsilon}^{\gamma/10}$ until time $T_{+}$, we indeed have, for ${\varepsilon}$ small enough:

[TABLE] 2. (2)

Similarly to (1) we construct a majorising process on $aA$. For the birth rate we use the result of (1):

[TABLE]

At the upper bound we get: $\dot{n}_{aA}\leq n_{AB}^{2}\left({\textstyle{9f+{\Theta}({\Delta})\over 2\bar{n}_{B}}}-{\textstyle{10f(f-2\eta\bar{n}_{B})\over 2\bar{n}_{B}(D-{\Delta})}}\right)$. We will show that

[TABLE]

For this we use the minorising process on $AB$ constructed in (4.257), given by $\dot{n}_{AB}\geq-{\textstyle{25f+{\Theta}({\Delta})\over 4\bar{n}_{B}}}n_{AB}^{2}$. Thus the rhs of (4.263) is larger or equal to $-{\textstyle{105f+{\Theta}({\Delta})\over\bar{n}_{B}(D-{\Delta})}}n_{AB}^{3}.$ Thus we have to ensure that $\left({\textstyle{9f+{\Theta}({\Delta})\over 2\bar{n}_{B}}}-{\textstyle{10f(f-2\eta\bar{n}_{B})\over 2\bar{n}_{B}(D-{\Delta})}}\right)<-{\textstyle{105f+{\Theta}({\Delta})\over\bar{n}_{B}(D-{\Delta})}}{\Delta}$. This yields the condition on $\eta$:

[TABLE]

∎

Remark*.*

Observe that the condition $\eta<{\textstyle{c\over 2}}\left(1-{\textstyle{{\Theta}({\Delta})\over f-D+{\Delta}}}\right)$ in Lemma 4.24 prevents $aA$ to grow exponentially fast. Inequality (4.261) shows that for larger values for $\eta$, $n_{aA}$ would start to grow out of itself and thus the system would converge towards the 6-point equilibrium as we checked numerically. Hence, the assumption $\eta<{\textstyle{c\over 2}}\left(1-{\textstyle{{\Theta}({\Delta})\over f-D+{\Delta}}}\right)$ is essential in this phase and propagates to the following lemmata since we need therein the $n_{AB}^{2}$-dependent bound on $aA$.

Using Lemma 4.24 we can also compute a lower bound for $AA$:

Lemma 4.25**.**

For $\eta\leq{\textstyle{c\over 2}}-{\Theta}({\Delta})$ and $t\in[T_{2},T_{+}]$, the $AA$ population is bounded from below by

[TABLE]

Proof.

By Lemma 4.2 we construct a minorising process on $AA$ using Lemma 4.23 for the death rate. The bound is satisfied at $t=T_{2}$ by Lemma 4.13 (2).

[TABLE]

At the lower bound we have $\dot{n}_{AA}\geq{\textstyle{D-{\Theta}({\Delta})\over 8\bar{n}_{B}}}n_{AB}^{2}={\Theta}(n_{AB}^{2})$. It is left to show that

[TABLE]

We construct a majorising process on $AB$:

[TABLE]

Using Lemma 4.24 we get that the rhs of (4.269) is smaller or equal to ${\textstyle{f\over 4\bar{n}_{B}}}n_{AB}n_{aA}+{\Theta}\left(n_{AB}^{3}\right)\leq{\Theta}\left(n_{AB}^{3}\right)<{\Theta}(n_{AB}^{2})$ for ${\varepsilon}$ small enough since $n_{AB}\leq{\varepsilon}^{\gamma/10}$ for $t\in[T_{2},T_{+}]$. ∎

With all these lemmata we are now able to show that $n_{AB}$ stays below ${\varepsilon}^{\gamma/10}$ until $T_{3}$ when $n_{aa}$ reaches the neighbourhood of its equilibrium.

Lemma 4.26**.**

For $t\in[T_{2},T_{+}]$, it holds

[TABLE]

Proof.

As in (4.270) we construct a majorising process on $AB$ using Lemma 4.24:

[TABLE]

With the initial condition $n_{AB}(T_{2})={\Theta}\left({\varepsilon}^{\gamma/2}\right)$, we get the upper bound

[TABLE]

Thus $T_{{\varepsilon}^{\gamma/10}}^{AB}\geq{\Theta}\left({\varepsilon}^{-{\gamma}/2}\right)$. Observe that it follows that $T_{{\varepsilon}^{\gamma/5}}^{AA}={\Theta}({\varepsilon}^{-{\gamma}/2})$. ∎

To ensure the exponential growth of $aa$ we need that the $aA$ population does not decay under the order ${\Theta}({\varepsilon}^{\gamma})$.

Lemma 4.27**.**

For $\eta\leq{\textstyle{c\over 2}}-{\Theta}({\Delta})$ and for all $t\in[T_{2},T_{2}+{\Theta}({\varepsilon}^{-\gamma/2})]$, if $n_{aa}\leq n_{AB}$, then

[TABLE]

Proof.

We construct a minorising process on $aA$ in a very rough way assuming that there is no birth. Since $n_{aa}\leq n_{AB}$ the death rate can be bounded by using Lemma 4.23:

[TABLE]

By Lemma 4.24 we get:

[TABLE]

At time $T_{2}$ we know that $n_{aA}=\delta n_{AA}={\Theta}\left(\delta n_{AB}^{2}\right)={\Theta}(\delta{\varepsilon}^{\gamma})$. Thus, using (4.277) to bound $n_{AB}$, we conclude that $n_{aA}$ satisfies until $T_{+}$:

[TABLE]

We have to ensure that the slope of a minorising process on $AB$ is of the same order. Thus we construct a minorising process on $AB$ using Lemma 4.24 and that $n_{aa}\leq n_{AB}$:

[TABLE]

where for the last inequality we used that $n_{aB}\leq n_{AB}$. Using Lemma 4.23 we get:

[TABLE]

Thus the minorising process on $AB$ is given by:

[TABLE]

The minorising process (4.285) stays of order ${\Theta}({\varepsilon}^{\gamma/2})$ during a time of order ${\Theta}({\varepsilon}^{-\gamma/2})$.

The minorising process (4.281) needs time of order ${\Theta}\left({\varepsilon}^{-\gamma/2}\right)$, to reach the order ${\Theta}\left(\delta{\varepsilon}^{\gamma+\alpha}\right)$ for any $\alpha>0$. The bound (4.278) is thus ensured for a time of order ${\Theta}({\varepsilon}^{-\gamma/2})$. ∎

Now we show that $n_{aa}$ increases to a neighbourhood of its equilibrium before time $T_{+}$.

Lemma 4.28**.**

For $\eta\leq{\textstyle{c\over 2}}-{\Theta}({\Delta})$ and all $t\in[T_{2},T_{+}]$ the $aa$ population increases to a ${\varepsilon}_{0}$-neighbourhood of its equilibrium $\bar{n}_{a}$ exponentially fast, and $T_{3}<T_{+}$.

Proof.

We construct a minorising process on $aa$ and distinguish some cases. We remember that $n_{aa}(T_{2})={\Theta}({\varepsilon}^{2})$ if we do not have a better starting value. First observe that by Lemma 4.24 it holds that $d_{aa}\leq n_{aa}\left(D+{\Delta}+cn_{aa}+{\Theta}\left(n_{AB}^{2}\right)\right)$.

1. (1)

If $n_{aa}\leq n_{aA}\leq n_{AA}$ or $n_{aA}\leq n_{aa}\leq n_{AA}$

In that case the birth-rate is given by $b_{aa}\geq fn_{aa}{\textstyle{n_{aA}\over 6n_{AA}}}$. With Lemma 4.24 and 4.27 we get $b_{aa}\geq{\Theta}({\delta})fn_{aa}$ and $\dot{n}_{aa}\geq n_{aa}\left({\Theta}({\delta})f-D-{\Delta}-{\Theta}\left(n_{AB}^{2}\right)\right)$.

Hence the time until $aa$ reaches $\bar{n}_{a}-{\varepsilon}_{0}$ satisfies $T_{3}\leq{\Theta}\left(\ln({\varepsilon}^{-2})\right)={\Theta}(\ln(1/{\varepsilon}))$. 2. (2)

If $n_{aa},n_{AA}\leq n_{aA}$

In that case, by Lemma 4.24 we get $b_{aa}\geq{\textstyle{f\over 6}}n_{aa}$ and $\dot{n}_{aa}\geq n_{aa}\left({\textstyle{f\over 6}}\!-\!D\!-\!{\Delta}\!-\!{\Theta}\left(n_{AB}^{2}\right)\!\right)$.

Hence $T_{3}\leq{\Theta}\left(\ln\left({\varepsilon}^{-2}\right)\right)={\Theta}(\ln(1/{\varepsilon}))$. 3. (3)

If $n_{aA}\leq n_{AA}\leq n_{aa}\leq n_{AB}$

In that case, by Lemma 4.24, $b_{aa}\geq{\textstyle{f\over 3}}n_{aa}$ and $\dot{n}_{aa}\geq n_{aa}\!\left({\textstyle{f\over 3}}\!-\!D\!-\!{\Delta}\!-\!cn_{aa}\!-\!{\Theta}\!\left(n_{AB}^{2}\right)\!\right)$. Hence $T_{3}\leq{\Theta}\left(\ln\left({\varepsilon}^{-{\gamma}/5}\right)\right)={\Theta}(\ln(1/{\varepsilon}))$. 4. (4)

If $n_{AA}\leq n_{aA}\leq n_{aa}\leq n_{AB}$

In that case, by Lemma 4.24, $b_{aa}\geq{\textstyle{f\over 3}}n_{aa}$ and $\dot{n}_{aa}\!\geq n_{aa}\!\left({\textstyle{f\over 3}}\!-\!D\!-\!{\Delta}\!-\!cn_{aa}\!-\!{\Theta}\left(n_{AB}^{2}\right)\!\right)$. Hence $T_{3}\leq{\Theta}\left(\ln\left({\varepsilon}^{-{\gamma}/5}\right)\right)={\Theta}(\ln(1/{\varepsilon}))$. 5. (5)

If $n_{aa}>n_{AB}$

In that case, by Lemma 4.24, $b_{aa}\geq n_{aa}(f-{\Theta}({n_{AB}^{2}}/{n_{aa}}))$ and $\dot{n}_{aa}\geq n_{aa}(f-D-{\Delta}-cn_{aa}-{\Theta}({n_{AB}^{2}}/{n_{aa}}))$. Hence $T_{3}\leq{\Theta}({\varepsilon}^{{\gamma}/10}\ln({\varepsilon}^{-{\gamma}/10}))={\Theta}({\varepsilon}^{\gamma/10}\ln(1/{\varepsilon}))$.

Thus, remembering Lemma 4.26, we proved that $T_{3}\leq{\Theta}(\ln(1/{\varepsilon}))\leq{\Theta}({\varepsilon}^{-\gamma/2})\leq T_{+}.$ ∎

4.5. Phase 4: Convergence to $p_{aB}=(\bar{n}_{a},0,0,0,0,\bar{n}_{B})$

**

The Jacobian matrix of the field (2.25) at the fixed point $p_{aB}$ has the 6 eigenvalues: [math] (double), and $-(2f-D),-(f-D+\Delta),-(f-D-\Delta),-((f-D)(5f-4D)+f\Delta)/(4(f-D)+\eta\bar{n}_{B})$ which are strictly negative under Assumptions (C). Because of the zero eigenvalues, $p_{aB}$ is a non-hyperbolic equilibrium point of the system and linearisation fails to determine its stability properties. Instead, we use the result of the center manifold theory [20, 28] that asserts that the qualitative behaviour of the dynamical system in a neighbourhood of the non-hyperbolic critical point $p_{aB}$ is determined by its behaviour on the center manifold near $p_{aB}$.

Theorem 4.29** (The Local Center Manifold Theorem 2.12.1 in [28]).**

Let $f\in C^{r}(E)$, where $E$ is an open subset of ${\mathbb{R}}^{n}$ containing the origin and $r\geq 1$. Suppose that $f(0)=0$ and $Df(0)$ has $c$ eigenvalues with zero real parts and $s$ eigenvalues with negative real parts, where $c+s=n$. Then the system $\dot{z}=f(z)$ can be written in diagonal form

[TABLE]

where $z=(x,y)\in{\mathbb{R}}^{c}\times{\mathbb{R}}^{s}$, $C$ is a $c\times c$-matrix with $c$ eigenvalues having zero real parts, $P$ is a $s\times s$-matrix with $s$ eigenvalues with negative real parts, and $F(0)=G(0)=0,DF(0)=DG(0)=0.$ Furthermore, there exists $\delta>0$ and a function, $h\in C^{r}(N_{\delta}(0))$, where $N_{\delta}(0)$ is the ${\delta}$-neighbourhood of [math], that defines the local center manifold and satisfies:

[TABLE]

for $|x|<{\delta}$. The flow on the center manifold $W^{c}(0)$ is defined by the system of differential equations

[TABLE]

for all $x\in{\mathbb{R}}^{c}$ with $|x|<{\delta}$.

Theorem 4.30**.**

The non-hyperbolic critical point $p_{aB}$ is a stable fixed point and the flow on the center manifold near the critical point approaches $p_{aB}$ with speed $\frac{1}{t}$.

Proof.

We apply the Local Center Manifold Theorem 4.29. All the calculations below were done with the program Mathematica 11.0.0.0 Student Edition. We do not write the results of all the intermediary calculations as they would take a few pages and bring no more information. Instead, we describe as precisely as possible the calculations we did so that they can be checked by the reader (using a similar computer program).

By the affine transformation $(n_{aa},n_{BB})\mapsto(n_{aa}-\bar{n}_{a},n_{BB}-\bar{n}_{B})$ we get a translated system $\tilde{F}(n)$ which has a critical point at the origin. We compute the two eigenvectors corresponding to 0 eigenvalues of the Jacobian matrix of $\tilde{F}$ at the fixed point $(0,0,0,0,0,0)$, which are

[TABLE]

We perform a new change of variable to work in the basis of eigenvectors of $\tilde{F}(n)$. Let us call the new coordinates $x_{1},\ldots,x_{6}$. Let $h(x_{1},x_{2})$ be the local center manifold (still unknown). We shall look at its local shape near $(0,0)$ and expand it up to second order. Let us write

[TABLE]

We then substitute the series expansions into the center manifold equation (4.288) which gives us 12 equations for the 12 unknowns $\lambda_{3},\ldots,\mu_{6}$. Substitution of the explicit second order approximation of the center manifold equation into (4.289) yields the (order 2 approximation of the) flow on the local center manifold:

[TABLE]

where we obtain

[TABLE]

and

[TABLE]

It is left to show that the above system flows toward the origin, at least for $\eta$ smaller than a certain constant. To do that, we perform another change of variables which allows us to work in the positive quadrant. We call the new coordinates (on the center manifold) $y_{1}$ and $y_{2}$, and the new field $\hat{F}$. Observe that it is sufficient to prove that the scalar product of the field with the position is negative. We thus consider the function

[TABLE]

which is a quadratic form in $y_{1}$ and $y_{2}$. As the field $\hat{F}$ is homogeneous of degree 2 in its variables, it is enough to consider any direction given by $y_{2}=\lambda y_{1}$, and prove that $s(y_{1},\lambda y_{1})<0$ for all $\lambda>0$. As the expressions are so ugly, we work perturbatively in $f$ and consider it as large as needed. Observe that the numerator and the denominator of $s(y_{1},\lambda y_{1})$ are polynomials of degree 5 in $f$. We thus compute the coefficient in front of $f^{5}$, and obtain by a series expansion:

[TABLE]

Observe that the denominator is always negative (because by the assumption that $\eta\leq c$). We finally compute the minimal value of the ratio

[TABLE]

and obtain $r_{max}\simeq 0.593644$. Thus, asymptotically as $f\to\infty$, the field is attractive for $\eta<c\cdot r_{max}$. Thus we see that $p_{aB}$ is a stable fixed point which is approached with speed $\frac{1}{t}$ as long as $\eta<c\cdot r_{max}$.

Figure 7 shows the flow in the center manifold of the fixed point $p_{aB}$, for two values of $\eta$, one below the threshold $cr_{max}$ and one above. We see that the flow is attractive in the first case and repulsive in the second one. ∎

5. Discussion

In the rigorous results we have presented in this paper, we made some particular assumptions on the parameters of our model in order to simplify the analysis of the (already difficult) dynamical system. In this section we discuss which of these assumptions can be relaxed, based on heuristic considerations and numerical simulations.

The no-reproduction-small-competition model.

In the model considered so far, we assume that the mutation to the $B$ allele produces a new species different to the one of phenotype $a$. This is done by the no reproduction assumption between individuals of phenotype $a$ and of phenotype $B$.

These requirements are not needed to observe the recovery of the $aa$ population. In fact, what we require is that the invasion fitness of the $aa$ population into a resident $BB$ population is positive. Therefore, we can relax the no-competition assumption and add a small competition, $c_{aB}$, between $aa$ individuals and $BB$ individuals. This additional competition increases the time until $aa$ can reinvade and also modify the two- population fixed point $p_{aB}$ (see Figure 8).

Adding the factor $\eta$ accelerates the process of recovery, and, consequently, allows to increase the competition rate $c_{aB}$ (see Figure 9).

For small $\eta$ we end up in a $aa$-$BB$ equilibrium, but by accelerating (increasing $\eta$ or decreasing $c_{aB}$) the process even more, we can also end up in a 6-point equilibrium (all six population coexist) (see Figure 10).

If there is competition between individuals of phenotype $a$ and of phenotype $B$ ($c_{aB}>0$) the $aa$ and $BB$ populations have smaller equilibria as the no-competition equilibria $\bar{n}_{a}$ and $\bar{n}_{B}$, obtained when $c_{aB}=0$. Thus the competition felt from $aA$ by $aa$ and $BB$ is lower and a smaller $\eta$ is enough to observe the 6-point equilibrium.

The all-with-all model.

The assumption of no reproduction between individuals of phenotype $a$ and of phenotype $B$ is also not really necessary for the recovery of the $aa$ population. Let us discuss the all-with-all model where all phenotypes can reproduce among themselves, that is, where the reproductive compatibility is $R_{i}(j)=1$, for all $i,j\in\mathcal{G}$.

If we analyse the invasion fitness of the $aa$ population in the macroscopic $BB$ population, it is positive (and we can observe the recovery of $aa$) if the fecundity $f$ scales with ${\varepsilon}$ in such a way that $f\cdot{\varepsilon}^{2}>D+\Delta$. Indeed, in this model the whole population acts as potential partner for each individual, and the birth rate of $aa$ scales with $\Sigma_{6}$ instead of $\Sigma_{3}$.

If this requirement on $f$ is fulfilled, then numerical simulations show that most of the results carry over to this model, but with the main difference that the 2-points equilibrium is replaced by a 3-points equilibrium. The reason for this is that the reproduction between $a$ and $B$ individuals will always give birth to $aB$ individuals and thus the $aB$ population also survives (see Figure 11 (up)).

We can also add a small competition between individuals of phenotype $a$ and individuals of phenotype $B$ and still get the 3-point equilibrium (see Figure 11(up)).

As in the no-reproduction model, adding the factor $\eta$ results in accelerating the process (see Figure 11 (middle-left)).

With reasonable choices for $\eta$ and $c_{aB}$, we end up in a 6-point-equilibrium where all populations coexist (see Figure 11 (down)). Observe, the $aB$ population can be bigger than the $AB$ population, because it gets an additional birth factor from the reproduction of individuals of genotype $aa$ with individuals of phenotype $B$ which outcompetes the birth of $AB$ individuals by reproduction of individuals of genotype $AA$ and of phenotype $B$.

A. Appendix

We collect in this appendix all the important definitions of times separating phases or subphases of the process, and provide a proof summary with all the implications. Recall Definition 4.1, and Figures 4 and 6. We write $N_{\delta}(x)$ for the $\delta$-neighbourhood of $x\in{\mathbb{R}}^{6}$. We consider

[TABLE]

and

where $\gamma:=2/(1+\eta\bar{n}_{B}-\Delta)$ is the order of magnitude of $n_{AA}$ at $T_{2}$ : $n_{AA}(T_{2})={\Theta}({\varepsilon}^{\gamma})$.

We summarise below the detailed structure of the proof (we abbreviate Lemma, Proposition and Theorem by L,P, and T respectively):

1. Phase 1:

$t\in[0,T_{1}]$

Initial conditions : $(n_{aa},n_{aA},n_{AA},n_{aB},n_{AB},n_{BB})=({\Theta}({\varepsilon}^{2}),{\varepsilon},\bar{n}_{A}\pm{\Theta}({\varepsilon}),0,{\varepsilon}^{3},0)$.

With those initial conditions the following bounds hold:

[TABLE]

Here are the main steps of the proof of Propostion 4.4 and Corollary 4.5. Until the time $T_{1}$,

[TABLE] 2. Phase 2:

$t\in[T_{1},T_{2}]$

Initial conditions :

$(n_{aa},n_{aA},n_{AA},n_{aB},n_{AB},n_{BB})=({\Theta}({\varepsilon}^{2-{\Theta}({\varepsilon}_{0})}),{\Theta}({\varepsilon}^{1-{\Theta}({\varepsilon}_{0})}),\bar{n}_{A}\pm{\Theta}({\varepsilon}_{0}),{\Theta}({\varepsilon}^{1-{\Theta}({\varepsilon}_{0})}{\varepsilon}_{0}),{\varepsilon}_{0},{\Theta}({\varepsilon}_{0}^{2})$.

As the process stays uniformly bounded in time, until $T_{2}$ we have

[TABLE]

[TABLE] 3. Phase 3:

$t\in[T_{2},T_{3}]$

Initial conditions :

$(n_{aa},n_{aA},n_{AA},n_{aB},n_{AB},n_{BB})=({\Theta}({\varepsilon}^{2}),{\Theta}({\delta}{\varepsilon}^{\gamma}),{\Theta}({\varepsilon}^{\gamma}),{\Theta}(\delta{\varepsilon}^{\gamma/2}),{\Theta}({\varepsilon}^{\gamma/2}),\bar{n}_{B}-{\Theta}({\varepsilon}^{\gamma/2}))$.

[TABLE] 4. Phase 4:

$t\in[T_{3},\infty]$

Initial conditions :

$(n_{aa},n_{aA},n_{AA},n_{aB},n_{AB},n_{BB})=(\bar{n}_{a}-{\varepsilon}_{0},{\Theta}({\varepsilon}^{\gamma/5}),{\Theta}({\varepsilon}^{\gamma/5}),{\Theta}({\varepsilon}^{\gamma/10}),{\Theta}({\varepsilon}^{\gamma/10}),\bar{n}_{B}\pm{\Theta}({\varepsilon}^{\gamma/10}))$.

T4.29$\Rightarrow$ For $\eta<c\cdot 0.593644$, the fixed point $p_{aB}=(\bar{n}_{a},0,0,0,0,\bar{n}_{B})$ is stable and for ${\varepsilon},{\varepsilon}_{0}$ small enough, the system converges to it at speed $1/t$ (Theorem 4.30).