Wright-Fisher diffusion bridges

TL;DR

This paper explores the genealogical structure of Wright-Fisher diffusion bridges, providing new interpretations, exact simulation algorithms, and connections to urn models, especially focusing on allele frequency trajectories with and without selection.

Contribution

It introduces a novel interpretation of coalescent genealogy in Wright-Fisher bridges, develops an exact simulation algorithm, and links the process to branching Polya urn models.

Findings

Genealogical structure involves two-direction coalescence from time t

Exact simulation algorithm for neutral Wright-Fisher bridges

Connection between Wright-Fisher bridges and Polya urn models

Abstract

{\bf Abstract} The trajectory of the frequency of an allele which begins at $x$ at time $0$ and is known to have frequency $z$ at time $T$ can be modelled by the bridge process of the Wright-Fisher diffusion. Bridges when $x=z=0$ are particularly interesting because they model the trajectory of the frequency of an allele which appears at a time, then is lost by random drift or mutation after a time $T$. The coalescent genealogy back in time of a population in a neutral Wright-Fisher diffusion process is well understood. In this paper we obtain a new interpretation of the coalescent genealogy of the population in a bridge from a time $t\in (0,T)$. In a bridge with allele frequencies of 0 at times 0 and $T$ the coalescence structure is that the population coalesces in two directions from $t$ to $0$ and $t$ to $T$ such that there is just one lineage of the allele under consideration at times $0$ and $T$. The genealogy in Wright-Fisher diffusion bridges with selection is more complex than in the neutral model, but still with the property of the population branching and coalescing in two directions from time $t\in (0,T)$. The density of the frequency of an allele at time $t$ is expressed in a way that shows coalescence in the two directions. A new algorithm for exact simulation of a neutral Wright-Fisher bridge is derived. This follows from knowing the density of the frequency in a bridge and exact simulation from the Wright-Fisher diffusion. The genealogy of the neutral Wright-Fisher bridge is also modelled by branching P\'olya urns, extending a representation in a Wright-Fisher diffusion. This is a new very interesting representation that relates Wright-Fisher bridges to classical urn models in a Bayesian setting.

Equations

$${\cal L}=\frac{1}{2}x(1-x)\frac{\partial^{2}}{\partial x^{2}}+\frac{1}{2}\Bigg{(}x(1-x)\gamma+(\theta_{1}-\theta x)\Bigg{)}\frac{\partial}{\partial x}.$$

$$f_{\theta_{1},\theta_{2}}(x,y;t)={\cal B}_{\theta_{1},\theta_{2}}(y)\Big{\{}1+\sum_{n=1}^{\infty}e^{-n(n+\theta-1)t/2}{\widehat{P}}^{(\theta_{1},\theta_{2})}_{n}(x){\widehat{P}}^{(\theta_{1},\theta_{2})}_{n}(y)\Big{\}},$$

$${\mathcal{B}}_{{\theta }_1,{\theta }_2}(y)=B({\theta }_1,{\theta }_2{)}^{-1}{y}^{{\theta }_1-1}(1-y{)}^{{\theta }_2-1},0<y<1$$

$$\widetilde{P}_{n}^{(a,b)}(r)=\frac{(a+1)_{(n)}}{n!}\>_{2}F_{1}\big{(}-n,a+b+n+1;a+1;(1-r)/2\big{)},$$

$$\displaystyle f_{\theta_{1},\theta_{2}}(x,y;t)=y^{\theta_{1}-1}(1-y)^{\theta_{2}-1}\Bigg{\{}\sum_{n=0}^{\infty}e^{-n(n+\theta-1)t/2}c_{n}(\theta_{1},\theta_{2})$$

$$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\times\widetilde{P}^{(\theta_{2}-1,\theta_{1}-1)}_{n}(r)\widetilde{P}^{(\theta_{2}-1,\theta_{1}-1)}_{n}(s)\Bigg{\}},$$

$$c_n({\theta }_1,{\theta }_2)=\frac{n!\Gamma (n+{\theta }_1+{\theta }_2-1)(2n+{\theta }_1+{\theta }_2-1)}{\Gamma (n+{\theta }_1)\Gamma (n+{\theta }_2)}.$$

$$f_{0,0}(x,y;t)=x(1-x)\sum _{i=1}^{\infty }{e}^{-i(i+1)t/2}(2i+1)i(i+1)R_{i-1}(r)R_{i-1}(s).$$

$$\underset{t\to \infty }{\lim }\frac{f_{0,0}(x,y;t)}{{\int }_0^1{f}_{0,0}(x,y;t)dy}$$

$$\text{constant}\times y\times y^{-1}(1-y{)}^{\theta -1}=\theta (1-y{)}^{\theta -1},0<y<1.$$

$$\text{constant}\times y(1-y)\times y^{-1}(1-y{)}^{-1}=1,0<y<1.$$

$$P(L^{\theta }(t)=l\mid L^{\theta }(0)=\infty )=q_l^{\theta }(t).$$

$$q_l^{\theta }(t)=\sum _{k=l}^{\infty }{\rho }_k^{\theta }(t)(-1{)}^{k-l}\frac{(2k+\theta -1)(l+\theta {)}_{(k-1)}}{l!(k-l)!},$$

$$\mathbb{E}\big{[}L^{\theta}(t)_{[k]}\big{]}=\sum_{l=k}^{\infty}\rho_{l}^{\theta}(t)(2l+\theta-1){l-1\choose k-1}(\theta+l)_{(k-1)}.$$

$${\phi }_l^{\theta }(t)=\frac{l(l+\theta -1)}{2}{q}^{\theta }(t),t>0.$$

$$f_{{\theta }_1,{\theta }_2}(x,y;t)=\sum _{l=0}^{\infty }{q}_l^{\theta }(t)\sum _{k=0}^l\left(\genfrac{}{}{0ex}{}{l}{k}\right)x^k(1-x{)}^{l-k}{\mathcal{B}}_{k+{\theta }_1,l-k+{\theta }_2}(y).$$

$$V\mathbf{U},(1-V)\{w_{(j)}\}$$

$$c(\mathbf{\theta }{)}^{-1}{e}^{\gamma x}B({\theta }_1,{\theta }_2{)}^{-1}{x}^{{\theta }_1-1}(1-x{)}^{{\theta }_2-1},0<x<1,$$

$$\frac{1}{2}x(1-x)w_{j}^{\prime\prime}+\frac{1}{2}\Bigg{(}\gamma x(1-x)+(\theta_{1}-\theta x)\Bigg{)}w_{j}^{\prime}+\lambda_{j}w_{j}=0.$$

$$f_{{\theta }_1,{\theta }_2,\gamma }(x,y,t)=\sum _{\mathbf{\alpha }\in {\mathbb{Z}}_{+}^2}{b}_{\mathbf{\alpha }}^{({\theta }_1,{\theta }_2)}(t,x)\pi [\mathbf{\alpha }+\mathbf{\theta }](y),$$

$$\pi [\mathbf{\theta }](y)=c(\mathbf{\theta }{)}^{-1}B({\theta }_1,{\theta }_2{)}^{-1}{e}^{\gamma y}{y}^{{\theta }_1-1}(1-y{)}^{{\theta }_2-1}.$$

$$q(\mathbf{\alpha },\mathbf{\alpha }-{\mathbf{e}}_i)$$

$$q(\mathbf{\alpha },\mathbf{\alpha }+{\mathbf{e}}_2)$$

$$q(\mathbf{\alpha },\mathbf{\alpha })$$

$$\frac{\partial }{\partial t}{b}_{\mathbf{\alpha }}^{({\theta }_1,{\theta }_2)}(t,x)=-q(\mathbf{\alpha },\mathbf{\alpha })b_{\mathbf{\alpha }}^{({\theta }_1,{\theta }_2)}(t,x)$$

