Notation

Recall that a rooted binary phylogenetic X-tree is a tree $\mathcal{T}=\left(V\right(\mathcal{T}),E(\mathcal{T}\left)\right)$ with the following properties: There is one vertex $\rho \in V\left(\mathcal{T}\right)$ with indegree 0 and outdegree 1, which is called the root of $\mathcal{T}$. All edges $e\in E\left(\mathcal{T}\right)$ are directed away from ρ, and all vertices $v\in V\left(\mathcal{T}\right)\setminus \left\{\rho \right\}$ have indegree 1 and outdegree 0 or 2. Vertices with outdegree 0 are usually referred to as leaves of $\mathcal{T}$. Remember that for an X-tree, there are exactly |X|=n leaves, which is why there is a bijection between the set of leaves of $\mathcal{T}$ and the taxon set X. Thus, when there is no ambiguity, we use the terms leaf and taxon synonymously. Moreover, we often just write “phylogenetic tree” or “tree” when referring to a rooted binary phylogenetic tree.

Furthermore, recall that a character f is a function $f:X\to \mathcal{C}$ for some set $\mathcal{C}:=\{c_1,c_2,c_3,\dots ,c_r\}$ of r character states ($r\in \mathbb{N}$). We denote by ${\mathcal{C}}^n$ the set of all r ⁿ possible characters on $\mathcal{C}$ and n taxa. For instance, for the four-state DNA alphabet, ${\mathcal{C}}_{\mathrm{DNA}}=\{A,G,C,T\}$ and the set ${\mathcal{C}}_{\mathrm{DNA}}^n$ consists of all 4 ⁿ possible characters.

An extension of f to $V\left(\mathcal{T}\right)$ is a map $g:V\left(\mathcal{T}\right)\to \mathcal{C}$ such that g(i)=f(i) for all i in X. For such an extension g of f, we denote by $l_{\mathcal{T}}\left(g\right)$ the number of edges e={u v} in $\mathcal{T}$ on which a substitution occurs, i.e. where g(u)≠g(v). The parsimony score of f on $\mathcal{T}$, denoted by $l_{\mathcal{T}}\left(f\right)$, is obtained by minimizing $l_{\mathcal{T}}\left(g\right)$ over all possible extensions g. Given a tree $\mathcal{T}$ and a character f on the same taxon set, one can easily calculate the parsimony score of f on $\mathcal{T}$ with the famous Fitch algorithm [7]. Moreover, when a character state changes along one edge of the tree, we refer to this state change as substitution or mutation. As for our purposes only so-called manifest mutations are relevant, i.e. those mutations that can be observed and are not reversed, we do not distinguish between mutations and substitutions, which is why we use these terms synonymously.

Construction of the OSM matrix

We now introduce the OSM framework in a stepwise fashion. The aim of the OSM approach is to determine the effects a single mutation occurring on a rooted tree $\mathcal{T}$ has on a character evolving on that tree.

The first task of this approach is to formalize the term mutation and its effects on a single character state in $\mathcal{C}$. A mutation is an operation $\sigma :\mathcal{C}\to \mathcal{C}$ which is bijective, i.e. it satisfies the following condition:

C1. For all $c_i\in \mathcal{C}$ there is a $c_j\in \mathcal{C}$ such that σ(c _i )=c _j , and if σ(c _i )=σ(c _j ), then c _i =c _j .

This guarantees that a mutation affects a character state in a unique fashion. It is well-known that any bijective function on a finite discrete state set is a permutation (e.g., [13]). Thus, a mutation is a specific instance of a permutation applied to a character.

The next step is to select the set Σ of admissible permutations acting on $\mathcal{C}$. It is mathematically convenient to select Σ such that it forms an Abelian group [9] with a regular (transitive and free) action on $\mathcal{C}$. Hence, Σ satisfies the following conditions:

C2. For every pair $c_i,c_j\in \mathcal{C}$ there is exactly one permutation σ ∈ Σ such that σ (c _i ) = c _j , i.e., the action of Σ on $\mathcal{C}$ is regular.

