This article includes a list ofgeneral references, butit lacks sufficient correspondinginline citations. Please help toimprove this article byintroducing more precise citations.(November 2010) (Learn how and when to remove this message)

A number of differentMarkovmodels of DNA sequence evolution have been proposed.^[1] Thesesubstitution models differ in terms of the parameters used to describe the rates at which onenucleotide replaces another during evolution. These models are frequently used inmolecular phylogenetic analyses. In particular, they are used during the calculation of likelihood of a tree (inBayesian andmaximum likelihood approaches to tree estimation) and they are used to estimate the evolutionary distance between sequences from the observed differences between the sequences.

These models are phenomenological descriptions of the evolution of DNA as a string of four discrete states. These Markov models do not explicitly depict the mechanism of mutation nor the action of natural selection. Rather they describe the relative rates of different changes. For example, mutational biases andpurifying selection favoring conservative changes are probably both responsible for the relatively high rate oftransitions compared totransversions in evolving sequences. However, the Kimura (K80) model described below only attempts to capture the effect of both forces in a parameter that reflects the relative rate of transitions to transversions.

Evolutionary analyses of sequences are conducted on a wide variety of time scales. Thus, it is convenient to express these models in terms of the instantaneous rates of change between different states (theQ matrices below). If we are given a starting (ancestral) state at one position, the model'sQ matrix and a branch length expressing the expected number of changes to have occurred since the ancestor, then we can derive the probability of the descendant sequence having each of the four states. The mathematical details of this transformation from rate-matrix to probability matrix are described inthe mathematics of substitution models section of thesubstitution model page. By expressing models in terms of the instantaneous rates of change we can avoid estimating a large numbers of parameters for each branch on a phylogenetic tree (or each comparison if the analysis involves many pairwise sequence comparisons).

The models described on this page describe the evolution of a single site within a set of sequences. They are often used for analyzing the evolution of an entirelocus by making the simplifying assumption that different sites evolveindependently and are identically distributed. This assumption may be justifiable if the sites can be assumed to be evolvingneutrally. If the primary effect of natural selection on the evolution of the sequences is to constrain some sites, then models of among-site rate-heterogeneity can be used. This approach allows one to estimate only one matrix of relative rates of substitution, and another set of parameters describing the variance in the total rate of substitution across sites.

Continuous-timeMarkov chains have the usual transition matrices which are, in addition, parameterized by time,${\displaystyle t}$. Specifically, if${\displaystyle E_1,E_2,E_3,E_4}$ are the states, then the transition matrix

${\displaystyle P(t)=(P_{ij}(t))}$ where each individual entry,${\displaystyle P_{ij}(t)}$ refers to the probability that state${\displaystyle E_i}$ will change to state${\displaystyle E_j}$ in time${\displaystyle t}$.

Example: We would like to model the substitution process in DNA sequences (i.e.Jukes–Cantor, Kimura,etc.) in a continuous-time fashion. The corresponding transition matrices will look like:

where the top-left and bottom-right 2 × 2 blocks correspond totransition probabilities and the top-right and bottom-left 2 × 2 blocks corresponds totransversion probabilities.

Assumption: If at some time${\displaystyle t_0}$, the Markov chain is in state${\displaystyle E_i}$, then the probability that at time${\displaystyle t_0+t}$, it will be in state${\displaystyle E_j}$ depends only upon${\displaystyle i}$,${\displaystyle j}$ and${\displaystyle t}$. This then allows us to write that probability as${\displaystyle p_{ij}(t)}$.

Theorem: Continuous-time transition matrices satisfy:

${\displaystyle P(t+\tau )=P(t)P(\tau )}$

Note: There is here a possible confusion between two meanings of the wordtransition. (i) In the context ofMarkov chains, transition is the general term for the change between two states. (ii) In the context ofnucleotide changes in DNA sequences, transition is a specific term for the exchange between either the two purines (A ↔ G) or the two pyrimidines (C ↔ T) (for additional details, see the article abouttransitions in genetics). By contrast, an exchange between one purine and one pyrimidine is called atransversion.