$$\displaystyle~{}~{}+\sum_{i=1}^{2}\Bigg{(}q(\bm{\alpha}+\bm{e}_{i},\bm{\alpha})b^{(\theta_{1},\theta_{2})}_{\bm{\alpha}+\bm{e}_{i}}(t,x)+q(\bm{\alpha}-\bm{e}_{i},\bm{\alpha})b^{(\theta_{1},\theta_{2})}_{\bm{\alpha}-\bm{e}_{i}}(t,x)\Bigg{)},$$

$$\frac{\partial }{\partial t}{b}_{\mathbf{\alpha }}^{({\theta }_1,{\theta }_2)}(t,x)=\mathcal{L}{b}_{\mathbf{\alpha }}^{({\theta }_1,{\theta }_2)}(t,x),$$

$$b_{\mathbf{\alpha }}^{({\theta }_1,{\theta }_2)}(0,x)=\left(\genfrac{}{}{0ex}{}{|\mathbf{\alpha }|}{\mathbf{\alpha }}\right)x^{{\alpha }_1}(1-x{)}^{{\alpha }_2},$$

$$f_{x,z,[0,T]}(u,v;s,t)=\frac{f(u,v;t-s)f(v,z;T-t)}{f(u,z;T-s)},$$

$$f(x,y;t)\frac{h(y,t)}{h(x,0)}.$$

$$\underset{{\tau }^{\ast }\to \infty }{\lim }\frac{P_x(\tau >{\tau }^{\ast })}{h(x,{\tau }^{\ast })}=\text{constant}.$$

$$f(x,y;t)\frac{P_y(\tau >{\tau }^{\ast }-t)}{P_x(\tau >{\tau }^{\ast })}$$

Full text

Wright-Fisher diffusion bridges

Robert Griffiths Department of Statistics, University of Oxford; email: griff@stats.ox.ac.uk   

Paul A. Jenkins Department of Statistics and Department of Computer Science, University of Warwick; email: P.Jenkins@warwick.ac.uk   

Dario Spanò Department of Statistics, University of Warwick; email: D.Spanò@warwick.ac.uk

∗ Corresponding author. Version 3.1: TPB 2nd revision

Abstract The trajectory of the frequency of an allele which begins at $x$ at time [math] and is known to have frequency $z$ at time $T$ can be modelled by the bridge process of the Wright-Fisher diffusion. Bridges when $x=z=0$ are particularly interesting because they model the trajectory of the frequency of an allele which appears at a time, then is lost by random drift or mutation after a time $T$. The coalescent genealogy back in time of a population in a neutral Wright-Fisher diffusion process is well understood. In this paper we obtain a new interpretation of the coalescent genealogy of the population in a bridge from a time $t\in(0,T)$. In a bridge with allele frequencies of 0 at times 0 and $T$ the coalescence structure is that the population coalesces in two directions from $t$ to [math] and $t$ to $T$ such that there is just one lineage of the allele under consideration at times [math] and $T$. The genealogy in Wright-Fisher diffusion bridges with selection is more complex than in the neutral model, but still with the property of the population branching and coalescing in two directions from time $t\in(0,T)$. The density of the frequency of an allele at time $t$ is expressed in a way that shows coalescence in the two directions. A new algorithm for exact simulation of a neutral Wright-Fisher bridge is derived. This follows from knowing the density of the frequency in a bridge and exact simulation from the Wright-Fisher diffusion. The genealogy of the neutral Wright-Fisher bridge is also modelled by branching Pólya urns, extending a representation in a Wright-Fisher diffusion. This is a new very interesting representation that relates Wright-Fisher bridges to classical urn models in a Bayesian setting.

This paper is dedicated to the memory of Paul Joyce.

Keywords: Wright-Fisher diffusion bridges, Coalescent processes in Wright-Fisher diffusion bridges.

1 Introduction

The trajectory of the frequency of an allele which begins at $x$ at time [math] and is known to have frequency $z$ at time $T$ can be modelled by the bridge process of the Wright-Fisher diffusion. Henceforth we call this process the Wright-Fisher bridge. Developing a theory of Wright-Fisher bridges is interesting from a population genetics point of view, particularly when an allele frequency is conditioned to begin at [math] and end at [math]. This interest is because, informally, in populations where there is no mutation there is a recurrent process where individual genes begin with [math] frequency and are either lost at [math] frequency or fix at frequency $1$ at a subsequent time. This also occurs in the infinitely-many-alleles model where new mutants arise at frequency [math] and are lost at frequency [math] because of mutation away from any particular type. Wright-Fisher bridges are also interesting from a theoretical point of view as to how they fit into the general theory of bridges in diffusion processes. The Kingman coalescent and variations, (Griffiths, 1980; Kingman, 1982; Tavaré, 1984; Barbour et al., 2000), look backward in time in Wright-Fisher diffusion processes and are important in describing genealogy in populations. They are formally dual processes to the diffusion processes. Schraiber et. al. (2013) investigated bridges in a Wright-Fisher process in a neutral model and a model with selection. The neutral model genealogy plays a role in obtaining theoretical results in Schraiber et. al. (2013). The main emphasis in their paper was a clever rejection sampling simulation of the bridge beginning and ending with zero frequency when there is genic selection acting on the allele in the model. In the current paper the genealogy in neutral bridges and bridges where there is genic selection is investigated. Emphasis in this paper is on a thorough understanding and description of coalescent processes in the bridge. The coalescent interpretation in bridges is new. The neutral model and the model with selection are both reversible before fixation occurs which implies a bridge from time [math] to $T$ looks the same probabilistically in reverse time from $T$ to [math]. The genealogy in a Wright-Fisher diffusion with selection is more complex than in a neutral model. It is a branching coalescing random graph. We consider theory analogous to the neutral model bridge theory as much as possible, using the transition probabilities of Barbour et al. (2000) for the branching coalescing genealogy.

Section 2 describes known results about Wright-Fisher diffusion processes with and without selection that are necessary to understand bridges. We are particularly concerned with transition density expansions which will allow us to make calculations in a bridge probability density and to explain the coalescent genealogy in the bridge. A spectral expansion of the transition density by Kimura (1964) is well known, however an expansion as a mixture in terms of transition functions in a dual coalescent process is important and less well known.

Section 3 contains the main results on Wright-Fisher bridges. We begin by explaining the use of Doob $h$-transforms in the derivation of the density of an allele frequency in a bridge. In Section 3.1 neutral Wright-Fisher bridges when there is no mutation are considered in detail with a new interpretation of coalescent genealogy in the bridge. The genealogical form of the transition density of the Wright-Fisher diffusion in Section 2.1 is important in showing the coalescent genealogy in Wright-Fisher bridges. An interpretation of Theorem 2 in Section 3.1 shows that in a bridge with allele frequency beginning at [math] and ending at [math] the coalescence structure of the population from time $t\in(0,T)$ is that the population coalesces in two directions from $t$ to [math] and $t$ to $T$ such that there is just one lineage of the allele under consideration at times [math] and $T$. This is illustrated in Figure 1. The new case of the infinitely-many-alleles model where mutation is away from the allele type is also considered by finding a similar coalescent interpretation where lineages are now lost by coalescence or mutation. In Section 3.2 there is a subsection on $h$-transformations that condition on non-absorption of the allele frequency and give rise to new interesting identities for the coalescent transition functions and the Jacobi polynomials in the spectral expansions for the unconditional transition density. In Section 3.3 we consider the genealogy of the neutral Wright-Fisher bridge modelled by branching Pólya urns, extending a representation in a Wright-Fisher diffusion by Griffiths and Spanò (2010). This is a new very interesting representation that relates Wright-Fisher bridges to classical urn models in a Bayesian setting. In Section 3.4 Wright-Fisher bridges when there is selection in the model are considered. The genealogy is a branching coalescing graph more complex than the tree in the neutral model, but still with the property of the population branching and coalescing in two directions from time $t\in(0,T)$. This is illustrated in Figure 4. The form of the bridge distribution and coalescent interpretation are new. We consider the $h$-transform of the Wright-Fisher transition density with selection conditioned on non-absorption in the future. The transform is related to the first eigenvalue-eigenfunction pair. The Yaglom density is calculated from the first eigenfunction. This is a spheroidal wave function which we show how to calculate with Mathematica using Wolfram Alpha.

