Coalescence in the diffusion limit of a Bienayme-Galton-Watson branching process

TL;DR

This paper develops a statistical framework for estimating the time since the most recent common ancestor in populations modeled by a supercritical Bienayme-Galton-Watson process, using limited information and diffusion approximations.

Contribution

It introduces a novel method for estimating the MRCA time and initial population size based on minimal data and diffusion limits, extending previous analyses.

Findings

Derived a probability formula for a sample descending from a single ancestor.

Provided confidence intervals for initial population size based on observed data.

Constructed a joint likelihood surface for estimating population size and MRCA time.

Abstract

We consider the problem of estimating the elapsed time since the most recent common ancestor of a finite random sample drawn from a population which has evolved through a Bienayme-Galton-Watson branching process. More specifically, we are interested in the diffusion limit appropriate to a supercritical process in the near-critical limit evolving over a large number of time steps. Our approach differs from earlier analyses in that we assume the only known information is the mean and variance of the number of offspring per parent, the observed total population size at the time of sampling, and the size of the sample. We obtain a formula for the probability that a finite random sample of the population is descended from a single ancestor in the initial population, and derive a confidence interval for the initial population size in terms of the final population size and the time since initiating the process. We also determine a joint likelihood surface from which confidence regions can be determined for simultaneously estimating two parameters, (1) the population size at the time of the most recent common ancestor, and (2) the time elapsed since the existence of the most recent common ancestor.

Figures

Equations

$$\mathbb{E}(S_{\alpha})=\lambda,\qquad\mbox{$\rm{Var}$\,}(S_{\alpha})=\sigma^{2},$$

$$M(t+1)=\sum _{\alpha =1}^{M(t)}{S}_{\alpha },i=1,2.$$

$${\kappa }_0=\frac{2m_0\log \lambda }{{\sigma }^2},$$

$$s=t\log \lambda ,K(s)=\frac{2M(t)\log \lambda }{{\sigma }^2}.$$

$$\frac{\partial f_K(\kappa ;s)}{\partial s}=-\frac{\partial (\kappa f_K(\kappa ;s))}{\partial \kappa }+\frac{{\partial }^2(\kappa f_K(\kappa ;s))}{\partial {\kappa }^2},$$

$$f_{K}(\kappa;s)\,d\kappa=\mbox{$\rm{Prob}$\,}(\kappa\leq K(s)<\kappa+d\kappa).$$

$$\begin{array}{rl}\delta s& =\log \lambda \\ \delta K(s)& =\frac{2\log \lambda }{{\sigma }^2}(M(t+1)-M(t)).\end{array}$$

$$\begin{array}{rl}\mathbb{E}(\delta K(s)\mid K(s)=\kappa )& =a(\kappa )\delta s+o(\delta s),\\ \mathbb{E}(\delta K(s{)}^2\mid K(s)=\kappa )& =b(\kappa )\delta s+o(\delta s),\\ \mathbb{E}(\delta K(s{)}^k\mid K(s)=\kappa )& =o(\delta s),k\ge 3,\end{array}$$

$$\frac{\partial f_K(\kappa ;s)}{\partial s}=-\frac{\partial }{\partial \kappa }\left(a(\kappa )f_K(\kappa ;s)\right)+\frac{1}{2}\frac{{\partial }^2}{\partial {\kappa }^2}\left(b(\kappa )f_K(\kappa ;s)\right).$$

$$\mathbb{E}(\delta K(s)\mid K(s)=\kappa )$$

$$\mathbb{E}(\delta K(s{)}^2\mid K(s)=\kappa )$$

$$f_K(\kappa ;s|{\kappa }_0)=\delta (\kappa )\exp \left(-\frac{{\kappa }_0}{1-e^{-s}}\right)+$$

$$f_K(\kappa ;0|{\kappa }_0)=\delta (\kappa -{\kappa }_0).$$

$$\varphi (\theta ,s)={\int }_0^{\infty }{e}^{-\theta \kappa }{f}_K(\kappa ;s|{\kappa }_0)d\kappa .$$

$$\frac{\partial \varphi (\theta ,s)}{\partial s}+\theta (\theta -1)\frac{\partial \varphi (\theta ,s)}{\partial \theta }=0.$$

$$\varphi (\theta ,0)=e^{-\theta {\kappa }_0}.$$

$$0=\frac{\partial \varphi (\theta ,s)}{\partial s}+\frac{d\theta }{ds}\frac{\partial \varphi (\theta ,s)}{\partial \theta }.$$

$$\frac{d\theta }{ds}=\theta (\theta -1),$$

$$\theta (s)=\frac{\theta (0)}{\theta (0)-(\theta (0)-1)e^s}.$$

$$\theta (0)=\frac{\theta e^s}{1+\theta (e^s-1)},$$

$$\varphi (\theta ,s)=\exp \left\{-\frac{{\kappa }_0\theta e^s}{1+\theta (e^s-1)}\right\}.$$

$$\varphi (\theta )=\exp \left\{-\frac{A\theta }{1+B\theta }\right\},$$