Consider a DNA sequence of fixed lengthm evolving in time by base replacement. Assume that the processes followed by them sites are Markovian independent, identically distributed and that the process is constant over time. For a particular site, let

${\displaystyle \mathcal{E}=\{A,G,C,T\}}$

be the set of possible states for the site, and

${\displaystyle \mathbf{p}(t)=(p_A(t),p_G(t),p_C(t),p_T(t))}$

their respective probabilities at time${\displaystyle t}$. For two distinct${\displaystyle x,y\in \mathcal{E}}$, let${\displaystyle {\mu }_{xy}}$ be the transition rate from state${\displaystyle x}$ to state${\displaystyle y}$. Similarly, for any${\displaystyle x}$, let the total rate of change from${\displaystyle x}$ be

${\displaystyle {\mu }_x=\sum _{y\ne x}{\mu }_{xy}.}$

The changes in the probability distribution${\displaystyle p_A(t)}$ for small increments of time${\displaystyle \mathrm{\Delta }t}$ are given by

${\displaystyle p_A(t+\mathrm{\Delta }t)=p_A(t)-p_A(t){\mu }_A\mathrm{\Delta }t+\sum _{x\ne A}{p}_x(t){\mu }_{xA}\mathrm{\Delta }t.}$

In other words, (in frequentist language), the frequency of${\displaystyle A}$'s at time${\displaystyle t+\mathrm{\Delta }t}$ is equal to the frequency at time${\displaystyle t}$ minus the frequency of thelost${\displaystyle A}$'s plus the frequency of thenewly created${\displaystyle A}$'s.

Similarly for the probabilities${\displaystyle p_G(t)}$,${\displaystyle p_C(t)}$ and${\displaystyle p_T(t)}$. These equations can be written compactly as

${\displaystyle \mathbf{p}(t+\mathrm{\Delta }t)=\mathbf{p}(t)+\mathbf{p}(t)Q\mathrm{\Delta }t,}$

where

is known as therate matrix. Note that, by definition, the sum of the entries in each row of${\displaystyle Q}$ is equal to zero. It follows that

${\displaystyle {\mathbf{p}}^{\prime }(t)=\mathbf{p}(t)Q.}$

For astationary process, where${\displaystyle Q}$ does not depend on timet, this differential equation can be solved. First,

${\displaystyle P(t)=\exp (tQ),}$

where${\displaystyle \exp (tQ)}$ denotes theexponential of the matrix${\displaystyle tQ}$. As a result,

${\displaystyle \mathbf{p}(t)=\mathbf{p}(0)P(t)=\mathbf{p}(0)\exp (tQ).}$

If the Markov chain isirreducible,i.e. if it is always possible to go from a state${\displaystyle x}$ to a state${\displaystyle y}$ (possibly in several steps), then it is alsoergodic. As a result, it has a uniquestationary distribution${\displaystyle \mathit{\pi }=\{{\pi }_x,x\in \mathcal{E}\}}$, where${\displaystyle {\pi }_x}$ corresponds to the proportion of time spent in state${\displaystyle x}$ after the Markov chain has run for an infinite amount of time. In DNA evolution, under the assumption of a common process for each site, the stationary frequencies${\displaystyle {\pi }_A,{\pi }_G,{\pi }_C,{\pi }_T}$ correspond to equilibrium base compositions. Indeed, note that since the stationary distribution${\displaystyle \mathit{\pi }}$ satisfies${\displaystyle \mathit{\pi }Q=0}$, we see that when the current distribution${\displaystyle \mathbf{p}(t)}$ is the stationary distribution${\displaystyle \mathit{\pi }}$ we have

${\displaystyle {\mathbf{p}}^{\prime }(t)=\mathbf{p}(t)Q=\mathit{\pi }Q=0.}$

In other words, the frequencies of${\displaystyle p_A(t),p_G(t),p_C(t),p_T(t)}$ do not change.

Definition: A stationary Markov process istime reversible if (in the steady state) the amount of change from state${\displaystyle x}$ to${\displaystyle y}$ is equal to the amount of change from${\displaystyle y}$ to${\displaystyle x}$, (although the two states may occur with different frequencies). This means that:

