Preliminaries

Let $$ \mathcal{T}$$ be a binary rooted time tree with N taxa, and data D associated with the tips, such as a multiple sequence alignment with L sites, morphological data, or geographic locations. The posterior density of a phylogenetic model is described by

Let t_i be the height (time) of node i. Each node i in $$ \mathcal{T}$$, except for the root, is associated with a parental branch length τ_i (the height difference between i and its parent) and a parental branch substitution rate $$ r_i=r\left({\mathcal{R}}_i\right)$$. In an uncorrelated relaxed clock model, each of the 2N − 2 elements in $$ {\overrightarrow{\mathcal{R}}}^{}$$ are independently distributed under the prior $$ p\left({\overrightarrow{\mathcal{R}}}^{}\right|\sigma )$$.

The posterior distribution is sampled by the Metropolis-Hastings-Green MCMC algorithm [7, 8, 34], under which the probability of accepting proposed state x′ from state x is equal to:

Branch rate parameterisations

In Bayesian inference, the way parameters are represented in the model can affect the mixing ability of the model and the meaning of the model itself [36]. Three methods for parameterising substitution rates are described below. Each parameterisation is associated with i) an abstraction of the branch rate vector $$ {\overrightarrow{\mathcal{R}}}^{}$$, ii) some function for transforming this parameter into unabstracted branch rates $$ r\left({\overrightarrow{\mathcal{R}}}^{}\right)$$, and iii) a prior density function of the abstraction $$ p\left({\overrightarrow{\mathcal{R}}}^{}\right|\sigma )$$. The three methods are summarised in Fig 1.

Fig 1

Branch rate parameterisations.

Top left: the prior density of a branch rate r under a Log-normal(−0.5σ², σ) distribution (with its mean fixed at 1). The function for transforming $$ \mathcal{R}$$ into branch rates $$ r\left(\mathcal{R}\right)$$ is depicted for real (top right), cat (bottom left), and quant (bottom right). For visualisation purposes, there are only 10 bins/pieces displayed, however in practice we use 2N − 2 bins for cat and 100 pieces for quant. The first and final quant pieces are equal to the underlying function (solid lines) however the pieces in between use linear approximations of this function (dashed lines).

1. Real rates

The natural (and unabstracted) parameterisation of a substitution rate is a real number $$ {\mathcal{R}}_i\in \mathbb{R},{\mathcal{R}}_i>0$$ which is equal to the rate itself. Under the real parameterisation:

Under a log-normal clock prior $$ p\left({\overrightarrow{\mathcal{R}}}^{}\right|\sigma )$$, rates are distributed with a mean of 1:

where μ = −0.5σ² is set such that the expected value is 1. In this article we only consider log-normal clock priors, however the methods discussed are general.

Zhang and Drummond 2020 introduced a series of tree operators which propose node heights and branch rates, such that the resulting genetic distances (r_i × τ_i) remain constant [33]. These operators account for correlations between branch rates and branch times. By keeping the genetic distance of each branch constant, the likelihood is unaltered by the proposal.

2. Categories

The category parameterisation cat is an abstraction of the real parameterisation. Each of the 2N − 2 branches are assigned an integer from 0 to n − 1:

These integers correspond to n rate categories (Fig 1). Let f(x|σ) be the probability density function (PDF) and let $$ F\left(x\right|\sigma )=\underset{0}{\overset{x}{\int }}f\left(t\right|\sigma )dt$$ be the cumulative distribution function (CDF) of the prior distribution used by the underlying real clock model (a log-normal distribution for the purposes of this article). In the cat parameterisation, f(x|σ) is discretised into n bins and each element within $$ {\overrightarrow{\mathcal{R}}}^{}$$ points to one such bin. The rate of each bin is equal to its median value:

where F⁻¹ is the inverse cumulative distribution function (i-CDF). The domain of $$ {\overrightarrow{\mathcal{R}}}^{}$$ is uniformly distributed under the prior:

The key advantage of the cat parameterisation is the removal of a term from the posterior density (Eq 1), or more accurately the replacement of a non-trivial $$ p\left({\overrightarrow{\mathcal{R}}}^{}\right|\sigma )$$ term with that of a uniform prior. This may facilitate efficient exploration of the posterior distribution by MCMC.

This parameterisation has been widely used in BEAST and BEAST 2 analyses [3]. However, the recently developed constant distance operators—which are incompatible with the cat parameterisation—can yield an increase in mixing rate under real by up to an order of magnitude over that of cat, depending on the dataset [33].

3. Quantiles

Finally, rates can be parameterised as real numbers describing the rate’s quantile with respect to some underlying clock model distribution [18]. Under the quant parameterisation, each of the 2N − 2 elements in $$ {\overrightarrow{\mathcal{R}}}^{}$$ are uniformly distributed.

Transforming these quantiles into rates invokes the i-CDF of the underlying real clock model distribution. Evaluation of the log-normal i-CDF is has high computational costs and therefore an approximation is used instead.

where $$ {\widehat{F}}^{-1}$$ is a linear piecewise approximation with 100 pieces. While this approach has clear similarities with cat, the domain of rates here is continuous instead of discrete. In this project we extended the family of constant distance operators [33] so that they are compatible with quant. Further details on the quant piecewise approximation and constant distance operators can be found in S1 Appendix.

Clock model operators

