How does geographical distance translate into genetic distance?

TL;DR

This paper models how geographic structure influences genetic differentiation by representing populations as a graph and analyzing genetic distances through random walk hitting times, revealing the importance of overall metapopulation structure.

Contribution

It introduces a graph-based model linking genetic distances to metapopulation structure using random walk hitting times under specific regimes.

Findings

Genetic distance converges to a deterministic value in large populations.

Distance depends on the entire metapopulation structure, not just direct migration.

Model simplifies to a population-based model under rare mutation and migration regimes.

Abstract

Geographic structure can affect patterns of genetic differentiation and speciation rates. In this article, we investigate the dynamics of genetic distances in a geographically structured metapopulation. We model the metapopulation as a weighted directed graph, with d vertices corresponding to d subpopulations that evolve according to an individual based model. The dynamics of the genetic distances is then controlled by two types of transitions -mutation and migration events. We show that, under a rare mutation - rare migration regime, intra subpopulation diversity can be neglected and our model can be approximated by a population based model. We show that under a large population - large number of loci limit, the genetic distance between two subpopulations converges to a deterministic quantity that can asymptotically be expressed in terms of the hitting time between two random walks in the metapopulation graph. Our result shows that the genetic distance between two subpopulations does not only depend on the direct migration rates between them but on the whole metapopulation structure.

Equations

$${\delta }_t^{\gamma ,ϵ}(x,y)=\frac{1}{l^{ϵ}}\#\{k\in \{1,\dots ,l^{ϵ}\}:\text{ x and y differ at locus k }\}.$$

$$d_t^{ϵ,\gamma }(i,j)=\left(\frac{1}{n_i^{ϵ}}\sum _{x\in \{i_1,\dots ,i_{n_i^{ϵ}}\}}\underset{y\in \{j_1,\dots ,j_{n_j^{ϵ}}\}}{\min }{\delta }_t^{\gamma ,ϵ}(x,y)\right)\vee \left(\frac{1}{n_j^{ϵ}}\sum _{y\in \{j_1,\dots ,j_{n_j^{ϵ}}\}}\underset{x\in \{i_1,\dots ,i_{n_i^{ϵ}}\}}{\min }{\delta }_t^{\gamma ,ϵ}(x,y)\right).$$

$$\forall i,j\in E,\frac{1}{b^{\gamma ,ϵ}{n}_i^{ϵ}},\frac{1}{m_{ij}^{\gamma }}\gg n_j^{ϵ},l^{ϵ}\gg 1,$$

$$m_{ij}^{\gamma }$$

$$b^{\gamma ,ϵ}$$

$$n_i^{ϵ}=[N_i/ϵ]\text{ where Ni>0 remains constant as ϵ→0}$$

$$l^{ϵ}\to \infty \text{ as ϵ→0}$$

$$\forall t\geq 0,\ \ D_{t}(i,j)\ =\ 1\ -\ \int_{0}^{t}e^{-2b_{\infty}s}{\mathbb{P}}(\tau_{ij}\in ds)\ -\ e^{-2b_{\infty}t}{\mathbb{P}}\left(\tau_{ij}>t\right),\$$

$$\underset{ϵ\to 0}{\lim }\underset{\gamma \to 0}{\lim }(d_{t/(\gamma ϵ)}^{\gamma ,ϵ}(i,j),t\ge 0)=(D_t(i,j),t\ge 0)\text{ in the sense of finite dimensional distributions (f.d.d.).}$$

$$\underset{t\to \infty }{\lim }{D}_t(i,j)=1-\mathbb{E}(e^{-2b_{\infty }{\tau }_{ij}}).$$

$$d^{\epsilon}_{t}({i,j})\ =\ \frac{1}{l^{\epsilon}}\#\{\ k\in\{1,\dots,l^{\epsilon}\}\ :\ \mbox{subpopulations $i$ and $j$ differ at locus $k$ at time $t$}\ \}$$

$$\underset{\gamma \to 0}{\lim }(d_{t_1/(\gamma ϵ)}^{\gamma ,ϵ},\dots ,d_{t_k/(\gamma ϵ)}^{\gamma ,ϵ})=(d_{t_1}^{ϵ},\dots ,d_{t_k}^{ϵ})\text{ in distribution.}$$

$$\{0\leq s\leq t\ :\ \forall i\in\ E\mbox{ pop. $i$ is homogeneous at time $s$}\}.$$

$$Z_t^{\gamma ,ϵ}=Y_{(T^{\gamma ,ϵ}{)}_t^{-1}}^{\gamma ,ϵ},\text{mbox}where(T^{\gamma ,ϵ}{)}_t^{-1}=\inf \{s\ge 0:T_s^{\gamma ,ϵ}\ge t\}.$$

$${\mathcal{F}}_j^{L^{ϵ},\lambda ,\Delta t/(\gamma ϵ)}\underset{\gamma \to 0}{⟹}{\mathcal{F}}_j^{L^{ϵ},\lambda },\text{mbox}and\underset{\gamma \to 0}{\lim }{p}_{\Delta t/(\gamma ϵ),i}=\frac{1}{n_j^{ϵ}}$$

$$\underset{\text{mbox}rateofmutationintheIBM}{\underbrace{b_{\infty }{l}^{ϵ}{n}_i^{ϵ}}}\times \underset{\text{mbox}probaoffixation}{\underbrace{\frac{1}{n_i^{ϵ}}}}=b_{\infty }{l}^{ϵ},$$

$$\frac{1}{ϵ}{M}_{ij}\mathbb{P}\left(k\in {\mathcal{F}}_j^{L^{ϵ},\lambda }\right)=\frac{1}{ϵ}{M}_{ij}/n_j^{ϵ}.$$

