{'fraction_non_alphanumeric': 0.03250404024309744, 'fraction_numerical': 0.012997063710650309, 'mean_word_length': 4.664517792676637, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 3, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fidelity of Hyperbolic Space for Bayesian Phylogenetic Inference

Bayesian inference for phylogenetics is a gold standard for computing distributions of phylogenies. It faces the challenging problem of. moving throughout the high-dimensional space of trees. However, hyperbolic space offers a low dimensional representation of tree-like data. In this paper, we embed genomic sequences into hyperbolic space and perform hyperbolic Markov Chain Monte Carlo for Bayesian inference. The posterior probability is computed by decoding a neighbour joining tree from proposed embedding locations. We empirically demonstrate the fidelity of this method on eight data sets. The sampled posterior distribution recovers the splits and branch lengths to a high degree. We investigated the effects of curvature and embedding dimension on the Markov Chain's performance. Finally, we discuss the prospects for adapting this method to navigate tree space with gradients.

Bayesian phylogenetics seeks the posterior distribution over discrete tree topologies and continuous tree branch lengths and evolutionary model parameters given an alignment of nucleotide sequences. Computing the posterior is analytically intractable so it is approximated using Markov chain Monte Carlo (MCMC) Rannala, 1997, Larget andSimon, 1999]. It is a workhorse algorithm of Bayesian phylogenetics that draws samples from the posterior distribution by proposing new samples that are "nearby" the current state of the Markov chain. However, the space of phylogenetic trees is super-exponential in the number of topologies and navigating tree space is difficult Matsen, 2015, Harrington et al., 2021]. Most methods to propose a new tree topology are too simplistic, e.g. nearest neighbour interchange (NNI) and subtree prune and re-graft (SPR), and improvements on these by guided proposals comes at the cost of additional computation [Höhna and Drummond, 2012]. This all poses the dilemma of how to best represent and navigate the space of phylogenies.
Hyperbolic space offers quality embeddings of tree-like data in low dimensions. Hyperbolic representations have successfully clustered data hierarchically into trees [Chami et al., 2020, Gu et al., 2018, Monath et al., 2019. They optimise a carefully devised objective function to find an optimal embedding. Unlike those methods, the objective function for Bayesian phylogenetics is prescribed by the model of evolution and prior distribution. Nonetheless, embedding genomic sequences as points in hyperbolic space and working in the embedding space could provide a way to move through tree space with an improved notion of"locality" -small changes to the embedding locations of sequences produce small changes in their distances on the tree. The embedding may enable a more natural representation of sequence divergences, potentially enabling a MCMC sample to make more complex changes to the model (topology, branch lengths etc) in a single move Recent works reviewed by Iuchi et al. [2021] have used representations, or specifically hyperbolic representations to learn a maximum likelihood phylogeny [Wilson, 2021, Matsumoto et al., 2021. However, their Bayesian counterpart is missing.
Bayesian phylogenetic practitioners increasingly need scalable methods to deal with larger sets of sequences, as highlighted by the current global pandemic of SARS-CoV-2 [Ki and Terhorst, 2022]. A demonstration that MCMC works in hyperbolic space could open up a wealth of possibilities for working with embeddings for Bayesian phylogenetics. As opposed to machine learning techniques, our work decodes a tree and uses the phylogenetic posterior probability as a cost function for hyperbolic embeddings. It also generates a distribution of phylogenetic trees rather than a point estimate.
To assess whether hyperbolic embeddings can represent a posterior distribution of trees, we perform MCMC on trees decoded from an embedding of nucleotide sequences. Trees are decoded using the neighbour-joining (NJ) algorithm [Saitou and Nei, 1987] before computing their likelihood and prior probability for each MCMC generation. We implement this MCMC in a python package called Dodonaphy, which is available on GitHub at https://github.com/mattapow/dodonaphy.
The goal of this paper is to empirically demonstrate the fidelity of hyperbolic embeddings for Bayesian phylogenetics. We begin by detailing the concepts from Bayesian phylogenetics, hyperbolic embeddings of trees and MCMC that are necessary to present the method. Once devised, we investigate the fidelity of the method and illustrate the embedding landscape. We then quantify the effects of the curvature and dimension of the embedding space on the MCMC performance. We also explore how a proposal distribution in the embedding space transfers to tree space. Our last result is the algorithm complexity and run-times. Finally, we discuss the numerous future research possibilities that this method opens.

Bayesian Phylogenetics
In phylogenetics, the posterior probability is the probability of a tree and evolutionary model parameters given an alignment of nucleotide sequences ψ. We refer to a phylogenetic tree T which includes both the topology and branch lengths of the tree. The posterior probability of a tree is Dodonaphy uses a simple model of evolution, the Jukes-Cantor (JC69) model [Jukes and Cantor, 1969] to compute the likelihood given a tree p(ψ|T ) using the Felsenstein pruning algorithm [Felsenstein, 1973]. For the prior p(T ), we utilise a uniform distribution over tree topologies and a Gamma-Dirichlet prior probability of a tree's branch lengths, as previously suggested by Rannala et al. [2012]. The Gamma-Dirichlet prior assigns a Gamma distribution Γ(α, β) with shape α = 1 and rate β = 0.1 on the total tree length (sum of branch lengths), before dividing this into individual branches with a Dirichlet distribution Dir(α k = 1.0).
The state-of-art phylogenetic software MrBayes [Ronquist and Huelsenbeck, 2003] offers both of these models for the prior and likelihood, which allows us to directly compare its results with Dodonaphy. We use two golden runs of MrBayes with 10 9 iterations and consider these the ground truth posterior distribution.
We approximate the posterior distribution through MCMC sampling, which does not require computing the intractable probability of the data p(ψ). As detailed later, it only computes the likelihood and prior probabilities of a tree. However, it does require generating proposals, for which, we move to a hyperbolic embedding.

The Hyperboloid Model
One common model of Hyperbolic space is the upper sheet of a hyperboloid H d = {x ∈ R d+1 : x, x = −1} where the Lorentz inner product is defined as The length of the geodesic between two points on the sheet is given by the metric In this model, for a fixed distance d(x, y), a more negative curvature κ < 0 stretches how far apart the points are on the curved sheet, imparting more curvature between the points. Conversely, to maintain a fixed distance as κ → 0, points move towards the origin where the space is flatter. The first coordinate is determined by the last d coordinates according to x, x = −1 to stay on the sheet. This rearranges to give Thus we need only store the last d points in R d and project up exactly onto the hyperboloid by φ : R d → H d ⊂ R d+1 . Explicitly, we use: projecting the first coordinate and leave the rest unchanged. The inverse mapping omits the first coordinate:

Hyperbolic Embeddings
To initialise an embedding of the taxa, Dodonaphy begins by minimising the stress of the embedding to a distance matrix. The distance matrix may either derive from a provided starting tree or the genetic distances between the aligned sequences. The stress of an embedding of points x i ∈ H d to a given set of pairwise distances D ij is Minimising this stress is a challenging non-convex problem, however, Hydra+ is a recent algorithm that provides a fast approximate solution [Keller-Ressel and Nargang, 2020]. It begins at a solution to the much simpler corresponding strain minimisation problem: finding a set of locations x to minimise ǫ though an Eigen-decomposition. Then it minimises the stress through gradient-based optimisation. We reimplemented Hydra+ as a Python package available on GitHub at https://github.com/mattapow/hydraPlus. We embed the n labelled taxa as a set of points in a hyperbolic space X = {x i }, i = 1, 2, ..., n using hydra+. By representing taxa as points on the hyperboloid x i ∈ H d we can work directly with their hyperbolic distances d(x i , x j ).

Tree Decoding
Dodonaphy then employs neighbour-joining to decode a tree T = NJ(X) from a set of embedding locations X. This decodes both branch lengths and a tree topology at once. The NJ algorithm takes pairwise distances as input, which is efficiently computed in the hyperboloid model of hyperbolic space using linear algebra [Chowdhary and Kolda, 2018]. The main advantage of using NJ is its consistency: it correctly constructs a tree when the given distances fit on that tree. However, instead of the sequence distances, Dodonaphy utilises the hyperbolic distances between taxa for NJ. This allows Dodonaphy to move through tree space by moving the set of n tip vectors X through hyperbolic space H n×d .

Unique tree decoding
The neighbour-joining algorithm fits a unique unrooted tree provided that the pairwise distances satisfy the four-point condition: ∀w, x, y, z When this occurs for all pairs of points, the distances are said to be additive and a tree consistent with the given distances can be recovered by NJ.
Hyperbolic space is often chosen for embedding trees because there exists a bound δ > 0 on how much hyperbolic distances violate the four-point condition: Such spaces are called δ-hyperbolic. Furthermore, this bound δ can be scaled to become arbitrarily small by making the curvature κ more negative [Wilson, 2021].

MCMC
Our interest is forming posterior distributions of trees, which we achieve through MCMC sampling of the posterior distribution. MCMC is a classical method for computing the posterior distribution p(T |ψ) given some data ψ. Starting at an initial location in the posterior X, a new point X * is proposed from a proposal distribution g(X * |X) and accepted or rejected according to the Metropolis-Hastings algorithm. A new proposal state X * decodes a tree T * = NJ(X * ) which is accepted with probability min(1, α), where the acceptance ratio is Importantly, it does not depend on the intractable component of the posterior p(ψ). Choosing a symmetric proposal g(X * |X), such as a Normal distribution, makes the last term simply cancel out. After a period of burn-in, points generated from this simple procedure compose samples from the targeted posterior distribution.

Improved MCMC
We make use of two modifications to the original MCMC algorithm. First, to improve mixing, Metropolis Coupled MCMC can use multiple chains that can swap intermittently through an MCMC move [Geyer, 1991, Altekar et al., 2004. The advantage of this is the ability to heat chains to move out of local optima more freely and it is routinely used in MrBayes. This is achieved by raising the acceptance ratio to a power given by a temperature parameter α τ , 0 ≤ τ ≤ 1. We use four chains with temperatures τ i = 1/(1 + λi), where i = 0, 1, 2, 3 is the chain index and λ = 0.1 is fixed, as is done in MrBayes. This ensures that the first chain remains cold (τ 0 = 1) and correctly samples from the posterior distribution.
The second idea is to adapt the covariance matrix of proposals using the Robust Adaptive Metropolis (RAM) algorithm [Vihola, 2012]. The RAM algorithm is a variant of the Metropolis-Hastings algorithm that uses a covariance matrix to adapt the proposal distribution. During a warm-up phase, we tune the initially diagonal covariance matrix simply by scaling it to achieve the target acceptance rate of 0.234, see appendix A. After the warm-up phase up we employ the RAM algorithm for the duration of the MCMC. Simultaneously, it tunes the acceptance rate to ensure computational efficiency.

MCMC Proposals
To propose a new MCMC point we draw from a multivariate Gaussian centred over the current taxa positions. We centre the Gaussian at the locations of each taxa x i ∈ H d , by omitting the first (prescribed) dimension and concatenating all locations into a single vector X ∈ R n×d . Note the distinction from the set of points denoted X = {x i } with i = 1, 2, ..., n indexing the taxa. The covariance matrix is initialised to a scalar multiple of the identity Σ = ζI n×d , with ζ = 0.1, before being tuned by the MCMC algorithms. Explicitly, each proposal point is drawn as X * ∼ N ( X, Σ). The proposal points are then projected onto the hyperboloid exactly using Eq.(1).

Experiments
In summary, Dodonaphy embeds pairwise tip distances D with hydra+ and performs phylogenetic MCMC on a hyperboloid, algorithm 1.6. An example embedding is illustrated in figure 1. The MCMC is run with a Gamma-Dirichlet prior on the branch lengths and a tree likelihood under a JC69 model of evolution. The MCMC runs four parallel chains and performs ten MCMC chain swap moves every 10 3 generations. MCMC proposals are multivariate normal distributions for all tip vectors. An initial warm-up period of 10 4 iterations serves to tune the covariance matrix before switching to the RAM algorithm until reaching 10 6 iterations. Dodonaphy draws 10 4 evenly spaced tree samples throughout the simulation. Unless otherwise stated, simulations are run in three dimensions with curvature κ = −1.
We looked at eight real datasets (DS1-DS8) compiled in Lakner et al. [2008] to evaluate the fidelity of hyperbolic space for Bayesian phylogenetics. The datasets contain between n = 27 and n = 64 taxa of aligned nucleotide sequences and are commonly used for phylogenetic benchmarking [Whidden et al., 2020]. They are freely available from treebase [Vos et al., 2012] with their identifier listed in table 1. The sequence lengths vary from 378 to 2520. We compress them into 256 ≤ L ≤ 1252 unique site patterns by weighting each pattern. Algorithm 1 Hyperbolic MCMC algorithm in Dodonaphy.

Results
To compare posterior distributions of phylogenetic trees we consider the distributions of splits, branch lengths, total tree lengths and the posterior tree probability. First, we present these in detail for DS1. Then we illustrate Dodonaphy's capacity over a range of datasets, curvatures and embedding dimensions.
To summarise the frequency of splits appearing in the posterior, we employ the average standard deviation of split frequencies (ASDSF) [Lakner et al., 2008]. It is a common statistic used to compare the splits between tree distributions, for example as used by MrBayes. A threshold of ASDSF < 0.05 is commonly used to indicate that the split frequencies from two distributions closely match each other.

Embedding Fidelity
We test Dodonaphy's capacity for MCMC starting directly from sequence data. Dodonaphy constructs the neighbour-joining tree directly from the multiple sequence alignment. The pairwise evolutionary distances (number of substitutions per site) between sequences are fed into hydra+ to embed the taxa as hyperbolic points. In this section, Dodonaphy runs for 10 7 generations. Figure 2a) demonstrates that Dodonaphy can move from this starting state into a similar region of tree space as the golden runs of MrBayes. The starting probability for the joint distribution of the data and parameters for Dodonpahy is about −8266, whereas for MrBayes it initialises to −9822. The trace plot illustrates that the burn-in period appears completed and the MCMC does not get stuck in one state for too long, pointing to good mixing.
Dodonaphy explores the posterior splits with considerable fidelity, figure. 2b). It recovered every split appearing in the golden run with support above 10 −3 . Of the 64 splits from MrBayes that Dodonaphy misses, only one split was visited by either golden run more than once. Indeed, the ASDSF between Dodonaphy and the first golden run is 0.003, well below the threshold for equivalence.
The mean length of splits appearing in both the golden runs and Dodonaphy are compared in figure 2c). They generally match well, however, Dodonaphy tends to slightly overestimate branch lengths. Consequently, the mean tree length posterior is overestimated by Dodonaphy: 0.453 compared to 0.408, a ratio of 1.11. The mean variance of the total tree length matches to a high degree: a ratio of 2.830 × 10 −4 /2.656 × 10 −4 = 1.007. Compared to the golden run, this overestimation reduces the overall log posterior, figure 2a) inset.
Repeated runs of Dodonpahy yield similar results on the tree length, figure 2d). The support, mean location and shape of the tree length distribution are self-consistent.
For the following sections, the MCMC chains are initialised from the consensus tree from MrBayes and run for a shorter 10 6 generations.