$$f_K(\kappa )=\delta (\kappa )\exp \left(-\frac{A}{B}\right)+\exp \left(-\frac{A+\kappa }{B}\right){\left(\frac{A}{\kappa B^2}\right)}^{\frac{1}{2}}{I}_1\left(\frac{2(A\kappa {)}^{\frac{1}{2}}}{B}\right).$$

$$f_{\mathrm{Feller}}(z;\xi )=\delta (z)e^{-\xi }+\xi z^{-\frac{1}{2}}{e}^{-\xi (1+z)}{I}_1\left(2\xi z^{\frac{1}{2}}\right),$$

$$f_K(\kappa ;s|{\kappa }_0)=\frac{1}{{\kappa }_0{e}^s}{f}_{\mathrm{Feller}}\left(\frac{\kappa }{{\kappa }_0{e}^s};\frac{{\kappa }_0}{1-e^{-s}}\right).$$

$$\displaystyle\mbox{$\rm{Prob}$\,}({\cal E}_{n}(\kappa_{0},s),K(s)\in[\kappa,\kappa+d\kappa)|K(0)=\kappa_{0})$$

$$w=\frac{(\kappa {\kappa }_0{e}^s{)}^{\frac{1}{2}}}{e^s-1}=\frac{(\kappa {\kappa }_0{)}^{\frac{1}{2}}}{2\sinh \frac{1}{2}s}.$$

$$K_1(0)=x_0{\kappa }_0,K_2(0)=(1-x_0){\kappa }_0,0\le x_0\le 1.$$

$$K(s)=K_1(s)+K_2(s),X(s)=\frac{K_1(s)}{K_1(s)+K_2(s)}.$$

$$f_{K,X}(\kappa ,x;s|{\kappa }_0,x_0)=\kappa f_K(x\kappa ;s|x_0{\kappa }_0)f_K((1-x)\kappa ;s|(1-x_0){\kappa }_0),$$

Full text

1 Introduction

Suppose one observes the current size of a population which has been evolving through a Bienaymé-Galton-Watson (BGW) branching process since some unknown time in the past. Is it possible to infer the time that has elapsed since the most recent common ancestor (MRCA) of the current population, and of a random sample taken from the current population?

This paper is a continuation of an earlier paper [5] which gives a partial answer to this question in the form of a joint estimate of the elapsed time since coalescence to the MRCA and the corresponding size of the entire population at the time of the MRCA. Our new results extend the previous calculation to the MRCA of finite random samples of the population, and to a proposed procedure for obtaining confidence regions in the joint parameter space of the time since coalescence and the population size at the time of coalescence.

There have been many studies of coalescence in BGW process [17, 20, 13], specifically for critical [9, 2] sub-critical [16, 18] and supercritical processes. In the current paper we concentrate on supercritical processes in the near-critical limit, or, equivalently, the continuous-time diffusion limit. Closest to our treatment is that of [19], who derives a formula (with a minor correction in [15]) for the coalescent time for sample of size $n=2$ as a fraction of the time since initialisation of the BGW process with a known initial population size. Numerical simulation of an ensemble of BGW processes initiated from a single founder by [8] is consistent with this formula. The ultimate aim of Cyran and Kimmel’s paper was to compare the effectiveness of various population genetics models in dating the time since mitochondrial Eve (mtE). More recently O’Connell’s formula has been generalised to samples of size $n>2$ by [14].

Our approach differs from previous treatments in that we initialise the BGW process at the (a priori unknown) time of coalescence of a random sample of the current population, treat the time since coalescence and initial population size as parameters to be estimated, and determine confidence regions in the two-dimensional parameter space defined by these two parameters. Central to our approach is the solution by [11] to the forward Kolmogorov equation defining the diffusion limit of a BGW process. This approach is closer in spirit to the intention of Cyran and Kimmel’s aim of determining the time since mtE without prior knowledge of the time at which the process was initiated from a single founding individual. It also enables a direct comparison with the analogous solution by [21] to the problem of determining the distribution of the time since coalescence of a sample of size 2 in a Wright-Fisher (WF) population with exponential growth.

The main results of the paper are stated in three theorems. Theorem 1 in Section 3 gives a formula for the probability that a uniform random sample is descended from a single individual in the initial population of a BGW process, given the initial and final population sizes and the time since the process was initiated. Theorem 2 in Section 4 gives confidence intervals for the unknown initial population size in a BGW process when the final population size and time since initiation are known. These two theorems enable us to map out the likely initial conditions in terms of two scaled parameters arising from the diffusion limit; (1) a scaled population size at the time of initialisation of the process, denoted $\kappa_{0}$, and (2) the scaled time elapsed from initialisation of the process to observation of the current population, denoted $s$. Theorem 3 in Section 6, which is contingent on a technical conjecture stated in the Section 5, gives a universal asymptotic density function for the likely time since coalescence of a sample in the limit of large scaled population sizes.