C3. For all σ₁,σ₂ ∈ Σ also the product σ₁ ∘ σ₂ ∈ Σ. Mathematically speaking, Σ is closed with respect to concatenation of its permutations.

C4. For all σ₁,σ₂∈Σ we have ${\sigma }_1\circ {\sigma }_2={\sigma }_2\circ {\sigma }_1$ Thus, Σ is commutative, and hence the order in which we assign permutations is irrelevant for the outcome.

C5. There is an element σ₀ ∈ Σ such that for all σ₁ ∈ Σ we have ${\sigma }_1\circ {\sigma }_0={\sigma }_0\circ {\sigma }_1={\sigma }_1$, i.e. there exists a so-called neutral element, namely the identity, in Σ. For all $c_i\in \mathcal{C}$ only σ₀(c _i ) = c _i , i.e. σ _i is fixed point free for all σ _i ≠σ₀.

C6. For every σ₁ ∈ Σ there exists a σ₂ ∈ Σ such that σ₁ ∘ σ₂ = σ₀. Mathematically speaking, for every element of Σ there exists an inverse element. This guarantees that every permutation can be reversed within a single step.

C7. For all σ₁,σ₂,σ₃ ∈ Σ we have σ₁ ∘ (σ₂ ∘ σ₃) = (σ₁ ∘ σ₂) ∘ σ₃ = σ₁ ∘ σ₂ ∘ σ₃, i.e. the associative law holds.

It should be noted that any set of permutations is associative, i.e. satifies C7. Thus, for a set of permutations Σ to be Abelian with a regular action on $\mathcal{C}$ it only needs to satisfy C1−C6.

In the following, we consider the matrix representation of permutations. A permutation matrix over $\mathcal{C}$ is an r×r matrix such that ${\sigma }_{c_i{c}_j}=1$ if σ(c _i )=c _j , and 0 otherwise. We consider it equivalent to discuss a permutation or its corresponding matrix. Therefore, concatenation “∘” is equivalent to the matrix multiplication “·”. We use σ to denote a permutation or a permutation matrix, depending on the context.

The set Σ_K3ST:={s₀s₁s₂s₃} coincides with the substitution matrices under the Kimura 3ST model[6]. In particular, s₁describes transitions within purines (A G) and pyrimidines (C T), s₂represents transversions within pairs (A C) and (G T), and s₃represents the remaining set of transversions within pairs (A T) and (C G).

Note that there are actually six different four-cycles for${\mathcal{C}}_{\mathrm{DNA}}$. These result in three distinguishable Abelian groups. Bryant[14]generates his cyclic group with the four-cycle A→C→G→T→A, and shows that the resulting set Σ_K2STunderlies the Kimura 2ST model[15], where${s^{\prime }}_2$corresponds to the transition within purines and pyrimidines, and${s^{\prime }}_1$and${s^{\prime }}_3$are the (not further distinguished) transversions.

This means that there are r ⁿ different operators in Σ ⁿ =Σ⊗⋯⊗Σ.

Remark 1. Therefore, for any pair of characters$f,g\in {\mathcal{C}}^n$we can find an operation σ ∈ Σ ⁿ such that σ (f) = g.

Another noteworthy consequence of using the Kronecker product is that the elements of Σ ⁿ are permutations over ${\mathcal{C}}^n$[16, 17], and in fact Σ ⁿ satisfies our Conditions C1−C7, i.e. Σ ⁿ is an Abelian group over ${\mathcal{C}}^n$.

where r is the number of states in Σ.

For every X-tree $\mathcal{T}$ we have ${\Sigma }^{\mathcal{T}}\supseteq {\Sigma }^X$, and therefore ${\Sigma }^{\mathcal{T}}$ is a generator for Σ ⁿ , too. An illustration of such a generator set ${\Sigma }^{\mathcal{T}}$ over the character set ${\mathcal{C}}^n$ is the so-called Cayley graph[18], which has as vertices the characters of ${\mathcal{C}}^n$, and two characters $f,g\in {\mathcal{C}}^n$ are connected if there is a permutation $\sigma \in {\Sigma }^{\mathcal{T}}$ such that σ(f)=g. In [5] Cayley graphs have been presented as alternative illustrations of the tree $\mathcal{T}$ over a binary state set $\mathcal{C}=\{0,1\}$.

