# Forward-backward algorithms with a biallelic mutation-drift model:
Orthogonal polynomials, and a coalescent/urn-model based approach

Claus Vogl<EMAIL_ADDRESS>Sandra Peer<EMAIL_ADDRESS>Lynette Caitlin Mikula<EMAIL_ADDRESS>Department of Biomedical
Sciences, Vetmeduni Vienna, Veterinärplatz 1, A-1210 Wien, Austria  Vienna
Graduate School of Population Genetics, A-1210 Wien, Austria  Centre for
Biological Diversity, School of Biology, University of St. Andrews, St Andrews
KY16 9TH, UK

###### Abstract

Inference of the marginal likelihood of sample allele configurations using
backward algorithms yields identical results with the Kingman coalescent, the
Moran model, and the diffusion model (up to a scaling of time). For inference
of probabilities of ancestral population allele frequencies at any given point
in the past - either of discrete ancestral allele configurations as in the
coalescent, or of ancestral allele proportions as in the backward diffusion -
backward approaches need to be combined with corresponding forward ones. This
is done in so-called forward-backward algorithms. In this article, we utilize
orthogonal polynomials in forward-backward algorithms. They enable efficient
calculation of past allele configurations of an extant sample and
probabilities of ancestral population allele frequencies in equilibrium and in
non-equilibrium. We show that the genealogy of a sample is fully described by
the backward polynomial expansion of the marginal likelihood of its allele
configuration.

###### keywords:

biallelic mutation-drift model , forward-backward algorithm , orthogonal
polynomials , coalescent model , urn model , boundary-mutation Moran model.

††journal: Theoretical Population Biology

##  1 Introduction

In population genetics, extant (present-time) data are usually used for
inference of past population genetic processes. The coalescent [16] is a
stochastic process that describes the genealogy of a sample from a single
locus back to the last common ancestor. It allows for convenient simulations
of genealogical trees, conditional on the current sample size and either the
current (effective) population size or past (effective) population sizes.
Along the branches of this tree, mutations may introduce new allelic states.
The coalescent has become a pillar of theoretical and empirical population
genetics (see [e.g., 13, 34] and e.g., the software ms [14]).

While the extant sample is assumed given, many intermediate configurations are
possible until the last common ancestor of the sample is reached. Mutations
increase this complexity further. Coalescent simulations have been used for
inferring parameters of relatively complex population genetic models assuming
high-throughput population genetic data. The combinatorial complexity of the
coalescent with mutation means that sufficient statistics are often not
available. Summary statistics and Approximate Bayesian Computation (ABC)
[e.g., 2, 8, 12] are then employed for inference of population genetic
parameters. Summary statistics are generally not sufficient and ABC is so
computationally demanding that often only subsets of the parameter space can
be investigated. Inference therefore becomes approximate.

The infinite sites model [15] is another classic population genetic model; one
of its assumptions is complete linkage disequilibrium, i.e., no recombination
[e.g., 35]. Using a forward algorithm (a dynamic programming method) with the
underlying logic of an urn model, Wu [37] showed that inference under the
infinite sites model is efficient in analysis of data sets from panmictic or
subdivided populations in equilibrium. Forward in time, the sample sizes
leading to the extant sample size are random variables. Therefore, the forward
approach cannot be extended to account for deviations from equilibrium, e.g.,
via changing population sizes. Faisal et al. [11] introduced a corresponding
backward algorithm; since it proceeds backward in time conditional on the
extant sample, modelling changing population sizes is possible.

Note that data from autosomes of some model species are further from complete
linkage disequilibrium, which is assumed by the infinite sites model, than
from linkage equilibrium: In flies of the genus Drosophila, deviations from
linkage equilibrium are barely noticeable in genomic regions of moderate to
high recombination rates [22] and weak enough to be negligible even within
short introns [7]. In great apes, recombination rates are considerably more
variable across the genome: while on average neither linkage equilibrium nor
complete disequilibrium can be assumed, either assumption may hold
approximately in certain genomic regions (see Supplementary Material, Figure
S3 in [21]).

In the protein coding genes of eukaryotes, the expected heterozygosity (which
is roughly equal to the scaled mutation rate) is approximately $\leq 10^{-2}$
or smaller [17]. This makes it unlikely to find more than one allele
segregating in a small to moderately sized sample. Thus bi-allelic mutation
models suffice to capture the true dynamics of a population. In genomic
regions where the recombination rate is much greater than the mutation rate,
polymorphic sites can additionally be considered independent: The distance
between polymorphic loci in small to moderately sized samples is great enough
for recombination to break their association. In other words, neighboring
polymorphic sites have different genealogies due to recombination. Data from
such populations can be represented in a site frequency spectrum (SFS) without
loss of information (and data from multiple populations in a joint site
frequency spectrum (jSFS)). The SFS records the frequency of the focal allele
in a sample. A bi-allelic mutation model that can be represented in this way
can also be re-parametrized to a parent-independent mutation model; this means
that both mutation and coalescent events are uninformative regarding the
immediately ancestral allelic state.

In this article, we assume data that can be represented as site frequency
spectra without loss of information. We further assume that data are generated
by a haploid bi-allelic Moran model of population size $N$ or by the
corresponding diffusion model. (Note that the Wright-Fisher model has the same
diffusion limit as the Moran model.) In Sec. (2), we first review the discrete
time and discrete space decoupled mutation-drift Moran model and the
derivation of its corresponding diffusion equation. Then we describe two known
forward-backward algorithms previously applied by Bergman et al. [3]. Such
algorithms enable inference of probabilities of ancestral population allele
frequencies at any time in the past in addition to the probabilities of extant
sample configurations.

The first of these algorithms is for the discrete case: Essentially, the
algorithm of Rabiner and Juang [23] can be directly applied [3]. Our
representation of this algorithm is based on matrix multiplications and the
population size is assumed constant over time. Changes in mutation-drift
parameters can however be immediately incorporated and changes in the
effective population size can be indirectly modeled through rescaling of time.

The second of these algorithm is for the continuous case: An analytical form
of the transition density function is found by decomposition of the diffusion
generator into eigenvectors of orthogonal polynomials (more precisely,
modified Jacobi polynomials) and their corresponding eigenvectors [25, 3]. The
prior distributions of both the ancestral population allele proportions and
the probabilities of extant sample allele configurations can then be extended
across time. This orthogonal polynomial approach is computationally efficient
compared to the matrix multiplications of the discrete case. Again, changes in
mutation-drift parameters can be accommodated, as well as demographic changes
through concurrent rescaling of time. These modifications necessitate a change
in the base of the orthogonal polynomials through linear transformation. Such
transformations can cause numerical inaccuracies.

In Sec. (3), we introduce our so-called particle model: Using coalescent
arguments, one can trace the probability of sample configurations backward
through the sample history similarly to the algorithm assuming complete
linkage and the infinite sites model [11]. By conditioning on the extant
sample, changing demography can be accounted for. Augmenting with temporal
dynamics, the joint probability of the number of focal alleles in the sample
and the past sample sizes can be determined for every point in time:
Essentially, the history of the extant sample can be inferred. We show that
the marginal distribution of the extant sample configuration derived using
this backward particle approach is equivalent to that determined using the
backward orthogonal polynomial approach. The sample genealogy can therefore be
considered embedded in the spectral decomposition of the backwards diffusion
generator similar to how it can be considered embedded within the discrete
Moran model [9, chapt. 2.8].

Assuming equilibrium and reversing time in the coalescent arguments, we obtain
a forward particle approach that can also be used to determine the
probabilities of sample configurations. It essentially corresponds to an urn
model [26] and is again similar to the algorithm assuming complete linkage
equilibrium and the infinite sites model [37]. Recall that running the
particle model forward in time means that past sample sizes become random
variables and modelling non-equilibrium is not feasible. For a full forward-
backward algorithm using the particle model, the backward particle approach
must therefore be combined with the forward orthogonal polynomial approach.

In Sec. (4), we discuss different forward-backward algorithms for the boundary
mutation-drift Moran model [30], which is a simplification of the general
mutation Moran model for small scaled mutation rates. Inference, especially in
non-equilibrium scenarios, becomes particularly efficient using orthogonal
polynomials because the change of base required in the general mutation model
can be avoided [29]. The corresponding boundary mutation particle model, which
we introduce in this article, allows for simple derivation of non-equilibrium
transition probabilities.

##  2 Population allele proportions in the Moran and diffusion models

We assume a haploid population of size $N$ evolving according to a discrete
space, bi-allelic, reversible, decoupled mutation-drift Moran model [18, 19,
20, 9, section 2.8]. Evolution proceeds stepwise, either by a birth-death
event during which a randomly chosen individual is replaced the randomly
chosen offspring of another or by an individual mutating. In Sec. (2.1) and
Sec. (2.2) we establish our notation for the time discrete and diffusion
versions of this Moran model respectively. The diffusion generator can be
decomposed into eigenvectors (spectral decomposition) of modified Jacobi
polynomials with corresponding eigenvalues. In other words, the transition
density function has an explicit spectral representation; we review this in
Sec. (2.2.1). To infer population genetic parameters under these models, the
marginal distributions of population allele frequencies have to be determined
conditional on an extant sample. In the discrete case, a classic forward-
backward algorithm on the transition rate matrix can be employed (Sec. 2.3.1).
In the diffusion limit, the Jacobi polynomials allow for efficient
representation of the result (Sec. 2.3.2).

###  2.1 The discrete decoupled Moran model

We consider a discrete time and discrete space, reversible, and decoupled
mutation-drift Moran model with haploid population size $N$ [9, section 2.8]
to describe the evolution of the proportions of a focal allele $1$ and a non-
focal allele $0$. We assume that birth-death events initially occur at rate
$1$, while mutations arise at rates $\mu_{1}$ towards and $\mu_{0}$ away from
the focal allele (so that the total mutation rate is $\mu=\mu_{0}+\mu_{1}$).
We re-parameterize by setting the mutation bias towards the focal allele to
$\alpha=\mu_{1}/(\mu_{0}+\mu_{1})$ with $0<\alpha<1$, and equivalently define
$\beta=1-\alpha$. Fig. (1) provides a visualisation for $N=6$ individuals, of
which $y=4$ are of the focal type at the present time $s=0$.

Figure 1: Schematic plot of a Moran model: The allelic type of the extant
individuals is given at the top of the plot; backward in time the arrows
indicate birth-death events and the ’x’ mark mutation events. Note that the
position and numbering of individuals on the x-axis is arbitrary.

Let $x_{s}$ ($0\leq x_{s}\leq 1$) denote the relative frequency of allele $0$
at a non-focal locus at time $s$. The numbering of individuals is arbitrary.
The transition rate matrix $\mathbf{T}_{i,j}$ is tridiagonal, aperiodic, and
right stochastic:

$\mathbf{T}_{i,j}=\Pr(x_{s+1}=\tfrac{j}{N}{\,|\,}x_{s}=\tfrac{i}{N})=\begin{cases}x_{s}(1-x_{s})+\beta\mu
x_{s}&\text{for j=i-1}\\\ 1-2x_{s}(1-x_{s})-\beta\mu x_{s}-\alpha\mu
x_{s}&\text{for j=i }\\\ x_{s}(1-x_{s})+\alpha\mu x_{s}&\text{for j=i+1}\\\
0\qquad&\text{otherwise}\,.\end{cases}$ (1)

Following convention, we re-scale the Moran events to bring the birth-death
rate to $N^{2}$, and equivalently time becomes $t=s/N$ (note that this differs
from the otherwise identical treatment in [9, formula 2.10], where time is
scaled by $\binom{N}{2}$ to match the rate of the coalescent in the Wright-
Fisher model). We also set the overall scaled mutation rate to
$\theta=N(\mu_{0}+\mu_{1})$. Then, the law of total probability allows us to
write the forward N-particle generator of the Moran process as [3, formula
75]:

$\begin{split}{\cal
L}_{N}\Pr(Nx_{t})&=\alpha\theta\big{(}(N-i+1)\Pr(Nx_{t}=i-1)-(N-i)\Pr(Nx_{t}=i)\big{)}\\\
\qquad+\beta\theta&\big{(}(i+1)\Pr(Nx_{t}=i+1)-i\Pr(Nx_{t}=i)\big{)}\\\
\qquad+&\big{(}(i-1)(N-i+1)\Pr(Nx_{t}=i-1)\\\
&+(i+1)(N-i-1)\Pr(Nx_{t}=i+1)-2i(N-i)\Pr(Nx_{t}=i)\big{)}\,.\end{split}$ (2)

Note that that the terms have been collected so that the first two summands
account for mutation events and the remaining for genetic drift.

Suppose a sample of size $K$, with $0\leq K\leq N$, is drawn from the
population at the current time $t=0$. The likelihood of observing $y$ of the
focal alleles in the sample follows a hypergeometric distribution:

$\Pr(y{\,|\,}K,N,Nx=i,t=0)=\frac{\binom{i}{y}\binom{N-i}{K-y}}{\binom{N}{K}}\,.$
(3)

Calculation of sample allele configurations at previous times requires the
backward N-particle generator (note the generator ${\cal L}_{N}^{{}^{\prime}}$
operating on the sample allele frequency in Eq. (3) below is defined
equivalently to the forward generator ${\cal L}_{N}$ in Eq. (2) acting on the
transition density from Eq. (1)):

$\begin{split}&{\cal L}_{N}^{{}^{\prime}}\Pr(y{\,|\,}K,Nx=i,t)=\\\
&\qquad\alpha\theta(N-i)\bigg{(}\Pr(y{\,|\,}K,Nx=i+1,t)-\Pr(y{\,|\,}K,Nx=i,t)\bigg{)}\\\
&\qquad+\beta\theta
i\bigg{(}\Pr(y{\,|\,}K,Nx=i-1,t)-\Pr(y{\,|\,}K,Nx=i,t)\bigg{)}\\\
&\qquad+i(N-i)\bigg{(}\Pr(y{\,|\,}K,Nx=i+1,t)+(\Pr(y{\,|\,}K,Nx=i-1,t)\\\
&\qquad\quad-2\Pr(y{\,|\,}K,Nx=i,t)\bigg{)}\,.\end{split}$ (4)

###  2.2 The diffusion limit

Passing to the diffusion limit $N\to\infty$ is comparatively straightforward
within the decoupled Moran model _vs_ the classic Moran model (see [10, chapt.
4] for the derivation assuming the classic Wright-Fisher model, rescaling
suffices to obtain equivalent results for the Moran model): Set $\delta x=1/N$
and denote the transition rate density of continuous allele proportions as
$\phi(x{\,|\,}t)$. Then, the continuous forward N-particle generator is
(compare: Eq. (2) [3, formula 76]:

$\begin{split}{\cal
L}_{N}\phi(x{\,|\,}t)&=\alpha\theta\bigg{(}\frac{(1-x+\delta x)\phi(x-\delta
x{\,|\,}t)-(1-x)\phi(x{\,|\,}t)}{\delta x}\bigg{)}\\\
&\qquad+\beta\theta\bigg{(}\frac{(x+\delta x)\phi(x+\delta
x{\,|\,}t)-x\phi(x{\,|\,}t)}{\delta x}\bigg{)}\\\
&\qquad+\bigg{(}\frac{(x-\delta x)(1-x+\delta x)\phi(x-\delta
x{\,|\,}t)}{\delta x^{2}}\\\ &\qquad+\frac{(x+\delta x)(1-x-\delta
x)\phi(x+\delta x{\,|\,}t)}{\delta x^{2}}-\frac{2x(1-x)\phi(x{\,|\,}t)}{\delta
x^{2}}\bigg{)}.\\\ \end{split}$ (5)

Taking $\lim_{N\to\infty}{\cal L}_{N}\phi(x{\,|\,}t)$ and recognizing that
each term defines a first or second derivative with respect to $x$ immediately
recovers the infinitesimal operator [9, formula 2.11]

${\cal L}\phi(x{\,|\,}t)=-\frac{\partial}{\partial
x}\theta(\alpha-x)\phi(x{\,|\,}t)+\frac{\partial^{2}}{\partial
x^{2}}x(1-x)\phi(x{\,|\,}t)\,$ (6)

of the Kolmogorov forward (Fokker-Planck) diffusion equation for a general
biallelic mutation-drift model:

$\frac{\partial}{\partial t}\phi(x{\,|\,}t)={\cal L}\phi(x{\,|\,}t).$ (7)

Again, let us take a population sample of size $K$, $0\leq K<\infty$, at the
current time $t=0$. The likelihood of a sample of size $K$ with $y$ alleles of
the focal type is binomial:

$\Pr(y{\,|\,}K,x=\tfrac{i}{N},t=0)=\binom{K}{y}x^{y}(1-x)^{K-y}\,.$ (8)

The trajectory of the sample allele proportions backward in time is described
by the Kolmogorov backward equation:

$-\frac{\partial}{\partial t}\Pr(y{\,|\,}K,x,t)={\cal
L}^{{}^{\prime}}\Pr(y{\,|\,}K,x,t),$ (9)

with the backwards operator derived from Eq. (4) [3, formula 80]:

${\cal
L}^{{}^{\prime}}\Pr(y{\,|\,}K,x,t)=\theta(\alpha-x)\frac{\partial}{\partial
x}\Pr(y{\,|\,}K,x,t)+x(1-x)\frac{\partial^{2}}{\partial
x^{2}}\Pr(y{\,|\,}K,x,t)\,.$ (10)

The discrete probability distribution $\Pr(y{\,|\,}K,x,t)$ in Eq. (9) is
interpreted as the probability of obtaining the extant sample configuration
$(y,K)$ conditional on an allele proportion $x$ at past times $t$ [3]. Note
that the negative sign on the left side of the backward diffusion equation Eq.
(9) (opposite to [10]), ensures compatibility of the direction of time between
the forward and backward Kolmogorov equations [38].

####  2.2.1 Modified Jacobi polynomials

Obtaining explicit analytical representations of the transition density at
different time steps is a non-trivial problem. We will use the method of Song
and Steinrücken [25] as adapted by Bergman et al. [3] and introduce the
(modified) Jacobi polynomials (compare also formula 22.3.2 in [1]):

$R_{n}^{(\alpha,\theta)}(x)=\sum_{l=0}^{n}(-1)^{l}\frac{\Gamma(n-1+l+\theta)\Gamma(n+\alpha\theta)}{\Gamma(n-1+\theta)\Gamma(l+\alpha\theta)l!(n-l)!}x^{l}\,,$
(11)

where $n$, with $0\leq n\leq\infty$, is the order of the polynomial. Note that
any regular polynomial of order $n$ can be represented as a weighted sum of
the above modified Jacobi polynomials.

The modified Jacobi polynomials for any orders $m,n$, fulfill the following
orthogonality relationship with respect to the weight function
$w(x,\alpha,\theta)=x^{\alpha\theta-1}(1-x)^{\beta\theta-1}$:

$\int_{0}^{1}R_{n}^{(\alpha,\theta)}(x)R_{m}^{(\alpha,\theta)}(x)\,w(x,\alpha,\theta)\,dx=\delta_{n,m}\Delta_{n}^{(\alpha,\theta)}\,,$
(12)

where $\delta_{n,m}$ is Kronecker’s delta, and

$\Delta_{n}^{(\alpha,\theta)}=\frac{\Gamma(n+\alpha\theta)\Gamma(n+\beta\theta)}{(2n+\theta-1)\Gamma(n+\theta-1)\Gamma(n+1)}\,$
(13)

is the proportionality constant.

The forward and backward Kolmogorov operators (Eqs. (7) and (9), respectively)
can be conveniently decomposed with the modified Jacobi polynomials as
eigenfunctions (for details see [25, 3]). The forward operator becomes:

$-\lambda_{n}w(x,\alpha,\theta)R_{n}^{(\alpha,\theta)}(x)={\cal
L}w(x,\alpha,\theta)R_{n}^{(\alpha,\theta)}(x)\,,$ (14)

and the backward operator:

$\lambda_{n}R_{n}^{(\alpha,\theta)}(x)={\cal
L}^{{}^{\prime}}R_{n}^{(\alpha,\theta)}(x)\,,$ (15)

with corresponding eigenvalues

$\lambda_{n}=n(n+\theta-1)\,.$ (16)

###  2.3 Forward-backward algorithm

Forward-backward algorithms are dynamic programming techniques that enable the
efficient calculation of model states at any time from a sequence of
observations. In our case, we have information from a population sample at the
current time and aim to infer the distribution of past population allele
frequencies. As shown in Bergman et al. [3], the forward-backward algorithm
classically used for hidden Markov models [23, 31] can be readily applied to
the discrete decoupled Moran model of Sec. (2.1): the population allele
proportions are considered ‘hidden’ states and the sample allele
configurations ‘emitted’. In Sec. (2.3.2), we again follow Bergman et al. [3]
in establishing that a forward-backward algorithm can be constructed for the
diffusion model of Sec. (2.2) if the diffusion model is represented using the
modified Jacobi polynomials.

####  2.3.1 Discrete Moran model

We here reproduce the outline of the forward-backward algorithm for the
discrete decoupled Moran model [3].

##### Forward in time

We start at time $s=S$ and assume that the allele proportions of the ancestral
population are distributed according to an arbitrary distribution
$\boldsymbol{\rho}(x)$. Recall that the beta distribution describes the allele
proportions of a bi-allelic, general mutation Moran model in equilibrium [36].
Multiplying with the binomial sampling likelihood and integrating over allele
proportions results a beta-binomial compound distribution as marginal
likelihood. For the biallelic Moran model, the following beta-binomial
distribution with arbitrary mutation-drift parameters is therefore the
standard prior $\boldsymbol{\rho}(x)$:

$\begin{split}\Pr(Nx=i{\,|\,}N,\alpha,\theta)\\\
&=\binom{N}{i}\,\frac{\Gamma(\theta)}{\Gamma(\alpha\theta)\Gamma(\beta\theta)}\frac{\Gamma(i+\alpha\theta)\Gamma(N-i+\beta\theta)}{\Gamma(N+\theta)}\,.\end{split}$
(17)

Any starting distribution can be represented as a row vector of probabilities
$\mathbf{f}_{S}=\boldsymbol{\rho}(x)$, where each entry corresponds to the
probability of allele proportions being $0$, $1$,…,$N$. The probabilities of
allele proportions at any time between $s=S$ and $s=0$ given our prior
distribution, $\mathbf{f}_{s}=\Pr(Nx_{s}{\,|\,}\boldsymbol{\rho})$, can be
determined via:

$\mathbf{f}_{s+1}=\mathbf{f}_{s}\mathbf{T}_{i,j}\quad(S\leq s<0),$ (18)

where ${T}_{i,j}$ is the transition matrix defined in Eq. (1).

At $s=0$, the entries of row vector $\mathbf{b}_{0_{i}}$ are given by the
hypergeometric sampling scheme from Eq. (3) for each possible extant focal
allele proportion between $0$ and $N$. The marginal likelihood of the observed
sample allele frequency is then:

$\begin{split}\Pr(y{\,|\,}K,x,\boldsymbol{\rho})&=\mathbf{f}_{0}\mathbf{b}_{0}^{\prime}\\\
&=\mathbf{f}_{S}\mathbf{T}_{i,j}^{|S|}\mathbf{b}_{0}^{\prime}\,.\end{split}$
(19)

##### Backward in time

The same marginal likelihood may be obtained by recursing backwards from our
sampling step at $s=0$ with initial probabilities $\mathbf{b}_{0_{i}}$ to the
ancestral population state at $s=S$. Define entries of the row vector
$\mathbf{b}_{s,i}=\Pr(y{\,|\,}K,N,Nx_{s}=i)$—they can be interpreted as the
probability of the data given the population allele proportion at time $s$. We
can recurse back in time by:

$\begin{split}\mathbf{b}_{s}^{\prime}=\mathbf{T}_{i,j}\mathbf{b}_{s+1}^{\prime}\quad(0\geq
s>S)\,,\end{split}$ (20)

At $s=S$, we again obtain the marginal likelihood Eq. (19).

##### Joint and conditional probabilities

At any time $s$, the joint probability of the population allele proportion
$x_{s}=\tfrac{i}{N}$, and the sample allele frequency $y$ conditional on the
starting distribution $\boldsymbol{\rho}$ is:

$\Pr(x_{s}=\tfrac{i}{N},y{\,|\,}\boldsymbol{\rho})=(\mathbf{f}_{s})_{i}(\mathbf{b}_{s})_{i}\,.$
(21)

Furthermore, the probability of the population allele proportions $x_{s}=i/N$
conditional on both the sample allele frequency and the starting distribution
is:

$\Pr(x_{s}=\tfrac{i}{N}{\,|\,}y,\boldsymbol{\rho})=\frac{(\mathbf{f}_{s})_{i}(\mathbf{b}_{s})_{i}}{\mathbf{f}_{s}\mathbf{b}_{s}^{\prime}}\,.$
(22)

##### Summary

The forward-backward algorithm with the Moran model conforms to the canonical
situation [23]: At each time point, the population allele proportions comprise
$N$ hidden states and the transition matrix $\mathbf{T}_{i,j}$ is of dimension
$N\times N$. Conditional on an observed sample at the current time $s=0$ and a
prior distribution on the ancestral allele configuration, the distribution of
past and current population allele proportions can be determined. In
population genetics, population demographic events are usually modeled to
occur at a specific time in the past. Note that driving changes in mutation
parameters can be incorporated into this approach by assuming different
parameters for the prior distribution at $s=S$ than for the transition rate
matrix and therefore the times $S<s\leq 0$. Further changes in population
demography may be modeled by time-dependent transition matrices.

####  2.3.2 Diffusion Model

We now adapt the forward-backward algorithm to the diffusion model described
in Sec. (2.2) using the modified Jacobi polynomials from Sec. (2.2.1).

##### Forward in time

Suppose that at time $t=S$, with $S\leq 0$, the distribution of ancestral
allele proportions is given by the arbitrary distribution
$\boldsymbol{\rho}(x)=\phi(x{\,|\,}t=S)$. This distribution is often assumed
to be the beta equilibrium distribution for the bi-allelic Moran model [36]:

$\boldsymbol{\rho}(x)=\frac{\Gamma(\theta)}{\Gamma(\alpha\theta)\Gamma(\beta\theta)}\,x^{\alpha\theta-1}(1-x)^{\beta\theta-1}\,.$
(23)

(Note that the weight function of the modified Jacobi polynomials is
proportional to this beta distribution.) The allele proportions further
forward in time, $\phi(x{\,|\,}t,\boldsymbol{\rho})$, are determined by the
forward diffusion equation Eq. (7). The solution to the forward equation can
be represented using modified Jacobi polynomials; the ancestral allele
proportion distribution is first expanded to:

$\boldsymbol{\rho}(x)=\sum_{n=0}^{\infty}\rho_{n}^{(\alpha,\theta)}R_{n}^{(\alpha,\theta)}(x),$
(24)

where the $\rho_{n}^{(\alpha,\theta)}$ are a possibly infinite number of
Jacobi coefficients that depend on $\boldsymbol{\rho}(x)$. More explicitly:

$\rho_{n}^{(\alpha,\theta)}=\frac{1}{\Delta_{n}^{(\alpha,\theta)}}\int_{0}^{1}w(x,\alpha,\theta)R_{n}^{(\alpha,\theta)}(x)\boldsymbol{\rho}(x)\,dx\,.$
(25)

We then incorporate temporal dynamics and obtain the full solution:

$\phi(x{\,|\,}t,\boldsymbol{\rho})=w(x,\alpha,\theta)\sum_{n=0}^{\infty}\rho_{n}^{(\alpha,\theta)}R_{n}^{(\alpha,\theta)}(x)e^{\lambda_{n}(S-t)}\,.$
(26)

##### Backward in time

Backward in time, we again start with our sample: The binomial likelihood of
the sampled allele proportions at time $t=0$ in Eq. (8) is expressed as a
regular polynomial up to order $K$ with coefficients

$a_{j=y+i}(K,y)=(-1)^{i}\binom{K}{y}\binom{K-y}{i}$

for $0\leq i\leq K-y$, and zero otherwise. Let $\mathbf{a}(K,y)$ be the vector
of coefficients $a_{j}(K,y)$ and $\mathbf{R}^{(\alpha,\theta)}$ be the matrix
of coefficients $R_{n}^{(\alpha,\theta)}(x)$. Note that this matrix is lower
triangular. Then the binomial distribution can be uniquely expanded into
Jacobi polynomials via the following linear algebraic equation:

$\mathbf{d}^{(\alpha,\theta)}(K,y)=\mathbf{a}(K,y)\mathbf{R}^{(\alpha,\theta)}$
(27)

Note that the triangular structure of $\mathbf{R}^{(\alpha,\theta)}$ obviates
matrix inversion. Now the binomial sampling distribution from Eq. 8 can be
rewritten:

$\Pr(y{\,|\,}K,x,\alpha,\theta,t=0)=\binom{K}{y}x^{y}(1-x)^{K-y}=\sum_{n=0}^{K}d_{n}^{(\alpha,\theta)}(K,y)R_{n}^{(\alpha,\theta)}(x)\,.$
(28)

Further back, at times $t$ ($S\leq t\leq 0$), the distribution of sample
proportions is given by:

$\Pr(y{\,|\,}K,x,\alpha,\theta,t)=\sum_{n=0}^{K}d_{n}^{(\alpha,\theta)}(K,y)R_{n}^{(\alpha,\theta)}(x)e^{\lambda_{n}t}\,.$
(29)

Using the orthogonality of the Jacobi polynomials (Eq. 12), the continuous
marginal likelihood becomes:

$\begin{split}\Pr(y{\,|\,}K,\alpha,\theta,S,\boldsymbol{\rho})&=\int_{0}^{1}\phi(x{\,|\,}t,\boldsymbol{\rho})\Pr(y{\,|\,}K,x,\alpha,\theta,t=0)\,dx\\\
&=\sum_{n=0}^{K}\rho_{n}^{(\alpha,\theta)}d_{n}^{(\alpha,\theta)}(K,y)\Delta_{n}^{(\alpha,\theta)}e^{\lambda_{n}S}\,.\end{split}$
(30)

Note that the expansion $\rho_{n}^{(\alpha,\theta)}$ may be infinite; however,
calculation of the marginal likelihood only requires expansion to the order of
the sample size $K$.

As briefly noted in the summary of Sec. (2.3.1), it is often convenient to be
able to account for population demographics. We will consider this possibility
for the diffusion approach: Assume the mutation parameters change from
$\alpha$ to $\alpha^{*}$ and from $\theta$ to $\theta^{*}$ at time $t=S$. The
coefficients of the new polynomial expansions can be obtained by linear
transformation. More explicitly, consider a simple model where:

  * 1.

i) at the present time $t=0$, the sample allele configuration $(y,K)$ is
given;

  * 2.

ii) between the times $S\leq t\leq 0$, the population genetic parameters
$\alpha$ and $\theta$ remain constant; and

  * 3.

