Neutral, independent, finite sites models

Most substitution models used to date are neutral, independent, finite sites models.

Neutral

Selection does not operate on the substitutions, and so they are unconstrained.

Independent

Changes in one site do not affect the probability of changes in another site.

Finite Sites

There are finitely many sites, and so over evolution, a single site can be changed multiple times. This means that, for example, if a character has value 0 at time 0 and at time t, it could be that no changes occurred, or that it changed to a 1 and back to a 0, or that it changed to a 1 and back to a 0 and then to a 1 and then back to a 0, and so on.

The molecular clock and the units of time

Different substitution models deal with time differently.

• • It is very common to measure time in substitutions. For example, if one was going to construct a phylogenetic tree while employing a substitution model, one could just measure the distance along the branches of the trees in substitutions. This is convenient because it avoids any question of whether the rate of substitution with respect to the unit of time has changed or not (because by definition the number of substitutions per substitution is one), and it does not need any information about timescales that could be called into question.

  • The molecular clock assumption is also very common, namely that the rate of substitutions with respect to time is constant. This is just a multiplying factor (usually called μ, the number of substitutions per unit time) different from measuring time in substitutions. To carry out this type of analysis, one needs to estimate μ first (which requires knowledge of at least one branch length ahead of time, often a difficult task, which can easily be disputed by others).

  • The assumption of a molecular clock is often unrealistic, especially across long periods of evolution. For example, even though rodents are genetically very similar to primates, they have undergone a much higher number of substitutions in the estimated time since divergence in some regions of the genome^[1]. This could be due to their shorter generation time^[2], higher metabolic rate, increased population structuring, increased rate of speciation, or smaller body size ^[3][4]. When studying events like the Cambrian explosion under a molecular clock assumption, poor concurrence between cladistic and phylogenetic data is often observed. There has been some work on models allowing variable rate of evolution (see for example ^[5] and ^[6]).

Time-reversible models

Most useful substitution models are time-reversible. In terms of substitution models, this simply means that over time, the relative frequencies of each character do not change.

For a time reversible model, there is no assumption that substitutions preferentially change in certain directions over time. For example A -> C -> G is the same as G -> C -> A.

The reason for this is because when an analysis of real biological data is performed, there is generally no access to the sequences of ancestral species, only to the species present today. However, when a model is time-reversible, which species was the ancestral species is irrelevant. Instead, the phylogenetic tree can be rooted along the branch leading to any arbitrary extant species, re-rooted later based on new knowledge, or left unrooted.

A time reversible model satisfies the following property π₁Q₁₂ = π₂Q₂₁ (see below for the definition of this notation)

The mathematics of substitution models

Neutral, independent, finite sites models (assuming a constant rate of evolution) have two parameters, Π, a vector of base (or character) frequencies at time zero (for a time reversible model, this vector is usually referred to as the equilibrium base frequencies, and applies at all times), and the rate matrix, Q, which describes the rate at which bases of one type change into bases of another type, Q_ij for

is the rate at which base i goes to base j. For convenience, the diagonals of the Q matrix are chosen so the rows sum to zero.

The transition matrix function is a function from the branch lengths (in some units of time, possibly in substitutions), to a matrix of conditional probabilities. It is denoted P(t). The entry in the i^th column and the j^th row (P_ij(t)) is the probability, after time t, that there is a base j at a given position, conditional on there being a base i in that position at time 0. When the model is time reversible, this can be performed between any two sequences, even if one is not the ancestor of the other, if you know the total branch length between them.

The asymptotic properties of P_ij(t) are such that

where δ_ij is the Kronecker delta function. That is, there is no change in base composition between a sequence and itself. At the other extreme, or, in other words, as time goes to infinity the probability of finding base j at a position given there was a base i at that position originally goes to the equilibrium probability that there is base j at that position, regardless of the original base.

The transition matrix can be computed from the rate matrix and the equilibrium base frequencies by P(t) = e^Qt. SinceQ is a matrix, this is a matrix exponential, and can be approximated by the Taylor series expansion

If Q is diagonalizable, the exponential can be computed directly: let be a diagonalization of Q, with

where {λ_i} are the eigenvalues of Q, each repeated according to its multiplicity. Then

where

