A tensorial approach to the inversion of group-based phylogenetic models

Background Hadamard conjugation is part of the standard mathematical armoury in the analysis of molecular phylogenetic methods. For group-based models, the approach provides a one-to-one correspondence between the so-called “edge length” and “sequence” spectrum on a phylogenetic tree. The Hadamard conjugation has been used in diverse phylogenetic applications not only for inference but also as an important conceptual tool for thinking about molecular data leading to generalizations beyond strictly tree-like evolutionary modelling. Results For general group-based models of phylogenetic branching processes, we reformulate the problem of constructing a one-one correspondence between pattern probabilities and edge parameters. This takes a classic result previously shown through use of Fourier analysis and presents it in the language of tensors and group representation theory. This derivation makes it clear why the inversion is possible, because, under their usual definition, group-based models are defined for abelian groups only. Conclusion We provide an inversion of group-based phylogenetic models that can implemented using matrix multiplication between rectangular matrices indexed by ordered-partitions of varying sizes. Our approach provides additional context for the construction of phylogenetic probability distributions on network structures, and highlights the potential limitations of restricting to group-based models in this setting.


Background
Fundamental to evolutionary biology is the development and implementation of molecular phylogenetic methods [1]. These methods provide the means to reconstruct the past evolutionary history of biological entities given present-day molecular data, such as DNA. Considering Kimura's neutral theory of molecular evolution, it is logical to apply a stochastic model at the level of DNA substitutions to construct probabilistic description of what molecular alignments are expected to be observed, given a proposed evolutionary history (tree topology and edge lengths). is commonly implemented assuming an IID (across sites in the alignment) and Markov process for DNA substitution, leading to a model that has a continuous-time Markov chain at its core (see Semple and Steel [2] for an introduction to the mathematics underlying modern phylogenetic methodology). In a series of papers, Hendy and colleagues introduced the Hadamard conjugation as a novel tool for phylogenetic analyses [3][4][5]. They found an invertible relationship between a phylogenetic tree, as characterized by its edge length spectrum, and the probability distribution of site patterns (referred to as the sequence spectrum). Originally introduced only for the 2-state symmetric model, the Hadamard conjugation was later extended to the K3ST model [6][7][8] and further to any of the so-called "groupbased" models [9]. Hadamard conjugation has been used as both a tool for simulation [10] and to look at statistical properties of methods, exploring the inconsistency of parsimony under a molecular clock [5,11]. For these sorts of applications, following the notation in Felsenstein [1], we can use the Hadamard transform H to start with an edge length spectrum γ and calculate the sequence spectrum s = H −1 log(Hγ). The beauty of Hadamard conjugations is that one can also begin with an observed sequence spectrumŝ and perform the inverse of the conjugation to empirically obtain an edge length spectrum http://www.biomedcentral.com/1471-2148/14/236 γ = H −1 log(Hŝ). Although it is not expected that thê γ spectrum will precisely match a tree, Hendy [12] proposed using an optimisation criterion to map fromγ to the "closest tree".
Several authors have commented that it is potentially a useful feature of Hadamard conjugation that data isn't forced onto a fixed tree. The conflicting information can be retained and interpreted in the form of a "lentoplot" [13] or a splits-graph [14], with both of these methods implemented in Spectronet [15]. Schliep [16] gives some more statistical justification for such an approach by making a link to modern statistical techniques such as the Lasso and Ridge regression.
von Haeseler and Churchill [17] seems to be the first paper that explicitly suggests using Hadamard conjugation to provide a likelihood framework for networks. The principle idea in this work was to start with an edge length spectrum that encodes a set of incompatible splits, use the Hadamard transformation to get site probabilities and use these to determine a likelihood. This idea was further explored by Bryant [18], and Bryant [19] followed this through defining the "n-taxon process" for group-based models. It should be noted that likelihoods calculated via Hadamard are not equivalent to likelihoods calculated by taking a mixture of trees. Indeed, Matsen and Steel [20], Matsen et al. [21] used Hadamard methods in combination with phylogenetic invariants to show that mixtures of trees with the same topology can exactly mimic another tree under the 2state model. Considering biological applications, thinking in terms of mixtures of trees or partitions where the data can be thought of as arising on a set of trees [22][23][24] seems more reasonable than the Hadamard conjugation. Strimmer and Moulton [25] suggested using split networks as a spring board to likelihood-based analyses on DAGs, but later identified several problems with the approach [26]; most notably, in split-networks internal nodes do not have a biological interpretation as an ancestor.
In Sumner et al. [27], we gave some additional insight into the interpretation of applying the Hadamard conjugation in a network setting. We showed that permutation group structure inherent to the Hadamard transformation -as for any group-based model -restricts the resulting process from being capable of reproducing truly convergent processes. This is a serious limitation, as one of the biological motivations for explicit network models is the ability to model convergent processes. We also presented an alternative algebraic formalism for the general Markov model, analogous to the n-taxon process, but capable of reproducing convergent processes.
From the point of view of group representation theory, the inversion of group-based models relies on the fact that the irreducible representations of an abelian group are one-dimensional, and the model structure essentially reduces to analysing group characters -hence the standard presentation of a Fourier inversion. In this article, we make this connection concrete. For the general Markov model, it is then immediately apparent that an analogous inversion is not possible because the algebraic structure underlying the model is not abelian and hence the irreducible representations are not one-dimensional. In fact, to obtain one-dimensional representations for the general Markov model, it is necessary to apply higher-degree polynomial maps (beyond the degree 1, linear case), and define "Markov invariants" [28]. These invariants present onedimensional representations but at the cost of the higher degree -degree 5 in the case of the general Markov model with four states on quartet trees [29,30]. This connection between Hadamard transformation and Markov invariants is an interesting one, but we do not discuss it further here.
In this paper we approach the inversion of group-based phylogenetic models by taking a representation-theoretic perspective and working explicitly with tensor indices. Our approach rests heavily on the formalism of "phylogenetic tensors", as presented in Bashford et al. [31], for the binary-symmetric and K3ST model, and Sumner et al. [27,28], for the general Markov model.
Although the main inversion results presented here are not more general than those in in Székely et al. [7], we think it is important to reformulate them using the language of tensors and representation theory. This viewpoint has already led to new approaches for modeling convergent evolution [27] and for studying non-groupbased models [28]. However, in none of our previous work was the link to Hadamard conjugation explicitly discussed. By presenting an old technique (Hadamard conjugation) in a new light we hope to introduce other researchers to the viewpoint of tensor analysis and representation theory.