In the Appendix a method of exact simulation is developed for neutral Wright-Fisher bridges. The algorithm is based on the coalescent genealogy in the bridge and the alternative $h$-transform bridge density. It is not easy to do exact simulation so this new algorithm is interesting.

2 The Wright-Fisher diffusion process

In this section we describe known results about Wright-Fisher diffusion processes with and without selection that are necessary to understand bridges. Genealogical interpretations which give rise to transition density expansions are important in calculating bridge transition functions.

Let $A$ and $a$ be two types of a gene in a population of individuals. A Wright-Fisher diffusion process modelling the relative frequency of $A$ genes over time has a generator

[TABLE]

The population is subject to genic selection, whose sign and strength are described by $\gamma$; and mutations $a\to A$ and $A\to a$ which occur respectively at rates $\theta_{1}/2$, $\theta_{2}/2$, with notation $\theta=\theta_{1}+\theta_{2}$.

2.1 Neutral Wright-Fisher diffusion

Let $f_{\theta_{1},\theta_{2}}(x,y;t)$ be the transition density of the diffusion process with generator (1) when $\gamma=0$ so there is no selection in the model. There are two forms of a transition density expansion of $f_{\theta_{1},\theta_{2}}(x,y;t)$. The spectral expansion was derived by Kimura (1964) as

[TABLE]

for $\theta_{1},\theta_{2}\geq 0$, where

[TABLE]

is the Beta $(\theta_{1},\theta_{2})$ density and $\{{\widehat{P}}^{(\theta_{1}\theta_{2})}_{n}(y)\}_{n=0}^{\infty}$ are orthonormal polynomials on the Beta density derived from scaling the usual Jacobi polynomials $\{{\widetilde{P}}^{(a,b)}_{n}(r)\}_{n=0}^{\infty}$ which are well defined for $a,b\geq-1$ and are orthogonal on $(1-r)^{a}(1+r)^{b}$, $-1<r<1$. See, for example Ismail (2005). An expression for these polynomials is

[TABLE]

where ${}_{2}F_{1}$ is a hypergeometric function. In terms of these polynomials the expansion (2) is

[TABLE]

where $r=2x-1$, $s=2y-1$, and

[TABLE]

The expansion (3) holds for $\theta_{1},\theta_{2}\geq 0$ taking care with the starting summation index $n$. If $\theta_{1},\theta_{2}>0$ then the starting index is $n=0$, $\widetilde{P}_{0}^{(\theta_{2}-1,\theta_{1}-1)}\equiv 1$ and the first term is $c_{0}(\theta_{1},\theta_{2})=B(\theta_{1},\theta_{2})^{-1}$. If one of $\theta_{1},\theta_{2}$ is zero and the other non-zero then the summation begins from $n=1$, and if $\theta_{1}=\theta_{2}=0$ then the summation begins from $n=2$. In the last case the expansion found by Kimura was the form

[TABLE]

The polynomials $\{R_{i}(r)\}_{i=0}^{\infty}$ are scaled Jacobi polynomials $\{\widetilde{P}_{i}^{(1,1)}(r)\}_{i=0}^{\infty}$ orthogonal on $x(1-x),\>0<x<1$ such that $R_{i}(1)=1$. There is an identity between these polynomials and the Jacobi polynomials with index $(-1,-1)$ shown later in the paper as (33). If $\theta_{1},\theta_{2}>0$ then $\{X(t)\}_{t\geq 0}$ is reversible with respect to the stationary Beta $(\theta_{1},\theta_{2})$ distribution; or if at least one of $\theta_{1},\theta_{2}$ is zero $\{X(t)\}_{t\geq 0}$ is reversible, before being fixed, with respect to the speed measure.

The Yaglom density

[TABLE]

is the limit density of the diffusion process conditioned on the diffusion not being fixed at [math] or $1$ in $(0,t)$. The Yaglom density is straightforward to obtain from (3) or (4) when either or both of $\theta_{1},\theta_{2}$ are zero as being proportional to the first orthogonal polynomial times the speed measure density. If $\theta_{1}=0,\theta_{2}=\theta$ the Yaglom density is

[TABLE]

If $\theta_{1}=\theta_{2}=0$ the density is

[TABLE]

A second form of the transition density depends on the coalescent genealogy of the population. The Kingman coalescent is a death process which counts the number of lineages back in time in a coalescent tree beginning with a finite or infinite number of leaves. The death rates are ${k\choose 2}$, $k=2,3,\ldots$. Let $L^{\theta}(t)$ be the number of non-mutant lineages at time $t$ back in a coalescent tree, beginning with $L^{\theta}(0)$ leaves, which can be finite or infinity. Mutations occur according to a Poisson process of rate $\theta/2$ on the edges of the coalescent tree, so the death rate of non-mutant lineages is ${k\choose 2}+k\theta/2$, $k=1,2,\ldots.$. The last lineage is lost by mutation. If there is no mutation $L^{0}(t)$ counts the number of lineages in the Kingman coalescent. $\{L^{\theta}(t)\}_{t\geq 0}$ is a dual process to the Wright-Fisher diffusion which describes the behaviour of the population back in time. The duality is described in Ethier and Griffiths (1993); Griffiths and Spanò (2013). Denote the transition functions of $L^{\theta}(t)$ when $L^{\theta}(0)=\infty$ by

[TABLE]

Notation in this paper will be that for an integer $k\geq 0$ and real number $a$, $a_{[k]}=a(a-1)\cdots(a-k+1)$ and $a_{(k)}=a(a+1)\cdots(a+k-1)$. An explicit expression (Griffiths, 1980; Tavaré, 1984) is

[TABLE]

where $\rho_{k}^{\theta}(t)=\exp\{-\frac{1}{2}k(k+\theta-1)t\}$. The falling factorial moments of $L^{\theta}(t)$ are known from Griffiths (1980); Tavaré (1984) as

[TABLE]

Let $\{T_{l}\}_{l=1}^{\infty}$ be times between events when there are $l$ and $l-1$ non-mutant lineages. The density of $\sum_{k=l}^{\infty}T_{k}$ is

[TABLE]

The transition density of $X(t)$ can be expressed as a mixture of Beta densities

[TABLE]

The expansion (6) is a special case of a more general model in Ethier and Griffiths (1993). It appeared first implicitly in Griffiths (1980) and then written formally in Griffiths and Li (1983). The range of the summation index $k$ depends on whether $\theta_{1},\theta_{2}$ are zero or positive. If $\theta_{1},\theta_{2}>0$ then $0\leq k\leq l$; if $\theta_{1}=0,\theta_{2}>0$ then $1\leq k\leq l$; if $\theta_{1}>0,\theta_{2}=0$ then $0\leq k\leq l-1$; and if $\theta_{1}=0,\theta_{2}=0$ then $1\leq k\leq l-1$. Griffiths and Spanò (2010) show the algebraic connection between the two forms of the transition density (2) and (6). The transition density is improper if either of $\theta_{1}$ or $\theta_{2}$ are zero because the population may be fixed for one of the allele types by time $t$ when $X(t)=0$ or $X(t)=1$ with positive probability. Then $\int_{0}^{1}f_{\theta_{1},\theta_{2}}(x,y;t)dy$ is the probability that the population is not fixed by time $t$. This probability has a coalescent interpretation found by integrating in (6).

