# Bayesian inference of species trees using diffusion models

Marnus Stoltz1 Boris Bauemer1 Remco Bouckaert2 Colin Fox3 Gordon Hiscott1
David Bryant1∗
$1.$ Department of Mathematics and Statistics, University of Otago
$2.$ Centre for Computational Evolution, University of Auckland
$3.$ Department of Physics, University of Otago.
$*$ Corresponding author<EMAIL_ADDRESS>
###### Abstract

We describe a new and computationally efficient Bayesian methodology for
inferring species trees and demographics from unlinked binary markers.
Likelihood calculations are carried out using diffusion models of allele
frequency dynamics combined with a new algorithm for numerically computing
likelihoods of quantitative traits. The diffusion approach allows for analysis
of datasets containing hundreds or thousands of individuals. The method, which
we call Snapper, has been implemented as part of the Beast2 package. We
introduce the models, the efficient algorithms, and report performance of
Snapper on simulated data sets and on SNP data from rattlesnakes and
freshwater turtles.

##  1 Introduction

Recent years have witnessed a proliferation in the number of methods for
inferring species trees from whole genomes (Bryant and Hahn, 2019). Some have
gone so far as to describe this as a paradigm shift in phylogenetics (Edwards,
2009). It is now widely accepted that (i) phylogenetic analysis needs to take
account of the varying evolutionary histories of different parts of the
genome, and (ii) estimation of evolutionary relationships between populations
(or species) should take account of evolutionary dynamics within populations
(or species).

In the systematics community, taking account of population dynamics in a
phylogenetic context has meant implementing some version of the multispecies
coalescent (Liu et al., 2008; Liu, 2008; Larget et al., 2010; Heled et al.,
2013). The coalescent models the distribution of gene trees within a
population; the multispecies coalescent is its natural extension to multiple
populations or species. One of the key advantages of the coalescent is that,
for neutral mutations, the gene tree can be decoupled from the mutation
process, a feature which forms the basis of many implementations of the model.
Nevertheless, the coalescent has its limitations. It is difficult to
incorporate selection, for example. Also, the running time of coalescent based
methods (e.g. BEAST, *BEAST, etc.) depends critically on the number of
individuals being sampled. In *BEAST a gene tree is sampled for each loci,
with the number of leaves in the gene tree given by the number of individuals
sampled from all populations.

Bryant et al. (2012) showed that explicit sampling of gene trees for each
locus could be avoided in the case of unlinked binary markers. This is
appropriate for the estimation of species trees from unlinked single
nucleotide polymorphisms (SNPs). Their algorithm was implemented in Snapp, and
can analyse data sets with hundreds of thousands of loci and up to 200
individuals (depending on computing resources and sampling difficulties).
Unfortunately Snapp does not scale that well as the number of individuals
increases. The running time of the algorithm algorithm is $O(n^{2}\log n)$
with $n$ individuals, which quickly becomes limiting. The far most common user
complaint about Snapp to developers has been running speed.

Rather than tinker with the Snapp algorithm, we have taken a completely new
approach with a different kind of model. Like Gutenkunst et al. (2010), Sirén
et al. (2010) and Lukic and Hey (2012) we use diffusion models, though we
apply them in a new way. The net result is that we can carry out Snapp-type
analyses but with essentially no limits on the number of individuals being
sampled.

Diffusion models, like the coalescent, are a convenient approximation of the
standard Wright-Fisher discrete models. In fact, diffusion models pre-date
coalescent models by a good fifty years (Wright, 1931; Kingman, 1982). Whereas
coalescent models describe the distribution of the genealogy or gene tree for
a sample of individuals, diffusion models describe the frequency of an allele
in the population as a whole. There are straightforward extensions
incorporating selection.

There are three challenges to overcome working with a model of allele
frequencies. First, we do not actually observe the population frequencies, we
just observe a sample drawn from that population. This problem is easily
solved by incorporating the sampling step explicitly into our likelihood.

Second, gene frequencies are continuous traits, so it is not possible to sum
over ancestral trait values as in Felsenstein’s pruning algorithm
(Felsenstein, 1981). For this, we follow a numerical approach developed in
Hiscott et al. (2016), though with some new twists to improve efficiency and
accuracy.

Third, diffusion models do not give explicit transition probabilities or
densities: these are only available via partial differential equations (PDEs).
We solve these differential equations numerically, extending a standard
spectral approach from a single population to the entire species tree. We note
that there is significant potential for extending our methods to other
diffusion models.

The outline of this paper is as follows: In Section 2 we briefly review the
population genetic diffusion models and some relevant theory, describe the
model on a species tree and discuss other related models. In Section 3 we
spend some time setting up the mathematics developed to compute the model
likelihood. Thereafter we present a numerical algorithm for computing the
likelihood and discuss the Bayesian software package Snapper. In Section 4 we
assess the perfomance of Snapper and compare running times to Snapp. In
Section 5 we showcase the capability of Snapper on a SNP dataset from
freshwater turtles. We conclude the paper with a summary of method
implementation and discuss some possible model extensions for future work.

##  2 Models and Methods

###  2.1 The Wright-Fisher diffusion

Our starting point is the Wright-Fisher discrete model of drift and mutation.
Suppose we have a population of $N$ diploid individuals, giving $2N$ copies of
each autosomal gene. We consider a gene with two alleles (say ‘red’ and
‘green’) and, to being with, a model with non-overlapping generations. With
random mating the number $X_{n+1}$ of red alleles in generation $n+1$ has a
binomial distribution with parameters $2N$ and
$(1-u)\frac{X_{n}}{2N}+v\left(1-\frac{X_{n}}{2N}\right)$, where $u$ is the
probability of mutating from red to green and $v$ the probability of mutating
from green to red. In this way the number of red alleles will follow a random
walk with discrete time (generation number) and states (from $0$ to $2N$).

The idea behind the diffusion models is to approximate this discrete random
walk by a continuous random walk that is easier to work with analytically. We
set up the approximation so that the larger $N$ gets, the better the
approximation fits. We describe the proportion of gene copies having the red
allele, rather than the total count. The state space of the process will be
the interval $[0,1]$. Instead of considering discrete generations, we
construct a random walk which is continuous in time.

Diffusion models for gene frequencies involve a rescaling of time. There are
simple, practical, reasons for this. As the population size $N$ gets larger,
the rate of genetic drift decreases, ultimately approaching zero drift in the
limit. However the effect of drift is something that we would like to model.
For this reason we change the units of time so that the rate of drift remains
approximately the same as $N$ increases, eventually converging to some non-
zero amount. The convention for diploid populations is to use a scale where
$1$ unit of time corresponds to $2N$ generations (one coalescent unit).

If we are to change the units of time, we need to adjust the rate of mutation
accordingly. Therefore the overall rate of change due to mutations from green
alleles to red alleles is $2Nu$. We adopt the standard notation and define
$\beta_{1}=2Nu$ and $\beta_{2}=2Nv$.

For each value of $t$ we let $f(x,t)$ denote the density of the allele
proportion at time $t$, noting that this allele proportion is a continuous
random variable with state space $[0,1]$. Surprisingly, there is very little
choice over what the function $f(x,t)$ might be, after a few basic assumptions
are made. See Etheridge (2011) for technical details and McKane and Waxman
(2007) for a discussion about how we need to incorporate point masses at $0$
and $1$ into $f(x,t)$.

As is often the case in mathematical modelling, we work with the function $f$
indirectly using a PDE. Stochastic process theory (Øksendal, 2003, chap. 5)
tells us that the function $f$ satisfies the PDE

$\displaystyle\frac{df(x,t)}{dt}$
$\displaystyle=-\frac{d}{dx}\left(\beta_{1}(1-x)-\beta_{2}x\right)f(x,t)+\frac{1}{2}\frac{d^{2}}{dx^{2}}x(1-x)f(x,t).$
(1)

This equation by itself is not enough to uniquely determine $f$. We also need
to specify what $f$ looks like at the boundaries. To specify that the
distribution of the initial state is given by some density $\pi$ we add the
initial condition

$f(x,0)=\pi(x)\mbox{ for $x\in[0,1]$ }.$ (2)

Note that to specify the initial value exactly, we need $\pi$ to be a Dirac-
delta function $\pi(x)=\delta(x-x_{0})$ which is essentially a infinitely thin
spike on one value $x_{0}$. Even with these initial conditions, the PDE (1)
does not uniquely determine the function $f$. We also need to add boundary
conditions at $x=0$ and $x=1$ to guarantee that the probability of going
outside the interval $[0,1]$ is always zero. Writing this as a condition on
integrals of $f$ and then plugging in the PDE, eventually leads to what in
physics is known as a zero-flux condition

$-(\beta_{1}(1-x)-\beta_{2}x)f(t,x)+\frac{1}{2}\frac{d}{dx}x(1-x)f(t,x)=0,$
(3)

when $x=0$ or $x=1$. See McKane and Waxman (2007) for details.

The PDE (1) together with the boundary conditions (2) and (3) are just a
mathematically convenient, if not particularly transparent, way to describe
the function $f$. The function $f$, in turn, is just a way to approximate the
probabilities in the original discrete process. The diffusion approximation
works well only if mutations rates are assumed small to begin with. This will
give a reasonably small error in the diffusion approximation. As shown in
Ethier and Norman (1977) we still have convergence in expectation with an
error bound on the diffusion approximation behaving like $O(u+v+1/N)$. Note
that if $u$ or $v$ are large (compared to, say, $\frac{1}{N}$) then the
diffusion approximation will fail miserably. In a recent simulation study,
Tataru et al. (2017) used simulations to quantify the error from the diffusion
approximation in small populations, and found that diffusion models gave
reasonable approximations even when the population has fewer than 100
individuals.

###  2.2 Modelling allele frequencies on a species tree

The diffusion model describes how allele frequencies change over time in a
single population. The model extends directly to multiple populations in a
species tree (Sirén et al., 2010). As in Bryant et al. (2012) we think of the
root of the species tree to be at the top of the tree and the leaves at the
bottom. The model describes evolution of the allele frequencies from the top
of the tree to the bottom.

The allele frequency at the root has a distribution given by the stationary
density of the diffusion model (in this case, a beta distribution). The allele
frequency at the bottom of a branch has a distribution given by the diffusion
model with an initial value equal to the allele frequency at the top of the
branch. At a speciation, the two daughter populations have the same allele
frequency as the parent population.

We now formalise these ideas. Suppose that the branches in the species tree
are numbered $i=1,2,\ldots,2n-2$ where, for convenience, the branches adjacent
to the leaves are numbered $1,2,\ldots,n$. Let $\mathcal{X}^{T}_{i}$ denote
the allele frequency immediately below the top of branch $i$. Let $X^{B}_{i}$
denote the allele frequency immediately above the bottom of branch $i$. Let
$X_{\rho}=X_{\rho}^{B}$ denote the allele frequency at the root.

At the root, $\rho$, the proportion of red alleles in the ancestral population
equals the stationary distribution of the diffusion model. The stationary
distribution of the diffusion equation (1) has a beta distribution (Ewens,
2004) with density

$\pi_{\rho}(x|\beta_{1},\beta_{2})=\frac{\Gamma\\{2\beta_{1}+2\beta_{2}\\}}{\Gamma\\{2\beta_{1}\\}\Gamma\\{2\beta_{2}\\}}x^{2\beta_{1}-1}(1-x)^{2\beta_{2}-1},\qquad
0<x<1$ (4)