Group-based models
We consider the continuous-time formulation of Markov processes, and show how to implement the inversion of a group-based phylogenetic model based on any abelian group G. We note that such an inversion requires a map from tensor product space (where elements are indexed by ordered-n-partitions) to phylogenetic splits (where elements are indexed by bipartitions). We achieve this by finding canonical maps from bipartitions to ordered-npartitions.
For a group G (not necessarily abelian) with order |G| = d, we write G = {σ 1 , σ 2 , . . . , σ d }, and, when necessary, write ∈ G to specify the identity element of G. We will discuss the "regular representation" of G shortly, but http://www.biomedcentral.com/1471-2148/14/236 skipping ahead we find that any rate matrix Q occurring in a group-based Markov model can be written in the form where each 0 ≤ α σ ∈ R, λ = =σ ∈G α σ and the K σ are the permutation matrices corresponding to the (nonidentity) group elements σ ∈ G.
For the reader interested in deriving this result, with scalar multiplication and vector addition defined via for all v, v ∈ G C and λ ∈ C. The regular representation, ρ reg : G → GL(d, C), is then defined by setting the group action for all v ∈ G C and σ ∈ G. If we fix {σ 1 , σ 2 , . . . , σ d } as an ordered basis for G C , it is then clear -via Cayley's theorem -that each group element σ gets mapped to a permutation matrix Consider the unit column vectors and identify each σ i ∈ G C with ξ i ∈ C d , so that the group action becomes σ : i have i as the column label and j as the row label.
A group-based Markov model is then obtained by taking a continuous-time Markov chain with state space G = {σ 1 , σ 2 , . . . , σ d } and using the group multiplication in G to assign a rate α σ to all substitutions σ 1 → σ 2 where σ σ 1 = σ 2 . Following this through (as is done in detail in [32]) we are led to the formula (1) for rate matrices in any group-based model.
The regular representation is one example of the general concept of a representation of G on a vector space V, defined as a homomorphism ρ : G → GL(V ) satisfying ρ(g 1 g 2 ) = ρ(g 1 )ρ(g 2 ) for all g 1 , g 2 ∈ G. A representation is said to be reducible if there exists a proper subspace U ⊂ V satisfying ρ(g)U ⊂ U, i.e. the set of matrices ρ(G) send vectors in U back to U. In this case, U is called an invariant subspace. The representation ρ is then called irreducible if V does not contain any invariant subspaces.
The reader should note that the usual construction of a "group-based" model [2] stipulates that G be abelian. Although the construction just given using the regular representation allows for non-abelian G, we will nonetheless only consider the abelian case in this paper, because, as discussed in the introduction, it is only in the abelian case that a (linear) inversion of phylogenetic models is possible. In this case the irreducible representations of G are all one-dimensional [33], and hence the analysis reduces to computations with group characters, as is exploited in the previous approaches using Fourier analysis [9,34].

