Transition in the ancestral reproduction rate and its implications for the site frequency spectrum

Cheek and Johnston [5] considered a continuous-time Galton-Watson process where each individual reproduces at rate $r$. In a reproduction event, an individual is replaced by $k\geq 0$ individuals with probability $p_{k}$. Denote by $N_{t}$ the population size at time $t$ and $F_{t}(s)=\mathbb{E}[s^{N_{t}}]$ for $s\in[0,1]$ its generating function. Conditional on survival up to time $T$, i.e., $N_{T}>0$, an individual is sampled uniformly at random from the population at time $T$. They studied the rate of reproduction events along the ancestral lineage of the sampled individual. The rate is biased towards large reproduction events and is inhomogeneous along the ancestral lineage. Specifically, they show the following theorem.

Igelbrink and Ischebeck [13] generalized this result to the case when the reproduction rates are age-dependent.

In general, one can also take a sample of size $n$ from the population. When we trace back the ancestral lineages in the sample, they will coalesce when individuals find their common ancestors, giving us a coalescent tree. A central question in population genetics is to understand this coalescent tree and the distribution of the mutations along the tree.

In this paper, we consider sampling $n$ individuals uniformly at random from a supercritical birth and death process with birth rate $\lambda$ and death rate $\mu$, where $\lambda>\mu$. We write $r=\lambda-\mu>0$ for the net growth rate. When there is a birth event, we assume the number of mutations that the child acquires has a Poisson distribution with mean $\nu$. The process starts with one individual and runs for $T$ units of time. As usual, we write $N_{t}$ for the population size at time $t$ and $F_{t}(s)=\mathbb{E}[s^{N_{t}}]$ for $s\in[0,1]$ for its generating function. Conditional on $N_{T}\geq n$, we take a uniform sample of size $n$.

We give the construction of the coalescent tree and the mutations in the following proposition, with the explanation given in Section 3. It is easier to construct the tree backward in time. That is, we trace back the lineages of the sampled individuals and track the time for individuals to find their common ancestor. If there are $n$ individuals in the sample, then there are $n-1$ coalescent times $H_{i,n,T}$ for $1\leq i\leq n-1$. We should emphasize that the construction of the coalescent tree without mutations is obtained by Harris, Johnston, and Roberts [12] and Lambert [17] using two different approaches. The approach described in steps 1 and 2 of construction below comes from Lambert [17].

We are interested in the case when the sample size is much smaller than the population size, so the sampling probability $Y_{n,T}$ is close to 0. To interpret the time-inhomogeneity, $F_{t}(1-y)$ is the probability conditional on $Y_{n,T}=y$, that an individual born $t$ time units before the sampling will have no descendants alive in the sample. As $t\rightarrow 0$, we have $F_{t}(1-y)\rightarrow 1-y\approx 1$, which is the probability that the newborn is not sampled. As $t\rightarrow\infty$, we have $F_{t}(1-y)\rightarrow\mu/\lambda<1$, which is the probability that the subtree generated by the newborn goes extinct. Therefore, we see more mutations along the branches near the bottom part of the tree.

To compare our result with that of Cheek and Johnston [5], observe that the density function of $Y_{n,T}$ with $n=1$ is

$$f_{Y_{1,T}}(y)=\frac{\delta_{T}}{(y+\delta_{T}-y\delta_{T})^{2}}=\frac{r(\lambda e^{rT}-\mu)}{(y\lambda e^{rT}-y\lambda+r)^{2}}.$$

Using the formula for $F_{t}(\cdot)$, one can compute the density function of the random variable $S$ in [5] to be

$$\frac{F_{T}^{\prime}(s)}{1-F_{T}(0)}=\frac{r(\lambda e^{rT}-\mu)}{\left((1-s)\lambda e^{rT}-(1-s)\lambda+r\right)^{2}}.$$

Therefore, the random variable $S$ agrees with $1-Y_{1,T}$. To interpret $r_{l}(s,t)$ in [5] correctly, note that the total rate of birth and death events is $\lambda+\mu$ and the probability of a size two reproduction event, i.e., a birth event, is $p_{2}=\lambda/(\lambda+\mu)$. Therefore, we have