The layout of this paper is as follows. The diffusion limit of a BGW process, notational conventions employed in subsequent sections, and a derivation of Feller’s solution for a specified initial population are detailed in Section 2. Arguments leading to Theorems 1 and 2 and the proofs of these theorems are in Sections 3 and 4 respectively. Section 3 is a generalisation of analogous derivations in [5], which concentrated on the entire population as opposed to a finite random sample. Section 5 is devoted to deriving a function ${\cal L}(s,\kappa_{0};\kappa)$ defined over the $(s,\kappa_{0})$, the contours of which are confidence regions for the scaled time since coalescence of a sample and the scaled population size at the time of coalescence. Here $\kappa$ is the scaled final population size at the time of sampling. By integrating out the initial population size $\kappa_{0}$, a comparison is made with the distributions derived by [21] for a sample of size $n=2$ taken from an exponentially growing WF population. The asymptotic $\kappa\rightarrow\infty$ limit is analysed in Section 6, and this enables us to compare our analytic formulae with a numerical simulation in Section 7. Results are discussed and conclusions drawn in Section 8. A technical appendix is devoted to details of calculations arising in Section 5.

2 The diffusion limit of a Bienaymé-Galton-Watson (BGW) branching process

Consider a Bienaymé-Galton-Watson (BGW) branching process consisting of a population of $M(t)$ haploid individuals reproducing in discrete, non-overlapping generations $t=0,1,2,\ldots$, with initial population size $M(0)=m_{0}$. The numbers of offspring per individual at each time step $t$ are identically and independently distributed (i.i.d.) random variables $S_{\alpha}$, $\alpha=1,\ldots,M(t)$, whose common distribution has mean and variance

[TABLE]

and finite moments to all higher orders. Thus

[TABLE]

We are interested in the supercritical case of the diffusion limit studied by Feller [11] in which the initial population $m_{0}$ becomes infinite while the growth rate simultaneously approaches the critical value $\lambda\rightarrow 1$ in such a way that the parameter

[TABLE]

remains fixed. As argued in Burden and Simon [5], this parameter has a straightforward physical meaning: Suppose the initial population is divided into two subpopulations of roughly equal size. If the initial population size $m_{0}$ is such that $\kappa_{0}>>1$, then with high probability descendant lineages of both subpopulations will survive as $t\rightarrow\infty$. On the other hand, if $\kappa_{0}<<1$ and the population does not go extinct, the entire surviving population as $t\rightarrow\infty$ will, with high probability, be descended from an individual in one of the two subpopulations. Thus $\kappa_{0}\approx 1$ sets a scale for inferring coalescent events. With this in mind, in the following lemma we introduce continuum quantities $K(s)$ and $s$ adapted from Burden and Simon [5, Eq. (36)].

Lemma 1**.**

Define a continuum time $s$ and continuous random population sizes $K(s)$ by

[TABLE]

The corresponding forward Kolmogorov equation for the density $f_{K}(\kappa;s)$ in the limit $\log\lambda\rightarrow 0$ and $m_{0}\rightarrow\infty$ such that $\kappa_{0}$ defined by Eq, (3) remains fixed is

[TABLE]

where

[TABLE]

Proof.

The derivation follows the formal method given in Ewens [10, Chapter 4] for obtaining a forward-Kolmogorov equation. Equating an increment of continuum to one time step in the discrete model, set

[TABLE]

In general, if

[TABLE]

for finite functions $a(\kappa)$ and $b(\kappa)$ as $\delta s\rightarrow 0$, then the forward Kolmogorov equation takes the form

[TABLE]

From Eqs. (1), (2) and (7), and setting $\kappa=2m\log\lambda/\sigma^{2}$ one obtains

[TABLE]

and

[TABLE]

Thus $a(\kappa)=\kappa$, $b(\kappa)=2\kappa$, and Eq. (5) follows. ∎

Note that the continuum scaling defined by Eqs. (3) and (4) is deliberately chosen so that all parameters defining the process are subsumed into a single parameter $\kappa_{0}$ which is present in the diffusion limit only via the initial conditions, and that the differential equation itself is parameter-free. However, while appropriate for current problem of exploring coalescence in the near-supercritical case $\lambda>1$, the scaling is not appropriate appropriate for the near-subcritical case $\lambda<1$ as the subsequent analysis assumes $\kappa_{0}>0$ and that the density $f_{K}(\kappa)$ is non-zero only for $\kappa\geq 0$.

The solution to the forward Kolmogorov equation for a BGW process was first solved, albeit with a different continuum scaling, by Feller [11, 12] by means of a Laplace transform. Summaries of Feller’s method of solution can be found in Cox and Miller [7, p236] and Burden and Wei [6]. For completeness, the following lemma derives the solution for a supercritical Feller diffusion using the continuum scaling and notation of the current paper.

Lemma 2**.**

The solution solution to the forward Kolmogorov equation, Eq. (5), corresponding to an initial scaled population $K(0)=\kappa_{0}>0$ is

[TABLE]

where $\delta(\cdot)$ is the Dirac delta function and $I_{1}(\cdot)$ is a modified Bessel function of the first kind. In particular,

[TABLE]

Proof.

Consider Laplace transform

[TABLE]

Applying the Laplace transform to both sides of Eq. (5) and carrying through standard manipulations [7, pp218, 219] gives

[TABLE]

This is a first order linear partial differential equation for $\phi(\theta;s)$ whose initial boundary condition, corresponding to Eq. (13), is

[TABLE]

The solution is found by observing that $\phi(\theta;s)$ is constant along characteristic curves in the $(s,\theta)$-plane defined by