The weight of an operator determines the probability of the operator being selected. Weights are typically fixed throughout MCMC. In BEAST 2, operators can have their own tunable parameter s, which determines the step size of the operator. This term is tuned over the course of MCMC by achieving a target acceptance rate, typically 0.234 [10, 37, 38]. We define clock model operators as those which generate proposals for either $$ {\overrightarrow{\mathcal{R}}}^{}$$ or σ. Pre-existing BEAST 2 clock model operators are summarised in Table 1, and further operators are introduced throughout this paper.

Table 1

Summary of pre-existing BEAST 2 operators.

Operator	Description	Parameters	Parameterisations
RandomWalk	Moves a single element by a tunable amount.	$$	cat, real, quant
Scale	Applies RandomWalk on the log-transformation (suitable for parameters with positive domains).	$$	cat, real, quant
Interval	Applies RandomWalk on the logit-transformation (suitable for parameters with upper and lower limits).	$$	quant
Swap	Swaps two random elements within the vector [3].	$$	cat, real, quant
Uniform	Resamples one element in the vector from a uniform distribution.	$$	cat, quant
ConstantDistance	Adjusts an internal node height and recalculates all incident branch rates such that the genetic distances remain constant [33].	$$	real, quant
SimpleDistance	Applies ConstantDistance to the root node [33].	$$	real, quant
SmallPulley	Proposes new branch rates incident to the root such that their combined genetic distance is constant [33].	$$	real, quant
CisScale	Applies Scale to σ. Then recomputes all rates such that their quantiles are constant (for real [33]) or recomputes all quantiles such that their rates are constant (quant).	$$	real, quant

Summary of pre-existing BEAST 2 operators, which apply to either branch rates $$ {\overrightarrow{\mathcal{R}}}^{}$$ or the clock standard deviation σ, and the substitution rate parameterisation they apply to. ConstantDistance and SimpleDistance also adjust node heights t.

The family of constant distance operators (ConstantDistance, SimpleDistance, and SmallPulley [33]) are best suited for larger datasets (or datasets with strong signal) where the likelihood distribution is peaked (Fig 1). While simple one dimensional operators such as RandomWalk or Scale must take small steps in order to stay “on the ridge” of the likelihood function, the constant distance operators “wander along the ridge” by ensuring that genetic distances are constant after the proposal.

Scale and CisScale both operate on the clock model standard deviation σ however they behave differently in the real and quant parameterisations (Fig 2). In real, large proposals of σ → σ′ made by Scale could incur large penalties in the clock model prior density $$ p\left({\overrightarrow{\mathcal{R}}}^{}\right|{\sigma }^{\prime })$$ and thus may be rejected quite often. This led to the development of the fast clock scaler [33] (herein referred to as CisScale). This operator recomputes all branch rates $$ {\overrightarrow{\mathcal{R}}}^{}\to {\overrightarrow{\mathcal{R}}}^{\prime }$$ such that their quantiles under the new clock model prior remain constant $$ p\left({\overrightarrow{\mathcal{R}}}^{}\right|\sigma )=p\left({\overrightarrow{\mathcal{R}}}^{\prime }\right|{\sigma }^{\prime })$$. In contrast, a proposal made by Scale σ → σ′ under the quant parameterisation implicitly alters all branch rates $$ r\left({\overrightarrow{\mathcal{R}}}^{}\right)$$ while leaving the quantiles $$ {\overrightarrow{\mathcal{R}}}^{}$$ themselves constant. Whereas, application of CisScale under quant results in all quantiles being recomputed $$ \overrightarrow{\mathcal{R}}\to {\overrightarrow{\mathcal{R}}}^{\prime }$$ such that their rates are constant, i.e. $$ r\left({\overrightarrow{\mathcal{R}}}^{}\right)=r\left({\overrightarrow{\mathcal{R}}}^{\prime }\right)$$. In summary, Scale and CisScale propose rates/quantiles in the opposite (trans) or same (cis) space that the clock model is parameterised under (Fig 2).

Fig 2

Clock standard deviation scale operators.

The two operators above propose a clock standard deviation σ → σ′. Then, either the new quantiles are such that the rates remain constant (“New quantiles”, above) or the new rates are such that the quantiles remain constant (“New rates”). In the real parameterisation, these two operators are known as Scale and CisScale, respectively. Whereas, in quant, they are known as CisScale and Scale.

Adaptive operator weighting

It is not always clear which operator weighting scheme is best for a given dataset. In this article we introduce AdaptiveOperatorSampler—a meta-operator which learns the weights of other operators during MCMC and then samples these operators according to their learned weights. This meta-operator undergoes three phases. In the first phase (burn-in), AdaptiveOperatorSampler samples from its set of sub-operators uniformly at random. In the second phase (learn-in), the meta-operator starts learning several terms detailed below whilst continuing to sample operators uniformly at random. In its final phase, AdaptiveOperatorSampler samples operators (denoted by ω) using the following distribution:

Under Eqs 11 and 12, operators which effect larger changes on the parameters of interest, in shorter runtime, are sampled with greater probabilities. Division of the squared distance by a parameter’s sample variance $$ {\sigma }_p^2$$ enables comparison between numerical parameters which exist in different spaces.

