EMBL-European Bioinformatics Institute, Hinxton, United Kingdom
Correspondence: E-mail: goldman{at}ebi.ac.uk.
![]() |
Abstract |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Key Words: amino acid replacement Dayhoff matrix Markov models phylogenetic inference protein evolution
![]() |
Introduction |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() | (1) |
In this article we summarize different methods that have been used to construct IRMs from the probability matrices of DSO78, identify their implementations in standard phylogenetic software packages, and indicate ways in which they lead to different and potentially inaccurate IRMs. We describe two simpler methods for deriving an IRM from information such as that published by Dayhoff and collaborators. We then compare the performance of all of these implementations using a test set of 200 protein domain families.
The Dayhoff model and other similar models (e.g., the JTT model of Jones, Taylor, and Thornton 1992) have been very influential in molecular phylogenetics, database searching, and other fields, and they continue to be widely and regularly used. It is important to have a complete understanding of how these models can be accurately derived from raw data. In the interests of (1) remaining faithful to the information collected and published by Dayhoff and colleagues and others, (2) keeping computations as simple and accurate as possible, and (3) facilitating agreement among scientists using different implementations of supposedly identical methods, we propose a standardization of these IRMs in phylogenetic software.
![]() |
Methods |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
The continuous-time Markov process is a stochastic model in which (P(t))ij gives the probability that amino acid i will change to amino acid j at any single site after any time t 0. Since there are 20 amino acids, i and j take the values 1, ... ,20. The 20 x 20 probability matrix is calculated as P(t) = exp(tQ) (eq. 1), where the matrix Q is independent of time in the Markov processes typically used in molecular phylogenetics. Q has off-diagonal entries qij,i
j equal to the instantaneous rates of replacement of i by j, and diagonal entries defined by the mathematical requirement that each row sum is 0.
Associated with a matrix Q are equilibrium frequencies of the 20 amino acids, denoted fi, and mutabilities, mi, defined as the rate at which amino acid i changes to any other amino acid: mi = j
i qij. Typically in phylogenetic applications, Q is normalized so that the mean rate of replacement at equilibrium (
i
j
i fiqij or
i fimi) is 1, meaning that times t are measured in units of expected numbers of changes per site. Here, we simplify our notation by omitting the trivial multiplicative constants that achieve this normalization. More detailed descriptions of Markov models in molecular phylogenetics are given by Liò and Goldman (1998), Felsenstein (2003), and Thorne and Goldman (2003).
Dayhoff and colleagues devised a method to estimate P(t) that relied on the comparison of closely related pairs of aligned sequences. They selected pairs of sequences that were 85% identical and counted the differences between them. In general these differences underestimate the actual numbers of changes because when multiple changes occur at a single site they are observed as at most single replacements. However, if sufficiently closely related sequences are used, then the probability of these "multiple hits" is reduced. Dayhoff and colleagues assumed that at 85% sequence identity there are no multiple hits (see also Jones, Taylor, and Thornton 1992; Goldman, Thorne, and Jones 1996, 1998). The cost of this 85% rule is some loss in accuracy and inefficient use of data, because divergent sequence pairs have to be discarded.
Advances in methodology and computer speed have made it possible to estimate Markov models of amino acid replacement without any constraints on the levels of divergence between the sequences (e.g., Adachi and Hasegawa 1996; Yang, Nielsen, and Hasegawa 1998; Adachi et al. 2000; Müller and Vingron 2000; Whelan and Goldman 2001; Devauchelle et al. 2001; Dimmic et al. 2002; Veerassamy, Smith, and Tillier 2003). While these methods may now give superior results, less sophisticated methods based on simple counts of changes under the 85% rule are still used (e.g., Goldman, Thorne, and Jones 1996, 1998) and were the basis for two of the most widely used evolutionary models in molecular phylogenetics (DSO78; Jones, Taylor, and Thornton 1992). The rest of our analysis concentrates solely on cases such as these, where differences between sequence pairs can be counted and assumed to accurately represent actual evolutionary events.
We write the observed number of occurrences of sites in all aligned sequence pairs with amino acids i in one sequence and j in the other as nij. The data collection technique of Dayhoff and colleagues leads to the relationships nij = nji, because the direction of an evolutionary change cannot be determined from the observation of two contemporary sequences. (As a consequence, equilibrium is assumed and resulting models must be reversible; see Liò and Goldman 1998.) Because we assume the nij accurately represent all evolutionary events, an intuitive estimate for the rate of change of i to j (j i) is given by the number of such events as a proportion of all observations of i:
![]() | (2) |
![]() | (3) |
![]() | (4) |
Equation (2) represents a simple way to estimate the IRM Q directly from counts of observed differences between closely related sequence pairs and was used by Goldman, Thorne, and Jones (1996, 1998). Applications of the models of Dayhoff and colleagues and Jones, Taylor, and Thornton (1992) in molecular phylogentics have, however, used a different approach, estimating Q not directly but via corresponding probability matrices P(t). We proceed by describing these methods.
The Dayhoff Model
Dayhoff and colleagues published probability matrices based on counts of sequence differences, frequencies, and mutabilities. These articles include the values of the counts nij, but in an incomplete manner: the numbers of positions containing amino acid i in both sequences, nii, are omitted.
The mutabilities mi were calculated slightly differently from equation (3) above. The data from individual sequence pairs s = 1, ..., S were considered separately, and the estimator
![]() | (5) |
Dayhoff and colleagues used much the same idea to estimate the frequencies fi of the amino acids, combining the data for each comparison s to derive the following estimator:
![]() | (6) |
Dayhoff and colleagues described an estimator PD = (pD)ij for probability matrices as follows:
![]() | (7) |
![]() | (8) |
To create the PAM1 matrix, Dayhoff and collaborators chose c so that is 1 observed mutation per 100 sites. Thus
![]() | (9) |
The "time" at PAM1 is defined by = 1 observed amino acid change per 100 amino acids, i.e., excluding changes that are unobservable because of multiple hits. Therefore, PAM1 does not correspond to a probability matrix P(0.01) calculated according to equation (1), but to some P(0.01 +
), where
is the extra time needed to account for unobserved substitutions. For small times t, the difference is negligible (e.g., PAM1 approximates P(0.01 + 6.85 x 105) when using the IRM derived with the DCMut method described below); even for larger times, the difference is still small (PAM250 is close to P(2.52) under the DCMut model).
Dayhoff and colleagues' estimator embodied in equations (7)(9) makes two assumptions: that the data from which the observed counts nij are derived are taken to be sufficiently closely related to exclude multiple hits (the 85% rule), and that there will be no multiple hits in the probability model defined by PD. The latter assumption is approximately true for = 0.01, and becomes more accurate as
decreases. For larger values of
it is less accurate; if
is large enough, the resulting matrix PD will not even be a valid probability matrix (see fig. 1).
|
![]() | (10) |
This method is appropriate for recovering Q if P is generated according to P(t) = exp(tQ) for any specific t. However, Dayhoff and colleagues' approach (equations 79) generates a matrix PD() which is only an approximation to P(t) because multiple hits occurring over the time period corresponding to
are neglected. (Indeed, the PAM1 matrix cannot be generated as exp(t*Q*) for any valid IRM Q* and time t*
0proof not shown.) This approximation is increasingly poor as
increases (fig. 1). Kishino, Miyata, and Hasegawa (1990) adopted
= 0.01, to generate PAM1 as in DSO78. It turns out that this value of
is not small enough. Using this eigen-decomposition method, we have calculated the IRMs Q for
in the range [3 x 107, 3 x 101], and figure 2 shows that these can suffer from two convergence problems. On the one-hand, if
is too large, elements qij may not yet be converged. On the other hand, choosing
too small may cause the numerical eigen-analysis that finds the
i and ui to become unstable. Thus it is necessary to check convergence for each of the 19 x 20 = 380 values of qij,i
j. We find that many of the qij computed by the method of Kishino, Miyata, and Hasegawa (1990) using
= 0.01 do not attain converged values. Therefore, the IRM derived this way does not give a fully accurate representation of the amino acid replacement data collected by DSO78.
|
We call the first method Direct Computation with Mutabilities (DCMut), because it uses only the observed changes and the mutabilities. Rearranging equations (2) and (3):
![]() | (11) |
The second method, Direct Computation with Frequencies (DCFreq), relies only on the observed changes and the frequencies and is defined by rearranging Equations (2) and (4):
![]() | (12) |
![]() | (13) |
Note that both the DCMut and DCFreq derivations of IRMs require neither the consideration of any limit
0 nor matrix eigen-analysis. They do, however, still require that the sequence pairs from which the nij are derived should be closely related (e.g., satisfying the 85% rule). Although there are many ways to calculate P(t) = exp(tQ) (Moler and Van Loan 2003), the most popular method in molecular phylogenetics uses eigen-decomposition (Liò and Goldman 1998) and thus even with the DCMut and DCFreq approaches to deriving Q we do not avoid eigen-analysis in its use.
Had Dayhoff and colleagues calculated mutabilities and frequencies according to equations (3) and (4), the IRMs given by equations (2), (11), and (13) would be identical. Likewise, if approaches zero then the Kishino, Miyata, and Hasegawa (1990) method applied to PD (eq. 10) would also converge to this rate matrix. But because Dayhoff and colleagues estimated the mutabilities and frequencies incorporating the weighting factors described above (eq. 5 and eq. 6), we expect slightly different IRMs. Figure 2 shows the element q21 (Arg
Asn) calculated according to equations (10), (11), and (13). Although very close, the values are not identical. To facilitate agreement among scientists wishing to use supposedly identical models, we would like to be able to suggest one standard implementation of the model of DSO78 for molecular phylogenetics. In the Results and Discussion we perform a comparison of the different IRMs' performance in practice and propose a new standard.
![]() |
Results and Discussion |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
KMH refers to the IRM calculated according to the eigen-decomposition method described by Kishino, Miyata, and Hasegawa (1990) applied to PAM1 of DSO78. No widely used software packages currently use exactly this implementation. We find that a number of off-diagonal elements of Q are < 0, meaning that this is not strictly a valid IRM (see the earlier section From Probability to Rate Matrix via Eigen-Decomposition). However, this implementation ran without problems in codeml, and so we made no further alterations.
The IRM for the Dayhoff model distributed with the PAML package, version 3.13d, (Yang 1997) is calculated using the Kishino, Miyata, and Hasegawa (1990) eigen-decomposition method applied to PAM1 (as for KMH above). The rate matrix Q = (qij) is then uniquely decomposed into qij = fj x sij (see, e.g., Whelan and Goldman 2001), where the fj are the equilibrium frequencies derived from this Q (and which thus differ minutely from those published by DSO78). The exchangeabilities sij are multiplied by approximately 100, rounded, and stored as integers. All except one of the sij < 0 become equal to 0 after this rounding. Following Yang (1997), we replace the rounded value of 1 for sGlu Arg with + 1. We refer to this version of the Dayhoff model as the Paml implemention; it is also used in the MrBayes v. 3.0 (Ronquist and Huelsenbeck 2003) and Phyml v. 2.4 (Guindon and Gascuel 2003) programs.
The proml program distributed with version 3.6 of Felsenstein's PHYLIP package (Felsenstein 2002) includes an implementation of the Dayhoff model based on an IRM QF defined by
![]() | (14) |
Early versions of the protml program of the MOLPHY package (Adachi and Hasegawa 1992) used the Kishino, Miyata, and Hasegawa (1990) eigen-decomposition of PAM1; later versions (e.g., v. 2.3b3) use the DCFreq approach, suggested to the authors of MOLPHY by Korbinian Strimmer and Arndt von Haeseler (Korbinian Strimmer, pers. comm.). However, the counts nij used were changed from those in DSO78: where zeros occurred, they were substituted by scaled values taken from Jones, Taylor, and Thornton (1992). This makes the later MOLPHY implementation, also used in the TREE-PUZZLE v. 5.2 package (Schmidt et al. 2002) and in the sequence simulation program PSeq-Gen v. 1.1 (Grassly, Adachi, and Rambaut 1997), a hybrid between the Dayhoff and JTT models.
The DCMut method has never been used before. We implemented the Dayhoff model according to this method, using data from DSO78. We have also calculated a DCFreq implementation, again using exactly the data from DSO78.
To get an idea of the impact of the different versions of the Dayhoff model on phylogenetic analysis, we performed a small test. We calculated the maximum likelihoods (Felsenstein 2003) for 200 protein families and their associated tree topologies taken from release 7.6 of the PANDIT database (Whelan, de Bakker, and Goldman 2003) using the implementations described above. Equilibrium frequencies were either taken from the appropriate implementation of the Dayhoff model (analyses labeled KMH, Paml, Proml, Molphy, DCMut, DCFreq), or were estimated separately for each protein family and incorporated using the "+F" method of Cao et al. (1994; see also Thorne and Goldman 2003; analyses KMH+F, Paml+F, Proml+F, Molphy+F, DCMut+F, DCFreq+F). (Note that it is not currently possible to use the Proml+F model in the proml program of Felsenstein's PHYLIP package. Note also the importance of re-normalizing Q so that the mean rate of change at equilibrium equals 1 when using the +F method, particularly if different software is used for simulating sequence data and then analyzing that simulated data.) Table 1 summarizes the results.
|
By modern standards not much data was available to Dayhoff and colleagues, and their published values of nij include a number of zeros. It seems likely that these give significant underestimates of the frequencies of replacement between certain amino acids, and we attribute the Molphy implementation's success to its hybrid nature, incorporating information from the counts published by Jones, Taylor, and Thornton (1992) to redress this underestimation. Among versions based solely on the information published by DSO78, the DCMut and DCFreq versions have advantages in computational ease and seem to give a small advantage in terms of maximum likelihood scores over the test data sets studied. The versions based on eigen-decomposition of PAM1 (i.e., with = 0.01, which does not guarantee convergence) give the worst performance and seem to have no advantages. It is interesting to note that the Paml implementation fares better than KMH; evidently the variations introduced by the rounding procedure in the Paml version, in combination with the usage of a non-converged eigen-decomposition method and the sparse data published by DSO78, give a model that is no worse, and may even be better, than the version without rounding.
![]() |
Conclusions |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
We have summarized different approaches (including two new ones) to calculating an instantaneous rate matrix from the (incomplete) information given by Dayhoff and colleagues. In practice, the differences are small but may be non-trivial. All the implementations studied are valid models of protein sequence evolution, and they may be applied to any protein sequence alignments. Nevertheless, while this and similar models remain of interest in molecular phylogenetics and other fields and for the sake of consistency, particularly among investigators developing models and software implementations, we suggest that it is of value to identify a "standard" implementation for the model of Dayhoff, Schwarz, and Orcutt (1978).
Although the Molphy implementation (Adachi and Hasegawa 1992) performs well, it uses a hybrid data set based on both Dayhoff, Schwarz, and Orcutt (1978) and Jones, Taylor, and Thornton (1992). The Proml implementation also uses the Dayhoff, Schwarz, and Orcutt (1978) data in a modified form. The Proml, Paml, and KMH implementations perform least well in our maximum likelihood implementation experiment, and the generation of their instantaneous rate matrices is based on relatively complex eigen-decomposition and, in the cases of Paml and KMH, exhibits convergence problems. Therefore, looking for standardization based only on the data published by Dayhoff, Schwarz, and Orcutt (1978), we propose the adoption of the DCMut matrix, which is straightforward to calculate and performed well in our test. Adopting DCMut also means that no decision has to be made regarding the fact that, because of rounding errors, the equilibrium frequencies published by Dayhoff, Schwarz, and Orcutt (1978) do not sum to 1. We do not suggest that DCMut is a better model than any other implementation described here, but only that it is a reasonable choice when a single model is required for comparing results across different analyses, computer programs, or data sets.
Jones, Taylor, and Thornton (1992) used much the same methodology as Dayhoff and colleagues, but based on a larger sequence database. Again, only a probability matrix, mutabilities, frequencies and incomplete counts are provided in their original article. The JTT probability matrix has also been used to compute an IRM for use in phylogenetics and has been implemented (again, in a variety of ways) in the MOLPHY, MrBayes, PAML, PHYLIP, Phyml, PSeq-Gen, and TREE-PUZZLE software. We again suggest that an implementation based on the DCMut method be adopted as standard.3
We hope to persuade software developers to agree to include our DCMut implementations of the Dayhoff and JTT models into their programs, and we are happy to discuss providing the necessary data in whatever forms required.4
![]() |
Acknowledgements |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() |
Footnotes |
---|
2 Files prepared following the pattern of the distributed files (e.g.) dayhoff.dat and jones.dat are available from http://www.ebi.ac.uk/goldman/dayhoff.
3 A file containing this matrix, suitable for use in the codeml program of the PAML package, is available from http://www.ebi.ac.uk/goldman/dayhoff.
4 The authors of MrBayes (Ronquist and Huelsenbeck 2003), PAML (Yang 1997), PAUP* (Swofford 2002), Phyml (Guindon and Gascuel 2003), Seq-Gen (Rambaut and Grassly 1997; the successor to PSeqGen) and TREE-PUZZLE (Schmidt et al. 2002) have agreed to incorporate the DCMut versions of the Dayhoff and JTT IRMs in future releases of their software.
Arndt von Haeseler, Associate Editor
![]() |
References |
---|
![]() ![]() ![]() ![]() ![]() ![]() ![]() ![]() |
---|
Adachi, J., and M. Hasegawa. 1992. MOLPHY version 2.3: Programs for Molecular Phylogenetics Based on Maximum Likelihood. Computer Science Monographs 28, Institute of Statistical Mathematics, Tokyo. http://www.ism.ac.jp/software/ismlib/softother.e.html#molphy
. 1996. Model of amino acid substitution in proteins encoded by mitochondrial DNA. J. Mol. Evol. 42:459468.[ISI][Medline]
Adachi, J., P. J. Waddell, W. Martin, and M. Hasegawa. 2000. Plastid genome phylogeny and a model of amino acid substitution for proteins encoded by chloroplast DNA. J. Mol. Evol. 50:348358.[ISI][Medline]
Cao, Y., J. Adachi, A. Janke, S. Pääbo, and M. Hasegawa. 1994. Phylogenetic relationships among eutherian orders estimated from inferred sequences of mitochondrial proteins: instability of a tree based on a single gene. J. Mol. Evol. 39:519527.[ISI][Medline]
Dayhoff, M. O., R. V. Eck, and C. M. Park. 1972. A model of evolutionary change in proteins. Pp. 8999 in M. O. Dayhoff, ed., Atlas of Protein Sequence and Structure Vol. 5. National Biomedical Research Foundation, Washington, D.C.
Dayhoff, M. O., R. M. Schwartz, and B. C. Orcutt. 1978. A model of evolutionary change in proteins. Pp. 345352 in M. O. Dayhoff, ed., Atlas of Protein Sequence and Structure Vol. 5, suppl. 3. National Biomedical Research Foundation, Washington, D.C.
Devauchelle, C., A. Grossmann, A. Hénaut, M. Holschneider, M. Monnerot, J. L. Riesler, and B. Torrésani. 2001. Rate matrices for analyzing large families of protein sequences. J. Comp. Biol. 8:381399.[CrossRef][ISI]
Dimmic, M. W., J. S. Rest, D. P. Mindell, and R. A. Goldstein. 2002. rtREV: an amino acid substitution matrix for inference of retrovirus and reverse transcriptase phylogeny. J. Mol. Evol. 55:6573.[CrossRef][ISI][Medline]
Felsenstein, J. 1996. Inferring phylogenies from protein sequences by parsimony, distance, and likelihood methods. Methods Enzymol. 266:418427.[ISI][Medline]
. 2002. PHYLIP (Phylogeny Inference Package) Version 3.6a. Department of Genome Sciences, University of Washington, Seattle, Wash. http://evolution.genetics.washington.edu/phylip.html
. 2003. Inferring phylogenies. Sinauer Associates, Sunderland, Mass.
Goldman, N., J. L. Thorne, and D. T. Jones. 1996. Using evolutionary trees in protein secondary structure prediction and other comparative sequence analysis. J. Mol. Biol. 263:196208.[CrossRef][ISI][Medline]
. 1998. Assessing the impact of secondary structure and solvent accessibility on protein evolution. Genetics 149:445458.
Grassly, N. C., J. Adachi, and A. Rambaut. 1997. PSeq-Gen: an application for the Monte Carlo simulation of protein sequence evolution along phylogenetic trees. CABIOS 13:559560. http://evolve.zoo.ox.ac.uk/software.html?id=pseqgen[Medline]
Guindon, S., and O. Gascuel. 2003. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst. Biol. 52:696704. http://atgc.lirmm.fr/phyml[CrossRef][ISI][Medline]
International Chimpanzee Chromosome 22 Consortium. 2004. DNA sequence and comparative analysis of chimpanzee chromosome 22. Nature 429:382388.[CrossRef][ISI][Medline]
International Human Genome Sequencing Consortium. 2001. Initial sequencing and analysis of the human genome. Nature 409:860921.[CrossRef][ISI][Medline]
International SNP Map Working Group. 2001. A map of human genome sequence variation containing 1.42 million single nucleotide polymorphisms. Nature 409:928933.[CrossRef][ISI][Medline]
Jones, D. T., W. R. Taylor, and J. M. Thornton. 1992. The rapid generation of mutation data matrices from protein sequences. CABIOS 8:275282.[Medline]
Kishino, H., T. Miyata, and M. Hasegawa. 1990. Maximum likelihood inference of protein phylogeny and the origin of chloroplasts. J. Mol. Evol. 31:151160.[ISI]
Liò, P., and N. Goldman. 1998. Models of molecular evolution and phylogeny. Genome Res. 8:12331244.
Moler, C., and C. Van Loan. 2003. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Rev. 45:349.[ISI]
Mouse Genome Sequencing Consortium. 2002. Initial sequencing and comparative analysis of the mouse genome. Nature 420:520562.[CrossRef][ISI][Medline]
Müller, T., and M. Vingron. 2000. Modeling amino acid replacement. J. Comp. Biol. 7:761776.[CrossRef][ISI]
Rambaut, A., and N. C. Grassly. 1997. Seq-Gen: an application for the Monte Carlo simulation of DNA sequence evolution along phylogenetic trees. CABIOS 13:235238. http://evolve.zoo.ox.ac.uk/software.html?id=seqgen[Medline]
Rat Genome Sequencing Project Consortium. 2004. Genome sequence of the Brown Norway rat yields insights into mammalian evolution. Nature 428:493521.[CrossRef][ISI]
Ronquist, F., and J. P. Huelsenbeck. 2003. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19:15721574. http://morphbank.ebc.uu.se/mrbayes3
Schmidt, H. A., K. Strimmer, M. Vingron, and A. von Haeseler. 2002. Tree-Puzzle: maximum likelihood phylogenetic analysis using quartets and parallel computing. Bioinformatics 18:502504. http://www.tree-puzzle.de
Swofford, D. L. 2002. PAUP*. *Phylogenetic analysis using parsimony (and other methods) version 4. Sinauer Associates, Sunderland, Mass. http://paup.csit.fsu.edu
Thorne, J. L., and N. Goldman. 2003. Probabilistic models for the study of protein evolution. Pp. 209226 in D. J. Balding, M. Bishop, and C. Cannings, eds. Handbook of Statistical Genetics, 2nd Ed. Wiley, Chichester.
Veerassamy, S., A. Smith, and E. R. M. Tillier. 2003. A transition probability model for amino acid substitutions from Blocks. J. Comp. Biol. 10:9971010.[CrossRef][ISI]
Whelan, S., P. I. W. de Bakker, and N. Goldman. 2003. Pandit: a database of protein and associated nucleotide domains with inferred trees. Bioinformatics 19:15561563. http://www.ebi.ac.uk/goldman-srv/pandit
Whelan, S., and N. Goldman. 2001. A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach. Mol. Biol. Evol. 18:691699.
Yang, Z. 1997. PAML: a program package for phylogenetic analysis by maximum likelihood. CABIOS 13:555556. http://abacus.gene.ucl.ac.uk/software/paml.html[Medline]
Yang, Z., R. Nielsen, and M. Hasegawa. 1998. Models of amino acid substitution and applications to mitochondrial protein evolution. Mol. Biol. Evol. 15:16001611.