The Common Ancestor Process for a Wright-Fisher Diffusion

Rates of molecular evolution along phylogenetic trees are influenced by mutation, selection and genetic drift. Provided that the branches of the tree correspond to lineages belonging to genetically isolated populations (e.g., multi-species phylogenies), the interplay between these three processes can be described by analyzing the process of substitutions to the common ancestor of each population. We characterize this process for a class of diffusion models from population genetics theory using the structured coalescent process introduced by Kaplan et al. (1988) and formalized in Barton et al. (2004). For two-allele models, this approach allows both the stationary distribution of the type of the common ancestor and the generator of the common ancestor process to be determined by solving a one-dimensional boundary value problem. In the case of a Wright-Fisher diffusion with genic selection, this solution can be found in closed form, and we show that our results complement those obtained by Fearnhead (2002) using the ancestral selection graph. We also observe that approximations which neglect recurrent mutation can significantly underestimate the exact substitution rates when selection is strong. Furthermore, although we are unable to find closed-form expressions for models with frequency-dependent selection, we can still solve the corresponding boundary value problem numerically and then use this solution to calculate the substitution rates to the common ancestor. We illustrate this approach by studying the effect of dominance on the common ancestor process in a diploid population. Finally, we show that the theory can be formally extended to diffusion models with more than two genetic backgrounds, but that it leads to systems of singular partial differential equations which we have been unable to solve.


