# Modeling Neutral Evolution Using an Infinite-Allele Markov Branching Process

Authors: Xiaowei Wu; Marek Kimmel
Journal: International Journal of Stochastic Analysis
 (2013)

---

## Abstract

We consider an infinite-allele Markov branching process (IAMBP). Our main focus is the frequency spectrum of this process, that is, the proportion of alleles having a given number of copies at a specified time point. We derive the variance of the frequency spectrum, which is useful for interval estimation and hypothesis testing for process parameters. In addition, for a class of special IAMBP with birth and death offspring distribution, we show that the mean of its limiting frequency spectrum has an explicit form in terms of the hypergeometric function. We also derive an asymptotic expression for convergence rate to the limit. Simulations are used to illustrate the results for the birth and death process.


---

## Body

## 1. Introduction

The infinite-allele branching process was first introduced by Griffiths and Pakes [1]. As a special type of branching process, this process allows individuals to mutate into infinitely many allelic variants, each of which is “new” in the sense of being different from all previously existing variants. This idealization is approximately correct for rare point mutations in long DNA sequences. Fundamental results for the discrete-time case (simple branching process) and for the continuous-time case (Markov branching process) have been obtained in [1, 2]. These include the number of alleles at a given generation or time, the generation number or time of the last mutation, and the limiting frequency spectrum. There exists an analogy between the results for the discrete-time and the continuous-time cases; however, the characteristics in the continuous case are relatively easier to derive [2]. Many evolutionary processes may be considered time continuous, and frequently we assume Markov property in modeling. A classical example is the discrete-time Wright-Fisher model, which is typically either approximated by a continuous-time diffusion or replaced by a continuous-time Markov chain, the so-called continuous-time Moran process [3]. Therefore, the time-continuous infinite-allele Markov branching process (TCIAMBP, or simply IAMBP) seems to be appropriate for modeling evolution in population genetics.Consider a Markov branching process with neutral mutations. Suppose that the process starts from a group of individuals carrying the same allele, and individuals can mutate into new allelic variants. We assume that the mutation is independent of the previous history of the process, and the offspring distribution is independent of the allelic type, that is, the selection is neutral for all alleles. The process can be described as an “infinitely-many-alleles” model (IAM). Whenever a mutation happens, it yields a new allele, which differs from all the previously existing ones. In this paper, we are interested in the frequency spectrum of the IAMBP, which may be defined as the number or proportion of alleles present in a given number of individuals at a specified time point. Frequency spectrum in this paper refers to random allele frequencies, not their expected value as in Griffiths and Pakes [1] and in Pakes [2], since we also consider the variance of the allele frequency later on. Unless specified otherwise, we will use terms “mean frequency spectrum” and “variance frequency spectrum” in the remainder of this paper to denote the expected value and variance of the allele frequency. The frequency spectrum plays an important role in many genetic processes, such as DNA sequence evolution. As an example, Kimmel and Mathaes [4] modeled the Alu sequence data using an infinite-allele simple branching process with linear-fractional offspring distribution, and the goodness of fit testing suggested that Alu sequences do not evolve neutrally and might be under selection. It has to be noted that the concept of the frequency spectrum is in some sense similar to the Ewens’ sampling formula [5] in population genetics. We will return to this subject in the discussion, although analysis of the analogies and differences transcends the scope of the present paper.The paper is organized as follows. In Section2, we rigorously define the IAMBP and the mean frequency spectrum of the IAMBP. Then, we provide explicit expressions for the special case of the birth and death process. In Section 3, we derive the variance frequency spectrum and discuss its use in interval estimation for process parameters. We perform simulations to illustrate the results using the birth and death process example in Section 4. Section 5 is a summary.


## 2. IAMBP and Its Limiting Mean Frequency Spectrum

### 2.1. Definition and Basic Properties in the Supercritical Case

