Richer models of sequence evolution

Rates, continuous time, and transition probabilities

Last time

The Jukes-Cantor (JC) model: Every change is equally likely

But in the yeast SSA3 gene, observed transitions outnumbered transversions about 3 to 1 (each transition > 6 times as fast as each transversion)

And when we simulated transition-biased evolution, the JC distance underestimated the true distance

A model that ignores how sequences really evolve gives biased answers

A second rate: K80

Kimura’s two-parameter model1

Two rates for the two kinds of change:

A C G T
A — \(\beta\) \(\alpha\) \(\beta\)
C \(\beta\) — \(\beta\) \(\alpha\)
G \(\alpha\) \(\beta\) — \(\beta\)
T \(\beta\) \(\alpha\) \(\beta\) —

\(\kappa = \alpha / \beta\): How many times faster each transition is than each transversion

\(\kappa = 1\) is the JC model

The K80 distance

\(S\) = proportion of sites that differ by a transition

\(V\) = proportion of sites that differ by a transversion

\[ d = -\frac{1}{2}\ln\left(1 - 2S - V\right) - \frac{1}{4}\ln\left(1 - 2V\right) \]

Yeast SSA3 JC K80 Estimated \(\kappa\)
Whole gene 0.182 0.187 6.2
3rd codon positions 0.570 0.655 9.7

The correction matters most where sequences have diverged the most

Continuous time

Do we need time steps?

Last time, we chopped time into tiny steps, each with a small chance of substitution

As the steps get smaller, we approach a continuous-time Markov chain

Instead of a probability per step, we describe the model by rates: How many changes we expect per unit of time

Question: If a site changes at rate \(\lambda\), how long until it changes?

Waiting times are exponential

Figure 1

Average wait \(= 1 / \text{rate}\)

The exponential is memoryless

  • A site that has gone a long time without changing is not “due” for a change
  • The time until the next change does not depend on how long it has already waited
  • This is the Markov property in continuous time

So to simulate, we just repeat:

  1. Calculate the total rate of leaving the current base
  2. Draw an exponential waiting time using the total rate
  3. Choose the new base in proportion to its rate

Unequal base frequencies

Not all genomes are 25% A, C, G, & T

  • Some genomes are GC-rich, others AT-rich
  • The HKY model1 adds base frequencies: \(\pi_A\), \(\pi_C\), \(\pi_G\), \(\pi_T\)
  • The rate of change to a base is proportional to its frequency

Stationary distribution

A site reaches equilibrium at \(\pi\), not \(\frac{1}{4}\) each

Saturation

Two random draws from \(\pi\) match with probability \(\sum_i \pi_i^2\)

With 40% C and G: sequences saturate at \(p \approx 0.66\), not 0.75

From rates to probabilities

The rate matrix, \(Q\)

The table of rates, with minus the total rate of leaving on the diagonal:

\[ Q_{\text{K80}, \kappa = 8} = \begin{bmatrix} -1 & 0.1 & 0.8 & 0.1 \\ 0.1 & -1 & 0.1 & 0.8 \\ 0.8 & 0.1 & -1 & 0.1 \\ 0.1 & 0.8 & 0.1 & -1 \end{bmatrix} \]

  • Rows and columns: A, C, G, T
  • Every row sums to zero
  • Rates scaled so time is in expected substitutions per site

Transition probabilities

The probability that a site that starts as base \(i\) is base \(j\) after time \(t\):

\[ P(t) = e^{Qt} \]

  • The matrix exponential: A version of \(e^x\) for matrices
  • We won’t calculate it by hand
    • We’ll use scipy.linalg.expm(Q * t)
  • For JC, \(e^{Qt}\) gives exactly last week’s formula
  • For every other model, \(e^{Qt}\) gives what we can’t easily write down

What \(P(t)\) looks like

K80 (\(\kappa = 8\)), \(t = 0.5\):

From \ To A C G T
A 0.658 0.045 0.251 0.045
C 0.045 0.658 0.045 0.251
G 0.251 0.045 0.658 0.045
T 0.045 0.251 0.045 0.658
  • Every row sums to 1; why?
  • Question: What will it look like when \(t = 0\)?
  • Question: What will it look like when \(t\) is very large?