A genealogical interpretation of the Wright-Fisher population when $\theta_{1}=\theta_{2}=0$ is that the coalescent tree backward in time of the individuals beginning at $t$ has $l$ founder edges at the origin. The joint distribution of the $l$ family sizes of the founders at time $t$ is Dirichlet $(1,\ldots,1)$. The types of the families at $t$ are determined by the identity of the $l$ edges at time 0. The probability that an edge is of type $A$ is $x$, and $a$ is $1-x$. The second summation in (6) takes account that there are $k$ edges of type $A$ and $l-k$ of type $a$. The population is not fixed for one of the allele types if and only if $l\geq 1$ and $l-k\geq 1$. Grouping the family sizes into $k$ and $l-k$ in the Dirichlet distribution gives the distribution of the frequency of $A$ alleles as Beta $(k,l-k)$. If $\theta>0$ then families at time $t$ can be labelled by those from non-mutant edges at time zero and those from each mutation occurring on the coalescent tree. The joint distribution of the family sizes if there are $l$ non-mutant edges at time zero is

[TABLE]

where $\bm{U}$ has a $l$-dimensional Dirichlet$(1,1,\ldots,1)$ distribution; $V$ has a Beta $(l,\theta)$ distribution; $\{w_{(j)}\}$ has a Poisson-Dirichlet $(\theta)$ distribution and the three random variables are independent. This genealogical interpretation occurs in Griffiths (1980); Ethier and Griffiths (1993). If there are two types and $\theta_{1},\theta_{2}>0$ then choosing the type of the non-mutant edges as $A$ with probability $x$ and mutant types to be of type $A$ with probability $\theta_{1}/\theta$ gives the density (6). If $\theta_{1}=0,\theta_{2}>0$ the non-mutant edges type is chosen in the same way and all mutations are of type $a$ instead.

2.2 Wright-Fisher diffusion with selection

If $\gamma\neq 0$ selection is present in the model. When $\theta_{1},\theta_{2}>0$ the stationary density of $\{X(t)\}_{t\geq 0}$ is a weighted Beta density

[TABLE]

where c(\bm{\theta})=\mathbb{E}\big{[}e^{\gamma X}\big{]}, with expectation in a Beta $(\theta_{1},\theta_{2})$ distribution. $\{X(t)\}_{t\geq 0}$ is reversible with respect to the stationary distribution if $\theta_{1},\theta_{2}>0$ and reversible with respect to the speed measure if one or both of $\theta_{1},\theta_{2}$ are zero. The eigenvalues and eigenfunctions $\{\lambda_{j},w_{j}\}$ satisfy ${\cal L}w_{j}=\lambda_{j}w_{j}$, or

[TABLE]

The eigenfunctions are orthogonal on the weighted Beta distribution or speed measure, but are not polynomials. If $\theta_{1},\theta_{2}>0$ then $\lambda_{0}=1$ and $w_{0}=1$. There is no explicit expression for the eigenvalues. Song and Steinrücken (2012) expand the eigenfunctions in terms of Jacobi polynomials and then use numerical methods to find a computational solution of the transition density. The Yaglom limit density when either or both of $\theta_{1},\theta_{2}$ are zero is still proportional to the speed measure times the first eigenfunction, which does not have a simple form. If $\theta_{1}=\theta_{2}=0$, $w$ can be expressed in terms of a spheroidal wave function. This was known by Kimura (1955). We consider the expression in more detail later in the paper. If $\theta>0$ then the eigenfunctions cannot be expressed as spheroidal wave functions multiplied by a function not depending on the eigenfunction.

A genealogical form of the transition density of the diffusion process with generator (1), derived in Barbour et al. (2000), is

[TABLE]

where

[TABLE]

The expansion is analogous to (2) in the neutral model. If $\theta_{1}>0,\theta_{2}>0$ then $\pi[\bm{\theta}]$ is the stationary density of $\{X(t)\}_{t\geq 0}$. $b^{(\theta_{1},\theta_{2})}_{\bm{\alpha}}(t,x)$ are transition functions in a two-dimensional birth and death process $\bm{\alpha}(t)$ beginning with an entrance boundary at infinity where the relative frequencies of the types are $x,1-x$. If $\theta_{1}>0,\theta_{2}>0$ summation is over $\alpha_{1}\geq 0,\alpha_{2}\geq 0$; if $\theta_{1}=0,\theta_{2}>0$ summation is over $\alpha_{1}\geq 1,\alpha_{2}\geq 0$; if $\theta_{1}>0,\theta_{2}=0$ summation is over $\alpha_{1}\geq 0,\alpha_{2}\geq 1$; and if $\theta_{1}=0,\theta_{2}=0$ summation is over $\alpha_{1}\geq 1,\alpha_{2}\geq 1$. $\{\bm{\alpha}(t)\}_{t\geq 0}$ is a dual process to the diffusion process describing coalescent genealogy back in time. The transition rates of $\{\bm{\alpha}(t)\}_{t\geq 0}$ when $\gamma>0$, $\theta\geq 0$, are

[TABLE]

$b^{(\theta_{1},\theta_{2})}_{\bm{\alpha}}(t,x)$ satisfies the usual forward equations

[TABLE]

and the backwards equation

[TABLE]

where ${\cal L}$ is the generator (1). The boundary condition is

[TABLE]

the same as in a neutral model. Krone and Neuhauser (1997); Neuhauser and Krone (1997) introduced an ancestral selection graph, which is a branching coalescing graph describing the coalescent genealogy in a Wright-Fisher model with selection. In their graph the types of ancestral lineages are not fully determined until the type of the ultimate ancestor is reached. The ancestral tree is then read off from the graph. The process here is related, but different, in that the lineages are all typed in the graph from time $t$ back to time [math]. A helpful detailed description of the genealogy in a pre-limit Moran model and in our process is in Etheridge and Griffiths (2009).

The frequencies $y_{1}=y,y_{2}=1-y$ can be decomposed into independent ${\cal PD}(\theta_{1}+\alpha_{1})$ and ${\cal PD}(\theta_{2}+\alpha_{2})$ processes $\{w_{1i}\}$, $\{w_{2i}\}$ which are then weighted as $V\{w_{1i}\},(1-V)\{w_{2i}\}$ where $V$ is independent of the two processes and has a density $\pi[\bm{\theta}+\bm{\alpha}](v)$. This representation is then at the level of individual family sizes either beginning from the origin, or new mutations.

3 Wright-Fisher diffusion bridges

Let $\{X^{x,z,[0,T]}(t)\}_{0\leq t\leq T}$ be a Wright-Fisher diffusion process $\{X(t)\}_{t\geq 0}$ conditioned on $X(0)=x$ and $X(T)=z$. Notation is similar to that used in Schraiber et. al. (2013). We dispense with subscripts $\theta_{1},\theta_{2},\gamma$ on $f$ in the explanation of $h$-transforms that follows to ease notation. The transition density of $X(t)^{x,z,[0,T]}$ at $v$ given $X(s)^{x,z,[0,T]}=u$ is clearly

[TABLE]

where $f(x,y;t)$ is the transition density in the diffusion process with generator (1). There is a very elegant theory of Doob $h$-transforms and the invariance of a bridge distribution under a transform. The reader is referred to Rogers and Williams (2000) for an introduction to $h$-transforms. We give a very brief introduction and explain the relevance to bridges. In a $h$-transform the transition density is mapped to

