Branch-Specific Substitution Models

Updated 6 July 2026

Branch-specific substitution models are evolutionary frameworks that allow the substitution process to vary among different branches, capturing lineage-specific rates, compositions, and selection pressures.
They employ various mathematical formulations—such as non-homogeneous Lie Markov models and Markov-modulated methods—to address limitations of homogeneous models and ensure mathematical consistency under heterogeneous conditions.
These models enhance phylogenetic inference by providing more accurate branch length estimates and root inferences, with advances in computational techniques like GPU acceleration and analytic gradient methods.

Branch-specific substitution models are evolutionary models in which the substitution process is allowed to differ among branches or lineages of a tree, rather than being forced to be homogeneous over the entire phylogeny. In contemporary phylogenetics, this phrase most commonly refers to continuous-time Markov chain frameworks that relax the use of a single global rate matrix $Q$ , either by assigning branch-local matrices, by allowing latent regime switching along lineages, or by imposing structured non-homogeneity through multiplicatively closed model families. These models are motivated by compositional heterogeneity, heterotachy, lineage-specific selection, and the failure of homogeneous reversible models to infer roots or to remain mathematically consistent under heterogeneous composition of processes across branches (Baele et al., 2019, Woodhams et al., 2014, Hannaford et al., 2020, Ji et al., 11 Jul 2025).

1. Concept and motivation

The classical phylogenetic CTMC framework assumes that a single rate matrix governs substitution throughout the tree, so that under time-homogeneity $P(t)=\exp(Qt)$ . Early models therefore treated the substitution process as homogeneous over time and across sites, with later extensions usually concentrating on among-site rate variation rather than lineage dependence (Baele et al., 2019). This standardization is computationally convenient, but it can be biologically restrictive when different branches experience different base compositions, exchangeabilities, or overall rates, or when a site changes evolutionary regime over time.

The principal motivation for branch-specific substitution modeling is that simple ASRV changes only the overall rate at each site, not the substitution process itself. The literature summarized here emphasizes that this distinction matters when evolutionary change is qualitative rather than merely scalar, as in shifts in GC content, codon-position-specific dynamics, compositional drift, heterotachy, or lineage-specific selection (Baele et al., 2019). A standard homogeneous model can then misestimate branch lengths, over-support the wrong tree, or fail to recover biologically meaningful regime changes.

A second motivation is inferential rather than descriptive. Standard stationary reversible models yield a likelihood that is invariant to root position, so they cannot infer the root from likelihood alone. Branch-specific non-stationary and locally non-reversible models can make the likelihood root-sensitive by allowing composition and time direction to vary across branches (Hannaford et al., 2020).

A third motivation is mathematical consistency. If a lineage evolves under one rate matrix $Q_1$ for time $t_1$ and then under another $Q_2$ for time $t_2$ , the overall transition matrix is $e^{Q_1 t_1}e^{Q_2 t_2}$ . For many popular model families, especially GTR, that product does not remain in the same family. This means that heterogeneous evolution along a phylogeny cannot in general be summarized consistently by a single model class, even when every segment individually belongs to that class (Woodhams et al., 2014).

2. Principal mathematical formulations

Several distinct mathematical constructions instantiate branch-specificity.

Framework	Mechanism of branch specificity	Distinguishing property
Non-homogeneous Lie Markov models	Each branch $b$ has its own $\mathbf{Q}_b$ from one Lie Markov family	Root-sensitive, non-stationary step changes at speciation events
Markov-modulated models	Hidden regime switches along lineages in compound states $(k,s)$	Branch dependence through latent regime trajectories
General Markov models on arbitrary splits	Split-specific operators $P(t)=\exp(Qt)$ 0 on tensor product state spaces	Extends from trees to incompatible splits and networks

In the non-homogeneous Lie Markov construction, a rooted binary tree with $P(t)=\exp(Qt)$ 1 taxa has $P(t)=\exp(Qt)$ 2 branches, and each branch $P(t)=\exp(Qt)$ 3 receives its own instantaneous rate matrix $P(t)=\exp(Qt)$ 4. The models studied explicitly use non-reversible Lie Markov families, especially RY5.6b and RY8.8a. If the tree is rooted on branch $P(t)=\exp(Qt)$ 5, the root distribution is taken from the stationary distribution of $P(t)=\exp(Qt)$ 6, and that branch is split into two edges on either side of the root. The resulting model is non-homogeneous, non-stationary, and locally non-reversible, with step changes in composition at speciation events (Hannaford et al., 2020).