Let us consider a continuous-time Markov branching process consisting of individuals with exponential life spans with meana-1. Let us assume that upon death, each individual produces a random number of offspring. As usually assumed, the offspring counts are identically distributed according to probability generating function (pgf) f(s), and they are independent conditional on the past process. The mean f′(1-) of the offspring distribution is m, regardless of the allelic type. We further assume that a newborn individual mutates into a new allelic type with probability μ independently of the previous history of the process. Let us denote by h(s)=f(μ+(1-μ)s) the offspring pgf in a clone, started by the overall ancestor or any of mutants, containing only the like-type individuals. The entire process is a union over all individual types of such clones. The theory of the IAMBP has been developed by Griffiths and Pakes [1] in the discrete-time case and then by Pakes [2] in the continuous-time case. We will assume m>1 and M=h′(1-)>1, although some results can be proved without this latter assumption.Letαt(j) be the number of alleles present in j individuals at time t and ϕi,t(j)=Ei[αt(j)], where subscript i indicates that the process begins with i individuals carrying the same allele. It has been shown that [2]
(1)ϕi,t(j)=qij(t)+iamμeλt∫0t‍e-λxq1j(x)dx,j≥0,
where λ=a(m-1) is the Malthusian parameter of the overall process and qij(t) is the probability of observing j individuals (j≥1) carrying the parental allele at time t when starting from i individuals with the parental allele at time t=0. Consequently, for the number Kt of alleles at time t, we have
(2)Ei[Kt]=∑j=1∞‍ϕi,t(j)=1-qi0(t)+iamμeλt∫0t‍e-λx[1-q10(x)]dx.LetGj=∫0∞‍e-λtq1j(t)dt, j≥0. If we define ψij(t)=ϕi,t(j)/Ei[Kt] and
(3)ψj=limt→∞ψij(t)=limt→∞qij(t)+iamμeλt∫0t‍e-λxq1j(x)dx1-qi0(t)+iamμeλt∫0t‍e-λx[1-q10(x)]dx
as the limiting mean frequency spectrum, that is, the expected proportion of alleles present in j individuals as t→∞, then we see that for the supercritical process such that λ>0,
(4)ψj=λGj1-λG0,j≥1.IfM>1, then the process of the like-type clones is supercritical, and as it is known [6], q10(t)↑q10(∞)<1 and q1j(t)→0, j≥1, as t→∞. Therefore, eλt|∫0t‍e-λxq10(x)dx-q10(∞)/λ|→0 and eλt∫0t‍e-λxq1j(x)dx→0 as t→∞, for j≥1. This yields the following asymptotic equivalence:
(5)ψij(t)-ψj~t→∞λGj[1-q10(∞)-(λ(1-qi0(∞))/iamμ)](1-λG0)2e-λt.
Details of the proof are omitted, since they appear elementary.


### 2.2. IAMBP with Birth and Death Offspring Distribution

