- Email: [email protected]
Reversibility and the age of an allele II. Two-allele models, with selection and mutation
THEORETICAL
POPULATIOK
12, 179-196 (1977)
BIOI.OCY
Reversibility and the Age of an Allele II. Two-Allele Models, with Selection and Mutation G. A. WATTERSOS Deportment
of Mathematics,
Monosh University,
Received October
Victorio
3168, Australia
29, 1976
A Markov process with absorbing boundaries may be made recurrent by returning the process to the interior whenever a boundary is reached. The age of such a process may be defined as the length of time since the last return event. Examples drawn from two-allele genetic models are discussed, in which reversibility of the return process means that the age of an allele, whose present frequency in the population is known, has the same probability distribution as its future extinction time. Some discrete models are not reversible, yet if approximated by diffusion processes, the (approximate) age distribution is the same as the future extinction time distribution. Various results in the literature arc unified by this viewpoint.
1. INTR~DLJCTI~S
In Watterson (1976), it was shown that time reversibility of an infinitelymany-alleles model led to a fairly easy study of the ages of the alleles. In that model, all mutations produced new alleles, and all alleles were selectively neutral. By reversibility, we mean that given the present state of a stochastic process, the statistical properties of its future behavior are the same as those for its past history treated as a stochastic process with time running in the reverse direction. (See Feller, 1968, p. 414; Kendall, 1975.) In the present paper, we restrict the number of allelic types to two, but allow a selective differential between them as well as recurrent mutations. We are able to make progress with reversible models, and nonreversible models approximated by diffusions, beyond that achieved by Kimura and Ohta (1973), Maruyama (1974), Li (1975), and Levikson (1977), who did not exploit reversibility. We point out that the considerable literature concerning first passage times may be reinterpreted to solve the ageof-alleles problem. Sawyer (1977) has a somewhat different approach to the problem. Levikson (1977) considered a discrete Markov chain {X, . tl = 0, 1, 2,...} having state space (0, 1, 2 ,... ), whose transition matrix I’ = (pi,) is aperiodic, has 0 as an absorbing state certain to be reached ultimately, and is otherwise process {x,,}, which has the same irreducible. He introduced the “return” 179 Gpyright - All rights
8? 1977 by Academic Press. Inc. of reproduction in any form reserved.
ISSS
OIMO-5809
180
G. A. WATTERSON
transition matrix as {X%} except that whenever state 0 is reached, the chain instantaneously returns to state 1. He let A(t) denote the “age” of the return process at time t, namely the length of time since X last returned from state 0. Levikson proved that
under the assumption that the return process is positive recurrent, having finite mean time between successive visits to state 0 and having a limiting stationary distribution. Pakes (1977) has discussed the generalization of (1) to cases when the return process is not positive recurrent. Levikson (1977) postulated that the corresponding formula for discrete states, continuous time processes (e.g., birth and death processes) is given in terms of the absorbing process transition probabilities by Is;_n;t Pr(A(s) E (t, t + dt) ! X(s) = j) = pJt>
dt/JoK p&t)
dt.
(2)
For diffusion processes for which the absorbing process has transition density p(x, s, y) dy = Pr(X(s) E (y, y + dy) / X(0) = x) and for which, whenever state 0 is reached a return to state x* occurs instantaneously, Levikson (1977) postulated that lii
Pr(A(s) E (t, t + dt) ’ k(s) = y) == p(x*, t, y) dt/p
p(x*,
t, y) dt.
(3)
In this section, we give two examples of the application of (1) and (2) respectively. We concentrate our attention in Section 2 on models for which (1) is strictly appropriate, but for which a diffusion approximation, as in (3) is employed. In Section 3, we consider diffusion approximations when there are two return boundaries. Examples of the diffusion approximations are given in Section 4. Moran (1958, 1959), introduced a population genetic model in which there are two possible alleles, and in which at each unit time one gene dies and is replaced. He allowed for the possibility of mutation, and for viability or fertility differences between the two alleles. The number of genes belonging to a particular allele evolved as a random walk process, changing by at most one gene per unit time. When the process has settled down to statistical stationarity, it may be proved to be reversible. Thus although the result (1) could be used for this process, it is equally possible to discuss the age distribution by equating it to the extinction time distribution. The time since an allele was last reintroduced to the population by mutation has the same distribution, conditional on the present number of genes of that allele, as the time from the present until the next occasion (in the future), when it will again be lost. We do not go into the details
AGE
OF
AN
ALLELE
181
of this example, because they are similar in some respects to the example presented in Watterson (1976), and in order to obtain numerical results, one would probably use a diffusion approximation similar to those discussed below in Sections 2, 3, and 4. There is perhaps need for a rigorous proof of (2), especially in the case of transient or null-recurrent return processes. Dr. A. Pakes (personal communication) has established some results in this direction. Levikson (1977) illustrated the use of (2) on a birth and death process, in which the infinitesimal birth and death rates are hn and pn, respectively, whenever the population is of size n. When 0 < h < p, Levikson evaluated the right side of (2) as dt . [I - CX(t)][l - B(t)] /3(t>j-‘j$/R’,
(4)
where a(t) = CL(exp@) - 1) A expbt) - P ’
p(t) = YexpW - 1) A expW - CL ’
Y==h-#LL.
Our purpose is simply to point out that the age density (4) may be obtained as being the density of the extinction time, the time to first passageto state 0 fromj as initial state. We have for j = 1,2, 3,...,
Pdt)lLm Pdt> dt = (d/dt)[P~,dt)l, t 3 0,
(5)
which may be verified from (4) above and using
The reason for the result (5) is that the return process (ignoring any instantaneous transitions from 1 to 0 and return) is a birth and death process on the positive integers, with stationary distribution being the logarithmic distribution mi
=
-W)i/(i
141
-
i = 1) 2, 3)...)
(VP.)),
and the process is reversible because of the connection between the stationary distribution and the infinitesimal parameters: m,Ai = mi+lp(i + 1)
for
i = 1, 2, 3,...;
see Kendall (1975, (12)) for the general theory. Even in the critical case when /I = p, whether or not (2) is itself valid, the
G. A. WATTERSOX
182
result (5) may be easily verified, using the birth and death transition p&t) forj
= (At)j-i/(1
+ xqi+i
and
probabilities
pj,, =: [At/(1 + ht)]j,
= 1,2,3 ,... .
2. DIFFUSION APPROXIMATIONS; ONE RETURN BOUNDARY Some return Markov chains, for which (1) is applicable exactly, are not strictly reversible but yet have the property that their stationary age distribution (1) is approximately the extinction time distribution pjyd - p$;‘) for the nonreturn process. We establish this result by considering only those processes which have good diffusion approximations. Such processes are particularly common in population genetics. Suppose therefore that, after suitable resealing of the state space and the time variables, the original absorbing process is approximated by a diffusion with drift and diffusion parameters denoted by p(x) and us(x), respectively, and with 0 as an accessible, absorbing state (see Goel and Richter-Dyn, 1974, p. 51, for terminology). Corresponding to the discrete process being returned to state 1 whenever state 0 is reached we suppose that the approximating “return diffusion” is returned to state x*, interior to the diffusion state space. Whether or not (3) may be rigorously proved for the age of a return diffusion process, it is certainly going to yield an approximation to the exact age distribution (1) in the cases under consideration. It6 and McKean (1965, p. 149, (2)) exhibit the diffusion transition densityp(x*, t, y) as the product of two functions, one of which is a function of y only and the other is symmetrical in x* and y. (I thank Frank Norman for the reference to the proof of this fact. He has shown (personal communication) that it implies reversibility of diffusion processes under certain sufficient conditions.) Defining the “natural scale”
G(Y) = exp ) -2 ly [libel and the “speed measure”
dr 1
(6)
density
4~) = [u2(y) G(~)l-l, we have p(x*, t> Y> = m(y) 4(x*, t> Y), where 4(x*, t, y) is a symmetrical
function
4(x*, t, Y)
in x* and y:
= 4(Y, 4 x*1.
(7)
AGE
The right-hand side of (3) (omitting symmetrical in x* and y also:
OF AN
183
ALLELE
the differential
dt) may now be seen to be
Y> dt=4~)4(x*> t,Y)/ jam 4~)dx*, t,Y)dt =m(x) dy,t,x*,/lx dt 04.4dv,t,"x*1 =P(Y,t, x*)/l= P(Y>t, x*>dt.
I+*, t,y)/jomP(x*>t>
Moreover, it is usually the case in population genetics that the return state xx is arbitrarily close to 0, being of the form x* = 1/(2N) as N + co, so that the right side of (3) is approximately
P(Y,t, O-t)/jomp(y, t, O+)dt.
(8)
This approximating age density is also the density of the extinction time, t, for the absorbing diSfusion process to jirst reach the absorbing state 0 from y as initial state (cf. Crow and Kimura, 1970, 8.3.21). It is this observation which is the main point of the paper. Most authors have studied the mean age, rather than the age distribution. The mean of the approximating age distribution (3) is, in Levikson’s notation,
= B(x*, Y)/G(x*Y)
say,
where B(x*, y) = lrn tP(x*, t,y) dt and G(x*, y) = jam P(x*, t, Maruyama
Y) dt.
(1974) showed that
Y)+,4x)(W) B(x, Y)=- G(x, Y),
W4(~zl~x2) B(x> and that
W4(W~z)
‘7x,
Y)+&)(VW G(x, y)=--s(x - y).
184 Levikson
G. A. WATTERSON
(1977) noted that
A(x*,y)
= joffi G(x*, 5) G(t, Y) G’G(x*,
Y).
We write
for the (approximating) mean age when the diffusion is returned to near the 0 state, conditional on a diffusion state y now. In view of the limiting relationship between (3) and (8), we find that M(Y)
= jm MY,
t, Of)
0
qj-=
0
p(y, t, O+) dt,
which is the mean extinction time for the diffusion conditional on state y now. This is given explicitly in terms of the approximating diffusion parameters CL(.) and CT’(.), by the standard formula for extinction times:
M(y) = 2 j’ G(x) j” [u”(z) G(z)]-l d.z dx, x 0
(9)
where G(z) is as in (6); see, for instance, Goel and Richter-Dyn (1974, Table 3.6). We have written co in (9) for the upper bound of the state space of the approximating diffusion, although this may well be replaced by a finite value in some applications. The mean extinction time and the mean age may be infinite in particular cases. An example is given in Section 4.1 below. Sometimes, population geneticists are interested in processes approximated by a diffusion for which there are two absorbing, accessible, boundaries 0 and 1, but conditional on boundary 1 not being reached. If the unconditional (absorbing) diffusion has parameters p(.), up(.), then the diffusion conditioned on ultimate absorption at 0 (rather than I) has parameters /Lo(x) = p(x) i- U2(X) @o’($@&), (10)
uo‘yx) = c?(x), where Go(x) is the probability and is given by
of ultimate
absorption
in 0 from the initial state X,
(11)
@o(x) x j1 G(Y) 41 jol G(Y) dy,
I
where G(.) is as in (6). (See Ewens, 1973; Narain, 1974.) The return process corresponding to the conditioned
process
(10) has 1
185
AGE OF AN ALLELE
as an inaccessible boundary. The mean age, and the mean first extinction time, is then as given by (5) but calculated with the conditional parameters (10): M,(y)
G(x)
= 2 j’
”
G,(x) j1 [uo2(z) G,(z)]-l dz dx, z
= exp 1-2 .\” iki~)I~~~(r)l
dyi
= G(M@,(412Integrating (12) by parts and also using (13) sh ows that the mean conditional age can be written
M,,(y)
= 2 1j’ @&)[l
0
- @&)][G(s)
U”(X)]-’ dx (14)
- j’
Y
[l - @oW%(r)l
@,MW
a2CW1 dx/ * j1 G(x) dx, 0
which can be shown to be equivalent to Li (1975, Eq. (lo)), which was derived from the work of Maruyama (1974). Ag ain, our main point is that (14) is both a mean age and a mean extinction time. In fact, (14) is equivalent to the mean extinction time formula studied by Kimura and Ohta (1969b, (1)). Admittedly, their emphasis was on takingy = 1/2Nas initial state (a new mutant’s frequency) whereas we are interested in larger values of y. The equivalence in the neutral allele case was noted by Kimura and Ohta (1973, p. 208). Examples of the use of (14) are given in Section 4.2 for genetic diffusions subject to mutation and selection.
3. PROCESSESWITH Two
RETURN BOUNDARIES
Levikson (1977) h as considered Markov chains having 0 and N as absorbing states. He modifies the process by the assumption that whenever 0 is reached, an instantaneous return to 1 occurs, and whenever N is reached, an instantaneous return to N - 1 occurs. We now consider the diffusion approximations to Levikson’s results. Consider a diffusion approximation scaled so that 0 and 1 are absorbing boundaries and a return process so that state 0 is instantaneously followed by state x* (x* + 0) and state 1 is instantaneously followed by state x** (x ** G 1). The formula corresponding to Levikson’s discrete age distribution is the stationary age density for the time, A(t), that has elapsed since the process was last returned from a boundary. Kimura and Ohta (1973, p. 208) called this
186
G.
A.
WATTERSON
time the “age of the polymorphism.” density is
Using the notation of Section 2, this
limPr(A(s) E (t, t + dt) X(s) = y)!dt S’XI
(1 - @&*))p(x**,
t, Y) + @,(x**)P(x*~
6 Y)
t > 0.
(1 - @0(X*)) s,- p(x**, t, y) dt + @I@**) s,- p(x*, t, Y) dt ’
Levikson’s formula for the probability that the last return, before time t, from a boundary was from boundary 1 becomes fim Pr(last return from 1 : T(t) = y)
(1 - @,,(x*))&(x**,
t, y) dt
(1 - C&(X*)) JR p(x**, t, y) dt + @a(~**) J,- p(x*, t, Y) dt To see the role played by reversibility in these expressions, we may deduce from (6) and (7) that
I% 4 Y)
= P(Y, t, 4 dYM4,
= P(Y, t, 4 44
G(~/[uYY)
G(Y)!
We find
v-2 Pr(A(s) E (t, t + dt) ~T(s) = y)/dt = ((1 - Q&X*)) 0*(x**) G(x**)p(y, + @,,(x**) u”(x*) G(x*)P(Y,
t, .**)
t, x*>>
x /(I - @a@*)) u*(x**) G(x**) lmp(y,
t, .**)
dt
+ Q&v**) u”(x*) G(x*) sump(y, t, x*) dtl-’ and iir Pr(last return from 1 ~X(t) = y) = /(I - CD&*)) u2(x**) G(x**) JT: p(y, t, .**) dt/ x /(I - Q&e*)) ,2(x**) G(x**) j-)(y, + Qo(x**) u”(x*) G(x*) j-)$y,
t, x**) dt
t, x*) dt/-‘.
187
AGE OF AN ALLELE
But now, if the return states are allowed to approach their respective boundaries, x* --f 0 and x** -+ I, we find using (6) and (11) that [I - tDo(x*)] u2(x**) G(x**) u’(x*) G(.r*) Bo(x*)
-
5 a0 ’
where a, =- -F Thus finally,
a, ~ - da2(4 dx i . x=0
and
1 , s=1
taking limits as x* + 0, x** ---f 1, and s -+ GO, hmPr(A(s)
E (t, t + dt) ( x(s)
= y),‘dt
a,P(y, t, 1-j + a,P(y, 6 O+> al J;P(r, t, 1-I dt + a0S,-P(Y, t, Of) dt
(15)
and Lm Pr(last return from 1 1X(t) = y)
al.frp(y, t,I-)dt aI J;P(y,
4 1-I
dt + a,
J"p(y, 0
t, 0+) dt
(16)
= 1 - Do(y). The result (15) has the interpretation of being the probability density of the time (in the future) at which one or other boundary, 0 or 1, is first reached from initial state y (see, e.g., Crow and Kimura, 1970, pp. 378,379). Thus the past age density equals the future absorption time density, when returns are made to states infinitesimally close to the boundaries and when stationarity is reached. Similarly, the probability of a last return from 1, (16), has the alternative interpretation of being the probability that absorption will first occur (in the future) into boundary 1 rather than boundary 0, given initial state y (see (11) above). These results explain why certain formulas in the literature of population genetics, which were derived as properties of ages of alleles, can also be found in papers concerning the extinction of alleles. We mention some examples in Section 4.3 below. 4. EXAMPLES 4.1. Diffusion Approximation
for a Branching Process
Levikson (1977) considered a branching tion with infinitesimal parameters p(x) = (m -
1)X,
process and its diffusion c?(x) = bx.
approxima-
188
G. A. WATTERSOlX
Here, m is the mean number of offspring per individual and b is the variance. Cox and Miller (1965, Sect. 5.11) show that p(x*, t, y) is a mixture of gamma densities with one term missing: $@*,
t, y)
=
f
[e~"*"iB(.~*~i~)"!n!]
py'~-le-ylfl,(n
-
I)!,
71=1
where (y =
e’m-lM,
/3 =- b(e(wL-l)t -
1)/2(m -
If m = 1, we have 01:= 1 and /3 =- bt/2. The probability time t, from x* as initial state, is the missing term Pr(X(t)
1). of absorption
= 0 / X(0) = .Y*) = e-r**16.
We assume from now on that m C< 1, so that ultimate the nonreturn process. But the transition density may also be written as
I+*, f, Y)
=
m(r)
4(.~*,
at 0 by
(17)
extinction
is certain for
t, Y),
where m(y)
I_
2$(n"-l)"/bj(by)
and q(x*, t,y)
:== $b i
e-‘“‘+““.‘6(~*y)n
anp-?n;n! (n -
I)!,
n=1
which is symmetrical written as
in x* and y. The approximating
P(Y, t, .x*1
age density (3) may be
dY1 t, x*1
J” P(Y, t, .x*) dt = s,- q(y, t, x*) dt 0
Letting
x* -
0 we find for the limiting
expression
as in (8)
But this density for the age, t, is exactly the extinction directly from (17): (djdt) Pr(X(t)
time density calculated
= 0 1X(0) = y) = ibe -~~@yajiF,
t > 0.
AGE
OF
AN
189
ALLELE
This means that (~-l = e(l-m)t has an extreme value (Type 1) distribution, (see Johnson and Katz, 1970, Chap. 21). The mean of the absorption time distribution, and the mean age, for a given state y now, is M(y) = IL (1 - e-vu/s) dt = (;I(1
- an)) I’
(1 - e-2(1-m)zlb) x-‘( 1 + x/y)-’
dx
(18)
using the substitution x = y~(l-“)~/(l - e(l-l”)t), assuming m < 1. This can also be derived from (9). The special case when m = 1 and b = 2 most interested Levikson. We find that the age density, given state y and as x* -+ 0, is e-gltyt-2
t
2
(19)
0.
This distribution is an extreme value (Type 2 or Cauchy) distribution (see Johnson and Kotz, 1970, Chap. 21). It does not have a finite mean, as can be seen directly, or by letting m t 1 in (18). But the modal age is t = y/2, which increases as y increases. These remarks agree with observations made by Levikson (1977) for a corresponding discrete branching process. For our diffusion, the reciprocal of the age has an exponential distribution, with mean 1/y, as x* + 0. As a check on the accuracy of the diffusion approximation (19) to the actual age distribution as tabulated by Levikson (1977) for the branching process with geometric offspring distribution pn = ($)n+l, n = 0, 1, 2, 3,..., we make the comparison in the form Pr(age E n / 2 = j) + In+' e-i/tjt-2 dt n-t = e-j/'n+:' . - e-j/(s-:,),
112 ,-jjtj-2 : 10
z
dt
=
e-2i
for
n =T 1, 2,...,
for
n = 0.
In Table I, we show the left side, as tabulated by Levikson (1977), with the approximating right side in parentheses. Clearly, the diffusion approximation is accurate for large values of j, and is still moderately accurate for small j but large values of n. 4.2. Genetic Da&ion,
with Mutation
and Selection
The most common population genetic diffusion models having 0 as an absorbing state have parameters of the form
190
G. A. WATTERSON
TABLE .4ge Probability
Distributions
I
for Branching
Process (Diffusion
Approximation)
Present frequency
Age
1
2 -~___
3
4
5
__-
6 ~~~ --
8
9
.608 (.135)
(.O! 8)
C.0002)
$0)
(.c&
co&
(.~O)
(.OiO)
(.OiO)
.I51 (.378)
.282 (.245)
.193 (.133)
.126 (.069)
.078 (.036)
.047 (.018)
.027 (.009)
.016 (.005)
(.ow
2
.067 (.157)
.167 (.186)
.152 (.166)
,132 (.132)
.lOl (.lOO)
,088 (.072)
.068 (.051)
.052 (.036)
.039 (.025)
3
,038 (.081)
.105 (.115)
,109 (.123)
.106 (.117)
.099 (.104)
.089 (.089)
,078 (.075)
.067 (.061)
.056 (.049)
4
.024 (.049)
.072 (.076)
,079 (.089)
.082 (.092)
.082 (.090)
.079 (.084)
.073 (.076)
,067 (.067)
,060 (.059)
5
.017 (.033)
.0.52 (.054)
.060 (.066)
.065 (.072)
.067 (.074)
.067 (.072)
.065 (.069)
.062 (.064)
.058 C.059)
6
.012 (.024)
.040 (.040)
.046 (.051)
.052 (.057)
.056 (.060)
.057 (.061)
.057 (.061)
.055 (.059)
.054 (.056)
7
.009 (.018)
.030 (.031)
.037 (.040)
.042 (.046)
.046 (.050)
.048 (.052)
.049 (.053)
.049 (.052)
,048 (.05 1)
8
.008 (.014)
.025 (.024)
.030 (.032)
.035 (.038)
.039 (0.42)
.041 (.044)
.043
.043
C.046) (.046)
.043 (.046)
9
.006 (.Oll)
,020 (.020)
.025 (.027)
,029 (.032)
.033 (.035)
.035 (.038)
.037 (.040)
,038 (.041)
.039 (.041)
10
305 (.009)
.017 (.016)
.021 (.022)
.025 (.027)
.028 (.030)
.031 (.033)
.033 (.035)
.034 (.036)
035 (.037)
15
.002 (.004)
.008
(.C’W
.OlO (.Oll)
.013 (.014)
.015 (.016)
.017 (.018)
.019 (.020)
,020 (.021)
.021 (.022)
.0014 (.002)
,005 (.@w
.006
(.oW
.008 (.008)
.009 (.OlO)
.Oll (.Oll)
.012 (.012)
,013 (.013)
.014 (.014)
.0009
.004 C.004)
.005 (035)
.006 (.007)
.007 (.008)
.008
.009
(.ocw
.003 (.003)
C.0’38) (.009)
.OlO (.OlO)
.0006 (.OOl)
,002 (.002)
003 C.003)
.004 (.004)
.005 (.005)
.005 C.005)
.006 (.006)
.007 (.007)
.007 (.007)
0
1
20
25
30
.009
AGE OF AN
191
ALLELE
p.(x) = 2s{h + (1 - 2h) X} X( I - X) - UX, (20)
U”(X) = X(1 - X)/2Ne ,
where 0 < x < 1, N, is the effective (diploid) population size, u is a mutation rate (0 < u < I), and s and h are selection and dominance parameters. Maruyama (1974) writes F = 1 - 4N,u, and S = 4N,s, so that p(x) = (2S[h + (1 - 2h)x](l - x) - (1 - F)}x/4Ne
.
We find from (6) that G(x) = (1 - ~)~-l exp{-S[2hx
+ (1 - 2h)x’]}.
(21)
We consider both the case when the age of an allele is calculated allowing for quasi-fixation occasions, and the case when the allele frequency is conditioned so as not to achieve the boundary value 1 during the lifetime of the allele. When u > 0 and hence F < 1, the back mutation is such as to make the boundary 1 nonabsorbing. The mean age of an allele of frequency y is, by (9), (2% and @I), M(y) = 4N, 1’ (1 - .~)~--lexp{ -S[2hx 0
+ (1 - 2h) a?])
1
X
sz
z-‘( I - z)-” exp(S[2hz + (1 - 2h) z2] dz dx.
(22)
The selectively neutral case provides a simpler expression; when S = 0 we have, integrating by parts, M(y)
= (4Ne.lF) 1-In y + J1 z-‘[( 1 - x)-’
- l] dz
0
- (1 - y)” j: z-l(l
- z+” dzj.
(23)
Priority for this result as an e&z&on time formula belongs to Ewens (1964, (2.6)) and as an average age formula, to Maruyama (1974, (12)). Two special cases deserve closer attention; when F = - 1, (23) reduces to WY)
= -4NAy
lnr)/(l
- Y),
(24)
while if F = 0, (24) is replaced by M(y) = -4N,
1’ (In z)/(l - z) dz, 0
(25)
192
C. A. WATTERSON
which approaches 4N$/6 as y T 1. Numerical values of M(y), as in (23) (24), (25) have been published by Maruyama (1974, Table 2) but these are sometimes incorrect in the third significant figure, no doubt due to numerical quadratures. If we wish to discuss the mean age (or mean extinction time) conditional on no intervening quasi-fixation occurrences, then if F < 0, the previous formulas (22)-(25) apply. We have Jfo(.?9 = M(Y)
if
F < 0,
because the boundary N = 1 is inaccessible in these cases. Hence column 4Nu = 1 in Maruyama’s Table 1 should agree with that column in his Table 2 and with (25) above. If, however, F > 0 so that x = 1 is an accessible boundary, the diffusion conditional on boundary 1 not being reached differs from the unconditional one. Li (1975) has tabulated the mean age M,,(y) as given by (14), for various s and h combinations, but in the absence of back mutation (u = 0, F = 1). Note that Li’s s corresponds to our -2s. There is some simplification of (14) possible in the case of no mutation and no dominance. With F = I, h -= Q , we have the conditional diffusion parameters corresponding to (20) as given by (10): pa(x) = (1/4N,) x( 1 -
.x) S(e+ + epss)/(e@ -
ecS”)
(26)
and uo2(x) =z (1/2N,) x(1 - x), which lead to the conditional
mean age, from (14) of
~ 1 (1 - e-s(l-s))(e-sY - e-S%)dr 4N, I. s7J e-S) x(17 Iy S(e-
(27)
This is identical with (13) in Marugama (1974), who derived it using a different argument. Kimura and Ohta (1969a, (17)) g ave essentially this result as an extinction mean time. If S = 0, we find the limiting (neutral allele) conditional age as MO(y) == --4N,y
lny/(l
- Y).
(28)
This is identical with (24); the reason will be discussed later. Maruyama (1974, Table 3) has tabulated values of (27) f or various S values. His column / 4Ns 1 = 0 should agree with (28), and hence with his column 4Nu = 0 in Table 1 and column 4Nu = 2 in Table 2; discrepancies must be due to numerical integration procedures.
AGE OF AN
193
ALLELE
As Maruyama noted, M,,(y) in (27) is an even function of S; that is, for a given allele frequency y, the conditional mean age is independent of whether selection is in favor of, or against, the allele. This appears at first sight to be paradoxical. The mathematical explanation lies in the fact that the conditional diffusion (26) is completely independent of the sign of S. The statisticaE explanation of the paradox is the following. Let X(t), 0 < t < T, be the realization of the unconditional diffusion with parameters p(x) =J Sx(I - x)/4Ne )
u’(x) = x(1 - x)/2Ne .
The stopping time, T, is that time when one or other boundary x = 0 or x = 1 is first reached. The likelihood of that realization, with respect to the measure of the diffusion having S = 0, is
= exp [$S[X( T) - X(O)] - & N:lS2 L’ X(t)( 1 - X(t)) dt/. This likelihood follows from the work of Brown and Hewitt (1975). The stochastic integral si X(t)( 1 - X(t)) dt is associated with the parameter F, and clearly contains information about the magnitude, but not the sign of S. The maximum likelihood estimator of S is, incidentally, L? = 4N,[X(T)
- X(O)]//rX(t)(l
- X(t)) dt.
0
By discussing the diffusion (26), which is conditional on the boundary 1 not being reached, we have conditioned on the sufficient statistic for the sign of S, X(T) - X(O), taking only the negative value -X(O), thus removing the dependence of the realization on the sign of S. In the case when there is mutation but there is no selection, then provided F > 0 so that boundary 1 is accessible, we have the conditional diffusion parameters
PO(X)= -4
+ F)/4N, ,
q,‘(x) = x(1 - x)/2Ne ,
which are similar to the unconditional diffusion parameters. Then the average age and the average extinction time may be derived from (14) to be
dx+
(1
Jy)'
i y
’ (’ ; X)F dsj, (29)
a result equivalent to Maruyama (1974, (10)) and Kimura and Ohta (1973,
194
G. A. WATTERSON
(13a)). Of course, (29) is just (23) with F replaced by -F, and so, in a sense, can be attributed to Ewens (1964) as an extinction time formula. Letting F + 1 yields the conditional mean age, as given by (28) in the absence of mutation and selection, and as given by (24) for the unconditional mean age with F = - 1. We have concentrated our attention on the mean age, either conditional or unconditional on no fixation. In many of the above examples, explicit formulas are known for the transition probability densities entering into the age density (8). Indeed, the time to extinction density (and hence the age density) for an allele of frequency y now, in the absence of mutation and selection and conditional on no fixation is ytgf’(2i+
l)(-l)i-1-lF(l
-i,i+2,2;
1 -y)hie-W,
t>O
where hi = i(i + 1)/4N, (see Kimura, 1970). It is this density which has mean i%!,(y) as in (28). M oreover, even though some mutation and selection cases have complicated transition densities, the approximation obtained by Voronka and Keller (1975, Sect. 7) may well be suitable, at least for relatively low ages. The higher order moments in some of the cases above can also be obtained explicitly. 4.3. Genetic Dzjfusion with Two Return Boundaries
We return to the diffusion whose parameters are given by (20), but specialized by putting u = 0: p(x) = 2s{h + (1 - 2h)x) x(1 - x),
and G(x) = x(1 - %)/2N, . Both boundaries x = 0 and x = 1 are now absorbing and accessible. In Section 4.2 we discussed, inter alia, the case when the process was conditional on 1 not being reached and the boundary 0 was a return boundary. We now consider the case, as discussed generally in Section 3, when both 0 and 1 are return boundaries. In the absence of selection, the mean age of a polymorphism was given by Kimura and Ohta (1973, p. 208) as -4N,[y
lny + (1 - A 141 - 391
when the return process is in state y. This is the mean of the distribution (15), and was first found by Watterson (1962) as the mean absorption time. The
AGE
OF
AN
195
ALLELE
probability that the last return was from state 1 is a special case of (16); Kimura and Ohta (1973) state this probability is y in the neutral alleles case. Of course, this is the well known result for the fixation probability of a neutral allele of frequency y. The full age density (15) may be obtained from Kimura’s (1955) work on fixation time densities: lim Pr(A(s) E (t, +dt) =y(l
-y)
f
1X(s) = y)/dt h,(2i+
l)(-l)i”[F(l
--,i
+2,2,y)
i=l
+ F(l
- i, i + 2, 2, 1 - y)] e&it,
where Ai = i(i + 1)/4Ne . The age density (and the absorption density) has a complicated form in the selection case, even when h = $ corresponding to no dominance effects, or when h = 0 or 1 corresponding to complete dominance of one or other allele; see Crow and Kimura (1970, 396-414). The approximation given by Voronka and Keller (1975) may prove useful here. The mean age may be equated to the mean absorption time from state y (see the Correction to W’atterson (1962)). The probability that the age corresponds to a return from boundary 1, (16), reduces to the fixation probability given by Kimura (1957, p. 896):
1-
@o(y) = 1 G(x) dx/(
where, for the process here, G(X) may be obtained G(x) = exp(-S[2hx
+ (1 -
G(x) dx, from (21) as 2h)x2]},
with S = 4N,s. ACKNOWLEDGMENTS I thank B. Brown paper.
and F. Norman
for helpful
suggestions
in the preparation
of this
REFERENCES B. M., AND HEWITT, J. I. 1975. Asymptotic likelihood theory cesses, J. Appl. Prob. 12, 228-238. Cox, D. R., AND MILLER, H. D. 1965. The Theory of Stochastic New York. CROW, J., AND KIMURA, M. 1970. An Introduction to Population Harper & Row, New York. EWENS, W. J. 1964. The pseudo-transient distribution and its uses in Probability 1, 141-156. BROWN,
^__653172124
for diffusion
pro-
Processes, Wiley, Genetics
Theory,
genetics, J. Appl.
196
G. A. WATTERSON
EWENS, W. J. 1973. Conditional diffusion processes in population genetics, Theor. Popul. Biol. 4, 21-30. FELLER, W. 1968. An Introduction to Probability Theory and Its Applications, Vol. I, 3rd Ed., Wiley, New York. GOEL, N. S., AND RICHTER-DYN, N. 1974. Stochastic Models in Biology, Academic Press, New York. 1~6, K., AND MCKEAN, H. 1965. Diffusion processes and their sample paths, Academic Press, New York and Springer-Verlag, Berlin. JOHNSON, N. L., AND KOTZ, S. 1970. Continuous Univariate Distributions I. Houghton Mifflin, Boston. KENDALL, D. G. 1975. Some problems in mathematical genealogy, in “Perspectives in Probability and Statistics: Papers in Honour of M. S. Bartlett” (J. Gani., Ed.), pp. 325-345, Applied Probability Trust and Academic Press, London. KIMKJRA, M. 1955. Solution of a process of random genetic drift with a continuous model, Proc. Nat. Acad. Sci. U.S.A. 41, 144-I 50. KIILIURA, M. 1957. Some problems of stochastic processes in genetics, Ann. Math. Statist. 28, 882-901. KIMURA, M. 1970. The length of time required for a selectively neutral mutant to reach fixation through random frequency drift in a finite population, &net. Res. Camb. 15 131-133. KIMURA, M., AND OHTA, T. 1969a. The average number of generations until fixation of a mutant gene in a finite population, Genetics 61, 763-771. KIMURA, M., OHTA, T. 1969b. The average number of generations until extinction of an individual mutant gene in a finite population, Genetics 63, 701-709. KIMURA, M., AND OHTA, T. 1973. The age of a neutral mutant persisting in a finite population, Genetics 75, 199-212. LEVIKSON, B. 1977. The age distribution of Markov processes. J. ilppl. Probability, to appear. LI, W. H. 1975. The first arrival time and mean age of a deleterious mutant gene in a finite population, Amer. J. Hum. Genet. 21, 274-286. MARUYAMA, T. 1974. The age of an allele in a finite population, Genet. Res. Camb. 23, 137-143. MORAN, P. A. P. 1958. Random processes in genetics, Proc. Camb. Phil. Sot. 54, 60-11. MORAN, P. A. P. 1959. The survival of a mutant gene under selection, J. Austral. Math. Sot. 1, 121-126. NARAIN, P. 1974. The conditional diffusion equation and its use in population genetics, J. Roy. Stat. Sot. B 36, 258-266. PAKES, A. G. 1977. On the age distribution of a Markov chain, J. Appl. Probability, to appear. SAWYER,S. 1977. On the past history of an allele now known to have frequencyp, to appear. VORONKA, R., AND KELLER, J. B. 1975. Asymptotic analysis of stochastic models in population genetics, Math. Biosci. 25, 331-362. WATTERSON, G. A. 1962. Some theoretical aspects of diffusion theory in population genetics, Anti. M&h. Statist. 33, 939-957. (1963) Correction 34, 352. WATTERSON, G. A. 1976. Reversibility and the age of an allele. I. Moran’s infinitely many neutral alleles model, Theor. Pop. Biol. 10, 239-253.
POPULATIOK
12, 179-196 (1977)
BIOI.OCY
Reversibility and the Age of an Allele II. Two-Allele Models, with Selection and Mutation G. A. WATTERSOS Deportment
of Mathematics,
Monosh University,
Received October
Victorio
3168, Australia
29, 1976
A Markov process with absorbing boundaries may be made recurrent by returning the process to the interior whenever a boundary is reached. The age of such a process may be defined as the length of time since the last return event. Examples drawn from two-allele genetic models are discussed, in which reversibility of the return process means that the age of an allele, whose present frequency in the population is known, has the same probability distribution as its future extinction time. Some discrete models are not reversible, yet if approximated by diffusion processes, the (approximate) age distribution is the same as the future extinction time distribution. Various results in the literature arc unified by this viewpoint.
1. INTR~DLJCTI~S
In Watterson (1976), it was shown that time reversibility of an infinitelymany-alleles model led to a fairly easy study of the ages of the alleles. In that model, all mutations produced new alleles, and all alleles were selectively neutral. By reversibility, we mean that given the present state of a stochastic process, the statistical properties of its future behavior are the same as those for its past history treated as a stochastic process with time running in the reverse direction. (See Feller, 1968, p. 414; Kendall, 1975.) In the present paper, we restrict the number of allelic types to two, but allow a selective differential between them as well as recurrent mutations. We are able to make progress with reversible models, and nonreversible models approximated by diffusions, beyond that achieved by Kimura and Ohta (1973), Maruyama (1974), Li (1975), and Levikson (1977), who did not exploit reversibility. We point out that the considerable literature concerning first passage times may be reinterpreted to solve the ageof-alleles problem. Sawyer (1977) has a somewhat different approach to the problem. Levikson (1977) considered a discrete Markov chain {X, . tl = 0, 1, 2,...} having state space (0, 1, 2 ,... ), whose transition matrix I’ = (pi,) is aperiodic, has 0 as an absorbing state certain to be reached ultimately, and is otherwise process {x,,}, which has the same irreducible. He introduced the “return” 179 Gpyright - All rights
8? 1977 by Academic Press. Inc. of reproduction in any form reserved.
ISSS
OIMO-5809
180
G. A. WATTERSON
transition matrix as {X%} except that whenever state 0 is reached, the chain instantaneously returns to state 1. He let A(t) denote the “age” of the return process at time t, namely the length of time since X last returned from state 0. Levikson proved that
under the assumption that the return process is positive recurrent, having finite mean time between successive visits to state 0 and having a limiting stationary distribution. Pakes (1977) has discussed the generalization of (1) to cases when the return process is not positive recurrent. Levikson (1977) postulated that the corresponding formula for discrete states, continuous time processes (e.g., birth and death processes) is given in terms of the absorbing process transition probabilities by Is;_n;t Pr(A(s) E (t, t + dt) ! X(s) = j) = pJt>
dt/JoK p&t)
dt.
(2)
For diffusion processes for which the absorbing process has transition density p(x, s, y) dy = Pr(X(s) E (y, y + dy) / X(0) = x) and for which, whenever state 0 is reached a return to state x* occurs instantaneously, Levikson (1977) postulated that lii
Pr(A(s) E (t, t + dt) ’ k(s) = y) == p(x*, t, y) dt/p
p(x*,
t, y) dt.
(3)
In this section, we give two examples of the application of (1) and (2) respectively. We concentrate our attention in Section 2 on models for which (1) is strictly appropriate, but for which a diffusion approximation, as in (3) is employed. In Section 3, we consider diffusion approximations when there are two return boundaries. Examples of the diffusion approximations are given in Section 4. Moran (1958, 1959), introduced a population genetic model in which there are two possible alleles, and in which at each unit time one gene dies and is replaced. He allowed for the possibility of mutation, and for viability or fertility differences between the two alleles. The number of genes belonging to a particular allele evolved as a random walk process, changing by at most one gene per unit time. When the process has settled down to statistical stationarity, it may be proved to be reversible. Thus although the result (1) could be used for this process, it is equally possible to discuss the age distribution by equating it to the extinction time distribution. The time since an allele was last reintroduced to the population by mutation has the same distribution, conditional on the present number of genes of that allele, as the time from the present until the next occasion (in the future), when it will again be lost. We do not go into the details
AGE
OF
AN
ALLELE
181
of this example, because they are similar in some respects to the example presented in Watterson (1976), and in order to obtain numerical results, one would probably use a diffusion approximation similar to those discussed below in Sections 2, 3, and 4. There is perhaps need for a rigorous proof of (2), especially in the case of transient or null-recurrent return processes. Dr. A. Pakes (personal communication) has established some results in this direction. Levikson (1977) illustrated the use of (2) on a birth and death process, in which the infinitesimal birth and death rates are hn and pn, respectively, whenever the population is of size n. When 0 < h < p, Levikson evaluated the right side of (2) as dt . [I - CX(t)][l - B(t)] /3(t>j-‘j$/R’,
(4)
where a(t) = CL(exp@) - 1) A expbt) - P ’
p(t) = YexpW - 1) A expW - CL ’
Y==h-#LL.
Our purpose is simply to point out that the age density (4) may be obtained as being the density of the extinction time, the time to first passageto state 0 fromj as initial state. We have for j = 1,2, 3,...,
Pdt)lLm Pdt> dt = (d/dt)[P~,dt)l, t 3 0,
(5)
which may be verified from (4) above and using
The reason for the result (5) is that the return process (ignoring any instantaneous transitions from 1 to 0 and return) is a birth and death process on the positive integers, with stationary distribution being the logarithmic distribution mi
=
-W)i/(i
141
-
i = 1) 2, 3)...)
(VP.)),
and the process is reversible because of the connection between the stationary distribution and the infinitesimal parameters: m,Ai = mi+lp(i + 1)
for
i = 1, 2, 3,...;
see Kendall (1975, (12)) for the general theory. Even in the critical case when /I = p, whether or not (2) is itself valid, the
G. A. WATTERSOX
182
result (5) may be easily verified, using the birth and death transition p&t) forj
= (At)j-i/(1
+ xqi+i
and
probabilities
pj,, =: [At/(1 + ht)]j,
= 1,2,3 ,... .
2. DIFFUSION APPROXIMATIONS; ONE RETURN BOUNDARY Some return Markov chains, for which (1) is applicable exactly, are not strictly reversible but yet have the property that their stationary age distribution (1) is approximately the extinction time distribution pjyd - p$;‘) for the nonreturn process. We establish this result by considering only those processes which have good diffusion approximations. Such processes are particularly common in population genetics. Suppose therefore that, after suitable resealing of the state space and the time variables, the original absorbing process is approximated by a diffusion with drift and diffusion parameters denoted by p(x) and us(x), respectively, and with 0 as an accessible, absorbing state (see Goel and Richter-Dyn, 1974, p. 51, for terminology). Corresponding to the discrete process being returned to state 1 whenever state 0 is reached we suppose that the approximating “return diffusion” is returned to state x*, interior to the diffusion state space. Whether or not (3) may be rigorously proved for the age of a return diffusion process, it is certainly going to yield an approximation to the exact age distribution (1) in the cases under consideration. It6 and McKean (1965, p. 149, (2)) exhibit the diffusion transition densityp(x*, t, y) as the product of two functions, one of which is a function of y only and the other is symmetrical in x* and y. (I thank Frank Norman for the reference to the proof of this fact. He has shown (personal communication) that it implies reversibility of diffusion processes under certain sufficient conditions.) Defining the “natural scale”
G(Y) = exp ) -2 ly [libel and the “speed measure”
dr 1
(6)
density
4~) = [u2(y) G(~)l-l, we have p(x*, t> Y> = m(y) 4(x*, t> Y), where 4(x*, t, y) is a symmetrical
function
4(x*, t, Y)
in x* and y:
= 4(Y, 4 x*1.
(7)
AGE
The right-hand side of (3) (omitting symmetrical in x* and y also:
OF AN
183
ALLELE
the differential
dt) may now be seen to be
Y> dt=4~)4(x*> t,Y)/ jam 4~)dx*, t,Y)dt =m(x) dy,t,x*,/lx dt 04.4dv,t,"x*1 =P(Y,t, x*)/l= P(Y>t, x*>dt.
I+*, t,y)/jomP(x*>t>
Moreover, it is usually the case in population genetics that the return state xx is arbitrarily close to 0, being of the form x* = 1/(2N) as N + co, so that the right side of (3) is approximately
P(Y,t, O-t)/jomp(y, t, O+)dt.
(8)
This approximating age density is also the density of the extinction time, t, for the absorbing diSfusion process to jirst reach the absorbing state 0 from y as initial state (cf. Crow and Kimura, 1970, 8.3.21). It is this observation which is the main point of the paper. Most authors have studied the mean age, rather than the age distribution. The mean of the approximating age distribution (3) is, in Levikson’s notation,
= B(x*, Y)/G(x*Y)
say,
where B(x*, y) = lrn tP(x*, t,y) dt and G(x*, y) = jam P(x*, t, Maruyama
Y) dt.
(1974) showed that
Y)+,4x)(W) B(x, Y)=- G(x, Y),
W4(~zl~x2) B(x> and that
W4(W~z)
‘7x,
Y)+&)(VW G(x, y)=--s(x - y).
184 Levikson
G. A. WATTERSON
(1977) noted that
A(x*,y)
= joffi G(x*, 5) G(t, Y) G’G(x*,
Y).
We write
for the (approximating) mean age when the diffusion is returned to near the 0 state, conditional on a diffusion state y now. In view of the limiting relationship between (3) and (8), we find that M(Y)
= jm MY,
t, Of)
0
qj-=
0
p(y, t, O+) dt,
which is the mean extinction time for the diffusion conditional on state y now. This is given explicitly in terms of the approximating diffusion parameters CL(.) and CT’(.), by the standard formula for extinction times:
M(y) = 2 j’ G(x) j” [u”(z) G(z)]-l d.z dx, x 0
(9)
where G(z) is as in (6); see, for instance, Goel and Richter-Dyn (1974, Table 3.6). We have written co in (9) for the upper bound of the state space of the approximating diffusion, although this may well be replaced by a finite value in some applications. The mean extinction time and the mean age may be infinite in particular cases. An example is given in Section 4.1 below. Sometimes, population geneticists are interested in processes approximated by a diffusion for which there are two absorbing, accessible, boundaries 0 and 1, but conditional on boundary 1 not being reached. If the unconditional (absorbing) diffusion has parameters p(.), up(.), then the diffusion conditioned on ultimate absorption at 0 (rather than I) has parameters /Lo(x) = p(x) i- U2(X) @o’($@&), (10)
uo‘yx) = c?(x), where Go(x) is the probability and is given by
of ultimate
absorption
in 0 from the initial state X,
(11)
@o(x) x j1 G(Y) 41 jol G(Y) dy,
I
where G(.) is as in (6). (See Ewens, 1973; Narain, 1974.) The return process corresponding to the conditioned
process
(10) has 1
185
AGE OF AN ALLELE
as an inaccessible boundary. The mean age, and the mean first extinction time, is then as given by (5) but calculated with the conditional parameters (10): M,(y)
G(x)
= 2 j’
”
G,(x) j1 [uo2(z) G,(z)]-l dz dx, z
= exp 1-2 .\” iki~)I~~~(r)l
dyi
= G(M@,(412Integrating (12) by parts and also using (13) sh ows that the mean conditional age can be written
M,,(y)
= 2 1j’ @&)[l
0
- @&)][G(s)
U”(X)]-’ dx (14)
- j’
Y
[l - @oW%(r)l
@,MW
a2CW1 dx/ * j1 G(x) dx, 0
which can be shown to be equivalent to Li (1975, Eq. (lo)), which was derived from the work of Maruyama (1974). Ag ain, our main point is that (14) is both a mean age and a mean extinction time. In fact, (14) is equivalent to the mean extinction time formula studied by Kimura and Ohta (1969b, (1)). Admittedly, their emphasis was on takingy = 1/2Nas initial state (a new mutant’s frequency) whereas we are interested in larger values of y. The equivalence in the neutral allele case was noted by Kimura and Ohta (1973, p. 208). Examples of the use of (14) are given in Section 4.2 for genetic diffusions subject to mutation and selection.
3. PROCESSESWITH Two
RETURN BOUNDARIES
Levikson (1977) h as considered Markov chains having 0 and N as absorbing states. He modifies the process by the assumption that whenever 0 is reached, an instantaneous return to 1 occurs, and whenever N is reached, an instantaneous return to N - 1 occurs. We now consider the diffusion approximations to Levikson’s results. Consider a diffusion approximation scaled so that 0 and 1 are absorbing boundaries and a return process so that state 0 is instantaneously followed by state x* (x* + 0) and state 1 is instantaneously followed by state x** (x ** G 1). The formula corresponding to Levikson’s discrete age distribution is the stationary age density for the time, A(t), that has elapsed since the process was last returned from a boundary. Kimura and Ohta (1973, p. 208) called this
186
G.
A.
WATTERSON
time the “age of the polymorphism.” density is
Using the notation of Section 2, this
limPr(A(s) E (t, t + dt) X(s) = y)!dt S’XI
(1 - @&*))p(x**,
t, Y) + @,(x**)P(x*~
6 Y)
t > 0.
(1 - @0(X*)) s,- p(x**, t, y) dt + @I@**) s,- p(x*, t, Y) dt ’
Levikson’s formula for the probability that the last return, before time t, from a boundary was from boundary 1 becomes fim Pr(last return from 1 : T(t) = y)
(1 - @,,(x*))&(x**,
t, y) dt
(1 - C&(X*)) JR p(x**, t, y) dt + @a(~**) J,- p(x*, t, Y) dt To see the role played by reversibility in these expressions, we may deduce from (6) and (7) that
I% 4 Y)
= P(Y, t, 4 dYM4,
= P(Y, t, 4 44
G(~/[uYY)
G(Y)!
We find
v-2 Pr(A(s) E (t, t + dt) ~T(s) = y)/dt = ((1 - Q&X*)) 0*(x**) G(x**)p(y, + @,,(x**) u”(x*) G(x*)P(Y,
t, .**)
t, x*>>
x /(I - @a@*)) u*(x**) G(x**) lmp(y,
t, .**)
dt
+ Q&v**) u”(x*) G(x*) sump(y, t, x*) dtl-’ and iir Pr(last return from 1 ~X(t) = y) = /(I - CD&*)) u2(x**) G(x**) JT: p(y, t, .**) dt/ x /(I - Q&e*)) ,2(x**) G(x**) j-)(y, + Qo(x**) u”(x*) G(x*) j-)$y,
t, x**) dt
t, x*) dt/-‘.
187
AGE OF AN ALLELE
But now, if the return states are allowed to approach their respective boundaries, x* --f 0 and x** -+ I, we find using (6) and (11) that [I - tDo(x*)] u2(x**) G(x**) u’(x*) G(.r*) Bo(x*)
-
5 a0 ’
where a, =- -F Thus finally,
a, ~ - da2(4 dx i . x=0
and
1 , s=1
taking limits as x* + 0, x** ---f 1, and s -+ GO, hmPr(A(s)
E (t, t + dt) ( x(s)
= y),‘dt
a,P(y, t, 1-j + a,P(y, 6 O+> al J;P(r, t, 1-I dt + a0S,-P(Y, t, Of) dt
(15)
and Lm Pr(last return from 1 1X(t) = y)
al.frp(y, t,I-)dt aI J;P(y,
4 1-I
dt + a,
J"p(y, 0
t, 0+) dt
(16)
= 1 - Do(y). The result (15) has the interpretation of being the probability density of the time (in the future) at which one or other boundary, 0 or 1, is first reached from initial state y (see, e.g., Crow and Kimura, 1970, pp. 378,379). Thus the past age density equals the future absorption time density, when returns are made to states infinitesimally close to the boundaries and when stationarity is reached. Similarly, the probability of a last return from 1, (16), has the alternative interpretation of being the probability that absorption will first occur (in the future) into boundary 1 rather than boundary 0, given initial state y (see (11) above). These results explain why certain formulas in the literature of population genetics, which were derived as properties of ages of alleles, can also be found in papers concerning the extinction of alleles. We mention some examples in Section 4.3 below. 4. EXAMPLES 4.1. Diffusion Approximation
for a Branching Process
Levikson (1977) considered a branching tion with infinitesimal parameters p(x) = (m -
1)X,
process and its diffusion c?(x) = bx.
approxima-
188
G. A. WATTERSOlX
Here, m is the mean number of offspring per individual and b is the variance. Cox and Miller (1965, Sect. 5.11) show that p(x*, t, y) is a mixture of gamma densities with one term missing: $@*,
t, y)
=
f
[e~"*"iB(.~*~i~)"!n!]
py'~-le-ylfl,(n
-
I)!,
71=1
where (y =
e’m-lM,
/3 =- b(e(wL-l)t -
1)/2(m -
If m = 1, we have 01:= 1 and /3 =- bt/2. The probability time t, from x* as initial state, is the missing term Pr(X(t)
1). of absorption
= 0 / X(0) = .Y*) = e-r**16.
We assume from now on that m C< 1, so that ultimate the nonreturn process. But the transition density may also be written as
I+*, f, Y)
=
m(r)
4(.~*,
at 0 by
(17)
extinction
is certain for
t, Y),
where m(y)
I_
2$(n"-l)"/bj(by)
and q(x*, t,y)
:== $b i
e-‘“‘+““.‘6(~*y)n
anp-?n;n! (n -
I)!,
n=1
which is symmetrical written as
in x* and y. The approximating
P(Y, t, .x*1
age density (3) may be
dY1 t, x*1
J” P(Y, t, .x*) dt = s,- q(y, t, x*) dt 0
Letting
x* -
0 we find for the limiting
expression
as in (8)
But this density for the age, t, is exactly the extinction directly from (17): (djdt) Pr(X(t)
time density calculated
= 0 1X(0) = y) = ibe -~~@yajiF,
t > 0.
AGE
OF
AN
189
ALLELE
This means that (~-l = e(l-m)t has an extreme value (Type 1) distribution, (see Johnson and Katz, 1970, Chap. 21). The mean of the absorption time distribution, and the mean age, for a given state y now, is M(y) = IL (1 - e-vu/s) dt = (;I(1
- an)) I’
(1 - e-2(1-m)zlb) x-‘( 1 + x/y)-’
dx
(18)
using the substitution x = y~(l-“)~/(l - e(l-l”)t), assuming m < 1. This can also be derived from (9). The special case when m = 1 and b = 2 most interested Levikson. We find that the age density, given state y and as x* -+ 0, is e-gltyt-2
t
2
(19)
0.
This distribution is an extreme value (Type 2 or Cauchy) distribution (see Johnson and Kotz, 1970, Chap. 21). It does not have a finite mean, as can be seen directly, or by letting m t 1 in (18). But the modal age is t = y/2, which increases as y increases. These remarks agree with observations made by Levikson (1977) for a corresponding discrete branching process. For our diffusion, the reciprocal of the age has an exponential distribution, with mean 1/y, as x* + 0. As a check on the accuracy of the diffusion approximation (19) to the actual age distribution as tabulated by Levikson (1977) for the branching process with geometric offspring distribution pn = ($)n+l, n = 0, 1, 2, 3,..., we make the comparison in the form Pr(age E n / 2 = j) + In+' e-i/tjt-2 dt n-t = e-j/'n+:' . - e-j/(s-:,),
112 ,-jjtj-2 : 10
z
dt
=
e-2i
for
n =T 1, 2,...,
for
n = 0.
In Table I, we show the left side, as tabulated by Levikson (1977), with the approximating right side in parentheses. Clearly, the diffusion approximation is accurate for large values of j, and is still moderately accurate for small j but large values of n. 4.2. Genetic Da&ion,
with Mutation
and Selection
The most common population genetic diffusion models having 0 as an absorbing state have parameters of the form
190
G. A. WATTERSON
TABLE .4ge Probability
Distributions
I
for Branching
Process (Diffusion
Approximation)
Present frequency
Age
1
2 -~___
3
4
5
__-
6 ~~~ --
8
9
.608 (.135)
(.O! 8)
C.0002)
$0)
(.c&
co&
(.~O)
(.OiO)
(.OiO)
.I51 (.378)
.282 (.245)
.193 (.133)
.126 (.069)
.078 (.036)
.047 (.018)
.027 (.009)
.016 (.005)
(.ow
2
.067 (.157)
.167 (.186)
.152 (.166)
,132 (.132)
.lOl (.lOO)
,088 (.072)
.068 (.051)
.052 (.036)
.039 (.025)
3
,038 (.081)
.105 (.115)
,109 (.123)
.106 (.117)
.099 (.104)
.089 (.089)
,078 (.075)
.067 (.061)
.056 (.049)
4
.024 (.049)
.072 (.076)
,079 (.089)
.082 (.092)
.082 (.090)
.079 (.084)
.073 (.076)
,067 (.067)
,060 (.059)
5
.017 (.033)
.0.52 (.054)
.060 (.066)
.065 (.072)
.067 (.074)
.067 (.072)
.065 (.069)
.062 (.064)
.058 C.059)
6
.012 (.024)
.040 (.040)
.046 (.051)
.052 (.057)
.056 (.060)
.057 (.061)
.057 (.061)
.055 (.059)
.054 (.056)
7
.009 (.018)
.030 (.031)
.037 (.040)
.042 (.046)
.046 (.050)
.048 (.052)
.049 (.053)
.049 (.052)
,048 (.05 1)
8
.008 (.014)
.025 (.024)
.030 (.032)
.035 (.038)
.039 (0.42)
.041 (.044)
.043
.043
C.046) (.046)
.043 (.046)
9
.006 (.Oll)
,020 (.020)
.025 (.027)
,029 (.032)
.033 (.035)
.035 (.038)
.037 (.040)
,038 (.041)
.039 (.041)
10
305 (.009)
.017 (.016)
.021 (.022)
.025 (.027)
.028 (.030)
.031 (.033)
.033 (.035)
.034 (.036)
035 (.037)
15
.002 (.004)
.008
(.C’W
.OlO (.Oll)
.013 (.014)
.015 (.016)
.017 (.018)
.019 (.020)
,020 (.021)
.021 (.022)
.0014 (.002)
,005 (.@w
.006
(.oW
.008 (.008)
.009 (.OlO)
.Oll (.Oll)
.012 (.012)
,013 (.013)
.014 (.014)
.0009
.004 C.004)
.005 (035)
.006 (.007)
.007 (.008)
.008
.009
(.ocw
.003 (.003)
C.0’38) (.009)
.OlO (.OlO)
.0006 (.OOl)
,002 (.002)
003 C.003)
.004 (.004)
.005 (.005)
.005 C.005)
.006 (.006)
.007 (.007)
.007 (.007)
0
1
20
25
30
.009
AGE OF AN
191
ALLELE
p.(x) = 2s{h + (1 - 2h) X} X( I - X) - UX, (20)
U”(X) = X(1 - X)/2Ne ,
where 0 < x < 1, N, is the effective (diploid) population size, u is a mutation rate (0 < u < I), and s and h are selection and dominance parameters. Maruyama (1974) writes F = 1 - 4N,u, and S = 4N,s, so that p(x) = (2S[h + (1 - 2h)x](l - x) - (1 - F)}x/4Ne
.
We find from (6) that G(x) = (1 - ~)~-l exp{-S[2hx
+ (1 - 2h)x’]}.
(21)
We consider both the case when the age of an allele is calculated allowing for quasi-fixation occasions, and the case when the allele frequency is conditioned so as not to achieve the boundary value 1 during the lifetime of the allele. When u > 0 and hence F < 1, the back mutation is such as to make the boundary 1 nonabsorbing. The mean age of an allele of frequency y is, by (9), (2% and @I), M(y) = 4N, 1’ (1 - .~)~--lexp{ -S[2hx 0
+ (1 - 2h) a?])
1
X
sz
z-‘( I - z)-” exp(S[2hz + (1 - 2h) z2] dz dx.
(22)
The selectively neutral case provides a simpler expression; when S = 0 we have, integrating by parts, M(y)
= (4Ne.lF) 1-In y + J1 z-‘[( 1 - x)-’
- l] dz
0
- (1 - y)” j: z-l(l
- z+” dzj.
(23)
Priority for this result as an e&z&on time formula belongs to Ewens (1964, (2.6)) and as an average age formula, to Maruyama (1974, (12)). Two special cases deserve closer attention; when F = - 1, (23) reduces to WY)
= -4NAy
lnr)/(l
- Y),
(24)
while if F = 0, (24) is replaced by M(y) = -4N,
1’ (In z)/(l - z) dz, 0
(25)
192
C. A. WATTERSON
which approaches 4N$/6 as y T 1. Numerical values of M(y), as in (23) (24), (25) have been published by Maruyama (1974, Table 2) but these are sometimes incorrect in the third significant figure, no doubt due to numerical quadratures. If we wish to discuss the mean age (or mean extinction time) conditional on no intervening quasi-fixation occurrences, then if F < 0, the previous formulas (22)-(25) apply. We have Jfo(.?9 = M(Y)
if
F < 0,
because the boundary N = 1 is inaccessible in these cases. Hence column 4Nu = 1 in Maruyama’s Table 1 should agree with that column in his Table 2 and with (25) above. If, however, F > 0 so that x = 1 is an accessible boundary, the diffusion conditional on boundary 1 not being reached differs from the unconditional one. Li (1975) has tabulated the mean age M,,(y) as given by (14), for various s and h combinations, but in the absence of back mutation (u = 0, F = 1). Note that Li’s s corresponds to our -2s. There is some simplification of (14) possible in the case of no mutation and no dominance. With F = I, h -= Q , we have the conditional diffusion parameters corresponding to (20) as given by (10): pa(x) = (1/4N,) x( 1 -
.x) S(e+ + epss)/(e@ -
ecS”)
(26)
and uo2(x) =z (1/2N,) x(1 - x), which lead to the conditional
mean age, from (14) of
~ 1 (1 - e-s(l-s))(e-sY - e-S%)dr 4N, I. s7J e-S) x(17 Iy S(e-
(27)
This is identical with (13) in Marugama (1974), who derived it using a different argument. Kimura and Ohta (1969a, (17)) g ave essentially this result as an extinction mean time. If S = 0, we find the limiting (neutral allele) conditional age as MO(y) == --4N,y
lny/(l
- Y).
(28)
This is identical with (24); the reason will be discussed later. Maruyama (1974, Table 3) has tabulated values of (27) f or various S values. His column / 4Ns 1 = 0 should agree with (28), and hence with his column 4Nu = 0 in Table 1 and column 4Nu = 2 in Table 2; discrepancies must be due to numerical integration procedures.
AGE OF AN
193
ALLELE
As Maruyama noted, M,,(y) in (27) is an even function of S; that is, for a given allele frequency y, the conditional mean age is independent of whether selection is in favor of, or against, the allele. This appears at first sight to be paradoxical. The mathematical explanation lies in the fact that the conditional diffusion (26) is completely independent of the sign of S. The statisticaE explanation of the paradox is the following. Let X(t), 0 < t < T, be the realization of the unconditional diffusion with parameters p(x) =J Sx(I - x)/4Ne )
u’(x) = x(1 - x)/2Ne .
The stopping time, T, is that time when one or other boundary x = 0 or x = 1 is first reached. The likelihood of that realization, with respect to the measure of the diffusion having S = 0, is
= exp [$S[X( T) - X(O)] - & N:lS2 L’ X(t)( 1 - X(t)) dt/. This likelihood follows from the work of Brown and Hewitt (1975). The stochastic integral si X(t)( 1 - X(t)) dt is associated with the parameter F, and clearly contains information about the magnitude, but not the sign of S. The maximum likelihood estimator of S is, incidentally, L? = 4N,[X(T)
- X(O)]//rX(t)(l
- X(t)) dt.
0
By discussing the diffusion (26), which is conditional on the boundary 1 not being reached, we have conditioned on the sufficient statistic for the sign of S, X(T) - X(O), taking only the negative value -X(O), thus removing the dependence of the realization on the sign of S. In the case when there is mutation but there is no selection, then provided F > 0 so that boundary 1 is accessible, we have the conditional diffusion parameters
PO(X)= -4
+ F)/4N, ,
q,‘(x) = x(1 - x)/2Ne ,
which are similar to the unconditional diffusion parameters. Then the average age and the average extinction time may be derived from (14) to be
dx+
(1
Jy)'
i y
’ (’ ; X)F dsj, (29)
a result equivalent to Maruyama (1974, (10)) and Kimura and Ohta (1973,
194
G. A. WATTERSON
(13a)). Of course, (29) is just (23) with F replaced by -F, and so, in a sense, can be attributed to Ewens (1964) as an extinction time formula. Letting F + 1 yields the conditional mean age, as given by (28) in the absence of mutation and selection, and as given by (24) for the unconditional mean age with F = - 1. We have concentrated our attention on the mean age, either conditional or unconditional on no fixation. In many of the above examples, explicit formulas are known for the transition probability densities entering into the age density (8). Indeed, the time to extinction density (and hence the age density) for an allele of frequency y now, in the absence of mutation and selection and conditional on no fixation is ytgf’(2i+
l)(-l)i-1-lF(l
-i,i+2,2;
1 -y)hie-W,
t>O
where hi = i(i + 1)/4N, (see Kimura, 1970). It is this density which has mean i%!,(y) as in (28). M oreover, even though some mutation and selection cases have complicated transition densities, the approximation obtained by Voronka and Keller (1975, Sect. 7) may well be suitable, at least for relatively low ages. The higher order moments in some of the cases above can also be obtained explicitly. 4.3. Genetic Dzjfusion with Two Return Boundaries
We return to the diffusion whose parameters are given by (20), but specialized by putting u = 0: p(x) = 2s{h + (1 - 2h)x) x(1 - x),
and G(x) = x(1 - %)/2N, . Both boundaries x = 0 and x = 1 are now absorbing and accessible. In Section 4.2 we discussed, inter alia, the case when the process was conditional on 1 not being reached and the boundary 0 was a return boundary. We now consider the case, as discussed generally in Section 3, when both 0 and 1 are return boundaries. In the absence of selection, the mean age of a polymorphism was given by Kimura and Ohta (1973, p. 208) as -4N,[y
lny + (1 - A 141 - 391
when the return process is in state y. This is the mean of the distribution (15), and was first found by Watterson (1962) as the mean absorption time. The
AGE
OF
AN
195
ALLELE
probability that the last return was from state 1 is a special case of (16); Kimura and Ohta (1973) state this probability is y in the neutral alleles case. Of course, this is the well known result for the fixation probability of a neutral allele of frequency y. The full age density (15) may be obtained from Kimura’s (1955) work on fixation time densities: lim Pr(A(s) E (t, +dt) =y(l
-y)
f
1X(s) = y)/dt h,(2i+
l)(-l)i”[F(l
--,i
+2,2,y)
i=l
+ F(l
- i, i + 2, 2, 1 - y)] e&it,
where Ai = i(i + 1)/4Ne . The age density (and the absorption density) has a complicated form in the selection case, even when h = $ corresponding to no dominance effects, or when h = 0 or 1 corresponding to complete dominance of one or other allele; see Crow and Kimura (1970, 396-414). The approximation given by Voronka and Keller (1975) may prove useful here. The mean age may be equated to the mean absorption time from state y (see the Correction to W’atterson (1962)). The probability that the age corresponds to a return from boundary 1, (16), reduces to the fixation probability given by Kimura (1957, p. 896):
1-
@o(y) = 1 G(x) dx/(
where, for the process here, G(X) may be obtained G(x) = exp(-S[2hx
+ (1 -
G(x) dx, from (21) as 2h)x2]},
with S = 4N,s. ACKNOWLEDGMENTS I thank B. Brown paper.
and F. Norman
for helpful
suggestions
in the preparation
of this
REFERENCES B. M., AND HEWITT, J. I. 1975. Asymptotic likelihood theory cesses, J. Appl. Prob. 12, 228-238. Cox, D. R., AND MILLER, H. D. 1965. The Theory of Stochastic New York. CROW, J., AND KIMURA, M. 1970. An Introduction to Population Harper & Row, New York. EWENS, W. J. 1964. The pseudo-transient distribution and its uses in Probability 1, 141-156. BROWN,
^__653172124
for diffusion
pro-
Processes, Wiley, Genetics
Theory,
genetics, J. Appl.
196
G. A. WATTERSON
EWENS, W. J. 1973. Conditional diffusion processes in population genetics, Theor. Popul. Biol. 4, 21-30. FELLER, W. 1968. An Introduction to Probability Theory and Its Applications, Vol. I, 3rd Ed., Wiley, New York. GOEL, N. S., AND RICHTER-DYN, N. 1974. Stochastic Models in Biology, Academic Press, New York. 1~6, K., AND MCKEAN, H. 1965. Diffusion processes and their sample paths, Academic Press, New York and Springer-Verlag, Berlin. JOHNSON, N. L., AND KOTZ, S. 1970. Continuous Univariate Distributions I. Houghton Mifflin, Boston. KENDALL, D. G. 1975. Some problems in mathematical genealogy, in “Perspectives in Probability and Statistics: Papers in Honour of M. S. Bartlett” (J. Gani., Ed.), pp. 325-345, Applied Probability Trust and Academic Press, London. KIMKJRA, M. 1955. Solution of a process of random genetic drift with a continuous model, Proc. Nat. Acad. Sci. U.S.A. 41, 144-I 50. KIILIURA, M. 1957. Some problems of stochastic processes in genetics, Ann. Math. Statist. 28, 882-901. KIMURA, M. 1970. The length of time required for a selectively neutral mutant to reach fixation through random frequency drift in a finite population, &net. Res. Camb. 15 131-133. KIMURA, M., AND OHTA, T. 1969a. The average number of generations until fixation of a mutant gene in a finite population, Genetics 61, 763-771. KIMURA, M., OHTA, T. 1969b. The average number of generations until extinction of an individual mutant gene in a finite population, Genetics 63, 701-709. KIMURA, M., AND OHTA, T. 1973. The age of a neutral mutant persisting in a finite population, Genetics 75, 199-212. LEVIKSON, B. 1977. The age distribution of Markov processes. J. ilppl. Probability, to appear. LI, W. H. 1975. The first arrival time and mean age of a deleterious mutant gene in a finite population, Amer. J. Hum. Genet. 21, 274-286. MARUYAMA, T. 1974. The age of an allele in a finite population, Genet. Res. Camb. 23, 137-143. MORAN, P. A. P. 1958. Random processes in genetics, Proc. Camb. Phil. Sot. 54, 60-11. MORAN, P. A. P. 1959. The survival of a mutant gene under selection, J. Austral. Math. Sot. 1, 121-126. NARAIN, P. 1974. The conditional diffusion equation and its use in population genetics, J. Roy. Stat. Sot. B 36, 258-266. PAKES, A. G. 1977. On the age distribution of a Markov chain, J. Appl. Probability, to appear. SAWYER,S. 1977. On the past history of an allele now known to have frequencyp, to appear. VORONKA, R., AND KELLER, J. B. 1975. Asymptotic analysis of stochastic models in population genetics, Math. Biosci. 25, 331-362. WATTERSON, G. A. 1962. Some theoretical aspects of diffusion theory in population genetics, Anti. M&h. Statist. 33, 939-957. (1963) Correction 34, 352. WATTERSON, G. A. 1976. Reversibility and the age of an allele. I. Moran’s infinitely many neutral alleles model, Theor. Pop. Biol. 10, 239-253.