Directed Acyclic Graph Expansion Navigates High Posterior Density Phylogenetic Tree Space

A systematic search framework designed to explore phylogenetic tree space by expanding a subsplit directed acyclic graph can successfully identify high posterior density evolutionary trees, yet aggregating trees from standard Markov chain Monte Carlo sampling remains computationally faster and captures a higher proportion of probable topologies. In Algorithms for Molecular Biology, researchers investigated whether graph expansion could replace traditional random-walk algorithms in Bayesian phylogenetics. Bayesian phylogenetics conventionally estimates tree posterior distributions using Markov chain Monte Carlo methods. These algorithms initialize a starting tree, generated randomly or via parsimony, and propose successive stochastic modifications accepted or rejected according to the Metropolis-Hastings ratio, biasing sampling toward higher likelihood and prior probability. However, because high-confidence topologies occupy a tiny fraction of a super-exponential search space, most modifications produce significantly worse trees and are discarded. This dynamic produces low acceptance rates, inherently constraining the efficiency of random walks across complex distributions.

In contrast to stochastic sampling, maximum likelihood approaches navigate phylogenetic space systematically by iteratively applying local rearrangements to improve the likelihood of the current top-scoring tree. While maximum likelihood methods operate faster than Markov chain Monte Carlo samplers, they focus on locating optimal trees rather than characterizing credible sets of plausible topologies. Previous attempts to bridge this divide systematically explored tree space in parallel by applying local rearrangements to visited topologies and retaining candidates exceeding a likelihood threshold. Normalizing the likelihoods of these topologies approximated a posterior distribution, but the strategy was defeated because the sheer volume of high-quality trees overwhelmed branch length optimization routines. To circumvent this bottleneck, recent work introduced the subsplit directed acyclic graph, or sDAG, which provides a compact mathematical structure capable of storing, manipulating, and evaluating vast collections of phylogenetic trees simultaneously without examining each individual topology in isolation.

Within the sDAG framework, each node represents a subsplit—a bipartition of a subset of taxa—whose partition components are termed subsplit-clades. Directed edges run leafward from the root, connecting parent subsplits to child subsplits that further divide constituent clades. Any bifurcating tree within the sDAG represents a valid topology, and nodes whose subsplit-clades possess an outdegree greater than one require a choice among descending paths to determine branching orders. Assigning probability distributions across edges emerging from each subsplit-clade yields a normalized probability distribution over all encoded topologies. This formulation strictly generalizes conditional clade distributions, or CCDs, because the sDAG allows splitting probabilities to depend directly on the sister clade rather than treating identical clades independently. When an sDAG is constructed from input topologies, it naturally encodes novel topological combinations beyond the initial inputs, combinatorially expanding support while risking the addition of non-credible topologies outside the target posterior distribution.

To systematically grow an sDAG, the researchers adapted the classical nearest neighbor interchange, or NNI, from an operation that alters a single topology into an operation that enlarges the graph itself. In traditional phylogenetics, an NNI swaps two subtopologies attached to adjacent nodes to generate an alternate tree. When applied to an sDAG containing subsplits t and s over clades X, Y, and Z, the interchange swaps clades Y and Z to construct new subsplit nodes t’ and s’, alongside a central edge connecting t’ to s’. Constructing the combined graph requires adding these new nodes alongside parent edges from ancestors of t to t’ and descendant edges from new subsplit-clades to child nodes of X, Y, and Z. Maintaining edges between all compatible subsplits guarantees that every newly introduced topology represents a genuine NNI modification of an existing topology and strictly contains the central edge.

Evaluating whether a candidate NNI should be incorporated into the sDAG requires formal optimality criteria that generalize tree likelihood calculations across multi-topology graphs. To maintain computational tractability, the researchers assigned a single fixed branch length to every edge within the sDAG, establishing an exact one-to-one correspondence between topologies and fully parameterized trees. The team adopted the Jukes-Cantor substitution model under site independence, computing the overall phylogenetic likelihood as the product of site likelihoods using Felsenstein’s pruning algorithm. Although Bayesian phylogenetic formulations typically integrate over continuous distributions of branch lengths and across-site rate variation, the researchers optimized fixed branch lengths using Brent optimization, an approximation strategy that preserves algorithmic efficiency while demonstrating reliable empirical performance. Based on these foundations, the investigators formulated two distinct objective criteria to guide systematic sDAG enlargement: top pruning and generalized pruning.

The first search method, designated as top pruning, aims to incorporate the highest-likelihood new tree introduced by an NNI without exhaustively evaluating all possible topologies generated by the operation. Because enumerating every tree containing the new central edge is computationally prohibitive, top pruning approximates the optimal tree using a recursive data structure termed a choice map. For each edge in the sDAG, a rootward choice map identifies a parent edge and sibling edge directed toward the root, while a leafward choice map identifies left and right child edges directed toward terminal taxa. Starting from any specified edge, applying these choice maps iteratively yields a deterministic tree called the best known tree. By caching partial likelihood vectors across graph nodes, this best known tree likelihood is calculated in constant time relative to sDAG leaf count.

