The PAM Model of Protein Sequence Evolution -Exam Notes

Just as models of DNA sequence evolution such as the Jukes–Cantor model are built around a 4 × 4 rate matrix describing substitutions between the four bases, a model of protein sequence evolution must account for substitutions among all twenty amino acids, and therefore requires a 20 × 20 substitution matrix.

The earliest and most widely used substitution matrices for amino acids are known as PAM matrices, where PAM stands for “Point Accepted Mutation.” An accepted mutation refers to a single amino acid change that has spread through a population and become fixed, so that it appears as a fixed difference when present-day sequences are compared.

Although the name “single amino acid substitution” would perhaps describe the concept more accurately, the term PAM is well established in the literature and is used throughout. The derivation of the PAM model presented here follows the original method of Dayhoff, Schwartz, and Orcutt (1978), which remains the conceptual basis for later, larger-scale versions of the model.

1. Counting Amino Acid Substitutions

To construct a PAM matrix, the first requirement is to count, from real sequence data, how often each amino acid has been substituted by every other amino acid. Dayhoff, Schwartz, and Orcutt based their analysis on alignments of 71 families of closely related proteins, where sequences within a family differed from one another by no more than 15%. For each family, an evolutionary tree was constructed using the parsimony method. The principle of parsimony is straightforward: among all possible trees, the one selected is the tree that requires the smallest possible total number of amino acid substitutions to explain the observed sequences.

This procedure can be illustrated using a short protein alignment containing only six different amino acids, chosen here purely to keep the example manageable. The known present-day sequences (labelled A to G) sit at the tips of the tree, while the sequences at the internal nodes (labelled 1 to 5) are not directly observed — they are inferred by a computer program so as to minimise the total number of substitutions required across the whole tree. Each substitution is then labelled on the branch of the tree where it is inferred to have occurred.

Once every substitution on the tree has been identified, the results are assembled into a matrix whose elements, A_ij, record the number of times amino acid i has been substituted by amino acid j. An important subtlety arises here: the parsimony method produces an unrooted tree, so it cannot tell us the direction in which a given substitution occurred. For example, a change between K and T on a particular branch could equally be interpreted as K→T or as T→K, depending on where the root of the tree is placed. For this reason, whenever a substitution between amino acids i and j is observed, one count is added to both A_ij and A_ji — that is, every substitution is counted in both directions. In the six-amino-acid example above, there are seven such substitutions in total, giving rise to the following matrix for the amino acids I, K, L, Q, T, and V:

When Dayhoff and colleagues applied this method to all twenty amino acids across their full dataset, the resulting A_ij matrix contained 1572 substitutions in total. A complication that arises in this method is that, at some internal nodes, more than one amino acid may be equally parsimonious that is, the ancestral sequence is ambiguous. In such cases, fractional substitutions are counted proportionally rather than assigning the full count to a single possibility.

A much larger-scale version of this analysis was later carried out by Jones, Taylor, and Thornton (1992), using all protein sequences available in Swiss-Prot at the time. Rather than working with small, hand-picked protein families, they automatically clustered sequences into groups sharing more than 85% similarity and applied a distance-matrix method to obtain phylogenetic trees for each cluster. Their resulting A_ij matrix contains 59,190 substitutions nearly forty times as many as the original Dayhoff dataset and every possible interchange between pairs of amino acids is observed at least once in this large dataset.

Why Substitutions Are Counted in Both Directions: Time Reversibility

There is a second, more theoretical reason for counting every substitution in both directions, beyond the fact that parsimony trees are unrooted. The aim is to construct a time-reversible model of evolution, exactly as discussed earlier for DNA substitution models. If the frequency of each amino acid, π_i, is assumed to remain constant on average over time, then we should expect to see, on average, an equal number of substitutions in each direction between any given pair of amino acids — that is, A_ij = A_ji on average. This follows because the number of observed substitutions from i to j is proportional to the probability π_i of being in state i, multiplied by the rate r_ij of change from i to j given that the site is currently in state i. By deliberately starting from a symmetric matrix (A_ij set equal to A_ji), the evolutionary model that is subsequently derived from it is guaranteed to satisfy the time-reversibility condition and behave as a proper reversible model.

