Models of Nucleic Acid Sequence Evolution and the Jukes Cantor Model

Molecular evolution deals with the changes that accumulate in DNA and protein sequences over time. When two sequences are compared, the number of differences between them tells us something about how much time has passed since they diverged from a common ancestor. However, a simple count of differences underestimates the true amount of evolutionary change, because multiple substitutions can occur at the same site without being detected. To overcome this problem, quantitative evolutionary models are used.

These models allow us to convert the observed difference between two sequences into a corrected evolutionary distance, which is essential for constructing phylogenetic trees, carrying out likelihood-based analyses, and designing scoring systems for sequence alignment. The Jukes–Cantor model is the simplest and most fundamental of these models and forms the basis for understanding more complex substitution models.

1. Why Evolutionary Models Are Needed

When two sequences have descended from a common ancestor, the most direct way to compare them is to align them and calculate the fraction of sites at which they differ. This quantity is denoted D, the observed proportion of different sites. For example, if 3 out of 10 aligned sites differ between two sequences, then D = 3/10 = 0.3.

As the time since divergence increases, more substitutions accumulate, and D increases. However, D does not behave in a simple manner with respect to time, for the following reasons:

• Multiple hits problem: When divergence time is short, each substitution occurs at a different site, so every change is visible when the two sequences are compared. As time increases, however, it becomes increasingly likely that the same site undergoes more than one substitution. If two substitutions occur at the same site, only one difference may be visible, or in some cases (for example, two changes that cancel out, such as G→A→G) no difference at all may be visible. Because of this, the number of visible differences is always less than or equal to the number of substitutions that have actually occurred.
• Saturation: When the time since divergence is extremely large, the sequences become effectively randomized with respect to one another. If all four bases occur with equal frequency, then on average three out of every four sites will differ between two completely randomized sequences. Hence D approaches an upper limit of 3/4, no matter how much additional time passes. This means D increases rapidly at first but levels off (saturates) at large times rather than increasing indefinitely.

1.1 Problems with Using D as a Distance Measure

Two major shortcomings of D make it unsuitable as a true evolutionary distance:

1. Non-linearity with time: D does not increase proportionally with time. A distance measure would ideally be twice as large if the divergence time were twice as long, but D fails to satisfy this because of the saturation effect described above.
2. Non-additivity: If the sequence of a common ancestor were known, one could measure the distance from the ancestor to each of the two descendant sequences, D₀₁ and D₀₂. Ideally, the distance between the two descendants should equal the sum D₁₂ = D₀₁ + D₀₂. This property is called additivity. In practice, D is notadditive — the observed D₁₂ is usually smaller than the sum of the two branch distances, again because of unseen multiple substitutions.

1.2 The Need for a Corrected Distance, d

Because of these limitations, molecular evolutionists define a corrected evolutionary distance, d, as the average number of substitutions that have actually occurred per site between two sequences (not merely the number of visible differences). If substitutions occur randomly at a constant rate, then:

• d is directly proportional to time, and
• d is additive, since the number of substitutions between two descendant sequences is, by definition, the sum of the substitutions that occurred along each branch from the common ancestor.

The difficulty is that d cannot be directly observed simply by comparing two sequences, unlike D. To calculate d from the observable quantity D, an explicit evolutionary model is required. This is precisely the purpose of models such as the Jukes–Cantor model.

2. The Jukes–Cantor (JC) Model

The Jukes–Cantor model, proposed by Jukes and Cantor in 1969, is the simplest possible mathematical model of nucleotide substitution.

2.1 Assumptions of the JC Model

• All four bases (A, C, G, T) occur with equal frequency (1/4 each).
• Every base has an equal rate of substitution, α, to each of the other three bases.
• Substitutions occur randomly and independently at each site, at a constant rate over time (this is the “molecular clock” assumption in its simplest form).
• Since there are three other bases that a given base can change into, the total rate at which any base leaves its current state is 3α.

2.2 Rate Matrix of the JC Model

The substitution process can be represented using a rate matrix, where the off-diagonal elements represent the rate of change from one base (row) to another base (column), and the diagonal elements are negative and equal to minus the sum of the other elements in that row (representing the total rate of leaving that state):

2.3 Derivation of the Distance Formula