iii) at time $t=S$ in the past the population allele proportion is beta
distributed according to:

$\boldsymbol{\rho}(x)=\frac{\Gamma(\theta^{*})}{\Gamma(\alpha^{*}\theta^{*})\Gamma(\beta^{*}\theta^{*})}\,x^{\alpha^{*}\theta^{*}-1}(1-x)^{\beta^{*}\theta^{*}-1}$
(31)

In fact, the allele proportion is beta distributed as above between
$-\infty\leq t\leq S$ and can be expressed as the series
$\rho_{n}^{(\alpha^{*},\theta^{*})}R_{n}^{(\alpha^{*},\theta^{*})}(x)$.

The continuous marginal distribution then becomes:

$\begin{split}\Pr(y{\,|\,}K,\alpha,\theta,S,\boldsymbol{\rho})&=\int_{0}^{1}\phi(x{\,|\,}t,\boldsymbol{\rho})\Pr(y{\,|\,}K,x,\alpha,\theta,t=0)\,dx\\\
&=\sum_{n=0}^{K}\rho_{n}^{(\alpha^{*},\theta^{*})}d_{n}^{(\alpha,\theta)}(K,y)\Delta_{n}^{(\alpha,\theta)}e^{\lambda_{n}S}\,.\end{split}$
(32)

##### Joint and conditional probabilities

At any time $t$, the joint probability of the population allele proportion
$x_{t}=\tfrac{i}{N}$, and the number of focal alleles in the sample $y$
conditional on the starting distribution $\boldsymbol{\rho}$ can be
determined:

$\Pr(x_{t}=\tfrac{i}{N},y{\,|\,}\boldsymbol{\rho})=\phi(x{\,|\,}t,\boldsymbol{\rho})\Pr(y{\,|\,}K,x,t=0).$

The probability of the population allele frequencies $x_{t}=\tfrac{i}{N}$
conditional on the both the sample allele proportions and the starting
distribution is:

$\Pr(x_{t}=\tfrac{i}{N}{\,|\,}y,\boldsymbol{\rho})=\frac{\Pr(x_{t}=\tfrac{i}{N},y{\,|\,}\boldsymbol{\rho})}{\Pr(y{\,|\,}K,\alpha,\theta,S,\boldsymbol{\rho})}\,.$

We will not detail these equations here.

##### Summary

The forward-backward algorithm with the continuous diffusion model represented
using modified orthogonal Jacobi polynomials deviates from the canonical
situation: A transition kernel for population allele proportions is employed,
which is expanded into an infinite-dimensional system of eigenfunctions and
corresponding eigenvalues. For representing sample allele frequencies,
however, only an expansion of the order of the sample size is needed. Indeed
most problems only require a polynomial expansion up to the order of the
sample size and the temporal system required is diagonal and thus extremely
simple. Furthermore, a change in the mutation parameters can also be
incorporated.

##  3 Particle models and orthogonal polynomials

In this section, we introduce a novel forward-backward algorithm that
harnesses together three components:

  * 1.

i) An approach based on orthogonal polynomials to describe the evolution of
the population allele frequencies forward in time (Sec. 3.8),

  * 2.

ii) a so-called particle model that yields the conditional probabilities of
the proportion of the focal allele at any point in the history of the sample
by running traditionally backwards-looking coalescent arguments not only
backward but also forward in time (Sec. 3.2), (Sec. 3.3); this is augmented by

  * 3.

iii) backward-in-time temporal dynamics accounting for the effect of changing
mutation parameters on the sample sizes (Sec. 3.5) to yield probabilities of
all past sample configurations.

In total, we will arrive at the following joint probability:

$\Pr(\upsilon,\kappa{\,|\,}K,y,\alpha,\theta,t)\,.$ (33)

This is the probability of seeing $\upsilon$ focal alleles in a sample of size
$\kappa$ at any time $t<0$ in the past, conditional on an extant sample of
size $K$ with $y$ focal alleles (and the underlying mutation-drift parameters
$\theta$ and $\alpha$) at time $t=0$.

We now begin by motivating the particle model.

##### Sample genealogy backward in time: coalescent

In the (decoupled) Moran model of population size $N$, time is conventionally
scaled so that the genealogies of all individuals in the population conform to
the Kingman coalescent. Furthermore, the genealogy of any sample of size $K<N$
is embedded within the (decoupled) Moran model: The sample probabilities and
transition rates of the coalescent remain unaffected by a change in sample
size. In other words, the classic (decoupled) Moran model is sample
consistent, and is dual to the Kingman coalescent in the sense that the
expected population allele frequencies are identical between the two models
[9, chapt. 2.8].

