Felsenstein's Tree-Pruning Algorithm

Felsenstein's tree-pruning algorithm is a foundational dynamic-programming algorithm introduced by Joseph Felsenstein in 1973 and popularized for nucleotide sequences in 1981 that computes the exact likelihood of a molecular sequence alignment on a fixed phylogenetic tree topology and branch-length vector. In maximum-likelihood phylogenetics and Bayesian inference, the characters observed at the tips (leaves) of an evolutionary tree are known from DNA sequencing, whereas the ancestral character states at the N − 1 internal nodes of a rooted binary tree of N taxa are unobserved latent variables.

Calculating the marginal probability P(x_h | τ, t) of observing alignment column x_h requires summing the joint probability over all A^N−1 possible combinations of ancestral character states at the internal nodes, where A is the alphabet size (A = 4 for nucleotides, A = 20 for amino acids, or A = 61 for sense triplets in codon substitution models). Naive brute-force summation over all A^N−1 ancestral state histories grows according to an exponential function exp((N − 1) ln A), reaching 4⁹⁹ ≈ 4.0 × 10⁵⁹ terms for a modest 100-taxon DNA tree. Felsenstein's pruning algorithm exploits the conditional independence of evolutionary lineages under a Markov chain model to move summation signs inward across the tree hierarchy, reducing the computational complexity per site from O(A^N) to linear time O(N · A²).

Conditional Subtree Likelihood Vectors L_v(s)

Consider a rooted binary phylogenetic tree τ with branch lengths t and a continuous-time Markov substitution model defined by equilibrium state frequencies π_s and branch transition probability matrices P(t) = exp(Qt), where P_sx(t) is the probability that state s ∈ {1, ..., A} at a parent node evolves into state x ∈ {1, ..., A} at a child node separated by branch length t. Assume that sites h = 1, ..., L evolve independently and identically, so that the total log-likelihood of the alignment X is the sum over unique site patterns:

ℓ(τ, t) = ln P(X | τ, t) = ∑_h=1^L ln P(x_h | τ, t)

For a single site h and any node v in the tree, Felsenstein defined the conditional subtree likelihood L_v^(h)(s) as the probability of observing all tip character states at or below node v (in the clade rooted at v), given that node v itself is in character state s ∈ {1, ..., A}.

Initialization at Tip Nodes and Post-Order Internal Recursion

Felsenstein's algorithm evaluates the conditional likelihood vector L_v^(h) = [L_v^(h)(1), ..., L_v^(h)(A)]^T from the tips down to the root via a bottom-up post-order tree traversal (visiting both child nodes u and w before visiting their parent node v):

1. Initialization at Leaf (Tip) Nodes: For any external tip node k ∈ {1, ..., N} where the observed character in alignment column h is x_k:

• Unambiguous Character: L_k^(h)(s) = 1 if s = x_k, and L_k^(h)(s) = 0 for all s ≠ x_k. (For example, if tip k has nucleotide A in state order {A, C, G, T}, its tip vector is [1, 0, 0, 0]).
• Missing Data, Gap (-), or IUPAC Ambiguity Code: Because missing data or ambiguity codes represent uncertainty about which state is observed at tip k, L_k^(h)(s) = 1 for every state s compatible with the symbol (e.g., [1, 1, 1, 1] for N or -, and [1, 0, 1, 0] for purine R).

1. Post-Order Recursion at Internal Nodes: For an internal node v with left child u (connected by branch length t_u) and right child w (connected by branch length t_w), evolution along branch t_u and its descendant clade is conditionally independent of evolution along branch t_w and its descendant clade once the ancestral state s at parent node v is specified. Therefore, for each possible parent state s ∈ {1, ..., A}:

L_v^(h)(s) = [ ∑_x=1^A P_sx(t_u) L_u^(h)(x) ] × [ ∑_y=1^A P_sy(t_w) L_w^(h)(y) ]

Once L_v^(h) is computed from L_u^(h) and L_w^(h), the descendant subtrees below u and w have been algebraically "pruned" into the A-element vector at node v.

1. Marginal Termination at the Root Node (r): When the post-order traversal reaches the root node r, the unconditional site likelihood P(x_h | τ, t) is obtained by weighting each possible root state s by its prior equilibrium frequency π_s:

P(x_h | τ, t) = ∑_s=1^A π_s L_r^(h)(s)

Felsenstein's Pulley Principle for Unrooted Trees

When the underlying substitution model Q is time-reversible—satisfying the detailed balance condition π_i q_ij = π_j q_ji, which in turn ensures π_i P_ij(t) = π_j P_ji(t)— combined with the Chapman-Kolmogorov semigroup property P(t₁ + t₂) = P(t₁)P(t₂), the site likelihood P(x_h | τ, t) is completely invariant to the placement of the root node r anywhere along any branch of the tree. Felsenstein termed this invariance the Pulley Principle.