For the IAMBP with birth and death offspring distributionf(s)=α+βs2, α+β=1, we are able to obtain an explicit form for Gj, j≥0; therefore, the limiting mean frequency spectrum ψj,j≥1 can be derived. The offspring pgf of the like-type individuals clone in the birth and death IAMBP is written as
(6)h(s)=f(μ+(1-μ)s)=α+β[μ+(1-μ)s]2,
where α,β and μ stand for the death, birth, and mutation probabilities for every individual and α+β=1. Note that under another parameterization where the two newborn individuals die, live, and mutate independently, this pgf may be formulated differently as h(s)=[α+βμ+β(1-μ)s]2. Under either parameterization, λ=a(2β-1). If, as assumed, M=m(1-μ)>1, then parameters α and μ are subject to a constraint
(7)(1-α)(1-μ)>12.Let us writeA2=α+βμ2 and B2=β(1-μ)2 (note, for the other formulation, A2=(α+βμ)2 and B2=β2(1-μ)2). The explicit form of Gj can be written as
(8)G0=1cA2B2Γ(λ/c)Γ(2)Γ(2+(λ/c))×F(1,λc;2+λc;A2B2),Gj=1c(1-A2B2)2Γ(1+(λ/c))Γ(j)Γ(j+1+(λ/c))×F(j+1,1+λc;j+1+λc;A2B2),j≥1,
where c=a(B2-A2)=a[2β(1-μ)-1] is the Malthusian parameter of the like-type clone and F(·,·;·;·) is the Gauss hypergeometric function [7], defined as
(9)F(a,b;c;z)=Γ(c)Γ(b)Γ(c-b)∫01‍tb-1(1-t)c-b-1(1-tz)-adt,c>b>0.
For a detailed derivation, see Appendix A. Note that the supercritical condition also guarantees that the argument of the hypergeometric function remains within its region of definiteness.It follows that(10)ψj=λGj1-λG0=(λc(1-A2B2)2×(Γ(1+λ/c)Γ(j)Γ(j+1+λ/c))×F(j+1,1+λc;j+1+λc;A2B2)(1-A2B2)2)×(1-λcA2B2Γ(λ/c)Γ(2)Γ(2+λ/c)×F(1,λc;2+λc;A2B2))-1,j≥1.Figure1 shows an example of the limiting mean frequency spectrum for the birth and death process with parameters a=1, α=0.25, and μ=10-4, based on formula (10). To see how the spectrum varies with different parameter settings, we plot in Figure 2(a), the 3-D surface of a major component of the spectrum, ψ1, for different α’s and μ’s. Figures 2(b) and 2(c) illustrate the effect of one parameter on ψ1 given a fixed value of the other parameter.Figure 1
Limiting mean frequency spectrum of the infinite-allele birth and death process witha=1, α=0.25, and μ=10-4.(a) Surface ofψ1 at different α and μ, for the infinite-allele birth and death process. (b) Relation between ψ1 and μ when fixing α. (c) Relation between ψ1 and α when fixing μ.
(a)
(b)
(c)We see that for fixedα, increasing μ causes an increase of ψ1. This can be intuitively explained by the offspring pgf h(s) of the like-type clone. From the pgf expression h(s)=α+β[μ+(1-μ)s]2, we see that the probability of obtaining one like-type individual in the offspring is 2(1-α)μ(1-μ), which is an increasing function of μ for a given α, under the constraint (1-α)(1-μ)>1/2. Therefore, increasing μ will finally lead to an increase of ψ1. The effect of α on ψ1 when fixing μ is not so obvious, but we notice that when fixing μ very close to 0, as α approaches 1/2, the process is approximately critical binary fission; therefore, ψ1 drops down because of almost sure extinction of the process, as seen from the tail behavior of the solid thick line in Figure 2(c).Arguably, the frequency spectrum can only be observed in finite time. The finite-time mean frequency spectrum can be obtained by computingGj(t)=∫0t‍e-λxq1j(x)dx, j≥0 numerically. For the birth and death process, this involves the computation of the incomplete hypergeometric function. The following is a valid question in this context. In order to safely use the limiting mean frequency spectrum, how long should the process history be? Figure 3(a) compares the limiting mean frequency spectrum with some long-term mean frequency spectra, for the birth and death process with parameters a=1, α=0.25, and μ=10-4. We see that under this setting, the long-term mean frequency spectrum is almost identical to the limiting mean frequency spectrum when t≥28. In general, this result depends strongly on parameters a, α, and μ, for example, small μ leads to longer t. This provides us with some intuitions concerning the sufficiently large t for approximating the limiting mean frequency spectrum. Figure 3(b) illustrates the difference between the finite-time mean frequency spectrum and the limiting mean frequency spectrum as a function of t, for large t, t∈[15,35] and for j=1,2, where lines represent the true difference and markers represent the asymptotic approximation by formula (5). To emphasize the agreement for t large, this figure is plotted in semilogarithmic scale. We see that the true difference drops exponentially fast, and the asymptotic approximation is good for large t.Comparison between the finite-time mean frequency spectrum and the limiting mean frequency spectrum of the infinite-allele birth and death process. (a)ψ1j(t) and ψj for 1≤j≤15, t=15,20,28 and ∞. (b) Difference between ψ1j(t) and ψj as a function of t, t∈[15,35] for j=1,2. Lines represent the true difference and markers represent asymptotic approximations.
(a)
(b)Given the observed long-term mean frequency spectrum, the parametersθ of the IAMBP, such as α,μ in the birth and death process, can be estimated by equating the observed long-term mean frequency spectrum ψobs from the sample to the expected limiting mean frequency spectrum ψexp from formula (3) and solving for the process parameters. In the case of the birth and death process, we may estimate α and μ for example by solving
(11)(A2B2Γ(j1)Γ(k1+1+λc)×F(A2B2j1+1,1+λc;j1+1+λc;A2B2))×(A2B2Γ(k1)Γ(j1+1+λc)×F(A2B2k1+1,1+λc;k1+1+λc;A2B2))-1=ψobs(j1)ψobs(k1)(A2B2Γ(j2)Γ(k2+1+λc)×F(A2B2j2+1,1+λc;j2+1+λc;A2B2))×(A2B2Γ(k2)Γ(j2+1+λc)×F(A2B2k2+1,1+λc;k2+1+λc;A2B2))-1=ψobs(j2)ψobs(k2)
for positive integers j1≠k1, j2≠k2, where λ/c and A2/B2 are both functions of α and μ.There is no explicit solution for such estimator, but numerical search according to some criteria is feasible. Another possibility is to minimize the distance (such as thel2 norm) between the observed long-term mean frequency spectrum and the expected limiting mean frequency spectrum, that is, θ^=argminθ∥ψobs-ψexp(θ)∥2.The estimated parameters can be used to check the goodness of fit of the IAMBP model. Another interesting problem is to test whether two sets of parameters are identical, given two observed mean frequency spectra. A simple approach is to use Pearson’sχ2 test, such as in Kimmel and Mathaes [4]. However, there may be restrictions to applying the χ2 test, such as small cell counts and inappropriateness due to the finite length of the observed spectrum. This motivates us to develop an interval estimator for the IAMBP parameters.