To derive the relationship between the observed difference D and time t, we consider a single site starting in state A and track the probability P_AA(t) that it is still in state A after time t.

Over a very small time interval δt, the site can be in state A at time (t + δt) in two ways: either it was already in A at time t and did not change (probability 1 − 3αδt), or it was in some other base (C, G, or T) at time t and changed to A during δt (probability αδt for each). This gives the relation:

P_AA(t + δt) = αδt [P_AC(t) + P_AG(t) + P_AT(t)] + (1 − 3αδt) P_AA(t)

Taking the limit as δt → 0 converts this into a differential equation:

dP_AA/dt = α (P_AC + P_AG + P_AT) − 3α P_AA

Since the probabilities of being in any of the four states must sum to 1, we know P_AC + P_AG + P_AT = 1 − P_AA. Substituting this in simplifies the equation to one involving a single unknown function:

dP_AA/dt = −4α P_AA + α

This is a standard first-order linear differential equation. Using the trial solution P_AA(t) = A·e^(−4αt) + B and applying the initial condition P_AA(0) = 1 (the site starts as A with certainty), the constants work out to B = 1/4 and A = 3/4. Hence:

P_AA(t) = (3/4) e^(−4αt) + 1/4

By the symmetry of the model, P_GG(t), P_CC(t), and P_TT(t) are all equal to P_AA(t). Similarly, P_AC(t), P_AG(t), and P_AT(t) are all equal to one another, and together they must equal 1 − P_AA(t). Therefore:

P_AC(t) = 1/4 − (1/4) e^(−4αt)

Both P_AA(t) and P_AC(t) tend to 1/4 as t becomes large — after a long time, a site is equally likely to be in any of the four states, regardless of its starting base.

Now, when comparing two present-day sequences that diverged from a common ancestor a time t ago, the total time separating the two sequences (added along both branches) is 2t. The observed proportion of different sites, D, is the probability that the base has ended up different from the ancestral base, considering there are three ways it could differ:

D = 3 × P_AC(2t) = 3/4 − 3/4 e^(−8αt)

This is the key relationship between the observed difference D and time:

D(t) = 3/4 − (3/4) e^(−8αt) … (Equation for D)

When t is small, expanding the exponential shows that D increases approximately linearly with time (D ≈ 6αt), confirming the expected behaviour for short divergence times. When t is large, the exponential term vanishes and D tends to its saturation value of 3/4, consistent with the random-sequence argument discussed earlier.

2.4 The Jukes–Cantor Distance Formula

The evolutionary distance d, defined as the mean number of substitutions per site, is obtained by multiplying the rate of change (3α) by the total time available for change (2t):

d = 2 × 3α × t = 6αt

Both the formula for d and the formula for D are functions of the product αt, so αt can be eliminated between the two equations. Rearranging the D(t) equation gives:

ln(1 − 4D/3) = −8αt

Comparing this with d = 6αt gives the final, and most important, result of this section — the Jukes–Cantor distance formula:

d = −(3/4) ln(1 − 4D/3)

This equation allows the evolutionary distance d to be calculated purely from the observable quantity D, without needing to know the actual values of α or t individually.

2.5 Interpretation of the JC Distance

• When D is small, d ≈ D, since there is little chance of multiple substitutions at the same site.
• When D is large (approaching 1/2 or more), d becomes substantially greater than D, because a high proportion of sites have undergone more than one substitution, most of which go undetected.
• As D approaches its theoretical maximum of 3/4, d tends to infinity — this reflects the fact that once sequences are fully randomized with respect to one another, no finite estimate of the number of substitutions is possible from comparing them.
• The Jukes–Cantor distance is thus always greater than or equal to the simple observed difference D, and the gap between them widens as sequences become more diverged.

3. More Complex Models of DNA Sequence Evolution

The JC model, while foundational, makes simplifying assumptions that are rarely true for real biological sequences — namely that all substitutions occur at equal rates and all bases occur with equal frequency. More realistic models relax these assumptions.

3.1 Kimura Two-Parameter (K2P) Model