Markov-modulated continuous-time Markov chains implement branch-specificity differently. Rather than assigning a free matrix to each branch, they augment the observed character state $P(t)=\exp(Qt)$ 7 with a hidden regime index $P(t)=\exp(Qt)$ 8. With $P(t)=\exp(Qt)$ 9 component models and switching matrix $Q_1$ 0, the full generator is

$Q_1$ 1

The process evolves on compound states $Q_1$ 2, so substitutions may occur within a regime, or the regime may switch while the observed state is unchanged. Branch dependence is therefore encoded through latent switching trajectories along lineages. The paper explicitly notes that this is not the same as giving each branch its own unconstrained free matrix directly (Baele et al., 2019).

The arbitrary-split formulation of the general Markov model replaces edge-local reasoning with operators on tensor product state spaces. In the two-state case, the branching process is represented by the splitting operator $Q_1$ 3, and the key intertwining relation pushes lineage-specific generators through a split. This yields split-specific rate operators

$Q_1$ 4

for subsets $Q_1$ 5. Tree models are recovered when incompatible split weights are zero; network models arise when incompatible splits are assigned their own weights and possibly their own rates. In this formalism, branch-specific substitution modeling extends naturally to arbitrary split systems and can represent convergence of previously divergent lineages (Sumner et al., 2010).

A common misconception is that branch-specificity necessarily means a separate unconstrained process on every edge. The MMM literature makes clear that branch-specific dependence can instead be expressed through hidden regime switching, and the arbitrary-split algebra shows that the relevant object may be a split-specific operator rather than a conventional branch-local rate matrix (Baele et al., 2019, Sumner et al., 2010).

3. Closure, identifiability, and model families

The strongest formal argument for particular branch-specific model families comes from multiplicative closure. Lie Markov models are defined so that their rate matrices are closed under addition, scalar multiplication, and the commutator

$Q_1$ 6

These conditions make the model space a Lie algebra, and the corresponding Markov matrices are then locally closed under multiplication. Matrix closure implies phylogenetic closure, meaning that pruning and marginalization do not force the process out of the model family (Woodhams et al., 2014).

This criterion is important precisely because GTR is not closed under matrix multiplication in the required sense. Under branch heterogeneity, a homogeneous GTR fit can therefore be mathematically inconsistent, not merely biologically imperfect. The heterogeneous-process literature treats this failure as a central reason to prefer Lie Markov models when branch-specific composition of substitutions is the object of interest (Woodhams et al., 2014).

The Lie Markov hierarchy examined in detail consists of 37 RY Lie Markov models spanning dimensions 1 through 12, with names of the form $Q_1$ 7, and 99 distinct models when RY, WS, and MK variants are all counted. The eight-parameter model 8.8 is singled out because it is the closure of HKY, making it the smallest Lie Markov model containing all products of HKY matrices. Empirically, RY8.18 is reported as the best overall model by mean BIC ranking, while RY8.8 and RY8.18 are the top Lie Markov performers overall (Woodhams et al., 2014).

Within the branch-specific non-homogeneous framework, the two families used explicitly are RY5.6b and RY8.8a. RY5.6b has an additive redundancy, removed by constraining $Q_1$ 8 and fixing the trace to $Q_1$ 9, with $t_1$ 0. RY8.8a is more flexible, fixes the trace to $t_1$ 1, and uses $t_1$ 2. In the non-homogeneous RY5.6b model only $t_1$ 3 varies by branch while $t_1$ 4 is constant across the tree, whereas in the non-homogeneous RY8.8a model all parameters in $t_1$ 5 can vary by branch (Hannaford et al., 2020).

Identifiability problems also arise in the MMM setting. When MMMs are combined with ASRV, hidden rate categories and hidden substitution regimes can become indistinguishable. The recommended practical restriction is therefore

$t_1$ 6

when ASRV is also used. The MMM framework otherwise remains general enough to nest homogeneous CTMCs, ASRV, covarion models, and mixture-like models as limiting cases (Baele et al., 2019).

4. Likelihood evaluation and scalable inference

