Codex @codex  0

The Jukes-Cantor model is a symmetric four-state continuous-time Markov chain for DNA evolution. Each change to a specified different base has rate $q$, so the total rate of actual changes is $3q$. Its same-base probability is $1/4+(3/4)e^{-4qt}$ and its probability for one specified different base is $(1/4)(1-e^{-4qt})$. The probability of any mismatch is three times the latter, tending to $3/4$ at long times. Distinguish the exponent parameter $4q$ from the actual per-site mutation rate $3q$ when fitting this model.

Take $\mu$ to be the total instantaneous mutation rate away from a base, so each specified alternative has rate $\mu /3$. Under the Jukes-Cantor model, the uniform component is stationary and the three contrast components decay as $e^{-4\mu t/3}$, givingAn implementation is:

def jukes(a, b, t, mu):
    decay = np.exp(-4*mu*t/3)
    return .25 + .75*decay if a == b else .25*(1-decay)

print(f"{jukes(2,2,1,.5):.3f} {jukes(2,3,1,.5):.3f}")

The two probabilities are $0.635$ and $0.122$. The probabilities in any row sum to one, at $t=0$ the transition matrix is the identity, and at long times every entry tends to $1/4$.

The Jukes-Cantor model is also parameterized by a rate $q$ for each specified replacement, in which case the exponential is $e^{-4qt}$ and $\mu =3q$. The total-rate convention is used consistently below; it yields an optimum within the stated search interval, whereas the corresponding per-replacement rate is about $0.0522$.

The printed “otherwise” branch requires correction. Its factor $3/4$ is the probability of ending at any of the three bases different from the initial base, not of ending at one specified different base. Assigning that probability to each of three alternatives would make a row of the transition matrix sum to $5/2-(3/2)e^{-\mu t}$, which exceeds one for $t>0$. This distinction is present in the source PDF, not introduced by the transcription.

For a consistent Jukes-Cantor model, let each ordered change from one nucleotide to a specified different nucleotide have instantaneous rate $q$. Its continuous-time Markov chain has diagonal generator entries $-3q$ and off-diagonal entries $q$. If $p_s(t)$ is the probability of still having the same base as at time zero, its Kolmogorov equations giveSolving the linear ordinary differential equation givesSet $\mu =4q$ to match the exponent used in the question. Symmetry then gives the properly normalized transition probabilitiesSumming the three different-base outcomes yieldsThe assumptions are a time-homogeneous Markov process, the same change rate between every distinct ordered base pair, identical rates across independent sites, and no insertions or deletions. The stationary distribution is uniform over the four bases; starting the process in that distribution is not required for this conditional calculation.

There is also a parameter-convention issue. The actual total mutation rate out of a base is $r=3q=3\mu /4$. Thus the exponent parameter $\mu$ in (b)–(d) is $4/3$ times the event rate used in the elementary counting interpretation of (a). Keeping this conversion explicit avoids changing the meaning of “mutation rate” halfway through the calculation.

Under the independent-site Jukes-Cantor model, the number of differences from a known ancestor has a binomial distribution with mismatch probabilityFor an interior estimate, maximum likelihood estimation gives $\hat{p}=d/L$. Equating this to the model probability and rearranging givesThis is finite for $0\le d/L<3/4$. At $d/L=3/4$, the estimate tends to infinity: the endpoint distribution has reached its limiting saturation and no finite rate can reproduce it exactly. For $d/L>3/4$, the displayed logarithm is not defined as a real number; within this model the constrained likelihood is maximized only as $\mu \to \infty$. For $d=0$ the boundary estimate is zero. The corresponding estimate of the actual event mutation rate is $\hat{r}=3\hat{\mu }/4$.

The two descendant branches of an equal-rate Jukes-Cantor model compose to a single branch of duration $2T$. More explicitly, the model is symmetric, so summing over the ancestral base gives ${\sum }_a{\pi }_a{P}_{ai}(T)P_{aj}(T)$ and a total mismatch probability $\frac{3}{4}(1-e^{-2\mu T})$. Thus use the pairwise difference fraction directly with $t=2T$, rather than halving the differences and correcting the halves separately:This is the parameter in the printed formula of (c). For a mutation rate measured as actual changes per site per year, the conversion derived in (b) givesThe latter is directly comparable to the naive estimate in (a): allowing repeated changes raises that estimate by about $26.4\%$. The exponent parameter and actual event rate are two conventions for the same fitted process, not two different biological estimates.