Datasets which contain very poor signal (or small L) are likely to mix better when more weight is placed on bold operators (Fig 3). We therefore introduce the SampleFromPrior($$ {\overrightarrow{x}}^{}$$) operator. This operator resamples ψ randomly selected elements within vector $$ {\overrightarrow{x}}^{}$$ from their prior distributions, where $$ \psi \sim \text{Binomial}(n=|{\overrightarrow{x}}^{}|,p=\frac{s}{|{\overrightarrow{x}}^{}|})$$ for tunable term s. SampleFromPrior is included among the set of operators under AdaptiveOperatorSampler and serves to make the boldest proposals for datasets with poor signal.

Fig 3

Traversing likelihood space.

The z-axes above are the log-likelihoods of the genetic distance r × τ between two simulated nucleic acid sequences of length L, under the Jukes-Cantor substitution model [40]. Two possible proposals from the current state (white circle) are depicted. These proposals are generated by the RandomWalk (RW) and ConstantDistance (CD) operators. In the low signal dataset (L = 0.1kb), both operators can traverse the likelihood space effectively. However, the exact same proposal by RandomWalk incurs a much larger likelihood penalty in the L = 0.5kb dataset by “falling off the ridge”, in contrast to ConstantDistance which “walks along the ridge”. This discrepancy is even stronger for larger datasets and thus necessitates the use of operators such as ConstantDistance which account for correlations between branch lengths and rates.

In this article we apply three instances of the AdaptiveOperatorSampler meta-operator to the real, cat, and quant parameterisations. These are summarised in Table 2.

Table 2

Summary of AdaptiveOperatorSampler operators and their parameters of interest (POI).

Meta-operator	POI	Operators
AdaptiveOperatorSampler(σ)	σ	CisScale($$)
		RandomWalk(σ)
		Scale(σ)
		SampleFromPrior(σ)
AdaptiveOperatorSampler($$)	$$	ConstantDistance($$)
		RandomWalk($$)
		Scale($$)
		Interval($$)
		Swap($$)
		SampleFromPrior($$)
AdaptiveOperatorSampler(root)	$$	SimpleDistance($$)
		SmallPulley($$)
AdaptiveOperatorSampler(leaf)	$$	ConstantDistance($$)
		LeafAVMVN($$)
		RandomWalk($$)
		Scale($$)
		Interval($$)
		Swap($$)
		SampleFromPrior($$)
AdaptiveOperatorSampler(internal)	$$	ConstantDistance($$)
		RandomWalk($$)
		Scale($$)
		Interval($$)
		Swap($$)
		SampleFromPrior($$)
AdaptiveOperatorSampler(NER)	$$	NER{}
		$$

Different operators are applicable to different substitution rate parameterisations (Table 1). Nodes are broken down into regions to enable operators to be weighted according to their dimension. AdaptiveOperatorSampler(root) applies the root-targeting constant distance operators only [33] while AdaptiveOperatorSampler($$ {\overrightarrow{\mathcal{R}}}^{}$$) targets all rates and all nodes heights t. Nodes are further broken down into leaf rate $$ {\overrightarrow{\mathcal{R}}}_{\text{leaf}}^{}$$ and internal node rate $$ {\overrightarrow{\mathcal{R}}}_{\text{int}}^{}$$ operators. This facilitates suitable weighting of the LeafAVMVN operator, which is only applicable to leaf nodes (4. Screening operators for acceptance rate using simulated data.). In this setup, the RandomWalk(x), Scale(x), and SampleFromPrior(x) operators apply to the corresponding set of branch rates x, whereas ConstantDistance(x, t) is only applicable to internal nodes which have at least one child of type $$ x\in \{{\overrightarrow{\mathcal{R}}}_{\text{leaf}}^{},{\overrightarrow{\mathcal{R}}}_{\text{int}}^{}\}$$. Finally, the Robinson-Foulds distance between trees before and after every proposal accept is used to train the weights behind AdaptiveOperatorSampler(NER) (see Narrow exchange rate). In the special case of NER proposals, the RF distance is always equal to 1.

Bactrian proposal kernel

The step size of a proposal kernel q(x′|x) should be such that the proposed state x′ is sufficiently far from the current state x to explore vast areas of parameter space, but not so large that the proposal is rejected too often [37]. Operators which attain an acceptance probability of 0.234 are often considered to have arrived at a suitable midpoint between these two extremes [10, 37]. The standard uniform distribution kernel has recently been challenged by the family of Bactrian kernels [31, 32]. The (Gaussian) Bactrian(m) distribution is defined as the sum of two normal distributions:

Fig 4

The Bactrian proposal kernel.

The step size made under a Bactrian proposal kernel is equal to sΣ where Σ is drawn from the above distribution and s is tunable.

In this article we compare the abilities of uniform and Bactrian(0.95) proposal kernels at estimating clock model parameters. The clock model operators which these proposal kernels apply to are described in S1 Appendix.

Narrow exchange rate

The NarrowExchange operator [41], used widely in BEAST [9, 42] and BEAST 2 [10], is similar to nearest neighbour interchange [43], and works as follows (Fig 5):

Step 1. Sample an internal/root node E from tree $$ \mathcal{T}$$, where E has grandchildren.

Step 2. Identify the child of E with the greater height. Denote this child as D and its sibling as C (i.e. t_D > t_C). If D is a leaf node, then reject the proposal.

Step 3. Randomly identify the two children of D as A and B.

Step 4. Relocate the B − D branch onto the C − E branch, so that B and C become siblings and their parent is D. All node heights remain constant.