2. Defining the Evolutionary Model: The PAM1 Matrix

As with the DNA substitution models, P_ij(t) denotes the probability that a site is in state j at time t, given that it was in state i at time 0. The PAM1 matrix is defined as:

M_ij = P_ij(δt)

where δt is a small unit of time referred to as 1 PAM unit. For sufficiently small time intervals, substitution probabilities can be assumed to be directly proportional to substitution rates — that is, the probabilities vary linearly with time when time is small. Our estimate of the substitution rate is therefore taken to be proportional to the number of observed substitutions, A_ij, divided by the total number of times amino acid i appears in the dataset, N_i:

M_ij = λ (A_ij / N_i)   for i ≠ j

where λ is a constant of proportionality still to be determined.

Note on convention: Here M_ij and r_ij represent substitutions from i to j, consistent with the convention used for DNA models, with the way Jones et al. (1992) present their data, and with several more recent authors. The original Dayhoff et al. (1978) paper used the opposite convention, so formulas quoted from that paper may appear with reversed indices here.

The frequency of amino acid i in the dataset is defined as π_i = N_i/N_tot, where N_tot is the total number of amino acids in the dataset. To pin down the value of λ, the convention adopted is that 1 PAM unit is the amount of time in which, on average, 1% of amino acids have changed. Writing out the total fraction of sites that have changed and simplifying gives:

Σ_i π_i Σ_j≠i M_ij = λ A_tot / N_tot = 0.01⇒ λ = 0.01 N_tot / A_tot

where A_tot is the sum of all the elements of the A_ij matrix. To complete the PAM1 matrix, the diagonal elements — the probability that amino acid i remains unchanged over 1 PAM unit — are obtained from:

M_ii = 1 − Σ_j≠i M_ij

In the PAM1 matrix obtained by Jones et al. (1992), values are conventionally multiplied by 100,000 for ease of presentation, so that every row sums to exactly 100,000. All diagonal elements are slightly less than one (i.e., slightly less than 100,000 on this scale), while off-diagonal elements are all very small. The two largest off-diagonal entries in each row are highlighted, since they represent the two fastest substitution partners for that amino acid — for instance, alanine (A) is more likely to be substituted by serine (S) and threonine (T) than by any of the other amino acids. Unlike the underlying A_ij matrix, the M_ij matrix itself is not symmetric, because amino acid frequencies are unequal. Since the model is time-reversible, the relationship π_iM_ij = π_jM_ji holds — for example, for alanine and valine, 0.077 × 0.00193 equals 0.066 × 0.00226, both giving 0.000149.

Relative Mutability

If every amino acid changed at exactly the same rate, all diagonal elements M_ii would equal 0.99, since 1 PAM unit is defined around an average change probability of 1%. In practice, amino acids that change faster than the average will have M_ii below 0.99, while those that change more slowly will have M_ii above 0.99. The total probability of substitution away from amino acid i in one PAM unit is (1 − M_ii); dividing this by the mean substitution probability of 0.01 gives the relative mutability, m_i, of amino acid i:

m_i = (1 − M_ii) / 0.01

A relative mutability of 1 indicates an amino acid that changes at exactly the average rate. Tryptophan (W) is the least mutable amino acid, with m_W = 0.316 — less than one-third of the average rate — while serine (S) is the most mutable, with m_S = 1.452. These mutability values reflect deeper structural and physicochemical properties of the amino acids. It is worth noting that Jones et al. (1992) and Dayhoff et al. (1978) originally used a different scale in which alanine’s mutability was fixed at 100; the scale used here, where the average mutability equals 1, is considered more logical, since there is nothing inherently special about alanine that would justify using it as the reference point.

3. Extrapolating the Model to Higher Evolutionary Distances

The PAM1 matrix is the protein-sequence equivalent of the continuous-time DNA evolutionary models discussed earlier, with one key difference: M_ij represents a transition probability over a small, discrete unit of time, whereas the rate matrices r_ij used for DNA describe rates of change in a continuous-time model. For two different amino acids, the probability of substitution over a small time δt is:

P_ij(δt) = δt × r_ij = M_ijP_ii(δt) = 1 − Σ δt × r_ij = M_ii

Choosing the constant λ described above is equivalent to choosing a value of δt such that there is a 1% probability of an amino acid changing within that time.

Combining Time Steps: Matrix Multiplication

An important point is that the transition probability over a doubled time interval, 2δt, is not simply 2 × M_ij. To calculate it correctly, every possible intermediate amino acid state that could have existed at the midpoint (time δt) must be considered and summed over:

P_ij(2δt) = Σ_k M_ik M_kj

For example, to go from alanine (A) to valine (V) in two steps, the path could be A→A→V, A→V→V, A→R→V, or any of the other seventeen possible intermediate routes through the remaining amino acids.

More generally, for n time steps (a total elapsed time of nδt), all possible amino acids occupying each of the n − 1 intermediate time points must be summed over. This is equivalent to multiplying the M matrix by itself n times and reading off the appropriate element of the resulting matrix, Mⁿ:

P(nδt) = Mⁿ_ij = Σ_k Σ_l … Σ_m M_ik M_kl … M_mj

This completes the derivation of the PAM model. The well-known substitution matrices PAM100, PAM250, and so on, are simply the PAM1 matrix raised to the appropriate power — that is, M multiplied by itself 100 or 250 times respectively. This is straightforward to compute using a computer, and the number of times the matrix has been multiplied by itself is referred to as the PAM distance.

Why the Model Is Built from Closely Related Sequences

The PAM model is deliberately derived starting from groups of closely related sequences, for several important reasons:

• Reliable alignment: when sequences are highly similar, they are easy to align because there are almost no gaps, so the alignment obtained is insensitive to the details of whatever scoring matrix is used at this stage.
• Avoiding circularity: since a proper scoring matrix has not yet been derived at this point in the analysis, it is important that the initial alignments do not depend on the quality of an as-yet-uncalculated scoring system — closely related sequences avoid this problem.
• Reliable tree-building: phylogenetic trees can generally be reconstructed more reliably for closely related sequences than for highly diverged ones.
• Equal weighting in parsimony: the parsimony method treats every type of amino acid substitution as equally likely when selecting the best tree, which again avoids relying on an evolutionary model that has not yet been established.
• Single-substitution assumption: when counting substitutions, it is assumed that no more than one substitution has occurred at any site on a given branch — for instance, a K↔T change is assumed to be a single, direct substitution rather than occurring indirectly via an intermediate route such as K↔Q↔T. This assumption is reasonable for closely related sequences but breaks down for highly diverged sequences, where there is no strong justification for using parsimony at all. Likelihood-based methods of estimating substitution matrices, which do not depend on this assumption, are preferred for more divergent sequence sets.

Log-odds scoring matrix: Alongside the observed substitution counts, a log-odds scoring matrix corresponding to PAM250 can be calculated from the same data as S_ij = 10 log₁₀ R_ij, rounded to the nearest integer. Positive scores indicate amino acid pairs that are more likely to interchange with one another than would be expected purely by chance, and such scores form the basis of substitution scoring systems used in sequence alignment.

4. Box: PAM Distances – Relating Observed Difference to Evolutionary Distance

For two proteins separated by a PAM distance of 1, there has been an average of 0.01 substitutions per site, so the probability that any two aligned sites differ is D = 0.01. For a PAM distance of n, there has been an average of 0.01 × n substitutions per site. To keep this parallel with the DNA models, where evolutionary distance is defined as the number of substitutions per site, the protein evolutionary distance is defined as:

d = 0.01 × n

so that PAM100 corresponds to an average of d = 1 substitution per site. This does not imply that every site has changed exactly once — in reality, some sites will have changed more than once while others will not have changed at all. The fraction of sites expected to differ at a PAM distance of n is given by:

D = Σ_i π_i (1 − Mⁿ_ii)