## 2.1. Definition and Basic Properties in the Supercritical Case

Let us consider a continuous-time Markov branching process consisting of individuals with exponential life spans with meana-1. Let us assume that upon death, each individual produces a random number of offspring. As usually assumed, the offspring counts are identically distributed according to probability generating function (pgf) f(s), and they are independent conditional on the past process. The mean f′(1-) of the offspring distribution is m, regardless of the allelic type. We further assume that a newborn individual mutates into a new allelic type with probability μ independently of the previous history of the process. Let us denote by h(s)=f(μ+(1-μ)s) the offspring pgf in a clone, started by the overall ancestor or any of mutants, containing only the like-type individuals. The entire process is a union over all individual types of such clones. The theory of the IAMBP has been developed by Griffiths and Pakes [1] in the discrete-time case and then by Pakes [2] in the continuous-time case. We will assume m>1 and M=h′(1-)>1, although some results can be proved without this latter assumption.Letαt(j) be the number of alleles present in j individuals at time t and ϕi,t(j)=Ei[αt(j)], where subscript i indicates that the process begins with i individuals carrying the same allele. It has been shown that [2]
(1)ϕi,t(j)=qij(t)+iamμeλt∫0t‍e-λxq1j(x)dx,j≥0,
where λ=a(m-1) is the Malthusian parameter of the overall process and qij(t) is the probability of observing j individuals (j≥1) carrying the parental allele at time t when starting from i individuals with the parental allele at time t=0. Consequently, for the number Kt of alleles at time t, we have
(2)Ei[Kt]=∑j=1∞‍ϕi,t(j)=1-qi0(t)+iamμeλt∫0t‍e-λx[1-q10(x)]dx.LetGj=∫0∞‍e-λtq1j(t)dt, j≥0. If we define ψij(t)=ϕi,t(j)/Ei[Kt] and
(3)ψj=limt→∞ψij(t)=limt→∞qij(t)+iamμeλt∫0t‍e-λxq1j(x)dx1-qi0(t)+iamμeλt∫0t‍e-λx[1-q10(x)]dx
as the limiting mean frequency spectrum, that is, the expected proportion of alleles present in j individuals as t→∞, then we see that for the supercritical process such that λ>0,
(4)ψj=λGj1-λG0,j≥1.IfM>1, then the process of the like-type clones is supercritical, and as it is known [6], q10(t)↑q10(∞)<1 and q1j(t)→0, j≥1, as t→∞. Therefore, eλt|∫0t‍e-λxq10(x)dx-q10(∞)/λ|→0 and eλt∫0t‍e-λxq1j(x)dx→0 as t→∞, for j≥1. This yields the following asymptotic equivalence:
(5)ψij(t)-ψj~t→∞λGj[1-q10(∞)-(λ(1-qi0(∞))/iamμ)](1-λG0)2e-λt.
Details of the proof are omitted, since they appear elementary.