Fig 5

Depiction of NarrowExchange and NarrowExchangeRate operators.

Proposals are denoted by $$ \mathcal{T}\to {\mathcal{T}}^{\prime }$$. The vertical axes correspond to node heights t. In the bottom figure, branch rates r are indicated by line width and therefore genetic distances are equal to the width of each branch multiplied by its length. In this example, the $$ {\mathcal{D}}_{AE}$$ and $$ {\mathcal{D}}_{CE}$$ constraints are satisfied.

We hypothesised that if NarrowExchange was adapted to the relaxed clock model by ensuring that genetic distances remain constant after the proposal (analogous to constant distance operators [33]), then its ability to traverse the state space may improve.

Here, we present the NarrowExchangeRate (NER) operator. Let r_A, r_B, r_C, and r_D be the substitution rates of nodes A, B, C, and D, respectively. In addition to the modest topological change applied by NarrowExchange, NER also proposes new branch rates $$ {r_A}^{\prime }$$, $$ {r_B}^{\prime }$$, $$ {r_C}^{\prime }$$, and $$ {r_D}^{\prime }$$. While NER does not alter t_D (i.e. $$ {t_D}^{\prime }\leftarrow t_D$$), we also consider NERw—a special case of the NER operator which embarks t_D on a random walk:

Unfortunately, there is no solution to all six $$ {\mathcal{D}}_{ij}$$ constraints unless non-positive rates or illegal trees are permitted. Therefore instead of conserving all six pairwise distances, NER conserves a subset of distances. It is not immediately clear which subset should be conserved.

2. Solving Jacobian determinants

The determinant of the Jacobian matrix J is required for computing the Green ratio of this proposal. J is defined as

Solving the determinant |J| invokes standard analytical differentiation and linear algebra libraries of MATLAB. 6 of the 54 solvable operators were found to have |J| = 0, corresponding to irreversible proposals, and were discarded.

4. Screening operators for acceptance rate using simulated data

Selecting the best NER variant to proceed to benchmarking on empirical data (Results) was determined by performing MCMC on simulated data, measuring the acceptance rates of each of the 96 NER/NERw variants, and comparing them with the null operator NER{} / NarrowExchange. In total, there were 300 simulated datasets each with N = 30 taxa and varying alignment lengths.

These experiments showed that NER variants which satisfied the genetic distances between nodes B and A (i.e. $$ {\mathcal{D}}_{AB}$$) or between B and C (i.e. $$ {\mathcal{D}}_{BC}$$) usually performed worse than the standard NarrowExchange operator (Fig 6). This is an intuitive result. If there is high uncertainty in the positioning of B with respect to A and C, then there is no value in respecting either of these distance constraints, and the proposals made to the rates may often be too extreme or the Green ratio |J| too small for the proposal to be accepted.

Fig 6

Screening of NER and NERw variants by acceptance rate.

Top left: comparison of NER variants with the null operator NER{} = NarrowExchange. Each operator is represented by a single point, uniquely encoded by the point stylings. The number of times each operator is proposed and accepted is compared with that of NER{}, and one-sided z-tests are performed to assess the statistical significance between the two acceptance rates (p = 0.001). This process is repeated across 300 simulated datasets. The axes of each plot are the proportion of these 300 simulations for which there is evidence that the operator is significantly better than NER{} (x-axis) or worse than NER{} (y-axis). Top right: comparison of NER and NERw acceptance rates. Each point is one NER/NERw variant from a single simulation. Bottom: relationship between the acceptance rates α of $$ \text{NER}\{{\mathcal{D}}_{AE},{\mathcal{D}}_{BE},{\mathcal{D}}_{CE}\}$$ and NER{} with the clock model standard deviation σ and the number of sites L. Each point is a single simulation.

Fig 6 also revealed a cluster of NER variants which—under the conditions of the simulation—performed better than the null operator NER{} around 25% of the time and performed worse around 10% of the time. One such operator was $$ \text{NER}\{{\mathcal{D}}_{AE},{\mathcal{D}}_{BE},{\mathcal{D}}_{CE}\}$$ (Algorithm 1). This variant conserves the genetic distance between nodes A, B, C and their grandparent E. This operator performed well when branch rates had a large variance (σ > 0.5), corresponding to non clock-like data. On the other hand, the null operator NER{} performed better on shorter sequences (L < 1kb) with weaker signal. Overall, $$ \text{NER}\{{\mathcal{D}}_{AE},{\mathcal{D}}_{BE},{\mathcal{D}}_{CE}\}$$ outperformed the standard NarrowExchange operator when the data was not clock-like and contained sufficient signal.

Finally, this initial screening showed that applying a (Bactrian) random walk to the node height t_D made the operator worse. This effect was most dominant for the NER variants which satisfied distance constraints (i.e. the operators which are not NER{}).

Although there were several operators which behaved equivalently during this initial screening process, we selected $$ \text{NER}\{{\mathcal{D}}_{AE},{\mathcal{D}}_{BE},{\mathcal{D}}_{CE}\}$$ to proceed to benchmarking (Results). Due to the apparent sensitivity of NER operators to the data, we introduce the adaptive operator AdaptiveOperatorSampler(NER) which allows the operator scheme to fall back on the standard NarrowExchange in the event of $$ \text{NER}\{{\mathcal{D}}_{AE},{\mathcal{D}}_{BE},{\mathcal{D}}_{CE}\}$$ performing poorly (Table 2).