${\displaystyle {\pi }_x{\mu }_{xy}={\pi }_y{\mu }_{yx}}$

Not all stationary processes are reversible, however, most commonly used DNA evolution models assume time reversibility, which is considered to be a reasonable assumption.

Under the time reversibility assumption, let${\displaystyle s_{xy}={\mu }_{xy}/{\pi }_y}$, then it is easy to see that:

${\displaystyle s_{xy}=s_{yx}}$

Definition The symmetric term${\displaystyle s_{xy}}$ is called theexchangeability between states${\displaystyle x}$ and${\displaystyle y}$. In other words,${\displaystyle s_{xy}}$ is the fraction of the frequency of state${\displaystyle x}$ that is the result of transitions from state${\displaystyle y}$ to state${\displaystyle x}$.

Corollary The 12 off-diagonal entries of the rate matrix,${\displaystyle Q}$ (note the off-diagonal entries determine the diagonal entries, since the rows of${\displaystyle Q}$ sum to zero) can be completely determined by 9 numbers; these are: 6 exchangeability terms and 3 stationary frequencies${\displaystyle {\pi }_x}$, (since the stationary frequencies sum to 1).

By comparing extant sequences, one can determine the amount of sequence divergence. This raw measurement of divergence provides information about the number of changes that have occurred along the path separating the sequences. The simple count of differences (theHamming distance) between sequences will often underestimate the number of substitution because of multiple hits (seehomoplasy). Trying to estimate the exact number of changes that have occurred is difficult, and usually not necessary. Instead, branch lengths (and path lengths) in phylogenetic analyses are usually expressed in the expected number of changes per site. The path length is the product of the duration of the path in time and the mean rate of substitutions. While their product can be estimated, the rate and time are not identifiable from sequence divergence.

The descriptions of rate matrices on this page accurately reflect the relative magnitude of different substitutions, but these rate matrices arenot scaled such that a branch length of 1 yields one expected change. This scaling can be accomplished by multiplying every element of the matrix by the same factor, or simply by scaling the branch lengths. If we use the β to denote the scaling factor, and ν to denote the branch length measured in the expected number of substitutions per site then βν is used in the transition probability formulae below in place of μt. Note that ν is a parameter to be estimated from data, and is referred to as the branch length, while β is simply a number that can be calculated from the rate matrix (it is not a separate free parameter).

The value of β can be found by forcing the expected rate of flux of states to 1. The diagonal entries of the rate-matrix (theQ matrix) represent -1 times the rate of leaving each state. Fortime-reversible models, we know the equilibrium state frequencies (these are simply the π_i parameter value for statei). Thus we can find the expected rate of change by calculating the sum of flux out of each state weighted by the proportion of sites that are expected to be in that class. Setting β to be the reciprocal of this sum will guarantee that scaled process has an expected flux of 1:

${\displaystyle \beta =1/\left(-\sum _i{\pi }_i{\mu }_{ii}\right)}$

For example, in the Jukes–Cantor, the scaling factor would be4/(3μ) because the rate of leaving each state is3μ/4.

JC69, theJukes andCantor 1969 model,^[2] is the simplestsubstitution model. There are several assumptions. It assumes equal base frequencies${\displaystyle \left({\pi }_A={\pi }_G={\pi }_C={\pi }_T=\frac{1}{4}\right)}$ and equalmutation rates. The only parameter of this model is therefore${\displaystyle \mu }$, the overall substitution rate. As previously mentioned, this variable becomes a constant when we normalize the mean-rate to 1.

When branch length,${\displaystyle \nu }$, is measured in the expected number of changes per site then:

It is worth noticing that${\displaystyle \nu =\frac{3}{4}t\mu =(\frac{\mu }{4}+\frac{\mu }{4}+\frac{\mu }{4})t}$ what stands for sum of any column (or row) of matrix${\displaystyle Q}$ multiplied by time and thus means expected number of substitutions in time${\displaystyle t}$ (branch duration) for each particular site (per site) when the rate of substitution equals${\displaystyle \mu }$.

Given the proportion${\displaystyle p}$ of sites that differ between the two sequences the Jukes–Cantor estimate of the evolutionary distance (in terms of the expected number of changes) between two sequences is given by

