A dual process for the coupled Wright-Fisher diffusion

TL;DR

This paper introduces a dual ancestral process for the coupled Wright-Fisher diffusion, capturing multi-locus, multi-allelic genetic interactions with explicit formulas for special cases.

Contribution

It derives a dual ancestral process for the coupled Wright-Fisher diffusion, incorporating pairwise interactions and providing explicit stationary and transition rates in special cases.

Findings

Derived the dual ancestral process with coalescence, mutation, and branching events.

Explicit stationary density and transition rates for two loci with selection and mutation.

Extended the ancestral selection graph framework to coupled multi-locus models.

Abstract

The coupled Wright-Fisher diffusion is a multi-dimensional Wright-Fisher diffusion for multi-locus and multi-allelic genetic frequencies, expressed as the strong solution to a system of stochastic differential equations that are coupled in the drift, where the pairwise interaction among loci is modelled by an inter-locus selection. In this paper, an ancestral process, which is dual to the coupled Wright-Fisher diffusion, is derived. The dual process corresponds to the block counting process of coupled ancestral selection graphs, one for each locus. Jumps of the dual process arise from coalescence, mutation, single-branching, which occur at one locus at the time, and double-branching, which occur simultaneously at two loci. The coalescence and mutation rates have the typical structure of the transition rates of the Kingman coalescent process. The single-branching rate not only contains the one-locus selection parameters in a form that generalises the rates of an ancestral selection graph, but it also contains the two-locus selection parameters to include the effect of the pairwise interaction on the single loci. The double-branching rate reflects the particular structure of pairwise selection interactions of the coupled Wright-Fisher diffusion. Moreover, in the special case of two loci, two alleles, with selection and parent independent mutation, the stationary density for the coupled Wright-Fisher diffusion and the transition rates of the dual process are obtained in an explicit form.

Equations

$${\mathbf{X}}^{(l)}(t)=(X_1^{(l)}(t),\dots ,X_{M_l}^{(l)}(t){)}^{\mathrm{T}}$$

$$\mathbf{X}(t)=({\mathbf{X}}^{(1)}(t),\dots ,{\mathbf{X}}^{(L)}(t){)}^{\mathrm{T}}.$$

$$\text{d}\mathbf{X}(t)=\mu (\mathbf{X}(t))\text{d}t+D(\mathbf{X}(t))\nabla V(\mathbf{X}(t))\text{d}t+D^{\text{nicefrac}12}(\mathbf{X}(t))\text{d}\mathbf{W}(t),$$

$$\mu (\mathbf{x})=\left(\begin{array}{c}{\mu }^{(1)}({\mathbf{x}}^{(1)})\\ ⋮\\ {\mu }^{(1)}({\mathbf{x}}^{(\mathbf{L})})\end{array}\right),D(\mathbf{x})=\left(\begin{array}{ccc}{D}^{(1)}({\mathbf{x}}^{(1)})& & \\ & \ddots & \\ & & D^{(L)}({\mathbf{x}}^{(\mathbf{L})})\end{array}\right),$$

$$\text{d}{\mathbf{X}}^{(l)}(t)={\mu }^{(l)}({\mathbf{X}}^{(l)}(t))\text{d}t+D^{(l)\text{nicefrac}12}({\mathbf{X}}^{(l)}(t))\text{d}{\mathbf{W}}^{(l)}(t),$$

$$p(\bar{\mathbf{x}})=\frac{1}{Z}\pi (\bar{\mathbf{x}})e^{2V(\bar{\mathbf{x}})},$$

$$\mathbb{E}\left[F(\mathbf{X}(t),\text{n})|\mathbf{X}(0)=\text{x}\right]=\mathbb{E}\left[F(\text{x},\mathbf{N}(t))|\mathbf{N}(0)=\text{n}\right],$$

$$\mathcal{L}F(\cdot ,\text{n})(\text{x})={\mathcal{L}}^{D}F(\text{x},\cdot )(\text{n}),$$

$$V(\mathbf{x})={\mathbf{x}}^{\mathrm{T}}\mathbf{h}+\frac{1}{2}{\mathbf{x}}^{\mathrm{T}}J\mathbf{x},$$

$$g_i^{(l)}(\mathbf{x})=\sum _{k=1}^{M_l}{d}_{ik}^{(l)}({\text{x}}^{(l)}){\tilde{h}}_k^{(l)}(\mathbf{x}),\text{with}{\tilde{h}}_k^{(l)}(\mathbf{x})=h_k^{(l)}+\sum _{\begin{array}{c}r=1\\ r\ne l\end{array}}^L\sum _{m=1}^{M_r}{J}_{km}^{(lr)}{x}_m^{(r)}.$$

$${\mu }_i^{(l)}({\text{x}}^{(l)})=\sum _{j=1}^{M_l}[u_{ji}^{(l)}{x}_j^{(l)}-u_{ij}^{(l)}{x}_i^{(l)}].$$