An adaptive leaf rate operator

The adaptable variance multivariate normal (AVMVN) kernel learns correlations between parameters during MCMC [30, 42]. Baele et al. 2017 observed a large increase (≈ 5 − 10×) in sampling efficiency from using the AVMVN kernel substitution model parameters [30]. Here, we consider application of the AVMVN kernel to the branch rates of leaf nodes. This operator, referred to as LeafAVMVN, is not readily applicable to internal node branch rates due to their dependencies on tree topology.

Leaf rate AVMVN kernel

The AVMVN kernel assumes its parameters live in $$ x\in {\mathbb{R}}^N$$ for taxon count N and that these parameters follow a multivariate normal distribution with covariance matrix Σ_N. Hence, the kernel operates on the logarithmic or logistic transformation of the N leaf branch rates, depending on the rate parameterisation:

where r is a real rate and q is a rate quantile. The AVMVN probability density is defined by

where $$ \mathcal{\text{MVN}}$$ is the multivariate normal probability density. β = 0.05 is a constant which determines the fraction of the proposal determined by the identity matrix $$ {\mathbb{I}}_N$$, as opposed to the covariance matrix Σ_D which is trained during MCMC. Our BEAST 2 implementation of the AVMVN kernel is adapted from that of BEAST [42].

LeafAVMVN has the advantage of operating on all N leaf rates simultaneously (as well as learning their correlations), as opposed to ConstantDistance which operates on at most 2, or Scale which operates on at most 1 leaf rate at a time. As the size of the covariance matrix Σ_N grows with the number of taxa N, LeafAVMVN is likely to be less efficient with larger taxon sets. Therefore, the weight behind this operator is learned by AdaptiveOperatorSampler.

To prevent the learned weight behind LeafAVMVN from dominating the AdaptiveOperatorSampler weighting scheme and therefore inhibiting the mixing of internal node rates, we introduce the AdaptiveOperatorSampler(leaf) and AdaptiveOperatorSampler(internal) meta-operators which operate exclusively on leaf node rates $$ {\overrightarrow{\mathcal{R}}}_{\text{leaf}}^{}$$ and internal node rates $$ {\overrightarrow{\mathcal{R}}}_{\text{int}}^{}$$ respectively (Table 2). The former employs the LeafAVMVN operator and learns its weight during MCMC (after providing it sufficient time to learn Σ_N).

Model specification and MCMC settings

In all phylogenetic analyses presented here, we use a Yule [45] tree prior $$ p\left(\mathcal{T}\right|\lambda )$$ with birth rate λ ∼ Log-normal(1, 1.25). Here and throughout the article, a Log-normal(a, b) distribution is parameterised such that a and b are the mean and standard deviation in log-space. The clock standard deviation has a σ ∼ Gamma(0.5396, 0.3819) prior. Datasets are partitioned into subsequences, where each partition is associated with a distinct HKY substitution model [46]. The transition-transversion ratio κ ∼ Log-normal(1, 1.25), the four nucleotide frequencies (f_A, f_C, f_G, f_T)∼Dirichlet(10, 10, 10, 10), and the relative clock rate μ_C ∼ Log-normal(−0.18, 0.6) are estimated independently for each partition. The operator scheme ensures that the clock rates μ_C have a mean of 1 across all partitions. This avoids non-identifiability with branch substitution rates. To enable rapid benchmarking of larger datasets we use BEAGLE for high-performance tree likelihood calculations [47] and coupled MCMC with four chains for efficient mixing [29]. The neighbour joining tree [48] is used as the initial state in each MCMC chain.

Throughout the article, we have introduced four new operators. These are summarised in Table 3.

Table 3

Summary of clock model operators introduced throughout this article.

Operator	Description	Parameters
AdaptiveOperatorSampler	Samples sub-operators proportionally to their weights, which are learned (see Adaptive operator weighting).	$$
SampleFromPrior	Resamples a random number of elements from their prior (see Adaptive operator weighting).	$$
NarrowExchangeRate	Moves a branch and recomputes branch rates so that their genetic distances are constant (see Narrow exchange rate).	$$
LeafAVMVN	Proposes new rates for all leaves in one move (see An adaptive leaf rate operator) [30].	$$

Pre-existing clock model operators are summarised in Table 1. t denotes node heights while $$ \mathcal{T}$$ denotes the whole tree, including topology.

In Table 4, we define all operator configurations which are benchmarked throughout Results.

Table 4

Operator configurations and the substitution rate parameterisations which each operator is applicable to.

Configuration	Operator	Weight	real	cat	quant
nocons	RandomWalk($$)	10	✓	✓
	Scale($$)	10	✓
	Uniform($$)	10		✓	✓
	Interval($$)	10			✓
	Swap($$)	10	✓	✓	✓
	Scale(σ)	10	✓	✓	✓
cons	ConstantDistance($$)	$$	✓		✓
	SimpleDistance($$)	$$	✓		✓
	SmallPulley($$)	$$	✓		✓
	RandomWalk($$)	5	✓
	Scale($$)	2.5	✓
	Uniform($$)	5			✓
	Interval($$)	2.5			✓
	Swap($$)	2.5	✓		✓
	CisScale($$)	10	✓
	Scale(σ)	10			✓