$${\mathcal{S}}_i^k(\Pi )=({\tilde{\Pi }}_1,\dots ,{\tilde{\Pi }}_{l^{ϵ}})\text{ with }\{\begin{array}{rlr}& {\tilde{\Pi }}_j={\Pi }_j& \forall j\ne k\\ & {\tilde{\Pi }}_k={\mathcal{S}}_i({\Pi }_j)& j=k\end{array}.$$

$${\mathcal{I}}_{ij}^S(\Pi )=({\tilde{\Pi }}_1,\dots ,{\tilde{\Pi }}_{l^{ϵ}})\text{ with }\{\begin{array}{rlr}& {\tilde{\Pi }}_j={\Pi }_j& \forall k\notin S\\ & {\tilde{\Pi }}_k={\mathcal{I}}_{ij}({\Pi }_k)& \forall k\in S\end{array}.$$

$${\mathbb{G}}^{ϵ,L^{ϵ}}h(\Pi )$$

$$<\cdot ,\cdot >:{\mathcal{M}}_d\times {\mathcal{M}}_d$$

$$m,v$$

$$X$$

$${\xi }_t^{ϵ,L^{ϵ}}:=X({\Pi }^{ϵ,L^{ϵ}}(t))$$

$${\xi }_t^{ϵ}\equiv {\xi }_t^{ϵ,{\mathcal{L}}^{ϵ}}=X({\Pi }^{ϵ,{\mathcal{L}}^{ϵ}}(t))\text{where Lϵ∼U([0,1]lϵ)}$$

$$d_t^{ϵ}(i,j)=1-{\xi }_t^{ϵ}(\{\pi \in {\mathcal{P}}_d:i{\sim }_{\pi }j)\}).$$

$$G^{ϵ}g(\pi )=\sum _{i,j=1}^d{M}_{ij}\frac{1}{ϵn_j^{ϵ}}\left(g({\mathcal{I}}_{ij}(\pi ))-g(\pi )\right)+b_{\infty }\sum _{i=1}^d\left(g({\mathcal{S}}_i(\pi ))-g(\pi )\right).$$

$$G^{\epsilon}g(\pi)\to Gg(\pi):=\sum_{i,j=1}^{d}\tilde{M}_{ji}\left(\ g({\cal I}_{ij}(\pi))-g(\pi)\ \right)\ +\ b_{\infty}\sum_{i=1}^{d}\left(\ g({\cal S}_{i}(\pi))-g(\pi)\ \right)\ \ \mbox{as $\epsilon\to 0$.}$$

$$\forall m\in {\mathcal{M}}_d,{}^{t}Gm=\sum _{i,j=1}^d{\tilde{M}}_{ji}({\mathcal{I}}_{ij}\ast m-m)+b_{\infty }\sum _{i=1}^d\left({\mathcal{S}}_i\ast m-m\right).$$

$${\xi }_0^{ϵ}\underset{ϵ\to 0}{⟶}{P}^0.$$

$$({\xi }_t^{ϵ};t\ge 0)\underset{\begin{array}{l}ϵ\to 0\end{array}}{⟹}(P_t;t\ge 0)\text{mbox}intheweaktopology,$$

$$\frac{d}{ds}{P}_s={}^{t}G{P}_s$$

Full text

††journal: Stochastic Processes and their Applications

1 Introduction

1.1 Genetic distances in structured populations. Speciation

In most species, the geographical range is much larger than the typical dispersal distance of its individuals. A species is usually structured into several local subpopulations with limited genetic contact. Because migration only connects neighbouring populations, more often than not, populations can only exchange genes indirectly, by reproducing with one or several intermediary populations. As a consequence, the geographical structure tends to buffer the homogenising effect of migration, and as such, it is considered to be one of the main drivers for the persistence of genetic variability within species (see Malécot (1948) or Karlin (1982)).

The aim of this article is to present some analytical results on the genetic composition of a species emerging from a given geographical structure. The main motivation behind this work is to study speciation. When two populations accumulate enough genetic differences, they may become reproductively isolated, and therefore considered as different species. As the geographic structure of a species is one of the main drivers for the genetic differentiation between subpopulations, this work should shed light on which are the geographic conditions under which new species can emerge.

Several authors have studied parapatric speciation, i.e. speciation in the presence of gene flow between subpopulations, for example Gavrilets et al. (1998, 2000a, 2000b) and Yamaguchi and Iwasa (2013, 2015). In their models, some loci on the chromosome are responsible for reproductive isolation. These loci may be involved in incompatibilities at any level of biological organisation (molecular, physiological, behavioural etc) and either prevent mating (pre-zygotic incompatibilities) or prevent the development of hybrids (post-zygotic incompatibilities). The number of segregating loci increases through the accumulation of mutations, and decreases after each migration event (creating the opportunity for some gene exchange between the migrants and the host population). When the number of segregating loci between two individuals reaches a certain threshold, they become reproductively incompatible. For example, Yamaguchi and Iwasa (2013, 2015) studied the case of a metapopulation containing two homogeneous subpopulations. The authors studied how the genetic distance, defined as the number of loci differing between the two subpopulations, evolves through time, using a continuous-time model. When considering metapopulations with more than two subpopulations, this kind of dynamics may translate into complex patterns of speciation. One particularly intriguing example is the case of ring species (Noest, 1997; Gavrilets et al., 1998), where two neighbouring subpopulations are too different to be able to reproduce with one another but can exchange genes indirectly, by reproducing with a series of intermediate subpopulations that form a geographic ‘ring’. How these patterns emerge and are maintained is still poorly understood, and we hope that our analytical result might shed some new light on the subject.

1.2 Population divergence and fitness landscapes

To study speciation by accumulation of genetic differences, we model the evolution of some loci on the chromosome, that are potentially involved in reproductive incompatibilities. To visualise these evolutionary dynamics, Wright (1932) suggested the metaphor of adaptive landscapes. Adaptive landscapes represent individual fitness as a function defined on the genotype space, which is a multi-dimensional space representing all possible genotypes. Wright emphasised the idea of ‘rugged’ adaptive landscapes, with peaks of fitness representing species and valleys representing unfit hybrids. Speciation, seen as a population moving from one peak to another, implies a temporary reduction in fitness, which is not very likely to occur in large populations, where genetic drift is not important enough to counterbalance the effect of selection (see Gavrilets (1997) for a more detailed discussion). However, Gavrilets (1997) suggested the idea of ‘holey’ adaptive landscapes, where local fitness maxima can be partitioned into connected sets (called evolutionary ridges). Speciation is therefore seen as a population diffusing across a ridge, by neutral mutation steps, until it stands at the other side of a hole. Theoretical models, such as Gavrilets and Gravner (1997), have shown, using percolation theory, that in high-dimensional genotype spaces, fit genotypes are typically connected by evolutionary ridges.

Our model (see Section 1.3) is built in this framework. In fact we will assume that, in large populations, deleterious mutations are washed away by selection at the micro-evolutionary timescale and describe the evolutionary dynamics for our set of incompatibility controlling loci as neutral (any genotype on the evolutionary ridge can be accessed by single mutation neutral steps). This is the idea behind the description of our model in Section 1.3.

Further, we consider that the evolutionary dynamics along the ridge are slow (as random mutations are very likely to be deleterious, mutations along the evolutionary ridge are assumed to be rarer than in the typical population genetics framework), which is why we study our model in a low mutation - low migration regime (see Section 1.4 for more details). This assumption is commonly made when studying speciation, for example in Gavrilets et al. (2000a) or Yamaguchi and Iwasa (2013).

1.3 An individual based model (IBM)

We model the metapopulation as a weighted directed graph with $d$ vertices, corresponding to the different subpopulations. Each directed edge $(i,j)$ is equipped with a migration rate in each direction. (In particular, if two subpopulations are not connected, we assume that the migration rates are equal to 0.) We assume the existence of two scaling parameters, $\gamma$ and $\epsilon$, that will converge to [math] successively (first $\gamma\to 0$ and then $\epsilon\to 0$, see Section 1.4 for more details).

Each subpopulation consists of $n_{i}^{\epsilon}$ individuals, $i\in E:=\{1,\dots,d\}$. Each individual carries a single chromosome of length 1, which contains $l^{\epsilon}$ loci of interest (that are involved in reproductive incompatibilities). We assume that the vector of positions for those loci – denoted by ${\cal L}^{\epsilon}=\{x_{1},\dots,x_{l^{\epsilon}}\}$ – is obtained by throwing $l^{\epsilon}$ uniform random variables on $[0,1]$. The positions are chosen randomly at time 0, but are the same for all individuals and do not change through time. Recall that the upper indices (such as in $l^{\epsilon}$ and ${\cal L}^{\epsilon}$) are indices and not exponents).

Conditioned on ${\cal L}^{\epsilon}$, each subpopulation then evolves according to an haploid neutral Moran model with recombination.

• •

Each individual $x$ reproduces at constant rate $1$ and chooses a random partner $y$ ($y\neq x$). Upon reproduction, their offspring replaces a randomly chosen individual in the population.

• •

The new individual inherits a chromosome which is a mixture of the parental chromosomes. Both parental chromosomes are cut into fragments in the following way: we assume a Poisson Point Process of intensity $\lambda$ on $[0,1]$. Two loci belong to the same fragment iff there is no atom of the Poisson Point Process between them. For each fragment, the offspring inherits the fragment of one of the two parents chosen randomly.

To our Moran model we add two other types of events:

• •

Mutation occurs at rate $b^{\gamma,\epsilon}$ per individual, per locus according to an infinite allele model.

• •

Migration from subpopulation $i$ to subpopulation $j$ occurs at rate $m^{\gamma}_{ij}$. At each migration event, one individual migrates from subpopulation $i$ to $j$, and replaces one individual chosen uniformly at random in the resident population. (We set $\forall i\in E,\ m^{\gamma}_{ii}=0$.)

We define the genetic distance between two individuals $x$ and $y$ at time $t$ as:

[TABLE]

Consider two subpopulations $i$ and $j$ and let $\{i_{1},\dots,i_{n^{\epsilon}_{i}}\}$ be the individuals in population $i$ and $\{j_{1},\dots,j_{n^{\epsilon}_{j}}\}$ the individuals in population $j$. The genetic distance between subpopulations $i$ and $j$ at time $t$ is defined as follows:

[TABLE]

This corresponds to the so-called modified Haussdorff distance between subpopulations, as introduced by Dubuisson and Jain (1994). (This distance has the advantage of averaging over the individuals in each subpopulation, so introducing a single mutant or migrant would produce a smooth variation in the genetic distances.)

1.4 Slow mutation–migration and large population–dense site regime.

In this section, we start by describing in more details the slow mutation–migration regime alluded to in Sections 1.2 and 1.3.

It is well known that in the absence of mutation and migration, the neutral Moran model describing the dynamics at the local level reaches fixation in finite time i.e. after a finite amount of time the population becomes homogeneous. The average time to fixation for a single locus is of the order of the size of the subpopulation (Kimura and Ohta, 1969; Kimura, 1968) (In our multi-locus model, it will also depend on the number of loci and on the recombination rate $\lambda$.) Heuristically, if we assume a low mutation - low migration regime, i.e. that

[TABLE]

the average time between two migration events ($1/{m_{ij}^{\gamma}}$), and the average time between two successive mutations at a given locus ($1/({b^{\gamma,\epsilon}n_{j}^{\epsilon}})$) are much larger than the average time to fixation. This ensures that the fixation process is fast compared to the time-scale of mutation and migration, and, as a result, when looking at a randomly chosen locus, subpopulations are homogeneous except for short periods of time right after a migration event or a mutation event. This suggests that if we accelerate time properly, we can neglect intra-subpopulation diversity and approximate our model by a population based model.

Inspired by these heuristics, we are going to take a low mutation - low migration regime, by making the mutation and migration rates depend on the scaling parameter $\gamma$ in the following way:

[TABLE]

Recall that we take a slow mutation-migration regime but the recombination rate is constant. This is consistent with the fact that in most species, mutation rates are very low (they vary from $10^{-6}$ to $10^{-8}$ per base per generation) compared to the recombination rates (for example, for the human genome there are approximatively 66 crossovers per generation).