[TABLE]

The function $h\geq 0$ is chosen so that (9) is the transition density of a stochastic process. Typical choices of $h$ provide an interpretation that the transformed process is conditioned on not being absorbed at boundary points in the future. In a Wright-Fisher diffusion without mutation [math] and $1$ are absorbing boundary points. In a model where absorption is certain let $\tau$ be the time when this occurs. Let $h(x)$ be the first eigenfunction of the generator satisfying ${\cal L}h(x)=-\lambda_{1}h(x)$, $h(x)\geq 0$, where $\lambda_{1}>0$ is the largest eigenvalue and set $h(x,t)=h(x)e^{-\lambda_{1}t}$. Then

[TABLE]

The transition functions (9) can be thought of as asymptotic transition functions of a process conditioned not to be lost or fixed. Let $\tau^{*}>t$, then the transition functions of a process conditioned on $\tau>\tau^{*}$ are

[TABLE]

which converges to (9) as $\tau^{*}\to\infty$. The generator of a process with infintesimal variance $\sigma^{2}(x)$ and mean $\mu(x)$ which is $h$-transformed by a $h$ not depending on time is

[TABLE]

and the stationary distribution is

[TABLE]

A more formal description of conditioning in general diffusion processes is the first theorem in Pinsky (1985).

The probability distribution of bridges is invariant under a $h$-transform. This is important for us in finding alternative forms for the density of the frequency of an allele in a bridge, particularly when the frequency is [math] at times [math] and $T$. As in Schraiber et. al. (2013) consider the joint density of the $h$-transformed bridge process at times $0<t_{1}<\cdots<t_{n}<T$. The density is

[TABLE]

This shows the invariance of the bridge process under $h$-transforms not depending on time. In a similar way the bridge density is invariant under a $h$-transform depending on time $h(x,t)=h(x)e^{\lambda t}$, where $\lambda$ is a constant. Reversibility of $\{X(t)\}_{t\geq 0}$ implies that the probabilistic behaviour of the bridge process should look the same from [math] to $T$ as it does backwards in time from $T$ to [math]. Let $\pi(x)$ be the density of the reversing measure so that $\pi(x)f(x,y;t)=\pi(y)f(y,x;t)$. Now consider the $h$-transform with $h(x)=\pi(x)^{-1}$. The transform reverses the path with (10) evaluating to

[TABLE]

which is the path density back in time from $T$ to [math] of a bridge from $y$ to $x$. The next theorem is general when a process is reversible. It appears for a general process in Fitzsimmons et. al. (1992) and in Schraiber et. al. (2013) for the Wright-Fisher diffusion.

Theorem 1**.**

The density of $X^{x,z,[0,T]}(t)$ and the reverse form is

[TABLE]

3.1 Neutral Wright-Fisher bridges

In this section the coalescent genealogy of bridges is considered. We are particularly interested in $\{X^{0,0,[0,T]}(t)\}_{0\leq t\leq T}$ bridges. Bridges where there is no mutation are studied, as well as bridges where mutation is away from the allele of interest at a rate $\theta$, as in the infinitely-many-alleles model. An approach is to let $x,z\to 0$ in the general bridge. Often beginning a diffusion process at [math] is analytically difficult, but using reversibility and a coalescent genealogy approach the treatment is relatively straightforward. The transition density as $x\to 0$ is

[TABLE]

The term of order $x$ in the first line of (12) is the joint density of $X(t)$ and the probability that there is exactly one coalescent lineage from time $t$ back to time zero that is of type $A$, and at least one other lineage of type $a$. Denote

[TABLE]

An interpretation is that $g(y;t)$ is the probability that there are at least two founder lineages with exactly one lineage of the same type as the allele with frequency $y$. This is because given $l$ founders at time $t$ back the family size at $t$ of a given lineage has density $(l-1)(1-y)^{l-2},\>0<y<1$ and there are $l$ choices for such a lineage.

Note the identity, from (6), that $f_{0,0}(y,0;t)=y(1-y)g(y;t)$. This is a consequence of reversibility with respect to the speed measure density $y^{-1}(1-y)^{-1}$ because

[TABLE]

Now dividing the numerator and denominator of the second line in (11) by $x$,

[TABLE]

Set $z=0$ in (14) to obtain that

[TABLE]

where

[TABLE]

Substituting in (5) with $\theta=0$ and $k=2$,

[TABLE]

Collecting results and substituting in (11) gives the next theorem.

Theorem 2**.**

(Schraiber et. al., 2013) The density of $X^{0,0,[0,T]}(t)$, $0\leq t\leq T$ in a neutral Wright-Fisher diffusion bridge with $\theta_{1}=\theta_{2}=0$ is

[TABLE]

There is a clear reversibility in the density (16) from both ends of the bridge.

A diagram of the genealogy of a (0,0) bridge is illustrated in Figure 1.

Figure 1: Genealogy in a neutral (0,0) bridge

x$$1-x$$z$$1-z$$y$$1-y[math]t$$T

There must be exactly 1 lineage from [math] to $t$ and from $T$ to $t$ in reverse time with leaves in the frequency $y$, conditional on $y$ being reached as $x,z\to 0$.

The description of the coalescent genealogy in the bridge is new in this paper and important. Informally the bridge density (16) times $dy$ is proportional to the proportion of sample paths in $(y,y+dy)$ at a given time where the allele of frequency $y$ has a single ancestor in both directions at times $t$ and $T-t$, divided by the proportion of sample paths in $(y,y+dy)$, which is the speed measure $y^{-1}(1-y)^{-1}$.

It is of interest to consider the coalescent genealogy at time points $s<t<T$ and calculate the number of lineages $b$, at $T$, which are non-mutant back to $t$, but not as far back as $s$; the number of non-mutant lineages $a$ back to time $s$; and the number of distinct non-mutant lineages $c$ back to time $s$. To begin condition on the number of non-mutant lineages from $s$ to $t$ $L^{\theta}(t-s)=l$. Then, from the representation at time $t$, considering the roots of the lineages from $T$ back to $t$, it must be that $a$ fall on frequencies from $V\bm{U}$ and $b$ fall on the frequencies $(1-V)\{w_{(j)}\}$. Therefore the conditional probability is

[TABLE]

The number of distinct non-mutant lineages $c$ back to $s$ is the number of distinct frequencies $V\bm{U}$ which are hit by the $a$ root lineages. The required probability is therefore

[TABLE]

Now

[TABLE]

so the sum (17) is equal to

[TABLE]

The marginal distribution of $c$ from (18) is, setting $m=a+b$,

[TABLE]

The inner sum of (19) is equal to

[TABLE]

which is the coefficient of $u^{c-1}$ in

[TABLE]

where $\xi$ has a Beta $(1,\theta)$ distribution. Now (21) evaluates to

[TABLE]

showing that (21) is equal to

[TABLE]

Finally

[TABLE]

From (22) an identity also follows that

[TABLE]

We now consider a bridge from 0 to 0 frequency when there is mutation away from a gene of type $A$. The importance of this bridge is that it models a new mutant in the infinitely-many-alleles model which arises at time 0, then is lost at time $T$ with mutation away from the new mutant to other types at rate $\theta/2$. The generator for the frequency of a given allele in a neutral model is

[TABLE]

The speed measure density is $y^{-1}(1-y)^{\theta-1}$ and the process is reversible with respect to this measure before absorption. Theorem 1 holds for this model. Again care is needed for a 0,0 bridge. As $x\to 0$

[TABLE]

and

[TABLE]

Denote

[TABLE]

A similar calculation to that in Theorem 2, substituting in the second line of (11) gives the bridge density at $t$.

Theorem 3**.**