Each permutation which acts on the characters is thus a symmetric 16×16 permutation matrix depicting a transition (s^e,1), transversion 1 (s^e,2), or transversion 2 (s^e,3) along edge $e\in E\left(\mathcal{T}\right)$. Figures 1b-d display the permutation matrices for a transition on branch e₁($s^{e_1,1}$), e₂($s^{e_2,1}$) and e₁₂($s^{e_{12},1}$), respectively. Figure 1e shows the Cayley graph associated with ${\Sigma }_{\text{K3ST}}^{\mathcal{T}}$.

$M_{\mathcal{T}}$ can be regarded as the weighted exchangeability matrix for all characters under the K3ST model assuming that a single substitution occurs on the tree $\mathcal{T}$. Figure 1f depicts the OSM matrix for the tree in Figure 1a. Here, colors indicate relative branch lengths p _e , and patterns denote permutation types α _i . E.g., a blue square with horizontal lines indicates the product $p_{e_2}{\alpha }_{e_2,1}$, i.e. the probability of observing a transition s₁on edge e₂.

The transformation problem

With the construction of ${\Sigma }^{\mathcal{T}}$ we have generated the tools needed to formally describe the computations in Step 4 of the Misfits algorithm [4]. Given a rooted tree $\mathcal{T}$ and two characters f and f ^d in ${\mathcal{C}}^n$, we want to compute the minimal number of substitutions required on the tree to convert f into f ^d . [4] presented an efficient procedure to compute this minimal number of substitutions.

Algorithm 1

INPUT: rooted binary phylogenetic tree $\mathcal{T}$ on leaf set X, characters f and f ^d on X, Abelian group Σ.

STEP 1: Using Remark 1, find the substitution type σ _i which translates f _j into $f_j^d$ for all positions j=1,…,|X|. Let σ∈Σ ⁿ be the resulting operation, i.e. σ(f)=f ^d .

STEP 2: Let c:=c₁…c₁be a constant character on X with $c_1\in \mathcal{C}$. Let h:=σ(c).

STEP 3: Calculate $m:=l_{\mathcal{T}}\left(h\right)$.

OUTPUT: m.

We prove the correctness of our algorithm. In our framework, m corresponds to the minimum number of permutations σ₁,…,σ _m ∈Σ such that σ₁⊗⋯⊗σ _m (f)=f ^d . In this form, m has multiple equivalent interpretations. It is the length of the shortest path between f and f ^d in the Cayley graph for ${\Sigma }^{\mathcal{T}}$, where this path corresponds to σ₁⊗⋯⊗σ _m . Further, m corresponds to the minimum power (k) of $M_{\mathcal{T}}$ such that $M_{\mathcal{T}}^j(f,f^d)=0$ for j<k and $M_{\mathcal{T}}^k(f,f^d)>0$, because a positive entry in $M_{\mathcal{T}}^k$ means that there is a concatenation of k permutations connecting the associated characters.

Example 3. Figure2demonstrates how Algorithm 4 works under the K3ST model, i.e. when the group is Σ = Σ _K3ST (Figure2a). Consider the rooted five-taxon tree in Figure2band the character GTAGA at the leaves. Assume that the character GTAGA is to be converted into character ACCTC. By comparing the two characters position-wise, we need a substitution s₁on the external branch leading to taxon 1 to convert G into A at the first position. Similarly, we need a substitution s₁on the external branch leading to taxon 2, and a substitution s₂on every external branch leading to taxa 3, 4, and 5. Thus, the operation s:=(s₁s₁s₂s₂s₂) transfers the character GTAGA into the character ACCTC. As the operation s also translates the constant character AAAAA into GGCCC, converting GTAGA into ACCTC is equivalent to evolving the character state A at the root along the tree to obtain the character GGCCC at the leaves. The Fitch algorithm[7]applied to the character GGCCC with the constraint that the character state at the root is A produces a unique most parsimonious solution of two substitutions as depicted by Figure2c.

The impact of parsimony on the estimation of substitutions