In Fig. (2), the history of a sample of size $K=6$ with $y=4$ focal alleles is
depicted. Both mutation and coalescent events are uninformative regarding the
immediately preceding allelic state. Starting at the present time $t=0$ and
looking back, the size of the sample $\kappa$ with $K\geq\kappa\geq 0$ is
reduced in discrete stages by these events [compare 11, for a similar argument
in the context of the infinite sites model]. With the usual scaling of time,
the rate of a coalescent is $\kappa(\kappa-1)$ and that of a mutation
$\kappa\theta$, hence the total rate of reduction events is
$\kappa(\kappa-1+\theta)$. At any time $t<0$ the number of alleles $\kappa$
remaining from the original sample of size $K$ at $t=0$ is a random variable,
as is the number of focal alleles remaining from the $y$ focal alleles at
$t=0$, these we denote $\upsilon$. The probability of a reduction in the
number of focal alleles from $\upsilon$ to $\upsilon-1$ in each backward time
step is proportional to $\upsilon(\upsilon-1+\alpha\theta)$ (the contribution
of the coalescent event is $\upsilon(\upsilon-1)$, that of the mutation event
$\alpha\theta\upsilon$). The analogous reduction probabilities for the non-
focal alleles are $(\kappa-\upsilon)(\kappa-\upsilon-1+\beta\theta)$. Note
that these reduction probabilities are for ordered events. However, we will
treat allele ’labels’ as interchangeable in each time step, and therefore use
unordered reduction probabilities. We obtain these by dividing each of the
previous ordered reduction probabilities by the number of potential alleles
selected for a reduction event.

Let us now reverse the direction of time: Consider the indicator variable
$z_{\kappa+1}$, which is one if the $\kappa+1$st allele is of the focal type
and zero otherwise. In each forward step, the probability of going from
$\kappa$ to $\kappa+1$ unordered focal alleles is:

$\begin{split}{\Pr}(z_{\kappa+1}{\,|\,}\kappa,\upsilon,\alpha,\theta)&=\frac{(\upsilon+\alpha\theta)^{1-z_{\kappa}}(\kappa-\upsilon+\beta\theta)^{z_{\kappa}}}{\kappa+\theta}\\\
&=\frac{\Gamma(\kappa+\theta)}{\Gamma(\kappa+1+\theta)}\frac{\Gamma(\upsilon+z_{\kappa}+\alpha\theta)}{\Gamma(\upsilon+\alpha\theta)}\frac{\Gamma(\kappa+1-\upsilon-
z_{\kappa}+\beta\theta)}{\Gamma(\kappa-\upsilon+\beta\theta)}\,.\end{split}$
(34)

Figure 2: Schematic plot of a Moran model and embedded sample genealogy: The
allelic type of the extant individuals is given at the top of the plot;
backwards in time the arrows indicate birth/death events and the ’x’ are
mutation events. Note that the position and numbering of individuals on the
lower x-axis is arbitrary. The bold lines and arrows correspond to the
coalescent and mutation history of the sample respectively.

##### Sample genealogy forward in time: urn models

Let us briefly view the probability of sample configurations running forward
in time as an urn model similar to the Polya- or Hoppe-urns or the urn-models
in [26]. Hereto, we introduce a new time index $t_{\kappa-1}$ that runs from
$t_{0}=0$, a time in the past at which only one individual from the sample is
present in the population, to the present time $t_{K}$. The urn model is then
given by the following algorithm:

  * 1.

Initiation. Start with a sample of size $\kappa=1$ at time $t_{0}=0$. The
probability of the focal allele is set to $\alpha$.

  * 2.

Recursion. Add a period of rate $\kappa(\kappa-1+\theta)$ to $t_{\kappa-1}$
obtain $t_{\kappa}$. Increase the sample size from $\kappa$ to $\kappa+1$ by a
focal allele with probability

$p_{\upsilon+1}=\frac{\upsilon+\alpha\theta}{\kappa+\theta}\,,$ (35)

or increase the sample size by a non-focal allele with probability
$1-p_{\upsilon+1}$. Together this corresponds to the probability of the
indicator variable in Eq. (34).

  * 3.

Stop. When the sample size $\kappa=K$ is reached, add a period of rate
$K(K-1+\theta)$ to $t_{K-1}$ obtain $t_{K}$.

In order to obtain our usual time index with an extant time of zero, $t_{K}$
must be subtracted from all $t_{\kappa-1}$. Thus running the coalescent
process forwards in time as in Eq. (34) yields an urn model.

###  3.1 Particle model: algorithms

We will now demonstrate how running coalescent arguments for the genealogy of
an unordered sample from a decoupled Moran model both forward and backward in
time can be used to determine sample allele frequencies for every past sample
size (Sec. 3.2), (Sec. 3.3). Note that the past sample sizes are initially
assumed given because they are determined by the constant rate of coalescence.
Clearly, these sample sizes depend on mutation-drift parameters that are not
necessarily constant across time. Note, however, that the number of focal
alleles is conditionally independent of time given the (past) sample sizes
$\kappa$ and the mutation-drift parameters. Hence we first treat the time
independent dynamics, and then augment this with differential backward
equations (Sec. 3.5) to account for changing mutation-drift parameters.

###  3.2 Particle model: algorithm for forward probabilities

For every $0\leq\kappa\leq K$, let $f(.,\kappa)$ denote a forward row vector
of length $\kappa+1$. Each entry $f(\upsilon,\kappa)$ can be interpreted as a
probability $\Pr(\upsilon{\,|\,}\kappa,\alpha,\theta)$, where
$0\leq\upsilon\leq\kappa$.

  * 1.

Initiation. Start with a sample of size $\kappa=0$. Trivially, the frequency
of the focal allele is $\upsilon=0$, so set $f(\upsilon=0,\kappa=0)=1$.

  * 2.

Recursion. Move from sample size $\kappa$ to $\kappa+1$ by calculating every
entry of $f(.,\kappa+1)$:

$\begin{split}f(\upsilon,\kappa+1)&=\frac{\kappa-\upsilon-1+\beta\theta}{\kappa+\theta}\,f(\upsilon,\kappa-1)+\frac{\upsilon-1+\alpha\theta}{\kappa+\theta}\,f(\upsilon-1,\kappa-1)\,.\end{split}$
(36)

  * 3.

Stop. Stop when the sample size $\kappa=K$ is reached.

With time independent mutation-drift parameters, $f(\upsilon,\kappa)$ is a
beta-binomial compound distribution:

$\begin{split}f(\upsilon,\kappa)&=\Pr(\upsilon{\,|\,}\kappa,\alpha,\theta)\\\
&=\binom{\kappa}{\upsilon}\,\frac{\Gamma(\theta)}{\Gamma(\alpha\theta)\Gamma(\beta\theta)}\frac{\Gamma(\upsilon+\alpha\theta)\Gamma(\kappa-\upsilon+\beta\theta)}{\Gamma(\kappa+\theta)}\,.\end{split}$
(37)

###  3.3 Particle model: algorithm for backward probabilities

For every $0\leq\kappa\leq K$, introduce the backward row vector $b(.,\kappa)$
of length $\kappa+1$. Every entry $b(\upsilon,\kappa)$ corresponds to
${\Pr}(y{\,|\,}K,\kappa,\upsilon,\alpha,\theta)$ for
$0\leq\upsilon\leq\kappa$, which is the probability of the observed number of
focal alleles given past allele configurations.

  * 1.

Initiation. Start with a sample of size $K$ with $y$ focal alleles. Set
$\upsilon=y$ and $\kappa=K$, and therefore $b(\upsilon=y,\kappa=K)=1$.

  * 2.

Recursion. Move from $\kappa+1$ to $\kappa$ by calculating the entries of
$b(.,\kappa)$:

$\begin{split}b(\upsilon,\kappa)&=\frac{\kappa-\upsilon+\beta\theta}{\kappa+\theta}\,b(\upsilon,\kappa+1)+\frac{\upsilon+\alpha\theta}{\kappa+\theta}\,b(\upsilon+1,\kappa+1)\,.\end{split}$
(38)

  * 3.

Stop. End the recursion, when the sample size $\kappa=0$ is reached. Note that
$b(\upsilon=0,\kappa=0)={\Pr}(y{\,|\,}K,\upsilon=0,\kappa=0,\alpha,\theta)$
corresponds to the likelihood ${\Pr}(y{\,|\,}K,\alpha,\theta)$.

With time independent mutation-drift parameters, $b(\upsilon,\kappa)$ is again
a beta-binomial compound distribution:

$\begin{split}b(\upsilon,\kappa)&=\Pr(y{\,|\,}K,\kappa,\upsilon,\alpha,\theta)\\\
&=\binom{K-\kappa}{y-\upsilon}\,\frac{\Gamma(\kappa+\theta)}{\Gamma(\upsilon+\alpha\theta)\Gamma(\kappa-\upsilon+\beta\theta)}\frac{\Gamma(y+\alpha\theta)\Gamma(K-y+\beta\theta)}{\Gamma(K+\theta)}\,.\end{split}$
(39)

###  3.4 Particle model: time independent conditional probabilities of focal
alleles

We can now combine the forward and backward algorithms and calculate the
likelihood of seeing $y$ focal alleles in a sample of size $K$ given time
independent mutation-drift parameters (note the use of Vandermonde’s identity
to obtain the final result):

$\begin{split}\Pr(y{\,|\,}K,\alpha,\theta)&=\sum_{\upsilon=0}^{y}f(\upsilon,\kappa)b(\upsilon,\kappa)\\\
&=\sum_{\upsilon=0}^{y}\binom{\kappa}{\upsilon}\,\frac{\Gamma(\theta)}{\Gamma(\alpha\theta)\Gamma(\beta\theta)}\frac{\Gamma(\upsilon+\alpha\theta)\Gamma(\kappa-\upsilon+\beta\theta)}{\Gamma(\kappa+\theta)}\\\
&\quad\times\binom{K-\kappa}{y-\upsilon}\,\frac{\Gamma(\kappa+\theta)}{\Gamma(\upsilon+\alpha\theta)\Gamma(\kappa-\upsilon+\beta\theta)}\frac{\Gamma(y+\alpha\theta)\Gamma(K-y+\beta\theta)}{\Gamma(K+\theta)}\\\
&=\sum_{\upsilon=0}^{y}\binom{\kappa}{\upsilon}\binom{K-\kappa}{y-\upsilon}\frac{\Gamma(\theta)}{\Gamma(\alpha\theta)\Gamma(\beta\theta)}\frac{\Gamma(y+\alpha\theta)\Gamma(K-y+\beta\theta)}{\Gamma(K+\theta)}\\\
&=\binom{K}{y}\frac{\Gamma(\theta)}{\Gamma(\alpha\theta)\Gamma(\beta\theta)}\frac{\Gamma(y+\alpha\theta)\Gamma(K-y+\beta\theta)}{\Gamma(K+\theta)}\,.\end{split}$
(40)

Similarly, the probability of seeing $\upsilon$ focal alleles in a given past
sample of size $\kappa$ assuming a current sample of size $K$ containing $y$
focal alleles and time independent mutation-drift parameters is:

$\Pr(\upsilon{\,|\,}K,y,\kappa,\alpha,\theta)=\frac{f(\upsilon,\kappa)b(\upsilon,\kappa)}{\sum_{\upsilon=0}^{\kappa}f(\upsilon,\kappa)b(\upsilon,\kappa)}\,.$
(41)

We provide a visual example of the likelihood of the focal locus traced
backward in time in Fig. (3).

Figure 3: Consider a sample of size $K=10$ with $y=5$ focal alleles, i.e., a
starting configuration of $(5,10)$. Assume a mutation bias of $\alpha=0.3$
towards the focal allele and an overall scaled mutation rate of $\theta=0.1$.
Here, we show the conditional probabilities of all possible particle
configurations between $(5,10)$ and $(0,0)$ as per Eq. (41). This is an
application of the forward- backward particle model in Secs. (3.2) and (3.3).

###  3.5 Particle model: temporal dynamics

So far, we have always assumed the past sample sizes $\kappa$ as given.
However, they can easily be modelled as a pure death process dependent on
mutation and drift: Both events simply reduce the sample size by one.

The following system of differential equations describes the evolution of the
sample size $\kappa$ backward in time:

$\begin{split}-\frac{d}{dt}\Pr(\kappa=K{\,|\,}K,\theta,t)&=\lambda_{K}\,\Pr(\kappa=K{\,|\,}K,\theta,t)\,,\quad\text{and}\\\
-\frac{d}{dt}\Pr(\kappa{\,|\,}K,\theta,t)&=-\lambda_{\kappa+1}\,\Pr(\kappa+1{\,|\,}K,\theta,t)\\\
&\qquad+\lambda_{\kappa}\,\Pr(\kappa{\,|\,}K,\theta,t)\,,\quad\text{for
$K>\kappa\geq 0$.}\end{split}$ (42)

with corresponding eigenvalues $\lambda_{\kappa}=\kappa(\kappa-1+\theta)$ (as
in Eq. (16)). The starting condition is: $\Pr(\kappa=K{\,|\,}K,\theta,t=0)=1$,
and $\Pr(K>\kappa\geq 0{\,|\,}K,\theta,t=0)=0$. Then the solution of the above
system of equations is [27, chapter 6, Eq. 2.2]

$\begin{split}\Pr(\kappa=K{\,|\,}K,\theta,t)&=e^{\lambda_{K}t}\,,\\\
\Pr(\kappa{\,|\,}K,\theta,t)&=\sum_{i=\kappa}^{K}c_{i,\kappa}\,e^{\lambda_{i}t}\,\text{,
for $K-1\geq\kappa\geq 0$}\\\ &\text{with
$c_{i,\kappa}=\frac{\prod_{j=\kappa+1}^{K}\lambda_{j}}{\prod_{j=\kappa,j\neq
i}^{K}(\lambda_{j}-\lambda_{i})}$.}\end{split}$ (43)

###  3.6 Particle model: total backward dynamics

Augmenting the time independent backward variables from Eq. (39) with the
temporal dynamics from the previous subsection, i.e., Eq. (43), we obtain the
joint probability of the number of focal alleles in the extant sample and the
past sample sizes:

$\Pr(y,\kappa{\,|\,}K,\upsilon,\alpha,\theta,t)=\Pr(y{\,|\,}K,\kappa,\upsilon,\alpha,\theta)\Pr(\kappa{\,|\,}K,\theta,t)\,.$
(44)

###  3.7 Particle model: joint and marginal probabilities

The aim of the forward-backward algorithm is to obtain the probability of
ancestral sample configurations $(\upsilon,\kappa)$ at arbitrary times $t$,
conditional on the sample configuration $(y,K)$ at time $t=0$, i.e.,
$\Pr(\upsilon,\kappa{\,|\,}K,y,\dots)$. Multiplying the probability in Eq.
(44), $\Pr(y,\kappa{\,|\,}K,\upsilon,\dots)$, with the probability of the
frequency of the focal allele $\upsilon$ given a past sample size $\kappa$ and
population proportion $x$ gives us the joint probability
$\Pr(y,\upsilon,\kappa{\,|\,}K,x,\dots)$. Recall that the probabilities in Eq.
(44) are for unordered samples, since we treat the allele ’labels’ as
interchangeable in each time step within the sample genealogy. Obtaining the
number of focal alleles $\upsilon$ in an ancestral sample from a given
population allele proportion $x$, however, requires ordered sampling (instead
of binomial i.e., unordered sampling):

$\Pr(\upsilon{\,|\,}\kappa,x)^{o}=x^{\upsilon}(1-x)^{\kappa-\upsilon}\,,$ (45)

to obtain:

$\Pr(y,\upsilon,\kappa{\,|\,}K,x,\alpha,\theta,t)=\Pr(y{\,|\,}K,\kappa,\upsilon,\alpha,\theta)\Pr(\upsilon{\,|\,}\kappa,x)^{o}\Pr(\kappa{\,|\,}K,\theta,t)\,.$
(46)

The marginal distribution of the number of focal alleles $y$ in the extant
sample of size $K$ is then determined by:

$\Pr(y{\,|\,}K,x,\alpha,\theta,t)=\sum_{\kappa=0}^{K}\sum_{\upsilon=\max(\kappa+y-K,0)}^{\min(\kappa,y)}\Pr(y,\upsilon,\kappa{\,|\,}K,x,\alpha,\theta,t)\,,$
(47)

for every $t<0$. This is a sum over all possible $(K+1)(y+1)$ allele
configurations. For $t\to-\infty$, the sum becomes dominated by the
probability for the case $\kappa=0$ since the likelihood in Eq. (45) is
trivially one for $\kappa=0$ irrespective of $x$. Hence, this can be used to
obtain the likelihood:

$\Pr(y{\,|\,}K,x,\theta,\alpha,t\to\infty)=\Pr(y{\,|\,}K,\theta,\alpha)\,.$
(48)

####  3.7.1 Relationship between particle probabilities and Jacobi polynomials

Recall that we previously obtained an expression for Eq. (47) in terms of
modified Jacobi polynomials (Eq. 29), which we reproduce here:

$\Pr(y{\,|\,}K,x,\alpha,\theta,t)=\sum_{n=0}^{K}e^{\lambda_{n}t}d_{n}(K,y)^{(\alpha,\theta)}R_{n}^{(\alpha,\theta)}(x)\,.$
(49)

Clearly, these two representations must be equivalent: The order of expansion
of the modified Jacobi polynomials is generally equivalent to the sample size,
so we can set $\kappa$ to $n$ in Eq. (47). Now, for each past sample size
$\kappa$ (or now $n$), the terms with the same temporal component
$e^{\lambda_{n}t}$ and the same power of $x$ from Eq. (47)(left below) and Eq.
(49)(right below) can be equated:

$\begin{split}e^{\lambda_{n}t}x^{n}c_{n,n}\sum_{\upsilon=\max(n+y-K,0)}^{\min(\kappa,y)}(-1)^{n-\upsilon}b(\upsilon,n)&=e^{\lambda_{n}t}x^{n}d_{n}(K,y)^{(\alpha,\theta)}r_{n,n}\\\
c_{n,n}\sum_{\upsilon=\max(n+y-K,0)}^{\min(\kappa,y)}(-1)^{n-\upsilon}b(\upsilon,n)&=d_{n}(K,y)^{(\alpha,\theta)}r_{n,n}\,\end{split}$
(50)

where the $R_{n}^{(\alpha,\theta)}(x)$ are written as the term $x^{n}$
multiplied by the corresponding coefficient $r_{n,l}$, and $c_{n,n}$ are the
$n$th coefficients solving the system of temporal differential equations in
Eq. (42). For $y=0$ or $y=K$, the above formula simplifies as there only a
single term remains in the summation.

###  3.8 Past allele configurations

In order to determine the distribution of allele configurations at any given
point in time conditional on a current sample, we must combine the total
backwards dynamics of a sample given in Eq. (44) with repeated sampling from
the general population forward in time. Note that the latter cannot be
modelled using a particle model: The sample sizes in the particle model - both
current and past - are assumed fixed until augmented backward in time by
temporal dynamics accounting for mutation-drift parameters; importantly, these
temporal components are conditional on the extant sample size. When sampling
from the population forward in time and agnostic to either past or future
sample sizes and configurations, however, the sample size $\kappa$ becomes a
random variable.

The population allele proportion can be modelled forward in time by expanding
the transition density $\phi(x{\,|\,}t,\boldsymbol{\rho})$ into orthogonal
polynomials as in Eq. (26); $\boldsymbol{\rho}$ here represented the beta
equilibrium distribution of the bi-allelic Moran model. Further, the
likelihood of observing $\upsilon$ alleles of the focal type in a sample of
size $\kappa$ conditional on the population allele frequency $x$ can also be
expanded into the orthogonal polynomials:

$\Pr(\upsilon{\,|\,}\kappa,x)=\binom{\kappa}{\upsilon}x^{\upsilon}(1-x)^{\kappa-\upsilon}=\sum_{n=0}^{\kappa}{d}_{n}(\kappa,\upsilon)R_{n}^{(\alpha,\theta)}(x)\,.$
(51)

Then the joint probability of the focal alleles $y$ in the extant sample and
the past particle configurations becomes:

$\begin{split}&\Pr(y,\upsilon,\kappa{\,|\,}K,\alpha,\theta,t,\boldsymbol{\rho})\\\
&\qquad=\Pr(y,\kappa{\,|\,}K,\upsilon,\alpha,\theta,t)\int_{0}^{1}\Pr(\upsilon{\,|\,}\kappa,x)\,\phi(x{\,|\,}t,\boldsymbol{\rho})\,dx\\\
&\qquad=\Pr(\upsilon{\,|\,}K,y,\kappa,\alpha,\theta)\,\Pr(\kappa{\,|\,}K,\theta,t)\\\
&\qquad\qquad\times\sum_{n=0}^{\kappa}\,e^{\lambda_{j}(S-t)}\int_{0}^{1}{d}_{n}(\kappa,\upsilon)R_{n}^{(\alpha,\theta)}(x)\rho_{j}^{(\alpha,\theta)}R_{n}^{(\alpha,\theta)}(x)x^{\alpha\theta-1}(1-x)^{\beta\theta-1}\,dx\\\
&\qquad=\Pr(\upsilon{\,|\,}K,y,\kappa,\alpha,\theta)\,\Pr(\kappa{\,|\,}K,\theta,t)\,\sum_{n=0}^{\kappa}e^{\lambda_{j}(S-t)}{d}_{n}(\kappa,\upsilon)\rho_{n}^{(\alpha,\theta)}\Delta_{n}^{(\alpha,\theta)}\,.\end{split}$
(52)

Note that the first term can be calculated using our forward-backward particle
algorithm in Eq. (41), the second is the solution of the temporal system in
Eq. (43).

Finally, the probability of ancestral allele configurations at any time,
conditional on an extant sample, becomes:

$\begin{split}\Pr(\upsilon,\kappa{\,|\,}K,y,\alpha,\theta,t,\boldsymbol{\rho})&=\frac{\Pr(y,\upsilon,\kappa{\,|\,}K,\alpha,\theta,t,\boldsymbol{\rho})}{\Pr(y{\,|\,}K,\alpha,\theta,t,\boldsymbol{\rho})}\,,\end{split}$
(53)

where the denominator can be determined using either

  * 1.

(i) the particle approach via

$\Pr(y{\,|\,}K,\alpha,\theta,t,\boldsymbol{\rho})=\Pr(y{\,|\,}K,x,\alpha,\theta,t)\int_{0}^{1}\phi(x{\,|\,}t,\boldsymbol{\rho})\,dx\\\
$

with the first term from Eq. (47), or

  * 2.

(ii) the equivalent polynomial approach in Eq. (30).

Changes in mutation parameters over time can also be incorporated in this
forward-backward algorithm: Once again consider the scenario in which the
mutation parameters $\alpha$ and $\theta$ remain constant between the times
$S\leq t\leq 0$, and at time $t=S$ in the past they change to $\alpha^{*}$ and
$\theta^{*}$. Recall that the ancestral allele configuration
$\boldsymbol{\rho}$ is then modelled as beta-binomial with parameters
$\alpha^{*}$ and $\theta^{*}$ instead of $\alpha$ and $\theta$ (Eq. (31)), and
this can be substituted into the equations of this subsection accordingly
(directly in Eq. (52) and (i), and see Eq. (32) for (ii)). An example is shown
in Fig. (4).

##### Summary

The forward-backward algorithm using a combination of the particle model and
orthogonal polynomials deviates considerably from the canonical situation:
Within the forward-backward recursions tracing the sample allele
configurations described by the particle model, each step leads to a change in
the sample size $\kappa$ and hence the dimension of the transition matrix,
which is no longer square. Furthermore, the backwards temporal system
accounting for the effect of changing mutation-drift parameters on the sample
size is not diagonal and thus moderately complex. Calculating the marginal
likelihood of the number of focal alleles in the extant sample involves a
summation over $(K+1)(y+1)$ allele configurations rather than an eigensystem
decomposition as with the orthogonal polynomials—this may sometimes be
convenient though it is generally less efficient. However, the two methods
yield equivalent results, which is of considerable theoretical interest. To
determine the likelihood of any past population allele frequencies, the
particle model is combined with forward population dynamics represented as
orthogonal polynomials.