Along any branch in the species tree, the changes in allele frequencies are
modelled using the diffusion process. The allele frequency at the start of the
branch gives the initial density, $f(x,0)$, for $x\in[0,1]$. The distribution
of the allele frequency $y$ at the end of the branch is then given by
$f(y,t)$, where $t$ is the length of the branch in units of $2N$ generations
and $y\in[0,1]$.

At a speciation, we make the assumption that there is no correlation between
allele type and speciation. Hence the two descendant species are assumed to
have identical allele frequencies to the parent node.

We have, therefore, a model for allele frequencies over the whole tree.
However allele frequencies at the leaves are not directly observed. Instead we
have a sample of $n_{i}$ individuals, giving $2n_{i}$ allele copies taken from
each species $i$. If $x_{i}$ is the red allele frequency for the population
then the observed number $r_{i}$ of red alleles in the the sample for this
gene has binomial distribution with parameters $2n_{i}$ and $x_{i}$, so that

$P[R_{i}=r|x_{i}]={2n_{i}\choose r}x_{i}^{r}(1-x_{i})^{2n_{i}-r},$

see Sirén et al. (2010).

In Algorithm 1 we describe how to simulate allele frequencies under this
model. Note that in a pre-order traversal we visit every node in order so that
the children of a node are always visited after the node itself.

Algorithm 1 Simulate allele frequencies on a tree

1:procedure TreeSim(S)

2: Generate $X^{B}_{\rho}$ from density $\pi(x|\beta_{1},\beta_{2})$ given in
(4)

3: for all nodes $i$ except the root in a pre-order traversal do

4: $X_{i}^{T}\leftarrow X_{p(i)}^{B}$, where $p(i)$ denotes the parent of $i$

5: Use diffusion to simulate $X^{B}_{i}$ given initial value $X^{T}_{i}$

6: if node $i$ is a leaf then

7: Generate $r_{i}$ from $Bin(2n_{i},X_{i}^{B})$

8: end if

9: end for

10:end procedure

We use the algorithm of Jenkins et al. (2017) to simulate diffusions along
each branch. Alternative methods for simulating diffusions include truncating
a spectral expansion of the transition function (Song and Steinrücken, 2012;
Steinrücken et al., 2014) or numerical solutions of the Kolmogorov equations
(Williamson et al., 2005; Bollback et al., 2008; Gutenkunst et al., 2009).
Care needs to be taken when implementing these methods since approximation
errors close to the boundary can become large, a problem that Jenkins et al.
(2017) explicitly avoid.

###  2.3 Analytical formula for partial likelihoods

We now describe how to describe the likelihood functions analytically. Overall
our approach is to use dynamic programming, as in Felsenstein’s pruning
algorithm (Felsenstein, 1981). For each node $i$ and each ancestral state $x$,
which is our case is an allele frequency, we define two partial likelihoods,
$\ell^{B}_{i}(x)$ and $\ell^{T}_{i}(x)$. The first is the probability of
observing all allele counts for taxa in the subtree rooted at $i$, conditional
on the state $X^{B}_{i}$ at the bottom of branch $i$ being equal to $x$. The
second is defined in the same way, but is instead conditioned on the state
$X^{T}_{i}$ at the top of the branch. A similar approach was used in Bryant et
al. (2012).

We now show how to describe these partial likelihoods recursively for an
individual site $m$.

####  2.3.1 Partial likelihoods at the leaves

Suppose that node $i$ is a leaf and that we have sampled $n_{i}$ diploid
individuals from this population. Let $r_{i}$ denote the number of observed
gene copies carrying the red allele for site $m$. If $x$ is the proportion of
red alleles in the population then $r_{i}$ has a binomial distribution with
parameters $n_{i}$ and $x$. Hence the partial likelihood at a leaf $i$ is
given by

$\ell_{i}^{B}(x)={2n_{i}\choose r_{i}}x^{r_{i}}(1-x)^{2n_{i}-r_{i}}.$ (5)

####  2.3.2 Partial likelihoods at a speciation

Let $j$ and $k$ be the children of node $i$. We assumed that the allele
frequencies for daughter populations after a speciation are exactly those of
the parent population before speciation. The partial likelihoods at the bottom
of the branch above $i$ is then the product of partial likelihoods at the tops
of the branches immediately below $i$, so

$\ell_{i}^{B}(x)=\ell_{j}^{T}(x)\ell_{k}^{T}(x)$ (6)

for all $x\in[0,1]$.

####  2.3.3 Partial likelihoods along a branch

Let $N_{i}$ be the effective population size for the branch directly above
node $i$, and let $\tau_{i}$ denote the length of the branch above node $i$,
measured in numbers of generations. Then $t_{i}=\frac{\tau_{i}}{2N_{i}}$
equals the length of the branch in coalescent units. To derive the partial
likelihood $\ell_{i}^{T}(x)$ at the top of this branch, we integrate over all
frequency values at the bottom of the branch.

$\ell_{i}^{T}(x)=\int_{0}^{1}f(y,t_{i})\ell_{i}^{B}(y)dy.$ (7)

In this equation, $f(y,t_{i})$ is the density for the allele frequency $y$ the
bottom of the branch conditional on an initial value $x$ at the top of the
branch, while $\ell_{i}^{B}(y)$ is the partial likelihood conditional on
allele frequency $y$ at the bottom of the branch.

To get a handle on (7) we define

$g(x,t)=\int_{0}^{1}f(y,t)\ell_{i}^{B}(y)dy,\text{ for $0<t<t_{i}$}.$ (8)

Then $g$ satisfies the PDE,

$\frac{dg(x,t)}{dt}=(\beta_{1}(1-x)-\beta_{2}x)\frac{d}{dx}g(x,t)+\frac{1}{2}x(1-x)\frac{d^{2}}{dx^{2}}g(x,t).$
(9)

with initial condition

$g(x,0)=\ell_{i}^{B}(x)\text{, for }x\in[0,1].$ (10)

Solutions to the PDE (9) are unique and well-behaved without specifying
boundary conditions for $t>0$ (Epstein and Mazzeo, 2010).

####  2.3.4 Likelihood at the root of a tree

Equations (5), (6) and (7) can be applied to the entire tree, starting with
the leaves and moving upwards towards the root. They define the root
likelihood $\ell_{\rho}^{B}(x)$, which is the probability of observing all
allele counts for all populations, conditioning on the allele frequency
$X_{\rho}=X_{\rho}^{B}$ being equal to $x$.

To define the likelihood of the tree we now integrate over the allele
frequency $x$ at the root, multiplied by the stationary density of $x$ given
in (4):

$L_{m}=\int_{0}^{1}\pi(x|\beta_{1},\beta_{2})\ell_{\rho}^{B}(x)dx.$ (11)

This is the probability for a single site $m$. The likelihood from all markers
is found by multiplying the probabilities from each site together (or, in
practice, summing the log-likelihoods)

$\log L=\sum_{m}\log(L_{m}).$ (12)

In the case of SNP data we need to compute the conditional likelihood,
$L_{m}/(1-L_{0})$, where $L_{0}$ is the probability of observing a constant
site (See Bryant et al. (2012) and Felsenstein (1992)).

###  2.4 Translation of parameters

One of the more confusing aspects of working across population genetics and
phylogenetics is the different ways that different communities parameterize
different variables. In this section we summarise the parameters and outputs
for our model, and show how they might be converted into other formulations.

There are four groups of input variables in our model. These are the variables
that could conceivably be inferred from data using the likelihood function.
They are

  1. 1.

The species tree.

  2. 2.

The branch lengths in the species tree. In our model, we measure the length of
a branch in number of generations. Let $\tau_{i}$ denote the length of the
branch above node $i$.

  3. 3.

The mutation rates $u$ and $v$ giving the probabilities of mutating from a red
to a green allele or a green to a red allele, each generation.

  4. 4.

The effective population size $N_{i}$ for each branch $i$ (that is, the branch
above node $i$).

Different methods for inferring species trees describe these parameters in
different ways. The convention in phylogenetics is to express branch lengths
in terms of expected substitutions per site. The use of the term
’substitution’ is confusing, as the whole premise behind the multispecies
coalescent is to model incomplete substitutions. However for neutral
mutations, the rate of substitutions equals the rate of mutations. Under our
model, the rate of mutations per generation is

$\displaystyle\mu$ $\displaystyle=P(\mbox{ lineage has red allele })u+P(\mbox{
lineage has green allele })v$ (13)
$\displaystyle=\frac{v}{u+v}u+\frac{u}{u+v}v,$ (14)

so a branch of length $\tau_{i}$ generations corresponds to a branch length of
$\mu\tau_{i}$ expected substitutions (mutations) per site.

Methods which ignore the branch lengths in gene trees are not able to infer
both population sizes and branch lengths, and instead use a single parameter
per branch, typically measured in coalescent units (e.g. (Vachaspati and
Warnow, 2015; Liu et al., 2010)). If $\tau_{i}$ is the length of a branch in
numbers of generations, the length of the branch in coalescent units is given
by $t_{i}=\frac{\tau_{i}}{2N_{i}}$.

In our model, we parameterize effective population size for branch $i$
directly as $N_{i}$. Snapp (Bryant et al., 2012), Best (Liu, 2008) and BPP
(Yang, 2015) instead use $\theta_{i}=4N_{i}\mu$ as the parameter for effective
population size. Under an infinite sites model, $\theta_{i}$ equals the
expected proportion of sequence differences between two individuals from the
same population. In a finite sites model, $4N_{i}\mu$ is an overestimate of
this expectation, due to backward mutations.

###  2.5 Dealing with confounded rates and rate matrices

An important issue with the diffusion model as we have described it, and
indeed with many multispecies coalescent models, is that there is an
identifiability problem with rates. If the mutation rates are multiplied by
some constant $c$ and at the same time branch lengths and population sizes are
multiplied by $\frac{1}{c}$, the probability of the data remains the same.

One solution is to estimate the mutation rate $\mu$ beforehand and use this
value, or a prior distribution around that value, to include $\mu$ as part of
our model. The prior distribution for $u$ and $v$ is then reformulated so that
(14) is satisfied. This strategy is used by StarBeast, where the average
substitution (mutation) rate $r_{\mu}$ is fixed ahead of time.

An alternative strategy is to express branch lengths in terms of expected
substitutions (mutations) per site and population sizes in terms of the
$\theta$ parameter, as both of these are invariant to a choice of $\mu$,
(Yang, 2015; Bryant et al., 2012). As a consequence, the effective population
sizes $N_{i}$ cannot be inferred without additional information. This approach
adapts well to phylogenetics where it is customary to describe mutation rates
using a normalised rate matrix, such as the Jukes-Cantor and HKY85 matrices

$\mathbf{Q}=\left[\begin{matrix}-&1/3&1/3&1/3\\\ 1/3&-&1/3&1/3\\\
1/3&1/3&-&1/3\\\ 1/3&1/3&1/3&-\end{matrix}\right]\mbox{ and
}\mathbf{Q}=\frac{1}{r}\left[\begin{matrix}-&\pi_{C}&\kappa\pi_{G}&\pi_{T}\\\
\pi_{A}&-&\pi_{G}&\kappa\pi_{T}\\\ \kappa\pi_{A}&\pi_{C}&-&\pi_{T}\\\
\pi_{A}&\kappa\pi_{C}&\pi_{G}&-\end{matrix}\right].$

In this matrix, the stationary probabilities for the nucleotides are
$\pi_{A},\pi_{C},\pi_{G},\pi_{T}$, the parameter $\kappa$ controls the ratio
of transitions to transversions, and the diagonal elements are chosen to make
rows sum to zero. The Jukes-Cantor matrix on the left is already normalised,
while in the HKY85 model $r$ would be chosen to make the overall substitution
rate equal to $1$.