Embedding Curvature
The embedding curvature is a freely chosen hyper-parameter of the embedding space. Figure 3 reveals how a wide range of curvatures are suitable for all datasets. Panel a) shows that, compared to a golden run, the ASDSF falls below the threshold of 0.05 for curvatures in the range −100 ≤ κ ≤ −1. Similarly, the relative difference in the median total posterior tree lengthl to a golden runl mb falls within 10% for −100 ≤ κ ≤ −1 range for all datasets.
As the space becomes flatter (κ → 0), the space becomes more Euclidean. For κ = 0, we replace the hyperbolic distance with the regular Euclidean metric. These embeddings were initialised by  feeding hydra+ a minuscule curvature of 10 −10 . Note that MCMC in Euclidean space differs from vanilla MCMC because the tips are still embedded. As the space becomes flatter, the variance of tree lengths is generally larger and the trees tend to be longer, but not monotonically. Euclidean embeddings had better tree length estimates than weakly hyperbolic embeddings. However, the ASDSF significantly worsens in this limit, indicating that the wrong trees were recovered in Euclidean space.
In the other extreme, when the curvature is below −100, the variance of tree lengths becomes smaller. MCMCs with low variance are consistent with being constrained to a local optimum in tree space. This is also evidenced by worsening ASDSFs in this limit, signifying that, in same data sets, the tree topologies of the full posterior are not fully explored.