## 2.2. IAMBP with Birth and Death Offspring Distribution

For the IAMBP with birth and death offspring distributionf(s)=α+βs2, α+β=1, we are able to obtain an explicit form for Gj, j≥0; therefore, the limiting mean frequency spectrum ψj,j≥1 can be derived. The offspring pgf of the like-type individuals clone in the birth and death IAMBP is written as
(6)h(s)=f(μ+(1-μ)s)=α+β[μ+(1-μ)s]2,
where α,β and μ stand for the death, birth, and mutation probabilities for every individual and α+β=1. Note that under another parameterization where the two newborn individuals die, live, and mutate independently, this pgf may be formulated differently as h(s)=[α+βμ+β(1-μ)s]2. Under either parameterization, λ=a(2β-1). If, as assumed, M=m(1-μ)>1, then parameters α and μ are subject to a constraint
(7)(1-α)(1-μ)>12.Let us writeA2=α+βμ2 and B2=β(1-μ)2 (note, for the other formulation, A2=(α+βμ)2 and B2=β2(1-μ)2). The explicit form of Gj can be written as
(8)G0=1cA2B2Γ(λ/c)Γ(2)Γ(2+(λ/c))×F(1,λc;2+λc;A2B2),Gj=1c(1-A2B2)2Γ(1+(λ/c))Γ(j)Γ(j+1+(λ/c))×F(j+1,1+λc;j+1+λc;A2B2),j≥1,
where c=a(B2-A2)=a[2β(1-μ)-1] is the Malthusian parameter of the like-type clone and F(·,·;·;·) is the Gauss hypergeometric function [7], defined as
(9)F(a,b;c;z)=Γ(c)Γ(b)Γ(c-b)∫01‍tb-1(1-t)c-b-1(1-tz)-adt,c>b>0.
For a detailed derivation, see Appendix A. Note that the supercritical condition also guarantees that the argument of the hypergeometric function remains within its region of definiteness.It follows that(10)ψj=λGj1-λG0=(λc(1-A2B2)2×(Γ(1+λ/c)Γ(j)Γ(j+1+λ/c))×F(j+1,1+λc;j+1+λc;A2B2)(1-A2B2)2)×(1-λcA2B2Γ(λ/c)Γ(2)Γ(2+λ/c)×F(1,λc;2+λc;A2B2))-1,j≥1.Figure1 shows an example of the limiting mean frequency spectrum for the birth and death process with parameters a=1, α=0.25, and μ=10-4, based on formula (10). To see how the spectrum varies with different parameter settings, we plot in Figure 2(a), the 3-D surface of a major component of the spectrum, ψ1, for different α’s and μ’s. Figures 2(b) and 2(c) illustrate the effect of one parameter on ψ1 given a fixed value of the other parameter.Figure 1
Limiting mean frequency spectrum of the infinite-allele birth and death process witha=1, α=0.25, and μ=10-4.(a) Surface ofψ1 at different α and μ, for the infinite-allele birth and death process. (b) Relation between ψ1 and μ when fixing α. (c) Relation between ψ1 and α when fixing μ.
(a)
(b)
(c)We see that for fixedα, increasing μ causes an increase of ψ1. This can be intuitively explained by the offspring pgf h(s) of the like-type clone. From the pgf expression h(s)=α+β[μ+(1-μ)s]2, we see that the probability of obtaining one like-type individual in the offspring is 2(1-α)μ(1-μ), which is an increasing function of μ for a given α, under the constraint (1-α)(1-μ)>1/2. Therefore, increasing μ will finally lead to an increase of ψ1. The effect of α on ψ1 when fixing μ is not so obvious, but we notice that when fixing μ very close to 0, as α approaches 1/2, the process is approximately critical binary fission; therefore, ψ1 drops down because of almost sure extinction of the process, as seen from the tail behavior of the solid thick line in Figure 2(c).Arguably, the frequency spectrum can only be observed in finite time. The finite-time mean frequency spectrum can be obtained by computingGj(t)=∫0t‍e-λxq1j(x)dx, j≥0 numerically. For the birth and death process, this involves the computation of the incomplete hypergeometric function. The following is a valid question in this context. In order to safely use the limiting mean frequency spectrum, how long should the process history be? Figure 3(a) compares the limiting mean frequency spectrum with some long-term mean frequency spectra, for the birth and death process with parameters a=1, α=0.25, and μ=10-4. We see that under this setting, the long-term mean frequency spectrum is almost identical to the limiting mean frequency spectrum when t≥28. In general, this result depends strongly on parameters a, α, and μ, for example, small μ leads to longer t. This provides us with some intuitions concerning the sufficiently large t for approximating the limiting mean frequency spectrum. Figure 3(b) illustrates the difference between the finite-time mean frequency spectrum and the limiting mean frequency spectrum as a function of t, for large t, t∈[15,35] and for j=1,2, where lines represent the true difference and markers represent the asymptotic approximation by formula (5). To emphasize the agreement for t large, this figure is plotted in semilogarithmic scale. We see that the true difference drops exponentially fast, and the asymptotic approximation is good for large t.Comparison between the finite-time mean frequency spectrum and the limiting mean frequency spectrum of the infinite-allele birth and death process. (a)ψ1j(t) and ψj for 1≤j≤15, t=15,20,28 and ∞. (b) Difference between ψ1j(t) and ψj as a function of t, t∈[15,35] for j=1,2. Lines represent the true difference and markers represent asymptotic approximations.
(a)
(b)Given the observed long-term mean frequency spectrum, the parametersθ of the IAMBP, such as α,μ in the birth and death process, can be estimated by equating the observed long-term mean frequency spectrum ψobs from the sample to the expected limiting mean frequency spectrum ψexp from formula (3) and solving for the process parameters. In the case of the birth and death process, we may estimate α and μ for example by solving
(11)(A2B2Γ(j1)Γ(k1+1+λc)×F(A2B2j1+1,1+λc;j1+1+λc;A2B2))×(A2B2Γ(k1)Γ(j1+1+λc)×F(A2B2k1+1,1+λc;k1+1+λc;A2B2))-1=ψobs(j1)ψobs(k1)(A2B2Γ(j2)Γ(k2+1+λc)×F(A2B2j2+1,1+λc;j2+1+λc;A2B2))×(A2B2Γ(k2)Γ(j2+1+λc)×F(A2B2k2+1,1+λc;k2+1+λc;A2B2))-1=ψobs(j2)ψobs(k2)
for positive integers j1≠k1, j2≠k2, where λ/c and A2/B2 are both functions of α and μ.There is no explicit solution for such estimator, but numerical search according to some criteria is feasible. Another possibility is to minimize the distance (such as thel2 norm) between the observed long-term mean frequency spectrum and the expected limiting mean frequency spectrum, that is, θ^=argminθ∥ψobs-ψexp(θ)∥2.The estimated parameters can be used to check the goodness of fit of the IAMBP model. Another interesting problem is to test whether two sets of parameters are identical, given two observed mean frequency spectra. A simple approach is to use Pearson’sχ2 test, such as in Kimmel and Mathaes [4]. However, there may be restrictions to applying the χ2 test, such as small cell counts and inappropriateness due to the finite length of the observed spectrum. This motivates us to develop an interval estimator for the IAMBP parameters.