These matrices can be incorporated directly into Snapper, either using prior
information for $\mu$, or when using the same branch length and population
parameter scheme as BPP. For any site, if $X$ and $Y$ are the two nucleotides
present and we (arbitrarily) assign red to $X$ and green to $Y$ then the
appropriate choices for $u$ and $v$ are simply $\mathbf{Q}_{XY}$ and
$\mathbf{Q}_{YX}$. For example, under the HKY85 model, if the red allele is
nucleotide A and green allele is nucleotide C then the corresponding rates are
$u=\mu\mathbf{Q}_{AC}=\pi_{C}$ and $v=\mathbf{Q}_{CA}\pi_{A}$. We are
approximating a four state model with a two state model, but the error
introduced should be minimal in the absence of divergences or highly variable
sites.

###  2.6 Comparison with allele frequency spectrum methods

We note that diffusion models are already used in population genetics to
approximate changes in the allele frequency spectrum (AFS) (Gutenkunst et al.,
2009; Lukic and Hey, 2012; Racimo et al., 2016). For each value of $x$ between
$0$ and $1$, the AFS gives the proportion of sites for which the derived
allele appears in the proportion $x$ of the sample. For example, we might
consider all sites where the the derived allele appears in 30 of the 100
sampled genomes. If $5\%$ of all sites fall into this category then the AFS
value for 0.3 will be $5\%$.

Methods based on the AFS to infer genetic parameters for a single population
are widespread in the literature (Boitard et al., 2016; Lapierre et al., 2017;
Nielsen, 2000). In all of these analyses, it is assumed that the population
dynamics have achieved stationarity when the AFS predicted by the model is
used.

When there are multiple populations divergence times become important. Suppose
that there are $s$ species. For every vector $(x_{1},x_{2},\ldots,x_{s})$ of
numbers between $0$ and $1$, the joint-AFS gives the proportion of sites where
the derived allele appears in a proportion $x_{1}$ of the individuals in the
first population, $x_{2}$ of the individuals in the second population, and so
on. While the AFS for each population by itself will be the stationary AFS,
the joint-AFS for multiple populations depend on the time since separation.
Fortunately the PDE determining the predicted joint-AFS is a straight-forward
extension of the PDE for a single population (Gutenkunst et al., 2009; Lukic
and Hey, 2012; Racimo et al., 2016).

There are two main advantages of an approach based on the joint-AFS when
compared to the method we describe here. First, the PDE for the joint-AFS only
needs to be solved once for all sites, whereas in our approach we end up
solving the PDEs once for each (distinct) site. Second, migration between
populations can be incorporated quite simply into AFS calculation, whereas it
would break down the dynamical programming strategy that we use.

The main disadvantage of the joint-AFS strategy is that the PDE has as many
dimensions as the number of species, and suffers from the curse of
dimensionality (Gutenkunst et al., 2009), leading to an exponential growth in
the size of grid used by standard numerical methods. This factor severely
limits the number of species which can be considered concurrently, whereas in
our approach the algorithm scales linearly with the number of species.

##  3 Computing the likelihood efficiently

In the previous section, we showed how the recursions for the partial
likelihoods can be derived analytically. We did not show how to actually
compute those likelihoods. Indeed, computing the likelihoods amounts to an
extremely high dimensional integration problem. Our approach is to use
numerical techniques combined with dynamic programming. While the likelihoods
we compute are approximate, we can still bound the error.

###  3.1 Shifted Chebyshev polynomials

Hiscott et al. (2016) describe a general strategy for computing likelihoods
numerically whereby partial likelihoods are evaluated on a mesh of values at
each node, and accurate quadrature methods used to carry out the actual
computation. We will extend the same strategy by using a set of basis
functions to express approximate partial likelihoods instead of a mesh of
values. The basis functions are cleverly chosen to help solve equation (9)
efficiently and accurately.

The basis functions we use are called the shifted Chebyshev polynomials of the
first kind, denoted

$T_{0}^{*}(x),T_{1}^{*}(x),T^{*}_{2}(x),\ldots.$

Shifted Chebyshev polynomials are defined on $[0,1]$ and have particularly
nice properties (Fox and Parker, 1968). They also have a lot of equivalent
definitions, but the simplest uses the following recursion,

$\displaystyle T^{*}_{0}(x)$ $\displaystyle=1$ $\displaystyle T^{*}_{1}(x)$
$\displaystyle=2x-1$ $\displaystyle T^{*}_{k}(x)$
$\displaystyle=2(2x-1)T^{*}_{k-1}(x)-T^{*}_{k-2}(x).$

The shifted Chebyshev polynomials are related to the better-known Chebyshev
polynomials $T_{0},T_{1},T_{2},\ldots$ by the identity:
$T^{*}_{k}(x)=T_{k}(2x-1)$. That is, they are obtained by shifting and scaling
the Chebyshev polynomials to have domain $[0,1]$, see Mason and Handscomb
(2002).

There are two main ways of using shifted Chebyshev polynomials to approximate
$\ell_{i}^{B}$ and $\ell_{i}^{T}$ as functions of $x$. The first is to
approximate the function directly as a linear combination of the shifted
Chebyshev polynomials, that is, finding sets of coefficients
$\lambda^{B}_{i,0},\lambda^{B}_{i,1},\lambda^{B}_{i,2},\ldots$ and
$\lambda^{T}_{i,0},\lambda^{T}_{i,1},\lambda^{T}_{i,2},\ldots$ so that for all
$x$ in $[0,1]$ we have

$\ell_{i}^{B}(x)\approx\sum_{k=0}^{K}\lambda_{i,k}^{B}T_{k}^{*}(x)\quad\text{and}\quad\ell_{i}^{T}(x)\approx\sum_{k=0}^{K}\lambda_{i,k}^{T}T_{k}^{*}(x).$

It can be shown that the error in this approximation drops exponentially
quickly as the number $K$ of basis functions increases. We therefore only need
to store a few coefficients in order to evaluate the partial likelihood
function at any $x$, with small error.

The second way of obtaining an approximation is by determining the values
$\ell_{i}^{B}(x)$ and $\ell_{i}^{T}(x)$ at a pre-specified set of points
$x_{0},x_{1},\ldots,x_{K}$ and then finding the unique combination of
coefficients $\lambda_{i,0}^{B},\dots,\lambda_{i,K}^{B}$ such that

$\sum_{k=0}^{K}\lambda_{i,k}^{B}(x_{j})=\ell_{i}^{B}(x_{j})$

for all $j=0,1,\ldots,K$. This is called polynomial interpolation. As it
happens, there is a particular choice of points $x_{0},\ldots,x_{K}$ for which
we can switch back and forth between function values

$\ell_{i}^{B}(x_{0}),\dots,\ell_{i}^{B}(x_{K})$

and interpolation coefficients

$\lambda_{i,0}^{B},\dots,\lambda_{i,K}^{B}$

and back with little numerical error and in $O(K\log K)$ time (Waldvogel,
2006; Trefethen, 2013). These are called the Chebyshev-Lobatto points, defined
by the rather opaque formula

$x_{j}=\left(1-\cos\left(\frac{2\pi j}{K}\right)\right)/2,\text{ for
$j=0,\dots,K$}.$

These points are all in the interval $[0,1]$ with a denser packing of points
nearer $0$ and $1$. They equal the $x$-coordinates of points spaced regularly
around a semi-circle.

We use both coefficient values $\lambda_{i,k}^{B},\lambda_{i,k}^{T}$ and
function values $\ell_{i}^{B}(x_{j}),\ell_{i}^{T}(x_{j})$ when approximating
the partial likelihoods $\ell_{i}^{B}$ and $\ell_{i}^{T}$.

###  3.2 Approximate partial likelihoods at a leaf

We compute our approximate likelihood at the bottom of a leaf branch by simply
evaluating the function

$\ell_{i}^{B}(x_{j})={2n_{i}\choose
r_{i}}x_{j}^{r_{i}}(1-x_{j})^{2n_{i}-r_{i}}\text{ for j = 1, \ldots, K}$ (15)

at the Gauss-Lobatto points. We then compute corresponding coefficients
$\lambda_{i,0}^{B},\lambda_{i,1}^{B},\ldots,\lambda_{i,K}^{B}$ using the FFT
algorithm (Trefethen, 2013).

###  3.3 Solving the diffusion model numerically

Shifted Chebyshev polynomials provide the foundation for the numerical methods
we use to solve the PDE (9). Note that, unlike the forward diffusion (1),
solutions to this backward diffusion are nice smooth functions, meaning that
we can avoid some of the numerical headaches encountered by those using the
forward diffusion directly (Lukic and Hey, 2012).

We approximate $g(x,t)$ in terms of shifted Chebyshev polynomials

$g(x,t)\approx\sum_{k=0}^{K}\lambda_{k}(t)T^{*}_{k}(x).$ (16)

Then

$\displaystyle\frac{\partial g(x,t)}{\partial t}$
$\displaystyle\approx\sum_{k=0}^{K}\frac{d\lambda_{k}(t)}{dt}T^{*}_{k}(x),$
(17) $\displaystyle\frac{\partial g(x,t)}{\partial x}$
$\displaystyle\approx\sum_{k=0}^{K}\lambda_{k}(t)\frac{dT^{*}_{k}(x)}{dx}$
(18) and $\displaystyle\frac{\partial^{2}g(x,t)}{\partial x^{2}}$
$\displaystyle\approx\sum_{k=0}^{K}\lambda_{k}(t)\frac{d^{2}T^{*}_{k}(x)}{dx^{2}}.$
(19)

Let ${\bm{\lambda}}(t)$ denote the vector of values
$\lambda_{0}(t),\ldots,\lambda_{K}(t)$. Formulas for the derivative of shifted
Chebyshev polynomials lead eventually to a $(K+1)\times(K+1)$ matrix $\bm{D}$
such that

$\sum_{k=0}^{K}\lambda_{k}(t)\frac{dT^{*}_{k}(x)}{dx}=\sum_{k=0}^{K}(\bm{D}\bm{\lambda}(t))_{k}T^{*}_{k}(x),$

and

$\sum_{k=0}^{K}\lambda_{k}(t)\frac{d^{2}T^{*}_{k}(x)}{dx^{2}}=\sum_{k=0}^{K}(\bm{D}^{2}\bm{\lambda}(t))_{k}T^{*}_{k}(x).$

Indeed after some tedious algebra we can derive a $(K+1)\times(K+1)$ $\bm{Q}$
matrix so that

$(\beta_{1}(1-x)-\beta_{2}x)\frac{d}{dx}\sum_{k=0}^{K}\lambda_{k}(t)T^{*}_{k}(x)+\frac{1}{2}x(1-x)\frac{d^{2}}{dx^{2}}\sum_{k=0}^{K}\lambda_{k}(t)T^{*}_{k}(x)=\sum_{k=0}^{K}(\bm{Q}\bm{\lambda}(t))_{k}T^{*}_{k}(x)$

for all $x\in[0,1]$. Plugging this and (17) into (9) we obtain an equation

$\sum_{k=0}^{K}\frac{d\lambda_{k}(t)}{dt}T^{*}_{k}(x)=\sum_{k=0}^{K}(\bm{Q}\bm{\lambda}(t))_{k}T^{*}_{k}(x)$

for an approximate solution to the PDE, giving the system

$\frac{d{\bm{\lambda}}(t)}{dt}={\bm{Q}}\cdot{\bm{\lambda}}(t).$ (20)