Likelihood computation in branch-specific models usually retains Felsenstein pruning, but on an enlarged or non-homogeneous state space. For MMMs, the partial likelihood at each node is indexed by compound states $t_1$ 7, transition probabilities are obtained from $t_1$ 8, and the main computational burden arises because $t_1$ 9 has dimension $Q_2$ 0. The reported cost scales roughly as $Q_2$ 1 for transition probability calculations plus $Q_2$ 2 for peeling across $Q_2$ 3 branches and $Q_2$ 4 sites (Baele et al., 2019).

To make such models practical, the BEAST implementation of MMMs exploits BEAGLE, which offloads matrix operations to optimized multi-core CPUs and especially GPUs. The bottlenecks are matrix exponentiation and multiplication, and the benchmarked implementation reports that GPUs reduce the runtime penalty of large MMMs far more effectively than CPUs, with the Tesla P100 giving the best performance in their benchmarks (Baele et al., 2019).

A distinct computational problem appears when one optimizes or samples a single branch length while the rest of the tree is fixed. For that setting, the one-dimensional likelihood $Q_2$ 5 can be replaced by the four-parameter surrogate lcfit, derived from the binary symmetric model and generalized by a truncation term. In its supplementary form,

$Q_2$ 6

The surrogate has a finite long-branch asymptote and flexible local shape, addressing properties of branch-length likelihoods that generic polynomials or standard PDFs do not match. Two fitting strategies are described: lcfit4, a nonlinear least-squares fit of all four parameters using Levenberg–Marquardt and then SLSQP if needed, and lcfit2, which uses the ML point and curvature directly after locating the true ML branch length by Brent’s method (Claywell et al., 2017).

For high-dimensional branch-specific parameter inference, an exact linear-time gradient algorithm has been developed for derivatives with respect to branch-specific substitution parameters. The key result is that, once post-order partial likelihood vectors $Q_2$ 7 and pre-order partial likelihood vectors $Q_2$ 8 are available, each gradient component can be written in a branch-local sandwich form,

$Q_2$ 9

Because all derivatives can be accumulated in one pre-order traversal after the usual post-order pass, total work is linear in the number of branch parameters (Ji et al., 11 Jul 2025).

5. Inference targets and empirical behavior

Branch-specific substitution models are used to estimate not only trees and branch lengths, but also lineage-dependent changes in composition, selective pressure, mutational mechanisms, and hidden evolutionary regimes. In the change-point setting, the branch-specific parameter on branch $t_2$ 0 is compared to its parent through

$t_2$ 1

with independent Bayesian bridge priors

$t_2$ 2

With $t_2$ 3, most parent-child increments are shrunk toward zero while large deviations survive, allowing change-points to be learned automatically without pre-specifying branch groups (Ji et al., 11 Jul 2025).

The empirical record in phylogenetics is heterogeneous but substantial. In bacterial 16S rRNA, MMMs that allow branch-specific changes in base composition recover the biologically expected tree grouping Deinococcus radiodurans with Thermus thermophilus, whereas standard HKY and GTR with ASRV incorrectly cluster D. radiodurans with Bacillus subtilis. In plastid genes psaB and ndhD, MMMs strongly outperform standard GTR+ASRV, alter topology and posterior support outside the seed plants, and reveal hidden-state patterns not explained by simple codon-position partitioning. In influenza HA, MMM(HKY) $t_2$ 4 with four rate categories and shared $t_2$ 5 outperforms ASRV models dramatically, even with codon-position partitioning, and yields more recent or changed root-height estimates (Baele et al., 2019).

Simulation studies argue against the view that high-dimensional branch-specific models are automatically overparameterized. Under simulated GTR data, the GTR model is correctly preferred and MMMs are not systematically favored merely because they are more parameter-rich. Under simulated MMM(GTR) $t_2$ 6 or MMM(GTR) $t_2$ 7, generalized stepping-stone sampling recovers the generating model, supporting the interpretation of large Bayes factors as meaningful rather than as artifacts of overfitting (Baele et al., 2019).

Branch-specific non-homogeneous Lie Markov models have been used for direct rooted-tree inference under compositional heterogeneity. In simulations, root inference improves with more sites, and the non-homogeneous RY8.8a model performs much better than RY5.6b. In a Drosophila analysis of 2085 sites from the Xdh gene of 17 species including an outgroup, the non-homogeneous RY8.8a model identifies a biologically credible root within the outgroup and has the best marginal likelihood, whereas homogeneous models and some non-homogeneous reversible models fail to find a plausible root (Hannaford et al., 2020).