Introduction
One of the key insights to emerge from population genetics theory is that the effectiveness of natural selection is reduced by random variation in individual survival and reproduction. Although the expected frequency of a mutation will either rise or fall according to its effect on fitness, evolution in finite populations also depends on numerous chance events which affect individual life histories in a manner independent of an individual's genotype. Collectively, these events give rise to a process of stochastic fluctuations in genotype frequencies known as genetic drift (Gillespie 2004). For example, a mutation which confers resistance to a lethal infection will still decline in frequency if, by chance, disproportionately many of the individuals carrying that mutation are killed in a severe storm. Moreover, if the mutation is initially carried by only a few individuals, then it may be lost altogether from the population following such a catastrophe. Because it is counterintuitive that populations may evolve to become less fit, there has been much interest in the consequences of stochasticity for other aspects of adaptive evolution, such as the origin of sex (Poon and Chao 2004; Barton and Otto 2005), genome composition (Lynch and Conery 2003), and speciation and extinction (Whitlock 2000;Gavrilets 2003).
Testing these theories requires quantifying genetic drift and selection in natural populations. Although selection and drift can sometimes be inferred from historical changes in the distribution of a trait (Lande 1976) or genotype frequencies (O'Hara 2005), population genetical processes are mainly investigated using sets of contemporaneously sampled DNA sequences. For our purposes, it is useful to distinguish two scenarios. On the one hand, sequences sampled from a single population will usually share a common history shaped by selection and drift, and must be analyzed using models which take that shared history into account. One approach is to reduce the data to a set of summary statistics whose distribution can be predicted using population genetical models (Sawyer and Hartl 1992;Akashi 1995;Bustamante et al. 2001). Alternatively, more powerful analyses can be designed by using coalescent models and Monte Carlo simulations to estimate the joint likelihood of the data and the unobserved genealogy under different assumptions about selection and drift (Stephens and Donnelly 2003;Coop and Griffiths 2004). In both cases, the selection coefficients estimated with these methods will reflect the combined effects of selection and genetic drift in the population from which the sample was collected.
In contrast, when the data consists of sequences sampled from different species, then the time elapsed since any of the ancestors last belonged to a common population may be so great that the genealogy of the sample is essentially unrelated to the population genetical processes of interest. In this case, the genealogy is usually inferred using purely phylogenetic methods, and evolutionary inferences are facilitated by making certain simplifying assumptions about the way in which natural selection influences the substitution process along branches of this tree, i.e., the process of mutations to the ancestral lineages of the members of the sample. It is usually assumed that the substitution process along each branch of the tree is a Markov process, and that substitutions by beneficial or deleterious mutations occur at rates which are either greater than or less than the neutral mutation rate (Yang 1996). While the first assumption is true only when evolution is neutral, i.e., mutations do not affect fitness, the latter assumption reflects the fact that mutations which either increase or decrease the likelihood of a lineage persisting into the future are likely to be over-or under-represented, respectively, on lineages which do in fact persist. For example, it is often possible to identify proteins which are under unusually strong selection simply by comparing the rates of substitutions which change the amino acid composition of the protein with those which do not (Nielsen and Yang 1998).
An important limitation of purely phylogenetic analyses of selection is that the relationship between the phylogenetic rate parameters and population genetical quantities is usually obscure. One exception is when less fit variants are in fact lethal, so that selection is fully efficient and certain substitutions are never observed in live individuals. Alternatively, if the mutation rates are small enough that each new mutation is either rapidly lost or fixed in the population, then under some circumstances the substitution rate can be approximated by the flux of mutations which go to fixation (Kimura 1964). This approach has been used by McVean and Vieira (2001) to estimate the strength of selection on so-called silent mutations (i.e, those which do not change amino acid sequences) in several Drosophila species.
The common ancestor process can be used to describe the relationship between phylogenetic substitution rates and population genetical processes when the preceding approximations do not hold. The common ancestor of a population is any individual which is ancestral to the entire population. For the models which will be studied in this paper, such an individual will be guaranteed to exist at some time sufficiently (but finitely) far into the past and will be unique at any time at which it does exist. Denoting the type of the common ancestor alive at time t by z t , we will define the substitution process to the common ancestor to be the stochastic process (z t : t ∈ R) and the common ancestor distribution to be the stationary distribution of z t . This process will be a good approximation to the substitution process along the branches of a phylogenetic tree provided that the time elapsed along each branch is large in comparison with the coalescent time scales of the populations containing the sampled individuals and their ancestors. In particular, the divergence between the sequences in the sample should be much greater than the polymorphism within the populations from which the sample was collected. As is customary in modeling molecular evolution (Zharkikh 1994), we will assume that these populations are at equilibrium and that evolutionary processes such as mutation and selection do not vary along ancestral lineages. Although common ancestor processes could also be defined for non-equilibrium and time-inhomogeneous models, characterization of such processes will be substantially more difficult than in the idealized cases considered here.
Common ancestor distributions were first described for supercritical multitype branching processes by Jagers (1989Jagers ( , 1992, who showed that the distribution of the type of an individual spawning a branching process which survives forever has a simple representation involving the leading left and right eigenvectors of the first moment generator of the branching process. Because such an individual gives rise to infinitely many lineages which survive forever, but which individually do not give rise to the entire future population, it is not meaningful to speak of the common ancestor process in this setting. Instead, we must study what Georgii and Baake (2003) call the retrospective process, which characterizes the substitution process along lineages which survive forever. This process was also first described by Jagers (1989Jagers ( , 1992, who showed it to be a stationary time-homogeneous Markov process having the common ancestor distribution as its stationary measure. Extensive results concerning the retrospective process and common ancestor distribution can be found in Georgii and Baake (2003) and Baake and Georgii (2007).
Much less is known about the common ancestor process for traditional population genetical models such as the Moran and Wright-Fisher processes in which the population size remains constant. For neutral models, the fact that the substitution process decouples from the genealogy of a sample can be used to deduce that the common ancestor process is simply the neutral mutation process and that the common ancestor distribution is the stationary measure of this process. That this also holds true in the diffusion limit can be shown using the look-down construction of Donnelly and Kurtz (1996), which provides a particle representation for the Wright-Fisher diffusion. The key idea behind this construction is to assign particles to levels and then introduce look-down events which differ from (and replace) the usual neutral two-particle birth-death events of the Moran model in the requirement that it is always the particle occupying the higher level which dies and is then replaced by an offspring of the particle occupying the lower level. In the absence of selection, the common ancestor is the particle occupying the lowest level, as this individual never dies and it can be shown that all particles occupying higher levels have ancestors which coalesce with this lowest level in finite time.
In contrast, when selection is incorporated into the look-down process, particles can jump to higher levels and the common ancestor is no longer confined to the lowest level (Donnelly and Kurtz 1999). Furthermore, because the effect of selection depends on the frequencies of the types segregating in the population, e.g., selection has no effect if the population is monomorphic, we do not expect the non-neutral common ancestor process to be a Markov process. However, the mathematical difficulties which this creates can be overcome with the same technique that is used to characterize the genealogical processes of such models, namely by enlarging the state space of the process of interest until we obtain a higher dimensional process which does satisfy the Markov property. One such enlargement is the ancestral selection graph of Krone and Neuhauser (1997), which augments the ancestral lineages of the genealogy with a random family of 'virtual' lineages which are allowed to both branch and coalesce backwards in time. Fearnhead (2002) uses a related process to identify the common ancestor process for the Wright-Fisher diffusion with genic selection. His treatment relies on the observation that when there is only a single ancestral lineage, certain classes of events can be omitted from the ancestral selection graph so that the accessible particle configurations consist of the common ancestor, which can be of either type, plus a random number of virtual particles, all of the less fit type. This allows the common ancestor process to be embedded within a relatively tractable bivariate Markov process (z t , n t ), where z t is the type of the common ancestor and n t is the number of virtual lineages.
In this article, we will use a different enlargement of the non-neutral coalescent. Our treatment relies on the structured coalescent introduced by Kaplan et al. (1988) and formalized by , which subdivides the population into groups of exchangeable individuals sharing the same genotype and records both the types of the lineages ancestral to a sample from the population and the past frequencies of those types. With this approach, the common ancestor process of a population segregating two alleles can be embedded within a bivariate process (z t , p t ), where p t is the frequency at time t of one of the two alleles. We will show that both the stationary distribution and the generator of this process can be expressed in terms of the solution to a simple boundary value problem (Eq. 9) which determines the distribution of the type of the common ancestor conditional on the frequency at which that type occurs within the population. In certain cases we can solve this problem exactly and obtain an analytical characterization of the common ancestor process. However, one advantage of the diffusion-theoretic approach described here is that even when we cannot write down an explicit solution, we can still solve the corresponding boundary problem numerically. This makes it possible to calculate the substitution rates to the common ancestor for a much more general set of population genetical models than can be dealt with using the ancestral selection graph, including models with frequency-dependent selection, which we illustrate in Section 5, as well as fluctuating selection and genetic hitchhiking which 812 will be described elsewhere.
The remainder of the article is structured as follows. In Section 2 we describe the class of diffusion processes to be studied and we briefly recall the construction of the structured coalescent in a fluctuating background as well as its restriction to a single ancestral lineage, which we call the structured retrospective process. Using calculations with generators, we describe the stationary distribution of the structured retrospective process and identify the common ancestor process by reversing the retrospective process with respect to this measure. We also give an alternative probabilistic representation for the conditional distribution of the type of the common ancestor, and in Section 3 we use this to derive asymptotic expressions for the substitution rates to the common ancestor when the mutation rates are vanishingly small. Sections 4 and 5 are concerned with applications of these methods to concrete examples, and we first consider the Wright-Fisher diffusion with genic (frequency-independent) selection. In this case we can write the density of the common ancestor distribution in closed form (Eq. 23), and we show that this quantity is related to the probability generating function of a distribution which arises in the graphical representation of Fearnhead (2002). Notably, these calculations also show that approximations which neglect recurrent mutation (e.g., the weak mutation limits) can underestimate the true substitution rates by an order of magnitude or more when selection is strong. In contrast, few explicit calculations are possible when we incorporate dominance into the model in Section 5, and we instead resort to numerically solving the associated boundary value problem to determine the substitution rates to the common ancestor. In the final section we show that some of these results can be formally extended to diffusion models with more than two genetic backgrounds, but that the usefulness of the theory is limited by the need to solve boundary value problems involving systems of singular PDE's.

Diffusions, coalescents and the common ancestor
We begin by recalling the structured coalescent process introduced by Kaplan et al. (1989) and more recently studied by  and . Consider a closed population, of constant size N , and let P and Q be two alleles which can occur at a particular locus. Suppose that the mutation rates from Q to P and from P to Q are µ 1 and µ 2 , respectively, where both rates are expressed in units of events per N generations. Suppose, in addition, that the relative fitnesses of P and Q are equal to 1 + σ(p)/N and 1, respectively, where p is the frequency of P . For technical reasons, we will assume that the selection coefficient σ : [0, 1] → ∞ is the restriction of a function which is smooth on a neighborhood of [0, 1], e.g., σ(p) could be a polynomial function of the frequency of P . If we let p t denote the frequency of P at time t and we measure time in units of N generations, then for sufficiently large N the time evolution of p t can be approximated by a Wright-Fisher diffusion with generator where φ ∈ C 2 ([0, 1]). If we instead consider a diploid population, then the time evolution of the frequency of P can be modeled by the same diffusion approximation if we replace N by 2N .
We note that because the drift and variance coefficients are smooth, with support contained in the interior of (0, 1) is a core for A. Furthermore, provided that both mutation rates µ 1 and µ 2 are positive, then the diffusion corresponding to (1) has a unique stationary measure π(dp) on [0, 1], with density (Shiga 1981, Theorem 3.1; Ewens 2004, Section 4.5), where C is a normalizing constant. Unless stated otherwise (i.e., when we consider weak mutation limits in Section 3), we will assume throughout this article that both mutation rates are positive.
Although the structured coalescent can be fully characterized for this diffusion model, for our purposes it will suffice to consider only the numbers of ancestral lineages of type P or Q, which we denoteñ 1 (t) andñ 2 (t), respectively. Here, and throughout the article, we will use the tilde, both on random variables and on generators, to indicate a stochastic process which is running from the present (usually the time of sampling) to the past. Then, as shown in , the generatorG of the structured coalescent process (ñ 1 (t),ñ 2 (t),p t ) can be written as where for each (n 1 , n 2 ) ∈ N ×N , we have φ(n 1 , n 2 , ·) ∈ C 2 ([0, 1]).  prove that a Markov process corresponding to this generator exists and is unique, and moreover that this process is the weak limit of a suitably rescaled sequence of Markov processes describing both the sample genealogy and the allele frequencies in a population of size N evolving according to a Moran model. One particularly convenient property of biallelic diffusion models is that the processp(t) governing the evolution of allele frequencies backwards in time in a stationary population has the same law as the original Wright-Fisher diffusion p(t) corresponding to the generator A. In fact, this property is shared by one-dimensional diffusions in general, which satisfy a detailed balance condition with respect to their stationary distributions (Nelson 1958). This will not be true (in general) of the multidimensional diffusion models considered in Section 6, where we will characterize the common ancestor process at a locus which can occur in more than two genetic backgrounds which can change either by mutation or by recombination.
Because we are only concerned with substitutions to single lineages, we need only consider sample configurations (n 1 , n 2 ) which are either (1, 0) or (0, 1), and so we can replace the trivariate process (ñ 1 (t),ñ 2 (t),p t ) with a bivariate process (z,p t ) taking values in the space E = ({1} × (0, 1]) ∪ ({2} × [0, 1)), wherez t = 1 if the lineage is of type P andz t = 2 if it is of type Q. We will refer to (z t ,p t ) as the structured retrospective process to emphasize the fact that it describes evolution backwards in time. (In contrast, Georgii and Baake (2003) define a retrospective process for a multitype branching process which runs forwards in time.) With this notation, the generator of the structured retrospective process can be written as for functions φ ∈ D(G) ≡ C 2 c (E) which are twice continuously differentiable on E and have compact support. For future reference we note that D(G) is dense in the spaceĈ(E) of continuous functions on E vanishing at infinity and that D(G) is an algebra. The key step in proving the existence and uniqueness of a Markov process corresponding to this generator is to show that the ancestral lineage is certain to jump away from a type before the frequency of that type vanishes, e.g., the ancestor will almost surely mutate from P to Q before the diffusionp t hits 0. This will guarantee that the jump terms appearing inG, which diverge at the boundaries of the state space, are in fact bounded along trajectories of the process over any finite time interval [0, T ]. That the jumps do happen in time is a consequence of Lemma 4.4 of , which we restate below as Lemma 2.1.
We also supply a new proof of this lemma to replace that given in , which contains two errors (Etheridge 2005). One is that the variance σ(W s ) appearing in the time change of the Wright-Fisher diffusion needs to be squared, so that the exponent α in the integral displayed in Eq. (16) of that paper is 2 rather than 1+ 1 2(1−2µ 2 ) . The second is that the divergence of this integral requires α ≥ 2 rather than α ≥ 1. Although this condition is (just barely) satisfied, we cannot deduce the divergence of the integral from the Engelbert-Schmidt 0-1 law (Karatzas and Shreve, 1991, Chapter 3, Proposition 6.27; see also Problem 1 of Ethier and Kurtz, 1986, Chapter 6) because this result applies to functionals of a Brownian path integrated for fixed periods of time rather than along sample paths which are stopped at a random time, as is the case in Eq. (16). Lemma 2.1. Let p t be the Wright-Fisher diffusion corresponding to the generator A shown in (1). Then, for any real number R < ∞, and observe that on this restricted set, Then, for 815 each k > p −1 0 , the stopped process is a continuous martingale with quadratic variation In particular, on the set {τ k < ∞}, we have which in turn implies that, for any R < ∞, the following three inequalities are satisfied on the set Now, because M ·∧τ k is a continuous, one-dimensional martingale, there is an enlargement Ω ′ of the probability space Ω on which the diffusion p t is defined and there is also a standard one-dimensional Brownian motion B t , defined on Ω ′ , such that [See Karatzas and Shreve (1991), Chapter 3, Theorem 4.6 and Problem 4.7.] Thus, in view of the conditions holding on Ω R,k , we obtain the following bound where C = 1 2 + ||b|| ∞ is independent of k. The first half of the proposition then follows from the fact that the probability on the right-hand side of the preceding inequality goes to 0 as k → ∞ with R fixed. The second half can be proved using a similar argument, with φ k (p) = − ln(1 − p) on [0, 1 − 1/k].
With Lemma 2.1 established, the next proposition is a special case of the existence and uniqueness results for structured coalescents proved in . Proposition 2.2. For any ν ∈ P(E), there exists a Markov process (z t ,p t ), which we call the structured retrospective process, which is the unique solution to the D E [0, ∞)-martingale problem for (G, ν).
Proof. Because the operatorG is a Feller generator when restricted to twice continuously differentiable functions on each of the sets show that a stopped version of the process exists on each of these sets and that this process is the unique solution of the corresponding stopped martingale problem. Then, using the Lemma 2.1 and noting that the diffusions p t andp t are identical in distribution, we can show that the sequence of hitting times of the boundaries of the sets E k is almost surely unbounded as k → ∞. Consequently, Theorem 4. Of course, as the name indicates, the process (z t ,p t ) describes the retrospective behavior of a lineage sampled at random from the population rather than forward-in-time evolution of the common ancestor of the entire population. However, because Kingman's coalescent comes down from infinity (Kingman 1982), we know that, with probability one, all extant lineages, including that ancestral to the sampled individual, will coalesce with the common ancestor within some finite time. That this is still true when we incorporate genetic structure into the coalescent is evident from the fact that the coalescent rates within a background are accelerated by the reciprocal of the frequency of that background; see Eq. (3). Furthermore, because lineages move between genetic backgrounds at rates which are bounded below by the (positive) mutation rates, lineages cannot be permanently trapped in different backgrounds.
These observations lead to the following strategy for identifying the common ancestor process in a stationary population. First, because the asymptotic properties of the retrospective process in the deep evolutionary past coincide with those of the common ancestor process itself, any stationary distribution ofG will also be a stationary distribution of the common ancestor process. Indeed, we will call this distribution (assuming uniqueness) the common ancestor distribution. Secondly, given such a distribution, it is clear that we can construct a stationary version of the retrospective process (z t ,p t ) which is defined for all times t ∈ R. However, because this lineage persists indefinitely, it is necessarily the common ancestor lineage for the whole population. Accordingly, we can characterize the joint law of the stationary process of substitutions to the common ancestor and the forward-in-time evolution of the allele frequencies by determining the law of the time reversal of the retrospective process with respect to its stationary distribution. (Observe that by time reversing the retrospective process, which runs from the present to the past, we obtain a process which runs from the past to the present.)

The common ancestor distribution
In this section we show that the common ancestor distribution, which we denote π(z, dp), can be found by solving a simple boundary value problem. We begin by observing that because D(G) = C 2 c (E) is an algebra which is dense inĈ(E) and because the martingale problem forG is well-posed, any distribution π(z, dp) which satisfies the condition, 1 0G φ(1, p) π(1, dp) + 1 0G φ(2, p) π(2, dp) = 0 (5) for all φ ∈ D(G), is a stationary distribution forG [Ethier and Kurtz (1986), Chapter 4, Proposition 9.17]. Assuming that we can write π(z, dp) = π(z, p)dp for z = 1, 2, where π(z, ·) ∈ C 2 ((0, 1)), integration-by-parts shows that this condition will be satisfied if Here A * is the formal adjoint of A with respect to Lebesgue measure on [0, 1] and is defined by the formula Because the marginal distribution over z ∈ {1, 2} of the stationary measure π(z, dp) is just the stationary measure π(p)dp of the diffusion process itself, it is convenient to write π(z, dp) in the form π(1, dp) = π(1, p)dp = h(p)π(p)dp π(2, dp) = π(2, p)dp = (1 − h(p))π(p)dp, where h(p) is the conditional probability that the common ancestor is of type P given that the frequency of P in the population is p. Substituting this expression into (6) leads to the following boundary value problem (BVP) for h(p), We show below that the smoothness of the selection coefficient σ(p) is sufficient to guarantee the existence of a solution h(p) to (9) which is smooth in (0, 1) and which has a derivative h ′ (p) that can be continuously extended to [0, 1], and that this implies that the common ancestor distribution can always be represented in the form (8), with h(p) the unique solution to (9). However, we first make two observations concerning equation (9) itself. First, if σ(p) ≡ 0, i.e., P and Q are selectively neutral, then h(p) = p solves (9) and the distribution of the common ancestor is the same as that of an individual sampled randomly from the population. Of course, this claim can also be deduced directly from the look-down formulation of Donnelly and Kurtz (1996): under neutrality, the common ancestor is the individual occupying the lowest level and, by exchangeability, the distribution of the type of this individual is given by the empirical measure carried by all of the particles, which is just pδ 1 + (1 − p)δ 0 for a biallelic model.
, then a simple calculation shows that h(p) will satisfy the BVP (9) if and only if ψ(p) satisfies the following BVP: This result is useful when numerically calculating h(p) because it replaces the divergent inhomogeneous term on the right-hand side of (9) with a term which is smooth on [0, 1]. Even so, because the inhomogeneous equation is singular at p = 0, 1, the usual shooting method (Press et al. (1992)) used to solve such two-point BVP's must be modified as we discuss briefly in the appendix. More importantly, we can use the BVP (10) to prove the existence and regularity of the conditional probability h(p).
Suppose that A is the generator of a Wright-Fisher diffusion as in (1). Then there exists a function ψ(p) satisfying the BVP (10) which is holomorphic on (0, 1) and its first derivative ψ ′ (p) can be continuously extended to [0, 1]. Furthermore, the function h(p) = p+ψ(p) is the unique solution to the BVP (9) sharing these regularity properties.
Proof. We begin by noting that p = 0 and p = 1 are regular singular points for the corresponding homogeneous equation and that the indicial equations have roots λ = 1, −2µ 1 at p = 0 and λ = 1, −2µ 2 at p = 1. Because the coefficients are smooth in (0, 1), Theorems 7 and 8 of Chapter 9 of Birkhoff and Rota (1989) can be used to deduce the existence of four functions, u 0,1 (·) and u 0,2 (·), analytic in a neighborhood of p = 0, and u 1,1 (·), and u 1,2 (·), analytic in a neighborhood of p = 1, as well as two constants C 0 and C 1 such that the following two pairs of functions, and each constitutes a set of linearly independent solutions to the homogeneous equation. Furthermore, because the diffusion operator A is uniformly elliptic on any interval (ǫ, 1 − ǫ) for any ǫ > 0, these solutions can be analytically continued to (0, 1).
Consequently, by taking suitable linear combinations of the ψ ij (·), we can construct a pair of linearly independent solutions, ψ 0 (·) and ψ 1 (·), analytic on (0, 1), such that for any ǫ > 0, A solution to the inhomogeneous equation can then be obtained by the method of variation of parameters, which gives: where W (p) is the Wronskian of the homogeneous equation and p 0 is an arbitrary point in (0, 1). Furthermore, in light of the boundary behavior of the functions ψ 1 (·) and ψ 0 (·), it is easy to check that ψ(·) is smooth in (0, 1), that ψ(0) = ψ(1) = 0, and that the limits exist and are finite. Clearly, these statements also hold for h(p) = p + ψ(p) and a simple calculation verifies that h(p) solves the BVP (9).
By combining this result with the formal calculations leading from Eq. (5) to (9), as well as Proposition 9.17 of Chapter 4 of Ethier and Kurtz (1986), we can deduce the existence of a stationary distribution forG. Our next proposition asserts that this distribution is also unique.
Proposition 2.4. The retrospective process (z t ,p t ) has a unique stationary distribution π(z, dp) of the form (8), where π(p) is the density (2) of the stationary distribution for the Wright-Fisher diffusion generated by A and h(p) is the unique solution to the BVP (9).
Proof. Since we have already demonstrated the existence of a stationary distribution corresponding to Eq. (8)- (9), we need only show that this measure is unique. To do so, we will prove that G is strongly connected (see Donnelly and Kurtz 1999, Section 9): if P ν (t) denotes the onedimensional distribution at time t of the solution to the martingale problem for (G, ν), then for any pair of distributions ν 1 , ν 2 ∈ P(E) and all times T > 0, P η 1 (T ) and P η 2 (T ) are not mutually singular. Uniqueness of the stationary distribution will then follow from Lemma 5.3 of Ethier and Kurtz (1993), which implies that if the embedded Markov chain, ((z nT ,p nT ) : n ≥ 1), has two distinct stationary distributions (as it will if the continuous time process has two distinct stationary distributions), then it also has two mutually singular stationary distributions.
t ) be solutions to the D E [0, ∞)-martingale problem forG with initial distributions ν 1 and ν 2 , respectively. Because the marginal processesp are Wright-Fisher diffusions corresponding to A, the positivity of the mutation rates µ 1 , µ 2 implies that for any t > 0, the one-dimensional distributions ofp (1) t andp (2) t are mutually absolutely continuous with respect to Lebesgue measure on [0, 1]. (In particular, these distributions do not have atoms at 0 or 1.) Furthermore, for every δ ∈ (0, 1/2) and every T > 0, there exists an ǫ > 0 such that the probabilities P{p Combining this observation with the fact that for fixed δ ∈ (0, 1/2), the jump rates of the componentz t are uniformly bounded above 0 and below ∞ whenever the frequency processp t is in [δ, 1 − δ], it follows that P ν 1 (T ) and P ν 2 (T ) are each mutually absolutely continuous with respect to the product measure (δ 1 (dz) + δ 2 (dz)) × m(dp) on E, where m(dp) is Lebesgue measure restricted to (0, 1). Since this implies that P ν 1 (T ) and P ν 2 (T ) are mutually absolutely continuous with respect to one another for every T > 0, the proposition follows.
We can also rewrite the inhomogeneous differential equation in (9) in a form which leads to an alternative probabilistic representation of h(p). Because h(0) = 0 and h(1) = 1, the solution h(p) to the BVP (9) is a harmonic function for the operatorÂ, defined aŝ Setting we see thatÂ is the generator of a jump-diffusion process,p t , which diffuses in (0, 1) according to the law of the Wright-Fisher diffusion corresponding to A until it jumps to one of the boundary points {0, 1} where it is absorbed. It follows from Lemma 2.1 that if the process does reach 0 or 1, then it is certain to have arrived there via a jump rather than by diffusing, even if that boundary is accessible to the pure diffusion process. Indeed, the existence of a unique Markov processp t corresponding toÂ can be deduced from Lemma 2.1 in precisely the same way that the existence and uniqueness of the structured coalescent was obtained, although it is now essential thatp t be absorbed once it hits the boundary. Furthermore, because the total rate of jumps to either boundary point from any point in the interior is bounded below by µ 1 ∧ µ 2 , the process is certain to jump to the boundary in finite time. Taken together these observations lead to the following representation for h(p).
Proposition 2.5. Letp(t) be the jump-diffusion process corresponding to the generatorÂ, and let τ = inf{t > 0 :p t = 0 or 1} be the time of the first (and only) jump to the boundary {0, 1}.
Then the solution h(p) to the BVP (9) is the probability thatp t is absorbed at 1 when starting from initial value p: and so we can use the optional sampling theorem and the fact that h(0) = 0 and h(1) = 1 to calculate Proposition 2.5 has several interesting consequences. First, by comparing the generatorÂ shown in (11) with the generator of the structured coalescent (3) for a sample of size two with one P allele and one Q allele, it is evident that the type of the common ancestor has the same distribution as the type of the sampled lineage which is of the more ancient mutant origin. In other words, the quantity h(p) is the probability that if we sample a P allele and a Q allele from a population in which P occurs at frequency p, then the Q allele has arisen from a mutation to an ancestral P individual more recently than the P allele in the sample has arisen from a mutation to an ancestral Q individual.
Secondly, because the rate at whichp t jumps to 1 is a strictly increasing function of p while the rate at whichp t jumps to 0 is strictly decreasing, (12) implies that h(p) is a strictly increasing function of p. While we would expect such a relationship to hold if the selection coefficient σ(p) is either non-decreasing or non-negative, it is noteworthy that h(p) is increasing even with negative frequency-dependent selection, e.g., under balancing selection, with σ(p) = s · (p 0 − p) for s > 0 and p 0 ∈ (0, 1).
Another consequence of Proposition 2.5 is that the probability that the common ancestor is of a particular genotype is an increasing function of the fitness of that genotype. To make this precise, suppose that A (1) and A (2) are a pair of Wright-Fisher generators as in (1) which differ only in their (smooth) fitness functions, σ 1 (p) and σ 2 (p), respectively, letÂ (i) , i = 1, 2 be the generators of the jump-diffusion processes obtained by taking A = A (i) in (11), and let h 1 (p) and h 2 (p) be the corresponding conditional probabilities that the common ancestor is of type P .
Proof. In view of the smoothness of the coefficients of the diffusion generators A (i) , i = 1, 2, there exists a probability space (Ω, F, P), a Brownian motion (W t , t ≥ 0), and diffusion processes (p (i) t , t ≥ 0), i = 1, 2 corresponding to these generators such that the stochastic differential equation is satisfied a.s. for i = 1, 2 and all p ∈ [0, 1]. (For example, such a coupling can be constructed using a sequence of coupled Markov chains which converge weakly to these diffusions.) Furthermore, because the drift coefficients satisfy the inequality b 1 (p) ≤ b 2 (p) for all p ∈ [0, 1], while the infinitesimal variance satisfies the regularity condition | p(1 − p) − q(1 − q)| < 2|p − q| 1/2 for all p, q ∈ [0, 1], we can use Lemma 3.4 in Shiga (1981) to conclude that To relate this inequality to the jump-diffusion processes generated by theÂ (i) , i = 1, 2, observe that because each process jumps exactly once, we can construct coupled versions of these processes, denotedp Here, Z 1 and Z 2 are unit mean exponential random variables which are independent of each other and of the diffusions p  Since the rate function governing jumps from (0, 1) to 1 is an increasing function of p ∈ [0, 1], while that governing jumps from (0, 1) to 0 is a decreasing function of p, the inequality in (13) implies that, with probability one, if the processp (1) t jumps to 1, then the processp (2) t must also have jumped to 1, possibly at an earlier time. Consequently, τ (2) = 1} = h 2 (p), and the proposition follows upon noting that the initial condition p ∈ [0, 1] is arbitrary.
In particular, taking σ 2 (p) ≥ σ 1 (p) ≡ 0, we can use Proposition 2.6 to conclude that h 2 (p) ≥ h 1 (p) = p. Furthermore, if σ 2 (p) is strictly positive on [0, 1], then the fact that the diffusions p t } has positive Lebesgue measure whenever the initial condition p ∈ (0, 1), and therefore that h 2 (p) > h 1 (p) for every p ∈ (0, 1). In other words, if P is unconditionally more fit than Q, then the common ancestor will be more likely to be of type P than an individual sampled at random from the population, both on average and when conditioned on the frequency p at which P is segregating in the population. Furthermore, this property implies that the mean fitness of the common ancestor is greater than the mean fitness of an individual chosen at random from the population, and generalizes Theorem 2 of Fearnhead (2002) which applies when P has a fixed (frequency-independent) advantage over Q.

The common ancestor process
Having found the common ancestor distribution, our next task is to identify the common ancestor process, which we will do by determining the time-reversal of the retrospective process (z t ,p t ) with respect to its stationary distribution. Because Proposition 2.4 asserts that this distribution is unique, the common ancestor process, at least in a stationary population, is also unique. We recall that time reversal preserves the Markov property, and that the generator G of the Markov process obtained by time reversal of the stationary process corresponding toG has the property that it is adjoint toG with respect to the measure π(z, dp) (Nelson 1958), i.e., z=1,2 1 0 G φ(z, p) ψ(z, p)π(z, p)dp = z=1,2 1 0 φ(z, p) Gψ(z, p) π(z, p)dp, for any ψ ∈ D(G) and φ ∈ D(G) (which is to be determined). A calculation making repeated use of the product rule and integration-by-parts, along with the characterization of the common ancestor distribution π(z, p)dp provided by Proposition 2.4 and the fact that A * π(p) = 0, with the density π(p) given by (2) and the adjoint operator A * given by (7), shows that this condition will be satisfied if Proposition 2.7. For any ν ∈ P({1, 2} × [0, 1]), there exists a Markov process (z t , p t ), which we call the common ancestor process, which is the unique solution to the martingale problem for (G, ν). Furthermore, (z t , p t ) ∈ E for all t > 0.
Proof. Since h(0) = 0 and h(1) = 1, the continuity of h ′ (p) on [0, 1] implies the existence of constants Consequently, all of the terms appearing on the right-hand side of (15) can be continuously extended (as functions of p) to [0, 1], and we define Gψ(z, p) accordingly if (z, p) = (1, 0) or (z, p) = (2, 1). In particular, the operators To prove that (z t , p t ) ∈ E for all t > 0, we observe that boundary points inconsistent with the type of the common ancestor are entrance boundaries for the frequency component of the common ancestor process. For example, if the type of the common ancestor is P , i.e., if z = 1, then because h(p) is continuously differentiable on [0, 1] (Lemma 2.3) and because p/h(p) ≈ 1/h ′ (0) when p ≈ 0 we can write the drift coefficient of where b(p) is the drift coefficient of A. That p = 0 is an entrance boundary for the diffusion corresponding to A 1 can then be shown using Feller's boundary classification (Ewens 2004, Section 4.7) and the fact that µ 1 + 1 > 1/2. Similar remarks apply to the boundary p = 1 and the diffusion corresponding to A 2 .
If A is the generator of a neutral Wright-Fisher diffusion (i.e., σ(p) ≡ 0), then h(p) = p and the generator of the common ancestor process is just As expected, the process governing the change of type of the common ancestor decouples from the frequency process and coincides with the mutation process itself. The only novel feature of the neutral common ancestor process is the presence of the additional drift terms in the diffusion which reflect the fact that because the common ancestor contributes more offspring to the population than an individual chosen at random, the population has a tendency to evolve towards the type of the common ancestor. Indeed, these extra births can be made explicit by formulating a finite population model (z (N ) , p (N ) ) which combines the usual Moran resampling with a neutral look-down process that operates only on the lowest level (i.e., birth-death events involving the lowest level always assign the birth to the lowest level, but all other birth-death events are resolved by choosing the parent at random from the two participating individuals). It is then straightforward to show that as N → ∞, suitably rescaled versions of (z (N ) , p (N ) ) converge weakly to the jump-diffusion process generated by G.
When there are fitness differences between the two alleles, in general h(p) = p and the substitution rates to the common ancestor depend on the allele frequency p, i.e., the substitution process to the common ancestor z t is not a Markov process. In this case, the substitution rates will differ from the corresponding mutation rates for most values of p. Moreover, because Proposition 2.6 shows that h(p) > p for all p ∈ (0, 1) whenever P is unconditionally more fit than Q (i.e., σ(p) > 0 for all p ∈ [0, 1]), Eq. (15) shows that the rate of substitutions from the less fit allele to the more fit allele is greater than the corresponding mutation rate, and vice versa. A less intuitive property of the generator of the common ancestor process is that for each value of p the geometric mean of the two substitution rates is the same as that of the two neutral mutation rates. While it is unclear what biological interpretation this invariant might have, one mathematical consequence is that for each fixed value of p only one of the two substitution rates can exceed the corresponding neutral mutation rate, while the other is necessarily less than it. However, the direction of these two inequalities may differ according to the frequency p if selection is frequency-dependent or fluctuates in time.

Weak mutation limits
Because single nucleotide mutation rates in DNA sequences are typically on the order of 10 −8 mutations per site per generation, while most effective population size estimates are less than 10 7 (Lynch and Connery 2003), the asymptotic properties of the common ancestor process in the limit of vanishing mutation rates are of special interest. (Here we temporarily relax our earlier assumption that the mutation rates are positive.) We first observe that if µ 1 and µ 2 are both zero, then the BVP (9) simplifies to the equation Ah(p) = 0, with h(0) = 0 and h(1) = 1, and the solution, is just the fixation probability of P when its initial frequency is p (Ewens 2004, Section 4.3). Furthermore, if we substitute this expression into the generator of the common ancestor process (15), then because both jump rates vanish, we are left with a pair of operators, which we recognize to be the generators of the diffusion process corresponding to A conditioned to absorb either at 1 (top line) or at 0 (lower line) (Ewens 2004, Section 4.6). That the limiting generator takes this form reflects the fact that in the absence of mutation, any population which is descended from the common ancestor will also be fixed for the type of that individual.
A more useful observation is that if the mutation rates are small enough that mutations occur rarely on the coalescent time scale, then we can approximate the non-Markovian substitution process to the common ancestor by a continuous time two state Markov chain. Although approximate, such a process would greatly simplify the numerical or Monte Carlo computations needed to infer selection coefficients and other model parameters from a set of DNA sequences. One possibility is to define the transition rates of the Markov chain to be equal to the mean substitution rates obtained by averaging the frequency-dependent substitution rates of the bivariate process (15) over the conditional distribution of the allele frequencies given the type of the common ancestor, e.g., µ CA 2 is the mean substitution rate to the common ancestor given that the type of that individual is P (which mutates to Q). Indeed, the ergodic properties of Wright-Fisher diffusions (Norman 1977) offer some justification for this approximation. Provided that the mutation rates are sufficiently small, the time elapsed between successive substitutions to the common ancestor will with very high probability be large enough for the allele frequencies to have relaxed to their stationary distribution well in advance of the next mutation. Moreover, because Lemma 2.3 guarantees that the jump rates appearing in the generator G in (15) are continuous functions of the frequency p, the time averages of the jump rates along paths of the diffusion will be approximately equal to the product of the time elapsed and the mean substitution rates shown above. Of course, for this approximation to be relevant to data, we will also need the phylogenetic tree describing the relationships among the sequences to be deep enough that the ergodic averages of the substitution rates are approached along each branch of the tree.
When the mutation rates are very small, the average substitution rates shown in (18) can be replaced by simpler expressions which depend only on the mutation rates and the fixation probabilities. Suppose that µ i = θν i , i = 1, 2, and write h θ (p), π θ (p), and µ CA i,θ to indicate the dependence of these quantities on θ. In Proposition 3.2 we evaluate the scaled, weak mutation limits µ CA i,low ≡ lim θ→0 θ −1 µ CA i,θ . However, we first state a technical lemma which will be needed in the proof of that proposition. Recall that we assume that the selection coefficient σ(p) is holomorphic on some neighborhood of [0, 1].
To make this argument precise, let us introduce the hitting times T q,θ = inf{t > 0 : p θ (t) = q} for q ∈ [0, 1], with T q,θ = ∞ if p θ (t) = q for all t > 0, and recall that P p {T b,θ < T a,θ } = (s θ (p) − s θ (a))/(s θ (b) − s θ (a)) for any 0 < a < b < 1, where the scale function s θ (p) for the Wright-Fisher diffusion p θ (t) is . Furthermore, if 2θν 1 and 2θν 2 are both less than 1, then the scale function is finite on [0, 1] and we can also allow a = 0 and b = 1 in the previous expression for the hitting probability. Consequently, for every p ∈ [0, 1], s θ (p) converges pointwise (in fact, uniformly on [0, 1]) to s 0 (p) as θ → 0, and for any fixed 0 ≤ a < b ≤ 1, the probabilities P p {T b,θ < T a,θ } converge to P p {T b,0 < T a,0 }. In particular, if we define u θ (p) = P p {T 1,θ < T 0,θ } to be the probability that the diffusion p θ (t) hits 1 before hitting 0, then u θ (p) converges uniformly to h 0 (p) on [0, 1]. We also observe that if we let T θ ≡ T 0,θ ∧ T 1,θ denote the first hitting time of 0 or 1 by the diffusion p θ (t), then the expectation E p [T θ ] is finite whenever 2θν 1 and 2θν 2 are both less than 1, in which case
The proof of Proposition 3.2 shows that the weak mutation limits µ CA i,low are closely related to an approximation commonly used to describe the 'flux of selected alleles' (Kimura 1964;Otto and Whitlock 1997) and incorporated into a phylogenetic framework by McVean and Vieira (2001).
, we see that the limiting substitution rate of P is approximately equal to the product of the number, N ν 1 , of new P mutants produced per generation and the fixation probability, h 0 (N −1 ), of a single such mutant in a population otherwise fixed for Q. In contrast, if we let u θ (p) denote the fixation probability of P when the mutation rates are θν 1 and θν 2 , then it is not true that ν 1 u ′ θ (0) converges to ν 1 h ′ 0 (0) as θ → 0 since u ′ θ (0) = ∞ whenever θ > 0 (see the scale function s θ (p) introduced in the proof of Lemma 3.1). Thus the additional regularization of h ′ θ (p) afforded by recurrent mutation to the common ancestor lineage (represented by the 'jump terms' in the BVP (9)) appears to be essential to the existence of the low mutation rate limit. We shall see in the next two sections that the approximation given by Proposition 3.2 is generally very good when selection and mutation are both weak, but tends to underestimate the substitution rates if either selection is strong or the mutation rates are high.

Purifying selection in a haploid population
We next show how the theory developed in the preceding section can be used to characterize the common ancestor process of a haploid population evolving according to a Wright-Fisher diffusion (1) with frequency-independent fitness differences between the alleles, i.e., σ(p) ≡ s = 0. Because we know from (2) that the density of the stationary distribution of this diffusion is π(p) = Cp 2µ 1 −1 (1 − p) 2µ 2 −1 e 2sp , our description of the common ancestor distribution will be complete if we can solve (9) for the conditional probability h(p).
To do so, we begin by supposing that we can expand h(p) in a power series in s Substituting this expansion into (9) and collecting all terms multiplying s n leads to a recursive series of BVP's for the functions h n (·), n ≥ 1, subject to the conditions h n (0) = h n (1) = 0. To solve these inhomogeneous equations, we first need to determine the general solution to the corresponding homogeneous equation. Some guesswork leads to one solution and a reduction of order calculation leads to a second linearly independent solution With these in hand, integration-by-parts and the method of variation of parameters can be used to find a recursive solution to the boundary value problem given in (22) h n (p) = 2 β ′ (p) Defining H n (p) = β ′ (p)h n (p) and H(p) = β ′ (p)h(p), we can rewrite this recursion as which, upon differentiating with respect to p, gives Term-by-term differentiation of the series expansion of H(p) itself leads to the following firstorder differential equation and to find h(p), we must divide the general solution to this equation by β ′ (p) and impose the original boundary conditions h(0) = 0 and h(1) = 1 (which can be jointly satisfied). These calculations lead to the following expression for the conditional probability that the genotype of the common ancestor is P , where the constantp is the expectation of the allele frequency p with respect to the variancebiased stationary distributionπ(p)dp ≡ Csp(1 − p)π(p)dp (where C is a normalizing constant): (Observe thatp is also the probability that a sample of three individuals from a stationary population contains two P and one Q individual conditional on it containing at least one individual of each genotype.) We can calculate the marginal probability, π 1 , that the common ancestor is of type P by integrating the density of the joint probability π(1, p) = h(p)π(p) over [0, 1]. Because this integral cannot be evaluated analytically, π 1 must be calculated by numerical integration, which can be done accurately using the method described in the appendix. Furthermore, by interchanging the order of integration in the resulting double integral, we arrive at the following intriguing expression for π 1 where Covπ (·, ·) denotes the covariance with respect to the variance-biased stationary measure defined above. Although this expression is reminiscent of Price's equation (Price 1970), which states that the change in a trait caused by selection is equal to the covariance between that trait and fitness, it is not clear how to interpret the terms appearing within the covariance in a way that would make this correspondence precise.
When s > 0, i.e., P is fitter than Q, it is clear that the integral on the right-hand side of (23) vanishes at p = 0, 1 and is strictly positive when p ∈ (0, 1). Consequently, h(p) ≥ p, as follows from Proposition 2.6, and thus the common ancestor is more likely to be of the fitter type than an individual chosen at random. Plots of h(p) for different values of the (symmetric) mutation rates and selection coefficients are shown in Figure 1A. For fixed values of the mutation rates, we see that h(p) is an increasing function of the selection coefficient s, which also follows from Proposition 2.6. On the other hand, for fixed positive values of s, h(p) is a decreasing function of the mutation rates, probably because mutation reduces the correlation between the type of an extant lineage and its probability of surviving into the future.
Although expressions (15) and (23) fully determine the generator G of the common ancestor process, none of the terms containing h(p) simplify and so we do not reproduce these here. Less cumbersome, approximate expressions for the substitution rates can be derived with the help of Proposition 3.2, which shows that the weak mutation limits are µ CA 2,low = µ 2 2s e 2s − 1 and µ CA 1,low = µ 1 2se 2s e 2s − 1 .
These are also derived in Corollary 3 of Fearnhead (2002) and have been used by McVean and Vieira (2001) to estimate the strength of selection on codon usage in several Drosophila species. Of course, we can also use expressions (15) and (23) to calculate the exact common ancestor substitution rates. Figure 1B shows how the relative deleterious substitution rate, , varies as a function of the frequency p of P . (As can be seen in (15), the relative beneficial substitution rate, µ CA 1 (p)/µ 1 , is always the reciprocal of this quantity and so is not shown.) Note that the mutation rates are symmetric in parts A-C of Figure 1, i.e., µ 1 = µ 2 ≡ µ, and that the substitution rates are scaled by the mutation rate. As expected, the relative deleterious substitution rate is always less than 1, i.e., the absolute substitution rate is less than the mutation rate, and this rate decreases as the selective advantage of P increases, but increases as the mutation rate µ increases. For comparison, we have also plotted the average deleterious substitution rates, µ CA 2 /µ 2 , calculated using (18) and scaled by µ, as bold horizontal line segments on the right side of Figure 1B. Examining this figure reveals that for each fixed pair of values of µ and s, the average deleterious substitution rate is nearly as small as the smallest frequency-dependent rate (i.e., the bold horizontal lines lie beneath the corresponding curve for most values of p). Presumably this is because the conditional distribution of p given that the common ancestor is of type P is concentrated in a small region abutting the boundary p = 1 whenever P is selectively advantageous and the mutation rate is not too large.  Figure 1C, we plot the average deleterious and favorable substitution rates, again scaled by the corresponding (symmetric) mutation rates, as functions of the selective advantage of P . (Because the conditional distributions used in (18) to define these two rates differ, the reciprocity property noted for the frequency-dependent rates no longer holds.) Also shown are the weak mutation limits, calculated using (21), which are shown as bold curves. It is notable that except when the mutation rate is very large (µ = 0.5), both the beneficial and the deleterious substitution rates are greater than their weak mutation limits, although the discrepancy is greatest for deleterious substitutions and grows as the selection coefficient s increases. For example, when s = 5 and either µ = 0.01 or µ = 0.1, the average deleterious substitution rates are approximately 7 and 59-times greater than the corresponding weak mutation limit, respectively, while the average beneficial substitution rates for these two cases are approximately 7 and 3 times as large. Furthermore, because Figure 1B shows that the average substitution rates are nearly as small as the minimum frequency-dependent rates, we can conclude that the discrepancy between the average substitution rates and their weak mutation limits arises primarily because the limiting values underestimate the true substitution rates and not because the average substitution rates are too large. While this difference is small when the mutation rates are small, it could be as much as an order of magnitude or more in organisms having either very large effective population sizes or high mutation rates, e.g., in HIV-1 the mutation rate is approximately 3 × 10 −5 mutations per nucleotide per viral generation, while the effective viral population size within infected hosts is usually estimated to be between 10 3 − 10 4 (Koyous et al. 2006), giving µ ≈ 3 × 10 −2 − 3 × 10 −1 . In such cases, the stationary averages given by (18) will be much more accurate numerical summaries of the true, frequency-dependent substitution rates than approximations which neglect recurrent mutation.
To disentangle the effects of the two mutation rates, Figure 1D shows how the two substitution rates change when one of the mutation rates is held fixed and the other is varied. The behavior of the deleterious substitution rate is easiest to understand: increasing either mutation rate increases this substitution rate, although the effect is greatest when the mutation rate to the favorable allele is increased, presumably because a Q lineage is then more likely to mutate to P before going extinct. In contrast, the beneficial substitution rate eventually decreases when either mutation rate is increased, but is initially an increasing function of µ 1 , possibly because of mutation-drift interactions, i.e., a larger mutation rate helps drive a rare favorable allele up to frequencies where selection can be effective compared with genetic drift.

Complementarity of the diffusion and graphical representations
The common ancestor process for a Wright-Fisher diffusion with frequency-independent selection has also been characterized by Fearnhead (2002) using the ancestral selection graph. This approach embeds the common ancestor process within a pure jump Markov process (z t , n t ) taking values in the state space E = {0, 1} × N , where n t denotes the number of virtual lineages, all of the less fit type. Fearnhead (2002) shows that the stationary distribution of this process is given by where E π [·] denotes an expectation with respect to the stationary measure π(dp) and λ n = lim k→∞ λ and λ (k) k+1 = 0, and we interpret empty brackets (n = 0) as being equal to 1. (This formula assumes that s ≥ 0; if P is less fit than Q, then we simply exchange indices. Also note that (26) and (27) have been rewritten to reflect the scalings of µ i and s used in this article rather than those in Fearnhead (2002).) The transition rates of the common ancestor process can be calculated by reversing the modified ancestral selection graph with respect to π(z, dn) and also depend on the λ n 's; we refer the reader to Corollary 2 of Fearnhead (2002) for their values.
Because the marginal laws of the genealogical processes embedded within the structured coalescent and the ancestral selection graph are identical, it is clear that this will also be true of the common ancestor processes identified by these two methods. However, deducing this equality directly from the generators of the bivariate processes appears to be difficult and here we merely show that the marginal stationary distributions of the type of the common ancestor are the same.
Lemma 4.1. Let h(p) be the conditional probability defined in (23) and let (λ n ) n≥1 be the sequence defined by (27). Then, Proof. Using recursion (27) and its initial condition, we can show inductively that, for every n ≥ 1 and every k ≥ 1, λ (k) n < s/(s + µ 1 ), and therefore that λ n ≤ s/(s + µ 1 ) and a n ≡ n i=1 λ i ≤ (s/(s + µ 1 )) n as well. It follows that the function g(p) = p n≥1 a n (1 − p) n is holomorphic in the open disk D = D(1; 1 + ǫ) for some ǫ > 0 and thus can be differentiated term-by-term on D. Substituting g(p) into the left-hand side of equation (10), we obtain the following expression and, using recursion (24), which holds with λ n in place of λ (k) n , we can show that this is equal to −sp(1 − p) for all p ∈ [0, 1]. Since g(0) = g(1) = 0, the uniqueness of solutions to second-order boundary value problems with smooth coefficients implies that h(p) = p + g(p).
The equality of the marginal stationary distributions follows upon integrating both sides of the identity asserted by the lemma with respect to the stationary measure π(dp), Another consequence of the lemma is that it provides an explicit formula for the constants λ n , where v(p) ≡ (h(p) − p)/p. Of course, the algebraic operations needed to analytically evaluate successive derivatives of v(p) essentially replicate the recursion satisfied by the λ n . However, the one advantage that (28) does have over (27) is that it gives an explicit formula for λ 1 , which allows the recursion to be solved from the bottom-up, starting with the calculated value of λ 1 , rather than from the top-down, with the approximation λ n ≈ 0 for some large value of n.

Selection and dominance in a diploid population
To illustrate the generality of the diffusion theoretic methods, in this section we will consider a diploid population and explore what effect the degree of dominance of fitness has on the common ancestor process. Several observations suggest that dominance plays an important role in molecular evolution. For example, it has long been known that deleterious mutations in coding sequences are usually recessive, possibly because of complementation by the fully functional allele or because of structural features of metabolic pathways (Kondrashov and Koonin 2004). Furthermore, even when the different alleles segregating at a locus are not individually advantageous, non-additive interactions between alleles can cause heterozygous genotypes to have higher or lower fitness than any of the possible homozygotes, leading to balancing or disruptive selection (Richman 2000). From classical population genetics theory we know that dominance relations affecting fitness can profoundly alter fixation probabilities and rates (Ewens 2004;Williamson et al. 2004) and so we would expect the same to be true of the substitution process of the common ancestor. We formulate our model by considering a diploid population of effective population size 2N e in which the relative fitnesses of the genotypes P P , P Q, and QQ are 1 + 2s : 1 + 2ds : 1, respectively, and d is a constant which quantifies the dominance (d > 0.5) or recessiveness (d < 0.5) of P relative to Q. Note that when d > 1, heterozygotes have higher fitness than either homozygote and are said to be overdominant, whereas when d < 0, heterozygotes have lower fitness and then are said to be underdominant. By rescaling both the selection coefficient s and the mutation rates by a factor of 1/2N e and speeding up time by a factor of 2N , we can again approximate the changes in the frequency of P by a Wright-Fisher diffusion with generator (1) where σ(p) = 2s(d − (2d − 1)p). When d = 0.5, σ(p) ≡ s is constant and the common ancestor process can be characterized using either the results in the preceding section or those of Fearnhead (2002). However, for any other value of d, σ(p) is frequency-dependent and neither set of results applies.
Because the ancestral selection graph has been identified for this diffusion model (Neuhauser 1999), one might try to identify the common ancestor process by generalizing the methods used by Fearnhead (2002). The main obstacle to implementing this approach is that trinary branchings are required to account for the frequency-dependence of fitness and it is unclear that S there is a pruned version of the ancestral selection graph of the kind found in Fearnhead (2002). However, without such a simplification, the common ancestor process will be embedded within a trivariate process (z t , n 1 (t), n 2 (t)), where n i (t) is the number of virtual lineages of type A i , and the stationary distribution will have to be determined by solving a multi-dimensional recursion.
In contrast, the diffusion theory developed in the first part of this article can be applied without modification to the new model. The price we pay for the added complexity of frequencydependent selection is that we can no longer write down an explicit formula for h(p): although we can still solve the BVP (9) recursively as in the previous section, we no longer obtain an integrating factor for the original equation. On the other hand, it is relatively easy to solve this problem numerically using the shooting method, even for large values of |s|. (See the appendix for a description of our numerical methods.) Furthermore, because the density π(p) of the stationary distribution is known explicitly, we can use our numerical estimates of h(p) to evaluate both the density of the common ancestor distribution and the substitution rates to the common ancestor. Of course, we can also use Proposition 3.2 to calculate the weak mutation limits to the substitution rates. Figure 2 shows the results of these calculations. In Figure 2A we have plotted the numerical solutions themselves to show how h(p) varies as a function of the dominance coefficient d when selection is moderately strong (s = 5). Qualitatively similar results are obtained both with weaker (s = 1) and stronger (s = 10) selection and so are not shown. We then substituted these numerical values into (15) to obtain the frequency-dependent 'deleterious' substitution rates (P → Q) shown in Figure 2B. (Note that these substitutions are unconditionally deleterious only when the dominance coefficient d lies between 0 and 1.) The patterns evident in both figures can be interpreted by considering how dominance affects the marginal fitnesses of the two alleles at high and low frequencies of P . With increasing dominance of P over Q, the difference between the marginal fitnesses of the two alleles is reduced at high frequencies of P , rendering selection less effective and causing both a small decline in h(p) but also a marked increase in the deleterious substitution rate at frequencies p close to 1. In contrast, because higher levels of dominance expose heterozygotes to stronger selection in populations which are nearly fixed for Q, these relationships are reversed when p is close to 0.
Similar considerations apply when heterozygotes are over-or under-dominant. One interesting feature of Figure 2B is that with disruptive selection (d = −0.5) the relative substitution rate of P → Q increases above 1 when the frequency of P is sufficiently close to 0. This is because the deterministic dynamics corresponding to this fitness scheme have an unstable internal equilibrium below which the marginal fitness of P is less than that of Q and selection favors Q substitutions. In contrast, when heterozygotes are over-dominant, the corresponding deterministic dynamics have a stable internal equilibrium and the marginal fitness of Q is an increasing function of the frequency of P , leading to the convex substitution rates seen in Figure 2B when d = 1.5.
An overview of how dominance affects the substitution process can be gleaned from Figure 2C, which shows plots of the (relative) average substitution rates for different values of d and µ, as well as the low mutation rate limits. We again see that the weak mutation rate limits (21) generally underestimate the average substitution rates (18), except when the mutation rate is so large that the long-term fitness of a lineage is partially decoupled from its current type. Also, whereas the deleterious (P → Q) substitution rates are increasing functions of the dominance coefficient d, the beneficial (Q → P ) substitution rates are either unimodal or decreasing. In fact, even when µ is as small as 0.01, the beneficial substitution rate is seen to decrease slightly as d exceeds 1.3, probably because heterozygote advantage favors Q when the frequency of P is very high.

Multiple genetic backgrounds: prospects and problems
In this article we have used the structured coalescent in a fluctuating background to characterize the common ancestor process associated with a class of diffusion models important in population genetics theory. In addition to the classical Wright-Fisher diffusion, which can model general forms of frequency-dependent selection in panmictic populations, one-dimensional diffusions arise as scaling limits of models incorporating population structure (Cherry and Wakeley 2003), group selection (Roze and Rousset 2004) and environmental variation (Gillespie 1991). Although a closed-form solution was found only for the model with genic selection, the theory can also be used to quantify the influence of selection and genetic drift on the rate of molecular evolution under more complicated scenarios by first solving the BVP (9) numerically and then substituting the results into the expressions for the substitution rates which were derived as part of (15). Furthermore, as the section on the weak mutation limits illustrates, we can also use the theory to obtain analytical approximations for the substitution rates when selection, mutation or genetic drift are either very strong or very weak.
The most serious limitation of the diffusion-theoretic approach is that it leads to a much less tractable description of the common ancestor process when there are more than two genetic backgrounds. To illustrate both the difficulties and the potential interest of this approach, consider a locus which we will call the focal locus and which can occur in m different genetic backgrounds, P 1 , · · · , P m , present at frequencies p 1 , · · · , p m , respectively. Because the frequencies sum to one, we can describe the genetic composition of the population using any m − 1 of these and so we will consider diffusion processes which take values in the (m − 1)-dimensional simplex K m−1 = {(p 1 , · · · , p m−1 ) : p 1 , · · · , p m−1 ≥ 0, p 1 + · · · + p m−1 ≤ 1}. As before, let N e denote the effective population size, let 1 + σ i (p)/N e denote the relative fitness of background P i , and suppose that mutations from P i to P k occur at rate µ ik /N e , and that recombinations (or gene conversion events) involving individuals of type P j change the background of the focal locus from type P i to type P k at rate ρ(i, j|k)/N e . It will also be convenient to define µ ii = ρ(i, j|i) = 0 for all i, j = 1, · · · , m. By rescaling the parameters and time in the usual manner and writinḡ σ(p) = m i=1 p i σ i (p) for the mean fitness of the population, we obtain a Wright-Fisher diffusion with generator (p k p j ρ(k, j|i) − p i p j ρ(i, j|k)) + p i (σ i (p) −σ(p))   ∂ i ψ(p), for ψ ∈ C 2 (K m−1 ). Although we allow for recombination in this model, we emphasize that we are considering the common ancestor process at the focal locus only, and that ρ(i, j|k) is the rate at which recombinations involving a non-ancestral lineage in background P j change the background of the ancestral lineage from P i to P k . We could also define a structured ancestral recombination graph (Griffiths and Marjoram 1997) and use this to characterize the type of the common ancestor at several recombining loci, but this process would be even more complicated than the one we do consider here.
Unfortunately, when we try to write down the generator for the coalescent process of a sample of n genes from a population evolving according to this model, we encounter several complications that do not occur in biallelic models. One is that because the diffusion corresponding to generator (29) need not be not time-reversible with respect to its stationary distribution, i.e., the detailed balanced conditions need not hold, the generatorÃ of the time-reversed process may differ from A. If we denote the density of the stationary distribution of A by π(p) (for which we will assume both existence and uniqueness), then we can use the adjoint condition (Nelson 1958) pzh k (p) p k hz(p) which account for the effects of selection. One new feature of (34) is that the rate at which recombination changes the type of the common ancestor is also influenced by selection, although the effect depends only on the types of the backgrounds A z and A k of the common ancestor before and after the recombination event, and not on the type A j of the individual with which the common ancestor recombines. This suggests that attempts to quantify recombination using phylogenetic methods (e.g., Patterson et al. 2006) could be confounded by selection.
There are a few situations in which the multidimensional Wright-Fisher diffusion is reversible with respect to its stationary distribution, allowing analytical expressions for the conditional probabilities h k (p) and the generator G to be found (Li et al. 1999). Under complete neutrality and parent-independent mutation, we have σ k (p) ≡ 0 and µ ik ≡ µ k for all i, k = 1, · · · , m, and direct substitution into equation (33) shows that h z (p) = p z as expected. In this case, the stationary distribution of background frequencies is known to be the Dirichlet distribution with parameters (2µ 1 , · · · , 2µ m ). Moreover, under complete neutrality, equation (33) shows that it is true that h z (p) = p z even when the mutation rates are parent-dependent, although we are then unable to write down an explicit formula for the density π(p).
If the genetic backgrounds can be partitioned into two fitness classes, say F and U, with fitnesses 1+s and 1, respectively, and mutation is parent-independent, then as in Fearnhead (2002) we can use the solution from the corresponding biallelic model to determine the stationary distribution and generator of the multi-allelic common ancestor process. Suppose that F = {P 1 , · · · , P l } and U = {P l+1 , · · · , P m }, and let µ F = µ 1 + · · · + µ l , µ U = µ l+1 + · · · + µ m , and p F = p 1 + · · · + p l . Then p F (t) evolves according to a Wright-Fisher diffusion with parameters µ F , µ U , and s, and so the probability that the common ancestor belongs to the fitness class F given that the frequency of that class is p is given by equation (23), where we set µ 1 = µ F and µ 2 = µ U . Furthermore, using equation (33), we can show that the multi-allelic conditional distribution of the type of the common ancestor is given by h z (p) = p z p F h(p F ) if z = 1, · · · , l and h z (p) = p z 1 − p F (1 − h(p F )) if z = l + 1, · · · , m.
(35) Since the density of the stationary distribution is given by Wright's formula, it follows that (35) determines both the common ancestor distribution and the generator of the common ancestor process.
Unfortunately, analytical solutions such as these are rarely available, and thus the difficulty of numerically solving the system of singular PDE's in (33) limits the usefulness of this theory. Extensions to models based on multidimensional diffusions are important for several reasons. On the one hand, while neutral substitutions at different sites will occur independently of one another (assuming that the mutation rates are not context-dependent), selection will lead to correlated substitution processes whenever fitness is determined epistatically or when there is genetic linkage between polymorphic loci. It is important to understand and to quantify these correlations not only because they may alter the marginal substitution rates, but also because of the significant role which they might play in processes such as speciation and the evolution of recombination. Furthermore, even if we could assume that the substitution processes at different This problem can be resolved by splitting the singular integrals into three parts, 1 0 F (p)π(p)dp ≡ 1 0 F (p)G(p)p 2µ 1 −1 (1 − p) 2µ 2 −1 dp = ǫ 0 F (p)G(p)p 2µ 1 −1 (1 − p) 2µ 2 −1 dp + 1−ǫ ǫ F (p)G(p)p 2µ 1 −1 (1 − p) 2µ 2 −1 dp where ǫ > 0 is chosen small enough that the locally smooth functions F (p)G(p)(1 − p) 2µ 2 −1 and F (p)G(p)p 2µ 1 −1 can be approximated by F (0)G(0) and F (1)G(1) in the first and third integral, respectively. With this approximation, the boundary integrals can be evaluated analytically, while the non-singular integral over (ǫ, 1 − ǫ) can be evaluated numerically. The accuracy of this scheme was tested by comparing the expected substitution rates for the Wright-Fisher model with genic selection obtained using the diffusion characterization with those reported in Fearnhead (2002) using an independent characterization, and the two sets of rates were seen to agree to within the number of digits reported in the latter paper. A Mathematica program implementing both the shooting and the integration methods described here is available from the author upon request.