adapt	AdaptiveOperatorSampler(σ)	10	✓	✓	✓
	AdaptiveOperatorSampler($$)	$$	✓		✓
	AdaptiveOperatorSampler($$)	30		✓
	AdaptiveOperatorSampler(root)	$$	✓		✓
AVMVN	AdaptiveOperatorSampler(σ)	10	✓		✓
	AdaptiveOperatorSampler(leaf)	$$	✓		✓
	AdaptiveOperatorSampler(internal)	$$	✓		✓
	AdaptiveOperatorSampler(root)	$$	✓		✓
NER{}	NarrowExchange	15	✓	✓	✓
NER	AdaptiveOperatorSampler(NER)	15	✓		✓

Within each configuration (and substitution rate parameterisation), the weight behind $$ {\overrightarrow{\mathcal{R}}}^{}$$ sums to 30, the weight of σ is equal to 10, and the weight of NER is equal to 15. Operators which apply to specific node sets (root, internal, leaf, or all) are weighted according to leaf count N. The adaptive operators are further broken down in Table 2. All other operators (i.e. those which apply to which apply to other terms in the state such as the nucleotide substitution model) are held constant within each dataset.

Round 1: A simple operator-weight learning algorithm greatly improved performance

We compared the nocons, cons, and adapt operator configurations (Table 4). nocons contained all of the standard BEAST 2 operator configurations and weightings for real, cat, and quant. cons additionally contained (cons)tant distance operators and employed the same operator weighting scheme used previously [33] (real and quant only). Finally, the adapt configuration combined all of the above applicable operators, as well as the simple-but-bold SampleFromPrior operator, and learned the weights of each operator using the AdaptiveOperatorSampler.

This experiment revealed that nocons usually performed better than cons on smaller datasets (i.e. small L) while cons consistently performed better on larger datasets (Fig 8 and S1 Fig). This result is unsurprising (Fig 3). Furthermore, the adapt setup dramatically improved mixing for real by finding the right balance between cons and nocons. This yielded an ESS/hr (averaged across all 9 datasets) 95% faster than cons and 520% faster than nocons, with respect to leaf branch rates, and 620% and 190% faster for σ. Similar results were observed with quant. However, adapt neither helped nor harmed cat, suggesting that the default operator weighting scheme was sufficient.

Fig 8

Round 1: Benchmarking the AdaptiveOperatorSampler operator.

Top left, top right, bottom left: each plot compares the ESS/hr (±1 standard error) across two operator configurations. Bottom right: the effect of sequence length L on operator weights learned by AdaptiveOperatorSampler. Both sets of observations are fit by logistic regression models. The benchmark datasets are displayed in Table 5. The cat and quant settings are evaluated in S1 Fig.

This experiment also revealed that the standard Scale operator was preferred over CisScale for the real configuration. Averaged across all datasets, the learned weights behind these two operators were 0.47 and 0.03. This was due to the computationally demanding nature of CisScale which invokes the i-CDF function. In contrast, the performance of Scale and CisScale were more similar under the quant configuration and were weighted at 0.28 and 0.45. For both real and quant, proposals which altered quantiles, while leaving the rates constant (Scale and CisScale respectively), were preferred.

Overall, the AdaptiveOperatorSampler operator was included in all subsequent rounds in the tournament.

Round 2: The real parameterisation yielded the fastest mixing

In Round 1 we selected the best configuration for each of the three rate parameterisations described in Branch rate parameterisations, and in Round 2 we compared the three to each other. adapt (real) and adapt (quant) both employed constant distance tree operators [33] and both used the AdaptiveOperatorSampler operator to learn clock model operator weights. Clock model operators weights were also learned in the adapt (cat) configuration.

This experiment showed that the real parameterisation greatly outperformed cat on most datasets and most parameters (Fig 9). This disparity was strongest for long alignments. In the most extreme case, leaf substitution rates r and clock standard deviation σ both mixed around 50× faster on the 15.5 kb seed plant dataset (Ran et al. 2018 [49]) for real than they did for cat. The advantages in using constant distance operators would likely be even stronger for larger L. Furthermore, real outperformed quant on most datasets, but this was mostly due to the slow computational performance of quant compared with real, as opposed to differences in mixing prowess (Fig 10). Irrespective of mixing ability, the adapt (real) configuration had the best computational performance and generated samples 40% faster than adapt (cat) and 60% faster than adapt (quant).

Fig 9

Round 2: Benchmarking substitution rate parameterisations.

Top left, top right, bottom left: the adapt (real), adapt (cat), and adapt (quant) configurations were compared. Bottom right: comparison of the mean tip substitution rate ESS/hr as a function of alignment length L.

Fig 10

Comparison of runtimes across methodologies.

The computational time required for a setting to sample a single state is divided by that of the nocons (cat) configuration. The geometric mean under each configuration, averaged across all 9 datasets, is displayed as a horizontal bar.

Overall, we determined that real, and its associated operators, made the best parameterisation covered here and it proceeded to the following rounds of benchmarking. This outcome allowed Rounds 3, 4, and 5 to commence. If cat dominated in this round, then the operators benchmarked in Rounds 4 and 5 would not be applicable to the parameterisation and the Bactrian proposal kernel benchmarked in Round 3 would not be applicable to branch rate categories in its present form.