In this section, we provide some mathematical insights into the role of maximum parsimony in the estimation of the number of substitutions needed to convert a character into another one as explained above. In particular, we deliver a proof for Algorithm 4.

Theorem 1. Let $\mathcal{T}$be a rooted binary phylogenetic tree on taxon set X and let f be a character that evolved on$\mathcal{T}$due to some evolutionary model and let f ^d be another character on X. Then, the minimum number of substitutions to be put on$\mathcal{T}$which change the evolution of f in such a way that f ^d is generated can be calculated with Algorithm 4.

The impact of different groups

For any alphabet $\mathcal{C}$, there might be more than one Abelian group. Different groups might result in different numbers of substitutions required to translate a character into another character. We illustrate this observation using the following example.

Example 4. Recall the starting point of Example 3, i.e. regard the five-taxon tree T from Figure3b, and the characters f = GTAGA and f ^d = ACCTC. Now, instead of using Σ_K3STwe use the permutations from the cyclic group${\Sigma }_{{\mathbb{Z}}_4}$. In this setting, we need a substitution${s^{\prime }}_3$(blue in Figure3a) on the external edge leading to taxon 1 to convert G into A at the first position, and so on. Thus, we get the operation${s^{\prime }}_2:=({s^{\prime }}_3,{s^{\prime }}_1,s_3^{\prime },s_1^{\prime },{s^{\prime }}_3)$such that s ^′ (f) = f ^d . We immediately see, that s ^′ transforms the constant character c = AAAAA into h = CGCGC. The Fitch algorithm applied to the character CGCGC with the constraint that the character state at the root is A produces a unique most parsimonious solution of three substitutions as depicted by Figure3c. Thus, under the Σ _c group we need one substitution more than under the Σ_K3STgroup (cf. Example 3).

Note that variation of the minimum number of substitutions needed to translate a character into another one between different groups is not surprising: As different substitution types are needed to translate one pattern into the other one, depending solely on the underlying group, one group might need the same substitution type for some neighboring branches in the tree and another group different ones. Informally speaking, this would imply that in the first case, the substitution could be “pulled up” by the Fitch algorithm to happen on an ancestral branch, whereas in the second case this would not be possible.

The link between substitution models and permutation matrices

In Examples 1 and 2 we have shown that the K3ST substitution model can be included into our framework. The connection between the Klein-Four-group and the K3ST model (as well as the one between the ${\mathbb{Z}}_2$ group and symmetric 2-state model) were described in-depth in [9]. This section aims at discussing alternative models and how to identify their use (or lack thereof) for our approach.

Most substitution models assume the independence of the different branches of a tree to compute the joint probability of the characters in ${\mathcal{C}}^n$. Therefore, they use the probabilities for substitutions among the character states in $\mathcal{C}$ along the edges of the tree $\mathcal{T}$. We now establish a probabilistic link between ${\Sigma }^{\mathcal{T}}$ and ${\mathcal{C}}^n$. This link is provided by Birkhoff’s theorem:

Theorem 2 (Birkhoff’s theorem, e.g.,[20], Theorem 8.7.1). A matrix M is doubly stochastic, i.e., each column and each row of M sum to 1, if and only if for some N < ∞ there are permutation matrices σ₁, …, σ _N and positive scalars α₁, …, α _N ∈ [0,1] such that α₁ + ⋯ + α _n = 1 and M = α₁σ₁ + ⋯ + α _N σ _N .

Therefore, the weighted sum of the permutation matrices in ${\Sigma }^{\mathcal{T}}$ yields a doubly stochastic matrix $M_{\mathcal{T}}$ as introduced above. $M_{\mathcal{T}}$ also describes a random walk on ${\mathcal{C}}^n$ governed by $\mathcal{T}$ where the single step in ${\mathcal{C}}^n$ is illustrated by the associated Cayley graph. Its stationary distribution is uniform, i.e. when we throw sufficiently many mutations on $\mathcal{T}$ then we expect to see each character with probability 1/r ⁿ .

Another, even more useful consequence of Birkhoff’s theorem is the fact that it tells us which substitution models are suited for the OSM approach. If the transition matrix associated with the substitution model is doubly stochastic, then we find a set of permutations which give rise to the model.