This expression is simply the probability that the amino acid in the first sequence is i, multiplied by the probability that the corresponding site in the second sequence has not remained i after n PAM units, summed over all twenty amino acids. This is the protein-model counterpart of the Jukes–Cantor distance relationship used for DNA sequences.

The relationship between D and d has the same overall shape as for the DNA models — approximately linear for small D, and diverging steeply as D becomes larger. For completely unrelated (random) proteins, D tends towards approximately 0.89 (this is higher than the 3/4 limit found for DNA, because amino acid frequencies are unequal and there are twenty possible states rather than four). At d = 1 (that is, PAM distance 100), D is approximately 0.57 — in other words, about 43% of residues remain identical even though there has been, on average, one substitution per site across the sequence.

A practical limitation is that, unlike the Jukes–Cantor formula for DNA, this distance relationship cannot easily be inverted — given an observed value of D from two aligned proteins, it is not straightforward to solve directly for d. One practical solution is to use the D-versus-d graph as a calibration curve: for example, an observed D of 0.5 corresponds to a distance of roughly d ≈ 0.8 (a PAM distance of about 80). To avoid relying on a graph, Kimura (1983) proposed a convenient closed-form approximation:

d = − ln(1 − D − 0.2D²)

This formula is extremely close to the curve derived from real protein data over most of the practical range of distances, and is virtually indistinguishable from it for values of D below about 0.75. Because it expresses d directly as a function of the observable quantity D, this is the more convenient “right way round” formula to use in practice.

Maximum-Likelihood Estimation of Protein Distances

Both of the distance formulas above strictly apply only to very long sequences in which amino acid frequencies are assumed to exactly match the average frequencies from the original reference dataset. For two real proteins of finite length, actual amino acid frequencies will deviate somewhat from this average, so a more precise estimate of evolutionary distance uses the maximum-likelihood method. For two sequences of length L, with a_k and b_k denoting the amino acids present at site k in sequences 1 and 2 respectively, the likelihood that sequence 1 evolved into sequence 2 over n PAM units is:

L(n) = Π_k π_ak Mⁿ_{ak, bk}

The maximum-likelihood PAM distance is the value of n that maximises L(n), and the corresponding number of substitutions per site is d = n/100. This approach has been used, for example, to calculate evolutionary distances between a set of hexokinase sequences using the “protdist” program from the Phylip software package, based on an alignment of 393 sites with all gapped regions excluded. Distances calculated this way for real sequence pairs generally lie close to, but not exactly on, the theoretical average curve, since the finite length and specific composition of any individual pair of sequences will cause some deviation from the idealised model.

It is worth noting that many different methods exist for estimating evolutionary distances between proteins, and the details of the specific method used in any software package matter. Distances calculated from evolutionary models based on different underlying sequence datasets will differ somewhat. Even when nominally using the same PAM model, different programs can give slightly different distance estimates — for instance, the Phylip protdist program keeps amino acid frequencies fixed at the values from the original reference dataset, whereas the Tree-Puzzle program re-estimates the amino acid frequencies directly from the sequences being analysed, leading to small discrepancies between the two.

Conclusion

The PAM model extends the logic used for DNA substitution models to the more complex case of protein sequences, where twenty amino acids must be related to one another through a 20 × 20 substitution matrix. The model is built empirically: substitutions are counted from parsimony-based trees of closely related sequence families, assembled into an observed substitution matrix (A_ij), and converted into a small-time-step transition probability matrix, the PAM1 matrix, using the assumption that 1 PAM unit corresponds to an average of 1% amino acid change. Repeated matrix multiplication then extrapolates this short-time model to any desired evolutionary distance, producing matrices such as PAM100 and PAM250 used widely in sequence alignment and phylogenetics. Distance formulas — either the direct D-versus-d relationship, Kimura’s closed-form approximation, or maximum-likelihood estimation — allow the observed similarity between two real protein sequences to be converted into a proper, additive measure of evolutionary distance, analogous to the Jukes–Cantor distance used for DNA.

Reference

Higgs, P. G., & Attwood, T. K. (2005). Bioinformatics and Molecular Evolution. Blackwell Publishing.