Figure 4: Let us again consider a sample of size $K=10$ with $y=5$ focal
individuals drawn at the present time $t=0$. We now consider a non-equilibrium
situation where the population size changes at $t=S=-0.5$ in the past; this
change is driven by a single change in the scaled mutation rate: The ancestral
rate $\theta^{*}=0.3$ changes to the current rate $\theta^{*}=0.1$,
$\alpha=0.3$ remains constant throughout. We show the probabilities of the
various possible allele configurations $\upsilon$ and $\kappa$ from Eq. (52);
the size of the dots is proportional to the probability of the configuration.

##  4 Transition probabilities of the boundary mutation-drift Moran model

Recall that in synonymous sites of protein coding genes of eukaryotes, the
expected heterozygosity, which is roughly equal to the scaled mutation rate,
is $\leq 10^{-2}$ [17]. Note that the expected heterozygosity corresponds to
the expected polymorphism in a sample of size two and is thus proportional to
the number of polymorphic sites in the sample; increasing the sample size
leads to a roughly logarithmic increase in the proportion of polymorphic
sites. Population genetic models derived in the limit of small scaled mutation
rates, $\theta\to 0$, often approximate general mutation dynamics with
sufficient accuracy and have clear numeric advantages [33]. A first order
Taylor expansion in $\theta$ of the beta-binomial equilibrium distribution of
a sample from the decoupled bi-allelic Moran model yields the equilibrium
distribution of the proportion of alleles $\mathbf{X}$ at any locus in the so-
called boundary mutation-drift Moran model [30]:

${\Pr}(\mathbf{X}=\frac{i}{N}{\,|\,}N,\alpha,\theta)=\begin{cases}\beta\big{(}1-\alpha\theta
H_{N-1}\big{)}&\text{for $i=0$,}\\\
\alpha\beta\theta\,\frac{N}{i(N-i)}&\text{for $1\leq i\leq(N-1)$,}\\\
\alpha\big{(}1-\beta\theta H_{N-1}\big{)}&\text{for $i=N$}\,,\end{cases}$ (54)

where $H_{N-1}=\sum_{i=1}^{N-1}\frac{i}{N}$ is the harmonic number. In
essence, the boundary mutation-drift Moran model assumes that only a single
mutation segregates in a population at any given time: Mutations arise
exclusively from the monomorphic boundaries rather than from near the
boundaries as in the general mutation model with low mutation rates; the
polymorphic interior is governed by drift (and selection in non-neutral
scenarios [32]). Simulations show that the boundary mutation-drift Moran model
is a good approximation for the general mutation Moran model if the expected
equilibrium heterozygosity $2\alpha\beta\theta\leq 10^{-2}$ [30]. Importantly,
the boundary mutation-drift Moran model simplifies inference considerably
[28].

Several approaches have been taken to derive transition probabilities
compatible with this equilibrium distribution [29, 6, 32]. Burden and
Griffiths [6] start their derivations from a Wright-Fisher diffusion with
$\theta\to 0$, and obtain the probability of the admissible sample
configurations under the constrained mutation rate at any time $t$ by
considering all coalescent sub-trees and their scaled branch lengths. The
resulting formulae are complex and not easily employed in an inferential
framework.

The transition rate matrix of the boundary-mutation Moran model given by Vogl
and Mikula [32] is comparatively tractable. To balance the constraint that
mutations are only allowed from the boundaries, mutation rates are re-scaled
to obtain an equilibrium distribution corresponding to the Taylor series
expansion of the general mutation Moran model. Importantly, this transition
matrix and in particular its eigenvalues are not consistent with varying $N$,
precisely because it is scaled so that mutations segregate at an identical
average rate in equilibrium regardless of the sample size. Therefore the
embedded genealogies of alleles depend on the sample size. For small scaled
mutation rates the deviation of the eigenvalues from those of the general
mutation Moran model (which are consistent for varying $N$) is negligible, as
shown in Fig. (5). For larger values such as $\theta=0.025$, which is close to
the limit of feasibility for the first order $\theta$ approximation, the first
non-zero eigenvalue $\lambda_{1}$, which determines the asymptotic speed of
approach to equilibrium, is increased by approximately $2.5\%$ to $5\%$. This
could become important in a phylogenetic context, as split times between
species may considerably exceed the (effective) population size $N$, which in
turn constrains the effect of drift.

The stationary distribution of the phylogenetic rendition of the boundary
mutation-drift Moran model has traditionally been scaled to maintain the
proportion of monomorphic sites across different sample sizes [24] (but note
that one version can be reparametrized to give the other), and the
corresponding transition matrix [5] is therefore similarly inconsistent.
(Note, however, that with a given population size $N$ this approach provides a
statistically consistent framework for inference [4].) Divergence estimators
based on branch lengths [24, 5] may be affected by this; note that this also
depends on whether the deviation actually exceeds the margin of numerical
accuracy of the implementation, which may not be the case. Determining
substitution rates by multiplying the mutation rates and fixation/hitting
probabilities derived from the continuous approximations avoids the problem
entirely [32].

Figure 5: Percentage bias of the eigenvalues of the boundary mutation-drift
Moran model vs. the eigenvalues of the general mutation Moran model in
dependence on the haploid population size $N$: Each plot shows a different
eigenvalue - A) $\lambda_{1}$, B) $\lambda_{2}$, C) $\lambda_{3}$, D)
$\lambda_{4}$ \- of the transition matrix with three different overall scaled
mutation rates: $\theta=c(0.1,0.05,0.01)$, whereby the largest overall scaled
mutation rate corresponds to the slope showing the highest percentage bias,
etc. In all plots, the reference line at 100% corresponds to the value of the
eigenvalue in the limit $\theta\to 0$.

An alternative approach to representing the transition rates of the boundary
mutation-drift Moran model has been proposed by Vogl and Bergman [29]: In the
limit of the small scaled mutation rate $\theta\to 0$, the interior of the
forward (Fokker-Planck) diffusion equation in Eq. (7) becomes a pure drift
diffusion model that can be decomposed into eigenvectors given by (modified)
Gegenbauer polynomials and corresponding eigenvalues $\lambda_{n}=n(n-1)$ for
$n\geq 2$, where $0\leq n\leq\infty$ is the order of the polynomials. The
(modified) Gegenbauer polynomials can be obtained from the (modified) Jacobi
polynomials by a Taylor series expansion up to zeroth order [29]. To
incorporate boundary dynamics similar to those in the boundary mutation-drift
Moran model, boundary terms are introduced to derive a system of inhomogeneous
differential equations yielding the first two eigenvalues $\lambda_{0}=0$, and
$\lambda_{1}=\theta$. Hence, the first non-zero eigenvalue is identical to
that of the transition density of the general mutation model represented by
orthogonal modified Jacobi polynomials (see Sec. 2.2.1). This representation
of the transition density of the boundary mutation-drift Moran model can be
used for inference with our forward-backward algorithm that combines the
particle model and the orthogonal polynomial approach: For this to work, we
must establish a continuous prior distribution $\boldsymbol{\rho(x)}$ for the
boundary mutation-drift Moran model (Sec. 4.1). We will then review the
forward-backward algorithm using orthogonal polynomials for the boundary
mutation-drift Moran model (Sec. 4.2.1), and introduce the corresponding
particle model (Sec. 4.2.2).

###  4.1 A functional as an improper prior

The beta equilibrium distribution Eq. (23) [36] is generally assumed as the
prior distribution $\phi(x{\,|\,}t=S)=\boldsymbol{\rho}(x)$ for the bi-allelic
Moran model in the continuous forward algorithms (Sec. 2.3.2). For the
boundary mutation-drift Moran model, we define the functional

$eq(x{\,|\,}\alpha,\theta)=\lim_{N\to\infty}\begin{cases}\beta\big{(}1-\alpha\theta\int_{1/N}^{1-1/N}\frac{1}{x}\,dx\big{)}&\text{for
$x=0$,}\\\ \alpha\beta\theta\,\frac{1}{x(1-x)}&\text{for $1/N\leq x\leq
1-1/N$,}\\\
\alpha\big{(}1-\beta\theta\int_{1/N}^{1-1/N}\frac{1}{x}\,dx\big{)}&\text{for
$x=1$,}\end{cases}$ (55)

to replace this beta distribution the prior. For the functional to be a valid
prior, taking a sample of size $K$ with $y$ focal alleles with replacement
from it should result in a marginal distribution for $y$ that is equal to the
stationary distribution in Eq. (54). We will prove this here:

  * 1.

(i) For polymorphic samples $1\leq y\leq(K-1)$, the function

$\begin{cases}\alpha\beta\theta\,x^{y-1}(1-x)^{K-y-1}&\text{within
$[1/N,1-1/N]$}\\\ 0&\text{otherwise}\end{cases}$ (56)

is bounded from above by $\alpha\beta\theta\,x^{y-1}(1-x)^{K-y-1}$ within the
interval $[0,1]$ (but converges monotonically towards it with increasing $N$).
By the monotone convergence theorem, the order of taking the limit and
integration can be reversed and we obtain:

$\begin{split}\Pr(1\leq y\leq
K-1{\,|\,}K,\alpha,\theta)&=\alpha\beta\theta\binom{K}{y}\lim_{N\to\infty}\int_{1/N}^{1-1/N}x^{y-1}(1-x)^{K-y-1}\,dx\\\
&=\alpha\beta\theta\frac{K}{y(K-y)}\,.\end{split}$ (57)

  * 2.

(ii) For monomorphic samples, analogous arguments apply. We show the case
$y=K$ (equivalent calculations hold for $y=0$):

$\begin{split}\Pr(y=K{\,|\,}K,\alpha,\theta)&=\lim_{N\to\infty}\int_{1/N}^{1-1/N}x^{K}\alpha\beta\theta\frac{1}{x(1-x)}\\\
&\qquad+\alpha\bigg{(}1-\beta\theta\int_{1/N}^{1-1/N}\frac{1}{x}\,dx\bigg{)}\,dx\\\
&=\lim_{N\to\infty}\int_{1/N}^{1-1/N}\alpha\beta\theta\frac{x^{K-1}}{1-x}+\alpha\bigg{(}1-\beta\theta\int_{1/N}^{1-1/N}\frac{1}{1-x}\,dx\bigg{)}\,dx\\\
&=\alpha+\lim_{N\to\infty}\alpha\beta\theta\int_{1/N}^{1-1/N}(x^{K-1}-1)(1+x+x^{2}+\dots)\,dx\\\
&=\alpha+\lim_{N\to\infty}\alpha\beta\theta\int_{1/N}^{1-1/N}1+x+x^{2}+\dots+x^{K-2}\,dx\\\
&=\alpha-\alpha\beta\theta H_{K-1}\,.\end{split}$ (58)

The functional Eq. (55) can thus be used as an improper prior distribution for
the boundary mutation-drift Moran model. As long as
$K<e^{\frac{1}{\max(\alpha,\beta)\theta}}$ approximately holds, a proper
marginal posterior distribution for the proportion of focal alleles in the
sample will result.

###  4.2 Transition probabilities with small scaled mutation rates

####  4.2.1 Forward-backward with augmented Gegenbauer polynomials

Recall that the forward and backward operators of the Kolmogorov diffusion
equation for the general mutation bi-allelic Moran model can be decomposed
into eigenvectors represented as (modified) Jacobi polynomials and
corresponding eigenvalues (Sec. 2.2.1, Sec. 2.3.2). A similar decomposition
can be achieved for the boundary mutation-drift Moran model; we will reproduce
an outline of the derivation here, details can be found in Vogl and Bergman
[29, Appendix]:

##### Pure drift model

In a first step, the eigenvectors of a pure drift model are obtained by a
zeroth order Taylor series expansion around $\theta$ in Sec. (2.2.1), Sec.
(2.3.2).

The backward eigenvectors of the pure drift model are [3]:

$\begin{cases}B_{0}^{(\alpha,\theta)}(x)&=B_{0}^{(\alpha,0)}(x)=1\\\
B_{1}^{(\alpha,\theta)}(x)&=B_{1}^{(\alpha,0)}(x)=x-\alpha\\\ B_{n\geq
2}^{(\alpha,0)}(x)&=w(x)U_{n}(x)\,,\end{cases}$ (59)

where for $n\geq 2$:

$\begin{split}U_{n}(x)&=\sum_{l=0}^{n-2}(-1)^{l}\binom{n+l}{l}\binom{n-1}{l+1}x^{l}\,.\end{split}$
(60)

are the modified Gegenbauer polynomials (compare [25] Eqs. (12),(13)) with
weight function $w(x)=x(1-x)$. For any $m,n$, the modified Gegenbauer
polynomials fulfill the following orthogonality relationship with respect to
the weight function:

$\int_{0}^{1}U_{n}(x)U_{n}(x)w(x)\,dx=\delta_{n,m}\Delta_{n}\,,$ (61)

where $\delta_{n,m}$ is Kronecker’s delta and
$\Delta_{n}=\frac{(n-1)}{(2n-1)n}$ for $n\geq 2$ are proportionality
constants. The forward eigenvectors of the pure drift model are:

$\begin{cases}F_{0}^{(\alpha,0)}(x)&=\beta\delta(x)+\alpha\delta(x-1)\\\
F_{1}^{(\alpha,0)}(x)&=-\delta(x)+\delta(x-1)\\\ F_{n\geq
2}^{(\alpha,0)}(x)&=-\frac{(-1)^{n}}{n}\delta(x)+U_{n}(x)-\frac{1}{n}\delta(x-1)\,.\end{cases}$
(62)

The corresponding eigenvalues of the system are $\lambda_{0}=0$,
$\lambda_{1}=\theta=0$, and $\lambda_{n}=n(n-1)$ for $n\geq 2$.

##### Boundary mutation-drift model

To introduce mutations arising exclusively at the monomorphic boundaries,
$\lambda_{1}$ is set to $\theta$ with $0<\theta<<1$. The discrete boundary
mutation-drift model is scaled so that mutations enter the polymorphic region
from the respective boundaries at average rates of $\frac{\alpha\theta}{N}$
and $\frac{\beta\theta}{N}$ per Moran drift event in equilibrium; this is
required to maintain a constant equilibrium mutation rate across generations.
In an orthogonal polynomial representation of the transition density of the
continuous boundary mutation-drift model in which the (effective) population
size and thus the generation lengths are not necessarily constant, mutations
enter the polymorphic region from the respective boundaries at rates
$\alpha\theta b_{0}(x,t)$ and $\beta\theta b_{1}(t)$: Here, $b_{0}(t)$ and
$b_{1}(t)$ represent the probability masses both already fixed at the boundary
at time $t$ or expected to fix there imminently by drift. In the polymorphic
region, drift operates at an exponential rate that depends on the sample size
$\kappa$, which is the order of the polynomial expansion and the sample size
$K$. Overall, the forward and backward eigenvectors of the pure drift model of
order $K$ are multiplied with the solution $\tau_{n}(t)$, with $0\leq n\leq
K$, of the following system of first order linear equations describing the
temporal dynamics induced by the boundary mutation ([29], but note the
different weighting there):

$\begin{split}&\frac{d}{dt}\tau_{0}(t)=0\\\
&\frac{d}{dt}\tau_{1}(t)=-\theta\tau_{1}(t)\\\
&\frac{d}{dt}\tau_{n}(t)=-\lambda_{n}\tau_{n}(t)+A_{n}+B_{n}e^{\theta t}\text{
for $n\geq 2$}\\\ \end{split}$ (63)

where

$A_{n}=-\alpha\beta\theta(2n-1)n((-1)^{n}+1)$

and

$B_{n}=-\theta(2n-1)n(b_{0}(t)-\beta)((-1)^{n}\alpha-\beta)\,.$

##### Diagonalized boundary mutation-drift model

The system (Eq. 63) corresponds to a triangular matrix with eigenvalues
$\lambda_{0}=0$, $\lambda_{1}=\theta$, $\lambda_{n}=n(n-1)$ on the main
diagonal. Hence the eigensystem can be easily diagonalized [3]:

The backward eigenfunctions of the boundary mutation-drift model become:

$\begin{cases}B_{0}^{(\alpha,\theta)}(x)&=B_{0}^{(\alpha,0)}(x)=1\\\
B_{1}^{(\alpha,\theta)}(x)&=B_{1}^{(\alpha,0)}(x)=x-\alpha\\\ B_{n\geq
2}^{(\alpha,\theta)}(x)&=B_{n}^{(\alpha,0)}(x)-\vartheta\frac{E_{n}\Delta_{n}}{\lambda_{n}}B_{0}^{(\alpha,0)}(x)-\theta\frac{B_{n}\Delta_{n}}{\lambda_{n}}B_{1}^{(\alpha,0)}(x)\,.\end{cases}$
(64)

with

$\begin{split}\vartheta&=\alpha\beta\theta\,,\\\
E_{n}&=-(n-1)\frac{((-1)^{n}+1)}{\Delta_{n}}\,,\\\
O_{n}&=-(n-1)\frac{(-1)^{n}\alpha-\beta}{\Delta_{n}}\,,\end{split}$ (65)

Equivalently, the forward eigensystem becomes [3]:

$\begin{cases}F_{0}^{(\alpha,\theta)}(x)&=eq(x{\,|\,}\alpha,\theta)\\\
F_{1}^{(\alpha,\theta)}(x)&=neq(x{\,|\,}\alpha,\theta)\\\ F_{n\geq
2}^{(\alpha,\theta)}(x)&=F_{n}^{(\alpha,0)}(x)\,,\end{cases}$ (66)

with the functional $eq(x{\,|\,}\alpha,\theta)$ from Eq. (55), and the similar
functional:

$neq(x{\,|\,}\alpha,\theta)=\lim_{N\to\infty}\begin{cases}-1+\alpha\theta\int_{1/N}^{1-1/N}\frac{1}{x}\,dx&\text{for
$x=0$,}\\\
\theta\,\bigg{(}-\frac{\alpha}{x}+\frac{\beta}{1-x}\bigg{)}&\text{for $1/N\leq
x\leq 1-1/N$,}\\\ 1-\beta\theta\int_{1/N}^{1-1/N}\frac{1}{x}\,dx&\text{for
$x=1$.}\end{cases}$ (67)

The weights are augmented by $\Delta_{0}=\Delta_{1}=1$.

Note that the first diagonalized eigenfunction $F_{0}^{(\alpha,\theta)}(x)$
corresponds to the equilibrium distribution. With this and the second
eigenfunction $F_{1}^{(\alpha,\theta)}(x)$ a quasi-equilibrium state of the
allele proportions in the polymorphic interior can be modelled:

$F_{0}^{(\alpha,\theta)}(x)+\alpha e^{-\theta t}F_{1}^{(\alpha,\theta)}(x)\,.$
(68)

Observe that the rate of approach to equilibrium depends on the mutation rate
$\theta$, which is considerably slower than the rate of drift.

##### Application of the forward-backward algorithm

If the ancestral allele configuration $\boldsymbol{\rho}$(x) is taken to be
distributed according to the functional in Sec. (4.1) with

$\boldsymbol{\rho}(x)=\sum_{n=0}^{\infty}\rho_{n}^{(\alpha,\theta)}F_{n}^{(\alpha,\theta)}(x),$
(69)

and $\rho_{n}^{(\alpha,\theta)}$ obtained according to:

$\rho_{n}^{(\alpha,\theta)}=\frac{1}{\Delta_{n}^{(\alpha,\theta)}}\int_{x=0}^{1}w(x)F_{n}^{(\alpha,\theta)}(x)\boldsymbol{\rho}(x)dx\,.$
(70)

The forward transition density of the population can be expanded to either:

$\phi(x{\,|\,}t,\boldsymbol{\rho})=w(x)\sum_{n=0}^{\infty}\rho_{n}^{(\alpha,\theta)}F_{n}^{(\alpha,0)}(x)\tau_{n}(t)\,,$
(71)

or

$\phi(x{\,|\,}t,\boldsymbol{\rho})=w(x)\sum_{n=0}^{\infty}\rho_{n}^{(\alpha,\theta)}F_{n}^{(\alpha,\theta)}(x)(\tau_{n}(t))^{-}\,,$
(72)

where $(\tau_{n}(t))^{-}$ is the solution to

$\begin{split}&\frac{d}{dt}\tau_{n}(t)=-\lambda_{n}\tau_{n}(t)\end{split}$
(73)

The marginal likelihood of the number of focal alleles $y$ in an extant sample
of size $K$ can similarly be determined via:

$\begin{split}\Pr(y{\,|\,}K,\alpha,\theta,t=0,\boldsymbol{\rho})=\sum_{n=0}^{K}\rho_{n}^{(\alpha,\theta)}d_{n}(K,y)\Delta_{n}\tau_{n}(t=0)\,.\end{split}$
(74)

or

$\begin{split}\Pr(y{\,|\,}K,\alpha,\theta,t=0,\boldsymbol{\rho})&=\sum_{n=0}^{K}\rho_{n}^{(\alpha,\theta)}d_{n}(K,y)\Delta_{n}(\tau_{n}(t=0))^{-}\,,\end{split}$
(75)

with

$d_{n}(K,y)=\frac{1}{\Delta_{n}}\int_{x=0}^{1}\binom{K}{y}x^{y}(1-x)^{K-y}\,F_{n}^{(\alpha,0)}(x)\,dx\,.$
(76)

##### Summary

A forward-backward algorithm approach using orthogonal polynomials can easily
be applied to the boundary mutation-drift Moran model. Note that using the un-
diagonalized version of the eigensystem is particularly efficient when
mutation parameters change over time, both with respect to time and numerical
accuracy, since this change does not affect the eigenvectors but only the
temporal system. Although the boundary mutation-drift Moran model is an
approximation to the general mutation bi-allelic Moran model, it may provide
more accurate results in an inference framework than direct use of the general
mutation model through this numerical advantage. The increased numerical
accuracy may allow the use of higher sample sizes, which is necessary for many
data analyses.

####  4.2.2 Particle transition probabilities

As with the general mutation model, Eq. (74) and Eq. (75) can be determined
with an appropriate particle model. This can again be constructed using the
ratios of the probabilities of the appropriate configurations of ancestral
focal alleles $\upsilon$ in a sample size of $\kappa$ (see Eq. 34).

Recall that the indicator variable $z_{\kappa+1}$ is one if the $\kappa+1$st
allele is of the focal type and zero otherwise. Note that the transition
probabilities from $z_{\kappa}$ to $z_{\kappa+1}$ must hold in equilibrium.
Recall that for the general mutation model, we have:

$\begin{split}{\Pr}(z_{\kappa+1}{\,|\,}\upsilon,\kappa,\alpha,\theta)&=\frac{{\Pr}(z_{\kappa+1}+\upsilon{\,|\,}\kappa+1,\alpha,\theta)}{{\Pr}(\upsilon{\,|\,}\kappa,\alpha,\theta)}\\\
&=\frac{(\upsilon+\alpha\theta)^{1-z_{\kappa}}(\kappa-\upsilon+\beta\theta)^{z_{\kappa}}}{\kappa+\theta}\,.\end{split}$
(77)

Note that the terms of
${\Pr}(z_{\kappa+1}+\upsilon{\,|\,}\kappa+1,\alpha,\theta)$ and
${\Pr}(\upsilon{\,|\,}\kappa,\alpha,\theta)$ that only hold in equilibrium
cancel out. Therefore, the transition probabilities for the general mutation
model as given in Eq. (34) hold in non-equilibrium situations.

For the boundary mutation-drift model, transition probabilities need to be
consistent with the equilibrium probability Eq. (55). For
$1\leq\upsilon\leq(\kappa-1)$, no mutations are possible and the transition
probabilities are independent of the mutation parameters $\alpha$ and
$\theta$:

$\begin{split}{\Pr}(z_{\kappa+1}{\,|\,}\upsilon,\alpha,\theta)&=\frac{{\Pr}(z_{\kappa+1}+\upsilon{\,|\,}\kappa+1,\alpha,\theta)}{{\Pr}(\upsilon{\,|\,}\kappa,\alpha,\theta)}\\\
&=\frac{(\kappa-\upsilon)^{z_{\kappa}}\upsilon^{1-z_{\kappa}}}{\kappa}\,.\end{split}$
(78)

For $\upsilon=0$, the transition probability consistent with that in Eq. (78)
is:

$\begin{split}{\Pr}(z_{\kappa+1}{\,|\,}\upsilon=0,\alpha,\theta)&=\frac{{\Pr}(z_{\kappa+1}{\,|\,}\kappa+1,\alpha,\theta)}{{\Pr}(\upsilon=0{\,|\,}\kappa,\alpha,\theta)}\\\
&=\frac{(\frac{\alpha\theta}{\kappa})^{z_{\kappa}}(1-\alpha\theta
H_{\kappa})^{1-z_{\kappa}}}{1-\alpha\theta H_{\kappa-1}}\,,\end{split}$ (79)

and the analogous transition rate holds for $\upsilon=\kappa$.

For the transition from $\kappa=1$ to $\kappa=0$ (or the reverse forward in
time), the transition probabilities of the boundary mutation-drift models are
identical to those of the general mutation model: Specifically, the
probabilities are $\alpha$ and $\beta$ to zeroth order in $\theta$.

##### Temporal dynamics

The backward temporal dynamics of the particle model are clearly equivalent to
those in Sec. (3.5) however with eigenvalues $\lambda_{0}=0$,
$\lambda_{1}=\theta$, and $\lambda_{n}=n(n-1)$ for $n\geq 2$.

##### Summary

A particle model can be used to trace sample allele configurations in the
boundary mutation-drift Moran model forward and backward in time. The
transition probabilities defined through it are sample consistent, and the
first non-zero eigenvalue (defining the speed of approach to equilibrium) is
identical to that of the general mutation Moran model and independent of the
sample size. In phylogenetic settings, where split times (in generations)
between species or populations are large compared to effective population
sizes, this eigenvalue dominates.

##  5 Conclusions

In population genetics, small to moderately sized extant samples are generally
used to infer parameters of past population processes. The evolutionary
trajectory of the population can, under certain assumptions, be described by
the bi-allelic decoupled Moran model in the diffusion limit.

To determine the marginal likelihood of sample allele proportions in the
extant sample, the backwards N-particle diffusion of the extant sampling
distribution can be decomposed into eigenvectors of orthogonal polynomials and
corresponding eigenvectors. The backward orthogonal polynomial approach is
computationally efficient and accommodates temporally changing effective
population sizes (and therefore mutation rates) or mutation biases. This is
even more true when scaled overall mutation rates are small and a boundary
mutation-drift Moran model is assumed to describe the evolutionary dynamics:
if the system is not diagonalized, time inhomogeneous dynamics allow for
effient numerics.

The probabilities of ancestral sample allele configurations can be inferred
with a so-called particle model augmented with backwards temporal dynamics: It
directly traces the size and allele proportions of a sample both forward (urn
model) and backward (coalescent) in time. Empirically, this approach may be
useful for inferring the timing of mutation-drift events at a specific locus.
However, its main value lies in its equivalence with the backwards orthogonal
polynomial method: The backward orthogonal polynomial expansion of an extant
sample to $n$th order has the genealogy of the sample embedded within it. It
therefore contains all the information traditionally modelled by the
coalescent, but is more tractable and numerically superior.

In order to infer the probability of ancestral population allele proportions
at any time in the past, the particle model needs to be complemented with a
forward in time orthogonal polynomial approach. The result is a full forward-
backward algorithm equivalent to the forward-backward algorithm using only
orthogonal polynomials.

## Acknowledgments

We thank the present and past members of the doctorate college population
genetics, especially Juraj Bergman, for stimulating discussions. CV and SP’s
research was supported by the Austrian Science Fund (FWF): DK W1225-B20; LCM’s
by the School of Biology at the University of St.Andrews and also partially
funded through Vienna Science and Technology Fund (WWTF) [MA016-061].

## References

  * Abramowitz and Stegun [1970] M Abramowitz and IA Stegun, editors.  _Handbook of Mathematical Functions_.  Dover, 9th ed. edition, 1970.
  * Beaumont et al. [2002] M Beaumont, W Zhang, and DJ Balding.  Approximate Bayesian Computation in population genetics.  _Genetics_ , 162:2025–2035, 2002.
  * Bergman et al. [2018] J Bergman, D Schrempf, C Kosiol, and C Vogl.  Inference in population genetics using forward and backward, discrete and continuous time processes.  _Journal of Theoretical Biology_ , 439:166–180, 2018.
  * Borges and Kosiol [2020] R Borges and C Kosiol.  Consistency and identifiability of the polymorphism-aware phylogenetic models.  _Journal of Theoretical Biology_ , 486:Letter to the Editor, 2020.
  * Borges et al. [2019] R Borges, G Szöllősi, and C Kosiol.  Quantifying GC-biased gene conversion in great ape genomes using polymorphism-aware models.  _Genetics_ , 212:1321–1336, 2019.
  * Burden and Griffiths [2019] C Burden and RC Griffiths.  The transition distribution of a sample from a Wright–Fisher diffusion with general small mutation rates.  _Journal of Mathematical Biology_ , 79:2315–2342, 2019\.
  * Clemente and Vogl [2012] F Clemente and C Vogl.  Unconstrained evolution in short introns?—an analysis of genome-wide polymorphism and divergence data from _Drosophila_.  _Journal of Evolutionary Biology_ , 25(10):1975–90, 2012.
  * Estoup and Guillemaud [2010] A Estoup and T Guillemaud.  Reconstructing routes of invasion using genetic data: why, how and so what?  _Molecular Ecology_ , 19:4113–4130, 2010.
  * Etheridge [2011] A Etheridge.  _Some Mathematical Models for Population Genetics_.  Springer, 1st edition, 2011.
  * Ewens [2004] WJ Ewens.  _Mathematical Population Genetics_.  Springer, N.Y., 2nd edition, 2004.
  * Faisal et al. [2015] M Faisal, A Futschik, and C Vogl.  Exact likelihood calculation under the infinite sites model.  _Computation_ , 3:701–713, 2015.
  * Fraimout et al. [2017] A Fraimout, V Debat, S Fellous, RA Hufbauer, J Foucaud, P Pudlo, JM5 Marin, DK Price, J Cattel, X Chen, M Deprá, P François Duyck, C Guedot, M Kenis, MT Kimura, G Loeb, A Loiseau, I Martinez-Sañudo, M Pascual, M Polihronakis Richmond, P Shearer, N Singh, K Tamura, A Xuéreb, J Zhang, and A Estoup.  Deciphering the routes of invasion of _Drosophila suzukii_ by means of ABC random forest.  _Molecular Biology and Evolution_ , 34:980–996, 2017.
  * Hein et al. [2005] J Hein, MH Schierup, and C Wiuf.  _Gene Genealogies, Variation, and Evolution: a Primer in Coalescent Theory._ Oxford University Press, 2005.
  * Hudson [1985] RR Hudson.  Gene genealogies and the coalescent process.  _Oxford Surveys in Evolutionary Biology_ , 7:1–44, 1985\.
  * Kimura [1971] M Kimura.  Theoretical foundation of population genetics at the molecular level.  _Theoretical Population Biology_ , 2:893–903, 1971.
  * Kingman [1982] JFC Kingman.  On the genealogy of large populations.  _Journal of Applied Probability_ , 19A:27–43, 1982.
  * Lynch et al. [2016] M Lynch, M Ackerman, J Gout, H Long, W Sung, W Thomas, and P Foster.  Genetic drift, selection and the evolution of the mutation rate.  _Nature_ , 17:704–714, 2016.
  * Moran [1958a] PAP Moran.  Random processes in genetics.  _Proceedings of the Cambridge Philosophical Society_ , 54:60–71, 1958a.
  * Moran [1958b] PAP Moran.  The effect of selection in a haploid genetic population.  _Proceedings of the Cambridge Philosophical Society_ , 54:463–467, 1958b.
  * Moran [1962] PAP Moran.  _Statistical processes of evolutionary theory._ Clarendon Press, Oxford, 1962.
  * Myers et al. [2005] S Myers, L Bottolo, C. Freeman, G McVean, and P Donnelly.  A fine-scale map of recombination rates and hotspots across the human genome.  _Science_ , 310:310–324, 2005.
  * Parsch et al. [2010] J Parsch, S Novozhilov, SS Saminadin-Peter, KM Wong, and P Andolfatto.  On the utility of short intron sequences as a reference for the detection of positive and negative selection in _Drosophila._ _Molecular Biology and Evolution_ , 27:1226–1234, 2010\.
  * Rabiner and Juang [1986] LR Rabiner and BH Juang.  An introduction to hidden Markov models.  _EEE ASSP Magazine_ , 3:4–15, 1986.
  * Schrempf et al. [2016] D Schrempf, BQ Minh, N De Maio, A von Haeseler, and C Kosiol.  Reversible polymorphism-aware phylogenetic models and their application to tree inference.  _Journal of Theoretical Biology_ , 407:362–370, 2016.
  * Song and Steinrücken [2012] Y Song and M Steinrücken.  A simple method for finding explicit analytic transition densities of diffusion processes with general diploid selection.  _Genetics_ , 190:1117–1129, 2012.
  * Stephens and Donnelly [2000] M Stephens and P Donnelly.  Inference in molecular population genetics.  _Journal of the Royal Statistical Society, Series B_ , 62:605–655, 2000.
  * Taylor and Karlin [1998] HM Taylor and S Karlin.  _An Introduction to Stochastic Modeling._ Academic Press, San Diego, 1998.
  * Vogl [2014] C Vogl.  Estimating the scaled mutation rate and mutation bias with site frequency data.  _Theoretical Population Biology_ , 98:19—27, 2014.
  * Vogl and Bergman [2016] C Vogl and J Bergman.  Computation of the likelihood of joint site frequency spectra using orthogonal polynomials.  _Computation_ , 4:6, 2016.
  * Vogl and Clemente [2012] C Vogl and F Clemente.  The allele-frequency spectrum in a decoupled moran model with mutation, drift, and directional selection, assuming small mutation rates.  _Theoretical Population Genetics_ , 81:197–209, 2012.
  * Vogl and Futschik [2010] C Vogl and A Futschik.  Hidden markov models in biology.  In O. Carugo and F. Eisenhaber, editors, _Biological Data Mining._ , Methods in Molecular Biology. Humana Press, 2010.
  * Vogl and Mikula [2021] C Vogl and LC Mikula.  A nearly-neutral biallelic Moran model with biased mutation and linear and quadratic selection.  _Theoretical Population Biology_ , 139:1–17, 2021.
  * Vogl et al. [2020] C Vogl, LC Mikula, and CJ Burden.  Maximum likelihood estimators for scaled mutation rates in an equilibrium mutation-drift model.  _Theoretical Population Biology_ , 134:106–118, 2020.
  * Wakeley [2009] J Wakeley.  _Coalescent Theory, an Introduction_.  Roberts and Co., 2009.
  * Watterson [1975] GA Watterson.  On the number of segregating sites in genetical models without recombination.  _Theoretical Population Biology_ , 7:256–276, 1975.
  * Wright [1931] S Wright.  Evolution in Mendelian populations.  _Genetics_ , 16:97–159, 1931.
  * Wu [2010] Y Wu.  Exact computation of coalescent likelihood for panmictic and subdivided populations under the infinite sites model.  _IEEE/ACM Transactions on Computational Biology and Bioinformatics_ , 7:611–618, 2010.
  * Zhao et al. [2013] L Zhao, M Lascoux, ADJ Overall, and D Waxman.  The characteristic trajectory of a fixing allele: a consequence of fictitious selection that arises from conditioning.  _Genetics_ , 195:993––1006, 2013.