## 3. Variance of the Frequency Spectrum

Moment estimators based on the mean frequency spectrum only give point estimates of the process parameters. In order to quantify the uncertainty of point estimates, an interval estimator is needed, which requires more information about the distribution of the statisticαt(j). First, it can be seen that [2]
(12)αt(j)=I0,j(t)+∑n=1Nt‍∑k=1Un‍In,k,j(t-Tn),
where T1,T2,… are the successive split times of the process, I0,j(t), In,k,j(t) are two indicators, and I0,j(t)=1 if there are j individuals alive at time t carrying the parental allele, and In,k,j(t)=1, for n,k≥1 if the kth individual born at time Tn (Tn<t) mutates to a novel allelic type and further produces j individuals carrying this allele t time units later. Nt is the number of split times in (0,t], and Un is the number of offspring produced at time Tn. Obtaining the distribution of αt(j) is not elementary. However, it may still be possible to define a confidence interval (CI) based on the first and second moments of αt(j).Letηi,t(j)=Vari(αt(j)) be the variance frequency spectrum; by the law of total variance and independence between the indicators in the expression of αt(j) (details in Appendix B), we have
(13)ηi,t(j)=qij(t)[1-qij(t)]+im2μ2[C(t)+(λ+a)eλt∫0t‍e-λxC(x)dx]+iamμeλt∫0t‍e-λxq1j(x)dx+ia(σ2-m)μ2eλt×∫0t‍e-λxq1j2(x)dx,
where
(14)C(t)=a∫0t‍[q1j2(x)+(σ2+m2)β12(x)]e-a(t-x)dx-[a∫0t‍q1j(x)e-a(t-x)dx]2-[am∫0t‍β1(x)e-a(t-x)dx]2.
In Expression (14), β1(x)=aeλx∫0x‍e-λuq1j(u)du, and σ2 is the variance of the offspring distribution, regardless of the allelic types.Similarly as in Expression (3), we may define a limiting variance frequency spectrum ξj=limt→∞ηi,t(j)/(Ei[Kt])2. Expression (13) is complicated and usually does not assume an explicit form, even for the special case of the birth and death process. Therefore, we will only give numerical solutions for the finite-time variance frequency spectrum. Figure 4 shows an example of the “2σ”-bands of the finite-time frequency spectrum for the infinite-allele birth and death process with i=100, α=0.25, μ=10-4, a=1, and t=28. To emphasize the tail probabilities, we draw this plot in semilogarithmic scale.Figure 4
“2σ”-bands in semi-logarithmic scale of the finite-time (t=28) frequency spectrum for the infinite-allele birth and death process with i=100, a=1, α=0.25, and μ=10-4.From the finite-time variance frequency spectrum, it is possible to define a CI[θl,θu] where the upper and lower bounds can be written as
(15)argminθ∥ψobs(t)-[ψexp(θ,t)±2ξexp(θ,t)]∥2.
This CI is useful for checking model validity and for testing whether two observed mean frequency spectra are from the same IAMBP model.