[TABLE]

Comparing Eqs. (14) and (17) we see that the characteristic curves are solutions to the differential equation

[TABLE]

namely

[TABLE]

Thus, $\phi(s,\theta)$ can be found by tracing the characteristic curve back to the boundary point

[TABLE]

obtained by inverting Eq. (19). Combining this with the initial condition Eq. (16) gives

[TABLE]

It remains to invert the Laplace transform. The inverse Laplace transform of

[TABLE]

is derived in Cox and Miller [7, p236, 250] as111The continuous part of the density stated in Cox and Miller [7, p250] is missing a factor $\exp(-A/B)$. This appears to be an error.

[TABLE]

Setting $A=\kappa_{0}e^{s}$ and $B=e^{s}-1$ gives Eq. (2). ∎

Remark 1**.**

In the first term of Eq. (2) the coefficient of the delta-function is the probability of the population becoming extinct up to time $s$. See, for instance, Bailey [3, p206].

Remark 2**.**

*Eq. (2) can be written in terms of a 1-parameter family of normalised density functions222Note that Eq. (31) of Burden and Wei [6] defining the function $f_{\rm Feller}(\cdot)$ contains an error: the argument of the modified Bessel function should be $2\kappa_{0}z^{\frac{1}{2}}$, not $2\kappa_{0}z^{-\frac{1}{2}}$. *

[TABLE]

where $z,\xi\geq 0$, as

[TABLE]

3 Coalescence of a sample of $n$ individuals

Our aim is to estimate the time of the most recent common ancestor (MRCA) of a sample of $n$ individuals chosen randomly and uniformly from a population of scaled size $\kappa=(2m\log\lambda)/\sigma^{2}$, where $m$ is the population size at the “current” time $s=t\log\lambda$ since the coalescent event, which we set to be at time [math]. The population is assumed to have evolved via a BGW process with constant values of the parameters $\lambda$ and $\sigma^{2}$, which are assumed to be given. Thus $\kappa$ is to be thought of as a known input parameter, while $s$ is an unknown parameter which seek to estimate. In the process, the scaled population size $\kappa_{0}=2m_{0}\log\lambda/\sigma^{2}$ at the time of the MRCA will also be estimated. The analysis follows the same reasoning as earlier work of Burden and Simon [5] on estimating the time of the MRCA of an entire population, known as mtE in the case of the existing human population.

Lemma 3**.**

Define ${\cal E}_{n}(\kappa_{0},s)$ to be the event that a random uniform sample of $n$ individuals taken at time $s$ are descended from a single individual in the original population at time [math], given that the initial scaled population size was $\kappa_{0}$. Then the probability of the joint event that

${\cal E}_{n}(\kappa_{0},s)$* happens; and* 2. 2.

the population size at time $s$ is in the range $[\kappa,\kappa+d\kappa)$

is

[TABLE]

where

[TABLE]

Proof.

Consider the population to be divided into two types with scaled population sizes $K_{1}(s)$ and $K_{2}(s)$, such that each type evolves as an independent BGW process with initial conditions

[TABLE]

Introduce two random variables

[TABLE]

Here $K(s)$ is the total population size and $X(s)$ is the fraction of the population which is of the first type at time $s$. Since $K_{1}(s)$ and $K_{2}(s)$ are independent, the joint density corresponding to the initial conditions $K(0)=\kappa_{0}$, $X(0)=x_{0}$ is

[TABLE]

defined over the range $0\leq\kappa<\infty$, $0\leq x\leq 1$, $0\leq s<\infty$. The function $f_{K}(\cdot)$ is defined by Eq. (2) and the factor of $\kappa$ arises from the Jacobian of the transformation Eq. (29).

In the event ${\cal E}_{n}(\kappa_{0},s)$, define descendants of the common ancestor of the sample of size $n$ to be of type-1, and members of the population who are not descended from this individual to be of type-2. Setting $x_{0}=1/m_{0}$ in the density $f_{K,X}$ we have

[TABLE]

The initial factor $m_{0}$ accounts for the $m_{0}$ possibilities for the common ancestor of the sample, the integral ranges over all possible fractions $x$ of the final population that can be descended from that ancestor, and the factor $x^{n}$ is the probability that all $n$ individuals in the sample are of type-1, given that $X(s)=x$.

From Eqs. (2) and (30), and using the identity $\kappa\delta(\kappa x)=\delta(x)$ and the property of the modified Bessel function that $I_{1}(2w)=w+O(w^{2})$ as $w\rightarrow 0$, we obtain

[TABLE]

Substituting back into Eq. (31), noting that the $\delta(x)$ term in the first factor of Eq. (32) does not contribute because of the factor $x^{n}$ in the integrand, integrating out the $\delta(1-x)$ in the second factor, and taking the limit $m_{0}\rightarrow\infty$ gives, after some straightforward algebra,

[TABLE]

where $w$ is defined by Eq. (27). The integral can be evaluated with the aid of Wolfram alpha [22] as

[TABLE]

from which the required result, Eq. (26), follows. ∎

Theorem 1**.**