Phylogenetic tensors
We denote [d] := {1, 2, . . . , d} as the state space for a continuous-time Markov chain. Consider an n-taxa phylogenetic tree and a d-state phylogenetic pattern distribution {p i 1 i 2 ...i n } i 1 ,i 2 ,...,i n ∈[d] with the interpretation that p i 1 i 2 ...i n is the probability that the observed state at the k th leaf on the tree is i k . As is shown in Sumner and Jarvis [35] and in more detail in Sumner et al. [27], such phylogenetic pattern distributions can be represented abstractly as tensors in the n-fold tensor product space ⊗ n C d := C d ⊗ C d ⊗ . . .⊗ C d , as follows. If we choose {ξ 1 , ξ 2 , . . . , ξ d } as an ordered basis for C d , and ordered basis for the tensor product space, a "phylogenetic tensor" P ∈ ⊗ n C d is then defined as For readers who are unfamiliar with tensor products, it is possible to understand the general concept via the definition of the "Kronecker" product of a n × m matrix A and a n × m matrix B as the nn × mm matrix given by We can index the matrix A ⊗ B with row indicies i 1 j 1 = 11, 12, . . . , nn and column indices j 1 j 2 = 11, 12, . . . , mm , i.e. generically (A ⊗ B) i 1 j 1 ,i 2 j 2 = A i 1 i 2 B j 1 j 2 and specifically (A ⊗ B) 12,32 = A 13 B 22 . This point of view is useful if one wants to write out specific matrix representations of tensors, however, in the development that follows will http://www.biomedcentral.com/1471-2148/14/236 focus heavily on the indexing of tensor components in the various cases discussed.
Suppose π = i∈[d] π i ξ i ∈ C d represents the state distribution of a single taxa, i.e. π i is the probability that a randomly chosen site in the sequence will be in state i. Now suppose a phylogenetic branching event occurs and the sequence is copied. The corresponding phylogenetic representing the joint distribution of the two-taxa just after the branching event then has the property that p i 1 i 2 = π i 1 if i 2 = i 1 and is zero otherwise. Thinking in terms of tensor operations, we find that phylogenetic branching events can be generated by a linear operator δ : C d → C d ⊗ C d determined by δ(π) = P and defined in general using our chosen basis as The remarkable fact for group-based models, central to the present article, is that the permutation matrices "intertwine" particularly simply with the branching operator: Thus, for any rate matrix Q arising from a group-based model, we have (via the linearity of δ): We also note that, since Q can be expressed a linear combination of permutation matrices representing elements in a group G, the matrix powers Q 2 , Q 3 , Q 4 . . . will also be expressible as linear combinations of the same permutation matrices (although precise expressions for the relevant coefficients may or may not be easily computable). Together with (3), this implies that, for any substitution matrix e Qt arising from matrix exponentiation, This relation shows that mathematically, and hence conceptually, "Markov evolution on a single followed by a branching event" can be replaced with "Branching event on a single taxon followed by (correlated) Markov evolution of two taxa. " This equivalence is illustrated in Figure 1, and should be compared to the equivalent discussion of the "n-taxa process" given in [18] and [19].