which is equal to P_sGTRif a + b + c + d + e + f=1. Thus, the set Σ_sGTRis to sGTR what Σ_K3ST is to K3ST. However, Σ_sGTR does not satisfy condition C5, because s _a ,⋯, s _f are not fixed point free. This can be seen as the main diagonal of s _a ,⋯, s _f does not only contain zeros. It is also not commutative (condition C4) as e.g. s _a ·s _c ≠s _c ·s _a . And it is not closed under matrix multiplication (condition C3), which means that a concatenation of permutations in Σ_sGTR might lead to a new permutation not in Σ_sGTR, e.g., s _a ·s _f ∉Σ_sGTR. Other complex models like Tamura-Nei [21] do not even permit the decomposition of its transition matrix into the convex sum of permutation matrices. All of this shows why the overall applicability of complex models to the OSM approach is rather limited.

There are other approaches to describe phylogenetic models based on the group structure of their substitution matrices. In particular, Sumner et al. [22] use Lie algebra to construct OSM type matrices for the general Markov model, and discuss shortcomings of the group structure for the general GTR model [23].

Application to other biologically interesting sets

As stated above, OSM-type models require an underlying Abelian group. Thus, the OSM setting is applicable not only to binary data or four-state (DNA or RNA) data, but also to alphabets of 16 (doublets), 64 (codons), and 20 characters (amino acids) respectively. We compare such extensions to existing biologically motivated binning approaches and discuss their relevance.

As we have shown in the previous sections, the symmetric form of the Klein-Four-Group ${\mathbb{Z}}_2\times {\mathbb{Z}}_2$ is mathematically beautiful, computationally convenient and biologically relevant. Similar statements can be made about all powers of ${\mathbb{Z}}^2$, including the biologically relevant alphabets of 16 (doublets) and 64 (codons) letters.

There are four Abelian groups for twenty-state alphabets, namely ${\mathbb{Z}}_2\times {\mathbb{Z}}_2\times {\mathbb{Z}}_5,{\mathbb{Z}}_4\times {\mathbb{Z}}_5,{\mathbb{Z}}_2\times {\mathbb{Z}}_{10}$, and the cyclic group ${\mathbb{Z}}_{20}$ (see e.g., [24] for a complete list of all groups with up to 35 elements). Their construction is analogous to the construction of the Klein-Four-group in Example 1. For example, the elements of ${\mathbb{Z}}_4\times {\mathbb{Z}}_5$ are Kronecker products of one of the four permutations in the cyclic group ${\mathbb{Z}}_4$ with one of the five permutations of the cyclic group ${\mathbb{Z}}_5$.

Figure 4 shows a heat-map type visualization of an OSM-type matrix on a single-leaf tree where the coloring of the cells corresponds to the weights given to the 20 permutations in the respective groups. We see that the coloring pattern nicely reflects the four cosets of the subgroup ${\mathbb{Z}}_5$ in ${\mathbb{Z}}_2\times {\mathbb{Z}}_2\times {\mathbb{Z}}_5$. This can also be interpreted as a binning of the 20 states in the underlying alphabet into four sets of five elements each. If the weighting corresponds to a convex combination of operations, then the visualized matrix is doubly stochastic.

Binnings are also done for amino acids, using either biochemical properties or evolutionary divergence. An example of a biochemical binning is the hydrophobic index, where the 20 amino acids are binned into four groups, very hydrophobic, hydrophobic, neutral, and hydrophilic. Unfortunately, this binning does not correspond to any of the proposed Abelian groups. Moreover, it is difficult to derive transitions between these groups just from the biochemical properties.

Transition matrices for evolutionary models for amino acid substitutions are usually generated by counting mutation types in the alignments (see, e.g., [25] for an overview). From these, optimal groupings can be obtained using clustering approaches [26]. The existence of estimates for the transition probability between all amino acids provides the possibility to get further information about between-group operations. These groupings could be forced to fit Abelian groups. However, as indicated in [26] a grouping into four groups of five amino acids each is rarely optimal.