Define $u_{n}(s,\kappa_{0}|\kappa)$ to be the probability that the initial population of a BGW process contains an individual who is the single common ancestor of a random uniform sample of $n$ individuals taken at a later time $s$, given the initial and final population sizes are $\kappa_{0}$ and $\kappa$ respectively. Then

[TABLE]

where $w$ is defined by Eq. (27).

Proof.

By definition,

[TABLE]

The numerator follows from Lemma 3. The denominator, obtained from the continuous part of Eq. (2), is

[TABLE]

Substituting Eqs.(26) and (37) into Eq. (36) yields Eq. (35). ∎

Corollary 1**.**

[5]** The probability that the entire population at time $s$ in a BGW process is descended from a single individual in the initial population at time 0, given the initial and final population sizes are $\kappa_{0}$ and $\kappa$ respectively, is

[TABLE]

Proof.

From the ascending series for modified Bessel function [p375 of 1] we have $I_{n}(2w)=(w^{n}/n!)(1+O(1/n))$ as $n\rightarrow\infty$. Taking the limit of Eq. (35) leads to Eq. (38). ∎

The interpretation of Theorem 1 is illustrated by the black contours in Figure 1. Figures 1(a) and 1(b) show a contour map of the density in Eq. (35) in the $s$-$\kappa_{0}$ plane for a random sample of size $n=2$ taken from an observed final scaled population of size $\kappa=3$ and $\kappa=2000$ respectively. The map grades from a probability near 0 on the left edge of the plot corresponding to recent initialisation times to a probability near 1 at earlier times. If the BGW process were initiated at a time $s$ in the past with a population of size $\kappa_{0}$, the contour passing through the point $(s,\kappa_{0})$ gives the probability that the ancestral coalescence of the sample occurred after the initialisation of the process. Note that for larger values of $\kappa$ the coalescent event is fixed within a relatively narrow range.

4 Confidence interval for $\kappa_{0}$ given an observed $K(s)$

Superimposed in blue in Figure 1 are contours of the likelihood function that the population at time $s$ will not exceed the observed population $\kappa$ for a given starting population $\kappa_{0}$,

[TABLE]

where the function $f_{K}$ is defined by Eq. (2). Note that the integrand includes the point mass at $\eta=0$ which accounts for the possibility of extinction of the entire population. The contour marked $0.1$ connects starting configurations for which there is a 10% probability the population will undershoot the observed population at time $s$, and the contour marked $0.9$ connects starting configurations for which there is a 10% chance of overshooting.

The function $v(s,\kappa_{0}|\kappa)$ enables us to establish confidence interval for estimating the initial population $K(0)=\kappa_{0}$ given that a population generated by a BGW process is observed to be of size $K(s)$ at a known time $s$. This confidence interval is define in the following theorem.

Theorem 2**.**

Consider the inverse function $v^{-1}(p|s,\kappa)$ mapping the probability $p\in[0,1]$ back to an initial condition $\kappa_{0}$ in the interval $[0,\infty)$, thus,

[TABLE]

Then if a BGW process initiated at time 0 with population $\kappa_{0}$ evolves to a population $K(s)$ at time $s$, for any $0\leq p_{1}<p_{2}\leq 1$, the random interval

[TABLE]

has a probability $p_{2}-p_{1}$ of containing the starting population $\kappa_{0}$.

Proof.

It is clear that $v^{-1}$ is in general a decreasing function of $p$, with $v^{-1}(0|s,\kappa)=\infty$ and $v^{-1}(1|s,\kappa)=0$. Then

[TABLE]

The second last line follows because $v^{-1}$ is a decreasing function of $p$, and the final line follows because the cumulative distribution $v(s,\kappa_{0}|K(s))$ is uniformly distributed on $[0,1]$. ∎

Thus, for instance, the blue contours marked $v=0.1$ and $v=0.9$ define the edges of an 80% confidence interval on starting configurations $(s,\kappa_{0})$ for the observed final population $\kappa$. The small circle in each plot is a “median estimate” of the coalescence point which we define as the (scaled) starting time and population size for which there is a 50% probability that the coalescent event is yet to occur, and a 50% probability that the final population will overshoot the observed current observed population size.

5 Likelihood surface given an observed $K(s)$

Naturally one would like to devise a more intuitive description of the likely location of the ancestral coalescent point of a random sample than the median estimate defined in the preceding section. Consider the following scenario. A BGW process is initiated with an initial scaled population $\kappa_{0}^{\rm true}$, and a random sample on $n$ individuals taken from the descendant population at some later time $s^{\rm true}$. The process is further conditioned on the event that the MRCA of the sample coincides with the initial time. This is in principle achieved by creating an arbitrarily large ensemble of processes, and rejecting those for which the required condition is not satisfied to within some small tolerance. Our aim is to determine a function ${\cal L}(s,\kappa_{0};\kappa)$ with the property that, for any subset $\Omega$ of the $(s,\kappa_{0})$-plane and observed final population $K(s^{\rm true})=\kappa$,

[TABLE]

We will refer to the function ${\cal L}(s,\kappa_{0};\kappa)$ as a “likelihood surface”. This surface can be considered as a tool for generating confidence regions in the sense that, by tuning $\Omega$ so that the integral in Eq. (43) has a value 0.95 for instance, one creates a $95\%$ confidence region for the point $(s^{\rm true},\kappa_{0}^{\rm true})$.