$$r_{2}(s,t)=(\lambda+\mu)\cdot 2\cdot\frac{\lambda}{\lambda+\mu}F_{T-t}(s)=2\lambda F_{T-t}(s).$$

The extra constant 2 comes from the fact that only the children get mutations in our model, and the discrepancy between $T-t$ and $t$ is because time goes forward in [5] and backward in Proposition 1.1.

A commonly used summary statistic for mutational data is the site frequency spectrum. The site frequency spectrum for $n$ individuals at time $T$ consists of $(M^{1}_{n,T},M^{2}_{n,T},\dots,M^{n-1}_{n,T})$, where $M^{i}_{n,T}$ is the number of mutations inherited by $i$ individuals. Durrett [9] considered a supercritical birth and death process where mutations occur at a constant rate along the branches with rate $\nu$. He considered sampling the entire population and obtained that as the sampling time $T\rightarrow\infty$, the number of mutations affecting at least a fraction $x$ of the population is approximately $\nu/\lambda x$. Similar results also appear in [2, 21]. Durrett [9] also considered taking a sample of size $n$ and showed that as the sampling time $T\rightarrow\infty$, for $2\leq k\leq n-1$

$$\lim_{T\rightarrow\infty}\mathbb{E}[M^{k}_{n,T}]=\frac{n\nu}{k(k-1)}.$$

(3)

Gunnarsson, Leder, and Foo [10] considered models where the mutations occur at the times of reproduction events or at a constant rate along the branches. In both models, they computed the exact site frequency spectrum when the sample consists of the entire population or the individuals with an infinite line of descendants, at both a fixed time and a random time when the population size reaches a fixed level. Lambert [16] and Lambert and Stadler [18] developed a framework to study the site frequency spectrum when each individual is sampled independently with probability $p$, conditioned to have $n$ individuals. Following this approach, Delaporte, Achaz, and Lambert [6] computed $\mathbb{E}[M^{k}_{n,T}]$ for critical birth and death process, for both fixed $T$ and for an unknown $T$ with some improper prior. Dinh et al. [8] computed $\mathbb{E}[M^{k}_{n,T}]$ for supercritical birth and death process. As noted by Ignatieva, Hein, and Jenkins [14], if we take the sampling probability $p\rightarrow 0$ in [8], then one can recover the $1/k(k-1)$ shape in the site frequency spectrum as in formula (3).

Beyond the expected value of $M^{k}_{n,T}$, there are other works on its asymptotic behavior. Gunnarsson, Leder, and Zhang [11] considered a supercritical Galton-Watson process and sampled the entire population. They proved a strong law of large numbers for the site frequency spectrum at deterministic times and a weak law of large numbers at random times when the population size reaches a deterministic level. Schweinsberg and Shuai [20] considered a sample of size $n$ from supercritical or critical birth and death processes. As the sample size and the population size go to infinity appropriately, they established the asymptotic normality for the length of the branches that support $k$ leaves. These results can then be translated into the asymptotic normality of the site frequency spectrum, assuming that mutations occur at a constant rate along the branches.

In this paper, since we assume mutations occur at the times of reproduction events, the site frequency spectrum is characterized by the number of reproduction events. Indeed, let $R^{k}_{n,T}$ (resp. $R^{\geq 2}_{n,T}$) be the number of reproduction events in which the child has $k$ (resp. at least 2) descendants in the sample. Then conditional on $(R^{1}_{n,T},R^{2}_{n,T},\dots,R^{n-1}_{n,T})$, $M^{i}_{n,T}$ has a Poisson distribution with mean $\nu R^{i}_{n,T}$, and the $M^{i}_{n,T}$ are independent of one another. Therefore, we will focus on $R^{k}_{n,T}$ and $R^{\geq 2}_{n,T}$. We show the following theorem.

Under a stronger assumption on the sample size and the population size, we can establish a central limit theorem for $R^{\geq 2}_{n,T_{n}}$. The corresponding $M^{\geq 2}_{n,T_{n}}=\sum_{i=1}^{n-1}M^{i}_{n,T_{n}}$ represents the number of mutations shared by at least two individuals in the sample. The quantity $M^{\geq 2}_{n,T_{n}}$ has been applied to estimate the growth rate of a tumor in [15]. We believe an analogous result should hold for $R^{k}_{n,T_{n}}$, but the covariance computation is more involved.

