## Sunday, February 28, 2010

### 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 TA -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*λ*tPXY(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 00 λ2 0 00 0 λ3 00 0 0 λ4`

And that means that `Λ Λ` is simply:

 `λ12 0 0 00 λ22 0 00 0 λ32 00 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 00 eλ2 0 00 0 eλ3 00 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).