with initial condition ${\bm{\lambda}}(0)$ given as the Chebyshev coefficients
of $g(x,0)$.

Like Bryant et al. (2012) we use rational approximations to $\exp({\bf
Q}t)\bm{\lambda}(0)$ on the negative real axis computed using a Caratheodory
Fejer-approximation. Additionally, we use techniques adapted from Fox and
Parker (1968) and outlined in the appendix (A.1) to take advantage of
structure in the matrix $\mathbf{Q}$, allowing an implementation of the
Caratheodory-Fejer approximation which runs in $\mathcal{O}(K)$ time.

To summarise, consider a node $i$ and let $t_{i}$ denote the length (in
coalescent units) of the branch connecting $i$ to its parent. Suppose
$\ell_{i}^{B}$ with corresponding coefficients $\lambda^{B}_{i}$ is already
computed at the bottom of the branch. We then compute partial likelihood
$\ell_{i}^{T}$ at the top of the branch by

  1. 1.

Setting ${\bm{\lambda}}(0)=\lambda^{B}_{i,0},\dots,\lambda^{B}_{i,K}$

  2. 2.

Computing a numerical approximation for
$\exp({\bm{Q}}t_{v}){\bm{\lambda}}(0)$.

  3. 3.

$\ell_{i}^{T}$ = $\lambda(t_{0})_{0},\dots,\lambda(t_{0})_{K}$

###  3.4 Approximate partial likelihoods at a speciation

Consider a node $i$ with two children $j$ and $k$. We suppose that the
approximations for the partial likelihoods
$\ell_{j}^{T}(x_{0}),\ell_{j}^{T}(x_{1}),\ldots\ell_{j}^{T}(x_{K})$ and
$\ell_{k}^{T}(x_{0}),\ell_{k}^{T}(x_{1}),\ldots\ell_{k}^{T}(x_{K})$. We then
have

$\displaystyle\ell_{i}^{B}(x_{p})=\ell_{j}^{T}(x_{p})\ell_{k}^{T}(x_{p})$ for
all $p=0,1,\ldots,K$. (21)

This computation takes $O(K)$ time.

###  3.5 Approximate likelihoods at the root

The final integration to carry out is over the partial likelihood at the root

$\int_{0}^{1}\pi(x|\beta_{1},\beta_{2})\ell_{\rho}^{B}(x)dx,$

see (11). The density for a beta distribution (4) can be infinite at the
boundaries, meaning that numerical integration techniques such as Clenshaw-
Curtis quadrature can give poor approximations. The solution is to separate
out those parts which are difficult to integrate numerically and determine
them analytically. Suppose we have an approximation of $\ell_{\rho}^{B}(x)$ as
a polynomial

$f(x)=\sum_{k=0}^{K}\lambda_{\rho,k}^{B}T^{*}_{k}(x)\approx\ell_{\rho}^{B}(x).$

We factor this polynomial as

$f(x)=x(1-x)g(x)+f(0)+(f(1)-f(0))x$

and then compute

$\int_{0}^{1}f(x)\pi(x|\beta_{1},\beta_{2})dx=\int_{0}^{1}g(x)x(1-x)\pi(x|\beta_{1},\beta_{2})dx+\int_{0}^{1}(f(0)+(f(1)-f(0))x)\pi(x|\beta_{1},\beta_{2})dx.$

Noting that the first integral is now well-behaved while the second integral
evaluates to $f(0)+(f(1)-f(0))\frac{\beta_{1}}{\beta_{1}+\beta_{2}}$ by
properties of the Beta distribution.

###  3.6 Computing the log-likelihood of a species trees

Algorithm 2 summarises the numerical calculation of the likelihood. This
algorithm takes $\mathcal{O}(sKlog(K))$ time per site where $\mathcal{O}(ns)$
pre-processing, where $s$ is the number of species, $K$ the number of
Chebyshev basis functions and $n$ is the number of individuals at a site. The
pre-processing step involves counting the frequencies of allele types in each
population. In practice this step could be carried out once per data set,
rather than once per tree evaluated.

Algorithm 2 Snapper: Computes tree log-likelihoods

1:procedure log-likelihood($X$, S)

2: for All sites $s$ do

3: for All nodes $i$ in a post-order traversal of the tree do

4: if $i$ is a leaf node then

5: Evaluate $\ell_{i}^{B}(x_{j})$ for each $j$ using (15)

6: else

7: Evaluate $\ell_{i}^{B}(x_{j})$ for each $j$ using (21)

8: end if

9: Compute coefficients $\lambda_{i,j}^{B}$ from values $\ell_{i}^{B}(x_{j})$

10: Compute coefficients $\lambda_{i,j}^{T}$ by solving the PDE (20) with
${\bm{\lambda}}(0)={\bm{\lambda}}_{i}^{B}$

11: Compute values $\ell_{i}^{T}(x_{j})$ from coefficients $\lambda_{i,j}^{T}$

12: end for

13: Compute $L_{s}$ using numerical integration.

14: end for

15: return $\sum_{s}log(L_{s})$

16:end procedure

The conversion to and from coefficients to function values in the Chebyshev
expansion in lines 7 and 9 each takes $O(K\log K)$ time, using the FFT
algorithm reviewed above. Solution of the PDE in line 10 takes $O(K)$ time.

###  3.7 Snapper

The likelihood algorithm forms the core of a Bayesian inference software
package, Snapper. The software is open-source and available to download at
https://github.com/rbouckaert/snapper. Like Snapp, it takes biallelic data at
multiple loci for multiple individuals, in a set of species and returns
samples from the joint posterior distribution of (i) species phylogenies, (ii)
species divergence times and (iii) effective population sizes. As in Snapp we
implemented a dynamic cache-based system to store partial likelihoods on
different subtrees and multithreading to take advantage of parallel
computation on multiple core machines or graphics processing units. The range
of prior distributions, and flexible prior specification, remains essentially
unchanged.

The MCMC proposal function implemented is almost the same subset of BEAST2
(Bouckaert et al., 2019) move proposals implemented for Snapp. We add one new
proposal selects a subtree of the species tree and scales all population sizes
within that subtree simultaneously, a refinement of an existing proposal.

The user can specify number $K$ of basis functions. Some recommendations are
made in the Snapper manual regarding the number of basis functions to use
given the size of a dataset. Also contained in the software package are
simulators and python scripts to integrate Snapper with the iPyRad data
pipeline to assist with streamlining data analysis.

##  4 Performance assessment

We have tested the likelihood algorithm and the Snapper software extensively.
Here we highlight some results related to the performance of Snapper. Details
of additional simulations, including tests checking that posterior
distributions recovered model parameters, can be found in the appendix (A.2).

We conducted the following experiments: (i) An assessment of error convergence
as the number $K$ of Chebyshev basis functions increase (ii) A comparison of
inferred parameters between Snapp and Snapper given a dataset (iii) An
assessment of how computational time estimates for Snapp and Snapper scale in
terms of size of the data set. All simulations and inference were run on a
laptop with an Intel i7-8565U CPU.

###  4.1 Numerical error convergence

We assessed the rate of error convergence for different tree heights, tree
sizes and tree topologies. In particular we wanted to verify that we have
exponential convergence of error.

We drew parameter sets from the prior distributions specified by Table 1, for
4 species trees and 16 species trees with ‘caterpillar’ and ‘well-balanced’
tree shapes. Furthermore for each tree we specified expected tree heights to
be either very small or very large. We then simulated a 1000 sites for each
parameter set and computed the relative error of the likelihood for number of
basis function, $K=5,\ldots,50$. Note that for this experiment we do not limit
$K$ to values of the form $2^{m}+1$. As there is no analytical expression for
the likelihood so to assess convergence we compared values calculated to those
for a large ($K=200$) number of basis functions. We summarize the results for
4 species trees in Figure 1 and for 16 species trees in Figure 2.

|  | Short trees | Tall trees
---|---|---|---
Parameter | Description | Distribution | Mean | Distribution | Mean
$\theta$ | Population size |  $\Gamma$(2, 200) | 0.01 |  $\Gamma$(2, 200) | 0.01
Tree | Shape and branch length | Yule(100) | 0.00111… | Yule(0.1) | 10
Table 1: Priors used to draw trees for error convergence study. Gamma priors
on population sizes and Yule prior on tree shape, height and branch length was
specified.

We see that in all the cases we have that the error decreases exponentially
with $K$, that is, it decreases like $\alpha^{K}$ for some $\alpha<1$. Error
is smaller for long branch lengths since approximate solutions are typically
lower degree polynomials. However when the branch lengths are very short the
error decreases more slowly. There are good reasons for this. Firstly, when
the partial likelihood function is spiked, the approximate solutions are high
degree polynomials requiring more basis functions. Secondly, population sizes
for short branches are intrinsically more difficulty to estimate, no matter
which model or method is used. The only information we have about population
sizes comes from the distribution of coalescent events, and on short branches
there are simply insufficient coalescent events to make sound inference.
Later, we address this by adopting a prior distribution on population sizes
which introduces correlation between neighbouring branches.

Apart from branch length distribution, tree shape does not seem to have a
noticeable effect on the error convergence.

Figure 1: Log scale plot of relative error for basis functions
$K=5,6,\dots,50$ of 4-taxa trees: (a) balanced short; (b) caterpillar short;
(c) balanced tall; (d) caterpillar tall. The sub-linear decrease in log-error
corresponds to an exponential decrease in error, until the limit of machine
precision is reached and no further improvements in error are possible.
Figure 2: Log scale plot of relative error for basis functions
$K=5,6,\dots,50$ of 16-taxa trees: (a) balanced short; (b) caterpillar short;
(c) balanced tall; (d) caterpillar tall. The sub-linear decrease in log-error
corresponds to an exponential decrease in error, until the limit of machine
precision is reached and no further improvements in error are possible.

###  4.2 Comparing inferences from Snapp and Snapper

Theory states that the multispecies coalescent model underlying Snapp and the
diffusion model behind Snapper are approximations of one another (Griffiths
and Spano‘, 2010), with differences decreasing as effective population size
increases. To demonstrate this we analyze the rattlesnake data set in Chifman
and Kubatko (2014). The data 3000 SNPs from a total of 59 individuals in 7
populations.

In Chifman and Kubatko (2014) it was reported that Snapp was not sampling from
the posterior distribution efficiently enough. It quickly became apparent that
the sole cause of difficulty was the inference of population sizes on branches
with close to zero length. We selected the CIR prior in Snapp and in Snapper,
which implements a model with correlated population sizes on nearby branches.
Convergence was also improved by implementing a proposal which scales all
population sizes within a randomly selected subtree.

We ran the MCMC chain for 750,000 iterations for both Snapp and Snapper. It
took Snapp approximately $\sim 56$ hours to run and Snapper approximately
$\sim 4$ hours to run. We give the summary statistics in the appendix (A.3).

Figure 3 uses Densitree (Bouckaert, 2010) to depict samples of the posterior
distributions of the most likely densitrees for Snapp and Snapper. We print
the posterior mean of effective population size above each branch. In both
cases we see that there is only one topology in the 95% highest posterior
density. Posterior distributions of branch lengths are mostly
indistinguishable. There are some small differences in posterior distributions
of populations sizes. However in all cases population sizes follow the same
apparent distributions.

Figure 3: A Snapp and Snapper inference of rattlesnake species trees displayed
using Densitree (Bouckaert, 2010). Rattlesnake populations are numbered from
1-7. Branch thickness is related to relative population sizes, tree height is
reported in expected number of mutations and population sizes is printed below
each related branch, as $\mu\pm 3\sigma$

