In biology, a substitution model, also called models of sequence evolution, are Markov models that describe changes over evolutionary time. These models describe evolutionary changes in macromolecules, such as DNA sequences or protein sequences, that can be represented as a sequence of symbols (e.g., A, C, G, and T in the case of DNA or the 20 "standard" proteinogenic amino acids in the case of proteins). Substitution models are 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 as input for distance methods such as neighbor joining). Substitution models are also central to phylogenetic invariants because they are necessary to predict site pattern frequencies given a tree topology. Substitution models are also necessary to simulate sequence data for a group of organisms related by a specific tree.
Theoretical foundations
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 Q i j {\displaystyle Q_{ij}} 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:
Q i i = − ∑ { j ∣ j ≠ i } Q i j , {\displaystyle Q_{ii}=-{\sum _{\lbrace j\mid j\neq i\rbrace }Q_{ij}}\,,}
The equilibrium row vector π must be annihilated by the rate matrix Q:
π Q = 0 . {\displaystyle \pi \,Q=0\,.}
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 ) {\displaystyle P(t)} . The entry in the ith column and the jth row, P i j ( t ) {\displaystyle 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 Pij(t) are such that Pij(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, lim t → ∞ P i j ( t ) = π j , {\displaystyle \lim _{t\rightarrow \infty }P_{ij}(t)=\pi _{j}\,,} 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 π P ( t ) = π {\displaystyle \pi P(t)=\pi } for all t. The transition matrix can be computed from the rate matrix via matrix exponentiation:
P ( t ) = e Q t = ∑ n = 0 ∞ Q n t n n ! , {\displaystyle P(t)=e^{Qt}=\sum _{n=0}^{\infty }Q^{n}{\frac {t^{n}}{n!}}\,,}
where Qn is the matrix Q multiplied by itself enough times to give its nth power. If Q is diagonalizable, the matrix exponential can be computed directly: let Q = U−1 Λ U be a diagonalization of Q, with
Λ = ( λ 1 … 0 ⋮ ⋱ ⋮ 0 … λ 4 ) , {\displaystyle \Lambda ={\begin{pmatrix}\lambda _{1}&\ldots &0\\\vdots &\ddots &\vdots \\0&\ldots &\lambda _{4}\end{pmatrix}}\,,}
where Λ is a diagonal matrix and where { λ i } {\displaystyle \lbrace \lambda _{i}\rbrace } are the eigenvalues of Q, each repeated according to its multiplicity. Then
P ( t ) = e Q t = e U − 1 ( Λ t ) U = U − 1 e Λ t U , {\displaystyle P(t)=e^{Qt}=e^{U^{-1}(\Lambda t)U}=U^{-1}e^{\Lambda t}\,U\,,}
where the diagonal matrix eΛt is given by
e Λ t = ( e λ 1 t … 0 ⋮ ⋱ ⋮ 0 … e λ 4 t ) . {\displaystyle e^{\Lambda t}={\begin{pmatrix}e^{\lambda _{1}t}&\ldots &0\\\vdots &\ddots &\vdots \\0&\ldots &e^{\lambda _{4}t}\end{pmatrix}}\,.}
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)
π i Q i j = π j Q j i {\displaystyle \pi _{i}Q_{ij}=\pi _{j}Q_{ji}}
or, equivalently, the detailed balance property,
π i P ( t ) i j = π j P ( t ) j i {\displaystyle \pi _{i}P(t)_{ij}=\pi _{j}P(t)_{ji}}
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.
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, π → = ( π 1 , π 2 , π 3 , π 4 ) {\displaystyle {\vec {\pi }}=(\pi _{1},\pi _{2},\pi _{3},\pi _{4})} , giving the frequency at which each base occurs at each site, and the rate matrix
Q = ( − ( x 1 + x 2 + x 3 ) x 1 x 2 x 3 π 1 x 1 π 2 − ( π 1 x 1 π 2 + x 4 + x 5 ) x 4 x 5 π 1 x 2 π 3 π 2 x 4 π 3 − ( π 1 x 2 π 3 + π 2 x 4 π 3 + x 6 ) x 6 π 1 x 3 π 4 π 2 x 5 π 4 π 3 x 6 π 4 − ( π 1 x 3 π 4 + π 2 x 5 π 4 + π 3 x 6 π 4 ) ) {\displaystyle Q={\begin{pmatrix}{-(x_{1}+x_{2}+x_{3})}&x_{1}&x_{2}&x_{3}\\{\pi _{1}x_{1} \over \pi _{2}}&{-({\pi _{1}x_{1} \over \pi _{2}}+x_{4}+x_{5})}&x_{4}&x_{5}\\{\pi _{1}x_{2} \over \pi _{3}}&{\pi _{2}x_{4} \over \pi _{3}}&{-({\pi _{1}x_{2} \over \pi _{3}}+{\pi _{2}x_{4} \over \pi _{3}}+x_{6})}&x_{6}\\{\pi _{1}x_{3} \over \pi _{4}}&{\pi _{2}x_{5} \over \pi _{4}}&{\pi _{3}x_{6} \over \pi _{4}}&{-({\pi _{1}x_{3} \over \pi _{4}}+{\pi _{2}x_{5} \over \pi _{4}}+{\pi _{3}x_{6} \over \pi _{4}})}\end{pmatrix}}}
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 μ {\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, you count the number of entries above the diagonal in the matrix, i.e. for n trait values per site n 2 − n 2 {\displaystyle {{n^{2}-n} \over 2}} , and then add n-1 for the equilibrium frequencies, and subtract 1 because μ {\displaystyle \mu } is fixed. You get
n 2 − n 2 + ( n − 1 ) − 1 = 1 2 n 2 + 1 2 n − 2. {\displaystyle {{n^{2}-n} \over 2}+(n-1)-1={1 \over 2}n^{2}+{1 \over 2}n-2.}
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 a codon substitution model (a codon is three bases and codes for one amino acid in a protein). There are 4 3 = 64 {\displaystyle 4^{3}=64} codons, resulting in 2078 free parameters. However, the rates for transitions between codons which differ by more than one base are often assumed to be zero, reducing the number of free parameters to only 20 × 19 × 3 2 + 63 − 1 = 632 {\displaystyle {{20\times 19\times 3} \over 2}+63-1=632} parameters. Another common practice is to reduce the number of codons by forbidding the stop (or nonsense) codons. This is a biologically reasonable assumption because including the stop codons would mean that one is calculating the probability of finding sense codon j {\displaystyle j} after time t {\displaystyle t} given that the ancestral codon is i {\displaystyle i} would involve the possibility of passing through a state with a premature stop codon. An alternative (and commonly used) way to write the instantaneous rate matrix ( Q {\displaystyle Q} matrix) for the nucleotide GTR model is:
Q = ( − ( a π C + b π G + c π T ) a π C b π G c π T a π A − ( a π A + d π G + e π T ) d π G e π T b π A d π C − ( b π A + d π C + f π T ) f π T c π A e π C f π G − ( c π A + e π C + f π G ) ) {\displaystyle Q={\begin{pmatrix}{-(a\pi _{C}+b\pi _{G}+c\pi _{T})}&a\pi _{C}&b\pi _{G}&c\pi _{T}\\a\pi _{A}&{-(a\pi _{A}+d\pi _{G}+e\pi _{T})}&d\pi _{G}&e\pi _{T}\\b\pi _{A}&d\pi _{C}&{-(b\pi _{A}+d\pi _{C}+f\pi _{T})}&f\pi _{T}\\c\pi _{A}&e\pi _{C}&f\pi _{G}&{-(c\pi _{A}+e\pi _{C}+f\pi _{G})}\end{pmatrix}}}
The Q {\displaystyle Q} matrix is normalized so − ∑ i = 1 4 π i Q i i = 1 {\displaystyle -\sum _{i=1}^{4}\pi _{i}Q_{ii}=1} . This notation is easier to understand than the notation originally used by Tavaré, because all model parameters correspond either to "exchangeability" parameters ( a {\displaystyle a} through f {\displaystyle f} , which can also be written using the notation r i j {\displaystyle r_{ij}} ) or to equilibrium nucleotide frequencies π → = ( π A , π C , π G , π T ) {\displaystyle {\vec {\pi }}=(\pi _{A},\pi _{C},\pi _{G},\pi _{T})} . Note that the nucleotides in the Q {\displaystyle Q} matrix have been written in alphabetical order. In other words, the transition probability matrix for the Q {\displaystyle Q} matrix above would be:
P ( t ) = e Q t = ( p A A ( t ) p A C ( t ) p A G ( t ) p A T ( t ) p C A ( t ) p C C ( t ) p C G ( t ) p C T ( t ) p G A ( t ) p G C ( t ) p G G ( t ) p G T ( t ) p T A ( t ) p T C ( t ) p T G ( t ) p T T ( t ) ) {\displaystyle P(t)=e^{Qt}={\begin{pmatrix}p_{\mathrm {AA} }(t)&p_{\mathrm {AC} }(t)&p_{\mathrm {AG} }(t)&p_{\mathrm {AT} }(t)\\p_{\mathrm {CA} }(t)&p_{\mathrm {CC} }(t)&p_{\mathrm {CG} }(t)&p_{\mathrm {CT} }(t)\\p_{\mathrm {GA} }(t)&p_{\mathrm {GC} }(t)&p_{\mathrm {GG} }(t)&p_{\mathrm {GT} }(t)\\p_{\mathrm {TA} }(t)&p_{\mathrm {TC} }(t)&p_{\mathrm {TG} }(t)&p_{\mathrm {TT} }(t)\end{pmatrix}}}
Some publications write the nucleotides in a different order (e.g., some authors choose to group two purines together and the two pyrimidines together; see also models of DNA evolution). These differences in notation make it important to be clear regarding the order of the states when writing the Q {\displaystyle Q} matrix. The value of this notation is that instantaneous rate of change from nucleotide i {\d