Because of the Pulley Principle, maximum-likelihood phylogenetics searches operate on unrooted trees (which have 2N − 3 branches for N taxa) and place a virtual root on any internal branch being optimized during Nearest-Neighbor Interchange or Subtree Pruning and Regrafting. By caching the partial likelihood vectors pointing into both ends of an edge, the branch length t_e of any single edge e can be optimized by Newton-Raphson iteration in O(L · A²) time without recomputing the conditional likelihoods of the rest of the tree.

Discrete Gamma (+Γ₄) Rate Heterogeneity Across Sites

In real biological sequences, functional constraints cause substitution rates to vary substantially across alignment positions (for example, between first, second, and third codon positions, or between catalytic active sites and solvent-exposed loops). Following Ziheng Yang (1994), rate heterogeneity across sites is modeled by multiplying all branch lengths at site h by a relative rate factor r_k drawn from a continuous Gamma(α, α) distribution of mean 1 and variance 1/α, approximated by K = 4 equal-probability discrete categories (+Γ₄) with rates r₁, r₂, r₃, r₄ (each having prior weight 1/K = 0.25).

To evaluate the likelihood under a +Γ_K model, Felsenstein's pruning algorithm is executed independently for each rate multiplier r_k (replacing P(t) with P(r_k t)), and the site likelihood is averaged across the K rate categories before taking the logarithm:

P(x_h | τ, t, α) = ∑_k=1^K (1 / K) [ ∑_s=1^A π_s L_r^(h,k)(s) ]

Site-Pattern Compression and Floating-Point Underflow Scaling

Modern phylogenetic engines (such as RAxML, IQ-TREE, and PhyML) accelerate Felsenstein's pruning algorithm via two essential engineering techniques:

1. Site-Pattern Compression (Aliasing): Any two alignment columns x_h and x_m that display the exact same pattern of characters across all N tips have identical conditional likelihood vectors at every node of the tree. Compressing L total alignment columns into L_u ≤ L unique site patterns with integer multiplicity weights w_p reduces the total runtime to O(N · L_u · K · A²).
2. Floating-Point Underflow Rescaling: On trees containing hundreds or thousands of taxa, multiplying probabilities P_sx(t) < 1 across hundreds of internal nodes causes the conditional likelihoods L_v(s) to decay exponentially toward zero as exp(−cN), eventually underflowing double-precision IEEE 754 floating-point representation (~2.2 × 10⁻³⁰⁸). Whenever max_s L_v(s) falls below a safety threshold (such as 2⁻²⁵⁶), the vector L_v is divided by its maximum element M_v = max_s L_v(s) (or scaled by a power of 2), and ln(M_v) is accumulated into a per-site log-scaling ledger added back to the final log-likelihood at the root.

Why is Felsenstein's pruning algorithm O(N · A²) per site instead of O(N · A)?

At each of the N − 1 internal nodes of a binary tree, the algorithm must evaluate the conditional likelihood L_v(s) for each of the A possible parent states s ∈ {1, ..., A}. For each parent state s, it computes the inner product ∑_x=1^A P_sx(t) L_u(x) over all A possible child states x on both the left and right branches, requiring 2A² floating-point multiplications per node and O(N · A²) operations across the entire tree.

How does Felsenstein's pruning algorithm handle gaps (-) and missing data (N)?

Under standard time-reversible Markov substitution models that do not explicitly model insertion-deletion events as an extra state, a gap (-) or unknown nucleotide (N) at tip k is treated as missing data (complete uncertainty about which character state is present). Its tip conditional likelihood vector is initialized to all ones, L_k(s) = 1 for all s ∈ {1, ..., A}, which algebraically prunes away a terminal branch Containing only missing data without biasing the remaining tree likelihood.

What is Felsenstein's Pulley Principle and when does it break down?

Felsenstein's Pulley Principle states that the likelihood of an alignment on a tree is unchanged when the root is moved along any branch, so long as the sum of the two sub-branch lengths attached to the root stays constant. It holds whenever the substitution model satisfies time-reversibility (π_i q_ij = π_j q_ji) and stationarity; it breaks down under non-reversible substitution models (such as general non-reversible 12-rate DNA matrices) or non-stationary base-composition models across lineages.

1. Felsenstein, J. (1981). "Evolutionary trees from DNA sequences: A maximum likelihood approach." Journal of Molecular Evolution, 17(6), 368–376. https://doi.org/10.1007/BF01734359
2. Felsenstein, J. (1973). "Maximum likelihood and minimum-steps methods for estimating evolutionary trees from data on discrete characters." Systematic Zoology, 22(3), 240–249. https://doi.org/10.2307/2412304
3. Yang, Z. (1994). "Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: Approximate methods." Journal of Molecular Evolution, 39(3), 306–314. https://doi.org/10.1007/BF00160154
4. Felsenstein, J. (2004). Inferring Phylogenies. Sinauer Associates. https://doi.org/10.1086/424804