###  4.3 Running time of Snapper

We compare the computation time between Snapper and Snapp on a 4- and 8-taxa
tree scaling the number of individuals at a leaf. We generated datasets of a
1000 SNPs given the number of species ($s=4,8$) and the number of individuals
at a leaf ($n_{\ell}=5,10,\dots,30$). The number of basis functions for
Snapper was fixed at K = 33. Likelihood computation without caching was done
on a computer with an Intel i3-7100 CPU. In Figure 4 we report the average
times (in seconds) to compute the likelihood of a 1000 sites on 4- and 8-taxa
trees, for Snapp and Snapper. As expected we can clearly see the advantage of
Snapper, in that computation time is hardly affected by the number of
individuals at a leaf. Note that there is some dependence on the number of
individuals. As having large samples can lead to a need for more basis
functions to ensure accuracy.

Figure 4:  Average times (in seconds) to compute the likelihood of a 1000
sites on 4- and 8-taxa trees, for Snapp and Snapper with $s=4,8$ (number of
taxa) and $n_{\ell}=5,10,\dots,30$ (number of individuals at each leaf).

##  5 Analysis of freshwater turtle; Emydura macquarii

To illustrate the application of Snapper we reanalyse SNP data of Georges et
al. (2018) from a group of freshwater turtles known collectively as Emydura.
The range of Emydura extends almost the full length of the Australian
continent from north to south. The group is currently recognized as a complex
of closely related and morphologically distinct allopatric forms, variously
regarded as species, subspecies or distinct morphological lineages (Georges et
al., 2018). The individual samples were taken from 57 distinct water bodies to
study and clarify some of the evolutionary relationships among different
Emydura populations found in Eastern Australia river basins and along the
coast. The sampling covers the coastal drainages of eastern Australia, from
the Hunter River in the south (New South Wales) to the Normanby River
(Queensland) in the north; the rivers of the Murray-Darling Basin (MDB),
including the Paroo drainage, and the Lake Eyre Basin (LEB), and the
intervening Bulloo Basin (Georges et al., 2018).

For the Snapper analysis we used 399 individuals divided into 41 populations
with a total of 5,186 filtered SNPs. We expect short branch lengths for some
of the populations due to geographical proximity of sampling sites. Therefore
we specify a CIR prior that allows for correlation between populations on
nearby branches. The analysis include two out-taxon populations therefore we
expect tree height to be large. Hence we use a Yule prior with appropriate
mean tree height. We have no prior information on mutation rates therefore we
fix $u=v=1$. See Table 2 for more details.

Parameter | Description | Distribution | Mean
---|---|---|---
$\theta$ | Population size | CIR(2, 200,1) | 0.01
Tree | Shape and branch length | Yule(2) | 0.2
Table 2: Priors used in Snapper analysis of freshwater turtle; Emydura
macquarii. CIR priors on population sizes and Yule prior on tree shape, height
and branch length was specified. Note that the mean given for CIR prior is the
mean of the stationary process. Figure 5:  41 population densitree of
freshwater turtle Emydura macquarii. On the x-axis, variation in the tree
represents uncertainty of branch lengths. Thickness of the branches represent
posterior mean of population sizes. Timescale grid at the top is given in
expected number of mutations per lineage per site.

We conducted a number of initial investigations to determine appropriate
proposals and proposal weights. We then ran the sampler for 2,000,000
iterations with convergence assessed in Tracer (Rambaut et al., 2018) using
estimated effective sample size. It took a total of $\sim 500$ hours for the
sampler to run on a computer with an Intel i3-7100. CPU. We provide a complete
list of summary statistics in the appendix (A.3).

Figure 5 depicts the tree with largest posterior probability with uncertainty.
The shape of the tree in Figure 5 agrees with the genetic distances and
SVDquartets trees computed in Georges et al. (2018). We extend the anaylsis in
Georges et al. (2018) by quantifying uncertainty around branch lengths and
include an additional parameter, namely population size. In Figure 5 we see
some uncertainty surrounding topology in the Fitzroy clade (samples 38 - 43).
The Cooper Creek clade (samples 96 - 102) is sister to the North-east Coast
clade (samples 10 - 19) in all three analyses, as is the sister relationship
between the Hunter R (sample 92) and the Murray-Darling Basin (samples 96 -
115) populations. The East Coast (samples 29 - 60) and South-east Coast
(samples 11 - 20) clade is sister to the Hunter-MDB clade (samples 20 - 31) in
the Snapper tree, in agreement with the Fitch-Margoliash distance tree,
whereas the arrangement is unresolved in the SVDquartets tree. Sensibly, the
Kolan-Burnett-Mary (sampels 48 - 58) clade is sister to the Fitzroy clade in
the Snapper tree, whereas the SVDquartets tree has the Mary-Burnett clade in
an unresolved trichotomy with the Cooper-NE Coast clade (samples 96 - 102;
samples 10 - 19) and the Fitzroy clade. The Snapper analysis supports
relationships within the NE Coast clade that reflect drainage proximity better
than does the SVD quartets analysis.

Figure 6:  11 population densitree of freshwater turtle Emydura macquarii
restricted to MDB clade to provide better resolution of branch lenghts and
topology. We show the fraction of trees in the tree set that contain a clade
as text on the graph. We also display the mean of a node as a marker on the
tree. Timescale grid is given in expected number of mutations per lineage per
site.

Also, we note that the bottom two population samples, i.e Emydura subglobosa
and Emydura victoriae populations, are considered as out-taxon groups.
Therefore we expect to see the divergence time from the rest of the
populations to be much earlier. This however leads to poor resolution for the
MDB clade (samples 112 - 131). Thus we present the MDB clade in Figure 6. The
figure also depicts the extent of uncertainty surrounding the tree due to
short branch lengths. As we discuss above, population size estimates here will
be mostly dependent on the prior due to little information available on such
short branch lengths.

##  6 Discussion

In this paper we present Snapper a computationally more efficient method to
supersede Snapp. Like Snapp the method takes biallelic markers sampled from
individuals in multiple populations and computes the likelihood of the species
tree topology together with branch lengths and population sizes. We achieved
this computational efficiency by computing the likelihood using diffusion
models rather than the multispecies coalescent. The Wright-Fisher diffusion
and Coalescent are dual processes and it is therefore not surprising that the
same inference can be made under these distinct but related model frameworks.

We utilize observed allele frequencies in the sample as initial conditions for
the backwards diffusion. The advantages of using the backward diffusion are
two-fold. Firstly the numerical solutions are stable and bounded. Secondly, it
is clear how to define the boundary conditions.

The Snapper sampler is based on the Snapp sampler and uses the same move
proposals, which are standard in the Beast2 software package. However to
improve sampling efficiency for large trees we implemented an additional move
to scale population sizes on subtrees.

We have reported some of the analyses performed in order to validate and more
importantly convince the reader of the ability of Snapper to infer population
genetic parameters. Further work needs to be done to establish regions in the
appropriate parameter space where specific parameters are irrecoverable.

Finally, we note that the general framework for defining likelihoods in terms
of backward diffusions can be extended to other difffusion models. These
models can include additional parameters such as selection, migration and
linkage disequilibrium. In some cases even the computational framework is
reusable.

## Acknowledgements

We thank Arthur Georges for making the SNP dataset of Emydura macquarii
available. Also for his help with the biogeographic interpretation of the
Emydura macquarii species phylogeny and comparison of Snapper analysis with
previous analysis done on the Emydura dataset. Marnus Stolz received a
doctoral scholarship from the NZ Marsden Fund (PIs David Bryant and Steven
Higgins).