Tree Length Continuity
The δ-hyperbolic property of hyperbolic space makes it suitable for distance-based phylogenetics. We build on work by Wilson [2021] to highlight the continuity of the decoded tree lengths. First, let δ d,κ denote the minimal such error δ in the four point condition in d dimensions with curvature κ, Eq. 2.
Theorem 1. [Wilson, 2021] For any d ≥ 2 and κ < 0, there exists δ > 0 such that H d with curvature κ is δ-hyperbolic. Furthermore, This scaling is significant for neighbour-joining because it leads to additive embedding distances. A distance matrix D is additive if there is a weighted tree such that the distances on the tree D T match the given distances D = D T . The distances are additive when the four point condition is satisfied δ = 0. When the input distances are additive, neighbour-joining decodes a consistent tree. That is, it decodes a unique tree with correct tree distances: D = D T .
Neighbour joining remains consistent when the distances are almost additive. For the following theorem, we need the l ∞ norm, written ||D − D T || ∞ , which is the maximum element-wise difference between D and D T . The radius of a method f is the maximum α in so that f (D) = T produces the same tree with edges E(T ).
Theorem 2. [Atteson, 1999] The l ∞ radius of neighbour joining is 1 2 . Corollary 3. In the limit κ → −∞, the length of a decoded tree l(NJ(X)) is continuous in the embedding locations X.
Proof. In the limit κ → −∞, the error becomes arbitrarily small δ d,κ → 0. In particular, there exists κ < 0 such that ||D − D T || ∞ < α min e∈E(T ) (l e ). When this occurs, theorem 2 ensures the correct topology is decoded and the distances are arbitrarily close to additive. Consequently, the hyperbolic distances match the tip-tip distances on the tree, even as the topology changes. Since the hyperbolic distances are continuous in the embedding locations, so are the tree lengths.
If κ is too large, NJ may not reconstruct a unique tree, making the MCMC dependent on the particular taxa selected to merge in the algorithm. This difficulty could partly explain the large variances observed in tree lengths when κ > −1.
In the next section, we explore how the continuity of decoded distances translates into the posterior landscape.