$$d_{ij}^{(l)}({\text{x}}^{(l)})=x_i^{(l)}({\delta }_{ij}-x_j^{(l)})\text{with}{\delta }_{ij}=\{\begin{array}{l}1\text{ if }i=j,\\ 0\text{ if }i\ne j,\end{array}$$

$$\mathcal{S}=\left\{\mathbf{x}\in [0,1{]}^M\text{ s.t. }\sum _{i=1}^{M_i}{x}_i^{(l)}=1\forall l=1,\dots ,L,\right\},$$

$$\mathcal{L}f(\text{x})=\sum _{l=1}^L\left[\sum _{i=1}^{M_l}\left({\mu }_i^{(l)}({\text{x}}^{(l)})+g_i^{(l)}(\text{x})\right)\frac{\partial f}{\partial x_i^{(l)}}(\text{x})+\frac{1}{2}\sum _{i,j=1}^{M_l}{d}_{ij}^{(l)}({\text{x}}^{(l)})\frac{{\partial }^{2}f}{\partial x_i^{(l)}\partial x_j^{(l)}}(\text{x})\right],$$

$$\bar{\mathcal{S}}=\left\{\bar{\mathbf{x}}\in [0,1{]}^{M-L}\text{ s.t. }\sum _{i=1}^{M_l-1}{\bar{x}}_i^{(l)}\le 1\forall l=1,\dots ,L.\right\},$$

$$\bar{\text{x}}=({\bar{x}}_1^{(1)},\dots ,{\bar{x}}_{M_1-1}^{(1)},\dots ,{\bar{x}}_1^{(L)},\dots ,{\bar{x}}_{M_L-1}^{(L)}{)}^{\mathrm{T}}\in \bar{\mathcal{S}}$$

$$\text{x}={\left({\bar{x}}_1^{(1)},\dots ,{\bar{x}}_{M_1-1}^{(1)},1-\sum _{i=1}^{M_1-1}{\bar{x}}_i^{(1)},\dots ,{\bar{x}}_1^{(L)},\dots ,{\bar{x}}_{M_L-1}^{(L)},1-\sum _{i=1}^{M_L-1}{\bar{x}}_i^{(L)}\right)}^{\mathrm{T}}\in \mathcal{S}.$$

$$\pi (\bar{\text{x}})=\prod _{l=1}^L{\pi }_l({\bar{\text{x}}}^{(l)}),\text{with}{\pi }_l(\bar{\text{x}})=\prod _{i=1}^{M_l-1}({\bar{x}}_i^{(l)}{)}^{2u_i^{(l)}-1}{\left(1-\sum _{i=1}^{M_l-1}{\bar{x}}_i^{(l)}\right)}^{2u_{M_l}^{(l)}-1}.$$

$$p(\bar{\mathbf{x}})=\frac{1}{Z}\pi (\bar{\mathbf{x}})e^{2V(\bar{\mathbf{x}})},$$

$$\mathcal{L}F(\cdot ,\text{n})(\text{x})={\mathcal{L}}^{D}F(\text{x},\cdot )(\text{n})\text{x}\in \mathcal{S},\text{n}\in {\mathbb{N}}^M,$$

$$\mathcal{L}F(\cdot ,\text{n})(\text{x})=\sum _{\hat{\text{n}}}q(\text{n},\hat{\text{n}})\left[F(\text{x},\hat{\text{n}})-F(\text{x},\text{n})\right],$$

$$F(\text{x},\text{n})=\frac{1}{k(\text{n})}\prod _{l=1}^L\prod _{i=1}^{M_l}(x_i^{(l)}{)}^{n_i^{(l)}},$$

$$\sum _{\hat{\mathbf{n}}}q(\text{n},\hat{\text{n}})\mathbb{E}\left[F(\tilde{\mathbf{X}},\hat{\text{n}})-F(\tilde{\mathbf{X}},\text{n})\right]=0,$$

$$k(\text{n})=\mathbb{E}\left[\prod _{l=1}^L\prod _{i=1}^{M_l}({\tilde{X}}_i^{(l)}{)}^{n_i^{(l)}}\right].$$

$$\mathcal{L}F(\cdot ,\text{n})(\text{x})=\sum _{\hat{n}\ne n}q(\text{n},\hat{\text{n}})F(\text{x},\hat{\text{n}})+q(\text{n},\text{n})F(\text{x},\text{n}),$$

$$\mathcal{L}F(\cdot,\textbf{n})(\textbf{x})=\sum_{\hat{n}\neq n}q(\textbf{n},\hat{\textbf{n}})\left[F(\mathbf{x},\mathbf{\hat{n}})-F(\mathbf{x},\mathbf{n})\right]+\Big{(}\sum_{\hat{n}\neq n}q(\textbf{n},\hat{\textbf{n}})+q(\textbf{n},\textbf{n})\Big{)}F(\textbf{x},\textbf{n}).$$

$$\sum _{\hat{n}\ne n}q(\text{n},\hat{\text{n}})+q(\text{n},\text{n})=0.$$

$$V(\text{x})=h_1{x}_1+h_2{x}_2.$$

$$j(i)=\{\begin{array}{ll}2& \text{if }i=1\\ 1& \text{if }i=2\end{array}.$$

$${\mu }_i(\text{x})=u_i{x}_{j(i)}-u_{j(i)}{x}_i,$$

Full text

In this paper, a dual process to the coupled Wright-Fisher diffusion is derived, which contains transition rates corresponding to coalescence and mutation as well as single-locus selection and double-locus selection. The coalescence and mutation rates correspond to the typical transition rates of Kingman’s coalescent process. The single-locus selection rate not only contains the single-locus selection parameters in a form that generalises the rates for an ancestral selection graph, but it also contains the double-selection parameters to include the effect of the pairwise interaction on the single locus. The double-locus selection rate reflects the particular structure of pairwise interactions of the coupled Wright-Fisher diffusion.

Moreover, in the special case of two loci, two alleles, with selection and parent independent mutation, the stationary density for the coupled Wright-Fisher diffusion and the transition rates of the dual process are obtained in an explicit form.

1. Introduction

The coupled Wright-Fisher diffusion was introduced by Aurell, Ekeberg and Koski [2] with the purpose of analysing networks of loci in recombining populations of bacteria, or more precisely, detecting couples of loci co-evolving under strong selective pressure when the linkage disequilibrium is low across the genome. The model includes parent dependent mutation, interlocus selection and free recombination. Mutation is assumed to occur independently at each locus, while selection consists of first and second order selective interaction among loci. The biological significance of this type of interaction is found in [9].

The model considers $L$ different loci where, at each locus, a number of variants (alleles) is possible. The allele types at locus $l$ are labelled by $1,\dots,M_{l}$, thus assuming that the type space at each locus is finite.

The coupled Wright-Fisher diffusion is obtained as the weak limit of a sequence of discrete Wright-Fisher models characterised by the assumption that the evolution of the population at one locus is conditionally independent on the other loci given that the previous generation at each locus is known. Here we state the definition of the diffusion as solution of a system of stochastic differential equations, without reference to the underlying discrete model.

The coupled Wright-Fisher diffusion, $\mathbf{X}=\{\mathbf{X}(t),t\geq 0\}$, represents the evolution of the vector of all frequencies of allele types at each locus. Let

[TABLE]

represent the vector of frequencies at locus $l$, with $X_{i}^{(l)}(t)$ being the frequency of allele type $i$ at locus $l$, then

[TABLE]

The process $\mathbf{X}$ is the strong solution to the system of stochastic differential equations

[TABLE]

where $V$ is a specific quadratic function encoding the structure of the interactions, while the mutation vector $\mu$ and the diffusion matrix $D$ have the following block structure

[TABLE]

with $\mu^{(l)}:\mathbb{R}^{M_{l}}\rightarrow\mathbb{R}^{M_{l}}$ and $D^{(l)}:\mathbb{R}^{M_{l}}\rightarrow\mathbb{R}^{M_{l}\times M_{l}}$. The functions $V$, $\mu$ and $D$ are described in detail in the next section. The process $\mathbf{W}=(\mathbf{W}^{(1)},\dots,\mathbf{W}^{(L)})$ is a multidimensional Brownian motion with $\mathbf{W}^{(l)}$ having the same dimension of $\mathbf{X}^{(l)}$.

The system of SDEs (1.1) consists of $L$ systems of equations for $\mathbf{X}^{(1)},\dots,\mathbf{X}^{(L)}$, coupled by the drift term $D\,\nabla V$. Note that, if $\nabla V=0$, there is no interaction among the loci and the coupled Wright-Fisher diffusion consists of $L$ independent Wright-Fisher diffusions, that is, each $\mathbf{X}^{(l)}$ solves

[TABLE]

which is the SDE for a single-locus, multi-type Wright-Fisher diffusion with mutations. In fact, the selective interaction among different loci is completely described by the drift term $D\,\nabla V$.

An interesting feature of the coupled Wright-Fisher diffusion, addressed by Aurell et al. [2] as one of the main motivations for its introduction, is its stationary density which appeared, in a more general form, as a conjecture by Kimura over half a century ago. In [6], Kimura suggests a Wright-Fisher model for multi-locus and multi-allelic genetic frequencies and conjectures that the stationary density is of the form $\pi(\bar{\mathbf{x}})e^{m(\bar{\mathbf{x}})}$, where $\pi$ is the product of Dirichlet densities and $m$ is a generic mean fitness term. The coupled Wright-Fisher diffusion is constructed so that the quadratic function $V$ could replace the generic $m$. Indeed, under the assumption of parent independent mutations, the stationary density of the coupled Wright-Fisher diffusion is known up to a normalising constant $Z$, and corresponds to the one conjectured by Kimura with $m=2V$,

[TABLE]

see Section 2 for the definition of $\pi$ and $V$. Note that, since the sum of the frequencies at each locus is equal to one, some components can be written as function of the others, which makes the distribution degenerate. Therefore, in order to deal with densities, a reduced process must be considered, that is, the last component at each locus is omitted and $\bar{\mathbf{x}}$ will be used, in place of $\mathbf{x}$, to denote the states of the reduced process.

In this paper the main result concerns the derivation of a dual process for the coupled Wright-Fisher diffusion. See [5] for a complete survey on duality for Markov processes. Markov duality has proven to be a useful technique in population genetics, several duality relationships have been established between different types of Wright-Fisher diffusions, modelling the evolution of allele frequencies forward in time, and reverse-time genealogical processes. The simplest and best known duality relationship in this framework is the one between the classic Wright-Fisher diffusion and the block counting process of the Kingman’s coalescent, which is a moment duality between a diffusion process and a jump process. Generalisations of this fundamental dual relation lead to the construction of many genealogical process, e.g. the ancestral selection graph [7, 8]. For an overview on the role of duality in population genetics see [4]. As pointed out in [4], the great value of the dual process is that it describes the history of a sample, more precisely, the posterior distribution given an allelic configuration at the present time which is needed for statistical inference. For example, knowing the backward dynamics, or a good approximation of them, turns out to be useful to construct good proposal distributions in importance sampling algorithms, see [10].

The main result in this paper is Theorem 5.1, which provides a description of the transition rates of a pure jump Markov process, $\mathbf{N}=\{\mathbf{N}(t),t\geq 0\}$, that is dual to the coupled Wright-Fisher diffusion, $\mathbf{X}$, through the duality relationship

[TABLE]

where $F$ is a moment duality function which will be determined. The derivation use a generator approach as in [4] and [3]. It is based on the duality relationship of the infinitesimal generators

[TABLE]

where $\mathcal{L}$ is the generator of the coupled Wright-Fisher diffusion and $\mathcal{L}^{D}$ the unknown generator of the dual process. By proposing an appropriate duality function $F$, the generator $\mathcal{L}^{D}$ of the dual process can be identified, from which transition rates of the dual process are obtained. It is known that the duality relationship (1.4) implies (1.3) under certain conditions on $F$, see [5] for details.

The paper is outlined as follows. In Section 2 a background on the coupled Wright-Fisher diffusions is provided. Section 3 outlines the general generator approach to derive the dual process. In Section 4 the case of one locus, two allele types and parent independent mutations is considered. In this case the dual process is related to the ancestral selection graph, moreover, explicit formulas for the stationary density of the diffusion and the transition rates of the dual process are obtained. The main result is provided in Section 5 where the dual process is derived in the general multi-locus setting. The final Section 6 provides additional details in the case of two loci, two alleles, selection and parent independent mutations, more precisely, the transition rates of the dual process are expressed in terms of beta and confluent hypergeometric functions.

2. Preliminaries on the coupled Wright-Fisher diffusion

In this section the coupled Wright-Fisher diffusion is introduced and the explicit expression for its infinitesimal generator is provided. The notation in this section differs slightly from that in [2]. In [2], the frequency of the last allele type at each locus is omitted, being function of the other frequencies, whereas in this paper an expanded version of the diffusion is considered, which includes all the frequencies. Since the frequencies sum up to one the descriptions are equivalent. For our purpose we find the expanded version more convenient to work with.

For a given integer $L\geq 1$, the number of loci, let $M_{1},\dots,M_{L}$ be positive integers representing the number of alleles at each loci. Put $M=\sum_{l=1}^{L}M_{l}$. A vector $\textbf{x}\in\mathbb{R}^{M}$ is interpreted as the concatenation of $L$ vectors with lengths $M_{1},\dots,M_{L}$, i.e. $\textbf{x}=(\textbf{x}^{(1)},\dots,\textbf{x}^{(L)})^{\mathrm{T}}$ with $\textbf{x}^{(l)}\in\mathbb{R}^{M_{l}}\quad l=1,\dots L$, and the coordinate $i$ in vector $\textbf{x}^{(l)}$ is denoted by $x_{i}^{(l)}$. Similarly, a matrix $A\in\mathbb{R}^{M\times M}$ consists of $L^{2}$ blocks with dimensions $(M_{l}\times M_{r})_{l,r=1,\dots,L}$. The block at position $(l,r)$ is denoted by $A^{(lr)}$ and its component at position $(i,j)$ is denoted by by $A_{ij}^{(lr)}$. Furthermore by $\mathbf{e_{i}^{(l)}}$ it is denoted the unit vector in $\mathbb{R}^{M}$ with the $i^{th}$ component of its $l^{th}$ building vector being equal to $1$.

In the following, each of the terms appearing in (1.1) will be described, starting from the interaction drift term. The quadratic function $V:\mathbb{R}^{M}\rightarrow\mathbb{R}$ is given by

[TABLE]

where $\mathbf{h}\in\mathbb{R}^{M}_{+}$ and $J\in\mathbb{R}^{M\times M}_{+}$ is a symmetric block matrix with the blocks on the diagonal equal to zero matrices, i.e. $J^{(ll)}=0\in\mathbb{R}^{M_{l}\times M_{l}}$ and $J^{(lr)}=(J^{(rl)})^{\mathrm{T}}$ for all $l,r=1,\dots,L$. The vector $\mathbf{h}$ and matrix $J$ contain the selection parameters, expressing, respectively, the single locus selection and the selective interaction among pairs of loci. Note that $\nabla V(\mathbf{x})=\mathbf{h}+J\mathbf{x}$, since the matrix $J$ is symmetric.

Let $g(\mathbf{x})=D(\mathbf{x})\nabla V(\mathbf{x})$. Then, the components of $g(\mathbf{x})$ are

[TABLE]

The drift function $\mu$ models the mutations. It is assumed that mutations occur independently at each locus, in particular, at the $l^{th}$ locus the mutation rate is $\frac{\theta_{l}}{2}$ and the probability matrix of mutations is $P^{(l)}=(P_{ij}^{(l)})_{i,j=1\dots,M_{l}}$. The transition rates of mutation from type $i$ to type $j$ at locus $l$ are thus $u_{ij}^{(l)}=\frac{\theta_{l}}{2}P_{ij}^{(l)}$. As in the standard Wright-Fisher model with parent dependent mutations, the components of the drift function are defined by

[TABLE]

Finally, the components of the diagonal block $D^{(l)}(\textbf{x}^{(l)})$ of the diffusion matrix $D(\textbf{x})$ are defined by

[TABLE]

which is characteristic for Wright-Fisher processes.

Having defined $\mu,D,$ and $V$, a compact definition of the coupled Wright-Fisher diffusion can be given, in terms of its infinitesimal generator. The coupled Wright-Fisher diffusion $\{\mathbf{X}(t)\}_{t\geq 0}$ is a $M$-dimensional diffusion process on the state space

[TABLE]

with generator

[TABLE]

where $\mu$, $g$ and $d$ are given by (2.2), (2.1) and (2.3), respectively. The generator $\mathcal{L}$ is defined on the domain $C^{2}(\mathcal{S})$.

Before proceeding with the derivation of the dual process, the stationary density (1.2) is considered. For existence of a density, the coupled Wright-Fisher diffusion must be restricted to the state space

[TABLE]

where

[TABLE]

corresponds to

[TABLE]

If there are no interactions among loci, the coupled Wright-Fisher diffusion consists of $L$ independent Wright-Fisher diffusions and the stationary density is well known when the mutations are parent independent. Wright himself proved that the stationary distribution of a single-locus, multi-type Wright-Fisher diffusion with parent independent mutations is Dirichlet [11]. Therefore, the stationary density of independent Wright-Fisher diffusions is the product of Dirichlet densities. More precisely, let

[TABLE]

$\pi(\bar{\textbf{x}})$ is the non-normalized stationary density of a coupled Wright-Fisher diffusion with no interaction among loci. In the presence of interaction, Aurell et al. [2] prove that there is an additional exponential factor in the stationary density, that is

[TABLE]

with $V$ defined on the restricted space by naturally defining the missing frequencies as one minus the sum of the other frequencies at the same locus. It is assumed that the mutations are parent independent, i.e. $u_{ij}^{(l)}=u_{j}^{(l)}$ $\forall l$. The form of the stationary density is explicit up to a normalising constant. In general, it is difficult to compute the normalising constant $Z$ explicitly, but under additional assumptions it can be computed numerically, as demonstrated in Section 4 and 6.

3. Outline of the derivation of a dual process

To derive a process that is dual to the coupled Wright-Fisher diffusion, a generator approach will be used as in [4], where the authors find a dual process for a multi-locus Wright-Fisher diffusion with recombination. In this section the method will be explained, in general terms.

Let $\mathcal{L}$ be the generator of the diffusion process (2.4) and $\mathcal{L}^{D}$ be the unknown generator of a dual process. Suppose that the following relationship holds

[TABLE]

for some duality function $F$ that needs to be determined. By using the relationship above the transition rates of a dual process can be identified from its generator. To pursue this approach, it is necessary to compute the left hand side of (3.1) by applying the generator $\mathcal{L}$ to the dual function $F$, considered as a function of x, and rewrite it on form

[TABLE]

for some non-negative coefficients $q(\textbf{n},\hat{\textbf{n}})$, $\hat{\textbf{n}}\in\mathbb{N}^{M}$, $\hat{\textbf{n}}\neq\textbf{n}$. In the light of the duality relationship, expression (3.2) can be interpreted as the generator $\mathcal{L}^{D}$ applied to the dual function $F$, considered as a function of n. Consequently, the dual process obtained this way is a pure jump process on the discrete space $\mathbb{N}^{M}$ with transition rate matrix $Q=(q(\cdot,\cdot))$, the off-diagonal elements being the non-negative coefficients in (3.2) and the diagonal elements being chosen so that the rows sums are equal to [math]. The alleged duality relationship is validated once the transition rates and the proper duality function are determined.

Consider the following proposal for the duality function, $F$. The inspiration for the proposal comes from the duality function for the single locus Wright-Fisher diffusion with mutations. It can be generalised to the multi-loci setting by taking

[TABLE]

for some function $k:\mathbb{N}^{M}\rightarrow\mathbb{R}$ that is determined as follows.

Let $\mathbf{\tilde{X}}$ be distributed according to the stationary distribution of the diffusion process $\{\mathbf{X}(t)\}_{t\geq 0}$, when such a density exists. Then $\mathbb{E}\left[\mathcal{L}F(\mathbf{\tilde{X}},\textbf{n})\right]=0$. Therefore, by taking expectation under the stationary distribution in (3.2), it follows that

[TABLE]

which implies that $\mathbb{E}\left[F(\mathbf{\tilde{X}},\cdot)\right]$ must be constant, the constant can be taken to be equal to $1$, and consequently,

[TABLE]

To find the transition rates of the dual process, it remains to obtain an expression of the form (3.2). In fact, it is sufficient to obtain an expression of the form

[TABLE]

with the requirement that $q(n,\hat{n})$ is positive for $\hat{n}\neq n$ (it will be soon clear that $q(\textbf{n},\textbf{n})$ is thus negative). In fact, once (3.5) is obtained, it is possible to derive expression (3.2) from it as follows. Rewriting (3.5) yields

[TABLE]

Keeping in mind that $\mathbb{E}\left[\mathcal{L}F(\mathbf{\tilde{X}},\textbf{n})\right]=0$ and that $k$ is chosen such that $\mathbb{E}\left[F(\mathbf{\tilde{X}},\mathbf{n})\right]$ is constant, one can apply expectation with respect to the stationary distribution to get

[TABLE]

Therefore (3.6) implies (3.2) and it remains to write $\mathcal{L}F$ as in (3.5) by finding the positive coefficients $q(n,\hat{n})$. Furthermore, (3.7) can be used to find a recursion formula for the function $k$, as demonstrated in the following section.

Throughout the rest of the paper, the emphasis will be on obtaining an expression of the type (3.5). This approach is first illustrated in a simpler case (single locus), in order to lighten the formulas and highlight the ideas, and is subsequently used in the general case of coupled Wright-Fisher diffusions. The simpler case turns out to be closely related to a well known model: the ancestral selection graph.

4. The ancestral selection graph

When only one locus is considered, the coupled Wright-Fisher diffusion is simply a Wright-Fisher diffusion with selection. Let $L=1$ , $M_{1}=2$ and assume that mutations are parent independent, i.e. $u_{ij}=u_{j}\;\text{for }i=1,2$. Then the matrix of pairwise selection parameters is the zero matrix and the quadratic function $V$ becomes linear

[TABLE]

Let $j(i)$ be the index opposite to $i$,

[TABLE]

Then, the drift terms can be written as follows

[TABLE]

The diffusion process solving (1.1) under the assumptions above is a two-types Wright-Fisher diffusion with selection and parent independent mutations. It is known that the genealogical process corresponding to this type of Wright-Fisher diffusion is embedded in a graph with coalescing and branching structure, the ancestral selection graph, studied by Krone and Neuhauser [7, 8]. A dual process for the same type of diffusion is derived using a generator approach also by Etheridge and Griffiths in [3], where the authors derive a dual process for the finite population size Moran model and use it to find the limiting transition rates of the dual process for the diffusion. Note that here is assumed, as in [3], that the types of individuals in the sample n are known, thus mutations are included in the dual process rather than superimposed afterwards on it, as in [7].

Following the outline in Section 3, a dual process is derived as follows. By applying the generator $\mathcal{L}$ to the duality function $F$ in (3.3), rewriting the derivatives of $F$, and rearranging the terms yields

[TABLE]

where $n=n_{1}+n_{2}$. To obtain an expression of the form (3.5) the expression in the last display can be rewritten as follows. First replace $x_{i}=1-x_{j(i)}$ to obtain positive coefficients for the selection terms, then use the identities, for $i=1,2$,

[TABLE]

where $\mathbf{e_{i}}$, $i=1,2$, are the unit vectors in $\mathbb{N}^{2}$. Finally it yields,

[TABLE]

which is the desired expression. As demonstrated in Section 3 the transition rates of a dual process can be identified directly from this expression. Therefore a dual process for the Wright-Fisher diffusion considered in this section is a pure jump process on the state space $\mathbb{N}^{2}$, with transition rates as follows. The dual process, in state n,

• •

jumps to state $\mathbf{n-e}_{i}$ $i,l=1,2$ at rate

[TABLE]

[coalescence]

• •

jumps to state $\mathbf{n+e}_{j(i)}-\mathbf{e}_{i}$ $i,l=1,2$ at rate

[TABLE]

[mutation]

• •

jumps to state $\mathbf{n+e}_{i}$ $i,l=1,2$ at rate

[TABLE]

[selection]

As anticipated, the dual process just described corresponds to the limiting process in [3] which mirrors the block counting process of the ancestral selection graph. From the transition rates $q$, it is observed that three types of events are possible for the dual process: mutation, coalescence and branching. The first two appear also in the Kingman’s coalescent, while the latter is a virtual addition to the true genealogical process which is characteristic of the ancestral selection graph. Seen forward in time, it corresponds to the event that two potential parents are chosen and only the one carrying the advantageous allele reproduces. Backward in time, when a branching happens, the individual splits into two individuals: its true parent and its virtual (potential) parent.

To complete the identification of the transition rates, $q(\textbf{n},\textbf{n})$ is defined as the coefficient of $F(\textbf{x},\textbf{n})$ in (4.4),

[TABLE]

Furthermore, (3.7) ensures that the rows sums of the transition matrix equal zero and provides a recursion formula for $k$

[TABLE]

with boundary conditions $k(\mathbf{e}_{1})=\rho_{1}$, $k(\mathbf{e}_{2})=1-\rho_{1}$, where $\rho_{1}$ is the probability that the ultimate ancestor is of type $1$, i.e. $\rho_{1}=\int_{0}^{1}x_{1}p(x_{1})dx_{1}$. This type of recursion formula is usually too large to be solved directly, importance sampling techniques have been developed to estimate $k$, see [10].

In general it is not possible to find a closed-form expression for $k$ and thus for the transition rates. However, when the mutations are parent independent, as in this example, the stationary density is explicitly known up to a normalizing constant $Z$ and thus $k$ can be written as an integral with respect to the stationary density

[TABLE]

The integral above cannot be computed analytically but it is related to the confluent hypergeometric function of the first kind, the Kummer’s function, which can be efficiently computed numerically. The idea of using Kummer’s functions is taken from [2] and [7], where this type of functions have been used to find respectively a series representation for the normalising constant in an example and a representation for the probability $\rho_{1}$. Let ${}_{1}F_{1}$ be the confluent hypergeometric function, then, using its integral representation, it yields

[TABLE]

see [1] for a complete collection of definitions and properties of confluent hypergeometric functions.

5. A multi-locus dual process

In this section the derivation in the previous section is extended to the general multi-locus setting, $L\geq 1$ and $M_{l}\geq 2$, $l=1,\dots,L$. The diffusion process is given by (1.1), with generator (2.4), where $\mu$, $g$ and $d$ are given by (2.2), (2.1) and (2.3), respectively. Let $k$ be defined by (3.4).

Theorem 5.1**.**

The dual process for the coupled Wright-Fisher diffusion is a pure jump process on the state space $\mathbb{N}^{M}$ with the following transition rates. From the current state, n, the dual process jumps to

• •

$\textbf{n}-\mathbf{e_{i}^{(l)}},i=1,\dots M_{l},l=1,\dots,L,$* at rate*

[TABLE]

[coalescence]

• •

$\textbf{n}-\mathbf{e_{i}^{(l)}+e_{j}^{(l)}},i,j=1,\dots M_{l},l=1,\dots,L,$* at rate*

[TABLE]

[mutation]

• •

$\textbf{n}+\mathbf{e_{j}^{(l)}},j=1,\dots M_{l},l=1,\dots,L,$* at rate*

[TABLE]

[single-locus selection]

• •

$\textbf{n}+\mathbf{e_{j}^{(l)}+e_{h}^{(r)}},j=1\dots M_{l},h=1\dots,M_{r},l,r=1,\dots,L,l\neq r,$* at rate*

[TABLE]

[double-locus selection]

where $n^{(l)}=\sum_{i=1}^{M_{l}}n_{i}^{(l)}$.

Note that the mutation and coalescence jumps involve one locus at the time. The coalescence and mutation rates are similar to the typical transition rates of the Kingman’s coalescent process with mutations, the only difference being the function $k$ which, despite having the same structure, is based on a different stationary density and depends on all the loci, not only on the one where the jump takes place. The single-locus selection rate not only contains the single-locus selection parameters in a form that generalises the rates in Section 4, but it also contains the double-selection parameters to include the effect of the pairwise interaction on the single locus. The single-locus selection, as the coalescence and mutation, jumps involve one locus at the time. Finally, the double-locus selection rate reflects the particular structure of pairwise interactions of the coupled Wright-Fisher diffusion and it is, to the best of our knowledge, a novel type of transition rate appearing in genealogical processes related to Wright-Fisher diffusions. The double-locus selection jumps correspond to simultaneous branchings at two different loci. The explicit parts of the transition rates (not depending on the function $k$) have a very natural interpretation. As in the simpler case studied in Section 4, the basic principle is that weak types branch at an higher rate. The difference is that, while in the simpler case there are only two types, a viable type and a weaker type, here there are many types and many loci all influencing each others branching rates. To understand this behaviour in greater detail, some terms will be investigated more thoroughly. The term

[TABLE]

arises purely from the single-locus selection and contributes to the rate of adding a gene of type $j$ at locus $l$. It depends on the single-locus viability of the other allele types (all except type $j$) at locus $l$, the higher their viability, the higher the rate of adding type $j$, and of course it is also directly proportional to the number of genes at locus $l$, $n^{(l)}$.

The rate of adding a couple of genes of type $j$ at locus $l$ and of type $h$ at locus $r$ is related to the term

[TABLE]

It depends on the viability of the other couples of allele types (all except couple $j,h$), the higher their viability, the higher the rate of adding type $j$ and $h$ at locus $l$ and $r$, respectively. Again the rate is directly proportional to the number of genes at loci $l$ and $r$.

Although the interpretation of some parts of the transition rates is straightforward, the function $k$ remains implicit, even for the simpler Kingman’s coalescent process with parent dependent mutations. When the mutations are parent independent, the stationary density is known up to a normalising constant and $k$ can be expressed as an integral that sometimes can be computed numerically, as shown in Section 4 and as it will be shown in Section 6, where a a series representation of $k$ involving Kummer’s and Beta functions will be given. However, in the general case, we have to settle for a recursion formula for $k$.

Proof of Theorem 5.1

Following the outline in Section 3, a dual process is derived as follows.

By applying the generator $\mathcal{L}$ to the duality function $F$ in (3.3) each term in the expression for $\mathcal{L}F$ is treated separately. As in the previous section, the terms corresponding to mutation and diffusion can be easily rewritten in the required form. Summing the mutation terms over allele types at locus $l$ yields

[TABLE]

By using identity (4.2) at locus $l$ the mutation terms can be rewritten in the desired form

[TABLE]

For the diffusion terms, the diagonal and off-diagonal terms are written separately as

[TABLE]

[TABLE]

Summing the diffusion terms at locus $l$ and rearranging yields

[TABLE]

Now use identity (4.1) at locus $l$ to obtain

[TABLE]

Next, consider the interaction terms. Using the definition of $g^{(l)}$ and rewriting the derivatives of $F$ yields,

[TABLE]

Note that the first group of sums, $S_{1}$, contains the single-locus selection parameters while the second, $S_{2}$, contains the pairwise selection parameters. Each of them will be treated separately. The single-locus selection term can be rearranged into

[TABLE]

As in Section 4, the fact that the sum of the frequencies at each locus equals one is used. Since $x_{k}^{(l)}=1-\sum_{\begin{subarray}{c}j=1\\ j\neq k\end{subarray}}^{M_{l}}x_{j}^{(l)}$, the terms can be rearranged to obtain

[TABLE]

The second part of (5.3) can be expressed as

[TABLE]

This time the equality

[TABLE]

will be used. To see that (5.4) holds, the fact that the frequencies sum up to one at each locus is used multiple times, as follows,

[TABLE]

Applying (5.4) in the expression for $S_{2}$ and rearranging the terms, leads to

[TABLE]

Summing over $l$ and putting similar terms together yields

[TABLE]

Use the identities (4.3) at locus $l$ and

[TABLE]

for the mixed terms involving loci $l$ and $r$, in order to rewrite the selection terms in the desired form

[TABLE]

The terms corresponding to mutation (5.1), diffusion (5.2) and selection (5.5) are now written in form (3.5). It is finally possible to identify the transition rates of the dual process and this completes the proof.

6. Two loci, two alleles, with selection and parent independent mutation

In this section a particular example will be considered, where there are two loci, $L=2$, and two allele types at each locus, $M_{1}=M_{2}=2$. The pairwise interactions are represented by the matrix

[TABLE]

and there is no single-locus selection, $\mathbf{h}=0$. Furthermore, parent independent mutations are assumed.

In this special case, the function $k$, in (3.4), and consequently the transition rates of the dual process can be computed rather efficiently. The main difficulty in the computation is that the normalising constant of the stationary density (2.5) is unknown. It may be noted that computing the normalising constant and the function $k$ are closely related problems. In fact, by defining

[TABLE]

the normalising constant can be written as

[TABLE]

and the function $k$ as

[TABLE]

The integral $I$ cannot be computed analytically, but it is possible to find a series representation of it in terms of Beta and Kummer functions, which can be truncated to numerically evaluate the function $k$. The following formula is derived by a straightforward, albeit cumbersome, application of definitions and properties of Kummer functions

[TABLE]

where $B$ is the Beta function, ${}_{1}F_{1}$ is the Kummer function and $[a]_{n}=a(a+1)\cdots(a+n-1)$, for $n>0$, $[a]_{0}=1$.

As an illustration, the stationary density of independent Wright-Fisher diffusions, with $J_{1}=J_{2}=0$, is compared to the stationary density of the coupled Wright-Fisher diffusion, with $J_{1}=J_{2}=2$, in Figure 6.1. Both distributions have mutation rates $u_{1}^{(1)}=u_{2}^{(1)}=u_{1}^{(2)}=u_{2}^{(2)}=0.8$.