In a second step, we will make an additional approximation: we will consider a large population - dense site limit. In fact, our second scaling parameter $\epsilon$, corresponds to the inverse of a typical subpopulation size. The parameters of the model depend on $\epsilon$ in the following way (corresponding to the second inequality in (2):

[TABLE]

In this article, we are going to take the limits successively: first $\gamma\to 0$ and then $\epsilon\to 0$, in order to be consistent with the informal inequality (2). We are now ready to state the main result of this paper.

Theorem 1.1**.**

For each pair of subpopulations $i,j\in E$, let $S^{i}$ and $S^{j}$ be two independent random walks on $E$ starting respectively from $i$ and $j$ and whose transition rate from $k$ to $p$ is equal to $\tilde{M}_{kp}:=M_{pk}/N_{k}$. Finally, define $D_{t}(i,j)$ as

[TABLE]

where $\tau_{ij}=\inf\{t\geq 0\ :\ S^{i}(t)=S^{j}(t)\}.$

If at time [math] the metapopulation is homogeneous (i.e. all the individuals in all subpopulations share the same genotype) then

[TABLE]

In particular,

[TABLE]

This result can be seen as a law of large numbers over the chromosome. Although the loci are linked and they do not fix independently (and the recombination rate is constant), when considering a large number of them, they become decorrelated, regardless of the value of $\lambda$. (Note that the limiting process does not depend on $\lambda$.) The model behaves as if infinitely many loci evolved independently according to a Moran model with inhomogeneous reproduction rates (see Remark 2.4). The expression of the genetic distances has then a natural genealogical interpretation. $S^{i}$ and $S^{j}$ can be interpreted as the ancestral lineages starting from $i$ and $j$, and our genetic distance is related to the probability that those lines meet before experiencing a mutation (or in other words, that $i$ and $j$ are Identical By Descent (IBD)).

Remark 1.2**.**

In Theorem 1.1, we considered a rather restrictive initial condition. In Section 5, we give a stronger version of this theorem, which works for a larger range of initial conditions, but that requires to introduce several cumbersome notations.

1.5 Consequences of our result

One interesting consequence of our result is that the genetic distance does not coincide with the classical graph distance, but instead it depends on all possible paths between $i$ and $j$ in the graph, and all the migration rates (and not only the shortest path and the direct migration rates $M_{ij}$ and $M_{ji}$), i.e., it does not only depend on the direct gene flow between $i$ and $j$ but on the whole metapopulation structure. In particular, this suggests that adding new subpopulations to the graph (which would correspond to colonisation of new demes), removing any edge (which could correspond to the emergence of a geographical or reproductive barrier between two subpopulations), or changing any migration rate (which could correspond to modifying the habitat structure, for example) can potentially modify the whole genetic structure of the population.

One striking illustration of the previous discussion is presented in Section 6, where we consider an example where a geographic bottleneck is dramatically amplified in our new metric. See Figure 1 and Theorem 6.1 for a more precise statement. If we consider, as in Yamaguchi and Iwasa (2013), that two populations are different species if their genetic distance reaches a certain threshold, that will mean that this metapopulation structure promotes the emergence of two different species, each one corresponding to the population in one subgraph. Very often, parapatric speciation is believed to occur only in the presence of reduced gene flow. Our example shows that in the presence of a geographic bottleneck, genetic differentiation is manly driven by the geographical structure of the population, i.e., even if the gene flow between two neighbouring subpopulations is approximatively identical in the graph, the genetic distance is dramatically amplified at the bottleneck (see Figure 1).

We note that using the hitting time of random walks as a metric on graphs is not new, and has been a popular tool in graph analysis (see Doyle and Snell (1984) and Klein and Randic (1993)). For example, the commute distance, which is the time it takes a random walk to travel from vertex $i$ to $j$ and back, is commonly used in many fields such as machine learning (Von Luxburg et al., 2014), clustering (Yen et al., 2005), social network analysis (Liben-Nowell and Kleinberg, 2003), image processing (Qiu and Hancock, 2005) or drug design (Ivanciuc, 2000; Roy, 2004). In our case the genetic distance is given by the Laplace transform of the hitting time between two random walks, which was already suggested as a metric on graphs by Hashimoto et al. (2015). In that paper the authors claimed that this metrics preserves the cluster structure of the graph. In the example alluded to above (Section 6), we found that our metric reinforces the cluster structure of the metapopulation graph. In other words, a clustered geographic structure tends to increase genetic differentiation.

1.6 Discussion and open problems

As already mentioned above, the main result is obtained by: (i) proving that, in a low mutation - low migration regime (i.e., when $\gamma\to 0$), subpopulations are monomorphic most of the time and our individual based model converges to a population based model, (ii) showing that, under a large population - dense site limit (i.e. taking $\epsilon\to 0$), the genetic distances between subpopulations (for the population based model) converge to a deterministic process (defined in Theorem 1.1). Taking these two limits successively gives no clue on how the parameters should be compared to ensure the approximation to be correct. It would be interesting to take the limits simultaneously but it is technically challenging (for example we would need to characterise the time to fixation for $l$ loci that do not fix independently, which is not easy).

As discussed in the previous paragraph, we can only show our results under some rather drastic constraints: subpopulations are asymptotically monomorphic. More generally, we believe that Theorem 1.1 should hold under relaxed assumptions, namely when the intra-subpopulation genetic diversity is low compared to the inter-subpopulation diversity (see Figure 2 for an example, where $\gamma=2e^{-6}$ and $\epsilon=5e^{-3}$ ). Technically, this would correspond to the condition that at a typical locus (i.e, a locus chosen uniformly at random) each subpopulations is monomorphic at that site with high probability (which is in essence (2)). Of course, proving such a result would be much more challenging, but would presumably correspond to a more realistic situation.

1.7 Outline

In Section 2, we show that in the rare mutation-rare migration regime (i.e. when $\gamma\to 0$ whereas $\epsilon$ remains constant), the individual based model (IBM) described above converges to a population based model (PBM) (see Theorem 3). This PBM is a generalization of the model proposed by Yamaguchi and Iwasa (2013, 2015) in three ways. First, it is an extension of their model from two to an arbitrary number of subpopulations, which is not trivial from a mathematical point of view. Second, in Yamaguchi and Iwasa (2013, 2015), the authors only assumed that the migrant alleles are fixed independently at every locus. To make the model more realistic, we took into account genetic linkage, which introduces a non-trivial spatial correlation between loci (along the chromosome). Finally, we suppose that the loci are distributed randomly along the chromosome (and not in a regular fashion). Section 2 is interesting on its own since it provides a theoretical justification of the model proposed by Yamaguchi and Iwasa (2013, 2015).

In Section 3 and 4, we study the PBM in the large population - dense site limit (i.e. when $\epsilon\to 0$). We properly introduce the main tool used to study the population based model – the genetic partition probability measure – and show an ergodic theorem related to this process (see Theorem 3.1).

Finally, in Section 5 we prove our main result (Theorem 5.1 which is an extension of Theorem 1.1) by combining the results of the previous sections.

Section 6 proves the result related to the geographical bottleneck alluded to in Section 1.5 (see Proposition 6).

2 Approximation by a population based model

We now describe a population based model (PBM) that can be seen as the limit of the IBM presented above, when $\gamma$ goes to [math] (whereas $\epsilon$ remains fixed) and time is rescaled by $1/(\gamma\epsilon)$. Consider a metapopulation where the individuals are characterised by a finite set of loci, whose positions are distributed as $l^{\epsilon}$ uniform random variables on $[0,1]$, and let ${\cal L}^{\epsilon}$ be the vector of the positions of the loci (as described in Section 1.3 for the IBM). We now describe the dynamics of the model, conditional on ${\cal L}^{\epsilon}=L^{\epsilon}$, with $L^{\epsilon}\in[0,1]^{l^{\epsilon}}$.

Before going into the description of our model, we start with a definition. It is well known that the Moran model reaches fixation in finite time, i.e., after a (random) finite time, every individual in the population carries the same genetic material, and from that time on, the system remains trapped in this configuration (see Kimura and Ohta (1969), Kimura (1968)).

Definition 2.1**.**

Consider a single population of size $n_{j}^{\epsilon}$ formed by a mutant individual (the migrant) and $n^{\epsilon}_{j}-1$ residents, that evolves according to a Moran model with recombination at rate $\lambda$ (as described in Section 1.3). We define ${\cal F}^{L^{\epsilon},\lambda}_{j}$ as the (random) set of loci carrying the mutant type when the population becomes homogeneous. (Note that ${\cal F}^{L^{\epsilon},\lambda}_{j}$ is potentially empty.)

We are now ready to describe our PBM. We represent each subpopulation as a single chromosome, which is itself represented by the set of loci $L^{\epsilon}$. The dynamics of the population can then be described as follows.

• •

For every $i\in E$: fix a new mutation in population $i$ at rate $b_{\infty}l^{\epsilon}$, the locus being chosen uniformly at random along the chromosome.

• •

For every $i,j\in E$ and every $S\subseteq\{1,\dots,l^{\epsilon}\}$: at every locus in $S$, fix simultaneously the alleles from population $i$ in population $j$ at rate $\frac{1}{\epsilon}M_{ij}{\mathbb{P}}\left({\cal F}^{L^{\epsilon},\lambda}_{j}=S\right)$.

In the PBM (parametrised by $\epsilon$), we define the genetic distance between subpopulations $i$ and $j$ at time $t$ as follows:

[TABLE]

as opposed to $d^{\gamma,\epsilon}$ which will refer to the genetic distances in the IBM as described in Section 1.3 (parametrised by $\gamma$ and $\epsilon$). We note that the definition of the genetic distance in the PBM is consistent with the one in the IBM (see (1)) in the sense that if the subpopulations are homogeneous in the IBM, (1) is equal to the RHS of the previous equation. We are now ready to state the main result of this section.

Theorem 2.2**.**

Assume that, at time [math], the subpopulations in the IBM are homogeneous and that $\forall i,j\in E$, $d^{\gamma,\epsilon}_{0}(i,j)=d^{\epsilon}_{0}(i,j)$. Then, for every $k\in\mathbb{N}$, $\forall\ 0\leq t_{1}<\dots<t_{k}$,

[TABLE]

Proof.

Recall that the loci are distributed randomly along the chromosome. In the proof, we assume that the vector of the positions of the loci ${\mathcal{L}}^{\epsilon}$ is fixed and equals to $L^{\epsilon}\in[0,1]^{l^{\epsilon}}$ (and is the same in the IBM and in the PBM). We also consider that IBM and the PBM start from the same deterministic initial condition. The unconditional extension of the proof can be easily deduced from there.

We define a coupling between the IBM and a new PBM that is close (in distribution) to the PBM defined at the beginning of this section. The idea behind the coupling is that, when time is accelerated by $1/(\gamma\epsilon)$, and $\gamma$ is small, in the IBM, the time to fixation after a mutation or migration event is short enough so that the population has become homogeneous before the next mutation or migration event takes place. Then, we can decompose the trajectories of the IBM into periods where the population is homogeneous (and waits for the next mutation or migration event to take place) and homogenization phases (where the dynamics of the population is described by a Moran model). See Figure 3 for an illustration of this concept.

More formally, let us consider $(Y_{t}^{\gamma,\epsilon};t\geq 0)$ the process recording the genetic composition in the IBM (i.e. a matrix containing the sequences of the chromosomes of all the individuals in the metapopulation) after rescaling time by $\gamma\epsilon$ so that

For $i\in E$, mutation events on the subpopulation $i$ occurs according to a Poisson Point Process (PPP) with intensity measure $b_{\infty}l^{\epsilon}n_{i}^{\epsilon}dt$. 2. 2.

For $i,j\in E$, migration events from $i$ to $j$ can be described in terms of a PPP with intensity measure $M_{ij}/\epsilon dt$.

Define ${\mathcal{E}}^{\gamma,\epsilon}$ the event that every time a subpopulation is affected by a mutation or a migration event on the interval $[0,T]$, the subpopulation is genetically homogeneous when the event occurs (as in Figure 3). In other words, there is no overlap between mutation and migration homogenization periods. The time to fixation in our (multi-locus) Moran model only depends on the number of individuals and the number of loci, so in our model it only depends on $\epsilon$ (but not on $\gamma$). As a consequence, ${\mathbb{P}}\left({\mathcal{E}}^{\gamma,\epsilon}\right)\to 1$ as $\gamma\to 0$.

Next, let us consider $T_{t}^{\gamma,\epsilon}$ as the Lebesgue measure of

[TABLE]

In words, $(T_{t}^{\gamma,\epsilon};t\geq 0)$ is the random clock which is obtained by skipping the homogenization period after a migration or mutation event (i.e. by skipping the red intervals in Figure 3). By arguing as in the previous paragraph, as $\gamma\to 0$, it is not hard to see that $(T^{\gamma,\epsilon}_{t};t\geq 0)$ converges to the identity in the Skorokhod topology on $[0,T]$ for every $T\geq 0$.

Let us now consider

[TABLE]

By construction, this process defines a PBM in the sense that at every time $t$, any subpopulation is composed by genetically homogeneous individuals. Further, since $(T^{\gamma,\epsilon}_{t};t\geq 0)$ converges to the identity and mutation and migration events occur at Poisson times, the finite dimensional distributions of $Z^{\gamma,\epsilon}$ are good approximations of the ones for the IBM.

Let us now show that $Z^{\gamma,\epsilon}$ (constructed from the IBM) is close in distribution to the PBM defined at the begining of this section. Conditioned on the event ${\cal E}^{\gamma,\epsilon}$ (whose probability goes to $1$) and on the PPP’s described in 1 and 2 above, the PBM $Z^{\gamma,\epsilon}$ can be described as follows. Define $p_{\Delta t,i}$ to be the probability for a mutant allele (at a given locus of a given individual) to fix in a population of size $n_{i}^{\epsilon}$, conditioned on the homogenization time to be smaller than $\Delta t$. Then the distribution of the conditioned PBM $Z^{\gamma,\epsilon}_{t}$ can be generated as follows.

1. (a)

At every mutation time $t$ in subpopulation $i\in E$, choose a locus $k$ uniformly at random and fix the mutation instantaneously with probability $p_{\frac{\Delta t}{\gamma\epsilon},i}$, where $\Delta t$ is the time between $t$ and the next mutation or migration event (in our new time scale). We note that if the mutation does not fix, then $Z^{\gamma,\epsilon}$ is not affected by the mutation event, and as a consequence “effective mutation” events in $Z^{\gamma,\epsilon}$ are obtained from the mutation events in the IBM after thinning each time with their respective probability $p_{\frac{\Delta t}{\gamma\epsilon},i}.$ 2. (b)

At every migration event $t$ on subpopulation $j$, fix a random set $S$ where $S$ is chosen according to ${\mathcal{F}}_{j}^{L^{\epsilon},\lambda,\Delta t/(\gamma\epsilon)}$, where $\Delta t$ is defined as in the previous point, and where ${\mathcal{F}}_{j}^{L^{\epsilon},\lambda,s}$ is the random variable ${\mathcal{F}}_{j}^{L^{\epsilon},\lambda}$ conditioned on the homogenization to occur in a time smaller smaller than $s$.

Since fixation (of one of the alleles) occurs in finite time almost surely, and the distribution of the homogenization time only depends on $\epsilon$, we have

[TABLE]

where the RHS of the second limit is the probability of fixation of a mutant allele in the absence of conditioning.

Putting all the previous observations together, one can easily show that the genetic distance in $Z^{\gamma,\epsilon}$ converges (in the finite dimensional distributions sense) to the ones of the PBM. In particular, we recover the mutation rate on subpopulation $i$ in the PBM

[TABLE]

which corresponds to the limiting “effective mutation rate” in the PBM $Z^{\gamma,\epsilon}$ (see (a) above) as $\gamma\to 0$. This completes the proof of Theorem 3. ∎

We note that Theorem 3 could be extended to the case where the subpopulations are not homogeneous in the IBM at time $t=0$. Indeed, arguing as in the proof of Proposition 3, if we start with some non-homogeneous initial condition, then each island becomes homogeneous before experiencing any mutation or mutation event with very high probability. In order to get an efficient coupling between the IBM and the PBM, we simply choose the initial condition of the PBM as the (random) state of the IBM after this initial homogenization period.

Remark 2.3**.**

Choose a locus $k\in\{1,\dots,l^{\epsilon}\}$. We let the reader convince herself that in the PBM the genetic composition at a given locus $k$ follows the following Moran-type dynamics:

1. (mutation)

“Individual” $i$ takes on a new type (or allele) at rate $b_{\infty}$. 2. (reproduction)

“Individual” $j$ inherits its type from “individual” $i$ at rate $(1/{\epsilon})M_{ij}{\mathbb{P}}\left(k\in{\cal F}^{L^{\epsilon},\lambda}_{j}\right).$ Further, in a neutral one-locus Moran model, the probability of fixation of a single allele in a resident population of size $n_{j}^{\epsilon}$ is equal to its initial frequency, which in our case is $1/n_{j}^{\epsilon}$. Thus

[TABLE]

This dynamics is not dependent on the position of the locus under consideration.

Remark 2.4**.**

Our model can be seen as a multi-locus Moran model with inhomogeneous reproduction rates. The main difficulty in analysing this model stems from the fact that there exists a non trivial correlation between loci. This correlation is induced by the fact that fixation of migrant alleles can occur simultaneously at several loci during a given migration event. In turn, the set of fixed alleles during a given migration event is determined by the local Moran dynamics described in the Introduction.

3 Large population - dense site limit

In this Section, we study the PBM described in the Section 2, in the large population - dene site limit. In particular, we study the dynamics of the genetic distances and we state Theorem 3.1, that together with Theorem 3 implies the main result of this article, namely Theorem 5.1, that is a stronger version of Theorem 1.1 (see Section 5).

3.1 The genetic partition measure

The main difficulty in dealing with the genetic distance is that it lacks the Markov property, and as a consequence, it is not directly amenable to analysis. In fact, when $d>2$, a migration event from $i$ to $j$ can potentially have an effect on the genetic distance between $j$ and another subpopulation $k$ (see Figure 4 for an example).

To circumvent this difficulty, we now introduce an auxiliary process – the genetic partition probability measure – from which one can easily recover the genetic distances (see (5) below), and whose asymptotical dynamics is explicitly characterised in Theorem 3.1 below.

Let ${\cal P}_{d}$ be the set of partitions of $\{1,\dots,d\}$. Fix $\pi\in{\cal P}_{d}$ and $i,j\in E$. Define $\mathcal{S}_{i}(\pi)$ as the element of ${\cal P}_{d}$ obtained from $\pi$ by making $i$ a singleton (e.g., ${\cal S}_{2}(\{1,2,3\})=\{1,3\}\{2\}$). Define ${\cal I}_{ij}(\pi)$ as the element of ${\cal P}_{d}$ obtained from $\pi$ by displacing $j$ into the block containing $i$ (e.g., ${\cal I}_{2,3}(\{1,3\}\{2\})=\{1\}\{2,3\}$).

At every locus $k\in\{1,\dots,l^{\epsilon}\}$ (ordered in increasing order along the chromosome) and every time $t$, the allele composition of the metapopulation induces a partition on $E$. More precisely, at locus $k$, two subpopulations are in the same block of the partition at time $t$ iff they share the same allele at locus $k$.

In the following, fix $L^{\epsilon}\in[0,1]^{l^{\epsilon}}$, the vector containing the positions of the loci. In the PBM parametrised by $\epsilon$, we condition on the loci being located at $L^{\epsilon}$ and for every $k\in\{1,\dots,l^{\epsilon}\}$ we let $\Pi^{\epsilon,L^{\epsilon}}_{k}(t)$ be the partition induced at locus $k$ (see Figure 4). The vector $\Pi^{\epsilon,L^{\epsilon}}(t)=\left(\Pi^{\epsilon,L^{\epsilon}}_{k}(t);\ k\in\{1,\dots,l^{\epsilon}\}\right)$ describes the genetic composition of the population at time $t$. According to the description of our dynamics, $\Pi^{\epsilon,L^{\epsilon}}(t)$ is a Markov chain with the following transition rates:

• •

(mutation) For every $i\in E$, $k\in\{1,\dots,l^{\epsilon}\}$, define ${\cal S}_{i}^{k}$ to be the operator on $({\mathcal{P}}_{d})^{l^{\epsilon}}$ such that $\forall\ \Pi\ \in({\mathcal{P}}_{d})^{l^{\epsilon}}$

[TABLE]

The transition rate of the process from state $\Pi$ to ${\cal S}_{i}^{k}(\Pi)$ is given by $b_{\infty}$.

• •

(migration from $i$ to $j$) For every $i,j\in E$, $S\subset\{1,\dots,l^{\epsilon}\}$, define ${\cal I}_{ij}^{S}$ the operator on $({\mathcal{P}}_{d})^{l^{\epsilon}}$ such that $\forall\ \Pi\ \in({\mathcal{P}}_{d})^{l^{\epsilon}}$

[TABLE]

The transition rate of the process from $\Pi$ to ${\cal I}_{ij}^{S}(\Pi)$ is given by $\frac{M_{ij}}{\epsilon}{\mathbb{P}}\left({\cal F}^{L^{\epsilon},\lambda}_{j}=S\right)$.

To summarise, for any function $h:{\cal P}_{d}^{l^{\epsilon}}\to\mathbb{R}$, the generator of $\Pi^{\epsilon,L^{\epsilon}}$ can be written as

[TABLE]

3.2 Some notation

Let ${\cal M}_{d}$ denote the space of signed finite measures on ${\cal P}_{d}$. Since ${\cal P}_{d}$ is finite, we can identify elements of ${\cal M}_{d}$ as vectors of $\mathbb{R}^{Bell_{d}}$, where $Bell_{d}$ is the Bell number, which counts the number of elements in ${\cal P}_{d}$ (the number of partitions of $d$ elements). In particular, if $\pi$ is a partition of $E$, and $\mu\in{\cal M}_{d}$, then $\mu(\pi)$ will correspond to the measure of the singleton $\{\pi\}$, or equivalently, to the “$\pi^{th}$ coordinate” of the vector $\mu$. We define the inner product $<\cdot,\cdot>$ as

[TABLE]

For every function $f:{\cal P}_{d}\rightarrow{\cal P}_{d}$, we define the operator $*$ s.t for every $m\in{\cal M}_{d}$, for every $\pi\in{\cal P}_{d}$, $f*m(\pi)=m(f^{-1}(\pi))$. In words, $f*m$ is the push-forward measure of $m$ by $f$. Further, we will also consider square matrices indexed by elements in ${\mathcal{P}}_{d}$. For such a matrix $K$ and an element $m\in{\cal M}_{d}$, we define $Km(\pi):=\sum_{\pi^{\prime}\in{\cal P}_{d}}K(\pi,\pi^{\prime})m(\pi^{\prime})$.

Define

[TABLE]

i.e., $X(\Pi)$ is the empirical measure associated to the “sample” $\Pi_{1},\dots,\Pi_{l^{\epsilon}}$. In the following, we define

[TABLE]

will be referred to as the (empirical) genetic partition probability measure of the population, conditional on the $l^{\epsilon}$ loci to be located at $L^{\epsilon}$. We also define

[TABLE]

will be referred to as the (empirical) genetic partition probability measure of the population.

The genetic distance between $i$ and $j$ at time $t$ can then be expressed in terms of $\xi^{\epsilon}_{t}$ as follows:

[TABLE]

In the following, we identify the process $(\xi^{\epsilon}_{t},t\geq 0)$ to a process in the set of the càdlàg functions from $\mathbb{R}^{+}$ to $\mathbb{R}^{Bell_{d}}$, equipped with the standard Skorokhod topology.

3.3 Convergence of the genetic partition probability measure

Following Remark 2.3, for every $k\in\{1,\dots,l^{\epsilon}\}$, the process $(\Pi_{k}^{\epsilon,{L}^{\epsilon}}(t);t\geq 0)$ – the partition at locus $k$ – obeys the following dynamics:

(reproduction event) $j$ is merged in the block containing $i$ at rate $M_{ij}\frac{1}{\epsilon n^{\epsilon}_{j}}$. 2. 2.

(mutation) Individual $i$ takes on a new type at rate $b_{\infty}$.

The generator associated to the allelic partition at locus $k$ is then given by

[TABLE]

Recall that the expression of the generator associated to the allelic partition at a given locus $k$ is independent on the position of locus $k$ and on $\lambda$. Also recall that $\epsilon n_{i}^{\epsilon}\to N_{i}$ and that, by definition, $\forall k,p\in E,\ \tilde{M}_{kp}=M_{pk}/N_{k}$. Thus

[TABLE]

Direct computations yield that ${}^{t}G$, the transpose of the matrix $G$ satisfies

[TABLE]

In the light of (7), the following theorem can be interpreted as an ergodic theorem. We show that the (dynamical) empirical measure constructed from the allelic partitions along the chromosome converges to the probability measure of a single locus. Although in the IBM the different loci are linked and do not fix independently (as already mentioned in Remark 2.4), as the number of loci tends to infinity, they become decorrelated. In the large population - dense site limit, the following result indicates that the model behaves as if infinitely many loci evolved independently according to the (one-locus) Moran model with generator $G$ provided in (7).

In the following $``\Longrightarrow^{\prime\prime}$ indicates the convergence in distribution. Also, we identify $(\xi_{t}^{\epsilon};\ t\geq 0)$ to a function from $\mathbb{R}^{+}$ to $\mathbb{R}^{Bell_{d}}$; and convergence in the weak topology means that for every $T>0$, the process $(\xi_{t}^{\epsilon};\ t\in[0,T])$ converges in the Skorohod topology $D([0,T],\mathbb{R}^{Bell_{d}})$.

Theorem 3.1** (Ergodic theorem along the chromosome).**

Assume that $\xi^{\epsilon}_{0}$ is deterministic and there exists a probability measure $P^{0}\in{\cal M}_{d}$ such that the following convergence holds:

[TABLE]

Then

[TABLE]

where $P$ solves the forward Kolmogorov equation associated to the aforementioned Moran model, i.e.,

[TABLE]

with initial condition $P_{0}=P^{0}$ and where ${}^{t}G$ denotes the transpose of $G$ (see (8)).

4 Proof of Theorem 3.1

The idea behind the proof is to condition on ${\cal L}^{\epsilon}=L^{\epsilon}$, and then decompose the Markov process $(<\xi^{\epsilon,L_{\epsilon}}_{t},v>;t\geq 0)$ into a drift part and a Martingale part. We show that the drift part converges to the solution of the Kolmogorov equation alluded to in Theorem 3.1 and that the Martingale part vanishes when $\epsilon\to 0$. The main steps of the computation are outlined in the next subsection. We leave technical details (tightness and second moment computations) until the end of the section.

4.1 Main steps of the proof

Fix $L^{\epsilon}\in[0,1]^{l^{\epsilon}}$. Recall the definition of ${\mathbb{G}}^{\epsilon,L^{\epsilon}}$, the generator of the process $(\Pi^{\epsilon,L^{\epsilon}}(t);\ t\geq 0)$, given in (4). Let $f:\mathbb{R}\to\mathbb{R}$ be a Borel bounded function and fix $v\in{\cal M}_{d}$. Define $h(\Pi)\ =\ f(\left<X(\Pi),v\right>)$. Then, it is straightforward to see from (4) that

[TABLE]

where in the first line $\mathbb{E}_{\lambda,L^{\epsilon},j}$ is the expected value taken with respect to the random variable $S$, distributed as ${\cal F}_{j}^{L^{\epsilon},\lambda}$ as defined in Definition 2.1. In the second line, $\mathbb{E}_{l^{\epsilon}}$ is the expected value is taken with respect to $K$, distributed as a uniform random variable on $\{1,\dots,l^{\epsilon}\}$.

Lemma 4.1**.**

Define $v\in{\cal M}_{d}$, $L^{\epsilon}\in[0,1]^{l^{\epsilon}}$, $g(\Pi):=\left<X(\Pi),v\right>$. Then ${\mathbb{G}}^{\epsilon,L^{\epsilon}}g(\Pi)\ =\ \left<{}^{t}G^{\epsilon}X(\Pi),v\right>$ where ${}^{t}G^{\epsilon}$ is the transpose of $G^{\epsilon}$ – the generator of the allelic partition at a single locus as defined in (6) – i.e.,

[TABLE]

Proof.

We define the two following signed measures:

[TABLE]

In words, $\partial_{ij}^{S}X(\Pi)$ is the change in the genetic partition measure $X(\Pi)$ if we merge $j$ in the block of $i$ at every locus in $S$, and $\partial_{i}^{K}X(\Pi)$ is the change in $X(\Pi)$ if we single out element $i$ at locus $K$. Using those notations, for our particular choice of $g$, (10) writes

[TABLE]

We now show that for every $v\in{\cal M}_{d}$,

[TABLE]

We only prove the first identity. The second one can be shown along the same lines. Again, we let $\Pi_{k}$ be the $k^{th}$ coordinate of $\Pi$. By definition, the vector ${\cal I}_{ij}^{S}(\Pi)$ is only modified at the coordinates belonging to $S$, and thus

[TABLE]

Secondly, for every $j\in E$,

[TABLE]

As $S$ is distributed as ${\cal F}_{j}^{L^{\epsilon},\lambda}$, we can use the fact that ${\mathbb{P}}\left(k\in{\cal F}_{j}^{L^{\epsilon},\lambda}\right)\ =\ M_{ij}\frac{1}{n^{\epsilon}_{j}}$ (see Remark 2.3), and then

[TABLE]

Furthermore, by applying (15) for every $\pi^{\prime}\in\mathcal{I}_{ij}^{-1}(\pi)$ and then taking the sum over every such partitions, we get

[TABLE]

This completes the proof of (13). From this result, we deduce that

[TABLE]

This completes the proof of Lemma 4.1. ∎

For every $L^{\epsilon}\in[0,1]^{l^{\epsilon}}$, for every $v\in{\cal M}_{d}$, define

[TABLE]

Since $\left<\xi_{t}^{\epsilon,L^{\epsilon}},v\right>$ is bounded, the previous result implies that $M^{\epsilon,L^{\epsilon},v}$ is a martingale with respect to $({\cal H}^{L^{\epsilon}}_{t})_{t\geq 0}$, the filtration generated by $(\Pi^{\epsilon,L^{\epsilon}}(t);t\geq 0)$. Further, the semi-martingale $\left<\xi_{t}^{\epsilon,L^{\epsilon}},v\right>$ admits the following decomposition:

[TABLE]

Lemma 4.2**.**

For every $v\in{\cal M}_{d}$, for every $L^{\epsilon}\in[0,1]^{l^{\epsilon}}$,

[TABLE]

with

[TABLE]

and where $\left<M^{\epsilon,L^{\epsilon},v}\right>_{t}$ denotes the quadratic variation of $M^{\epsilon,L^{\epsilon},v}$.

Proof.

For $v\in{\cal M}_{d}$, define $h(\Pi)=\left<X(\Pi),v\right>^{2}$. Then, by (10)

[TABLE]

Since

[TABLE]

the previous identities yield

[TABLE]

where $g(\Pi)=\left<X(\Pi),v\right>$. As a consequence

[TABLE]

is a martingale. Further using Itô’s formula, the process

[TABLE]

is also a martingale. Combining the two previous results completes the proof of Lemma 4.2. ∎

Proposition 4.3**.**

[TABLE]

where the expected value is taken with respect to the random variable ${\cal L}^{\epsilon}$.

Proposition 4.4**.**

For $T>0$, the family of random variables $(\xi^{\epsilon};\ {\epsilon>0})$ is tight in the weak topology $D([0,T],\mathbb{R}^{{Bell}_{d}})$.

We postpone the proof of Propositions 4.3 and 4.4 until Sections 4.2 and 4.3 respectively.

Proof of Theorem 3.1 based on Proposition 4.3 and 4.4.

Since $(\xi^{\epsilon};\ {\epsilon>0})$ is tight, we can always extract a subsequence converging in distribution (for the weak topology) to a limiting random measure process $\xi$. We will now show that $\xi$ can only be the solution of the Kolmogorov equation alluded to in Theorem 3.1. From (11), for every probability measure $m$ on ${\cal P}_{d}$, for every $v\in{\cal M}_{d}$,

[TABLE]

where $||v||_{\infty}:=\max_{\pi\in{\cal P}_{N}}v(\pi)$. Since as $\epsilon\to 0$, $n_{i}^{\epsilon}\epsilon\to N_{i}$, the term between parentheses also converges, and thus the RHS is uniformly bounded in $\epsilon$. Finally, the bounded convergence theorem implies that for every $v\in{\cal M}_{d}$,

[TABLE]

where we used the fact that ${}^{t}G^{\epsilon}m\ \to\ ^{t}Gm$ for every $m\in{\cal M}_{d}$ (where $G$ is defined as in Theorem 3.1). On the other hand, since

[TABLE]

Lemma 4.2 and Proposition 4.3 imply that

[TABLE]

which ends the proof of Theorem 3.1. ∎

4.2 Proof of Proposition 4.3

Our first step in proving Proposition 4.3 is to prove the following result.

Lemma 4.5**.**

$\forall\ j\ \in E$, $\forall\lambda>0$,

[TABLE]

Before turning to the proof of this result, we recall the definition of the ancestral recombination graph (ARG) (see also Hudson (1983), Griffiths (1981), Griffiths (1991)) for the case of two loci. Fix ${L}^{\epsilon}=\{\ell_{1},\dots\ell_{l^{\epsilon}}\}$ the positions of the loci in the chromosome, $\lambda$ the recombination rate, and choose two loci $k$ and ${k^{\prime}}$ among the $l^{\epsilon}$ loci. In order to compute the probability that for both loci the allele from the migrant is fixed in the host population – $\mathbb{P}(k,k^{\prime}\in{\cal F}_{j}^{L^{\epsilon},\lambda})$ – we follow backwards in time the genealogy of the corresponding alleles carried by a reference individual in the present population, assuming that a migration event occurred in the past (sufficiently many generations ago, so that we can assume that the present population is homogeneous).

More precisely, at locus $k$, we consider the ancestral lineage of a reference individual (chosen uniformly at random) in the extant population. We envision this lineage as a particle moving in $\{1,\dots,n^{\epsilon}_{j}\}$: time $t=0$ corresponds to the present, and the position of the particle at time $t$ – denoted by $A^{L^{\epsilon},\lambda,j}_{k}(t)$ – identifies the ancestor of locus $k$, $t$ units of time in the past (i.e., at locus $k$, the reference individual in the extant population inherits its genetic material from individual $A^{L^{\epsilon},\lambda,j}_{k}(t)$ at time $-t$) (see Figure 5).

The recombination rate between the two loci, $k$ and $k^{\prime}$, $r^{L^{\epsilon},\lambda}_{k,k^{\prime}}$ corresponds to the probability that there is at least one Poisson point between $\ell_{k}$ and $\ell_{k^{\prime}}$ and the two fragments are inherited from different parents and is given by

[TABLE]

${\bf A}^{L^{\epsilon},\lambda,j}=(A^{L^{\epsilon},\lambda,j}_{k},\ A^{L^{\epsilon},\lambda,j}_{k^{\prime}})$ defines a $2$-dimensional stochastic process on $\{1,\dots,n_{j}^{\epsilon}\}$. At time [math], the two particles have the same position (they coincide at a randomly chosen individual as in Figure 5) and then evolve according to the following dynamics:

• •

When both particles are occupying the same location $z$, the group splits into two at rate $r^{L^{\epsilon},\lambda}_{k,k^{\prime}}$ (see (17)). Forward in time this corresponds to a reproduction event where $z$ is replaced by the offspring of $x$ and $y$. Each individual $x$ reproduces at rate 1 (chooses a random partner $y$), and with probability $1/n_{j}^{\epsilon}$ his offspring replaces individual $z$. There are $n_{j}^{\epsilon}$ possible choices for $x$. Following (17), the probability that both loci are inherited from different parents is $r^{L^{\epsilon},\lambda}_{k,k^{\prime}}$, so the rate of fragmentation for loci $k,k^{\prime}$ is given by $n_{j}^{\epsilon}\ .\ \frac{1}{n_{j}^{\epsilon}}\ .\ r^{L^{\epsilon},\lambda}_{k,k^{\prime}}$.

• •

When the two particles are occupying different positions, they jump to the same position at rate $2/n_{j}^{\epsilon}$. Forwards in time, this corresponds to a reproduction event where the individual located at $A^{L^{\epsilon},\lambda,j}_{k}$ (resp. $A^{L^{\epsilon},\lambda,j}_{k^{\prime}}$) replaces the one at $A^{L^{\epsilon},\lambda,j}_{k^{\prime}}$ (resp. $A^{L^{\epsilon},\lambda,j}_{k}$), and the offspring inherits the allele at locus $k^{\prime}$ (resp. $k$) from this parent. A reproduction event where the individual located at $A^{L^{\epsilon},\lambda,j}_{k}$ (resp. $A^{L^{\epsilon},\lambda,j}_{k^{\prime}}$) replaces the one at $A^{L^{\epsilon},\lambda,j}_{k^{\prime}}$ (resp. $A^{L^{\epsilon},\lambda,j}_{k}$) occurs at rate $2/n_{j}^{\epsilon}$ (as the individual at $A^{L^{\epsilon},\lambda,j}_{k}$ –resp. $A^{L^{\epsilon},\lambda,j}_{k^{\prime}}$– can be the mother or the father); and the probability that the offspring inherits the locus $k^{\prime}$ (resp. $k$) from this parent is $1/2$. The total rate of coalescence is $2\ .\ \frac{2}{n_{j}^{\epsilon}}\ .\ \frac{1}{2}$.

Since we assume that the migration event occurred far back in the past, the following duality relation holds:

[TABLE]

In other words, assuming that the migrant is labelled 1, the set on the RHS corresponds to the set of loci inheriting their genetic material from the migrant.

Proof of Lemma 4.5.

Define $(Y^{L^{\epsilon},\lambda,j}(t):=1_{A^{L^{\epsilon},\lambda,j}_{k}(t)=A^{L^{\epsilon},\lambda,j}_{k^{\prime}}(t)};t\geq 0)$. It is easy to see from the previous description of the dynamics that $Y$ is a Markov chain on $\{0,1\}$ with the following transition rates:

[TABLE]

and further

• •

conditional on $Y^{L^{\epsilon},\lambda,j}(t)=1$, the two lineages $(A^{L^{\epsilon},\lambda,j}_{k}(t),A^{L^{\epsilon},\lambda,j}_{k^{\prime}}(t))$ occupy a common position that is distributed as a uniform random variable on $\{1,\dots,n_{j}^{\epsilon}\}$.

• •

conditional on $Y^{L^{\epsilon},\lambda,j}(t)=0$, $(A^{L^{\epsilon},\lambda,j}_{k}(t),A^{L^{\epsilon},\lambda,j}_{k^{\prime}}(t))$ are distinct and are distributed as a two uniformly sampled random variables (without replacement) on $\{1,\dots,n_{j}^{\epsilon}\}$.

We have:

[TABLE]

Furthermore, it is straightforward to show that

[TABLE]

From (18), we get that,

[TABLE]

One can then easily check that, $\exists\ \alpha>0$ such that, for every ${L}^{\epsilon}=\{\ell_{1},\dots\ell_{l^{\epsilon}}\}$,

[TABLE]

Thus,

[TABLE]

In addition, using the fact that $n_{j}^{\epsilon}=[N_{j}/\epsilon]$,

[TABLE]

∎

We are now ready to prove Proposition 4.3.

Proof of Proposition 4.3.

Using the definition given in Lemma 4.2,

[TABLE]

To bound the second term on the RHS, we note that, by definition, ${\cal S}_{i}^{k}(\Pi)$ and $\Pi$ only differ in one component, so from the definition of $\partial_{i}^{K}X(\Pi)$ (see (12)), it is not hard to see that

[TABLE]

It follows that,

[TABLE]

Since $l^{\epsilon}\to\infty$ as $\epsilon\to 0$, this term converges and can be bounded from above, uniformly in $\Pi$ and $\epsilon\in(0,1)$. Note that this bound does not depend on the choice of ${L}^{\epsilon}$.

For the second term on the RHS, we simply note that expanding $\frac{1}{\epsilon}\mathbb{E}_{\lambda,L^{\epsilon},j}\left(\left<\partial_{ij}^{S}X(\Pi),v\right>^{2}\right)$ (see (14)), yields a sum of four terms that can be upper bounded by

[TABLE]

where $p_{1}$ and $p_{2}$ are alternatively replaced by $\{\pi\},{\cal I}_{ij}^{-1}(\pi)$ with $\pi\in{\cal P}_{d}$. Finally, $\forall L^{\epsilon}\in[0,1]^{l^{\epsilon}}$,

[TABLE]

randomising the positions of the loci and using Lemma 4.5 the term on the RHS also converges and can also be bounded from above, which completes the proof. ∎

Remark 4.6** (Magnitude of the stochastic fluctuations).**

Lemma 4.1 and the proof Proposition 4.3 entail that:

[TABLE]

where $C,C^{\prime}$ are constants. This suggests that the order of magnitude of the fluctuations should be of the order of $\max(\sqrt{\epsilon\log(1/\epsilon)},\sqrt{(}1/{l^{\epsilon}}))$.

In Yamaguchi and Iwasa (2013), the authors proposed a diffusion approximation (only for the case of two subpopulations). Their approximation is based on the simplifying hypothesis that loci are fixed independently on each other – the number of fixed loci (after each migration event) follows a binomial distribution–, and the hypothesis that the number of loci $l$ is s.t. $l>>\frac{1}{\epsilon}$. They found that the magnitude of the stochastic fluctuations was $\sqrt{\epsilon}$.

In summary, the previous heuristics suggest that taking into account correlations between loci increases the magnitude of the stochastic fluctuations.

4.3 Tightness: Proof of Proposition 4.4

We follow closely Fournier and Méléard (2004). It is sufficient to prove that for every $v\in{\cal M}_{d}$, the projected process $\left(\left<\xi^{\epsilon},v\right>;\ \epsilon>0\right)$ is tight. To this end, we use Aldous criterium (see Aldous (1989)). In the following, we define

[TABLE]

We first note that

[TABLE]

which implies that that for every deterministic $t\in[0,T]$, the sequence of random variables $(\left<\xi^{\epsilon}_{t},v\right>;\epsilon>0)$ is tight. Thus, the first part of Aldous criterium is sastified. Next, fix $\delta>0$, and take two stopping times $\tau^{\epsilon}$ and $\sigma^{\epsilon}$ with respect to $({\cal H}^{\epsilon}_{t})_{t\geq 0}$ the filtration generated by $(\Pi^{\epsilon,{\cal L}^{\epsilon}}_{t},t\geq 0)$, such that $0\leq\tau^{\epsilon}\leq\sigma^{\epsilon}\leq\tau^{\epsilon}+\delta\leq T$. Since $\left<\xi_{t}^{\epsilon},v\right>\ =\ M^{\epsilon,v}_{t}\ +\ B^{\epsilon,v}_{t},$ it is enough to show that the quantities

[TABLE]

are bounded from above by two functions in $\delta$ (uniformly in the choice of $\tau^{\epsilon},\sigma^{\epsilon}$ and $\epsilon$) going to [math] as $\delta$ go to [math].

The rest of the proof is dedicated to proving those two inequalities. We start with the martingale part. First,

[TABLE]

Recall that $\forall\ L^{\epsilon}\in[0,1]^{l^{\epsilon}}$, $M^{\epsilon,L^{\epsilon},v}$ is a martingale. Thus, $M^{\epsilon,v}$ is a martingale with respect to $({\cal G}_{t}^{\epsilon})_{t\geq 0}=({\cal H}^{\epsilon}_{t})_{t\geq 0}\vee\sigma({\cal L}^{\epsilon})$, where $({\cal H}^{\epsilon}_{t})_{t\geq 0}$ is the filtration generated by $(\Pi^{\epsilon,{\cal L}^{\epsilon}}_{t})$. As $({\cal H}^{\epsilon}_{t})_{t\geq 0}\subset({\cal G}^{\epsilon}_{t})_{t\geq 0}$, $\tau^{\epsilon}$ and $\sigma^{\epsilon}$ are also stopping times for the filtration $({\cal G}^{\epsilon}_{t})_{t\geq 0}$, so that

[TABLE]

where $m^{\epsilon,L^{\epsilon},v}(\Pi)$ was defined in Lemma 4.2 and where the second line follows from the fact that $\tau^{\epsilon}$ and $\sigma^{\epsilon}$ are stopping times for the filtration $({\cal G}^{\epsilon}_{t})_{t\geq 0}$.

If there exists $C_{1}$ such that

[TABLE]

then,

[TABLE]

thus showing the desired inequality for the martingale part $M^{\epsilon,v}$. To prove (22), we recall the definition of $m^{\epsilon,L^{\epsilon},v}(\Pi)$,

[TABLE]

The second term on the RHS can be bounded as in the proof of Proposition 4.3 (see (20)). For the first term on the RHS, we use the bound given by (LABEL:majoration-rhs). We only need to prove that $\frac{1}{(l^{\epsilon})^{2}\epsilon}\ \mathbb{E}_{\lambda,L^{\epsilon},j}(\sum\limits_{k=1}^{l^{\epsilon}}\sum\limits_{k^{\prime}=1}^{l^{\epsilon}}1_{k\in S}1_{k^{\prime}\in S})$ is bounded. Using (19),

[TABLE]

so (22) is proved.

We now turn to the drift part. First, for every $L^{\epsilon}\in[0,1]^{l^{\epsilon}}$,

[TABLE]

We already showed in (16), that the integrand on the RHS is uniformly bounded in $\epsilon$. Thus, there exists $C_{2}$ such that, for every $L^{\epsilon}\in[0,1]^{l^{\epsilon}}$:

[TABLE]

So,

[TABLE]

which is the desired inequality. This completes the proof of Proposition 4.4.

Remark 4.7**.**

Notice that the tightness (and the convergence) does not depend on the recombination rate. However, for small values of $\lambda$, or if $L^{\epsilon}$ is such that the positions of the loci are all very close to each other, correlations between loci are very high. This means that, when a migration event takes place, either no locus will be fixed (with high probability), or almost all loci from the migrant will be fixed. Therefore, if we let $\lambda\to 0$ be the process of the genetics distances converges to a process that increases continuously (due to mutation) and has negative jumps (due to migration events). See Figure 6 for a numerical simulation.

5 Proof of Theorem 1.1 and more

In this section we state and prove a stronger version of Theorem 1.1 (see Theorem 5.1 below). As in Theorem 1.1, we consider, for each pair of subpopulations $i,j\in E$, $S^{i}$ and $S^{j}$, two independent random walks on $E$ starting respectively from $i$ and $j$ and whose transition rate from $k$ to $p$ is equal to $\tilde{M}_{kp}$. We have the following generalization of of Theorem 1.1.

Theorem 5.1**.**

Assume that

• •

At time [math], in the IBM, subpopulations are homogeneous and that, the genetic partition measure of the population (in the associated PBM) is given by $\xi^{\epsilon}_{0}$, a deterministic probability measure in ${\cal P}_{d}$.

• •

There exists a probability measure $P^{0}\in{\cal M}_{d}$ such that the following convergence holds:

[TABLE]

For every $t\geq 0$, define

[TABLE]

where $\tau_{ij}=\inf\{t\geq 0\ :\ S^{i}(t)=S^{j}(t)\}.$ Then,

[TABLE]

In particular,

[TABLE]

Proof.

We start by proving that,

[TABLE]

where $(D_{t}(i,j);t\geq 0)$ is the deterministic process defined in Theorem 1.1.

From equation (5) and Theorem 3.1 we get that $\forall i,j\in E$, $(d_{t}^{\epsilon}(i,j);\ t\geq 0)$ converges in distribution in the weak topology to $(1-P_{t}(\pi\in{\cal P}_{d},i\sim_{\pi}j);\ t\geq 0)$. It remains to show that this expression is identical to the one provided in Theorem 1.1. This is done in a standard way by using the graphical representation associated to the one-locus Moran model whose generator is specified by $G$ (defined in (7)). It is well known that such a Moran model is encoded by a graphical representation that is generated by a sequence of independent Poisson Point Processes as follows:

• •

$B^{i}$, with intensity measure $b_{\infty}dt$, that corresponds to mutation events at site $i$. If $b^{i}\in B^{i}$, at $(i,b^{i})$ we draw a $\star$ in the graphical representation (Figure 7(a)).

• •

$T^{i,j}$, with intensity measure $\tilde{M}_{ji}dt$, that corresponds to reproduction events, where $j$ is replaced by $i$. If $t^{i,j}\in T^{i,j}$, we draw an arrow from $(i,t^{i,j})$ to $(j,t^{i,j})$ in the graphical representation to indicate that lineage $j$ inherits the type of lineage $i$.

We now give a characterisation of the dual process starting at $t$. We define $(S^{1}_{t},S^{2}_{t},\dots,S^{d}_{t})$ a sequence of piecewise continuous functions $[0,t]\to E$, where $\forall i\in E,\ S^{i}_{t}$ represents the ancestral lineage of individual $i$ (sampled at time $t$). $S^{i}_{t}(t)=i$ and as time proceeds backwards, each time $S^{i}_{t}$ encounters the tip of an arrow it jumps to the origin of the arrow. It is not hard to see that $S^{i}_{t}$ is distributed as a random walk started at $i$ and with transition rates from $k$ to $p$ equal to $\tilde{M}_{kp}$ and that $(S^{1}_{t},S^{2}_{t},\dots,S^{d}_{t})$ are distributed as coalescing random walks running backwards in time, i.e. they are independent when appart and become perfectly correlated when meeting each other. In Figure 7(c), $S^{1}_{t},S^{2}_{t},\dots,S^{d}_{t}$ are represented in different colours.

Define $\tau^{t}_{ij}=\inf\{s\geq 0,\ S^{i}_{t}(t-s)=S^{j}_{t}(t-s)\}$. By looking carefully at Figures 7(b) and 7(c), we let the reader convince herself that two individuals $i$ and $j$ have the same type at time $t$ iff:

• (i)

$\tau^{t}_{ij}\leq t$ and there are no $\star$ in the paths of $S^{i}_{t}$ and $S^{j}_{t}$ before $\tau^{t}_{ij}$, or

• (ii)

$\tau^{t}_{ij}\geq t$ and $S^{i}(0)\sim_{\pi_{0}}S^{j}(0)$ and $S^{i}_{t}$ and $S^{j}_{t}$ their is no $\star$ in their paths, where $\pi_{0}$ is the initial genetic partition of the metapopulation, that is random partition of law $P^{0}$.

From here, it is easy to check that:

[TABLE]

As $\forall i,j\in E$, $(D_{t}(i,j))$ is continuous, the fact that $(d_{t}^{\epsilon}(i,j))$ converges in distribution (in the weak topology) to $(D_{t}(i,j))$ (24) implies (by the continuous mapping theorem) that $(d_{t}^{\epsilon}(i,j))$ converges to $(D_{t}(i,j))$ in the sense of finite dimensional distributions, as $\epsilon\to 0$.

This result, combined with Theorem 3, also implies that:

[TABLE]

The fact that $\lim_{t\to\infty}D_{t}(i,j)\ =\ 1-\mathbb{E}(e^{-2b_{\infty}\tau_{ij}})$ is a direct consequence of the definition of $(D_{t}(i,j);t\geq 0)$ and the dominated convergence theorem.

This completes the proof of Theorem 1.1. ∎

6 An example: a population with a geographic bottleneck

Fix $d\in\mathbb{N}\setminus\{0\}$. We let ${\cal G}_{1}$ and ${\cal G}_{2}$ be two complete graphs of $d$ vertices. We link the two graphs ${\cal G}_{1}$ and ${\cal G}_{2}$ by adding an extra edge $(v_{1},v_{2})$, where $v_{k},k=1,2$ is a given vertex in ${\cal G}_{k}$. We call ${\cal G}$ the resulting graph. We equip ${\cal G}$ with the following migration rates: if $i$ is connected to $j$, then $M_{ij}=1/d$ (so that the emigration rate from any vertex $i$ is $1$ if $i\neq v_{1},v_{2}$ and $1+\frac{1}{d}$ otherwise). We also assume that $N_{i}=1$, so that $\tilde{M}_{ij}=1/d$.

We think of ${\cal G}$ as two well-mixed populations connected by a single geographic bottleneck.

Theorem 6.1**.**

Fix $c>0$, $b_{\infty}=\frac{c}{d}$. Then for any two neighbours $i,j\in{\cal G}$

[TABLE]

Before going into the details of the proof, we present some heuristics for the formulae. First, if $i$ and $j$ belong to the same subgraph, we can assume that most of the time the two random walks hit each other before hitting the other subgraph. So we can consider $\tilde{S}^{i}$ and $\tilde{S}^{j}$, two random walks on a complete graph with $d$ vertices. If $J$ is the number of jumps made by $\tilde{S}^{i}$ and $\tilde{S}^{j}$ before hitting each other, $J$ follows a geometric distribution of parameter $1/d$, which properly renormalized converges to an exponential distribution of parameter $1$ (when $d\gg 1$). In addition, the mean time between two consecutive jumps (of $S^{i}$ or $S^{j}$) is $2\times 1/d$ so the distribution of $\tau_{ij}$ can be approximated by an exponential distribution of parameter $2/d$. If $e$ is an exponentially distributed random variable, with parameter $2/d$, then $\mathbb{E}\left(\exp(-2b_{\infty}e)\right)=1/(1+2b_{\infty}/d)=1/(1+c)$ which gives the desired result. Second, $i=v_{1}$ and $j=v_{2}$, with probability $1/(d+1)\simeq 1/d$, the first jump is from $v_{1}$ to $v_{2}$ (or $v_{2}$ to $v_{1}$), so the two random walks hit very fast, and the genetic distance is close to [math]. Otherwise, each random walk “gets lost” in its subgraph and the hitting time becomes very large. In that case, the genetic distance is approximatively $1$.

Proof.

We give a brief sketch of the computations since the method is rather standard. We start with some general considerations. Consider a general meta-population with $\bar{d}$ subpopulations. Define $a(i,j)\ =\ \mathbb{E}\left(\exp(-2b_{\infty}\tau_{ij})\right)$. By conditioning on every possible move of the two walks on the small time interval $[0,dt]$, it is not hard to show that the $a(i,j)$’s satisfy the following system of linear equations: $\forall i\in\{1,\dots,\bar{d}\},\ a(i,i)=1$ and $\forall i,j\in\{1,\dots,\bar{d}\}$ with $i\neq j$:

[TABLE]

Let us now go back to our specific case (in particular $\bar{d}=2d$). We distinguish between two types of points: the boundary points (either $v_{1}$ or $v_{2}$), and the interior points of the subgraphs ${\cal G}_{1}$ and ${\cal G}_{2}$ (points that are distinct from $v_{1}$ and $v_{2}$). For $(i,j)$, with $i\neq j$, we say that $(i,j)$ is of type

• •

$(II)$ if the vertices belong to the interior of the same subgraph (either ${\cal G}_{1}$ or ${\cal G}_{2}$).

• •

$(I\bar{I})$ if the vertices belong to the interior of distinct subgraphs.

• •

$(IB)$ if one of the vertex is in the interior of a subgraph, and the other vertex belongs to the boundary point of the same subgraph.

• •

$(I\bar{B})$, $(B\bar{B})$ are defined analogously.

By symmetry, $a(i,j)$ is invariant in each of those classes of pairs of points. We denote by $a(II)$ the value of $a(i,j)$ for $(i,j)$ in $(II)$. $a(I\bar{I}),a(IB),a(I\bar{B}),a(B\bar{B})$ are defined analogously. From this observation, we can inject those quantities in (26): this reduces the dimension of the linear problem from $\bar{d}(\bar{d}-1)$ to only $5$. The system can then be solved explicitly and straightforward asymptotics yield Theorem 6.1.

∎

Acknowledgements

We would like to thank the Editor and the reviewers for their useful comments on the previous version of the manuscript. We also thank Florence Débarre and Amaury Lambert for helpful discussions.