Adaptive dating and fast proposals: Revisiting the phylogenetic relaxed clock model

Relaxed clock models enable estimation of molecular substitution rates across lineages and are widely used in phylogenetics for dating evolutionary divergence times. Under the (uncorrelated) relaxed clock model, tree branches are associated with molecular substitution rates which are independently and identically distributed. In this article we delved into the internal complexities of the relaxed clock model in order to develop efficient MCMC operators for Bayesian phylogenetic inference. We compared three substitution rate parameterisations, introduced an adaptive operator which learns the weights of other operators during MCMC, and we explored how relaxed clock model estimation can benefit from two cutting-edge proposal kernels: the AVMVN and Bactrian kernels. This work has produced an operator scheme that is up to 65 times more efficient at exploring continuous relaxed clock parameters compared with previous setups, depending on the dataset. Finally, we explored variants of the standard narrow exchange operator which are specifically designed for the relaxed clock model. In the most extreme case, this new operator traversed tree space 40% more efficiently than narrow exchange. The methodologies introduced are adaptive and highly effective on short as well as long alignments. The results are available via the open source optimised relaxed clock (ORC) package for BEAST 2 under a GNU licence (https://github.com/jordandouglas/ORC).

where v = q(n − 1) indexes quantile q into piece number v .Values from the underlying function F −1 are cached, enabling rapid computation.

Tree operators for rate quantiles
Zhang and Drummond 2020 introduced several tree operators for the real parameterisation -including ConstantDistance, SimpleDistance, and SmallPulley [1].In this appendix, these three operators are extended to the quant parameterisation.Following the notation presented in the main article, let t i be the time of node i, let 0 < q i < 1 be the rate quantile of node i, and let r i = F −1 (q i ) be the real rate of node i where F −1 is the linear approximation of the i-CDF.

Constant Distance
Let X be a uniformly-at-random sampled internal node on tree T .Let L and R be the left and right child of X , respectively, and let P be the parent of X .Under the quant parameterisation, the ConstantDistance operator works as follows: Step 1 .Propose a new height for t X : where Σ is drawn from a proposal transition distribution (Uniform or Bactrian), and s is a tunable step size.Ensure that max{t L , t R } < t X < t P , and if the constraint is broken then reject the proposal.
Step 2 .Recalculate q X as: This ensures that the genetic distance between X and P remains constant after the operation by enforcing the constraint: (4) Step 3 .Similarly, propose new rate quantiles for the two children C ∈ {L, R}: Ensure that 0 < q i < 1 for all proposed nodes i ∈ {X , L, R}, and if the constraint is broken then reject the proposal.This constraint can only be broken from numerical issues.
Step 4 .Finally, in order to calculate the Metropolis-Hastings-Green ratio, return the determinant of the Jacobian matrix: As J is triangular, its determinant |J| is equal to the product of diagonal elements: The derivatives D F and D F −1 are computed using numerical approximations for the first and last pieces, or as the gradient of the linear approximation for internal pieces.As its final step, the operator returns ln |J|.

Simple Distance
While ConstantDistance proposes internal node heights, SimpleDistance operates on the root.Let X be the root node and let L and R be its two children.
Step 1 .Propose a new height for t X : Ensure that max{t L , t R } < t X , and if the constraint is broken then reject the proposal.
Step 2 .Propose new rate quantiles for the two children C ∈ {L, R}: These proposals ensure that the genetic distance between X and its children C remain constant after the operation by enforcing the constraint: Ensure that 0 < q C < 1, and if the constraint is broken then reject the proposal.
Step 3 .Finally, in order to calculate the Metropolis-Hastings-Green ratio, return the determinant of the Jacobian matrix: As J is triangular, its determinant |J| is equal to the product of diagonal elements: As its final step, the operator returns ln |J|.

Small Pulley
Just like the previous operator, SmallPulley operates on the root.Let X be the root node and let L and R be its two children.However, unlike SimpleDistance, this operator alters the two genetic distances Step 1 .Propose new genetic distances for d L and d R : and if the constraint is broken then reject the proposal.
Step 2 .Propose new rate quantiles for the two children L and R: Step 3 .Return the determinant of the Jacobian matrix: As J is triangular/diagonal, its determinant |J| is equal to the product of diagonal elements: Thus, as its final step, the operator returns ln |J|.

CisScale operator
CisScale was originally introduced by Zhang and Drummond 2020 for the real parameterisation (therein named ucldstdevScaleOperator).Under the quant configuration, the CisScale operator works as follows.
Step 1 .Propose a new value for the relaxed clock standard deviation σ Step 2 .Recalculate all branch substitution rate quantiles q such that their rates r remain constant Step 3 .Return the log Hastings-Green ratio of this proposal.If Σ was drawn from a symmetric proposal kernel (such as the Bactrian distribution), this is equal to: where derivatives D F and D F −1 can be approximated using either the piecewise linear model or standard numerical libraries.

Narrow exchange rates
The NarrowExchangeRate operator is also compatible with rate quantiles.This operator behaves the same as presented in the main article however the Hastings-Green ratio requires further augmentation due to changes in dimension throughout the proposal.
Step 1 .Apply NarrowExchange to the current tree topology as described in the main article.This will return a Hastings ratio H due to the asymmetry of this proposal.
Step 2 .Compute the relevant branch rates r i for r ∈ {A, B, C, D} of the current state from their respective quantile parameters.
Step 3 .Propose new rates and node heights and compute the Hastings-Green ratio of the real-space component of the proposal (e.g.Algorithms 1-2 of the main article).
Step 4 .Transform the rates back into quantiles.
Step 5 .Compute the log Hastings-Green ratio of the interconversion between rates and quantiles.log |J q | = log F (q) + log F −1 (r ). (28) Step 6 .Return the total log Hastings-Green ratio of this proposal: log Operators whose proposal kernels are affected by the decision to use a Bactrian kernel, as opposed to a uniform kernel, are specified below.

Table 1 :
Proposal kernels q(x |x) of clock model operators.In each operator, Σ is drawn from either a Bactrian(m) or uniform distribution.The scale size s is tunable.ConstantDistance and SimpleDistance propose tree heights t.The Interval operator applies to rate quantiles and respects its domain i.e. 0 < x, x < 1.A second NER algorithm is presented below.This operator was rejected by the screening protocol on simulated data.Algorithm 1 The NER{D BC , D CE } operator.1: procedure proposal(t A , t B , t C , t D , t E , r A , r B , r C , r D ) B (t D −t B )+r D (t E −t D )+r D (t E −t D ) t D −t B C (t E −t C )−r D (t E −t D ) t D −t C |J| ← (t D −t B )(t E −t C ) (t D −t B )(t D −t C ) r r return (r A , r B , r C , r D , t D , |J|)