Codon Substitution Models

Codon substitution models are continuous-time Markov chain models of molecular evolution that treat the sense triplets of the genetic code—rather than individual nucleotides or translated amino acids—as the fundamental discrete states evolving along the branches of a phylogenetic tree. Introduced independently in 1994 by Nick Goldman and Ziheng Yang and by Spencer V. Muse and Brandon S. Gaut, codon models bridge the gap between mutational processes acting at the DNA level (such as transition/transversion bias and uneven nucleotide composition) and natural selection acting on protein structure and function at the amino acid level.

In the universal genetic code, 3 of the 64 possible nucleotide triplets (TAA, TAG, and TGA) are nonsense stop codons whose premature appearance within an open reading frame typically truncates the polypeptide and is eliminated by purifying selection. The remaining 61 sense codons define the state space of a 61 × 61 instantaneous substitution rate matrix Q = {q_ij} (or 60 × 60 in vertebrate mitochondrial DNA, where AGA and AGG are also stop codons). By parameterizing the instantaneous rates of synonymous (silent) and nonsynonymous (amino-acid-altering) substitutions within a unified maximum-likelihood phylogenetics framework, codon substitution models provide the standard statistical machinery for estimating the selective pressure parameter ω = d_N/d_S across sites and lineages.

Structure of the 61 × 61 Goldman-Yang (GY94) Matrix

In the Goldman-Yang (GY94) formulation (implemented in the PAML codeml package), the instantaneous rate q_ij of substitution from sense codon i = (i₁, i₂, i₃) to sense codon j = (j₁, j₂, j₃) (for i ≠ j) assumes that mutations occur one nucleotide at a time. Any pair of codons differing at two or three nucleotide positions has an instantaneous transition rate of zero (although multiple substitutions can occur over a finite branch length t via the matrix exponential distribution P(t) = exp(Qt)). For codons differing at a single nucleotide position, q_ij is factored into mutational and selective terms:

• q_ij = 0, if codons i and j differ at more than one of the three codon positions;
• q_ij = μ π_j, if i → j is a synonymous transversion;
• q_ij = μ κ π_j, if i → j is a synonymous transition (A ↔ G or C ↔ T);
• q_ij = μ ω π_j, if i → j is a nonsynonymous transversion;
• q_ij = μ κ ω π_j, if i → j is a nonsynonymous transition.

Here, κ is the transition/transversion rate ratio, π_j is the equilibrium frequency of target codon j (satisfying ∑_j=1⁶¹ π_j = 1), and ω = d_N/d_S is the acceptance ratio of nonsynonymous to synonymous mutations:

1. ω < 1: Indicates purifying selection (negative selection), where deleterious amino acid replacements are eliminated and only a fraction ω of nonsynonymous mutations reach fixation.
2. ω = 1: Indicates neutral evolution (as in pseudogenes), where nonsynonymous mutations fix at the same rate as synonymous mutations.
3. ω > 1: Indicates diversifying positive selection (adaptive evolution), where advantageous amino acid replacements fix at a higher rate than neutral synonymous drift.

The diagonal elements of Q are set to q_ii = −∑_j≠i q_ij so that every row sums to zero, and the scaling factor μ is chosen so that the mean substitution rate across all sense codons equals one expected nucleotide substitution per codon per unit branch length: −∑_i=1⁶¹ π_i q_ii = 1.

Muse-Gaut (MG94) Target Nucleotide Parameterization and F1×4, F3×4, F61

Whereas the Goldman-Yang (GY94) model multiplies every instantaneous rate q_ij by the equilibrium frequency π_j of the entire target codon j, the Muse-Gaut (MG94) model (widely used in HyPhy) multiplies q_ij only by the equilibrium frequency π_nj of the single target nucleotide n_j ∈ {A, C, G, T} at the mutated position k ∈ {1, 2, 3}. Because GY94 uses the full codon frequency π_j, changes in nucleotides at non-mutating positions of the codon can subtly influence q_ij; MG94 isolates single-nucleotide mutational biases from codon-level context.

In both frameworks, the 61 equilibrium codon frequencies π_j are estimated using one of three nested empirical models:

• F1×4 (3 free parameters): Assumes equal nucleotide frequencies (π_A, π_C, π_G, π_T) across all three codon positions, rescaled after removing the 3 stop codons.
• F3×4 (9 free parameters): Estimates distinct nucleotide frequencies (π_A^(k), π_C^(k), π_G^(k), π_T^(k)) at the 1st, 2nd, and 3rd codon positions to capture strongly biased third-position GC content (GC₃).
• F61 (60 free parameters): Estimates all 61 empirical sense codon frequencies directly from the observed sequence alignment, fully accounting for translational codon usage bias.

Spectral Decomposition of P(t) = exp(Qt)

To compute the probability P_ij(t) that codon i evolves into codon j over a phylogenetic branch of length t (measured in expected nucleotide substitutions per codon), the instantaneous rate matrix Q is exponentiated via the matrix exponential distribution:

P(t) = exp(Qt) = ∑_m=0^∞ (Qt)^m / m!

Because both GY94 and MG94 satisfy the detailed balance (time-reversibility) condition π_i q_ij = π_j q_ji for all i and j, the asymmetric matrix Q can be transformed into a real symmetric matrix S = Π^1/2 Q Π^−1/2, where Π = diag(π₁, ..., π₆₁). Diagonalizing S = U Λ U^T yields real eigenvalues Λ = diag(λ₁, ..., λ₆₁) (with λ₁ = 0 and λ_k < 0 for k ≥ 2) and orthonormal eigenvectors U, allowing P(t) to be evaluated rapidly for any branch length t:

P(t) = Π^−1/2 U exp(Λt) U^T Π^1/2

Once P(t) is computed for every branch and ω category, the log-likelihood of the codon alignment is evaluated over the tree topology via Felsenstein's tree-pruning algorithm.

Heterogeneous Selection Across Amino Acid Sites (M0, M1a, M2a, M7, M8)

In almost all functional proteins, the vast majority of structural core residues are under strong purifying selection (ω ≪ 1), while adaptive positive selection (ω > 1) acts on only a small subset of surface or active-site residues. A single average ratio across the whole gene (model M0, one-ratio) averages over all sites and rarely exceeds 1.0 even during intense adaptive evolution. To detect site-specific adaptation, Ziheng Yang and Rasmus Nielsen (1998, 2000) developed random-sites codon models, which treat the selective regime ω_h at each codon column h as a random variable drawn from a discrete mixture distribution across K site classes with proportions p₀, p₁, ..., p_K−1 (where ∑ p_c = 1):

To test statistically for the presence of positively selected sites, a Likelihood Ratio Test (LRT) compares twice the log-likelihood difference, 2Δℓ = 2(ℓ₁ − ℓ₀), between a nested null model that forbids ω > 1 (such as M1a or M7) and an alternative model that adds an extra site class with ω > 1 (such as M2a or M8) against a χ² distribution with 2 degrees of freedom (critical values 5.99 at α = 0.05 and 9.21 at α = 0.01) or via the Akaike information criterion. When the global LRT rejects the null hypothesis, individual positively selected codon sites are pinpointed by computing posterior probabilities via Bayes Empirical Bayes site selection.

Why do codon substitution models have 61 states instead of 64?

Three of the 64 triplets in the standard genetic code (TAA, TAG, and TGA) are nonsense stop codons that terminate translation. Because a mutation from a sense codon to a premature stop codon within a protein-coding gene almost always destroys protein function and is rapidly purged by purifying selection, standard codon models assign zero probability to internal stop codons and restrict the continuous-time Markov chain to the 61 sense codons.

Why were the original M1 (Neutral) and M2 (Selection) models replaced by M1a and M2a?

In the original Yang et al. (2000) M1 and M2 models, the purifying-selection class was fixed strictly at ω₀ = 0, forcing any moderately conserved site with 0 < ω < 1 (such as ω = 0.25) to be absorbed into the neutral class (ω₁ = 1) or the positive-selection class (ω₂). In 2005, Yang, Wong, and Nielsen replaced ω₀ = 0 with a freely estimated parameter 0 < ω₀ < 1 in M1a (NearlyNeutral) and M2a (PositiveSelection), eliminating false positives caused by weakly constrained residues.

What is the difference between site models, branch models, and branch-site models?

Site models (such as M1a/M2a and M7/M8) allow ω to vary among codon columns while holding ω constant across all branches of the phylogenetic tree. Branch models hold ω uniform across sites but allow a designated foreground lineage to have a different ω_fg from background lineages ω_bg. Branch-site models (such as Model A) allow ω to vary simultaneously across both codon sites and specific foreground branches, detecting episodic bursts of positive selection along a single evolutionary lineage.

1. Goldman, N., & Yang, Z. (1994). "A codon-based model of nucleotide substitution for protein-coding DNA sequences." Molecular Biology and Evolution, 11(5), 725–736. https://doi.org/10.1093/oxfordjournals.molbev.a040153
2. Muse, S. V., & Gaut, B. S. (1994). "A likelihood approach for comparing synonymous and nonsynonymous nucleotide substitution rates, with application to the chloroplast genome." Molecular Biology and Evolution, 11(5), 715–724. https://doi.org/10.1093/oxfordjournals.molbev.a040152
3. Yang, Z., Nielsen, R., Goldman, N., & Pedersen, A.-M. K. (2000). "Codon-substitution models for heterogeneous selection pressure at amino acid sites." Genetics, 155(1), 431–449. https://doi.org/10.1093/genetics/155.1.431
4. Yang, Z. (2007). "PAML 4: Phylogenetic Analysis by Maximum Likelihood." Molecular Biology and Evolution, 24(8), 1586–1591. https://doi.org/10.1093/molbev/msm088