The density of $X^{0,0,[0,T]}(t)$, $0\leq t\leq T$ in a neutral Wright-Fisher diffusion bridge with mutation away from $A$ at rate $\theta/2$ is

[TABLE]

The distribution of the frequency in the bridge in Theorem 3 is the same from both ends as it is when there is no mutation. There is a similar diagram and interpretation to that in Theorem 2, Figure 1 for the genealogy, except that only non-mutant lines from the frequency $y$ are considered and there can only be single non-mutant founder lineages from each of the frequencies at [math] and $T$ which have families in the $y$ frequency. Theorem 2 and the coalescent interpretation are new.

The mean frequency in a bridge is found by multiplying (24) by $y$ and integrating over $(0,1)$.

Corollary 1**.**

[TABLE]

An algorithm to simulate exactly from the density in a neutral Wright-Fisher bridge is given in the Appendix.

3.2 $h$-transform identities for bridges

There are useful and interesting identities for bridge densities that follow from $h$-transforms. The transforms relate neutral Wright-Fisher transition functions when one or both of $\theta_{1},\theta_{2}$ are zero to Wright-Fisher transition functions with recurrent mutation when $\theta_{1},\theta_{2}>0$. The bridge distribution is invariant under the $h$-transform. Bridges beginning at 0 are easier to deal with after the transformation when beginning at zero, because the Wright-Fisher transition function is then well defined.

Suppose $\theta_{1}=0,\theta_{2}=0$, then an identity is

[TABLE]

To show the identity (25) consider the following calculation with the backward generator.

[TABLE]

This shows that the generator in the $h$-transformed process is

[TABLE]

which corresponds to a Wright-Fisher diffusion with $\theta_{1}=\theta_{2}=2$. The bridge density (11) is invariant under the $h$-transform $h(x)=x(1-x)$ and also under $h(x,t)=e^{t}x(1-x)$. This choice of $h(x,t)$ is made because $x(1-x)$ is the first eigenfunction and the first eigenvalue is 1.

There is a $h$-transform identity similar to (25) when $\theta_{1}=0,\theta_{2}=\theta>0$ where

[TABLE]

A straightforward calculation is that

[TABLE]

This shows that the generator in the $h$-transformed process is

[TABLE]

which corresponds to a Wright-Fisher diffusion with $\theta_{1}=2$, $\theta_{2}=\theta$. The first eigenfunction in the process is $x$ and first eigenvalue $\theta/2$.

Corollary 2**.**

An alternative form for the density of $X^{x,z,[0,T]}(t)$ in a neutral model when $\theta_{1}=\theta_{2}=0$, for $0\leq x,z\leq 1$ is

[TABLE]

An alternative form when $\theta_{1}=0,\theta_{2}=\theta>0$ $0\leq x,z\leq 1$ is

[TABLE]

The densities on the right of these expressions are well defined at the boundaries when $x,z\in\{0,1\}$ because of non-zero recurrent mutation rates. Notice that the two alternative forms of the bridge density (27) and (28) are not continuous in $\theta$ at $\theta=0$.

The alternative forms become important in numerical calculations for the bridge density in a later section.

Corollary 3**.**

The limit density of $X^{x,z,[0,T]}(T/2+t)$, $\theta_{1}=0,\theta_{2}=0$ as $T\to\infty$ is independent of $t$ and is

[TABLE]

for $0\leq x,z\leq 1$.

Corollary 4**.**

The limit density of $X^{x,z,[0,T]}(T/2+t)$, $\theta>0$ as $T\to\infty$ is independent of $t$ and is

[TABLE]

for $0\leq x,z\leq 1$.

The density (27) in Corollary 2 gives an alternative expression to (16) by setting $x=z=0$, noting that terms on the right side are well defined in this case.

Other identities for the coalescent transition functions follow.

Letting $x\to 0$ in the identity (25), using (6)

[TABLE]

and equating coefficients of $y(1-y)^{l-1}$ shows the identity

[TABLE]

Another identity from the spectral expansion (3), and (25), equating coefficients of $e^{-\frac{1}{2}n(n+4-1)t}$, is that for $n\geq 0$,

[TABLE]

Letting $x\to 0$ in the identity (26) with $\theta_{1}=0,\theta_{2}=\theta>0$

[TABLE]

and equating coefficients of $y(1-y)^{l+\theta-2}$ shows the identity

[TABLE]

Also from (3) and (26)

[TABLE]

There is interest in the way the $h$-transforms are used to show the pairs of identities (32), (34) and (33) and (35). The first set is new, but probably not the second as so much is already known about Jacobi polynomials. The $h$-transform method of proof is new.

3.3 Wright-Fisher bridges and branching Pólya urns

A new representation of neutral Wright-Fisher bridges that come from Polya urns is now detailed. This has interest as a very classical probabilistic connection.

Genealogy in the Wright-Fisher diffusion process can be modelled by branching Pólya urns (Griffiths and Spanò, 2010). The joint density of $X=X(0)$ and $Y=X(t)$ when the marginal density of $X$ is Beta $(\theta_{1},\theta_{2})$, from the conditional density of $Y$ given $X$, (6), can be written as

[TABLE]

There is a Pólya urn interpretation of (36), explained in Griffiths and Spanò (2010, 2013). Consider two Pólya urns with an initial configuration of $\theta_{1}$ red and $\theta_{2}$ blue balls. The urns are identical until there are $l$ extra balls, then branch to two urns with probability $q^{\theta}_{l}(t)$. The composition of both urns at the moment of splitting has a Beta binomial distribution

[TABLE]

Draws continue independently in each urn after branching. The unconditional relative frequencies of red balls in the urns after an infinity of draws are Beta $(\theta_{1},\theta_{2})$ and the relative frequencies conditional on a configuration of $(k,l-k)$ are Beta $(\theta_{1}+k,\theta_{2}+l-k)$. In both cases the relative frequency distributions are the de Finetti measures for the sequences of draws. The transition density $f_{\theta_{1},\theta_{2}}(x,y;t)$ is the density of the relative frequency $Y$ in the second urn given the relative frequency $x$ in the first urn. Instead of arguing to obtain a joint distribution of $X$ and $Y$ it is possible to argue directly to obtain the conditional distribution of $Y$ given $X=x$. If $X=x$, then since the distribution of $X$ is the de Finetti measure for whether the draws are red or blue, the probability of a configuration of $(k,l-k)$ balls given the urns branch after $l$ draws and $X=x$ is ${l\choose k}x^{k}(1-x)^{l-k}$. the distribution of $Y$ given the $(k,l-k)$ configuration is then Beta $(\theta_{1}+l,\theta_{2}+l-k)$. Therefore the conditional distribution of $Y$ given $X=x$ is (6). The conditional density $f_{\theta_{1},\theta_{2}}(x,y;t)$ converges as $\theta_{1},\theta_{2}\to 0$ to

[TABLE]

Relating the limit to convergence in an urn model: if $l>1$ the conditional distribution of the urn configuration at $l$ given $0<x<1$ as $\theta_{1},\theta_{2}\to 0$ is

[TABLE]

as a consequence of the de Finetti representation. The joint probability that $k=0$ (and similarly for $k=l$) at trial $l$ and $0<x<1$ tends to zero as $\theta_{1},\theta_{2}\to 0$, even though it formally appears to be $(1-x)^{l}$.

We now give a new interpretation of Wright-Fisher bridges in terms of an urn model with three urns. In a neutral Wright-Fisher diffusion with mutation the density of the frequency at time $t$ in a bridge from $x$ to $z$ at time $t$ is (8), where $f_{\theta_{1},\theta_{2}}(x,y;t)$ is the density (2) or (6). An urn representation is found by considering the joint distribution of $X,Y,Z$ when $X$ and $Z$ have marginal Beta $(\theta_{1},\theta_{2})$ distributions. The joint distribution of the three random variables is then