## 4. Simulation Study

We perform a simulation study of the birth and death process to illustrate the finite-time mean and variance frequency spectra. First we generate samples (genealogical trees) from an IAMBP with birth and death offspring distribution starting from100 individuals carrying the same parental allele. Due to memory restrictions caused by forward simulation, we limit our simulations to 12 generations and a relatively large mutation probability μ=0.01. The other parameters of the process are set to be a=1 and α=0.25. At time t=2, we record the number of alleles αt(j) represented by j copies, for j=1,2,…. Repeating the simulation 1000 times, we then obtain the simulated finite-time mean and variance frequency spectra from the replicates.Figure5 shows side-by-side bar plots of the simulated and expected finite-time mean and variance frequency spectra, for j=1,…,10. We plot the mean frequency spectrum and the variance frequency spectrum in semi-logarithmic scale to emphasize the tail probabilities. In each bar plot, the first black bar represents the expected finite-time mean frequency spectrum ψj or variance frequency spectrum ξj. The remaining ten white bars represent ten replicates of the simulated finite-time mean or variance frequency spectrum as described above. We see that for some classes, the expected mean or variance frequency spectrum is slightly different from the simulated spectrum. Beside sampling bias, this may be caused by the small scale of the simulations. For small mutation probability μ, we have to set large initial population size i and a long time t to obtain acceptable values of ψi,t(j) and ξi,t(j) from the simulated genealogical trees. We note that if one tries to use a naive method to calculate the variance frequency spectrum, that is, assume the proportions ψj of alleles having j representatives to be the mean of some independent Bernoulli random variables (they are not independent) and employ ξj=ψj(1-ψj), such method performs much worse than ξj based on the derivation of the variance frequency spectrum.Comparison between the simulated and expected finite-time (t=2) frequency spectra of the infinite-allele birth and death process with i=100, a=1, α=0.25, and μ=0.01. In each class j, 1≤j≤10, the first black bar represents the expected frequency spectrum and the rest 10 white bars represent 10 replicates of the simulated frequency spectrum. (a) mean frequency spectrum ψj in semi-logarithmic scale; (b) variance frequency spectrum ξj in semi-logarithmic scale.
(a)
(b)