Posterior Landscape
To gain an intuition for the effect of curvature on the posterior landscape, we mapped it for one taxon's position. We selected a set of embedding locations from an MCMC run in two dimensions H 2 . We then fixed all taxa locations except one, for which we performed a grid search. At each embedding location of this node, an NJ tree is decoded and the joint probability of this tree and the data is recorded. To plot in two dimensions, we projected the locations on the hyperboloid onto R 2 by φ −1 . A heat map of the posterior in H d (d = 2) reveals that the posterior distribution is multifaceted and depends on the locations of the other nodes, figure 4a). It is possible to move between different topologies, however the posterior appears multi-modal and is only smooth piecewise, with distinct regions corresponding to different tree topologies. This correspondence is apparent when compared with the middle panels, which shows the symmetric difference between the decoded tree to the optimal tree. Coarsely, locations leading to a larger symmetric difference are further in the embedding space from the optimal location. Consistent with corollary 3, the tree length surface is smoother with lower curvature, figure 4cf). With more negative curvature κ = −1000, changes in the posterior surface are less sharp, figure 4d). There are fewer jumps in the posterior and the regions have more uniform boundaries between them, figure 4e). This reflects the more local nature of the very curved space; moving one node has a smaller and more local effect on the tree. Whereas distal nodes are almost unaffected by perturbing this node. Sampling from a Normal distribution about the node could match the posterior surface better with this low curvature.