Subject to a conjecture which will be elucidated in the proof below, we have the following proposition:

Proposition 1**.**

The function

[TABLE]

where $K=\kappa$ is the observed current population, the function $u_{n}(s,\kappa_{0}|\kappa)$, defined in Theorem 1, and $v(s,\kappa_{0}|\kappa)$ is defined by Eq. (39) satisfies Eq. (43).

Proof.

Consider two neighbouring contours in the $(s,\kappa_{0})$-plane of the function $u_{n}$ given by Eq. (35) or (38), taking values $u$ and $u+du$ respectively, as shown in the left panel of Fig. 2. Recall that $u_{n}(s,\kappa_{0}|\kappa)$ is the probability that a BGW process initiated at any point on the first contour will result in a sample of size $n$ at time $s$ having a single ancestor in the initial population, conditional on $K(0)=\kappa_{0}$ and $K(s)=\kappa$. Similarly $u(s,\kappa_{0}|\kappa)+du(s,\kappa_{0}|\kappa)$ is the corresponding probability for the second contour. Then for fixed values of $\kappa_{0}$ and $\kappa$, the interval

[TABLE]

has a probability $du$ of containing the time $s^{\rm true}$ since the MRCA of the sample.

Furthermore, by Theorem 2, at a given value of $s$ the interval delimited by the contours $v$ and $v+dv$ shown in the left panel of Fig. 2, namely

[TABLE]

has a probability $dv$ of containing the initial population $\kappa_{0}^{\rm true}$, given an observed final population size $\kappa$.

These two intervals are random inasmuch as they depend on the final population size $K(s^{\rm true})|(K(0)=\kappa_{0}^{\rm true})$ of the conditioned BGW process described above. We conjecture that the two intervals are independent, and that therefore the probability that the intersection of the two corresponding regions in the $(s,\kappa_{0})$-plane contains the point $(s^{\rm true},\kappa_{0}^{\rm true})$ is $du\,dv$. Equivalently, we conjecture that if the $(s,\kappa_{0})$-plane is mapped onto the $(u,v)$-plane, as in Fig. 2, then the required function ${\cal L}(s,\kappa_{0};\kappa)$ maps to a uniform density over $[0,1]\times[0,1]$. The inverse mapping to the $(s,\kappa_{0})$-plane is effected by the Jacobian of the transformation, leading to Eq. (44). ∎

Contours of the likelihood surface ${\cal L}(s,\kappa_{0};\kappa)$ are marked in red in Figure 1 for $\kappa=3$ and $\kappa=2000$. Details of the calculation of the required derivatives are set out in A. These contour lines have the property that they bound a minimal area confidence region in the $(s,\kappa_{0})$-plane for a given integrated likelihood. Immediately noticeable is that in both cases the confidence regions are extremely skewed, and the median estimate defined at the end of the previous section is well separated from the maximum of the likelihood surface.

The marginal likelihood of the scaled time since coalescence, which we define from Eq. (44) as

[TABLE]

is plotted in Fig. 3 for a sample of $n=2$ individuals and various values of the currently observed scaled population $\kappa$. The numerical integration is tricky to evaluate for large values of $\kappa$ and small values of $s$ because of the difficulty in numerically evaluating $\partial v/\partial\kappa_{0}$ along the ridge in the joint density ${\cal L}(s,\kappa_{0};\kappa)$ running up to the point $(s,\kappa_{0})=(0,\kappa)$ (see Fig. 1(b), for instance). To produce the plots in Fig. 3(b) we avoided this problem by integrating by parts to obtain

[TABLE]

where the second partial derivative of $u_{n}$ is given in A, Eq. (63).

Since the marginal likelihood integrates to 1, a comparison can be made with the probability distribution of times since coalescence for a sample of size $n=2$ taken from an exponentially growing WF population, as determined by [21]. Their distribution is quoted in terms of a single parameter $\alpha=N_{0}r$ where $N_{0}$ is the final observed haploid population size and $r$ is the exponential growth rate. In terms of our notation, the expected number of offspring per individual per generation is $\lambda=e^{r}$, and the variance of the number of offspring per individual per generation for a WF model in the diffusion limit is $\sigma^{2}=1+O(1/M(t))$. Thus the condition $K(s)=2M(t)\log\lambda/\sigma^{2}=\kappa$ translates to $\kappa=2\alpha$, and Slatkin and Hudson’s Eq. (6) translates to

[TABLE]

This distribution is plotted as dashed curves in Fig. 3. Note that the parameters chosen in Fig. 3(b) match those of Fig. 4 in [21].

Figure 4 shows the marginal likelihood of the time since coalescence for samples of size $n=2,3,\ldots$. The choice of scaled current population, $\kappa=4500$, matches the naive model of human population growth during the upper Paleolithic period employed by Burden and Simon [5] to illustrate an application of BGW processes to mtE. In common with the predictions of the WF model for an exponentially growing population [21], we observe that MRCAs for samples of any size are likely to be restricted to a limited time range, consistent with a “star” shaped phylogeny.

6 Asymptotic behaviour for large $\kappa$