${\displaystyle \hat{d}=-\frac{3}{4}\ln (1-\frac{4}{3}p)=\hat{\nu }}$

The${\displaystyle p}$ in this formula is frequently referred to as the${\displaystyle p}$-distance. It is asufficient statistic for calculating the Jukes–Cantor distance correction, but is not sufficient for the calculation of the evolutionary distance under the more complex models that follow (also note that${\displaystyle p}$ used in subsequent formulae is not identical to the "${\displaystyle p}$-distance").

K80, theKimura 1980 model,^[3] often referred to asKimura's two parameter model (or theK2P model), distinguishes betweentransitions (${\displaystyle A\leftrightarrow G}$, i.e. from purine to purine, or${\displaystyle C\leftrightarrow T}$, i.e. from pyrimidine to pyrimidine) andtransversions (from purine to pyrimidine or vice versa). In Kimura's original description of the model the α and β were used to denote the rates of these types of substitutions, but it is now more common to set the rate of transversions to 1 and use κ to denote the transition/transversion rate ratio (as is done below). The K80 model assumes that all of the bases are equally frequent (${\displaystyle {\pi }_A={\pi }_G={\pi }_C={\pi }_T=\frac{1}{4}}$).

Rate matrix with columns corresponding to${\displaystyle A}$,${\displaystyle G}$,${\displaystyle C}$, and${\displaystyle T}$, respectively.

The Kimura two-parameter distance is given by:

${\displaystyle K=-\frac{1}{2}\ln ((1-2p-q)\sqrt{1-2q})}$

wherep is the proportion of sites that show transitional differences andq is the proportion of sites that show transversional differences.

K81, theKimura 1981 model,^[4] often calledKimura's three parameter model (K3P model) or the Kimura three substitution type (K3ST) model, has distinct rates fortransitions and two distinct types oftransversions. The twotransversion types are those that conserve the weak/strong properties of the nucleotides (i.e.,${\displaystyle A\leftrightarrow T}$ and${\displaystyle C\leftrightarrow G}$, denoted by symbol${\displaystyle \gamma }$^[4]) and those that conserve the amino/keto properties of the nucleotides (i.e.,${\displaystyle A\leftrightarrow C}$ and${\displaystyle G\leftrightarrow T}$, denoted by symbol${\displaystyle \beta }$^[4]). The K81 model assumes that all equilibrium base frequencies are equal (i.e.,${\displaystyle {\pi }_A={\pi }_G={\pi }_C={\pi }_T=0.25}$).

Rate matrix with columns corresponding to${\displaystyle A}$,${\displaystyle G}$,${\displaystyle C}$, and${\displaystyle T}$, respectively.

The K81 model is used much less often than the K80 (K2P) model for distance estimation and it is seldom the best-fitting model in maximum likelihood phylogenetics. Despite these facts, the K81 model has continued to be studied in the context of mathematical phylogenetics.^[5][6][7] One important property is the ability to perform aHadamard transform assuming the site patterns were generated on a tree with nucleotides evolving under the K81 model.^[8][9][10]

When used in the context of phylogenetics the Hadamard transform provides an elegant and fully invertible means to calculate expected site pattern frequencies given a set of branch lengths (or vice versa). Unlike many maximum likelihood calculations, the relative values for${\displaystyle \alpha }$,${\displaystyle \beta }$, and${\displaystyle \gamma }$ can vary across branches and the Hadamard transform can even provide evidence that the data do not fit a tree. The Hadamard transform can also be combined with a wide variety of methods to accommodate among-sites rate heterogeneity,^[11] using continuous distributions rather than the discrete approximations typically used in maximum likelihood phylogenetics^[12] (although one must sacrifice the invertibility of the Hadamard transform to use certain among-sites rate heterogeneity distributions^[11]).

F81, theFelsenstein's 1981 model,^[13] is an extension of the JC69 model in which base frequencies are allowed to vary from 0.25 (${\displaystyle {\pi }_A\ne {\pi }_G\ne {\pi }_C\ne {\pi }_T}$)