The hypothesis is similar to (4), and the $1/k(k-1)$ site frequency spectrum is consistent with (3). Conditional on $L^{k}_{n,T_{n}}$, the number of mutations inherited by $k$ individuals has a Poisson distribution with mean $\widetilde{\nu}L^{k}_{n,T_{n}}\approx n\widetilde{\nu}/(r_{n}k(k-1))=n\lambda\nu/(r_{n}k(k-1))$, consistent with Theorem 1.2.

Johnson et al. showed a result similar to Corollary 1.2. The quantity $M^{in}_{n}$ therein corresponds to $M^{\geq 2}_{n,T_{n}}$. In the special case where the rates $\lambda,\mu$ and $\widetilde{\nu}$ are constants, Corollary 2 of [15] reduces to the following theorem.

The hypothesis is the same as (5), and the result is the same when we plug in $\widetilde{\nu}=\lambda\nu$. Therefore, our results establish that the assumption that mutations occur at a constant rate along the branches is not essential to the results in [15, 20], and similar results hold when mutations occur only at the birth times, which may be a more natural assumption for some applications.

Another closely related summary statistic for the mutational data is the allele frequency spectrum. For the allele frequency spectrum, the individuals are partitioned into families where each family consists of individuals carrying the same mutations. Then the allele frequency spectrum consists of $(A^{1}_{n,T},A^{2}_{n,T},\dots,A^{n}_{n,T})$ where $A^{i}_{n,T}$ is the number of families of size $i$. Richard [19] considered the allele frequency spectrum of the population in a birth and death process with age-dependent death rates and mutations at birth times. He computed the expected allele frequency spectrum and proved a law of large numbers for the allele frequency spectrum. Champagnat and Lambert considered the case where the death rate is age-dependent and the mutations occur at a constant rate along the branches. They computed the expected allele frequency spectrum when the sample consists of the entire population in [3] and studied the size of the largest family and the age of the oldest family in [4].

In this section, we describe the coalescent tree, together with the red and blue mutations in the Bernoulli sampling scheme. Lambert and Stadler [18] obtained the construction of the subtree. We outline the construction here for readers’ convenience. Recall that the contour process $X^{(T)}$ makes jumps at rate $\lambda$. That is, birth events occur at rate $\lambda$ as we trace back the lineage. Consider a birth event occurring $t$ time units before the sampling time. Recalling the definition of $N_{t}$ from the introduction, the probability that this birth event does not give rise to a new leaf in the coalescent tree is

$$q(y,t)=\mathbb{P}\left(N_{t}=0\right)+\sum_{k=1}^{\infty}\mathbb{P}\left(N_{t}=k\right)(1-y)^{k}=\mathbb{E}\left[(1-y)^{N_{t}}\right]=F_{t}(1-y).$$

(6)

Therefore, red birth events occur with rate $\lambda(1-q(y,t))$, and blue birth events occur with rate $\lambda q(y,t)$, independently of each other. For a coupling argument in Section 4.2, red and blue birth events are represented using a Poisson point process. To be more precise, for any positive continuous functions $L\leq U$ on $[a,b]$, we let $\mathcal{A}(L,U,[a,b])=\{(t,x)\in\mathbb{R}^{2}:a\leq t\leq b,L(t)\leq x\leq U(t)\}$ be the area between $L$ and $U$ on $[a,b]$. In the special case where $L\equiv 0$ on $[a,b]$, we abbreviate the notation as $\mathcal{A}(U,[a,b])$. We also let $\{\Pi_{i}\}_{i\geq 0}$ be independent Poisson point process on $\mathbb{R}^{2}$ with unit intensity. Then, for the Bernoulli$(y)$ sampling, one can obtain the coalescent tree and the mutations in the following way:

1. 1.

   We get an empty sample with probability $\mathbb{E}[(1-y)^{N_{T}}]=F_{T}(1-y)$. If the sample is non-empty, then we have an initial branch $H_{0,y,T}$ of length $T$.

