Substitution model

In biology, a substitution model, also called models of DNA sequence evolution, are Markov models that describe changes over evolutionary time. These models describe evolutionary changes in macromolecules (e.g., DNA sequences) represented as sequence of symbols (A, C, G, and T in the case of DNA). Substitution models is used to calculate the likelihood of phylogenetic trees using multiple sequence alignment data. Thus, substitution models are central to maximum likelihood estimation of phylogeny as well as Bayesian inference in phylogeny. Estimates of evolutionary distances (numbers of substitutions that have occurred since a pair of sequences diverged from a common ancestor) are typically calculated using substitution models (evolutionary distances are used input for distance methods such as neighbor joining). Substitution models are also central to phylogenetic invariants since they can be used predict the frequencies of site pattern frequencies given a tree topology. Substitution models are necessary to simulate sequence data for a group of organisms related by a specific tree. 
Phylogenetic tree topologies and other parameters.
Phylogenetic tree topologies are often the parameter of interest; thus, branch lengths and any other parameters describing the substitution process are often viewed as nuisance parameters. However, biologists are sometimes interested in the other aspects of the model. For example, branch lengths, especially when those branch lengths are combined with information from the fossil record and a model to estimate the timeframe for evolution. Other model parameters have been used to gain insights into various aspects of the process of evolution. The Ka/Ks ratio (also called ω in codon substitution models) is a parameter of interest in many studies. The Ka/Ks ratio can be used to examine the action of natural selection on protein-coding regions; it provides information about the relative rates of nucleotide substitutions that change amino acids (non-synonymous substitutions) to those that do not change the encoded amino acid (synonymous substitutions).
Application to sequence data.
Most of the work on substitution models has focused on DNA/RNA and protein sequence evolution. Models of DNA sequence evolution, where the alphabet corresponds to the four nucleotides (A, C, G, and T), are probably the easiest models to understand. DNA models can also be used to examine RNA virus evolution; this reflects the fact that RNA also has a four nucleotide alphabet (A, C, G, and U). However, substitution models can be used for alphabets of any size; the alphabet is the 20 proteinogenic amino acids for proteins and the sense codons (i.e., the 61 codons that encode amino acids in the standard genetic code) for aligned protein-coding gene sequences. In fact, substitution models can be developed for any biological characters that can be encoded using a specific alphabet (e.g., amino acid sequences combined with information about the conformation of those amino acids in three-dimensional protein structures).
The majority of substitution models used for evolutionary research assume independence among sites (i.e., the probability of observing any specific site pattern is identical regardless of where the site pattern is in the sequence alignment). This simplifies likelihood calculations because it is only necessary to calculate the probability of all site patterns that appear in the alignment then use those values to calculate the overall likelihood of the alignment (e.g., the probability of three "GGGG" site patterns given some model of DNA sequence evolution is simply the probability of a single "GGGG" site pattern raised to the third power). This means that substitution models can be viewed as implying a specific multinomial distribution for site pattern frequencies. If we consider a multiple sequence alignment of with four DNA sequences there are 256 possible site pattens so there are 255 degrees of freedom for the site pattern frequencies. However, it is possible to specify. the expected site pattern frequencies using five degrees of freedom if using the Jukes-Cantor model of DNA evolution, which is a simple substitution model that allows one to calculate the expected site pattern frequencies only the tree topology and the branch lengths (given four taxa an unrooted bifurcating tree has five branch lengths).
Substitution models also make it possible to simulate sequence data using Monte Carlo methods. Simulated multiple sequence alignments can be used to assess the performance of phylogenetic methods and generate the null distribution for certain statistical tests in the fields of molecular evolution and molecular phylogenetics. Examples of these tests include tests of model fit and the "SOWH test" that can be used to examine tree topologies.
Application to morphological data.
The fact that substitution models can be used to analyze any biological alphabet has made it possible to develop models of evolution for phenotypic datasets (e.g., morphological and behavioural traits). Typically, "0" is. used to indicate the absence of a trait and "1" is used to indicate the presence of a trait, although it is also possible to score characters using multiple states. Using this framework, we might encode a set of phenotypes as binary strings (this could be generalized to "k"-state strings for characters with more than two states) before analyses using an appropriate mode. This can be illustrated using a "toy" example: we can use a binary alphabet to score the following phenotypic traits "has feathers", "lays eggs", "has fur", "is warm-blooded", and "capable of powered flight". In this toy example hummingbirds would have sequence 11011 (most other birds would have the same string), ostriches would have the sequence 11010, cattle (and most other land mammals) would have 00110, and bats would have 00111. The likelihood of a phylogenetic tree can then be calculated using those binary sequences and an appropriate substitution model. The existence of these morphological models make it possible to analyze data matrices with fossil taxa, either using the morphological data alone or a combination of morphological and molecular data (with the latter scored as missing data for the fossil taxa).
There is an obvious similarity between use of molecular or phenotypic data in the field of cladistics and analyses of morphological characters using a substitution model. However, there has been a vociferous debate in the systematics community regarding the question of whether or not cladistic analyses should be viewed as "model-free". The field of cladistics (defined in the strictest sense) favor the use of the maximum parsimony criterion for phylogenetic inference. Many cladists reject the position that maximum parsimony is based on a substitution model and (in many cases) they justify the use of parsimony using the philosophy of Karl Popper. However, the existence of "parsimony-equivalent" models (i.e., substitution models that yield the maximum parsimony tree when used for analyses) makes it possible to view parsimony as a substitution model.
The molecular clock and the units of time.
Typically, a branch length of a phylogenetic tree is expressed as the expected number of substitutions per site; if the evolutionary model indicates that each site within an ancestral sequence will typically experience "x" substitutions by the time it evolves to a particular descendant's sequence then the ancestor and descendant are considered to be separated by branch length "x".
Sometimes a branch length is measured in terms of geological years. For example, a fossil record may make it possible to determine the number of years between an ancestral species and a descendant species. Because some species evolve at faster rates than others, these two measures of branch length are not always in direct proportion. The expected number of substitutions per site per year is often indicated with the Greek letter mu (μ).
A model is said to have a strict molecular clock if the expected number of substitutions per year μ is constant regardless of which species' evolution is being examined. An important implication of a strict molecular clock is that the number of expected substitutions between an ancestral species and any of its present-day descendants must be independent of which descendant species is examined.
Note that the assumption of a strict 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. This could be due to their shorter generation time, higher metabolic rate, increased population structuring, increased rate of speciation, or smaller body size. When studying ancient 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.
Models that can take into account variability of the rate of the molecular clock between different evolutionary lineages in the phylogeny are called “relaxed” in opposition to “strict”. In such models the rate can be assumed to be correlated or not between ancestors and descendants and rate variation among lineages can be drawn from many distributions but usually exponential and lognormal distributions are applied. There is a special case, called “local molecular clock” when a phylogeny is divided into at least two partitions (sets of lineages) and a strict molecular clock is applied in each, but with different rates.
Time-reversible and stationary models.
Many useful substitution models are time-reversible; in terms of the mathematics, the model does not care which sequence is the ancestor and which is the descendant so long as all other parameters (such as the number of substitutions per site that is expected between the two sequences) are held constant.
When an analysis of real biological data is performed, there is generally no access to the sequences of ancestral species, only to the present-day species. However, when a model is time-reversible, which species was the ancestral species is irrelevant. Instead, the phylogenetic tree can be rooted using any of the species, re-rooted later based on new knowledge, or left unrooted. This is because there is no 'special' species, all species will eventually derive from one another with the same probability.
A model is time reversible if and only if it satisfies the property (the notation is explained below)
or, equivalently, the detailed balance property,
for every "i", "j", and "t".
Time-reversibility should not be confused with stationarity. A model is stationary if "Q" does not change with time. The analysis below assumes a stationary model.
The mathematics of substitution models.
Stationary, neutral, independent, finite sites models (assuming a constant rate of evolution) have two parameters, "π", an equilibrium vector of base (or character) frequencies and a rate matrix, "Q", which describes the rate at which bases of one type change into bases of another type; element formula_3 for "i" ≠ "j" is the rate at which base "i" goes to base "j". The diagonals of the "Q" matrix are chosen so that the rows sum to zero:
The equilibrium row vector "π" must be annihilated by the rate matrix "Q":
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 formula_6. The entry in the "i"th column and the "j"th row, formula_7, 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 "P""ij"(0) = δ"ij", 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, formula_8 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. Furthermore, it follows that formula_9 for all "t".
The transition matrix can be computed from the rate matrix via matrix exponentiation:
where "Q""n" is the matrix "Q" multiplied by itself enough times to give its "n"th power.
If "Q" is diagonalizable, the matrix exponential can be computed directly: let "Q" = "U"−1 Λ "U" be a diagonalization of "Q", with
where Λ is a diagonal matrix and where formula_12 are the eigenvalues of "Q", each repeated according to its multiplicity. Then
where the diagonal matrix "e"Λt is given by
Generalised time reversible.
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. The GTR model is often called the general time reversible model in publications; it has also been called the REV model.
The GTR parameters for nucleotides consist of an equilibrium base frequency vector, formula_15, giving the frequency at which each base occurs at each site, and the rate matrix
Because the model must be time reversible and must approach the equilibrium nucleotide (base) frequencies at long times, each rate below the diagonal equals the reciprocal rate above the diagonal multiplied by the equilibrium ratio of the two bases. As such, the nucleotide GTR requires 6 substitution rate parameters and 4 equilibrium base frequency parameters. Since the 4 frequency parameters must sum to 1, there are only 3 free frequency parameters. The total of 9 free parameters is often further reduced to 8 parameters plus formula_17, the overall number of substitutions per unit time. When measuring time in substitutions (formula_17=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 <math>