Maximum-Likelihood Phylogenetics

Maximum-likelihood phylogenetics is a parametric statistical methodology that estimates the evolutionary relationships among homologous biological sequences by identifying the tree topology τ, branch-length vector t, and substitution model parameters θ that maximize the probability of generating the observed multiple sequence alignment D. Originally formulated for discrete gene-frequency characters by Ronald Fisher and Luigi Luca Cavalli-Sforza and extended to nucleotide sequences by Joseph Felsenstein in 1981, the likelihood function is defined as L(τ, t, θ | D) = P(D | τ, t, θ). Unlike distance-matrix clustering (such as Neighbor-Joining) or maximum parsimony, maximum likelihood incorporates an explicit continuous-time Markov chain model of character state transitions along every branch of the phylogeny, accounting for unobserved multiple hits, back-mutations, compositional bias, and site-to-site rate heterogeneity.

Under the standard assumption that alignment columns D₁, D₂, ..., D_L evolve independently according to the same tree and model parameters, the joint likelihood across L sites factors into the product of individual site likelihoods, and numerical optimization is performed on the log-likelihood surface ℓ(τ, t, θ) = ∑_h=1^L ln P(D_h | τ, t, θ). Because evaluating a single site likelihood requires marginalizing over all A^N−1 unobserved ancestral character assignments at the internal nodes of an N-taxon tree (where A = 4 for DNA or RNA, A = 20 for protein sequences, and A = 61 for sense triplets in a codon substitution model), modern phylogenetic engines rely on Felsenstein's tree-pruning algorithm to evaluate each site in O(N · A²) operations while heuristically searching the combinatorial space of (2N − 5)!! tree topologies.

Markov Generator Matrices and Felsenstein's Pulley Principle

Along any branch of length t (measured in expected substitutions per site), character evolution is governed by a transition probability matrix P(t) = exp(Qt), where Q is an A × A instantaneous rate matrix whose off-diagonal entries q_ij ≥ 0 specify the instantaneous rate of substitution from state i to state j, and whose diagonal entries are set to q_ii = −∑_j≠i q_ij so that rows sum to zero. Most nucleotide and amino-acid models—including Jukes-Cantor (JC69), Hasegawa-Kishino-Yano (HKY85), General Time-Reversible (GTR), and empirical protein matrices such as LG and WAG—enforce time-reversibility, parameterized as q_ij = s_ij π_j with symmetric exchangeabilities s_ij = s_ji and stationary state frequencies π = (π₁, ..., π_A). Time-reversibility satisfies the detailed balance identity π_i P_ij(t) = π_j P_ji(t) for all branch lengths t ≥ 0.

Through Felsenstein's tree-pruning algorithm, a post-order traversal computes the partial likelihood vector L_v^(h) at every internal node v from its descendant subtrees. Detailed balance gives rise to Felsenstein's Pulley Principle: under a stationary, time-reversible substitution model, the likelihood P(D_h | τ, t, θ) is invariant to the position of the root node along any branch of the tree, provided the sum of the two sub-branch lengths connected to the root remains constant. Consequently, maximum-likelihood inference evaluates unrooted binary trees containing 2N − 3 branches for N taxa, placing a virtual root on whichever internal edge is undergoing branch-length optimization via Newton-Raphson or Brent's method without recomputing conditional likelihood vectors outside the focal branch.

Discrete Gamma (+G4) and FreeRate (+R) Site Rate Heterogeneity

Biochemical and structural constraints cause evolutionary rates to vary by several orders of magnitude across positions in a gene or protein: catalytic core residues and ribosomal stem helices evolve slowly under strong purifying selection, whereas solvent-exposed surface loops and third-codon positions accumulate substitutions rapidly. Ignoring this variation leads to systematic underestimation of evolutionary distances on long branches and triggers long-branch attraction (LBA), a statistical inconsistency in which rapidly evolving lineages are erroneously grouped together regardless of true ancestry.

To model rate variation across sites, Ziheng Yang introduced the discrete Gamma (+Γ₄ or +G4) model, which assumes that the relative substitution rate r at each alignment site is drawn from a continuous Gamma(α, α) probability distribution with unit mean and variance 1/α. When the shape parameter α < 1 (typically α ∈ [0.15, 0.80] for empirical alignments), the distribution is strongly right-skewed, indicating that most sites evolve slowly while a minority experience rapid turnover. For computational tractability, the continuous Gamma density is approximated by K = 4 equiprobable discrete categories with representative rates r₁, r₂, r₃, r₄ (each weighted by 1/K = 0.25), scaling the transition matrix to P(r_k t) and averaging the site likelihood across the K categories:

P(D_h | τ, t, θ, α) = ∑_k=1^K (1 / K) P(D_h | τ, r_k t, θ)

When empirical rate distributions are multimodal or deviate from Gamma parametric constraints, the FreeRate (+R_K) model (Soubrier et al., 2012) generalizes +G4 by estimating both the category rates r_k and their mixture weights w_k directly from the alignment subject to ∑_k=1^K w_k = 1 and ∑_k=1^K w_k r_k = 1, routinely combined with an explicit proportion of invariant sites (+I).

Topological Rearrangement Operators: NNI, SPR, and TBR

Whereas branch lengths and substitution parameters are continuous variables optimized on a smooth surface, the number of unrooted binary tree topologies for N ≥ 3 taxa is given by the double factorial (2N − 5)!! = 1 × 3 × 5 × ... × (2N − 5), which grows faster than an exponential function exp(N ln N). For N = 10 taxa there are 2,027,025 unrooted trees; for N = 50 taxa there are ~2.75 × 10⁷⁴ unrooted trees, rendering exhaustive enumeration intractable beyond ~12–14 taxa (and establishing maximum-likelihood tree search as NP-hard).