In practice, the top pruning algorithm operates through a structured ten-step iterative cycle. The procedure begins by initializing an sDAG from a list of starting trees ordered by likelihood, assigning edge branch lengths from the first tree containing each edge, and configuring initial choice maps. Edges are subsequently introduced between all compatible subsplits, with choice maps assigned via an argmax strategy and branch lengths optimized for resulting best known trees. The algorithm compiles a ranked list of candidate NNIs that enlarge the graph—evaluating two possible interchanges per edge—ordered by top pruning likelihood. In each iteration, the algorithm applies the top NNI, introduces compatible edges, optimizes branch lengths of newly added edges using Brent optimization, and inserts new candidate interchanges into the sorted list. This cycle repeats for a designated number of iterations or until remaining interchanges fall below a likelihood threshold.

The alternative search strategy, termed generalized pruning or GP, evaluates candidate NNIs based on their impact on a topology-marginal composite likelihood rather than the likelihood of a single representative tree. For an interchange producing central edge e, generalized pruning computes the product across all alignment sites of the sum of site likelihoods over all trees containing e, weighted by a uniform prior distribution across those topologies. This per-edge composite marginal likelihood differs fundamentally from the true marginal likelihood: whereas the true marginal likelihood requires separate partial likelihood vectors for every combination of sequence site and topology, the composite approximation shares partial likelihood vectors across topologies at each node. This sharing enables generalized pruning calculations to scale linearly with the number of edges in the sDAG. Both candidate interchange evaluation and branch length optimization on newly added edges are governed by maximizing this composite metric.

The research team implemented both search algorithms within the bito C++ library, providing a Python interface to manage systematic graph operations, and benchmarked performance against empirical datasets executed via the SLURM job scheduler. Benchmarks utilized seven standard reference datasets, designated DS1 through DS8 (excluding DS2 due to negligible posterior diversity), alongside a 100-sequence influenza A dataset known as flu100. Reference standards were established using long MrBayes golden runs to delineate empirical 95% credible sets, defined as minimal topology collections accounting for at least 95% of cumulative posterior density. The datasets spanned a broad gradient of posterior diffuseness, ordered from least to most diffuse as DS3, DS1, DS4, DS7, DS8, DS6, and DS5. Performance was measured by tracking cumulative posterior density captured per search iteration and against clock run time, comparing algorithms to short MrBayes runs evaluated both as raw topology lists and as aggregated sDAGs.

Across all benchmarked datasets, generalized pruning exhibited poor performance and failed to capture meaningful posterior density within practical timeframes. Even on DS3, the least diffuse dataset, generalized pruning consumed numerous iterations adding edges and subsplits that contributed nothing to the credible set, proving uncompetitive even against raw, unaggregated MCMC topology lists. In contrast, top pruning demonstrated substantially greater efficacy. It consistently captured higher cumulative posterior density than raw MCMC lists and performed competitively against sDAGs constructed by aggregating short MCMC runs, successfully navigating through prolonged plateaus to recover primary posterior density regions. However, evaluating the credibility of subsplits introduced during search operations revealed clear limitations: on diffuse datasets such as DS5 and DS6, top pruning incorporated a substantial fraction of non-credible subsplits alongside credible ones, whereas short MCMC runs maintained markedly higher true-positive subsplit discovery rates.

Seeding top pruning with distinct topologies from 200 maximum likelihood RAxML runs improved early interchange selections by exploring separate distributional peaks, whereas seeding with the top ten credible trees offered negligible benefit. Runtime profiling confirmed that top pruning operates at worst quadratically rather than linearly because later iterations introduce progressively more compatible edges requiring likelihood evaluations; nevertheless, consistent results on the 100-sequence flu100 dataset indicated stable scaling across moderate taxa numbers. Across datasets, the likelihood of an edge’s best known tree correlated strongly with its empirical posterior probability, demonstrating accurate ranking within credible regions. The researchers also examined a variant restricting additions to five edges per interchange, but demonstrated that this approach leaves certain valid topologies permanently unreachable. Ultimately, the authors concluded that aggregating trees from short MCMC runs into an sDAG remains the faster, superior strategy for recovering high posterior density regions prior to downstream variational inference.

Subject of Research: Systematic exploration of Bayesian phylogenetic tree space using subsplit directed acyclic graphs

Article Title: Finding high posterior density phylogenies by systematically extending a directed acyclic graph

Article References: Jennings-Shaffer, C., Rich, D. H., Macaulay, M., Karcher, M. D., Ganapathy, T., Kiami, S., Kooperberg, A., Zhang, C., Suchard, M. A., & Matsen, F. A., IV (2025). Finding high posterior density phylogenies by systematically extending a directed acyclic graph. Algorithms for Molecular Biology, 20(1), Article 2. https://doi.org/10.1186/s13015-025-00273-x

Image Credits: AI Generated

DOI: 10.1186/s13015-025-00273-x

Keywords: Bayesian phylogenetics, subsplit directed acyclic graph, Markov chain Monte Carlo, phylogenetic tree space, nearest neighbor interchange, top pruning, generalized pruning, computational biology, Finding, high, posterior, density