Figure 1 Markov evolution on a single followed by a branching event (illustrated on the left), is equivalent to a branching event on a single taxon followed by correlated Markov evolution of two taxa (illustrated on the right).
Mathematically, this equivalence can be implemented by exploiting the equality given in (4).
In Sumner et al. [27] we showed how to generalise this intertwining action to the case of the general Markov model. Interestingly, for the general Markov model the appropriate intertwining has quite a different structure from what occurs in group-based models, and hence the simplicity of (4) is somewhat misleading in general. We refer the reader to Sumner et al. [27] for more discussion on this point.
Returning to the case of group-based models, for each subset A ⊆ [n], we define a linear map on ⊗ n C d as the tensor product K (A) σ := K a 1 σ ⊗K a 2 σ ⊗. . .⊗K a n σ where a i = 1 if i ∈ A and 0 otherwise. For example, if n = 5, we have To develop a phylogenetic tensor on a tree, we root the phylogenetic tree at taxon n, and label edges by subsets ∅ = e ⊆ [n − 1], where i ∈ e if the path from taxon n to taxon i crosses the edge labelled by e. A five taxon tree with this labelling, is presented in Figure 2. To each edge labelled by ∅ = e ⊆ [n − 1], we assign the rate matrix where each α σ e ≥ 0 is the rate of substitution for all states σ 1 to σ 2 satisfying σ = σ 2 σ −1 1 , and λ e = σ ∈G α σ e . Each edge is then assigned substitution matrix M e = e Q e , so http://www.biomedcentral.com/1471-2148/14/236 that the time parameter for each edge is absorbed into the definition of Q e . Now iterating (4) multiple times, Bashford et al. [27,31] show that any phylogenetic tensor can be written as where λ = ∅ =e⊆[n−1] λ e = ∅ =e⊆[n−1], =σ ∈G α σ e , and δ n−1 π is the d × d × . . . × d tensor that represents the "zero edge-length star tree" distribution on n taxa. It is this form of phylogenetic tensors that will do a lot of the heavy lifting in the discussion that follows. The reader should note that under this representation, there is no need for the edge parameters α σ e : ∅ = e ⊆ [n − 1], σ ∈ G to be chosen to be compatible with a particular tree, hence the possibilities for generalising to non-tree-like or network models, as discussed in the introduction.
The stationary distribution for group-based models is uniform (because the rate matrices are doubly stochastic). In this paper we always assume a stationary distribution, so that: and δ n−1 π has tensor components This concludes our discussion of the tensor presentation of phylogenetic probability distributions under group-based models. It is important to note that everything discussed so far works for any group-based model, with no requirement that the underlying group G be abelian.
In what follows, we discuss the inversion of abelian group-based models. We present the simplest case with G = Z 2 ; the G = Z 3 case; the G = Z 2 × Z 2 case; the general G = Z r case; and finally we discuss the case of any abelian group.

The binary-symmetric case
We begin with the inversion of the so-called "binarysymmetric" model. Consider C 2 with standard basis As a group-based model, the binary-symmetric model arises by taking the group with a generic rate matrix given by where K = 0 1 1 0 is the permutation matrix representing σ in the standard basis. Now ρ reg : Z 2 → M 2 (C), with σ → K, is the regular representation of Z 2 , and the character table of Z 2 given in Table 1 is easily recognised to be the Hadamard matrix As Z 2 is an abelian group, the irreducible representations are one-dimensional.
The corresponding projection operators can be read off from the columns of the character table. That is, the operators project ρ reg = id ⊕sgn onto the id and sgn representations of Z 2 , respectively. This observation prompts us to work in the alternative basis: In this basis the permutation matrix is diagonal: The representation-theoretic perspective on K is to observe that id(σ ) = 1 and sgn(σ ) = −1.
Referring to (5), we know that we can write a generic phylogenetic tensor as We index matrix and tensor indices by using i, j, k = 0, 1 ∈ Z 2 and allow multiplication × in the ring of integers Z. The Hadamard matrix then has matrix elements [h] j i = (−1) i×j where j is the row index and i is the column index. Observe that in the diagonal basis, the permutation matrix has elements Thus we have expressions such as As we are dealing with tensors of arbitrary size, it is convenient to represent a string such as i 1 i 2 . . . i n as an ordered-bipartition μ = μ 0 :μ 1 (1) := h, in the diagonal basis P := h (n) · P and using our notation h (n) has tensor components The zero edge-length star-tree initial distribution has tensor components (where, although it seems we have given preference to taxon 1 in this expression, there are many ways that this distribution can be expressed using the δ ij ). In the diagonal basis with δ n−1 π := h (n) · δ n−1 π, we have components which is exactly the statement δ n−1 π μ = δ n−1 π μ 0 :μ 1 Since K is diagonal in the transformed basis, we can conclude that [2,n] α e K (e) μ 0 :μ 1 Of course many of these tensor components will be zero and we would like to ignore these.
Take u = u 0 :u 1 as an ordered bipartition of the reduced set [n − 1], so that u ≡ i 1 i 2 . . . i n−1 where j ∈ u k if and only if i j = k, and define then we can write the non-zero components as This is the first part of the inversion. We would like to go further and actually recover the individual edge weights α e . To do this we define the (square) 2 n−1 × 2 n−1 matrix F with components with e a subset and u an ordered-bipartition of [n − 1]. As h (n−1) 2 = 1 2 n−1 1, we see that F provides its own inverse F −1 with components Defining the column vectors α = {α e } and η = {η u }, we can write the matrix equations Together with the first part of the inversion (6), these equations give a one-one map between pattern probabilities and edge weights for the binary-symmetric model.

