Skip to content

Review: Evolutionary Trees from DNA Sequences: A Maximum Likelihood Approach

Citation

  • Felsenstein, J. (1981). Evolutionary trees from DNA sequences: a maximum likelihood approach. Journal of Molecular Evolution, 17(6), 368–376.
  • DOI
  • PubMed

Abstract

The application of maximum likelihood techniques to the estimation of evolutionary trees from nucleic acid sequence data is discussed. A computationally feasible method for finding such maximum likelihood estimates is developed, and a computer program is available. This method has advantages over the traditional parsimony algorithms, which can give misleading results if rates of evolution differ in different lineages.


The central problem in phylogenetics is to infer which tree, among the astronomically many possible branching patterns, best explains a given alignment of DNA or protein sequences. Before Felsenstein (1981), maximum likelihood was considered computationally out of reach for this problem. The key insight of the paper is that the tree structure can be exploited via dynamic programming: the likelihood factors over the branches, and partial results computed from the leaves toward the root never need to be recomputed.

The algorithm processes the tree in postorder (leaves first, root last). For each internal node, the conditional likelihood \( \mathcal{L}_k(x) \) gives the probability of all observed data in the subtree rooted at that node, given that the ancestral character state at the node is \( x \). This is computed from the children via the pruning recursion (Eq. 4, p. 370):

\[\mathcal{L}_k(x) = \left(\sum_{s_i} P_{x s_i}(v_i)\,\mathcal{L}_i(s_i)\right) \left(\sum_{s_j} P_{x s_j}(v_j)\,\mathcal{L}_j(s_j)\right)\]

where \( v_i \) and \( v_j \) are the branch lengths to the two children and \( P_{xy}(t) \) is the transition probability. For tip nodes, the conditional likelihood is 1 for the observed character and 0 for all others. At the root (Eq. 5, p. 370), the full likelihood is recovered by weighting over stationary frequencies \( \pi_{s_0} \):

\[L = \sum_{s_0} \pi_{s_0}\,\mathcal{L}_{s_0}^{(0)}\]

The paper uses a specific substitution model in which each nucleotide is replaced by any other at rate \( u \), with the target nucleotide drawn from a frequency distribution \( \pi \). This yields a closed-form transition probability (Eq. 7, p. 371):

\[P_{ij}(t) = e^{-ut}\,\delta_{ij} + (1 - e^{-ut})\,\pi_j\]

This is the forerunner of all GTR models: the Jukes-Cantor model with heterogeneous stationary frequencies. Tavaré (1986) later generalized it to the full six-parameter GTR matrix.

A key theoretical result is the Pulley Principle (Appendix, p. 376): when the model satisfies detailed balance \( \pi_i P_{ij}(t) = \pi_j P_{ji}(t) \), the likelihood is identical regardless of where the root is placed on the unrooted tree. This makes the root an irrelevant nuisance parameter, and Hifuku follows standard practice by working with unrooted trees throughout. All of Hifuku's substitution models satisfy this condition (see Section 4.1 of the theory document).

Hifuku's pruning recursion in log_likelihood and log_likelihood_gpu is a direct implementation of Eq. (4) and Eq. (5), generalized to arbitrary GTR matrices and extended with per-node float32 rescaling (see Section 3.2) to prevent underflow on long branches.

In-document navigation

Equation Location
Pruning recursion Eq. (4), p. 370
Root likelihood Eq. (5), p. 370
Infinitesimal substitution Eq. (6), p. 371
Transition probability Eq. (7), p. 371
Detailed balance Eq. (8), p. 371
Pulley Principle proof Appendix, p. 376, Eqs. (A1)–(A2)