The time reversibility (or stationarity) constraint is ΠQ = 0 because the rows were defined to sum to zero, and the overall base frequencies must not systematically change from Π. This is equivalent to saying ΠP(t) = Π for all t.

GTR: Generalised time reversible

GTR is the most general neutral, independent, finite-sites, time-reversible model possible. It was first described in a general form by Simon Tavaré in 1986.^[7]

The GTR parameters consist of an equilibrium base frequency vector, Π = (π₁,π₂,π₃,π₄), giving the frequency at which each base occurs at each site, and the rate matrix

Therefore, GTR (for four character states, as is often the case in phylogenetics, namely when working with DNA or RNA data) requires 6 substitution rate parameters, as well as 4 equilibrium base frequency parameters. Since the 4 frequency parameters have to sum to 1, there are really only 3 free frequency parameters. This is usually further reduced to 8 parameters plus μ, the overall number of substitutions per unit time. When measuring time in substitutions (μ=1) only 8 free parameters remain.

In general, to compute the number of parameters, you count the number of entries above the diagonal in the matrix, i.e. for n trait values per site , and then add n-1 for the equilibrium base frequencies, and subtract 1 because μ is fixed. You get

For example, for an amino acid sequence (there are 20 "standard" amino acids that make up proteins), you would find there are 208 parameters. However, when studying coding regions of the genome, it is more common to work with acodon substitution model (a codon is three bases and codes for one amino acid in a protein). There are 4³ = 64codons, resulting in 2078 free parameters, but when the rates for transitions between codons which differ by more than one base are assumed to be zero, then there are only

parameters.

Mechanistic vs. empirical models

A main difference in evolutionary models is how many parameters are estimated every time for the data set under consideration and how many of them are estimated once on a large data set. Mechanistic models describe all substitution as a function of a number of parameters which are estimated for every data set analyzed, preferably usingmaximum likelihood. This has the advantage that the model can be adjusted to the particularities of a specific data set (e.g. different composition biases in DNA). Problems can arise when too many parameters are used, particularly if they can compensate for each other. Then it is often the case that the data set is too small to yield enough information to estimate all parameters accurately.

Empirical models are created by estimating many parameters (typically all entries of the rate matrix and the character frequencies, see the GTR model above) from a large data set. These parameters are then fixed and will be reused for every data set. This has the advantage that those parameters can be estimated more accurately. Normally, it is not possible to estimate all entries of the substitution matrix from the current data set only. On the downside, the estimated parameters might be too generic and do not fit a particular data set well enough.

With the large-scale genome sequencing still producing very large amounts of DNA and protein sequences, there is enough data available to create empirical models with any number of parameters. Because of the problems mentioned above, the two approaches are often combined, by estimating most of the parameters once on large-scale data, while a few remaining parameters are then adjusted to the data set under consideration. The following sections give an overview of the different approaches taken for DNA, protein or codon-based models.

Models of DNA substitution

See main article: Models of DNA evolution for more formal descriptions of the DNA models.

Models of DNA evolution were first proposed in 1969 by Jukes and Cantor ^[8], assuming equal transition rates as well as equal equilibrium frequencies for all bases. In 1980 Kimura^[9] introduced a model with two parameters: one for the transition and one for the transversion rate and in 1981, Felsenstein^[10] made a model in which the substitution rate corresponds to the equilibrium frequency of the target nucleotide. Hasegawa, Kishino and Yano (HKY)^[11] unified the two last models to a six parameter model. In the 1990s, models similar to HKY have been developed and refined by several researchers (e.g. ^[12] and ^[13]).

For DNA substitution models, mainly mechanistic models (as described above) are employed. The small number of parameters to estimate makes this feasible, but also DNA is often highly optimized for specific purposes (e.g. fastexpression or stability) depending on the organism and the type of gene, making it necessary to adjust the model to these circumstances.

Models of amino acid substitutions

For many analyses, particularly for longer evolutionary distances, the evolution is modeled on the amino acid level. Since not all DNA substitution also alter the encoded amino acid, information is lost when looking at amino acids instead of nucleotide bases. However, several advantages speak in favor of using the amino acid information: DNA is much more inclined to show compositional bias than amino acids, not all positions in the DNA evolve at the same speed (non-synonymous mutations are more likely to become fixed in the population than synonymous ones), but probably most important, because of those fast evolving positions and the limited alphabet size (only four possible states), the DNA suffers much more from back substitutions, making it difficult to accurately estimate longer distances.