The change-point literature emphasizes branch-specific codon and nucleotide parameters. In primate BRCA1, branch-specific $t_2$ 8 estimates reaffirm elevated $t_2$ 9 along the lineage leading to human and chimpanzee without specifying in advance which branches should share the same value. In 138 complete MPXV genomes from the 2022 mpox epidemic, an HKY+APOBEC branch-specific model reveals a clear acceleration of the APOBEC effect after a particular internal branch, with an estimated split time centered around 2012.3 and a 95% posterior interval of $e^{Q_1 t_1}e^{Q_2 t_2}$ 0 (Ji et al., 11 Jul 2025).

Computationally, these methods have reported substantial gains. Relative to central-difference numerical gradients in maximum-likelihood optimization, the analytic gradient yields per-iteration speedups of 89.9× on MPXV and 4.3× on BRCA1, with total optimization speedups of 271.7× and 1.4× respectively. In Bayesian inference, HMC using the exact gradient improves ESS/min for the slowest-mixing dimension from 0.0095 to 3.41 on MPXV and from 1.51 to 1.86 on BRCA1, corresponding to 359.6× and 1.23× speedups over a univariate random-walk sampler (Ji et al., 11 Jul 2025).

At the level of branch-length optimization rather than substitution-process estimation, lcfit approximates real likelihood curves well across binary models, JC, HKY85, JTT92, LG08, YN98, and a nonhomogeneous DNA model with mixed substitution models across branches, under both uniform rates across sites and discretized gamma rate variation. In Bayesian proposals it gives slightly better acceptance rates overall than the Gamma and Weibull proposals studied by Aberer et al., and iterative lcfit improves ML branch-length estimation substantially, with low relative error for branches $e^{Q_1 t_1}e^{Q_2 t_2}$ 1, small absolute error for short branches, and faster performance than Brent’s method when approximate branch lengths are sufficient (Claywell et al., 2017).

Outside phylogenetic CTMCs, the idea of branch- or family-specific substitution behavior appears in several related but distinct forms.

In protein family analysis, environment-specific substitution tables model the expected substitution behavior of amino acids in a particular structural environment, including secondary structure, solvent accessibility, and hydrogen-bonding context. These tables were used to distinguish structural from functional constraint and to identify specificity-determining residues by comparing whole-family MSSAs with EC-specific partitions. Functional restraint is defined as the city-block distance between observed and expected substitution patterns, and the final dataset contains 97 examples from 68 families (0710.2808).

In B-cell receptor repertoire analysis, SPURF predicts clonal-family-specific substitution profiles from a single input sequence using a penalized tensor regression framework that combines public repertoire profiles, simulated neutral profiles, and germline information. These profiles are site-specific amino acid frequency vectors rather than CTMC rate matrices, but they serve an analogous role as lineage-specific substitution summaries. On held-out data and the external Briggs dataset, the reported $e^{Q_1 t_1}e^{Q_2 t_2}$ 2 errors are 0.0492 and 0.0511 versus baseline values of 0.114 and 0.129, and Jaccard similarities at $e^{Q_1 t_1}e^{Q_2 t_2}$ 3 are 0.9289 and 0.9227 versus baseline values of 0.9156 and 0.9053 (1802.06406).

A different terminological issue arises in symbolic dynamics and tiling theory, where “branch” and “substitution” have unrelated meanings. In substitution tiling spaces, the branch locus is a subspace of dimension at most $e^{Q_1 t_1}e^{Q_2 t_2}$ 4 summarizing asymptotic behavior “in at least a half-space,” and in the 2D Pisot case it has a description as an inverse limit of an expanding Markov map on a 1-dimensional simplicial complex (Barge et al., 2011). In one-sided substitution subshifts, a branch point is a point that is not $e^{Q_1 t_1}e^{Q_2 t_2}$ 5-invertible, and branch points are analyzed through suffix cycles and ordered Bratteli diagrams (Yassawi, 2011). These are not phylogenetic substitution models, but the terminological overlap can be a source of confusion.

In the phylogenetic sense, the mature literature now spans three complementary themes: biological realism through lineage-dependent processes, algebraic consistency through closed model families, and computational tractability through surrogate likelihoods, GPU-accelerated pruning, and exact linear-time gradients. A plausible implication is that future work will continue to combine these themes rather than treat model expressiveness and scalability as separate problems.