Round 3: Bactrian proposal kernels were around 15% more efficient than uniform kernels

We benchmarked the adapt (real) configuration with a) standard uniform proposal kernels, and b) Bactrian(0.95) kernels [31]. These kernels applied to all clock model operators. These results confirmed that the Bactrian kernel yields faster mixing than the standard uniform kernel (Fig 11). All relevant continuous parameters considered had an ESS/hr, averaged across the 9 datasets, between 15% and 20% faster compared with the standard uniform kernel. Although the Bactrian proposal made little-to-no difference to the caterpillar dataset (Kawahara and Rubinoff 2013 [54]), every other dataset did in fact benefit. Bactrian proposal kernels proceeded to round 4 of the relaxed clock model optimisation protocol.

Fig 11

Round 3: Benchmarking the Bactrian kernel.

The ESS/hr (±1 s.e.) under the Bactrian configuration, divided by that under the uniform kernel, is shown in the y-axis for each dataset and relevant parameter. Horizontal bars show the geometric mean under each parameter.

Round 4: NER operators outperformed on larger datasets

Our initial screening of the NarrowExchangeRate (NER) operators revealed that the $$ \text{NER}\{{\mathcal{D}}_{AE},{\mathcal{D}}_{BE},{\mathcal{D}}_{CE}\}$$ operator outperformed the standard NarrowExchange / NER{} operator about 25% of the time on simulated data, however it was also very sensitive to the dataset. Therefore we wrapped up the two operators (NER{} and $$ \text{NER}\{{\mathcal{D}}_{AE},{\mathcal{D}}_{BE},{\mathcal{D}}_{CE}\}$$) within an AdaptiveOperatorSampler operator so that the appropriate weights could be learned. In this round we benchmarked the Bactrian + adapt (real) setting with the adaptive NER operator (Table 4). The benchmark datasets are fairly non clock-like and therefore could potentially benefit from NER (Table 5).

Our experiments confirmed that $$ \text{NER}\{{\mathcal{D}}_{AE},{\mathcal{D}}_{BE},{\mathcal{D}}_{CE}\}$$ was indeed superior on larger datasets (where L > 5kb; Fig 12). While there was no significant difference in the ESS/hr of continuous parameters, $$ \text{NER}\{{\mathcal{D}}_{AE},{\mathcal{D}}_{BE},{\mathcal{D}}_{CE}\}$$ did have an acceptance rate 41% higher than that of the standard NarrowExchange operator in the most extreme case (the bony fish alignment by Broughton et al. [53]). The moderate variance in branch substitution rates ($$ \widehat{\sigma }=0.5$$), coupled with a long alignment (7kb), and high topological uncertainty (Fig 12) made this dataset the perfect target. Every acceptance of a branch rearrangement proposal yields a new topology and thus facilitates traversal of tree space.

Fig 12

Round 4: Benchmarking the NER operators.

Top: the learned weights (left) behind the two NER operators (NER{} and $$ \text{NER}\{{\mathcal{D}}_{AE},{\mathcal{D}}_{BE},{\mathcal{D}}_{CE}\}$$), and the relative difference between their acceptance rates α (right), are presented as functions of sequence length. Logistic and logarithmic regression models are shown, respectively. Bottom: maximum clade credibility tree of the bony fish dataset by Broughton et al. 2013 [53]. This alignment received the strongest boost from NER, likely due to its high topological uncertainty and branch rate variance. Branches are coloured by substitution rate, the y-axis shows time units, and internal nodes are labelled with posterior clade support. Tree visualised using UglyTrees [59].

In contrast, the standard NarrowExchange operator outperformed on smaller datasets. The new operator was not always helpful and sometimes it even hindered performance. Use of an adaptive operator (AdaptiveOperatorSampler) removes the burden from the user in making the decision of which operator to use. The AdaptiveOperatorSampler(NER) operator proceeded into the final round of the tournament.

Round 5: The AVMVN leaf rate operator was computationally demanding and improved mixing very slightly

We tested the applicability of the AVMVN kernel to leaf rate proposals. This operator exploits any correlations which exist between leaf branch substitution rates. To do this, we wrapped the LeafAVMVN operator within an AdaptiveOperatorSampler (Table 2). The two configurations compared here were a) adapt + Bactrian + NER (real) and b) AVMVN + NER + Bactrian + adapt (real) (Table 4).

These results showed that the AVMVN operator yielded slightly better mixing (around 6% faster) for the tree likelihood, the tree length, and the mean branch rate (Fig 13). However, it also produced slightly slower mixing for κ, reflecting the high computational costs associated with the LeafAVMVN operator (Fig 10). The learned weight of the LeafAVMVN operator was quite small (ranging from 1 to 8% across all datasets), again reflecting its costly nature, but also reinforcing the value in having an adaptive weight operator which penalises slow operators. The LeafAVMVN operator provided some, but not much, benefit in its current form.

Fig 13

Round 5: Benchmarking the LeafAVMVN operator.

See Fig 11 caption for figure notation.

Overall, we determined that the AVMVN operator configuration was the final winner of the tournament, however its performance benefits were minor and therefore the computational complexities introduced by the LeafAVMVN operator may not be worth the trouble.

Tournament conclusion