Plots of the marginal likelihood and of the Slatkin-Hudson density in Fig. 3 suggest that these functions each converge to an asymptotic shape as $\kappa\rightarrow\infty$. Specifically, by defining a shifted time scale

[TABLE]

one easily checks from Eq. (49) that

[TABLE]

as $\kappa\rightarrow\infty$. Thus the Slatkin-Hudson density rapidly takes on a shape related to the Gumbel distribution for large $\kappa$. As we next show, it turns out that each of the likelihood functions defined above for a BGW process also converges to an asymptotic shape depending only on $s^{\prime}$ (and, where relevant, $\kappa_{0}$) for large $\kappa$.

Theorem 3**.**

The functions $u_{n}$ and $v$ whose contours are plotted in Fig. 1, and hence the likelihood surface Eq,(44) and marginal likelihood Eq. (47), each take on a universal form as functions of $s^{\prime}=s-\log\kappa$ and $\kappa_{0}$ as $\kappa\rightarrow\infty$.

Proof.

From Eqs. (35), (38) and (27) we have that $u_{n}$ and $u_{\infty}$ are functions only of $w$ which behaves asymptotically as $(\kappa_{0}e^{-s^{\prime}})^{\frac{1}{2}}(1+O(\kappa^{-1}))$. It immediately follows that for large $\kappa$,

[TABLE]

where

[TABLE]

The asymptotic form of $v$ is obtained from Eqs. (39) and (25) by substituting Eq. (50) and discarding terms of $O(\kappa^{-1})$ to obtain

[TABLE]

where

[TABLE]

Here $F_{\rm Feller}(z,\kappa_{0})=\int_{0-}^{z}f_{\rm Feller}(\zeta,\kappa_{0})d\zeta$ is the cumulative distribution corresponding to the density defined by Eq. (24). ∎

Figure 5 shows contours of the functions $u_{2}$, $u_{\infty}$, $v$, and contours of associated likelihood surfaces, ${\cal L}(s^{\prime}+\log\kappa,\kappa_{0};\kappa)$ defined by Eq. (44), for relatively large values $\kappa=50$ and $1000$ of the observed current population. On the scale shown, the $\kappa=1000$ contours are indistinguishable from those corresponding to the universal functions $\tilde{u}_{n}$, $\tilde{u}_{\infty}$, and $\tilde{v}$ (not shown). Note that the universal functions are approached more rapidly for larger values of $n$. Contours of the limiting likelihood $\tilde{\cal L}(s^{\prime},\kappa_{0}):=\lim_{\kappa\rightarrow\infty}{\cal L}(s^{\prime}+\log\kappa,\kappa_{0};\kappa)$ and the limiting marginal likelihood $\tilde{\cal L}(s^{\prime}):=\lim_{\kappa\rightarrow\infty}{\cal L}(s^{\prime}+\log\kappa;\kappa)$, calculated from $\tilde{u}_{n}(s^{\prime},\kappa_{0})$ and $\tilde{v}(s^{\prime},\kappa_{0})$ are plotted in Fig. 6(c) to (f). We have observed that the numerically calculated limiting marginal likelihood for $n=\infty$, that is, the red curve in Figures 6(d) and (f), is a very close fit to a shifted Gumbel distribution. So far we have been unable to verify this analytically.

7 Numerical simulation

By exploiting the asymptotic universal functions arising from the $\kappa\rightarrow\infty$ limit we can compare the above theory with simulated data without needing to deal with explicit values of the final scaled population size $\kappa$.

Plotted in Fig. 6 are simulated data produced as follows. A set of trees was generated, each starting with a single ancestor and evolving as a BGW process, with each parent in the process independently producing a Poisson number of offspring with $\lambda=\sigma^{2}=1.01$. Trees were generated until a dataset of 40000 family trees, each surviving for $t_{\rm sim}=1000$ generations, was accumulated. For each tree $n=2$ individuals were randomly chosen from the final generation and their ancestry traced to locate the time of their MRCA, and to record the population size at that time. Each sample was taken from a separate tree to avoid correlations between the MRCAs of samples with a common history [4, 21]. For each family tree the time and population size corresponding to the MRCA of the entire population in the final generation was also recorded. These coalescence times and population sizes were then transformed to the scaled quantities $s^{\prime}$ and $\kappa_{0}$.

For consistency with the $\kappa\rightarrow\infty$ limit, a small number of trees with final scaled population sizes $\kappa<50$ were discarded to produce the scatter plots in Figures 6(a) and (b) and the histograms of $s^{\prime}$ values in Figures 6(c) and (d). It is clear that the right hand tails of the histograms fall short of the theoretical marginal likelihood curves. To appreciate the cause of this discrepancy, note that the marginal likelihoods in Figures 6(c) and (d) drop to almost zero outside a finite range, $s^{\prime}_{\rm lower}<s^{\prime}<s^{\prime}_{\rm upper}$. However, for practical reasons, a numerical simulation is restricted to be of finite time, namely $t_{\rm sim}$ generations, and because of this most of the trees in the sample have not existed sufficiently far in the past since the initial founder to cover this entire range. In an attempt to remedy this in Figures 6(e) and (f) we have also culled from the data all trees for which

[TABLE]