The family of models, revisited

Model Rates Base frequencies Free parameters
JC691 All equal Equal 0
K802 Transitions vs. transversions (\(\kappa\)) Equal 1
F813 All equal Unequal (\(\pi\)) 3
HKY854 \(\kappa\) Unequal (\(\pi\)) 4
GTR5 Every pair its own rate Unequal (\(\pi\)) 8

Each model is a rate matrix \(Q\), which enables us to get \(P(t) = e^{Qt}\)

Which model is better?

Likelihood under any model

For two aligned sequences, product over sites (base \(x\) in one sequence, base \(y\) in the other):

\[ L(t) = \prod_{\text{sites}} \pi_x P_{xy}(t) \]

For our computer’s sake, we turn the product into a sum of logs:

\[ \ln L(t) = \sum_{\text{sites}} \Big[\ln(\pi_x) + \ln\big(P_{xy}(t)\big)\Big] \]

Likelihood of yeast SSA3 sequences

3rd codon positions:

Model Best \(t\) Log-likelihood
JC 0.570 −1563.6
K80 (\(\kappa = 9.66\)) 0.655 −1473.8

The \(t\) that maximizes the likelihood under each model is that model’s distance formula!

More parameters always help

A model with an extra parameter always fits at least as well; why?

So we need to ask: Is the improvement worth the extra parameter?

  • Likelihood ratio test: Is \(2 \times\) the difference in log-likelihood bigger than chance alone would produce?
  • AIC1: “Charge” each parameter 1 log-likelihood unit (2 AIC units), and pick the model with the best score

K80 vs. JC for the yeast gene: About 90 log-likelihood units better, for 1 extra parameter

Summary

  • K80 gives transitions their own rate; its distance corrects JC’s bias when transitions are faster/slower than transversions
  • In continuous time, the wait until the next substitution is exponentially distributed and memoryless
  • With unequal base frequencies, sites equilibrate at \(\pi\) and saturate below a p-distance of 0.75
  • Any model is a rate matrix \(Q\); \(P(t) = e^{Qt}\) gives its transition probabilities
  • With \(P(t)\), we can calculate the likelihood under any model, and compare models

In the exercise

  1. count_transitions_transversions and k80_distance: Fix last week’s distance bias
  2. simulate_site: Simulate evolution in continuous time
  3. Watch a GC-rich model “settle” at its base frequencies
  4. get_rate_matrix and get_transition_probabilities: Build \(Q\) and \(P(t)\) with NumPy and SciPy
  5. log_likelihood_pair: Compare JC and K80 on the yeast gene

Watch for

Your simulation and \(P(t)\) are two completely different calculations. Do they agree?

Exercise

  1. Log into molab
  2. Open the exercise in molab
  3. Click “Fork” button

References

Akaike, Hirotugu. 1974. “A New Look at the Statistical Model Identification.” IEEE Transactions on Automatic Control 19 (6): 716–23.
Felsenstein, Joseph. 1981. “Evolutionary Trees from DNA Sequences: A Maximum Likelihood Approach.” Journal of Molecular Evolution 17: 368–76.
Hasegawa, Masami, Hirohisa Kishino, and Taka-aki Yano. 1985. “Dating of the Human-Ape Splitting by a Molecular Clock of Mitochondrial DNA.” Journal of Molecular Evolution 22 (2): 160–74.
Jukes, Thomas. H., and Charles. R. Cantor. 1969. “Evolution of Protein Molecules.” Chap. 24 in Mammalian Protein Metabolism, edited by H. N. Munro, III. Academic Press.
Kimura, Motoo. 1980. “A Simple Method for Estimating Evolutionary Rates of Base Substitutions Through Comparative Studies of Nucleotide Sequences.” Journal of Molecular Evolution 16: 111–20. https://doi.org/10.1007/BF01731581.
Tavaré, Simon. 1986. “Some Probabilistic and Statistical Problems in the Analysis of DNA Sequences.” In Some Mathematical Questions in Biology: DNA Sequence Analysis, edited by Robert M. Miura. American Mathematical Society.

Acknowledgements

I drafted and iterated on these slides in collaboration with Anthropic’s Claude (Opus 5.5 model).