# Base Substitution Probabilities Derivation

In Felsenstein's 1981 Paper Evolutionary Trees from DNA Sequences: A Maximum Likelihood Approach, he describes (in equations (6) and (7)) a derivation of the equation of the probability of a particular base substitution. I am having trouble making the connection between equations 6 and 7. That is, going from

$$P_{ij}(dt) = (1- u\ dt)\ \delta_{ij} + u\ dt \ \pi_j \quad\quad [6]$$

to

$$P_{ij}(t) = e^{-ut}\ \delta_{ij} + (1- e^{-ut}) \pi_j \quad\quad [7]$$

where $u$ is the rate of base substitution per unit time, $\pi_j$ is the probabiltiy of observing a particular base and $\delta_{ij}$ is the Kronecker delta function (1 if i=j, 0 if i != j).

In the paragraph that follows these equations there is an explanation as to why the exponential term makes sense (if you already know thats what it is), but I'm not understanding how to get from the first equation to the second. Any insights?

In case anyone is ever interested (or perhaps if I eventually forget), it appears that Felsenstein is using a Poisson Process approximation for equation [7], although even with this new information I still do no see the connection between equations [6] and [7].

Page 9 of this thesis had a pretty good explanation of this.

-

(1) There are two cases. For the first equation (6), the case $i \neq j,$

$$(1) \hspace{13mm} P_{ij}(dt) = u~dt ~\pi_j$$ if we have a transition from one base to another (for example G to C). In the case $i = j,$

$$(2)\hspace{10mm} P_{jj}(dt) = (1 - u~dt) + u~dt~\pi_j$$ holds for a same-base transition (G to G) or a non-transition (say, G and no change). Intuitively, we expect that, all things being equal, the probability of a base to same-base situation is slightly higher: it is the sum of probability that the base changes to itself plus the probability the base sits and does nothing in time $dt$ (with the same net result).

(2) Now we look at the equation(s) (7). Again there are two cases.

$$(1a) \hspace{10mm} P_{ij}(t) = (1 - e^{-ut})\pi_j$$

For case 2 (base to same base transition, $i = j$):

$$(2a)\hspace{13mm} P_{jj}(t) = e^{-ut} + (1 - e^{-ut})\pi_j.$$

(3) The author uses $P(dt),$ $P$ at a very small (infinitesimal) quantity $dt.$ He tells us the leap from (1) to (1a) follows "almost immediately" when we see that $e^{-ut}$ is the probability that no base transition occurs. Well, then $(1 - e^{-ut})$ is the probability that a transition occurs.

Now $e^{-ut} \approx (1 - ut)$ (a truncated Taylor series) or $ut\approx (1 - e^{ut}).$ So for short times $dt$ we have $u~dt \approx (1 - e^{-u dt})$ and this explains his first equality $P(dt) = u~dt~\pi_j.$

So we have to recognize when he tells us about $e^{-ut}$ that $u t$ is the [truncated] Taylor series representation of $(1-e^{-ut})$ for $t$ very near 0. That is how he gets from 1 to 1a and (6) to (7).

(4) $(1 - u)$ is the probability that no transition from one base to another occurs in some short time $t_i.$ The probability that this happens over n consecutive time intervals is $(1 - u)^{(nt_i)}.$ Let $nt_i = T.$

Note that $u \ll 1.$ We can write $(1 - u)^T = (1 - \frac{uT}{T})^T$ and $$\lim_{T \to \infty} (1 - \frac{uT}{T})^T= e^{-uT}.$$ This reasoning comes directly from that used to derive the Poisson distribution which is used to study rare events. See, for example, Deonier, Computational Genome Analysis, p.74.

-
It seems that you are giving a description of why equation [7] refers to the same probability in equation [6]. But this is only possible because you are given equation [7] already. How would one arrive at equation [7] given only equation [6]. Is there a mathematical trick short of reversing the process you showed at the end of your post? – fermats_last_nerve May 14 '14 at 15:56
@fermats_last_nerve: The bottom line is that he only includes (6) because he begins with the probability $u$ and the statement of (6) is an intuitive re-statement of the English. He (or Kaplan before him) already modeled the substitution as a rare event, so he is just trying to justify (7) in the simplest way possible. The power series for $e^{-ut}$ is $1 - ut + \frac{(ut)^2}{2!}-\frac{(ut)^3}{3!}+...$ Someone working in this setting is probably expecting this trick. When I see $dt$ I think integration but his hint is generous enough. – daniel May 17 '14 at 14:02