Embedding Dimension
The dimension of the hyperbolic manifold is also not prescribed and may affect the quality of the MCMC. Figure 5 highlights that just three dimensions is optimal to perform Bayesian inference. Increasing beyond three dimensions does not significantly affect the splits (ASDSF). Similarly, the variance of the tree length is generally well estimated in three or more dimensions. Two dimensions appears insufficient on both fronts to capture the phylogenetic posterior. Higher dimensions appeared to produce longer trees. This effect appears linear in the embedding dimension and depends on the data set. DS1, DS6 and DS8 are most affected, whereas DS2, DS3, DS4 and DS7 are least affected. These groups correspond to the shortest and longest trees according to the median length from the golden runs respectively. We checked that similar results are obtained with a lower curvature (Fig. 5 bottom row).

Locality in Tree Space
As Dodonaphy proposes new hyperbolic embeddings, the decoded tree topology freely changes according to NJ : {x i } → T . To quantify this change we looked at the symmetric difference from the current tree state to the proposal tree for each generation. This gives an idea of how topologically different proposal moves are from the current tree. The symmetric difference measures how many splits are in one tree but not the other. Looking at the cold chain in one MCMC run with 10 6 proposals, most proposals make small or no changes to the tree topology, figure 6. These include proposals from one heated chain, which are more likely to be in distant parts of the embedding space and could lead to larger symmetric differences. Longer proposals occur less and less frequently, decaying roughly exponentially.
In addition to the topology one might expect that branch lengths within a topology could be updated far more efficiently with this approach than in classical phylogenetic MCMC where they are generally updated one-at-a-time. This may matter more as trees grow large and practitioners start stacking in other continuous model parameters, as is done for viral phylodynamics.

Algorithm Run Time
The three most time-consuming algorithms of Dodonaphy's MCMC are the likelihood calculation O(nL), neighbour-joining O(n 3 ) and the pairwise distance computation O(n 2 d) [Chowdhary and Kolda, 2018]. We run Dodonaphy on Dell PowerEdge R640 or R740 servers with dual Intel(R) Xeon(R) Gold 6240 CPUs running at 2.60GHz. For each run, we record the duration per chain at each generation t over a range of taxa sets n with varying unique site patterns L and various embedding dimensions, figure 7. Adding the three main algorithms together, the best fit with standard deviations is an 3 + bLn + cn 2 d with a = 5.086 ± 0.766 × 10 −05 ms, b = 2.529 ± 0.205 × 10 −04 ms and c = 6.884 ± 0.661 × 10 −04 ms. Assuming three dimensions and taxa with L = 1000 unique sites, this means that the majority of the time is spent on neighbour joining when there are more than 237 taxa. Overall, the expected performance complexity is between quadratic and cubic in the (finite) number of taxa. Improvements to the NJ algorithm would yield the most significant speed-up as the number of taxa grows. For example, there are fast NJ implementations suitable for hundreds of thousands of sequences [Wheeler, 2009]. Similarly, BEAGLE provides a fast implementation of the likelihood computation [Ayres et al., 2012[Ayres et al., , 2019. Nevertheless, asymptotically, NJ remains the algorithmic bottleneck of Dodonaphy with O(n 3 ).