Inversion of the Z 3 model
Taking confidence from the previous case we now discuss the inversion of the group-based phylogenetic model with G = Z 3 . We take and, by analogy to the Z 2 case, index tensors with indices i, j = 0, 1, 2 and allow multiplication × by extending Z 3 to the ring F 3 = {0, 1, 2} +,×(mod 3) .
In this case a generic rate matrix is given by We define ω = e 2πi/3 , and present the character table of Z 3 in Table 2. The decomposition of the regular representation is ρ reg = id ⊕ ω ⊕ ω 2 , and the columns of the character table give the projection operators onto the (one-dimensional) irreducible subspaces: Therefore, the matrix diagonalizes the generic rate matrix for this model: or, equivalently, We recall our basic result (5) that for group-based models, a generic phylogenetic tensor can be expressed as where λ = ∅ =e⊆[n−1] (α e + β e ). We take the stationary distribution as initial distribution, so π = 1 3 , 1 3 , 1 3 T .
If |μ 1 | + 2|μ 2 | = 1 (mod 3), then and if |μ 1 | + 2|μ 2 | = 2 (mod 3), then Thus we have found a basis where all the elements of the initial star-tree tensor are zero unless the tripartion μ satisfies |μ 1 |+2|μ 2 | = 0 (mod 3). Crucially, this statement also holds for the phylogenetic tensor P because in this basis the rate matrices of this model are diagonal: 1 +β e K (e) 2 μ 0 :μ 1 :μ 2 μ 0 :μ 1 :μ 2 We deal with this condition on μ by taking u = u 0 :u 1 :u 2 as an ordered-tripartion of the reduced set [n − 1] and setting μ = u · γ (u) (considered as the concatenation of strings) where If we make the definitions P u := P u·γ (u) , , we then have the first part of the inversion As in the Z 2 case, we would like to use η u to recover the rate parameters α e , β e for all ∅ = e ⊆ [n − 1] and thus complete the full inversion for this model. Of course, it is little bit more difficult this time.
Recall that μ = μ 0 :μ 1  Considering for all ordered-triparitions u, w of [n − 1], we have the matrix products Thus the second part of the inversion for this model is Together with (7), these equations give a one-one map between pattern probabilities and edge weights for the group-based model with G = Z 3 .
In this model a generic rate matrix is given by We already know that the 2 × 2 Hadamard matrix h diagonalizes K, so we see immediately that H = h ⊗ h diagonalizes this model: Of course H is the character table of Z 2 × Z 2 and the permutation matrices (8), together with K 00 := 1, give the regular representation ρ reg ∼ = id⊗id⊕id⊗sgn⊕sgn⊗id⊕ sgn ⊗ sgn, where we recall the basic result that the tensor product of two irreducible representations of a group G gives an irreducible representation of G × G. http://www.biomedcentral.com/1471-2148/14/236 Simplifying notation, for this model we index tensors with indices given as pairs: i, j = 00, 01, 10, 11 ∈ Z 2 × Z 2 ; and we express the individual parts using lower case Roman characters. For example, we write i := ab = 01, with a = 0 and b = 1. This gives matrix elements: and more complicated tensor products such as Again we interpret strings such as μ ≡ a 1 a 2 . . . a n and ν ≡ b 1 b 2 . . . b n as ordered-bipartitions μ = μ 0 :μ 1 and ν = ν 0 :ν 1 of the set [n]. We can then write matrix elements of tensor products as Taking the stationary distribution π = 1 4 (1, 1, 1, 1) T as initial distribution, the zero edge-length star-tree distribution is given by which in the finer index representation is δ n−1 π a 1 b 1 a 2 b 2 ...a n b n = 1 4 δ a 1 a 2 δ a 1 a 3 . . . δ a 1 a n δ b Recall that elements of the Hadamard matrix can be written as [h] a b = (−1) a×b , where a, b ∈ Z 2 and we allow multiplication × by extending to the ring of integers Z. In the transformed basis, we have a 1 a 2 δ a 1 a 3 . . . δ a 1 a n δ b 1 b 2 δ b 1 b 3 . . . δ b 1 ..a n )×c 1 +(b 1 +b 2 +...+b n )×d 1 = 1 4 1 + (−1) a 1 +...+a n + (−1) b 1 +...+b n +(−1) a 1 +...+a n +b 1 +...+b n = 0, if either |μ 1 | or |ν 1 | is odd; 1, if |μ 1 | and |ν 1 | are both even.
To exclude the vanishing components we define, for all ordered bipartitions u = u 0 :u 1 of the reduced set [n − 1], and intepret u · γ (u) as the string u · γ (u) = a 1 a 2 . . . a n−1 γ (u). Then, for each pair u, v of ordered-bipartitions of [n − 1], we define This gives the inversion Consider the 2 n ×2 n−1 rectangular matrices F 01 , F 10 and and completes the inversion for the K3ST model.