It is commonly observed that transitions (purine↔purine: A↔G, or pyrimidine↔pyrimidine: C↔T) occur more frequently in real sequences than transversions (purine↔pyrimidine changes). The Kimura two-parameter model, proposed by Kimura in 1980/1983, accounts for this by introducing two separate rate parameters:

• α — the rate of transitions
• β — the rate of transversions

Typically α > β, reflecting the empirical observation that transitions are favoured over transversions. The corresponding rate matrix has α for transition pairs and β for the (twice as numerous) transversion pairs, with diagonal elements equal to −(α + 2β).

Distance formula for K2P: If S is the observed fraction of sites differing by a transition and V is the fraction differing by a transversion (so D = S + V), the estimated number of substitutions per site is:

d = −(1/2) ln(1 − 2S − V) − (1/4) ln(1 − 2V)

Worked example: For a pair of sequences with S = 0.2 and V = 0.1 (so D = 0.3), the K2P distance works out to d ≈ 0.402, whereas the simpler JC distance for the same D gives d ≈ 0.383. This illustrates an important point: the estimated number of substitutions depends on which evolutionary model is used, and since tree-building methods can be sensitive to small changes in input distances, choosing a model appropriate to the data is important for reliable results.

3.2 HKY Model (Hasegawa, Kishino, and Yano, 1985)

Real sequences frequently show unequal base frequencies (denoted π_A, π_G, π_C, π_T), which neither the JC nor the K2P model accounts for. The HKY model incorporates these base frequencies directly into the rate matrix, so that the rate of substitution towards a given base is proportional to that base’s equilibrium frequency, in addition to the transition rate α and transversion rate β.

• A key property of the HKY model is that the probability P_ij(t) of being in state j after a long time tends to the equilibrium frequency π_j of base j, regardless of the starting base i.
• When all four base frequencies are equal (1/4 each), the HKY model mathematically reduces to the K2P model, showing that HKY is a generalisation of K2P that additionally allows for base composition bias.

3.3 General Reversible (GR) Model

The general reversible model is the most parameter-rich model in this family that still obeys time reversibility (see below). It has:

• Four base frequency parameters (as in HKY), and
• Six independent rate parameters, one for each distinct pair of bases (α_AG, α_AC, α_AT, α_GC, α_GT, α_CT).

Because it has the largest number of free parameters, the GR model is, in principle, capable of describing the substitution process in real sequences more accurately than the simpler models. However, this comes at the cost of requiring more data to estimate the additional parameters reliably. Model parameters such as these are typically estimated by fitting the model to real sequence data using the maximum-likelihood criterion.

3.4 Time Reversibility

All of the models described above (JC, K2P, HKY, GR) share an important mathematical property called time reversibility:

π_i r_ij = π_j r_ji

This states that, for any pair of bases i and j, the number of substitutions per site per unit time occurring in the forward direction (i → j) equals the number occurring in the reverse direction (j → i). A consequence of time reversibility is that the base frequencies remain constant, on average, over time, and that the equilibrium probability relationships hold irrespective of the direction of comparison. Time reversibility is an important simplifying assumption used almost universally in phylogenetic methods, because it means the likelihood of a tree does not depend on where the root is placed.

4. Variation of Substitution Rates Between Sites

All the models discussed so far assume that every site in a sequence evolves at the same rate. In reality, this is almost never true, because sites that are structurally or functionally important (for example, in a protein-coding region or an active site) tend to evolve more slowly than less constrained sites, due to purifying selection.

Figure. Relationship between evolutionary distance (d) and the observed fraction of differing sites (D) based on the Jukes–Cantor model. The graph compares three evolutionary models: (i) a uniform substitution rate across all sites, (ii) a model with 25% invariant sites (f = 0.25), and (iii) a model with rate variation among sites described by a gamma distribution (α = 1). The dashed line represents the uncorrected evolutionary distance, where d = D.

4.1 Invariant Sites Model

In the simplest extension, it is assumed that a certain fraction, f, of sites are completely invariant (unable to change at all due to very strong selective constraint), while the remaining (1 − f) fraction evolve according to the standard JC model. Under this assumption, the corrected evolutionary distance becomes:

d = −[3(1 − f)/4] × ln[1 − 4D / (3(1 − f))]