2. 2.

   Let $\{\pi_{i,y}\}_{i\geq 0}$ be independent Poisson processes with rate $\lambda(1-q(y,t))$ on $[0,T]$ defined by $\pi_{i,y}(t)=\Pi_{i}(\mathcal{A}(\lambda q(y,\cdot),\lambda,[0,t]))$ for $0\leq t\leq T$.

3. 3.

   Having constructed $H_{0,y,T},\dots,H_{i,y,T}$, if $\pi_{i,y}(T)>0$, then the $i$th branch length is the time of the first jump in $\pi_{i,y}$, i.e. $H_{i,y,T}=\inf\{t\geq 0:\pi_{i,y}(t)>0\}$. Otherwise, we have found all the individuals in the sample. Note that $H_{i,y,T}$ is also the time for the red reproduction event on the $i$th individual.

4. 4.

   Blue birth events occur along the branches at rate $\lambda q(y,t)$. That is, the number of mutations along the $i$th branch up to time $t$ is $\Pi_{i}(\mathcal{A}(\lambda q(y,\cdot),[0,t]))$.

We compute $q(y,t)=F_{t}(1-y)$ here and refer the reader to equation (7) of [17] for the formula for the density of $H_{i,y,T}$. We define

$$\delta_{t}=\mathbb{P}(\pi_{i,1}(t)=0)=\exp\left(-\int_{0}^{t}\lambda\mathbb{P}(N_{s}>0)ds\right).$$

(7)

By Remark 3.1, conditional on $N_{t}>0$, $N_{t}$ has a geometric distribution with

$$\mathbb{P}(N_{t}\geq k+1|N_{t}\geq k)=\mathbb{P}(\pi_{i,1}(t)>0)=1-\delta_{t},\qquad\text{ for }k\geq 1.$$

Therefore, we have

$$e^{rt}=\mathbb{E}[N_{t}]=\mathbb{P}(N_{t}>0)/\delta_{t}.$$

(8)

Combining (7) and (8), we get

$$\delta_{t}=\exp\left(-\int_{0}^{t}\lambda e^{rs}\delta_{s}ds\right).$$

Solving this equation for $\delta_{t}$ and using the initial condition $\delta_{0}=1$, we get the following formula for $\delta_{t}$.

$$\delta_{t}=\frac{r}{\lambda e^{rt}-\mu},$$

which agrees with the formula of $\delta_{t}$ in (1). This formula is also found in equation (4) of [17]. To compute $q(y,t)$, by (8) and Remark 3.1 again, we have

	$\displaystyle q(y,t)$	$\displaystyle=F_{t}(1-y)$
		$\displaystyle=\mathbb{P}(N_{t}=0)+\mathbb{P}(N_{t}\geq 1)\mathbb{E}\left[(1-y)^{N_{t}}|N_{t}\geq 1\right]$
		$\displaystyle=1-\delta_{t}e^{rt}+\delta_{t}e^{rt}\cdot\frac{\delta_{t}(1-y)}{\delta_{t}+y-y\delta_{t}}.$		(9)

Lambert [17] showed that the uniform$(n)$ sampling can be realized by taking $y$ from the density (1) and then performing the Bernoulli($y$) sampling conditioned to have size $n$. Therefore, for the uniform$(n)$ sampling, one can obtain the coalescent tree and the mutations in the following way:

1. 1.

   Take the sampling probability $Y_{n,T}$ from the density (1).

2. 2.

   Conditional on $Y_{n,T}=y$, for $i=1,2,\dots,n-1$, let $\pi_{i,y}$ be the Poisson process be defined in Section 3.1, conditional on $\pi_{i,y}(T)>0$.

3. 3.

   For $i=1,2,\dots,n-1$, the $i$th branch length $H_{i,n,T}$ is the time of the first jump in $\pi_{i,y}$, conditional on $\pi_{i,y}(T)>0$. That is,

   $$\mathbb{P}(H_{i,n,T}>t|Y_{n,T}=y)=\mathbb{P}(\pi_{i,y}(t)=0|\pi_{i,y}(T)>0).$$

   The density for $H_{i,n,T}$ is given by (2).

4. 4.

   Blue birth events occur along the branches at rate $\lambda q(y,t)$.