## References

  * Boitard et al. (2016) Boitard, S., W. Rodríguez, F. Jay, S. Mona, and F. Austerlitz. 2016. Inferring population size history from large samples of genome-wide molecular data-an approximate Bayesian computation approach. PLoS genetics 12:e1005877.
  * Bollback et al. (2008) Bollback, J. P., T. L. York, and R. Nielsen. 2008. Estimation of 2Nes from temporal allele frequency data. Genetics 179:497–502.
  * Bouckaert (2010) Bouckaert, R. R. 2010. Densitree: making sense of sets of phylogenetic trees. Bioinformatics 26:1372–1373.
  * Bouckaert et al. (2019) Bouckaert, R., T. G. Vaughan, J. Barido-Sottani, S. Duchêne, M. Fourment, A. Gavryushkina, J. Heled, G. Jones, D. Kühnert, N. De Maio and others. 2019\. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS computational biology 15, 4:e1006650.
  * Bryant et al. (2012) Bryant, D., R. Bouckaert, J. Felsenstein, N. A. Rosenberg, and A. RoyChoudhury. 2012\. Inferring Species Trees Directly from Biallelic Genetic Markers: Bypassing Gene Trees in a Full Coalescent Analysis. Molecular Biology and Evolution 29:1917–1932.
  * Bryant and Hahn (2019) Bryant, D. and M. Hahn. 2019. The concatenation question. _in_ Phylogenetics in the Genomics Era (F. Delsuc, N. Galtier, and C. Scornavacca, eds.). (in press).
  * Chifman and Kubatko (2014) Chifman, J. and L. Kubatko. 2014. Quartet inference from SNP data under the coalescent model. Bioinformatics 30:3317–3324.
  * Drummond and Rambaut (2007) Drummond, A. J. and A. Rambaut. 2007. Beast: Bayesian evolutionary analysis by sampling trees. BMC Evolutionary Biology 7:214.
  * Edwards (2009) Edwards, S. V. 2009. Is a new and general theory of molecular systematics emerging? Evolution: International Journal of Organic Evolution 63:1–19.
  * Epstein and Mazzeo (2010) Epstein, C. L. and R. Mazzeo. 2010. Wright–Fisher diffusion in one dimension. SIAM Journal on Mathematical Analysis 42:568–608.
  * Etheridge (2011) Etheridge, A. 2011. Some mathematical models from population genetics : École d’éte de probabilités De Saint-flour XXXIX-2009. Springer.
  * Ethier and Norman (1977) Ethier, S. N. and M. F. Norman. 1977. Error estimate for the diffusion approximation of the Wright–Fisher model. Proceedings of the National Academy of Sciences 74:5096–5098.
  * Ewens (2004) Ewens, W. 2004. Mathematial population genetics. I. Theoretical introduction. Interdisciplinary Applied Mathematics 2 ed. Springer, New York.
  * Felsenstein (1981) Felsenstein, J. 1981. Evolutionary trees from dna sequences: a maximum likelihood approach. Journal of Molecular Evolution 17:368–376.
  * Felsenstein (1992) Felsenstein, J. 1992. Phylogenies from restriction sites: a maximum-likelihood approach. Evolution 46:159–173.
  * Fox and Parker (1968) Fox, L. and I. Parker. 1968. Chebyshev Polynomials in Numerical Analysis. Oxford Mathematical Handbooks Oxford University Press, London.
  * Georges et al. (2018) Georges, A., B. Gruber, G. B. Pauly, D. White, M. Adams, M. J. Young, A. Kilian, X. Zhang, H. B. Shaffer, and P. J. Unmack. 2018. Genomewide SNP markers breathe new life into phylogeography and species delimitation for the problematic short-necked turtles (Chelidae: Emydura) of eastern Australia. Molecular Ecology 27:5195–5213.
  * Griffiths and Spano‘ (2010) Griffiths, R. C. and D. Spano‘. 2010. Diffusion processes and the coalescent. Pages 358–379 _in_ Probability and Mathematical Genetics: Papers in Honour of Sir John Kingman (N. H. Bingham and C. M. Goldie, eds.) London Mathematical Society Lecture Note Series. Cambridge University Press.
  * Gutenkunst et al. (2010) Gutenkunst, R., R. Hernandez, S. Williamson, and C. Bustamante. 2010. Diffusion approximations for demographic inference: DaDi. Nature Precedings 5:1–1.
  * Gutenkunst et al. (2009) Gutenkunst, R. N., R. D. Hernandez, S. H. Williamson, and C. D. Bustamante. 2009\. Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. PLoS Genetics 5.
  * Heled et al. (2013) Heled, J., D. Bryant, and A. J. Drummond. 2013. Simulating gene trees under the multispecies coalescent and time-dependent migration. BMC evolutionary biology 13:44.
  * Hiscott et al. (2016) Hiscott, G., C. Fox, M. Parry, and D. Bryant. 2016. Efficient Recycled Algorithms for Quantitative Trait Models on Phylogenies. Genome biology and evolution 8:1338–1350.
  * Jenkins et al. (2017) Jenkins, P. A., D. Spano, et al. 2017. Exact simulation of the Wright–Fisher diffusion. The Annals of Applied Probability 27:1478–1509.
  * Kingman (1982) Kingman, J. 1982. The coalescent. Stochastic Processes and their Applications 13:235–248.
  * Lapierre et al. (2017) Lapierre, M., A. Lambert, and G. Achaz. 2017. Accuracy of demographic inferences from the site frequency spectrum: the case of the Yoruba population. Genetics 206:439–449.
  * Larget et al. (2010) Larget, B. R., S. K. Kotha, C. N. Dewey, and C. Ané. 2010. BUCKy: gene tree/species tree reconciliation with Bayesian concordance analysis. Bioinformatics 26:2910–2911.
  * Liu (2008) Liu, L. 2008. Best: Bayesian estimation of species trees under the coalescent model. Bioinformatics 24:2542–2543.
  * Liu et al. (2008) Liu, L., D. K. Pearl, R. T. Brumfield, and S. V. Edwards. 2008. Estimating species trees using multiple-allele DNA sequence data. Evolution: International Journal of Organic Evolution 62:2080–2091.
  * Liu et al. (2010) Liu, L., L. Yu, and S. V. Edwards. 2010. A maximum pseudo-likelihood approach for estimating species trees under the coalescent model. BMC Evolutionary Biology 10:302.
  * Lukic and Hey (2012) Lukic, S. and J. Hey. 2012. Demographic inference using spectral methods on SNP data, with an analysis of the human out-of-Africa expansion. Genetics 192:619–639.
  * Mason and Handscomb (2002) Mason, J. C. and D. C. Handscomb. 2002. Chebyshev polynomials. Chapman and Hall/CRC.
  * McKane and Waxman (2007) McKane, A. J. and D. Waxman. 2007. Singular solutions of the diffusion equation of population genetics. Journal of Theoretical Biology 247:849–858.
  * Nielsen (2000) Nielsen, R. 2000. Estimation of population parameters and recombination rates from single nucleotide polymorphisms. Genetics 154:931–942.
  * Øksendal (2003) Øksendal, B. 2003. Stochastic differential equations, an Introduction with Applications. Universitext Springer.
  * Racimo et al. (2016) Racimo, F., G. Renaud, and M. Slatkin. 2016. Joint estimation of contamination, error and demography for nuclear dna from ancient humans. PLoS genetics 12:e1005972.
  * Rambaut et al. (2018) Rambaut, A., A. J. Drummond, D. Xie, G. Baele, and M. A. Suchard. 2018. Posterior summarization in Bayesian phylogenetics using tracer 1.7. Systematic Biology 67:901–904.
  * Sirén et al. (2010) Sirén, J., P. Marttinen, and J. Corander. 2010. Reconstructing population histories from single nucleotide polymorphism data. Molecular Biology and Evolution 28:673–683.
  * Song and Steinrücken (2012) Song, Y. S. and M. Steinrücken. 2012. A simple method for finding explicit analytic transition densities of diffusion processes with general diploid selection. Genetics 190:1117–29.
  * Steinrücken et al. (2014) Steinrücken, M., A. Bhaskar, and Y. S. Song. 2014. A novel spectral method for inferring general diploid selection from time series genetic data. The Annals of Applied Statistics 8:2203–2222.
  * Tataru et al. (2017) Tataru, P., M. Simonsen, T. Bataillon, and A. Hobolth. 2017. Statistical inference in the Wright–Fisher model using allele frequency data. Systematic biology 66:e30–e46.
  * Trefethen (2013) Trefethen, L. N. 2013. Approximation theory and approximation practice. SIAM, Philadelphia.
  * Vachaspati and Warnow (2015) Vachaspati, P. and T. Warnow. 2015. ASTRID: accurate species trees from internode distances. BMC Genomics 16:S3.
  * Waldvogel (2006) Waldvogel, J. 2006. Fast construction of the Fejér and Clenshaw–Curtis quadrature rules. BIT Numerical Mathematics 46:195–202.
  * Williamson et al. (2005) Williamson, S. H., R. Hernandez, A. Fledel-Alon, L. Zhu, R. Nielsen, and C. D. Bustamante. 2005. Simultaneous inference of selection and population growth from patterns of variation in the human genome. Proceedings of the National Academy of Sciences 102:7882–7887.
  * Wright (1931) Wright, S. 1931. Evolution in mendelian populations. Genetics 16:97.
  * Yang (2015) Yang, Z. 2015. The BPP program for species tree estimation and species delimitation. Current Zoology 61:854–865.

##  Appendix A Appendix

###  A.1 The $O(K)$ algorithm for solving diffusions

In the main text, we introduced a matrix $\mathbf{Q}$ which plays a key role
in the numerical solution of diffusions. The matrix is chosen so that if

$g(x,t)=\sum_{k=0}^{K}\bm{\lambda}(t)_{k}T_{k}^{*}(x)$ (22)

then

$(\beta_{1}(1-x)-\beta_{2}x)\frac{d}{dx}g(x,t)+\mbox{$\frac{1}{2}$}x(1-x)\frac{d^{2}}{dx^{2}}g(x,t)\approx\sum_{k=0}^{K}(\mathbf{Q}\bm{\lambda}(t))_{k}T_{k}^{*}(x).$
(23)

The bottleneck in the numerical method we use is the repeated solution of
linear systems that look like

$(\mathbf{Q}-z\mathbf{I})\mathbf{x}=\mathbf{v}$

for difference complex values $z$ and vectors $\mathbf{v}$. Using a direct
method, these take $O(K^{2})$ time each as $\mathbf{Q}$ is lower triangular.
However we can do better, using a trick described in Fox and Parker (1968).
Here we give a very high level description of the approach.

The key idea is to apply integration twice to the LHS of (23). Using
integration by parts multiple times we have

$\displaystyle\int_{0}^{x}\int_{0}^{y}\left[(\beta_{1}(1-z)-\beta_{2}z)\frac{d}{dz}g(z,t)+\mbox{$\frac{1}{2}$}z(1-z)\frac{d^{2}}{dz^{2}}g(z,t)\right]\,dz\,dy$
$\displaystyle=-\int_{0}^{x}\int_{0}^{y}g(z,t)dz+\left(\beta_{1}(1-x)-\beta_{2}x-1+2x\right)\int_{0}^{x}g(z,t)dz$
$\displaystyle\quad\quad+\mbox{$\frac{1}{2}$}x(1-x)g(x,t)+(\mbox{$\frac{1}{2}$}-\beta_{1})xg(0,t).$
(24)

We introduce two new matrices $\mathbf{X}$ and $\mathbf{Y}$. The matrix
$\mathbf{X}$ is derived from properties of shifted Chebyshev polynomials, and
is defined so that if $g$ is expanded as in (22) then

$\displaystyle\int_{0}^{x}\int_{0}^{y}g(z)dz\,dy$
$\displaystyle\approx\sum_{k=0}^{K}(\mathbf{X}\bm{\lambda}(t))_{k}T^{*}_{k}(x)$

The matrix $\mathbf{Y}$ comes from (24) and has the property that if $g$ is
expanded as in (22) then the RHS of (24) equals

$\sum_{k=0}^{K}(\mathbf{Y}\bm{\lambda}(t))_{k}T^{*}_{k}(x).$

We then obtain

$\mathbf{X}\mathbf{Q}\approx\mathbf{Y}.$

The usefulness of this follows from that fact that, with the exception of two
rows, all the non-zero entries in $\mathbf{X}$ and $\mathbf{Y}$ are on or near
the digaonal: both matrices are almost banded. To solve
$(\mathbf{Q}-z\mathbf{I})\mathbf{x}=\mathbf{v}$ we multiply both sides by
$\mathbf{X}$ and solve the sparse system that results. Overall, this now takes
$O(K)$ time.

###  A.2 Further details on the simulation studies

Simulated from | Priors | $\theta_{1}$ | $\theta_{2}$ | $\theta_{3}$ | $\theta_{4}$ | $\theta_{5}$ | $\theta_{6}$ | $\theta_{root}$ | height
---|---|---|---|---|---|---|---|---|---
$\Gamma$(2,200) | Correct | 0.91 | 0.91 | 0.94 | 0.93 | 0.94 | 0.97 | 0.95 | 1.00
$\Gamma$(2,200) | Incorrect | 0.84 | 0.88 | 0.87 | 0.84 | 0.75 | 0.71 | 0.06 | 0.99
$\Gamma$(4,400) | Correct | 0.93 | 0.97 | 0.92 | 0.97 | 0.9 | 0.98 | 0.95 | 1.00
$\Gamma$(4,400) | Incorrect | 0.91 | 0.95 | 0.93 | 0.94 | 0.78 | 0.68 | 0.32 | 0.99
$\Gamma$(80,8000) | Correct | 0.98 | 0.97 | 0.99 | 0.99 | 0.97 | 0.98 | 0.97 | 1.00
$\Gamma$(80,8000) | Incorrect | 0.94 | 0.95 | 0.96 | 0.93 | 0.84 | 0.83 | 0.51 | 1.00
Table 3: Frequency that parameters fall within the $95\%$ highest posterior
density of the MCMC chain for a 100 simulated datasets varying the prior on
the population sizes for the trees from which the data is simulated. Note that
we display the frequencies above as percentages.

We conducted a simulation study to check whether the model posterior
distributions recover model parameters. We simulate data from three different
4 species tree distributions. Parameter means are equal but we change the
variance of the 4 species tree distribution. We then draw a 100 parameter sets
from each 4 species tree distribution. Thereafter we simulate a 1000 SNPs for
32 individuals distributed over 4 species for each parameter set. For every
simulated SNP dataset we ran two chains for 100,000 iterations, one with
correct prior, i.e 4 species tree distribution used for simulation, and one
with incorrect prior. See Table 4 for details on priors. We then count the
number of parameters that were recovered, i.e fall within the 95% highest
posterior density of the MCMC chain. MCMC chains took on average $\sim$200
seconds to run. We report these counts in Table 3. We see that correct priors
on the model has a high rate of parameter recovery. However, when incorrect
priors are used recovery rates suffer for population size parameters close to
the root. This does not come as a surprise due to information loss on branches
as we go up the tree.

|  | Correct prior | Incorrect prior
---|---|---|---
Parameter | Description | Distribution | Mean | Distribution | Mean
$\theta$ | Population size |  $\Gamma$(2, 200) | 0.01 |  $\Gamma$(2, 20) | 0.1
Tree | Topology and branch length | Yule(10) | 0.0111… | Yule(1) | 1
Table 4: We used the priors listed here to conduct the simulation study. In
the case of incorrect priors, we changed the magnitude of both the Gamma and
Yule priors.