Fidelity Of Embedded MCMC
The results indicate that Hyperbolic space offers fruitful embeddings for Bayesian phylogenetics. A standard MCMC algorithm can use hyperbolic space to navigate the posterior phylogenetic landscape and produce the posterior tree distribution. It can start from an embedding initialised by the NJ tree on the sequence alignment to arrive at the posterior. Once converged, it recovers similar splits frequencies, split lengths and tree lengths to the ground-truth.
As expected, hyperbolic embeddings performed significantly better than Euclidean embeddings. In Euclidean space, phylogenies can be embedded by taking the square root of their distances [de Vienne et al., 2011]. However, embedding n taxa could require as many as n − 1 dimensions [Layer and Rhodes, 2017]. This makes it less scalable and less suitable for the large datasets being produced by infectious disease surveillance programs.
Previous work on maximum likelihood phylogenetics via hyperbolic embeddings has shown that the fidelity of the embedding improves as the embedding curvature decreases [Wilson, 2021]. In line with this existing work but in a Bayesian context, we find that decreasing the curvature from zero improves the MCMC quality. This improvement saturates at about κ = −1, indicating that the error to the four-point condition δ is negligible compared to the branch lengths. Decreasing beyond this to a wide range of curvature [−100, −1] provided good results, indicating the suitability of hyperbolic space for phylogenetic embeddings.
Lower curvatures tended to produce smaller tree length variances than the golden run. Underestimated variances could signify that the MCMC chains are contained in local optima, which is consistent with the higher split deviation (ASDSF) from the golden run. This suggests that there is an optimal range of curvature for embedding phylogenies, balancing between competing factors. Decreasing the curvature makes the tree length and, to some degree, the posterior smoother. However, a very curved embedding space becomes too localised and it is difficult to make significant changes to the tree topology.
Overestimated tree length variances indicate the MCMC is proposing trees outside the credible set, and sampled trees are either too long or too short compared to the credible set. Generally this occurred for embeddings with an almost Euclidean curvature (close to zero), where the error on the four-point condition is larger. This could indicate that the proposal distribution is not matching the posterior landscape well. A "neighbouring" topology in the embedding space with low posterior probability could be over-sampled in these finite MCMC runs. However, with a more curved space (κ = −100) this did not occur as significantly for three or more dimensions.
Just three hyperbolic dimensions appears sufficient to continuously represent phylogenetic trees. Higher dimensions led to overestimated tree lengths, especially for short trees. This overestimation occurred regardless of the embedding curvature and appeared linear in the embedding dimension. This result highlights the suitability of low dimensional hyperbolic embeddings, regardless of the number of taxa.
One possible explanation of the tendency to overestimate branch lengths is the Jacobian involved in transforming proposal samples in hyperbolic space to samples in tree space, via neighbour joining. The Jacobian comes from each of a series of three transformations, starting from the sample space R d and ending in tree space. Points are sampled in R nd and 1) projected onto the hyperboloid H nd ⊂ R n(d+1) via φ, 2) then these points are transformed into distances in R ( n 2 ) ≥0 before 3) finally being transformed into a tree by neighbour joining. Closed forms for the Jacobians of the first two are given in appendix D.
However, the determinant of the Jacobian of the second transformation is necessarily zero because for each of the n 2 the distances, only 2(d + 1) of the inputs (the two embedding points) contribute to the distance. The derivative of each distance is zero with respect to the other (n − 2)(d + 1) embedding locations. As we discuss below, multiple embedding configurations may produce the same set of distances, yielding a non-invertible transformation, consistent with a null Jacobian determinant. Further, computing the Jacobian of neighbour joining is non-trivial because it requires the min function. These Jacobian adjustments are omitted from Dodonaphy's algorithm.
One remedy for this issue is to use a prior on the embedding locations rather than on trees. We demonstrate that using a standard multivariate Normal as a prior does alleviate this trend in appendix C. A drawback of the multivariate Normal prior is that it is an uncommon prior for phylogenetics, making it difficult to compare to state-of-art methods.

Tree Considerations
Embedding the internal nodes of a tree could lead to more degrees of freedom to explore tree space faster. This is simple for a fixed, given topology. However, connecting the internal nodes into a consistent bifurcating tree is a challenging problem.
Multi-furcating trees could be decoded using a (minimum) spanning tree or the approach used by SLANTIS [Koptagel et al., 2022], which connects internal nodes using Bernoulli trials. The efforts of Friedman et al. [2001] yielded a way to convert a multi-furcating tree into a bifurcating tree with the same likelihood. One could make use of this algorithm after forming a minimum spanning tree of the embedded points. However, the method does not maintain consistency: distances in hyperbolic space between internal nodes wouldn't translate to distances on the tree.
Dodonaphy constructs unrooted NJ trees, which is useful for exploring the relationship between taxa. If a rooted time tree is desired, other clustering algorithms such as the UPGMA or WPGMA would produce rooted trees from the taxa embeddings. However, these algorithms require ultrametric distances to output a unique rooted tree, a property that hyperbolic space doesn't guarantee. If the distances are ultra-metric, one could also add an MCMC parameter for the rate of evolution to map from evolutionary distance to time. In this case, distance on the hyperboloid would represent time rather than substitutions per site.

Beyond JC69
The JC69 model has no free parameters, so on a given topology, the branch lengths alone determine the model likelihood. Substitution models with free parameters, such as the generalised time reversible substitution model [Tavare, 1986], are not directly encoded by the embedding. Incorporating these additional model parameters could be done as for a regular MCMC, rather than in an embedding space.

