## Tuesday, December 7, 2010

### Models of sequence evolution 3

I wrote a couple of long posts about models of sequence evolution (here and here), in which (and it may not be obvious) we solved the problem of how to obtain:

`P = eQt`

for any t, by getting eigenvalues and eigenectors, using both R and Python. That's outstanding.

If you look carefully, you might notice something. There are just two terms in Q for the JC69 model:

`-0.01 (on diagonal) and 0.003333 (off diagonal)`

P has only two terms as well and they turn out to be (for t = 0.1):

`0.999 and 0.000333`

That is, the off-diagonal terms are apparently just multiplied by t, and naturally the diagonal is 1 - 3 times the other term. There was really no need for all the exponentiation of matrices, etc.

How does this happen? The formulas for the two terms in JC69 are:

 `PNN(t) = 1/4 + 3/4 exp(-4λt)PNM(t) = 1/4 - 1/4 exp(-4λt)`

But going back to the infinite series expansion for e (or just remembering the approximation for small x):

 `ex = 1 + x1/1! + x2/2! + ..ex = 1 + x`

So the exponential is just `1 - 4λt` and we obtain:

 `PNN(t) = 1/4 + 3/4 (1 - 4λt) = 1 - 3λtPNM(t) = 1/4 - 1/4 (1 - 4λt) = λt`

Let's take another look at the formulas for these terms in the TN93 model. Look at only the first row of P. We first consider the two exponential terms:

 `e2 = exp(-βt)e4 = exp[-(πYα1 + πRβ)t]`

Approximating for small t:

 `e2 = 1 - βte4 = 1 - (πYα1 + πRβ)t`

Look what happens to `P13` and `P14`:

 `P13 = &piA(1 - e2) = &piAβt`

This is t times `Q13`. Exactly the same thing happens for `P14`. `P12` is a bit more complicated, but it's just algebra:

 `P12 = πC + (πCπR/πY)e2 - (πC/πY)e4`

The second term is:

 `πCπR/πY (1 - βt)πCπR/πY - (πCπR/πY)βt`

The third term is:

 `πC/πY (-1 + πYα1t + πRβt)`

We can see that (considering both the second and the third term) there are two copies of: `(πCπR/πY)βt` with opposite sign, leaving:

 `P12 = πC + πCπR/πY - πC/πY + πCα1t`

The last term is t times `Q12`. So I'm betting that the rest will go away. We have:

 `πC + πCπR/πY - πC/πY`

but

 `πY + πR = 1πR - 1 = -πY`

so

 `πC + πCπR/πY - πC/πY = πC + πC/πY(πR - 1) = πC - πC = 0`

 `P12 = πCα1t = t Q12`