Approximations for the sampling probability $Y_{n,T}$ and the branch lengths $\{H_{i,n,T}\}_{i=1}^{n-1}$ when $n$ and $T$ are large are obtained in [15].

1. 1.

   Choose $W$ from the exponential density

   $$f_{W}(w)=e^{-w},\qquad w\in(0,\infty).$$

2. 2.

   Choose $U_{i}$ for $i=1,\dots,n-1$ independently from the logistic distribution, with density

   $$f_{U_{i}}(u)=\frac{e^{u}}{\left(1+e^{u}\right)^{2}},\quad u\in(-\infty,\infty).$$

3. 3.

   Let $Y=n\delta_{T}/W$.

4. 4.

   Let ${H}_{i}=T-(\log n+\log(1/W)+U_{i})/r$ for $i=1,\dots,n-1$.

Note that the $H_{i}$’s (hence the quantities defined below) depend on $n$, although this dependence is not explicit in the notation. Using this approximation, we can, therefore, define an approximation for the number of reproduction events with at least two descendants in the sample on the $i$th branch as

$$R^{\geq 2}_{0}=\Pi_{0}(\mathcal{A}(\lambda q(n\delta_{T}/W,\cdot),[H_{1},T])),$$

and

$$R^{\geq 2}_{i}=\Pi_{i}(\mathcal{A}(\lambda q(n\delta_{T}/W,\cdot),[H_{i+1},H_{i}]))+\mathbbm{1}_{\{H_{i+1}\leq H_{i}\}},\qquad\text{ for }1\leq i\leq n-2.$$

An approximation for the total number of reproduction events with at least two descendants in the sample, excluding the $0$th branch, is

$$R^{\geq 2}=\sum_{i=1}^{n-2}R^{\geq 2}_{i}.$$

Similarly, for $k\geq 2$, we define

$$R^{k}_{0}=\Pi_{0}\left(\mathcal{A}\left(\lambda q(n\delta_{T}/W,\cdot),\left[\max_{1\leq j\leq k-1}H_{j},H_{k}\right]\right)\right),$$

and for $1\leq i\leq n-k-1$

$\displaystyle R^{k}_{i}=\Pi_{i}\left(\mathcal{A}\left(\lambda q(n\delta_{T}/W,\cdot),\left[\max_{i+1\leq j\leq i+k-1}H_{j},H_{i}\wedge H_{i+k}\right]\right)\right)+\mathbbm{1}_{\{\max_{i+1\leq j\leq i+k-1}H_{j}\leq H_{i}\leq H_{i+k}\}},$

and

$\displaystyle R^{k}_{n-k}=\Pi_{i}\left(\mathcal{A}\left(\lambda q(n\delta_{T}/W,\cdot),\left[\max_{n-k-1\leq j\leq n-1}H_{j},H_{n-k}\right]\right)\right)+\mathbbm{1}_{\{\max_{n-k-1\leq j\leq n-1}H_{j}\leq H_{n-k}\}},$

An approximation for the total number of reproduction events inherited by $k$ individuals, excluding the $0$th and the $(n-k)th$ branches, is

$$R^{k}=\sum_{i=1}^{n-k-1}R^{k}_{i}.$$

In the rest part of the paper, the sampling time $T=T_{n}$ depends on $n$. The errors in the approximation resulting from taking the limit as $n\rightarrow\infty$ are bounded by Lemmas 4 and 5 in the supplementary material of [15]. We combine these lemmas into Lemma 4.1. The event $A_{n}$ in Lemma 4.1 is the intersection of the events $Q_{n,T_{n}}=Q_{n,\infty}\leq h(n):=n^{1/2}e^{rT_{n}/2}$ and $\widetilde{Q}_{n,\infty}=1/W$ in [15]. The event $B_{n}$ is the intersection of $Q_{n,T_{n}}=Q_{n,\infty}\leq h(n):=n^{1/2}(\log n)^{-1/3}e^{r_{n}T_{n}/3}$ and $\widetilde{Q}_{n,\infty}=1/W$ in [15]. The choice of $h(n)$ comes from the proof of Lemma 6 of [15].

These error bounds allow us to bound the error in the approximations for $R^{\geq 2}_{n,T_{n}}$ and $R^{k}_{n,T_{n}}$, with the aid of another technical lemma.