Rate matrix:

When branch length, ν, is measured in the expected number of changes per site then:

${\displaystyle \beta =1/(1-{\pi }_A^2-{\pi }_C^2-{\pi }_G^2-{\pi }_T^2)}$

HKY85, the Hasegawa, Kishino and Yano 1985 model,^[14] can be thought of as combining the extensions made in the Kimura80 and Felsenstein81 models. Namely, it distinguishes between the rate oftransitions andtransversions (using the κ parameter), and it allows unequal base frequencies (${\displaystyle {\pi }_A\ne {\pi }_G\ne {\pi }_C\ne {\pi }_T}$). [ Felsenstein described a similar (but not equivalent) model in 1984 using a different parameterization;^[15] that latter model is referred to as the F84 model.^[16] ]

Rate matrix

If we express the branch length,ν in terms of the expected number of changes per site then:

${\displaystyle \beta =\frac{1}{2({\pi }_A+{\pi }_G)({\pi }_C+{\pi }_T)+2\kappa [({\pi }_A{\pi }_G)+({\pi }_C{\pi }_T)]}}$

${\displaystyle P_{AA}(\nu ,\kappa ,\pi )=\left[{\pi }_A\left({\pi }_A+{\pi }_G+({\pi }_C+{\pi }_T)e^{-\beta \nu }\right)+{\pi }_G{e}^{-(1+({\pi }_A+{\pi }_G)(\kappa -1.0))\beta \nu }\right]/({\pi }_A+{\pi }_G)}$

${\displaystyle P_{AC}(\nu ,\kappa ,\pi )={\pi }_C\left(1.0-e^{-\beta \nu }\right)}$

${\displaystyle P_{AG}(\nu ,\kappa ,\pi )=\left[{\pi }_G\left({\pi }_A+{\pi }_G+({\pi }_C+{\pi }_T)e^{-\beta \nu }\right)-{\pi }_G{e}^{-(1+({\pi }_A+{\pi }_G)(\kappa -1.0))\beta \nu }\right]/\left({\pi }_A+{\pi }_G\right)}$

${\displaystyle P_{AT}(\nu ,\kappa ,\pi )={\pi }_T\left(1.0-e^{-\beta \nu }\right)}$

and formula for the other combinations of states can be obtained by substituting in the appropriate base frequencies.

T92, the Tamura 1992 model,^[17] is a mathematical method developed to estimate the number of nucleotide substitutions per site between two DNA sequences, by extendingKimura's (1980) two-parameter method to the case where aG+C content bias exists. This method will be useful when there are strong transition-transversion and G+C-content biases, as in the case ofDrosophila mitochondrial DNA.^[17]

T92 involves a single, compound base frequency parameter${\displaystyle \theta \in (0,1)}$ (also noted${\displaystyle {\pi }_{GC}}$)${\displaystyle ={\pi }_G+{\pi }_C=1-({\pi }_A+{\pi }_T)}$

As T92 echoes theChargaff's second parity rule — pairing nucleotides do have the same frequency on a single DNA strand, G and C on the one hand, and A and T on the other hand — it follows that the four base frequences can be expressed as a function of${\displaystyle {\pi }_{GC}}$

${\displaystyle {\pi }_G={\pi }_C=\frac{{\pi }_{GC}}{2}}$ and${\displaystyle {\pi }_A={\pi }_T=\frac{(1-{\pi }_{GC})}{2}}$

Rate matrix

The evolutionary distance between two DNA sequences according to this model is given by

${\displaystyle d=-h\ln (1-\frac{p}{h}-q)-\frac{1}{2}(1-h)\ln (1-2q)}$

where${\displaystyle h=2\theta (1-\theta )}$ and${\displaystyle \theta }$ is the G+C content (${\displaystyle {\pi }_{GC}={\pi }_G+{\pi }_C}$).

TN93, the Tamura andNei 1993 model,^[18] distinguishes between the two different types oftransition; i.e. (${\displaystyle A\leftrightarrow G}$) is allowed to have a different rate to (${\displaystyle C\leftrightarrow T}$).Transversions are all assumed to occur at the same rate, but that rate is allowed to be different from both of the rates for transitions.