Maximum-likelihood programs therefore construct one or more starting trees via rapid distance methods (such as BioNJ) or randomized stepwise addition parsimony, and then iteratively Climb the log-likelihood landscape using three hierarchical topological rearrangement operators:

• Nearest-Neighbor Interchange (NNI): Swaps two of the four subtrees incident to a single internal branch, generating 2 candidate topologies per internal branch and 2(N − 3) neighbors across the entire tree. NNI is fast—requiring only local partial-likelihood updates—but can become trapped on local likelihood peaks separated by multiple branch swaps.
• Subtree Pruning and Regrafting (SPR): Detaches a subtree by severing an edge, suppresses the resulting degree-2 node, and re-inserts the pruned subtree into any of the remaining O(N) edges of the residual tree, yielding O(N²) candidate topologies per sweep. Bounded-radius lazy SPR moves form the primary search engine in RAxML-NG and PhyML.
• Tree Bisection and Reconnection (TBR): Cuts an internal branch to split the tree into two separate subtrees and reconnects them by joining any branch of the first subtree to any branch of the second, exploring O(N³) candidate topologies and escaping deep local optima during broad topological exploration.

Ultrafast Bootstrap (UFBoot), SH-aLRT, and IQ-TREE 2 / RAxML-NG Benchmarks

To quantify statistical confidence in inferred clades without running hundreds of independent full-likelihood tree searches (Felsenstein's standard non-parametric bootstrap), modern phylogenomic workflows combine two rapid branch-support metrics:

1. Ultrafast Bootstrap Approximation (UFBoot / UFBoot2): Developed by Minh et al. (2013) and Hoang et al. (2018) in IQ-TREE 2, UFBoot uses the Resampling Estimated Log-Likelihoods (RELL) method to evaluate thousands of bootstrap replicate weight vectors simultaneously against candidate trees encountered during the stochastic NNI/parsimony search. Unlike standard bootstrap proportions (which are conservative), UFBoot values are unbiased probabilities when ≥ 95%, meaning a clade with UFBoot ≥ 95% has a ~95% probability of being correct.
2. Shimodaira-Hasegawa-like Approximate Likelihood Ratio Test (SH-aLRT): Introduced by Guindon et al. (2010), SH-aLRT compares the log-likelihood of the optimal four-subtree arrangement around an internal branch against the second-best of the three possible NNI configurations, using RELL bootstrap resampling of site log-likelihood differences to account for model misspecification. Empirical benchmarks show that clades supported by both SH-aLRT ≥ 80% and UFBoot ≥ 95% achieve robust false-positive control across genome-scale datasets.

In comparative benchmarks across multi-gene supermatrices, IQ-TREE 2 (combining ModelFinder BIC/AICc model selection, stochastic NNI perturbation, and SIMD-vectorized pruning kernels) and RAxML-NG (employing fine-grained MPI/pthreads parallelization, site-pattern compression, and adaptive SPR radius tuning) converge to virtually identical log-likelihood maxima within 0.01 log-likelihood units while scaling to alignments of > 10,000 taxa and > 10⁶ nucleotide or amino-acid sites.

How does maximum likelihood overcome long-branch attraction compared to maximum parsimony?

Maximum parsimony counts only the minimum number of character state changes needed to explain an alignment column and assumes every branch has an equal probability of change, causing two long, rapidly evolving branches to share convergent identical states by chance. Maximum likelihood explicitly models branch length t in expected substitutions per site through the matrix exponential P(t) = exp(Qt) and site rate variation (+G4/+R), discounting convergent changes on long branches and remaining statistically consistent when the substitution model matches the evolutionary process.

How should researchers choose among nucleotide or protein substitution models?

Model selection is performed using information-theoretic criteria that penalize excess free parameters k to prevent overfitting—specifically the Bayesian Information Criterion (BIC = −2 ℓ_max + k ln L, where L is alignment length) or the sample-size-corrected Akaike Information Criterion (AICc = −2 ℓ_max + 2k + 2k(k + 1)/(L − k − 1)). Automated tools such as ModelFinder in IQ-TREE 2 evaluate candidate matrices, base frequencies (+F), invariant sites (+I), and rate heterogeneity (+G4 or +R) on a fixed baseline tree before initiating the full topological search.

What is phylogenetic terraces in partitioned supermatrix analyses?

In phylogenomic supermatrices where multiple genes are concatenated across N taxa but several species lack sequence data for subsets of loci (missing gene partitions), distinct unrooted tree topologies can display identical induced subtrees for every individual gene locus. When branch lengths are unlinked across partitions, these distinct topologies yield the exact same composite log-likelihood, forming a flat "phylogenetic terrace" of equally optimal trees that can only be resolved by adding overlapping locus coverage or linking branch-length proportionality across partitions.

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. Guindon, S., Dufayard, J.-F., Lefort, V., Anisimova, M., Hordijk, W., & Gascuel, O. (2010). "New algorithms and methods to estimate maximum-likelihood phylogenies: Assessing the performance of PhyML 3.0." Systematic Biology, 59(3), 307–321. https://doi.org/10.1093/sysbio/syq010
3. Kozlov, A. M., Darriba, D., Flouri, T., Morel, B., & Stamatakis, A. (2019). "RAxML-NG: A fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference." Bioinformatics, 35(21), 4453–4455. https://doi.org/10.1093/bioinformatics/btz305
4. Minh, B. Q., Schmidt, H. A., Chernomor, O., Schrempf, D., Woodhams, M. D., von Haeseler, A., & Lanfear, R. (2020). "IQ-TREE 2: New models and efficient methods for phylogenetic inference in the genomic era." Molecular Biology and Evolution, 37(5), 1530–1534. https://doi.org/10.1093/molbev/msaa015