Unlike the DNA models, amino acid models traditionally are empirical models. They were pioneered in the 1970s by Dayhoff and co-workers ^[14], by estimating replacement rates from protein alignments with at least 85% identity. This minimized the chances of observing multiple substitutions at a site. From the estimated rate matrix, a series of replacement probability matrices were derived, known under names such as PAM250. The Dayhoff model was used to assess the significance of homology search results, but also for phylogenetic analyses. The Dayhoff PAM matrices were based on relatively few alignments (since not more were available at that time), but in the 1990s, new matrices were estimated using almost the same methodology, but based on the large protein databases available then (^[15][16], the latter being known as "JTT" matrices).

References

• 1. ^ X Gu and W Li (1992): Higher rates of amino acid substitution in rodents than in man. Molecular Phylogenetics and Evolution, 1:211-214. DOI

  2. ^ WH Li, J Ellsworth, BH Krushkal, J Chang, and D Hewett-Emmett (1996): Rates of nucleotide substitution in primates and rodents and the generation-time effect hypothesis. Molecular Phylogenetics and Evolution, 5:182-187. DOI

  3. ^ AP Martin and SR Palumbi (1993): Body size, metabolic rate, generation time, and the molecular clock. Proceedings of the National Academy of Science, USA, 90:4087-4091. PNAS

  4. ^ Z Yang and R Nielsen (1998): Synonymous and nonsynonymous rate variation in nuclear genes of mammals. Journal of Molecular Evolution, 46:409-418. DOI

  5. ^ H Kishino, JL Thorne and WJ Bruno (2001): Performance of a divergence estimation method under a probabilistic model of rate evolution. Molecular Biology and Evolution 18: 352-361.PubMed

  6. ^ JL Thorne, H Kishino and IS Painter (1998): Estimating the rate of evolution of the rate of molecular evolution: Molecular Biology and Evolution 15: 1647-1657. MBE

  7. ^ Tavaré S. "Some Probabilistic and Statistical Problems in the Analysis of DNA Sequences". American Mathematical Society: Lectures on Mathematics in the Life Sciences 17: 57–86.

  8. ^ Jukes, TH and Cantor, CR. 1969. Evolution of protein molecules. Pp. 21-123 in H. N. Munro, ed. Mammalian protein metabolism. Academic Press, New York.

  9. ^ Kimura, M. 1980. A simple method for estimating evolutionary rate of base substitution through comparative studies of nucleotide sequences. Journal of Molecular Evolution 16:111-120. DOI

  10. ^ Felsenstein, J. 1981. Evolutionary trees from DNA sequences: a maximum likelihood approach. Journal of Molecular Evolution 17:368-376. DOI

  11. ^ Hasegawa, M., Kishino, H and Yano, T. 1985. Dating the human-ape splitting by a molecular clock of mitochondrial DNA. Journal of Molecular Evolution 22:160-174. DOI

  12. ^ Tamura, K. 1992. Estimation of the number of nucleotide substitutions when there are strong transition-transversion and G+C content biases. Molecular Biology and Evolution 9:678-687. MBE

  13. ^ Tamura, K., and M. Nei. 1993. Estimation of the number of nucleotide substitutions in the control region of mitochondrial DNA in humans and chimpanzees. Molecular Biology and Evolution 10:512-526. MBE

  14. ^ Dayhoff MO, Schwartz RM, Orcutt BC (1978). "A model for evolutionary change in proteins". Atlas of Protein Sequence and Structure5: 345–352.

  15. ^ Gonnet GH, Cohen MA, Benner SA (1992). "Exhaustive matching of the entire protein sequence database". Science 256: 1443–1445.doi:10.1126/science.1604319. PMID 1604319.

  16. ^ Jones DT, Taylor WR, Thornton JM (1992). "The rapid generation of mutation data matrices from protein sequences". Comput Applic Biosci 8: 275–282.

Categories: Bioinformatics | Stochastic processes | Computational phylogenetics | Statistical genetics