n a l r u o o f J P c r i o n b E l e c t r o a b i l i t y Vol.12 (2007), Paper no. 28, pages 808–847. Journal URL http://www.math.washington.edu/~ejpecp/ The Common Ancestor Process for a Wright-Fisher Diffusion Jesse E. Taylor ∗ Department of Statistics, University of Oxford 1 South Parks Road, Oxford, OX1 3TG, United Kingdom [email protected] Abstract Rates of molecular evolution along phylogenetic trees are influenced by mutation, selection andgeneticdrift. Providedthatthe branchesofthe treecorrespondtolineagesbelongingto geneticallyisolatedpopulations(e.g.,multi-speciesphylogenies),theinterplaybetweenthese three processes can be described by analyzing the process of substitutions to the common ancestorofeachpopulation. Wecharacterizethisprocessforaclassofdiffusionmodelsfrom population genetics theory using the structured coalescent process introduced by Kaplan et al.(1988)andformalizedinBartonetal.(2004). Fortwo-allelemodels,thisapproachallows both the stationary distribution of the type of the common ancestor and the generator of thecommonancestorprocesstobedeterminedbysolvingaone-dimensionalboundaryvalue problem. In the case of a Wright-Fisher diffusion with genic selection, this solution can be foundinclosedform,andweshowthatourresultscomplementthoseobtainedbyFearnhead (2002) using the ancestral selection graph. We also observe that approximations which ne- glect 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 ∗Supportedby EPSRCgrant EP/E010989/1. 808 it leads to systems of singular partial differential equations which we have been unable to solve. Key words: Common-ancestor process, diffusion process, structured coalescent, substitu- tion rates, selection, genetic drift. AMS 2000 Subject Classification: Primary 60J70,92D10, 92D20. Submitted to EJP on November 14, 2006, final version accepted May 22, 2007. 809 1 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, thenit may be lost altogether from the population following sucha 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. Althoughselectionanddriftcansometimesbeinferredfromhistoricalchangesinthedistribution 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). Inboth cases, the selection coefficients estimated withthese methods willreflect 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 firstassumption 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 810 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 betweenthephylogeneticrateparametersandpopulationgenetical quantitiesisusuallyobscure. 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 (Kimura1964). This approach has been used by McVean and Vieira (2001) to estimate thestrengthof selection on so-called silentmutations(i.e, those whichdonotchange 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 individualwill 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 , we will define the substitution process to the common ancestor to be the stochastic t process (z : t ∈ R) and the common ancestor distribution to be the stationary distribution of t z . This process will be a good approximation to the substitution process along the branches t 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 pro- cesses by Jagers (1989, 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. Be- cause 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 inthissetting. Instead, wemuststudywhatGeorgii andBaake(2003) call the retrospective process, which characterizes the substitution process along lineages which survive forever. This process was also first described by Jagers (1989, 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. Forneutralmodels,thefactthatthesubstitutionprocessdecouplesfromthegenealogy of a sample can be used to deduce that the common ancestor process is simply the neutral 811 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. Thekeyideabehindthisconstruction istoassign particles tolevelsand then introduce look-down events which differ from (and replace) the usual neutral two-particle birth-deatheventsoftheMoranmodelintherequirementthatitisalwaystheparticleoccupying the higher level which dies and is then replaced by an offspring of the particle occupying the lowerlevel. Intheabsenceof selection, thecommon ancestor is theparticle occupyingthelowest 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 Kurtz1999). Furthermore,becausetheeffectofselection dependsonthefrequenciesofthetypes 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 mathematicaldifficultieswhichthiscreatescanbeovercomewiththesametechniquethatisused 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 Markovproperty. Onesuchenlargementistheancestralselection graphofKroneandNeuhauser (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 thattheaccessible particleconfigurationsconsistof thecommonancestor, whichcan beof 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 ,n ), t t where z is the type of the common ancestor and n is the number of virtual lineages. t t 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 Barton et al. (2004), 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 processofapopulationsegregatingtwoallelescanbeembeddedwithinabivariateprocess(z ,p ), t t wherep isthefrequencyattimetofoneofthetwoalleles. Wewillshowthatboththestationary t distributionandthegeneratorofthisprocesscanbeexpressedintermsofthesolutiontoasimple boundary value problem (Eq. 9) which determines the distribution of the type of the common ancestorconditionalonthefrequencyatwhichthattypeoccurswithinthepopulation. Incertain cases wecan solvethis problemexactly andobtain ananalytical 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. Theremainderofthearticleisstructuredasfollows. InSection2wedescribetheclassofdiffusion 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 whenthemutation ratesarevanishinglysmall. Sections 4and5areconcerned withapplicationsofthesemethodstoconcreteexamples, andwefirstconsidertheWright-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, andweinsteadresorttonumericallysolvingtheassociated boundaryvalueproblemtodetermine 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. 2 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 Barton et al. (2004) and Barton and Etheridge (2004). 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 µ and 1 µ , respectively, where both rates are expressed in units of events per N generations. Suppose, 2 in addition, that the relative fitnesses of P and Q are equal to 1+σ(p)/N and 1, respectively, wherep is thefrequencyof P. For technicalreasons, wewillassumethat theselection 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 denote the frequency t of P at time t and we measure time in units of N generations, then for sufficiently large N the time evolution of p can be approximated by a Wright-Fisher diffusion with generator t 1 Aφ(p) = p(1−p)φ′′(p)+ µ (1−p)−µ p+σ(p)p(1−p) φ′(p), (1) 1 2 2 (cid:0) (cid:1) where φ∈ C2([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, Theorem 2.1 of Ethier and Kurtz [(1986), Chapter 8] tells us that the set C∞([0,1]) of infinitely differentiable functions 0 813 with support contained in the interior of (0,1) is a core for A. Furthermore, provided that both mutation rates µ and µ are positive, then the diffusion corresponding to (1) has a unique 1 2 stationary measure π(dp) on [0,1], with density (Shiga 1981, Theorem 3.1; Ewens 2004, Section 4.5), p π(p) = Cp2µ1−1(1−p)2µ2−1exp 2 σ(q)dq , (2) (cid:18) Z0 (cid:19) whereC isanormalizingconstant. Unlessstatedotherwise(i.e.,whenweconsiderweakmutation limits inSection 3), wewillassumethroughoutthisarticlethatbothmutation ratesarepositive. 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 n˜ (t) and n˜ (t), respectively. Here, and throughout the article, we will use the tilde, 1 2 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 Barton et al. (2004), the generator G˜ of the structured coalescent process (n˜ (t),n˜ (t),p˜) can be written as 1 2 t n 1 G˜φ(n ,n ,p) = 1 [φ(n −1,n ,p)−φ(n ,n ,p)]+ (3) 1 2 1 2 1 2 2 p (cid:18) (cid:19)(cid:18) (cid:19) n 1 2 [φ(n ,n −1,p)−φ(n ,n ,p)]+ 1 2 1 2 2 1−p (cid:18) (cid:19)(cid:18) (cid:19) 1−p n µ [φ(n −1,n +1,p)−φ(n ,n ,p)]+ 1 1 1 2 1 2 p (cid:18) (cid:19) p n µ [φ(n +1,n −1,p)−φ(n ,n ,p)]+Aφ(n ,n ,p), 2 2 1 2 1 2 1 2 1−p (cid:18) (cid:19) whereforeach(n ,n ) ∈N×N,wehaveφ(n ,n ,·) ∈ C2([0,1]). Bartonetal.(2004) provethat 1 2 1 2 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 process p˜(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. Becauseweareonlyconcernedwithsubstitutionstosinglelineages,weneedonlyconsidersample configurations(n ,n )whichareeither(1,0)or(0,1),andsowecanreplacethetrivariateprocess 1 2 (n˜ (t),n˜ (t),p˜) with a bivariate process (z˜,p˜) taking values in the space E = ({1}×(0,1])∪ 1 2 t t ({2}×[0,1)),wherez˜ =1ifthelineageisoftypeP andz˜ =2ifitisoftypeQ. Wewillreferto t t (z˜,p˜) as the structured retrospective process to emphasize the fact that it describes evolution t t 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 814 the structured retrospective process can be written as 1−p G˜φ(1,p) = µ [φ(2,p)−φ(1,p)]+Aφ(1,p) 1 p (cid:18) (cid:19) p G˜φ(2,p) = µ [φ(1,p)−φ(2,p)]+Aφ(2,p), (4) 2 1−p (cid:18) (cid:19) for functionsφ ∈ D(G˜) ≡ C2(E) which are twice continuously differentiable on E and have com- c pact support. For future reference we note that D(G˜) is dense in the space Cˆ(E) of continuous functions on E vanishing at infinity and that D(G˜) is an algebra. The key step in proving the existenceanduniquenessofaMarkovprocess correspondingtothisgeneratoristoshowthatthe ancestral lineage is certain to jumpaway from a type before the frequency of that type vanishes, e.g., the ancestor will almost surely mutate from P to Q before the diffusion p˜ hits 0. This will t guarantee that the jump terms appearing in G˜, 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 Barton et al. (2004), which we restate below as Lemma 2.1. We also supply a new proof of this lemma to replace that given in Barton et al. (2004), which contains two errors (Etheridge 2005). One is that the variance σ(W ) appearing in the time s changeof theWright-Fisher diffusionneedsto be squared,so thatthe exponentα in theintegral displayedinEq.(16)ofthatpaperis2ratherthan1+ 1 . Thesecondisthatthedivergence 2(1−2µ2) 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 be the Wright-Fisher diffusion corresponding to the generator A shown in t (1). Then, for any real number R < ∞, τk 1 lim P ds >R = 1 p k→∞ (cid:26)Z0 (cid:18)ps(cid:19) (cid:27) τk′ 1 lim P ds >R = 1, p k→∞ (cid:26)Z0 (cid:18)1−ps(cid:19) (cid:27) where τk = inf{t > 0 :pt = 1/k} and τk′ = inf{t > 0: pt = 1−1/k}. Proof. For each positive integer k choose φ ∈ C(2)([0,1]) such that φ (p) = −ln(p) on [1/k,1] k k and observe that on this restricted set, 1 1 b(p) Aφ(p) = − − , 2p 2 p where b(p) = µ (1−p)−µ p+σ(p)p(1−p) is the infinitesimal drift coefficient in A. Then, for 1 2 815 each k > p−1, the stopped process 0 t∧τk M = φ (p )−φ (p )− Aφ (p )ds t∧τk k t∧τk k 0 k s Z0 1 t∧τk 1−2b(p ) 1 s = −ln(p )+ln(p )− ds+ (t∧τ ) t∧τk 0 2 p 2 k Z0 s is a continuous martingale with quadratic variation t∧τk t∧τk ds hMi = p (1−p )(φ′(p ))2ds = − (t∧τ ). t∧τk s s k s p k Z0 Z0 s In particular, on the set {τ < ∞}, we have k 1 τk 1−2b(p ) 1 s M = ln(k)+ln(p )− ds− τ τk 0 2 p 2 k Z0 s τk ds hMi = − τ , τk p k Z0 s which in turn implies that, for any R < ∞, the following three inequalities τ < R k hMi < R τk 1 M > ln(k)+ln(p )− +||b|| R τk 0 2 ∞ (cid:18) (cid:19) are satisfied on the set τk ds Ω = < R . R,k p (cid:26)Z0 s (cid:27) Now, because M is a continuous, one-dimensional martingale, there is an enlargement Ω′ ·∧τk of the probability space Ω on which the diffusion p is defined and there is also a standard t one-dimensional Brownian motion B , defined on Ω′, such that t M = B . t∧τk hMit∧τk [See Karatzas and Shreve (1991), Chapter 3, Theorem 4.6 and Problem 4.7.] Thus, in view of the conditions holding on Ω , we obtain the following bound R,k P{Ω } ≤ P supB > ln(k)+ln(p )−CR , R,k t 0 (t≤R ) where C = 1+||b|| is independent of k. The first half of the proposition then follows from the 2 ∞ 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 φ (p) = −ln(1−p) k on [0,1−1/k]. With Lemma 2.1 established, the next proposition is a special case of the existence and unique- ness results for structured coalescents proved in Barton et al. (2004). 816 Proposition 2.2. For any ν ∈ P(E), there exists a Markov process (z˜,p˜), which we call the t t structured retrospective process, which is the unique solution to the D [0,∞)-martingale problem E for (G˜,ν). Proof. Because the operator G˜ is a Feller generator when restricted to twice continuously dif- ferentiable functions on each of the sets E = {1}×[k−1,1] ∪ {2}×[0,1−k−1] , we can k show that a stopped version of the process exists on each of these sets and that this process is (cid:0) (cid:1) (cid:0) (cid:1) the unique solution of the corresponding stopped martingale problem. Then, using the Lemma 2.1 and noting that the diffusions p and p˜ are identical in distribution, we can show that t t the sequence of hitting times of the boundaries of the sets E is almost surely unbounded as k k → ∞. Consequently, Theorem 4.2 and Proposition 4.3 of Barton et al. (2004) imply the existence of a Markov process (z˜,p˜) defined on all of E which is the unique solution to the t t D [0,∞)-martingale problem for G˜. E Of course, as the name indicates, the process (z˜,p˜) describes the retrospective behavior of a t t 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 thefrequencyof that background;see Eq. (3). Furthermore, because lineages move betweengenetic backgroundsatrates whichareboundedbelowbythe(positive) mutationrates, 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 stationarydistributionofG˜ willalsobeastationarydistributionofthecommonancestorprocess. 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˜,p˜) which is defined for all times t ∈ R. However, because this lineage t t 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.) 2.1 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˜) = C2(E) is an algebra which is dense in Cˆ(E) and because the martingale problem for G˜ c is well-posed, any distribution π(z,dp) which satisfies the condition, 1 1 G˜φ(1,p)π(1,dp)+ G˜φ(2,p)π(2,dp) = 0 (5) Z0 Z0 817