Jukes-Cantor (6)

I'm exploring the derivation of equations for the probability that a nucleotide either stays the same over a period of evolutionary time, or changes to be one of the other three. The Jukes-Cantor model is that all of these rates (X to Y) are the same, described by a parameter α that is the instantaneous rate.

One of my sources is Ziheng Yang's Computational Molecular Evolution. According to him (and wikipedia), this model can be represented by a substitution-rate matrix which is usually called Q:

A C G T A -3λ λ λ λ C λ -3λ λ λ G λ λ -3λ λ T λ λ λ -3λ

Yang uses λ in place of the often-used α. For this post we'll follow along with him.

As I posted about the other day, we need to have equations that will give us the probabilities for any time t. These should be recognizable:

PXX(t) = 1/4 + 3/4*e-4*λ*t PXY(t) = 1/4 - 1/4*e-4*λ*t

Yang calls these p0(t) and p1(t).

The corresponding transition-probability matrix is designated P(t) and it is:

P(t): p0(t) p1(t) p1(t) p1(t) p1(t) p0(t) p1(t) p1(t) p1(t) p1(t) p0(t) p1(t) p1(t) p1(t) p1(t) p0(t)

And if we had P for one unit of evolutionary time, then we could do matrix multiplication (see here) to generate P for any time.

Now, Yang (and other sources) also says that:

P(t) = eQ(t)

And I've been trying to wrap my head around that, i.e. matrix exponentiation. Luckily, I've been getting into Gilbert Strang's MIT lectures on linear algebra (here). So I can appreciate what the inverse of a matrix is, such that:

Q Q-1 = I

where I is just the identity matrix of the same size and shape as Q. Now, suppose we can diagonalize Q. That is, we can find an invertible matrix U and a diagonal matrix Λ such that:

Q = U Λ U-1

A diagonal matrix means that the matrix has non-zero values only on the diagonal, e.g.:

λ1 0 0 0 0 λ2 0 0 0 0 λ3 0 0 0 0 λ4

And that means that Λ Λ is simply:

λ12 0 0 0 0 λ22 0 0 0 0 λ32 0 0 0 0 λ42

Look what happens when we do something like:

Q Q = U Λ U-1 U Λ U-1 = U Λ2 U-1

We only have to deal with Λ, and it's straightforward. Now for the real magic. I am reliably informed that this trick works for any function. For example:

eQ(t) = U e&Lambda(t) U-1

where we have:

e&Lambda(t) = eλ1 0 0 0 0 eλ2 0 0 0 0 eλ3 0 0 0 0 eλ4

And, most impressive of all, U and U^-1 are the right and left eigenvectors of Q, and the coefficients of Λ are its eigenvalues! We've looked at that in detail previously(here).