TN93 also allows unequal base frequencies (${\displaystyle {\pi }_A\ne {\pi }_G\ne {\pi }_C\ne {\pi }_T}$).

Rate matrix

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

GTR parameters consist of an equilibrium base frequency vector,${\displaystyle \mathrm{\Pi }=({\pi }_A,{\pi }_G,{\pi }_C,{\pi }_T)}$, giving the frequency at which each base occurs at each site, and the rate matrix

Where

${\displaystyle \begin{array}{r}\alpha =r(A\to G)=r(G\to A)\\ \beta =r(A\to C)=r(C\to A)\\ \gamma =r(A\to T)=r(T\to A)\\ \delta =r(G\to C)=r(C\to G)\\ ϵ=r(G\to T)=r(T\to G)\\ \eta =r(C\to T)=r(T\to C)\end{array}}$

are the transition rate parameters.

Therefore, GTR (for four characters, as is often the case in phylogenetics) requires 6 substitution rate parameters, as well as 4 equilibrium base frequency parameters. However, this is usually eliminated down to 9 parameters plus${\displaystyle \mu }$, the overall number of substitutions per unit time. When measuring time in substitutions (${\displaystyle \mu }$=1) only 8 free parameters remain.

In general, to compute the number of parameters, one must count the number of entries above the diagonal in the matrix, i.e. for n trait values per site${\displaystyle \frac{n^2-n}{2}}$, and then addn for the equilibrium base frequencies, and subtract 1 because${\displaystyle \mu }$ is fixed. One gets

${\displaystyle \frac{n^2-n}{2}+n-1=\frac{1}{2}{n}^2+\frac{1}{2}n-1.}$

For example, for an amino acid sequence (there are 20 "standard" amino acids that make upproteins), one would find there are 209 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${\displaystyle 4^3=64}$ codons, but the rates for transitions between codons which differ by more than one base is assumed to be zero. Hence, there are${\displaystyle \frac{20\times 19\times 3}{2}+64-1=633}$ parameters.

An alternative way to analyze DNA sequence data is to recode the nucleotides as purines (R) and pyrimidines (Y);^[20][21] this practice is often called RY-coding.^[22] Insertions and deletions in multiple sequence alignments can also be encoded as binary data^[23] and analyzed in using a two-state model.^[24][25]

The simplest two-state model of sequence evolution is called the Cavender-Farris model or the Cavender-Farris-Neyman (CFN) model; the name of this model reflects the fact that it was described independently in several different publications.^[26][27][28] The CFN model is identical to the Jukes-Cantor model adapted to two states and it has even been implemented as the "JC2" model in the popularIQ-TREE software package (using this model in IQ-TREE requires coding the data as 0 and 1 rather than R and Y; the popularPAUP* software package can interpret a data matrix comprising only R and Y as data to be analyzed using the CFN model). It is also straightforward to analyze binary data using the phylogeneticHadamard transform.^[29] The alternative two-state model allows the equilibrium frequency parameters of R and Y (or 0 and 1) to take on values other than 0.5 by adding single free parameter; this model is variously called CFu^[20] or GTR2 (in IQ-TREE).

Other recoding methods are WS (weak-strong) and MK (amino-keto).

Lie Markov models are, from a mathematical point of view, Markov models that form aLie algebra.^[30] For the mathematician, this makes them closed undermatrix multiplication. From a phylogeneticist's point of view, these models have the benefit of being able to add or remove taxa without affecting the site patterns that the model can generate over the remaining taxa. There is also a natural hierarchy of models based on how many parameters can be changed. Some existing models such as JC and F81 are already Lie Markov models, while GTR is not.^[31] Lie Markov models (with RY, WS, or MK) are available in IQ-TREE.^[32]

• Molecular evolution
• Molecular clock
• UPGMA

• DAWG: DNA Assembly With Gaps — free software for simulating sequence evolution

• Bioinformatics
• Phylogenetics
• Computational phylogenetics
• Markov models

• Articles with short description
• Short description is different from Wikidata
• Articles lacking in-text citations from November 2010
• All articles lacking in-text citations