Embedding Isometries
Isometries of hyperbolic space can move an embedding to numerous locations and produce the same tree. These form an equivalence class on trees [T ]. The isometries of hyperbolic space are given by the orthochronous Lorentz group O + (1, d). These are characterised by the d + 1 × d + 1 orthogonal matrices g such that g T Hg = H.
Rotating all the tip locations would leave the tip distances unaltered and decoded to the same tree. In addition to the symmetries of hyperbolic space, the tree structure may possess symmetries. For example, swapping the position of two taxa in a cherry could leave the decoded tree unaltered. Accounting for symmetries in hyperbolic space and a decoded tree determines the equivalence classes of tree embeddings and is, to our knowledge, an open question.

Local Tree Moves
When the curvature is sufficiently small, neighbour joining decodes trees where the taxon-taxon distances in the decoded tree are arbitrarily close to the embedding distances, corollary 3. As a result, small perturbations of the locations of embedded sequences produce small changes to the taxon-taxon distances in the decoded tree. This provides a natural way to move "locally" in treespace: comparing the taxon-taxon distances, rather than the discrete number of NNI moves.
Dodonaphy proposes new MCMC states by continuous probability distributions. Here it uses a multivariate Normal distribution for proposals, however other distributions could equally be used.
We discuss an alternative proposal with similar results in appendix B. Matching the proposal distribution to the posterior is sought for the efficiency of MCMC chains. Exploring additional proposal distributions could lead to a better match.
This representation of the tree space opens up the possibility of exploring the neighbourhood of trees. Dodonaphy's representation of tree space bypasses classical topology changes such as NNI and SPR. Relating new moves found by Dodonaphy to existing moves (NNI, SPR) could improve methods that need to propose trees, such as MCMC.

Adding Taxa To an Existing Phylogeny
Online viral phylogenetics takes genetic sequences as they arrive, i.e. sequentially. Including these new sequences into trees in the posterior space is a matter of active research [Matsen et al., 2010, Fourment et al., 2018. Dodonaphy could offer the ability to place a new taxon onto a tree in constant time via the embedding space.
One simple way to approximate the embedding location z ′ of a new taxon labelled k is to use its distances to a subset of the taxa. Select a subset Z ⊆ X of taxa with cardinality one more than the embedding dimension |Z| = d + 1. Compose the matrix S ij = z i containing the taxa locations on the hyperboloid z i ∈ H d ⊂ R d+1 . Then, the genetic distances from the new taxa k to the selected subset are D i,k ∈ R d+1 ≥0 for i = 0, 1, ..., d + 1. These distances satisfy which has solution z ′ = S −1 cosh(D) in R d+1 . This solution could be mapped onto the hyperboloid sheet by selecting the closest point, or simply adjusting the first coordinate using eq. 1. Incorporating new taxa in this way could lead to a fast method of adding taxa to a phylogeny.
For more accuracy, S could be extended with additional taxa to become an over-determined system and solved approximately by its least squares solution z ′ = S T S −1 S T cosh(D). Alternatively, the subset of taxa Z could be reselected multiple times, effectively working as mini-batches. Then the most likely placement could be used based on the reconstructed super tree's probability.

Differentiable Embeddings
An exciting avenue for future research is the potential for gradient-based phylogenetic methods in the embedding space. Differentiating the posterior probability as a function of the embedding locations would allow this. However, gradients cannot flow through tree decoding by NJ as it is not differentiable. The difficulty lies in selecting which two taxa to connect during NJ, which requires using the (non-differentiable) minimum function.
Decoding a differentiable tree would provide a gradient-based framework for exploring tree space. This would be useful for point estimates of the maximum a posteriori tree or the maximum likelihood tree. It could also allow gradient-based Bayesian methods such as Hamiltonian Monte Carlo (HMC) and variational inference, which we discuss below.
Hamiltonian Monte Carlo assigns a mass to each tip position in the embedding space and draws on Hamiltonian mechanics to cover more of the posterior space in fewer iterations. For tree space, moving through the discrete barrier between topologies with momentum is still challenging. Dinh et al. [2017] propose a high dimensional d = (2n − 3)!! solution using a surrogate function to smoothly bridge the potential energy between topologies. Dodonaphy offers a low dimensional embedding space for HMC. Whilst the posterior may only be piecewise smooth in the embedding space, some methods cope with these transitions [Mohasel Afshar andDomke, 2015, Dinh et al., 2017].
Variational inference approximates the posterior distribution with much simpler analytical distributions. The goal is to minimise the error of the approximation via optimisation. Trees are discrete objects so this is usually considered a categorical distribution. For example, recent works appeal to graphical models to encode such categorical distributions Matsen, 2018, 2019]. In contrast, continuous distributions for the location of each node in the embedding space provide continuous distributions in tree space. Making use of distributions in hyperbolic space could open up an alternative technique for gradient-based variational inference.