[TABLE]

where $n=l+m$. An urn representation has three urns $U,V,W$ which share draws before branching after draw $n$ in $W$. The probability of branching after $n$ draws and that there are $l$, $m$ draws respectively in $U,V$ is $q_{l}(t)q_{m}(T-t)$. The $l,m$ draws are chosen at random from the $n$ draws. After branching the urn draws continue independently in $U,V,W$. For fixed $l,m$ the probability of configurations $(j,l-j)$ and $(k,m-k)$ of (red,blue) balls in urns $U,V$ and $(j+k,n-j-k)$ in urn $W$ is

[TABLE]

Conditional on $l,m,j,k$ the limit relative frequencies of red balls in $U,V,W$ have the joint density in the last line of (38). The unconditional joint density of the relative frequencies in urns $U,V,W$ is the full expression (38).

The configuration in the urns at branching, conditional on $X=x,Z=z$ has a probability distribution

[TABLE]

Then for $l,m\geq 1$ as $\theta_{1},\theta_{2}\to 0$ the distribution (39) converges to

[TABLE]

There is a similar limit if $\theta_{1}\to 0$ and $\theta_{2}>0$. The relative frequency of red balls in urn $W$ after an infinite number of draws when $\theta_{1},\theta_{2}>0$ has a Beta $(\theta_{1}+j+k,\theta_{2}+n-j-k)$ distribution so the density of $Y$ conditional on $X=x,Z=z$ is

[TABLE]

with a similar expression when $\theta_{1}=0,\theta_{2}=0$ where summation is over $l,m\geq 1$, $1\leq j\leq l-1$ and $1\leq k\leq m-1$. The density (41) is the same as the neutral bridge densities in Section 3.1 with the different cases of mutation. The limit distribution when $\theta_{1}=0,\theta_{2}=0$ and $x,z\to 0$ is now considered. The denominator of (40) converges to

[TABLE]

which simplifies as (15). The conditional probability distribution of the urn configurations (40) converges to zero unless $j=1,k=1$ when the limit is

[TABLE]

The conditional distribution of the relative frequency of red balls in urn $W$ is then

[TABLE]

3.4 Wright-Fisher bridges with genic selection

In this section we consider the genealogy of a Wright-Fisher bridge when there is selection in the model. This is a new approach different from that in Schraiber et. al. (2013). The genealogy of the Wright-Fisher diffusion with selection is more complex that in a neutral model and the transition functions for the coalescent genealogy do not have an explict form.

Theorem 1 holds when $\gamma\neq 0$, if $x\neq 0$, however care is again needed when beginning from $x=0$. For definiteness let $\gamma>0$, then the speed measure density when $\theta_{1}=0$ and $\theta_{2}=\theta\geq 0$ is

[TABLE]

The density is reversible before fixation and with $\bm{\theta}=(0,\theta)$,

[TABLE]

$\alpha_{i}\geq 1$ in the summation, $i=1,2$. In a neutral model

[TABLE]

An asymptotic form as $x\to 0$ is

[TABLE]

since $\pi[\bm{\alpha}+\bm{\theta}](x)=o(x^{\alpha_{1}-1})$ as $x\to 0$ for $\alpha_{1}>1$ and

[TABLE]

Note that

[TABLE]

Denote

[TABLE]

We have shown that as $x\to 0$

[TABLE]

and that

[TABLE]

The next theorem is an analogue of Theorem 2 which follows by substitution in the second line of the bridge density (11) and the asymptotic forms in this section for the transition density.

Theorem 4**.**

The density of $X^{0,0,[0,T]}(t)$, $0\leq t\leq T$ in a Wright-Fisher diffusion bridge when $\gamma>0$ with mutation away from $A$ at rate $\theta/2\geq 0$ is

[TABLE]

where prime denotes differentiation with respect to $y$.

The theorem shows that in a non-neutral $(0,0)$ Wright-Fisher bridge there is exactly one lineage from $t$ to [math] and from $t$ to $T$ and that there is a reversibility from either end of the bridge. This is a simlar interpretation to in the neutral bridge, Theorems 2 and 3. An illustration of the genealogy is in Figure 4. Theorem 4 and the coalescent interpretation are new.

Figure 4: Genealogy in a non-neutral (0,0) bridge

y$$1-y[math]t$$T

We are interested in finding the Yaglom density and the limit distribution of $X^{x,z,[0,T]}(T/2+t)$ as $T\to\infty$. These calculations depend on the first eigenfunction $w(y)$. The Yaglom density is proportional to $y^{-1}(1-y)^{\theta-1}e^{\gamma y}w(y)$. An analogy of the $h$-transform (25) in the non-neutral model is

[TABLE]

where $(\lambda,w)$ are the first eigenvalue, eigenfunction pair satisfying (7) when either or both of $\theta_{1},\theta_{2}$ are zero. A calculation similar to that in the neutral model shows that the generator associated with the transition functions $f_{w}$ when $\theta_{1}=0,\theta_{2}=\theta\geq 0$ is

[TABLE]

This is a generator of a Wright-Fisher diffusion with frequency dependent selection. The stationary distribution of the process with generator (43) is

[TABLE]

An alternative form for the density of $X^{x,z,[0,T]}$ is, for $0\leq t\leq T$

[TABLE]

The limit density as $T\to\infty$ of $X^{x,z,[0,T]}(T/2+t)$, from (44) is

[TABLE]

Computing $\psi(y)$ (when $\theta=0$) rests on finding a numerical solution of the eigenfunction equation

[TABLE]

for a maximal $\lambda$. There is a unique solution with $\lambda>0,w(x)\geq 0$. Clearly $w(0)=w(1)=0$ but the maximal $\lambda$ is unknown. Parameterizing $\lambda\equiv\lambda(\gamma)$, then $\lambda(0)=1$, the first eigenvalue in a neutral model. Kimura (1955) recognizes a connection with Spheroidal wave functions, tabulates $\lambda(\gamma)$, and expresses $w(x)$ as a series in Gegenbauer polynomials. The spheroidal wave equations are

[TABLE]

see NIST Handbook of Mathematical Functions (2010). If $c^{2}>0$ the equations are prolate or $c^{2}<0$ oblate. Solutions to (45) are the Spheroidal wave functions $P_{S_{n}^{m}}(z,c^{2})$ with eigenvalues $\mu_{n}^{m}(c^{2})$, $n=m,m+1,\ldots$. Letting

[TABLE]

then

[TABLE]

This shows that $e^{\gamma/2y}y^{-1/2}(1-y)^{-1/2}w(y)$ is proportional to $P_{S_{1}^{1}}(2y-1),-\gamma^{2}/16)$. Parameters are $m=1$, $c^{2}=-\gamma^{2}/16$, and $\mu-c^{2}=2\lambda$ and $z=1-2y$. Scaling $w(y)$ to be a probability distribution

[TABLE]

Spheroidal wave functions are implemented in Mathematica and it is straightforward to calculate and plot $w(y)$ via WolframAlpha on a web interface. The basic command for solutions to (45) is SpheroidalPS[n,m,c,z]. In our application $c^{2}\leq 0$ so $c$ is a multiple of the complex number $i$. For example if $\gamma=3$, then $c=0.75i$, and $n=m=1$,

Plot[(-0.166599)^(-1)exp(-1.5y)(y(1-y))^(1/2) SpheroidalPS[1,1,0.75i,2y-1], {y,0,1}]

is the appropriate command. The constant $-0.166599$ is the scaling factor returned by