Inversion of the Z r model
We now consider the group based model for Z r = {0, 1, 2, . . . (r − 1)} +(mod r) ∼ = σ : σ r = e . For this model the generic rate matrix has the form Lemma 1.
The regular representation contains exactly one copy of every irreducible representation and the irreducible representations of Z r are given by the powers of ω: Thus the change of basis  Translating this result using the ordered-r-partitions for indices, we have Lemma 3. In the diagonal basis, the uniform initial distribution on the star tree has components where μ = μ 0 :μ 1 :μ 2 : . . . :μ r−1 is an ordered-r-partition of the set [n].
Again recall that for this model a generic phylogenetic tensor can be written as where π = 1 r (1, 1, . . . , 1) T . In the diagonal basis P := f (n) · P and as a consequence of Lemma 3 P will have many vanishing components. To avoid these we take u = u 0 :u 1 :u 2 : . . . :u r−1 as an ordered-r-partition of [n − 1] and set γ (u) = r − (0|u 0 | + 1|u 1 | + 2|u 2 | + . . . If we define P u := P u·γ (u) and , we then have the first part of the inversion for the Z r model: , so we have the vector equation We claim that .
Of course G i F j = δ ij 1, so we now have the second part of the inversion:

Inversion of any abelian group-based model
Lemma 5. Any (finitely generated) abelian group G is isomorphic to a direct product of cyclic groups of primepower order, ie. G ∼ = Z r 1 × Z r 2 × . . . × Z r q where each r i = p n i i where p i is prime and n i is a positive integer. http://www.biomedcentral.com/1471-2148/14/236 Lemma 6. The group-based model arising from the G is defined only up to group isomorphisms of G.
Proof. A generic rate matrix for the group-based model arsing from G is given by Recall (2), so that the matrix elements [K σ ] j i is set via the action σ i → σ σ i = σ j . If we consider the regular representation of G we then have and, because φ is a group isomorphism, this occurs if and only if σ This means that we can restrict attention to a single representitive in the isomorphism class of G. Of course, for this purpose we choose the representative guaranteed by Lemma 5.
Thus, for any abelian group G, with generators σ 1 , σ 2 , . . . , σ q the corresponding group-based model has rate generators given by for all e = σ = σ m 1 1 , σ m 2 2 , . . . , σ m q q ∈ G, where K σ i is the permutation matrix representing the generator σ i ∈ Z r i .
The character table f of G is simply the tensor product of the individual character tables of the Z r i : In the diagonal basis we have matrix elements where ω k is a k th root of unity. Thus

Conclusion
In this article we have given an alternative derivation of the inversion of group-based phylogenetic models. Primarily our method relies on the remarkable intertwining relation between branching events and Markov evolution (4), and the resulting simplified expression of phylogenetic tensors given in (5). From there we took a representation theoretic approach concentrating on the structure of tensor indices.