Conclusions
Hyperbolic space offers a fruitful embedding space for Bayesian phylogenetics. Performing MCMC on eight data sets captured the splits (ASDSF < 0.05) and median tree lengths within 20% on all datasets. Embedding in three dimensions with a curvature between −100 ≤ κ ≤ −1 is sufficient to attain these results. The MCMC can start from an embedding of a NJ tree on the sequence data and move throughout hyperbolic space to explore the posterior space of trees.
Phylogenetic embeddings allow MCMC algorithms to propose new states (topology and branch lengths) from continuous probability distributions, which provides a novel way to make "local" changes in tree space.

Funding
This work was supported by the Australian Government through the Australian Research Council (project number LP180100593).

Acknowledgements
Computational facilities were provided by the UTS eResearch High Performance Computer Cluster.

Conflict of Interest
All authors declare that they have no conflicts of interest.

A A: Tuning MCMC
A simple method for tuning the covariance matrix in MCMC to achieve a target acceptance rate a * can be achieved as follows. Given an initial covariance matrix Σ, we scale it by a factor s that is tuned. Assume that increasing the step size decreases the acceptance rate and vice versa. We solve the following ordinary differential equation for the step size s: using the Euler method: s n+1 = s n + η(a n − a * ).
The learning rate η = (1 + n) −λ decays to zero. This ensures that as n → ∞ the step size is constant and the target distribution is sampled. We make the simple choice of the small multiple ζ = 0.1 of identity matrix as the initial covariance and choose λ = 0.5. This method is similar to the adaptive scaling metropolis (ASM) algorithm [Atchadé and Rosenthal, 2005] with a small difference. In Eq. 3, replacing the acceptance rate over all past states a n by the acceptance probability of the current proposal α reproduces ASM. We find that using a n leads to better performance for the warm up phase, giving faster convergence and an acceptance rate closer to the target acceptance of 0.234. B B: Wrapping onto the Hyperboloid Nagano et al. [2019] present a method to wrap Euclidean vectors onto the hyperboloid. This allows Dodonaphy to sample vectors drawn from arbitrary distributions in hyperbolic space. We provide a brief outline of this method and use it for proposals within Dodonaphy to compare results.
The method has two steps, parallel transport of vectors tangent to the hyperboloid and then exponentially mapping from the tangent space onto the hyperboloid. The tangent space of a point µ on the hyperboloid is the set of vectors on the tangent plane at µ. These are the points orthogonal to µ in R d+1 : T µ H d := {x ∈ R d+1 : x, µ = 0}.
Parallel transport takes a vector in one tangent plane x ∈ T ν H d and places it in another tangent plane T µ H d whilst conserving the vector's metric. It is given by: PT ν→µ (x) = x + µ − αν, x α + 1 (ν + µ) where α = − ν, µ Then, the exponential map wraps this vector in the tangent space onto the hyperboloid. It is constructed to preserve the vector's norm ||x|| L = x, x as follows exp ν (x) = cosh(||x|| L )ν + x ||x|| L sinh(||x|| L ).
As an example, Dodonaphy samples Euclidean vectors from a Normal proposal distribution and wraps them onto the hyperboloid. This is instead of projecting the first coordinate of the hyperbolic vector. Like before, MCMC is run for 10 7 generations and all else is held equal with κ = −1, d = 3.
As expected from MCMC, the results of using these proposals is reasonably similar to projecting proposal vectors via with φ. Both the split frequencies (ASDSF of 0.004) and mean lengths are captured well, figure 8. The similar results of this method demonstrate the possibility of using alternative proposals.

C C: Normal Prior on Locations
Placing a prior on embedding locations, rather than tree samples avoids the Jacobian involved in MCMC for transforming from embeddings to trees. Here, we use a standard multivariate normal distribution N ( 0, Σ) with covariance given by the d dimensional identity matrix. As before, Dodonaphy runs for 10 6 iterations including a 10 4 warm-up period using the tuning method presented above. It starts from the distances on consensus tree from the golden run of MrBayes.
As expected, the tree length estimates vary from those given by the golden run, merely because they have different priors, figure 9. More to the point, it demonstrates that there is no clear trend in trees becoming longer with higher dimensions. This suggests that when a Gamma-Dirichlet prior on trees is used, the trend in increasing tree length (in the main text) derives from the non-trivial Jacobian involved in the MCMC.

D D: Closed form of Jacobians
The Jacobian of the first projection φ : R d → R d+1 is non-square because it attaches an extra dimension to each taxa. For a taxon with stored location x ∈ R d , it takes the form For the second transformation from the hyperboloid to the pairwise distances, the Jacobian the distance between two points x, y ∈ H d is adapted from Chowdhary and Kolda [2018]: ∂d(x, y, κ) ∂x = 1 −κ((x * y) 2 − 1) ||y|| 2 2 + 1 ||x|| 2 2 + 1 x − y .
It is clear that the distance between x and y does not depend on a third point z: ∂d(x, y, κ) ∂z = 0.