###  A.3 Tabulated posterior statistics for Snapp and Snapper analyses

(Note: these could go online or in supplementary data)

Table 5: Snapp parameter summary after 750,000 MCMC iterations for rattlensnake dataset Snapp | mean | variance | HPD | ACT | ESS
---|---|---|---|---|---
theta0 | 0.1121 | 1.5601E-05 | [0.1042,0.1192] | 1122.248 | 590.7785
theta1 | 0.0145 | 5.4621E-07 | [0.0129,0.0158] | 1515.4263 | 437.5007
theta2 | 0.0186 | 9.4307E-07 | [0.017,0.0206] | 3186.245 | 208.0819
theta3 | 0.023 | 1.286E-06 | [0.0208,0.0253] | 2561.4008 | 258.8427
theta4 | 0.0103 | 5.7183E-07 | [8.8272E-3,0.0119] | 3089.9697 | 214.5652
theta5 | 0.0228 | 1.2949E-06 | [0.0209,0.0253] | 2843.4933 | 233.1639
theta6 | 0.01 | 5.3459E-07 | [8.8004E-3,0.0115] | 2988.3364 | 221.8626
theta7 | 0.0196 | 5.7842E-07 | [0.0183,0.0213] | 3545.1482 | 187.0162
theta8 | 0.0248 | 1.0101E-06 | [0.023,0.0269] | 3654.5756 | 181.4164
theta9 | 0.0107 | 3.8089E-07 | [9.4059E-3,0.0118] | 3322.3719 | 199.5562
theta10 | 0.0225 | 7.0178E-07 | [0.0209,0.0241] | 2743.8447 | 241.6317
theta11 | 0.0258 | 7.1985E-07 | [0.0241,0.0274] | 1528.5784 | 433.7363
theta12 | 0.0669 | 1.2663E-05 | [0.0601,0.0735] | 1035.9779 | 639.975
tree.height | 0.0489 | 3.8258E-06 | [0.0452,0.0529] | 1822.8776 | 363.7106
Table 6: Snapper parameter summary after 750,000 MCMC iterations for rattlensnake dataset Snapper | mean | variance | HPD | ACT | ESS
---|---|---|---|---|---
theta0 | 0.1178 | 2.3175E-05 | [0.1095,0.1272] | 1076.1399 | 469.2698
theta1 | 0.0117 | 4.4551E-07 | [0.0104,0.0131] | 1833.7797 | 275.3875
theta2 | 0.0176 | 1.1315E-06 | [0.0156,0.0195] | 1830.8958 | 275.8213
theta3 | 0.0213 | 1.4143E-06 | [0.019,0.0236] | 1664.7158 | 303.3551
theta4 | 0.0068167 | 4.6294E-07 | [5.311E-3,7.9448E-3] | 1783.8144 | 283.1012
theta5 | 0.0208 | 1.4231E-06 | [0.0185,0.023] | 2257.7471 | 223.6743
theta7 | 0.0193 | 6.8242E-07 | [0.0179,0.0211] | 2945.9165 | 171.4237
theta8 | 0.0244 | 1.1309E-06 | [0.0222,0.0263] | 2591.8036 | 194.845
theta9 | 0.0077304 | 2.493E-07 | [6.8474E,3-8.7642E-3] | 2657.2438 | 190.0465
theta10 | 0.0228 | 9.5682E-07 | [0.0208,0.0247] | 2374.8237 | 212.6474
theta11 | 0.028 | 1.0207E-06 | [0.0261,0.03] | 1396.2835 | 361.6744
theta12 | 0.08 | 1.8343E-05 | [0.0713,0.0873] | 1000 | 505
tree.height | 0.0449 | 1.4012E-05 | [0.0413,0.0486] | 1341.1584 | 376.5402
Table 7: Snapper parameter summary after 2,000,000 MCMC iterations for freshwater turtle dataset Snapper | mean | variance | HPD | ACT | ESS
---|---|---|---|---|---
pdist | -32288 | 58.555 | [-32289,-32286] | 6921.9588 | 237.4934
theta0 | 0.0892 | 1.09E-05 | [0.0829,0.0956] | 4559.0819 | 360.5814
theta1 | 0.0326 | 2.03E-06 | [0.0299,0.0354] | 3334.7123 | 492.972
theta2 | 0.0326 | 2.31E-06 | [0.0297,0.0355] | 4882.2185 | 336.7158
theta3 | 0.0315 | 2.28E-06 | [0.0288,0.0347] | 4594.6906 | 357.7868
theta4 | 0.089 | 1.16E-05 | [0.0823,0.0955] | 5834.7779 | 281.745
theta5 | 0.0209 | 7.38E-07 | [0.0195,0.0228] | 5032.2875 | 326.6746
theta6 | 0.0209 | 7.52E-07 | [0.0193,0.0227] | 5336.7619 | 308.037
theta7 | 0.0211 | 7.42E-07 | [0.0195,0.0227] | 5739.448 | 286.4248
theta8 | 0.0568 | 5.41E-06 | [0.0525,0.0613] | 7592.2756 | 216.5254
theta9 | 0.0576 | 5.36E-06 | [0.0533,0.0622] | 5374.5587 | 305.8706
theta10 | 0.0576 | 5.40E-06 | [0.0535,0.0621] | 5700.7353 | 288.3698
theta11 | 0.0398 | 3.51E-06 | [0.0362,0.0437] | 4824.2314 | 340.763
theta12 | 0.0475 | 2.82E-06 | [0.0439,0.0506] | 5591.2287 | 294.0176
theta13 | 0.0474 | 3.08E-06 | [0.0441,0.0508] | 7868.6887 | 208.9192
theta14 | 0.0473 | 3.13E-06 | [0.0441,0.0509] | 6942.6134 | 236.787
theta15 | 0.0474 | 2.99E-06 | [0.0441,0.0507] | 7085.1136 | 232.0246
theta16 | 0.0473 | 2.84E-06 | [0.0442,0.0503] | 4927.4417 | 333.6254
theta17 | 0.0563 | 5.91E-06 | [0.0512,0.0608] | 6188.3147 | 265.649
theta18 | 0.0254 | 1.90E-06 | [0.0229,0.0283] | 4461.407 | 368.4756
theta19 | 0.021 | 1.58E-06 | [0.0186,0.0233] | 3920.8397 | 419.2776
theta20 | 0.0351 | 2.47E-06 | [0.032,0.0381] | 5149.6604 | 319.2288
theta21 | 0.035 | 2.45E-06 | [0.032,0.0382] | 5103.4682 | 322.1182
theta22 | 0.0312 | 1.98E-06 | [0.0286, 0.034] | 5843.4713 | 281.326
theta23 | 0.0275 | 1.68E-06 | [0.0249,0.0298] | 4633.1697 | 354.8154
theta24 | 0.0137 | 5.95E-07 | [0.0121,0.0151] | 6877.4641 | 239.03
theta25 | 0.0137 | 6.94E-07 | [0.0121,0.0152] | 7605.5597 | 216.1472
theta26 | 0.0105 | 3.66E-07 | [9.3707E-3,0.0116] | 5984.3981 | 274.701
theta27 | 0.035 | 1.34E-06 | [0.0328,0.0371] | 4350.192 | 377.896
theta28 | 0.0324 | 1.52E-06 | [0.0301,0.0349] | 6937.3698 | 236.9658
theta29 | 0.0324 | 1.51E-06 | [0.03,0.0347] | 6434.1179 | 255.5004
theta30 | 0.033 | 1.42E-06 | [0.0305,0.0352] | 5917.0715 | 277.8266
theta31 | 0.033 | 1.34E-06 | [0.0304,0.035] | 5047.0492 | 325.719
theta32 | 0.0346 | 1.40E-06 | [0.0325,0.0371] | 5272.0715 | 311.8168
theta33 | 0.0345 | 1.55E-06 | [0.0321,0.0367] | 4934.2925 | 333.1622
theta34 | 0.0346 | 1.52E-06 | [0.0322,0.037] | 5797.7463 | 283.5446
theta35 | 0.0345 | 1.39E-06 | [0.0319,0.0366] | 5713.2156 | 287.7398
theta36 | 0.0343 | 1.49E-06 | [0.0321,0.0367] | 5738.9424 | 286.45
theta37 | 0.032 | 1.73E-06 | [0.0296,0.0346] | 5566.4597 | 295.326
theta38 | 8.24E-03 | 1.83E-07 | [7.3532E-3,9.0356E-3] | 4660.0198 | 352.771
theta39 | 0.0285 | 2.06E-06 | [0.0254,0.031] | 2754.3786 | 596.8388
theta40 | 0.045 | 3.94E-06 | [0.0416,0.0493] | 2735.1968 | 601.0244
theta41 | 0.0891 | 1.03E-05 | [0.0823,0.0946] | 4363.7001 | 376.7262
theta42 | 0.0328 | 1.74E-06 | [0.0302,0.0352] | 4046.0439 | 406.303
theta43 | 0.0344 | 1.63E-06 | [0.0321,0.037] | 4026.7212 | 408.2528
theta44 | 0.0553 | 2.70E-06 | [0.0522,0.0585] | 7152.4991 | 229.8386
theta45 | 0.0209 | 7.09E-07 | [0.0196,0.023] | 4865.2479 | 337.8902
theta46 | 0.021 | 6.24E-07 | [0.0195,0.0226] | 4733.2039 | 347.3166
theta47 | 0.0518 | 4.84E-06 | [0.0481,0.0555] | 26723.5617 | 61.5158
theta48 | 0.0575 | 4.88E-06 | [0.0538,0.0621] | 5530.8802 | 297.2258
theta49 | 0.0573 | 4.28E-06 | [0.0535,0.0612] | 6291.947 | 261.2736
theta50 | 0.0567 | 3.64E-06 | [0.0531,0.0607] | 7435.9237 | 221.0782
theta51 | 0.0475 | 2.81E-06 | [0.0441,0.0506] | 7190.8831 | 228.6116
theta52 | 0.0474 | 2.92E-06 | [0.0442,0.0506] | 6880.0099 | 238.9416
theta53 | 0.0473 | 2.65E-06 | [0.0444,0.0507] | 5266.9057 | 312.1226
theta54 | 0.0473 | 2.64E-06 | [0.044,0.0502] | 5930.9334 | 277.1772
theta55 | 0.0452 | 2.35E-06 | [0.0423,0.0481] | 6519.7251 | 252.1456
theta56 | 0.0531 | 2.40E-06 | [0.0502,0.0562] | 10545.795 | 155.884
theta57 | 0.0246 | 1.12E-06 | [0.0225,0.0265] | 6463.5958 | 254.3352
theta58 | 0.0351 | 2.26E-06 | [0.0324,0.0382] | 4395.8089 | 373.9744
theta59 | 0.0348 | 1.46E-06 | [0.0326,0.0372] | 12531.3457 | 131.1846
theta60 | 0.0387 | 1.12E-06 | [0.0366,0.0405] | 16803.2068 | 97.8338
theta61 | 0.0276 | 1.16E-06 | [0.0254,0.0294] | 11382.032 | 144.4312
theta62 | 0.0136 | 5.75E-07 | [0.0122,0.015] | 7378.0828 | 222.8112
theta63 | 0.0132 | 3.93E-07 | [0.012,0.0144] | 9066.0779 | 181.3264
theta64 | 0.025 | 7.92E-07 | [0.0232,0.0266] | 13972.7813 | 117.6516
theta65 | 0.0492 | 2.12E-06 | [0.0465,0.052] | 24561.4047 | 66.931
theta66 | 0.0345 | 1.30E-06 | [0.0324,0.0367] | 4752.1445 | 345.9322
theta67 | 0.0344 | 1.22E-06 | [0.0323,0.0366] | 4528.4154 | 363.0232
theta68 | 0.0345 | 1.28E-06 | [0.0324,0.0368] | 5549.8598 | 296.2092
theta69 | 0.0344 | 1.21E-06 | [0.0323,0.0364] | 4893.2052 | 335.9598
theta70 | 0.0343 | 1.13E-06 | [0.0323,0.0363] | 4619.4067 | 355.8726
theta71 | 0.0324 | 1.44E-06 | [0.0301,0.0346] | 6038.6503 | 272.233
theta72 | 0.033 | 1.33E-06 | [0.0304,0.035] | 5481.3538 | 299.9114
theta73 | 0.0329 | 1.25E-06 | [0.0305,0.0349] | 5176.6725 | 117.563
theta74 | 0.0336 | 1.09E-06 | [0.0317,0.0357] | 4940.5417 | 332.7408
theta75 | 0.0334 | 1.06E-06 | [0.0314,0.0354] | 5603.2564 | 293.3866
theta76 | 0.037 | 3.01E-06 | [0.0336,0.0401] | 22711.2203 | 72.3836
theta77 | 0.0492 | 4.41E-06 | [0.0456,0.0533] | 33200.5358 | 49.5148
theta78 | 0.06 | 5.72E-06 | [0.0558,0.064] | 34594.3404 | 47.52
theta79 | 0.0383 | 9.33E-06 | [0.0328,0.0438] | 25502.3797 | 64.4614
theta80 | 0.0531 | 3.43E-06 | [0.0495,0.0565] | 33045.9805 | 49.7464
tree.height | 0.1982 | 3.96E-05 | [0.1863,0.2107] | 1811.7552 | 907.3632
tree.length | 0.754 | 2.41E-04 | [0.7242,0.7845] | 2135.0573 | 769.9652