In conjunction with all settings which came before it, the tournament winner outperformed both the historical cat configuration [3] as well as the recently developed cons (real) scheme [33]. Averaged across all datasets, this configuration yielded a relaxed clock mixing rate between 1.2 and 13 times as fast as cat and between 1.8 and 7.8 times as fast as cons (real), depending on the parameter. For the largest dataset considered (seed plants by Ran et al. 2018 [49]), the new settings were up to 66 and 37 times as fast respectively. This is likely to be even more extreme for larger alignments.

Modern operator design

Adaptability and advanced proposal kernels, such as Bactrian kernels, are increasingly prevalent in MCMC operator design [60–63]. Adaptive operators undergo training to improve their efficiency over time [64]. In previous work, the conditional clade probabilities of neighbouring trees have served as the basis of adaptive tree operators [26, 27]. Proposal step sizes can be tuned during MCMC [38]. The mirror kernel learns a target distribution which acts as a “mirror image” of the current point [32]. The AVMVN operator learns correlations between numerical parameters in order to traverse the joint posterior distribution efficiently [30].

Here, we introduced an adaptive operator which learns the weights of other operators, by using a target function that rewards operators which bring about large changes to the state and penalises operators which exhibit poor computational runtime (Eq 11). We demonstrated how learning the operator weightings, on a dataset-by-dataset basis, can improve mixing by up to an order of magnitude. We also demonstrated the versatility of this operator by applying it to a variety of settings. Assigning operator weights is an important task in Bayesian MCMC inference and the use of such an operator can relieve some of the burden from the person making this decision. However, this operator is no silver bullet and it must be used in a way that maintains ergodicity within the search space [64].

We also found that a Bactrian proposal kernel quite reliably increased mixing efficiency by 15–20% (Fig 11). Similar observations were made by Yang et al. [31]. While this may only be a modest improvement, incorporation of the Bactrian kernels into pre-existing operators is a computationally straightforward task and we recommend implementing them in Bayesian MCMC software packages.

Traversing tree space

In this article we introduced the family of narrow exchange rate operators (Fig 5). These operators are built on top of the narrow exchange operator and are specifically designed for the (uncorrelated) relaxed clock model, by accounting for the correlation which exists between branch lengths and branch substitution rates. This family consists of 48 variants, each of which conserves a unique subset of genetic distances before and after the proposal. While most of these operators turned out to be worse than narrow exchange, a small subset were more efficient, but only on large datasets.

Lakner et al. 2008 categorised tree operators into two classes. “Branch-rearrangement” operators relocate a branch and thus alter tree topology. Members of this class include narrow exchange, nearest neighbour interchange, and subtree-prune-and-regraft [43]. Whereas “branch-length” operators propose branch lengths, but can potentially alter the tree topology as a side-effect. Such operators include subtree slide [65], LOCAL [66], and continuous change [67]. Lakner et al. 2008 observed that topological proposals made by the former class consistently outperformed topological changes invoked by the latter [68].

We hypothesised that the increased efficiency behind narrow exchange rate operators could facilitate proposing internal node heights in conjunction with branch rearrangements. This would enable the efficient exploration of both topology and branch length spaces with a single proposal. Unfortunately, by incorporating a random walk on the height of the node being relocated, the acceptance rate of the operator declined dramatically (Fig 6). This decline was greater when more genetic distances were conserved.

These findings support Lakner’s hypothesis. The design of operators which are able to efficiently traverse topological and branch length spaces simultaneously remains an open problem.

Larger datasets require smarter operators

As signal within the dataset becomes stronger, the posterior distribution becomes increasingly peaked (Fig 3). This change in the posterior topology necessitates the use of operators which exploit known correlations in the posterior density; operators such as AVMVN [30], constant distance [33], and the narrow exchange rate operators introduced in this article.

We have shown that while the latter two operators are efficient on large alignments, they are also quite frequently outperformed by simple random walk operators on small alignments. For instance, we found that constant distance operators outperformed standard operator configurations by up to two orders of magnitude on larger datasets but they were up to three times slower on smaller ones (Fig 8). Similarly, our narrow exchange rate operators were up to 40% more efficient on large datasets but up to 10% less efficient on smaller ones (Fig 12).

This emphasises the value in our adaptive operator weighting scheme, which can ensure that operator weights are suitable for the size of the alignment. Given the overwhelming availability of sequence data, high performance on large datasets is more important than ever.

Future outlook

Based on these experiments, we have corroborated the results of Zhang and Drummond 2020 [33] by showing that the relaxed clock rate parameterisation can be greatly superior to the categories parameterisation. However this is conditional on the use of operators. Standard scale operators are ineffective in the real rate space on large datasets and constant distance operators tend to be the same on smaller ones. As the categories configuration is the current relaxed clock model default in both BEAST as well as BEAST 2, we therefore recommend researchers make the transition to the real rate parameterisation and employ an operator setup similar to that described in this article.

However, there remain several avenues for improvement. First, the use of a Weibull clock model prior, as opposed to a log-normal, could further decrease calculation time due to its closed-form inverse cumulative distribution function. Second, the leaf rate AVMVN operator presented in this article could be improved by proposing rates for only a small subset of correlated branch rates per proposal, as opposed to all at once. This would mitigate the computational burden associated with a covariance matrix which grows with the tree size (An adaptive leaf rate operator). Finally, while our adaptive weighting operator can learn operator weights within parameter spaces, it is not clear how to best assign weights between spaces; the relative weighting between clock model operators and tree topology operators for instance.