Integrate[exp(-1.5x)(x(1-x))^(1/2) SpheroidalPS[1,1,0.75i,2x -1],{x,0,1}]

It is not necessary to find the first eigenvalue for the above commands, however SpheroidalEigenvalue[n,m,c] will return it. In this example Eigenvalue[1,1,0.75i] returns $2.44853$, so $\lambda=1.505515$. The full expansion of the transition function is

[TABLE]

This expansion is analogous to the neutral expansion (4). Mathematica also includes a function for the derivative of a Spheroidal wave function, SpheroidalPSPrime[n,m,c,z]. The drift coefficient in the $h$-transformed generator (43) when $\theta=0$ can thus be easily evaluated or plotted. The coefficient is easily seen to be

[TABLE]

Carrying through the example above, if $\gamma=3$, the command to plot the drift coefficient is

Plot[0.5-x +2.0x(1-x)SpheroidalPSPrime[1,1,0.75i,2x-1] /SpheroidalPS[1,1,0.75i,2*x-1],{x,0,1}]

In this particular example the the plot shows that the drift coefficient is very close to $1-2x$, which would be exact if $\gamma=0$ and $w(x)=x(1-x)$.

Appendix: Simulating the neutral Wright-Fisher bridge

It is possible to simulate exactly from the density of a Wright-Fisher bridge in a neutral model. We develop a new algorithm here when one or both of $x,z$ are zero. Jenkins and Spanò (2017) gave a sophisticated algorithm for this in the special case that both $\theta_{1},\theta_{2}>0$ and $x,z>0$. The identities in this paper provide a new direct method for other cases considered here:

$\theta_{1}=0$, $\theta_{2}=\theta>0$, and $x=z=0$. 2. 2.

$\theta_{1}=\theta_{2}=0$ and $x=z=0$. 3. 3.

$\theta_{1}=0$, $\theta_{2}=\theta>0$, and $x,z>0$. 4. 4.

$\theta_{1}=\theta_{2}=0$ and $x,z>0$.

1. $\theta_{1}=0$, $\theta_{2}=\theta>0$ and $x=z=0$

The density of $X^{0,0,[0,T]}(t)$ is given in Theorem 3. It can be written as

[TABLE]

where

[TABLE]

We recognise (46) as the density of an infinite mixture of Beta random variables, so that $(p_{l_{1},l_{2}})_{l_{1},l_{2}\in\mathbb{N}}$ is a probability mass function on $\mathbb{N}^{2}$ (for convenience, set $p_{l_{1},l_{2}}=0$ if $l_{1}=0$ or $l_{2}=0$). The $(l_{1},l_{2})$th component corresponds to a $\text{Beta}(2,l_{1}+l_{2}+\theta-2)$ variate. Therefore in order to simulate from $f^{0,0,[0,T]}(t)$, it suffices to simulate from the discrete distribution $(p_{l_{1},l_{2}})$. This is complicated by the infinite series representations for $q_{l_{1}}^{\theta}(t)$, $q_{l_{2}}^{\theta}(T-t)$, and $h^{\theta}(T)$, but a solution is possible via the series method (Devroye, 1986, Ch. IV.5); we summarise the strategy as follows.

Suppose $(p_{l})_{l\in\mathbb{N}}$ is a probability mass function whose masses may not be computable with finite resource but for which we have available a pair of sequences $(p^{-}_{l}(k))_{k\in\mathbb{N}}$, $(p^{+}_{l}(k))_{k\in\mathbb{N}}$ for each $l$ with the following properties:

$p^{-}_{l}(k)\leq p_{l}\leq p^{+}_{l}(k)$ for all $l,k\in\mathbb{N}$. 2. 2.

For each $l$, $p^{-}_{l}(k)\uparrow p_{l}$ as $k\to\infty$. 3. 3.

For each $l$, $p^{+}_{l}(k)\downarrow p_{l}$ as $k\to\infty$.

For $U\sim\text{Uniform}[0,1]$, standard inversion sampling implies that $\inf\{L\in\mathbb{N}:\sum_{l=0}^{L}p_{l}>U\}$ is distributed according to $(p_{l})$, but this is not computable, as we have noted. However, it is easily verified that if

[TABLE]

then

[TABLE]

is also distributed as $(p_{l})$, and this can be computed from finitely many terms in the double arrays $(p^{-}_{l}(k))$, $(p^{+}_{l}(k))$.

To employ this strategy for $(p_{l_{1},l_{2}})$ in (46) we must find monotonically converging upper and lower bounds on each of $q_{l_{1}}^{\theta}(t)$, $q_{l_{2}}^{\theta}(T-t)$, and $h^{\theta}(T)$. The first two expressions are covered by Jenkins and Spanò (2017, Proposition 1), who showed that

[TABLE]

provided $2k+l\geq C^{(t,\theta)}_{l}:=\inf\left\{i\geq l:\frac{\theta+i+l-1}{i-l+1}\frac{\theta+2i+1}{\theta+2i-1}e^{-(i+\theta/2)t}<1\right\}$, a constant beyond which convergence of these bounds becomes monotonic in $k$. Jenkins and Spanò (2017) assumed that $\theta_{1},\theta_{2}>0$, but their result follows without change when either or both of $\theta_{1}=0$, $\theta_{2}=0$.] It remains to find similar converging bounds on $h^{\theta}(T)$. Writing

[TABLE]

we find for $l\geq 2$ that

[TABLE]

with the final inequality holding provided $l>D^{(\theta,T)}:=\frac{\log 5}{T}-\frac{\theta}{2}$. It follows immediately that, for $k>\max(2,D^{(\theta,T)})$,

[TABLE]

and hence that

[TABLE]

serve as monotonically converging lower and upper bounds on $h^{\theta}(T)$ as $k\to\infty$. We summarise the above argument in the following result.

Theorem 5**.**

Let $\Sigma:\mathbb{N}\to\mathbb{N}^{2}$ be any bijective pairing function denoted by $\Sigma(l)=(l_{1},l_{2})$, and let

[TABLE]

If $k_{l}>\max\left((C_{l_{1}}^{(t)}-l_{1})/2,(C_{l_{2}}^{(T-t)}-l_{2})/2,D^{(T,\theta)},2\right)$ for each $l=0,1,\dots,L$, then $S^{-}_{\bm{k}}(L)$ and $S^{+}_{\bm{k}}(L)$ are respectively lower and upper bounds on the distribution function $\sum_{l=\Sigma(0)}^{\Sigma(L)}p_{l_{1},l_{2}}$ which are monotonically converging to this limit as ${\bm{k}}\to\infty$ (that is, as $k_{l}\to\infty$ for every component of ${\bm{k}}=(k_{0},k_{1},\dots,k_{L})$). Thus the following algorithm returns exact samples from the density of $X^{0,0,[0,T]}(t)$.

2. $\theta_{1}=\theta_{2}=0$ and $x=z=0$

This case is very similar to that of the previous subsection and so is omitted. (In fact, one can go through the previous subsection and replace $\theta$ with [math] throughout.) The density of the bridge when $\theta_{1}=\theta_{2}=0$ simplifies slightly (see Theorem 2) to

[TABLE]

3. $\theta_{1}=0$, $\theta_{2}=\theta>0$ and $x,z>0$

By Corollary 2, the density of $X^{x,z,[0,T]}(t)$ in this case is the same as for $\theta_{1}=2$, $\theta_{2}=\theta$, which is covered by Algorithm 4 of Jenkins and Spanò (2017).

4. $\theta_{1}=\theta_{2}=0$ and $x,z>0$

By Corollary 2, the density of $X^{x,z,[0,T]}(t)$ in this case is the same as for $\theta_{1}=\theta_{2}=2$, which is covered by Algorithm 4 of Jenkins and Spanò (2017).