# Solution of the Basel problem using the Feynman integral trick

F. L. Freitas<EMAIL_ADDRESS>
###### Abstract.

Euler’s solution in 1734 of the Basel problem, which asks for a closed form
expression for the sum of the reciprocals of all perfect squares, is one of
the most celebrated results of mathematical analysis. In the modern era,
numerous proofs of it have been produced, each emphasizing a different style
of calculation, as a way of testing the power of some demonstration method.
It’s often thought that solutions using calculus need to involve advanced
contour integration techniques or geometric coordinate transformations. We
show that this is not the case, as the result can be derived by analyzing
basic properties of a particular one-dimensional integral, and as such, can be
obtained with techniques typical of regular calculus tests and math
competitions.

##  1\. Introduction

Solving the Basel problem amounts to proving the following identity:

(1)
$\sum_{n=1}^{\infty}\frac{1}{n^{2}}=1+\frac{1}{2^{2}}+\frac{1}{3^{2}}+\frac{1}{4^{2}}+\dots=\frac{\pi^{2}}{6}$

The problem was first solved by Euler in 1734[1]. Ever since his famous
original proof, which relied on a generalized Vieta formula for the roots of
the function $\frac{\sin x}{x}$, numerous other proofs have appeared[2], most
emphasizing a different method to attack the problem.

Of all of these, the most elementary is perhaps that of Apostol[3], where the
problem is reduced to evaluating a two-dimensional integral:

(2)
$\sum_{n=1}^{\infty}\frac{1}{n^{2}}=\int_{0}^{1}\int_{0}^{1}\frac{\mathrm{d}x\mathrm{d}y}{1-xy}$

Apostol noted that, if we perform the 2D coordinate transformation
$(u,v)=((x+y)/2,(x-y)/2)$, the integral can be evaluated exactly, and the
$\pi^{2}/6$ answer emerges.

We give a proof similar in spirit to the one of Apostol, but arguably simpler,
because it uses only a 1D integral and instead of an _ad hoc_ coordinate
transformation, it employs differentiation under the integral sign, a
technique popularized by physicist Richard Feynman[4].

##  2\. Result

We attempt to evaluate the integral:

(3) $\int_{0}^{\infty}\ln\left(1+e^{-2x}\right)\mathrm{d}x$

However, instead of evaluating it as is, we solve a more general problem,
where the integral has an extra parameter $\alpha$:

(4) $I(\alpha)=\int_{0}^{\infty}\ln\left(1+\alpha
e^{-x}+e^{-2x}\right)\mathrm{d}x$

Although it looks like an extra complication, the idea is that the dependency
of the integral on $\alpha$ can be very simple, and can be deduced from a
differential equation in the parameter $\alpha$, which might be easier to
solve than the original integral. If we differentiate, we obtain:

(5)
$\frac{\mathrm{d}I}{\mathrm{d}\alpha}=\int_{0}^{\infty}\frac{e^{-x}}{1+\alpha
e^{-x}+e^{-2x}}\mathrm{d}x$

This integral can be simplified by changing variables to $u=e^{-x}$, where we
get $\mathrm{d}u=-e^{-x}\mathrm{d}x$. We can get rid of the minus sign by
swapping the integration limits. Renaming $u\to x$, we can write:

(6) $\frac{\mathrm{d}I}{\mathrm{d}\alpha}=\int_{0}^{1}\frac{1}{1+\alpha
x+x^{2}}\mathrm{d}x$

Since are interested in the range $-2<\alpha<2$, the integrand can be
evaluated with standard techniques, by completing the square in the
denominator. The answer is:

(7)
$\frac{\mathrm{d}I}{\mathrm{d}\alpha}=\left.\left[\frac{2}{\sqrt{4-\alpha^{2}}}\arctan\left(\frac{\alpha+2x}{\sqrt{4-\alpha^{2}}}\right)\right]\right|_{x=0}^{1}$

Substituting the limits:

(8)
$\frac{\mathrm{d}I}{\mathrm{d}\alpha}=\frac{2}{\sqrt{4-\alpha^{2}}}\left(\arctan\left(\frac{\alpha+2}{\sqrt{4-\alpha^{2}}}\right)-\arctan\left(\frac{\alpha}{\sqrt{4-\alpha^{2}}}\right)\right)$

The expression can be simplified with the $\arctan$ identity:

(9) $\arctan(x)-\arctan(y)=\arctan\left(\frac{x-y}{1+xy}\right)$

We obtain:

(10)
$\frac{\mathrm{d}I}{\mathrm{d}\alpha}=\frac{2}{\sqrt{4-\alpha^{2}}}\arctan\left(\sqrt{\frac{2-\alpha}{2+\alpha}}\right)$

Although the expression looks formidable, it can be cast into an elementary
form by substituting $\alpha=2\cos u$. The argument of $\arctan$ simplifies
with the following trigonometric identity:

(11) $\sqrt{\frac{1-\cos u}{1+\cos
u}}=\sqrt{\frac{2\sin^{2}(u/2)}{2\cos^{2}(u/2)}}=\tan\left(\frac{u}{2}\right)$

And the end result is:

(12) $I(\alpha)=-\frac{1}{2}\arccos\left(\frac{\alpha}{2}\right)^{2}+c$

We only need to determine the integration constant $c$. To do that, notice
that we have a factorization in the integrand when $\alpha=2$:

(13)
$I(2)=\int_{0}^{\infty}\ln\left(1+2e^{-x}+e^{-2x}\right)\mathrm{d}x=2\int_{0}^{\infty}\ln\left(1+e^{-x}\right)\mathrm{d}x$

By applying the transformation $x\to 2x$, we can reduce it to the integral at
$\alpha=0$:

(14)
$2\int_{0}^{\infty}\ln\left(1+e^{-x}\right)\mathrm{d}x=4\int_{0}^{\infty}\ln\left(1+e^{-2x}\right)\mathrm{d}x=4I(0)$

The relation $I(2)=4I(0)$ allows us to determine $c$. We have:

(15) $I(2)=c\qquad I(0)=-\frac{\pi^{2}}{8}+c$

Therefore:

(16) $I(2)=4I(0)\Rightarrow c=-\frac{\pi^{2}}{2}+4c\Rightarrow
c=\frac{\pi^{2}}{6}$

Now we are able to determine (4) for any value of $\alpha$. In particular, we
can calculate it for $\alpha=-2$. The corresponding integral is:

(17)
$I(-2)=\int_{0}^{\infty}\ln\left(1-2e^{-x}+e^{-2x}\right)\mathrm{d}x=2\int_{0}^{\infty}\ln\left(1-e^{-x}\right)\mathrm{d}x$

Substituting $\alpha=-2$ into (12), we obtain:

(18)
$I(-2)=-\frac{\pi^{2}}{3}\Rightarrow\int_{0}^{\infty}\ln\left(1-e^{-x}\right)\mathrm{d}x=-\frac{\pi^{2}}{6}$

From the geometric series expansion, we have:

(19) $\frac{1}{1-x}=1+x+x^{2}+x^{3}+\dots=\sum_{n=0}^{\infty}x^{n}$

Integrating (19) term by term we find:

(20)
$-\ln(1-x)=x+\frac{x^{2}}{2}+\frac{x^{3}}{3}+\dots=\sum_{n=1}^{\infty}\frac{x^{n}}{n}$

Substituting $e^{-x}$ and integrating we find:

(21)
$-\int_{0}^{\infty}\ln(1-e^{-x})\mathrm{d}x=\sum_{n=1}^{\infty}\int_{0}^{\infty}\frac{e^{-nx}}{n}\mathrm{d}x=\sum_{n=1}^{\infty}\frac{1}{n^{2}}$

Together, equations (18) and (21) imply our desired result:

(22) $\sum_{n=1}^{\infty}\frac{1}{n^{2}}=\frac{\pi^{2}}{6}$

##  3\. Conclusion

It has been shown that the solution to the classic Basel problem can be found
in a straightforward manner using the Feynman integral trick. The method
applied here is quite unusual, because in most applications of differentiation
under the integral sign, one obtains the integration constant by going to a
limit where the integral is trivial.

Here, the determination of the constant is done by looking at two
configurations of the parametrized integral related by a symmetry. This is a
quite interesting technique, which might have applications in other integrals.
The difficulty of applying the Feynman trick rests, as always, on how to
introduce the parameter to simplify the differentiated integral. However, like
many of these tricks, when it’s seen once, it’s much easier to apply it in new
settings.

In summary, we have provided a very simple proof of the Basel problem which
relies on the Feynman trick in a one-dimensional integral. Although hard to
find, after the parameter is introduced, the proof follows naturally, which
means such a technique might be useful in situations where a solution has to
be produced in tight time constraints, as is often the case in math
competitions.

## References

  * [1] R. Ayoub, “Euler and the zeta function,” The American Mathematical Monthly, vol. 81, no. 10, pp. 1067–1086, 1974.
  * [2] R. Chapman, “Evaluating $\zeta$(2),” Preprint. Available at https://empslocal.ex.ac.uk/people/staff/rjchapma/etc/zeta2.pdf, 2003.
  * [3] T. M. Apostol, “A proof that euler missed: evaluating $\zeta$ (2) the easy way,” The Mathematical Intelligencer, vol. 5, no. 3, pp. 59–60, 1983.
  * [4] R. Feynman and R. Leighton, "Surely You’re Joking, Mr. Feynman!": Adventures of a Curious Character.  W. W. Norton, 2010.