This second culling is more severe, and removes the majority of generated trees. Setting $s^{\prime}_{\rm upper}=5$, the dataset is reduced to trees satisfying $50<\kappa<141.2$, which reduces the original set to 192 trees. Nevertheless, the 2 final histograms are in closer agreement with the right hand tails of the theoretical marginal likelihoods.

8 Conclusions

We have addressed the problem of establishing confidence regions for the scaled time $s=t\log\lambda$ since coalescence to a MRCA of a random sample and the scaled population population size $\kappa_{0}=(2m_{0}\log\lambda)/\sigma^{2}$ at the time of the MRCA under the assumption that the current observed population is the result of a BGW process that has been evolving since some unknown time in the past. Here $t$ is the number of time steps since the MRCA and $m_{0}$ is the corresponding population size. The scalings are defined in terms of the mean $\lambda$ and variance $\sigma^{2}$ of the number of offspring per parent. The mean $\lambda$ is assumed to be close to and slightly larger than 1, and $t$ is assumed to be large, so that the diffusion limit can be applied. After scaling via Eq. (4) to implement the diffusion limit, the only free parameter in the problem of estimating $s$ and $\kappa_{0}$ is the currently observed scaled population $\kappa=2m\log\lambda/\sigma^{2}$, where $m$ is the currently observed physical population size.

The approach differs markedly from earlier analyses of the MRCA in a BGW diffusion [19, 14] in that the BGW process is not initiated from a pre-specified point in the past. Lacking information about the starting conditions is not a problem for tracing ancestry in WF-like models [21] which are inherently backward-looking in time. In such models Kingman’s coalescent is immediately applicable, and the concept of a probability density for the time since a MRCA is meaningful. For a BGW process, however, causality runs forward in time, and to be on firm ground we have couched our results in terms of confidence regions. Nevertheless some interpretation is in order.

To return to the question posed in the opening paragraph of this paper, suppose we are confronted with a population which we are told has evolved through a BGW process. Although the process is stochastic and its history unknown to us, one maintains an objective belief in the existence of a single realised path through the $(s,\kappa_{0})$-plane leading to the current observed scaled population $\kappa$, and at some point that unknown path will have passed through co-ordinates corresponding to the MRCA of the current population. Given that a BGW process is Markovian, the current population can be considered to be the result of a BGW process initiated from any point along the path, and in particular, from those co-ordinates corresponding to the MRCA. Our claim is that Eq. (44) calculated from the observed value of $\kappa$ defines a “likelihood surface” in the sense that any subset $\Omega$ of the $(s,\kappa_{0})$-plane has a probability $\int_{\Omega}{\cal L}(s,\kappa_{0};\kappa)dsd\kappa_{0}$ of containing the unknown but extant co-ordinates corresponding to the MRCA of the current population.

Having made this claim, we mention two final caveats. Firstly, note that Proposition 1 is contingent on a conjecture that the likelihood surface over the $(u,v)$-plane is a uniform joint density. Although the marginal densities in $u$ and $v$ are uniform, this is not necessarily the case for the joint density. Indeed any density of the form $1+\psi_{1}(u-v)+\psi_{2}(u+v)$ where $\psi_{1}$ and $\psi_{2}$ are periodic functions of a single variable with period 1 and which integrate to zero over $[0,1]$ will have uniform marginal densities. Without further mathematical proof, uniformity of the likelihood surface over the $(u,v)$-plane remains a conjecture.

Secondly, if the likelihood surface ${\cal L}(s,\kappa_{0};\kappa)$ is well defined, one feels it should be possible, at least in principle if not in practice, to devise a numerical simulation which will confirm its interpretation. The simulation in Section 7 attempts to do this, but as we have seen, there is some ambiguity about how to match this particular simulation with theory. Devising an appropriate simulation which is guaranteed to include the MRCA of the final population but does not impose any prior constraint on the population size at the time of the MRCA remains an open problem.

Acknowledgements

CB would like to thank Robert Griffiths, Robert Clark and Francis Hui for interesting and helpful discussions.

Appendix A Derivatives of $u_{n}(s,\kappa_{0}|\kappa)$ and $v(s,\kappa_{0}|\kappa)$

Derivatives of $u_{n}(s,\kappa_{0}|\kappa)$ and $v(s,\kappa_{0}|\kappa)$ are needed for the numerical determination of the likelihood surface defined by Eq. (44) and marginal likelihood calculation Eq. (48). Derivatives of $u_{n}(s,\kappa_{0}|\kappa)$ are straightforward to calculate from Eqs. (35) and (27). Using the identity [1, p376]

[TABLE]

we obtain

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

From Eq. (38) the analogous derivatives for $n=\infty$ are

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

Derivatives of $v(s,\kappa_{0}|\kappa)$ are less straightforward. The derivative with respect to $s$ can be calculated from the forward-Kolmogorov equation, Eq. (5), as follows;

[TABLE]

where Eqs. (2), (27) and (57) have been used in the last line. Evaluation of the derivative of $v$ with respect to $\kappa_{0}$ involves a numerical integration. Some straightforward but lengthy algebra gives

[TABLE]

Then from Eq. (39)

[TABLE]

where

[TABLE]