This reduces to the standard JC distance formula when f tends to zero, as expected. Importantly, even a relatively small fraction of invariant sites can make a considerable difference to the estimated distance — the presence of invariant sites means that all the observed change is concentrated in a smaller proportion of the sequence, so the true number of substitutions per variable site is higher than the whole-sequence average would suggest. Note that in this formula, D and d still refer to the whole molecule (invariant and variable sites combined), while the average number of substitutions occurring specifically at the variable sites is given by d/(1 − f).

4.2 Gamma-Distributed Rates Model

A more general and more commonly used approach is to model the distribution of relative rates across sites using a gamma distribution:

f(r, a) = constant × r^(a−1) × e^(−ar)

Here, r is the rate of substitution at a given site relative to the average rate across the whole sequence, and a is a shape parameter that controls the degree of rate variability:

• When a is large, the distribution becomes sharply peaked around r = 1, meaning nearly all sites evolve at essentially the same (average) rate — i.e., little rate variation.
• When a = 1, the distribution takes the form of a simple exponential.
• When a is small (less than 1), the distribution has a high probability near r = 0 (many nearly invariant sites) combined with a long tail at high r (a smaller number of rapidly evolving sites) — i.e., strong rate heterogeneity.

Distance formula with gamma-distributed rates (JC + Γ model):

d = (3a/4) × [ (1 − 4D/3)^(−1/a) − 1 ]

Both the invariant-sites correction and the gamma-distribution correction give a higher estimate of the true number of substitutions per site, for any given observed D, compared with the basic JC model that assumes uniform rates. This means that ignoring rate variation among sites leads to a systematic underestimation of evolutionary distance. The effect is small when overall divergence is low but becomes very large at higher divergence levels — a point of considerable practical importance when using molecular data to estimate divergence times between species (for instance, human–ape divergence dates are sensitive to how rate variation is modelled).

The gamma distribution is popular because it introduces only one additional parameter (a) while being flexible enough to represent a wide range of rate-variation patterns, from almost uniform to highly heterogeneous. The best-fitting value of a for a real dataset is usually estimated using maximum likelihood.

5. Applications and Significance of Evolutionary Models

Quantitative evolutionary distance measures derived from these models are essential in several major areas of molecular evolution and bioinformatics:

• Phylogenetic tree construction: Distance-matrix methods of building trees require pairwise evolutionary distances between all sequences in a dataset; branch lengths on such trees are typically drawn proportional to evolutionary distance, allowing visual identification of lineages with the greatest amount of change.
• Maximum-likelihood phylogenetics: With an explicit quantitative model of substitution, it becomes possible to calculate the likelihood that a given set of sequences would have evolved along any particular proposed tree topology, allowing alternative trees to be statistically compared and the most likely one selected.
• Sequence alignment scoring systems: Evolutionary models of amino acid substitution are used to derive scoring matrices (such as those used in pairwise and multiple sequence alignment), assigning higher scores to pairs of residues that are expected to substitute for one another frequently during evolution, so that the optimally scoring alignment corresponds to the evolutionarily most probable alignment.
• Estimating divergence times and rates: Because the substitution rate depends on mutation, natural selection, and random genetic drift acting together (not simply the mutation rate), evolutionary distances calculated from these models reflect the outcome of population-level fixation processes, and can be used, along with a calibrated molecular clock, to estimate approximate divergence times between species.

Conclusion

The simple observed difference between two sequences, D, is not a reliable measure of evolutionary distance because it is neither linear with time nor additive across branches of a tree, owing to the possibility of multiple, superimposed, or reversed substitutions at the same site. The Jukes–Cantor model resolves this by defining a corrected distance, d, based on an explicit model of random substitution at a constant rate, giving the well-known relationship d = −(3/4) ln(1 − 4D/3).

More elaborate models — Kimura’s two-parameter model, the HKY model, and the general reversible model — progressively relax the restrictive assumptions of equal substitution rates and equal base frequencies, providing a more realistic description of real sequence data at the cost of additional parameters.

Further realism is added by allowing substitution rates to vary between sites, through either an invariant-sites correction or a gamma-distributed rates model, both of which correct for the systematic underestimation of evolutionary distance that would otherwise occur. Together, these models form the mathematical foundation of modern molecular phylogenetics and sequence alignment methodology.