## 5.  Summary

In this paper, we consider the frequency spectrum of the IAMBP of Pakes [2]. We develop an explicit expression for the limiting mean frequency spectrum for the special case of the birth and death process, which can be stated in terms of the hypergeometric function. We also derive an asymptotic expression for the rate of convergence of the finite-time mean frequency spectrum to the limiting mean frequency spectrum and illustrate the convergence using the birth and death process. We further state and prove a theorem concerning the variance frequency spectrum of the IAMBP, which helps to quantify uncertainty in parameter estimation and hypothesis testing. We illustrate the results using simulations of the birth and death process case.As noted in the introduction, the frequency spectrum is similar to the Ewens’ sampling formula [5] in population genetics, since they both concern the count or frequency spectrum based on an infinite-allele model under neutral selection. However, they differ in several aspects. (1) Our frequency spectrum describes population property under a branching process, whereas the Ewens’ sampling formula describes allelic class count probabilities caused by a sampling procedure and further requires the sample size n to be small compared to the size of the whole population which is assumed constant. (2) Our results concerning the frequency spectrum only provide the first and second moments and not the distribution function of the proportion of alleles having a given number of copies at a specified time point, whereas the Ewens’ sampling formula gives the joint probability of all allelic classes. We also note that the Poisson-Dirichlet process [8] is usually used to describe the equilibrium behavior of the neutral infinite-allele model. Study of the relation between variations of the frequency spectrum under different models is of our future interest.The question of validity of the Wright-Fisher and Moran models of population genetics [3], as compared to stochastic population processes such as the IAMBP or O’Connell process [9], has importance for estimation of parameters based on genetic data. As an example, Cyran and Kimmel [10] compared estimates of the age of the Mitochondrial Eve based on the various versions of the Wright-Fisher model with those based on various branching process models. The outcomes showed differences of about 10–15%.



---

*Source: 963831-2013-03-17.xml*