A representation-theoretic approach to the calculation of
evolutionary distance in bacteria
Jeremy G Sumner, Peter D Jarvis, and Andrew R Francis

arXiv:1612.06035v1 [q-bio.PE] 19 Dec 2016

December 20, 2016
Abstract
In the context of bacteria and models of their evolution under genome rearrangement, we
explore a novel application of group representation theory to the inference of evolutionary
history. Our contribution is to show, in a very general maximum likelihood setting, how
to use elementary matrix algebra to sidestep intractable combinatorial computations and
convert the problem into one of eigenvalue estimation amenable to standard numerical
approximation techniques.

1

Phylogenetics and bacterial evolution

Phylogenetics is the suite of mathematical and computational methods which provide biologists
with the means to infer past evolutionary relationships between present-day observed species.
The usual computational output is an inferred evolutionary tree, together with branch lengths
and dates of divergences. Input to these methods most typically consist of molecular data such
as sequence alignments (for eukaryote DNA, the basic datum is the pattern of nucleotides at
each site in the alignment). The state of the art is to model molecular evolution as a continuoustime Markov chain on a tree, with the content of each site in the alignment presumed to have
evolved under the statistical assumption of being independent and identically distributed (IID).
Proposed evolutionary histories are inferred using standard statistical frameworks, such as
maximum likelihood [14, 20], or a Bayesian approach [9, 39] (in the latter the output is not just
a single tree proposal, but rather a posterior distribution on all model parameters, including
the evolutionary tree itself). These methods are only feasible because of significant modern
developments in both stochastic modelling and the almost weekly rise of available computer
power.
Historically, more elementary approaches to phylogenetic estimation have included maximum parsimony [12, 16] and distance-based methods [29, 33]. Maximum parsimony proceeds by
asking, given a sequence alignment and a candidate evolutionary tree, what is the least number
of changes required to fit the observed nucelotide patterns to the leaves of the tree. This method
is clearly blind to multiple changes and has long been known to be statistically biased (most
notably via the well-studied “long-branch attraction” effect [13]). For this reason, parsimony
has somewhat fallen out of favour and has been replaced by model-based approaches. On the
other hand, distance-based methods are a sensible intermediary, as they can be implemented in
a statistically consistent manner (with the assumption of a Markov model on a tree), and computationally are very efficient in comparison to maximum likelihood or Bayesian approaches.
Indeed, the main difficulty for any phylogenetic method is the huge number of possible (binary) trees: for L taxa there are ∼ 2L L! trees1 , and the attractive feature of distance-based
1

There are (2L − 3)!! = (2L − 3)(2L − 5) . . . 3 · 1 rooted binary trees with L leaves.

1

methods is that they proceed by clustering and hence can create a sensible candidate tree in
O(L3 ) time [15]. These include the ever-popular “Neighbour-Joining” algorithm [29], and its
generalizations “Neighbour-Net” [5] and “Splits-Tree” [22]. The latter methods return a network structure which is close to a tree but allows for the visual representation of signals caused
by evolutionary events which deviate from strict vertical descent (for example, “lateral gene
transfer” events2 ) and/or random noise in the sequence data. Importantly, current standard
practice for both likelihood and Bayesian approaches is to first use a distance-based method
to construct a reasonable starting tree and then use hill-climbing on the likelihood function to
attempt to find a better solution using small perturbations to the tree under standard set of
moves [4].
At their core, distance-based methods rely on finding a statistically consistent method for
calculating the distance between a pair of observed sequences (sometimes referred to as a
“genetic” distance). Solving this problem credibly consists of three parts:
(i) Specify a stochastic model of how molecular sequences evolve randomly, for example a
fixed Markov model on a tree;
(ii) Define mathematically what is meant by distance relative to the model. For example, in a
continuous-time formulation the distance is most sensibly made to be time itself (although
this is not the only choice and may depend on extra additional unknowns such as relative
rates of molecular substitutions);
(iii) Analyse the model to produce a robust and bias free statistical estimator of the distance
from observed data.
A simple example illustrating this process is the well-known “Jukes-Cantor correction” [23]
which (i) assumes a Markov model where every nucleotide substitution occurs with the equal
rate r; (ii) defines distance to be the time passed; and (iii) corrects for unobserved multiple
changes in the raw number of observed changes ∆ between two sequences of length N (also
known as the Hamming distance). This is achieved with a statistically consistent estimate of
time elapsed via the inversion formula Tb = −(3/4r) log(1 − 43 ∆/N), which can be shown is the
analytic expression for the maximum likelihood estimate of time under the Jukes-Cantor model
[18]. Analogous to maximum parsimony, failing to use this correction, for example applying
the Hamming distance ∆ directly, can lead to incorrect evolutionary tree inference.
The Jukes-Cantor correction can be generalized to distances consistent with more complicated models such as the Kimura 2ST model [24], or even the general Markov model [2] using
the so-called “log-det” [26]. However, as always in any statistical estimation problem, there
is a compromise between assuming a complicated model to increase model realism (and hence
decreasing bias), with the inevitable increase in variance of the additional parameters requiring
estimation in more complex models [6]. This so-called “bias-variance” tradeoff is an extremely
important consideration in all approaches of phylogenetic estimation methods.
Phylogenetic modelling in the context of bacterial genomes is radically different. At the
DNA level, evolution in bacteria is far more dynamic than that in eukaryotes, including evolutionary events that violate the assumption of vertical descent. These include genome rearrangements (relocation of genes within the genome), deletions of large segments of the genome, and
lateral gene transfer between species. This fluidity makes the methods described above – the
observation of substitutions at single nucleotide sites to infer evolutionary distances – highly
problematic for bacteria. A major practical obstruction is that accurate sequence alignment
of bacterial genomes at individual nucleotide sites is unattainable. For instance, it is often the
case that genomes within the same bacterial species share as little as 40% of the same genes
(the situation with E.coli is discussed in [37]).
2

Where genes from one species are inserted into the genome of another.

2

The occurrence of large-scale rearrangements deems the previously described classic phylogenetic approaches (where individual nucleotide sites are modelled under an IID assumption)
rather inappropriate for bacteria, and a more coarse-grained approach is required. In this context, it is sensible to divide bacterial genomes into “regions” (these can be genes or clusters
of genes; essentially any component which is sufficiently conserved such that identification of
the region is plausible), and then compare similarity between genomes by comparing firstly
how many regions the bacteria have in common, and secondly the relative arrangement of the
regions on the genome.
A major drawback for existing models of bacterial evolution is the absence of a statistically consistent distance estimator analogous to the Jukes-Cantor correction. In principle there
is no difficulty for, as we will describe in the next section, it is straightforward to develop a
stochastic model of genome rearrangement and produce an estimator for evolutionary distance
which is statistically consistent with the model. However, as we will also discuss, in practice
the calculations involved are of factorial complexity (in the number of regions) and, to date
only provisional approaches have been explored. These approaches either compromise on model
realism (making strong assumptions regarding the types of rearrangements that can occur), or
compromise on statistical consistency, or both. For instance, the means to compute, in polynomial time, the minimal number of “signed inversions” required to convert one genome into
another was provided in [21] (this was improved to linear time in [1]). Comparing to the DNA
case, this is analogous to computing a Hamming (or edit) distance, and is usually referred to
as “minimal distance” in the context of genome rearrangement models. The results given in
[21] are however limited by strong assumptions on the genome rearrangements that may occur
(specifically, each signed inversion is assumed to be equally likely — an assumption that is
known to be violated, see [8] for example). Earlier, [30] presented the so-called breakpoint distance. This provides a minimal distance between two genomes that is very easy to calculate but
has the disadvantage of not being closely aligned to a specific model of genome rearrangement.
Using simulation, an empirically-derived correction to the breakpoint distance was presented in
[36], and in [32] we explored the importance of specifying a stochastic model of rearrangements,
explicitly showing there are cases where the maximum likelihood estimate of time elapsed and
the minimal distance produce contradictory answers.
In this paper we work with a model based on a simple stochastic process that allows for
an arbitrary set of genome rearrangements. Our original contribution is to explore the efficient
calculation of estimates of elapsed evolutionary time under this model using methods from group
representation theory. We discuss how to convert the problem into the calculation of group
characters and hence show the core computational difficulty reduces to an eigenvalue problem.
This approach has an immediate pay off in the reduction of computational complexity from
factorial complexity to the square root of factorial complexity (as a conservative estimate).
Although this is a useful theoretical observation it of course does not solve the computational
issues inherent in the problem. Our broader goal is to use these ideas to produce the means to
find numerical approximations to maximum likelihood estimates of elapsed time, without the
need to make unrealistic modelling assumptions. Ultimately we hope our approach will lead
to computationally efficient, usable software for inferring bacterial phylogenetic trees under
rearrangement models.
Our mathematical methods are nothing more than undergraduate level linear algebra and
standard results from the representation theory of finite groups (we consider [28] an excellent
reference for the later). It is worth noting that the results on rearrangement distances presented
in [11] feature application of the irreducible characters of the symmetric group, and hence we
are not the first authors to apply character theory to the problem of genome rearrangement
models. However, the results provided in that work are again highly specialized as it is assumed
that each transposition of genome regions is equally likely. To the best of our knowledge, the
3

work presented in this paper is the first to apply more general methods from representation
theory, and certainly the first to do so in the context of likelihood calculations under a general
model of genome rearrangement.
In Section 2 we discuss the details of the basic mathematical framework for setting up a
genome rearrangement model. In particular, but without loss of generality, we discuss how to
deal with the case of “unsigned” circular genomes. In Section 3 we use the representation theory
of the symmetric group to find a novel means of counting passages between genomes under a
given rearrangement model. We explore this idea in detail and show how — through character
theory — the relevant computations reduce to computing the eigenvalues of an operator on each
irreducible representation. In Section 4, we return to consideration of the likelihood function
for elapsed time under a stochastic model and discuss possible future directions to efficiently
evaluate this function using approximate, numeric means based solely on eigenvalue algorithms.
We close with a discussion of our plans for future explorations of the key ideas presented in
this paper.

2

Rearrangement models for bacterial genomes

There is a diverse range of evolutionary events that cause rearrangements of regions on the
bacterial genome. The rearrangements that are responsible for the majority of large-scale
changes include inversion, translocation, duplication, insertion (through lateral gene transfer),
and deletion (for example, see [25] or the reviews [19, 38]). The techniques used in this paper
are suited to those that yield rearrangement models based on permutation groups [10, 17], in
which each event is invertible (able to be undone). These include inversions (the most frequent
rearrangement event), as well as translocation. To simplify the discussion, we will focus on
inversions. An inversion refers to the excision of a segment of DNA, followed by its reinsertion
in the same place but with the opposite orientation. These events are commonly performed by
the action of “site-specific recombinases”. These enzymes act by cutting two strands that cross
(at specific short sequences that the enzyme recognizes), and rejoining them to each other.
A common modelling simplification for these processes is to envisage the circular bacterial
genome as consisting of a cyclic ordering of N genes (or, more generally, any identifiable set
of homologous3 regions of DNA). We depict the canonical cyclic ordering in Figure 1. In
this picture, each individual gene/region can be considered as either (i) orientated, or (ii)
unorientated (with directionality defined by the chemical orientation of a single strand of nucleic
acid). For the unoriented case the relevant group is the symmetric group SN , while in the
oriented case the group is the hyperoctahedral group (these are the Coxeter groups of types
A and B respectively, in the context of Lie theory). In this paper we will focus on the the
simpler unorientated formulation because we are primarily interested in addressing the core
computational obstructions inherent in all rearrangement models, and there is no advantage
in further obscuring the discussion with more complicated (albeit, more realistic) models. In
any case, the algebraic ideas we present are easily adaptable to any group-based rearrangement
model one cares to use. See [10, 17] for descriptions of a range of possible models amenable to
algebraic approaches.
Given a genome as an (unorientated) cyclic order of N genes, the effect of a rearrangement
is described as the application of a permutation α ∈ SN understood as moving the gene in
location i to location α(i) (this is the “positions” paradigm as described in [3]). Taking the
canonical cyclic ordering as reference, each genome is determined by a permutation σ ∈ SN
where j = σ(i) is understood as indicating that gene i is in location j = σ(i). We can then realise
3

Evolutionary related.

4

10

1

9

10
2

8

5
3

7

2

8

4
6

6

3
7

5

1
4

9

Figure 1: Two circular genomes, each represented by an arrangement of N = 10 regions. The
“canonical” cyclic ordering is on the left.

a sequence of rearrangements in mathematical terms as a composition (in the symmetric group4 )
of a sequence of permutations. That is, the rearrangement α converts the genome determined
by σ to the genome determined by α ◦ σ ≡ ασ.
Without a fixed frame of reference for gene locations, we must treat as equivalent any two
genomes that can be identified after rotation and/or reflection. Given we only observe the
relative location of each gene in a genome, it is appropriate to describe these equivalences
using dihedral permutations DN , that is, the group of symmetries generated by reflection and
cyclic rotation on N elements. To implement this, we declare two permutations σ, σ ′ ∈ SN
to determine equivalent genomes if there is a permutation d ∈ DN such that σ ′ = dσ. This is
reasonable, since k = d(j) = dσ(i) = σ ′(i) is understood as indicting gene i is in location k = d(j).
When a frame of reference is not fixed, this means we may identify each genome with a right
coset DN σ, which is consistent with the orbit-stabilizer theorem in identifying the number of
distinct genomes as |SN /DN | = N!/2N.
Defining a biological model of rearrangements M amounts to fixing a subset M ⊂ SN , most
commonly consisting entirely of involutions (self inverse permutations5 ). For example, perhaps
the most simple yet biologically plausible model emerges by restricting to inversion of adjacent
regions only, so
M = {(1, 2), (2, 3), (3, 4), . . . , (N − 1, N), (N, 1)}.

(1)

On this point however we keep the discussion general and suppose the biological model consists
of R involutions:
M = {s1 , s2 , . . . , sR } ⊂ SN
so s2i = e for each i = 1, 2, . . . , R.
Given a genome G0 and a fixed frame of reference, without loss of generality we may label
the regions of G0 canonically such that G0 ≡ e, and the biological model is implemented by
composing members of M with e. For instance the sequence of rearrangements si1 , si2 , . . . , sik
produces the new genome G ≡ σ where
σ = s ik . . . s i2 s i1 e = s ik . . . s i2 s i1 .
On top of the chosen biological model M, we choose a simple stochastic model of genome
evolution where rearrangements occur randomly in time with:
(i) Events governed by a Poisson distribution with fixed rate r;
(ii) Each event is a rearrangement chosen uniformly at random from M.
4

We compose permutations from right to left, so (1, 2) ◦ (2, 3) ≡ (1, 2)(2, 3) = (1, 2, 3), for instance.
All inversions are described by involutions, however the converse is false. For example, the permutation
(1, 2)(3, 4) is an involution but does not describe an inversion.
5

5

We note that when M is closed under inverses (which is the case if M consists entirely of
involutions), we have a reversible stochastic process.
As is quite normal in phylogenetic modelling, we cannot infer “molecular rate” r independently from time elapsed T without further calibration [27] (from, for example, the fossil
record). Thus, for simplicity we set r = 1 throughout and interpret T in units of mean number
of events.
Given this setting, the primary computational obstruction we explore — and ultimately
bypass — is the calculation of the number of ways of rearranging one genome into another
using k rearrangements taken from the fixed, but arbitrary, biological model M. For each
σ ∈ SN , this amounts to the calculation of the number βk (σ) of length k words with letters sij
taken from the alphabet M such that
σ = sik sik−1 . . . si1 .
From the biological point of view, when considering two present-day genomes, we rely on
homology to identify common regions and consider one as the canonically ordered reference
G0 ≡ e and the other as the target G ≡ σ. Suppose, further, we would like to obtain the
evolutionary distance between these two genomes. As discussed above, since we are assuming
the biological model M consists entirely of involutions, our stochastic model is reversible and
we may, without loss of generality, consider G to be the descendant of G0 . However, under this
assumption, we cannot be certain how to orientate G with respect to G0 under the dihedral
symmetries DN mentioned above. These observations imply we should take
X
βk (dσ),
αk (G0 → G) ≡ αk (σ) :=
(2)
d∈DN

as the number of passages from genome G0 to the descendant G (or equivalents) under k
rearrangements si ∈ M.
The likelihood for the passage of evolutionary time T is defined as the conditional probability
L(T |G0 → G) = prob(G0 → G|T ) =

X αk (G0 → G) T k e−T
,
k
R
k!
k≥0

(3)

where, consistent with our assumption that each rearrangement is equally likely, Rk is recognized as the total number of passages of length k. One would then like to find the optimum
of this function and hence obtain the maximum likelihood estimate Tb of elapsed time. As we
will see, the key computational obstruction to this ideal is the efficient calculation of the counts
βk (σ). In [32] we found that the iterative algorithm described in [7] can be used to effectively
compute (in a matter of hours) counts for genomes with a maximum of N = 9 regions.

3

Computation of path counts using group characters

In this section we show how to compute the counts βk (σ) using the irreducible representations
of the symmetric group, ultimately showing the computation reduces to an eigenvalue problem.
Consider the symmetric group SN and a rearrangement model M = {s1 , s2 , . . . , sR } ⊂ SN .
Given a target genome G ≡ σ, the basic data we need to compute
βk (σ), the number of words
Pis
R
of length k equal to σ using the letters si ∈ M. We define s := i=1 si , where the sum is taken
in the group algebra CSN (formal linear combinations of group elements over the complex field)
and observe that
X
X
βk (σ)σ.
sk = (s1 + s2 + . . . + sR )k =
(words of length k in the si ) :=
(4)
σ∈SN

6

We take χ : SN → C as the group character of the associated regular representation θ :
SN → GL(N!, C) (obtained by letting SN act on itself on the left). Since there are no fixed
points for this action, χ(σ) = N! if σ = e and is zero otherwise. Thus we have (extending χ
linearly to all of the group algebra):
βk (σ) =

1
1
tr(θ(σ −1 sk )) =
χ(σ −1 sk ).
N!
N!

(5)

Immediately it follows that, for each f ∈ SN satisfying f −1 sf = s, we have equality of path
counts βk (σ) = βk (f −1 σf ) since (using the cyclic property of the trace):




χ((f −1 σf )−1 sk ) = χ f σ −1 f −1 sk = χ σ −1 f −1 sk f = χ σ −1 (f −1 sf )k = χ σ −1 sk .

Consistent with [32, Prop 5.1], in the specific case of the biological model (1) this equality of
path counts holds for all f ∈ DN , and, in general, the equality holds for all f in the stabilizer
of s under the conjugation action s 7→ f −1 sf .
Now suppose θ = ⊕µ dµ θµ is the decomposition of θ into irreducible representations θµ :
SN → GL(dµ , C) indexed by integer partitions µ of N, with each irreducible appearing with
multiplicity equal to its dimension dµ [28]. This reduces the calculation of βk (σ) to computing,
for each integer partition µ, the trace of the matrix θµ (σ −1 sk ) = θµ (σ −1 )θµ (s)k . That is,
βk (σ) =


1 X
1 X
dµ χµ (σ −1 sk ) =
dµ tr θµ (σ)−1 θµ (s)k ,
N! µ
N! µ

(6)

where χµ is the irreducible character corresponding to the irreducible representation labelled
by µ.
Although this reduces the computation of βk (σ) to the calculation of some matrix products
in each irreducible representation of the symmetric group, this procedure remains unfeasible
for N of reasonable size, since the dimension of the irreducible representations can still be very
large. (For example, when N = 20 — which is hardly a biologically impressive case — the
maximum dimension over all irreducible representations is ∼ 109 .)
We further our analysis by noticing that if the biological model M consists entirely of involutions, that is s2i = e, then, in the regular
PN representation, each θ(si ) is a symmetric permutation
matrix. This in turn implies θ(s) = i θ(si ) is symmetric and hence diagonalizable (with real
eigenvalues) and hence each irreducible representation θi (s) is also diagonalizable.
Recalling that a linear operator is diagonalizable if and only if its minimal polynomial q is
a product of distinct linear factors, the minimal polynomial of s must have the form
q = (x − λ1 )(x − λ2 ) . . . (x − λD )
for distinct eigenvalues λi ∈ R. The degree D of the minimal polynomial of θ(s) will of course
depend on the particular model chosen but we can at least say (via the Cayley-Hamilton
theorem) that it is less than N!, since |SN | = N! is the dimension of the regular representation.
For each irreducible θµ it is also the case that q(θµ (s)) = 0, which in turn implies that the
minimal polynomial qµ of each irreducible θµ (s) must divide q. We again conclude (via the
Cayley-Hamilton theorem) that Dµ := deg(qµ ) ≤ dµ .
(µ)
(µ)
(µ)
For each µ, we label the distinct eigenvalues of θµ (s) as λ1 , λ2 , . . . , λDµ ∈ R and define,
(µ)

for each a = 1, 2 . . . , Dµ , the projection operators Ea onto the eigenspaces corresponding to
(µ)
eigenvalue λa via
Y θµ (s) − λ(µ)
b
Ea(µ) =
.
(µ)
(µ)
λ
−
λ
a
b6=a
b
7

(4)
(31)
(22 )
(212)
(14 )

[e]
1
3
2
3
1

[(1,2)]
1
1
0
-1
-1

[(1,2)(3,4)] [(1,2,3)]
1
1
-1
0
2
-1
-1
0
1
1

[(1,2,3,4)]
1
-1
0
1
-1

Table 1: The character table of S4 . The rows are labeled by the irreducible representations, the
columns by conjugacy classes. Note that (4) and (14 ) label the trivial and sign representations,
respectively.
(µ)

By construction, for each irreducible µ, the Ea provide a set of orthogonal projector operators:
Dµ
X

(µ)

Ea(µ) = 1,

Ea(µ) Eb

= δab 1,

a=1

where δab is the Kronecker delta (equal to 1 if i = j and 0 otherwise).
These observations allow us to write, for each µ:
k

θµ (s) =

Dµ
X

(λa(µ) )k Ea(µ)

(7)

a=1

and hence
−1

tr θµ (σ) θµ (s)

k



=

Dµ
X

(λa(µ) )k tr

θµ (σ)

−1

a=1

Ea(µ)



=

Dµ
X
a=1


(λa(µ) )k tr Ea(µ) θµ (σ)−1 Ea(µ) ,

where in the last equality we use the cyclic property of the trace and orthogonality of the
(µ)
projection operators Ea .
Recall that each irreducible representation θµ can be associated with the action of SN on an
irreducible SN -module V µ ∼
= CDµ [28]. Our results thus far reduce the computation of βk (σ) to
(µ)

(i) Computation of the eigenvalues {λa } of s on each irreducible module V µ under this
action;
(µ)

(µ)

(ii) For each eigenvalue λa and eigenspace Eig(s; µ, a) = Ea V µ , computation of the eigen(µ)
(µ)
values of θµ (σ)−1 restricted to Eig(s; µ, a) under θµ (σ)−1 7→ Ea θµ (σ)−1 Ea .
As a simple first example, we use these ideas to compute exact, closed forms of αk (σ)
for N = 4 under the biological model M = {(1, 2), (2, 3), (3, 4), (1, 4)}. In this case, we have
s = (1, 2) + (2, 3) + (3, 4) + (1, 4) ∈ S4 and the regular representation θ(s) can be computed by
hand without too much trouble as an explicit 24 × 24 matrix. We find that θ(s) has minimal
polynomial
q = x(x − 2)(x + 2)(x − 4)(x + 4).
The irreducible representations of S4 are labelled by the integer partitions {4, 31, 22, 212, 14 },
and the character table is given in Table 1.
The corresponding minimal polynomials are given by
q4 = (x − 4),

q31 = x(x − 2),

q22 = (x − 2)(x + 2),

8

q212 = x(x + 2),

q14 = (x + 4).

(The answers for the trivial and sign representations are immediate given sgn(si ) = −1 for each
si ∈ M.) Using the general formula (7) we have, for k > 0:
θ4 (sk ) = 4k ;

θ14 (sk ) = (−4)k ;

θ22 (sk ) = 2k 41 θ22 (s + 2) + (−2)k (− 41 )θ22 (s − 2);

θ31 (sk ) = 0k (− 12 )θ31 (s − 2) + 2k 21 θ31 (s);

θ212 (sk ) = 0k ( 12 )θ212 (s + 2) + (−2)k (− 21 )θ212 (s).

After simplifying and with reference to the character table, for k > 0 odd we obtain the
(constant time in k) expressions
χ(σ −1 sk ) = 4k (1 − sgn(σ)) + 3 · 2k−1 χ31 (σ −1 s) + 2 · 2k−1 χ22 (σ −1 s) + 3 · 2k−1χ212 (σ −1 s),
where we have used the decomposition of the regular representation into irreducibles
χ = χ4 + 3χ31 + 2χ22 + 3χ212 + χ14 .
For k > 0 even, we similarly obtain
χ(σ −1 sk ) = 4k (1 + sgn(σ)) + 3 · 2k−1χ31 (σ −1 s) + 2 · 2k χ22 (σ −1 ) − 3 · 2k−1χ212 (σ −1 s).
Dividing each of these by 4! and summing over the dihedral group as in (2), we find explicit
expressions (constant time in k) for the number of passages G0 → G ≡ σ using rearrangement
moves belonging to M = {(1, 2), (2, 3), (3, 4), (4, 1)}. In particular we see that we have the
limits (as a proportion of total number of passage counts):
 1
αk (σ)
if sgn(σ) = ±1 and k is even/odd, respectively,
12
→
0, otherwise;
4k
simply reflecting that this particular model M consists entirely of odd involutions.

4

Calculation of likelihood function

Here we use the results of the previous section to further our goal of finding an efficient (possibly approximate) method for evaluating the likelihood function (3). Firstly, we continue
our description using the group algebra (and regular representation) to identity our model as a
continuous-time Markov chain with associated rate matrices arising from a ‘group-based’ model
(as coined in the phylogenetics literature [31]). This discussion furthers our understanding of
the model and guides the way to significant generalization in which different rearrangements
may occur at different (Poisson distributed) rates. Secondly we discuss the decomposition of
the likelihood function (as a function on the symmetric group) into contributions from each
irreducible part.
As discussed in Section 2, we suppose mutation events are Poisson distributed in time with
rate r and are chosen uniformly at random from a model M, consisting of R = |M| allowed
rearrangements (as above, we again set r = 1 and interpret T in units of mean number of
events). Interpreting genome G ≡ σ as the descendant G0 ≡ e, using (2) and (5) we can express
the likelihood function (3) as

X αk (σ) T k e−T
1 X X χ σ −1 d−1 sk T k e−T
L(T |σ) ≡ L(T |G0 → G) =
=
Rk
k!
N! k≥0 d∈D
Rk
k!
k≥0
N

1 X X χ σ −1 dsk T k e−T
.
=
N! k≥0 d∈D
Rk
k!
N

9

Rescaling time as T → T /R, we find
L(T |σ) = e−RT

Tk
1 X X
χ(σ −1 dsk )
N! k≥0 d∈D
k!
N

X Tk
X
−1
−RT 1
χ σ d
=e
sk
N! d∈D
k!
k≥0
N

1 X
χ σ −1 desT
= e−RT
N! d∈D
N

1 X
=
χ σ −1 de(s−R)T ,
N! d∈D

!

N

where R ≡ Re ∈ CSN is understood as a scaling of the identity element in the group algebra.
Since each si ∈ M occurs as a permutation matrix in the regular representation, we observe
that Q := θ(s − Re) = θ(s1 + s2 + . . . + sN − Re) is a “rate matrix” (zero-column sums)
associated to a continuous-time Markov chain. Thus our likelihood function can be written in
an elegant form as the character of an analytic function on group algebra:
L(T |G0 → G) =

1 X
χ(dσ −1 eQT ).
N! d∈D
N

From this point of view eQT is the homogeneous, time-dependent probability transition
matrix whose ij entry gives the probability of a transition from the j th group element to the
ith group element in time T . This shows that our model is a special case of a so-called “groupbased” model [31] (usually formulated for abelian groups only6 ), where a more general, rate
matrix would naturally be given as
X
γf θ(f ),
Q = −b
γ1 +
f ∈SN

P
with b
γ = f ∈SN γf for some parameters γf ≥ 0 giving the rate of transition of group elements
σ → f −1 σ.
From this point of view, one can think of our stochastic model as following from the biological
model M together with the special choice γf = 1 when f ∈ M and γf = 0 otherwise. Treating
each γf as an independent modelling parameter is a nice generalization, but we do not explore
it further beyond making the point it seems unlikely one would be able infer values of the rates
γf from observed data given only two genomes and their relative gene arrangement as input
data.
Returning to the likelihood function, we see that
X
L(T |σ) =
dµ Lµ (T |σ),
µ

where each Lµ (T |σ) is a contribution to the likelihood function corresponding to each irreducible
representation of the symmetric group:
Lµ (T |σ) :=

1 X
χµ (σ −1 deQT ).
N! d∈D
N

6

See [34] for construction of a group-based model using the regular representation of any finite group.

10

In turn, each summand in this expression can be computed, using the orthogonal projection
(µ)
operators Ea , as
Dµ
X
(µ)
−1
QT
−RT
eλa T tr(θ(σ −1 d)Eaµ ).
χµ (σ de ) = e
a=1

This illustrates that, in general terms, each eigenvalue λ of the operator s contributes
to the likelihood function via the substitution λ → eλT . Hence, by far the most significant
contributions to the likelihood function are going to come from the largest eigenvalues. We hope
that this observation will lead to development of an efficient, approximate method for computing
the maximum likelihood estimate of time elapsed under general genome rearrangement models.

5

Discussion

In this paper we have explored the application of the representation theory of finite groups to the
problem of computing evolutionary distance between bacterial genomes. In Section 2 we showed
how to construct genome rearrangement models in the context of permutation groups and gave
the likelihood function for time elapsed under a simple stochastic model based on Poisson
distributed events (3). In Section 3 we reformulated the problem of counting the number
of passages from one genome to another (under a selected set of possible rearrangements)
into the calculation of coefficients in the group algebra of the symmetric group (4). We then
interpreted these coefficients through a trace calculation on each irreducible representation of
the symmetric group (6). Finally, in Section 4 we showed that, if one is primarily interested in
maximum likelihood estimates of elapsed time, it is possible to sidestep computation of these
counts and evaluate the likelihood function using the eigenvalues of the matrices representing
the group algebra element s under each irreducible representation.
Our future directions will focus on this final observation, particular in the context of numerical approximations. For instance, we intend to couple combinatorial constructions of the
irreducible representations of the symmetric group (such as the Specht modules [28]) with standard eigenvalue algorithms (such as Lanczos iteration [35]). The outcomes of these explorations
will be the subject of future work.

References
[1] D. Bader, B. Moret, and M. Yan. A linear-time algorithm for computing inversion distance between signed permutations with an experimental study. Journal of Computational
Biology, 8(5):483–491, 2001.
[2] D. Barry and J. A. Hartigan. Asynchronous distance between homologous DNA sequences.
Biometrics, 43:261–276, 1987.
[3] Sangeeta Bhatia, Pedro Feijão, and Andrew R Francis. Position and content paradigms in
genome rearrangements: the wild and crazy world of permutations in genomics. . arXiv,
pages 1–16, October 2016.
[4] D. Bryant, N. Galtier, and M.-A. Poursat. Likelihood calculation in molecular phylogenetics. In Olivier Gascuel, editor, Mathematics of Evolution and Phylogenetics, pages 33–62.
Oxford University Press, 2005.
[5] David Bryant and Vincent Moulton. Neighbor-Net: An agglomerative method for the
construction of phylogenetic networks. Mol. Biol. Evol., 21:255–265, 2004.
11

[6] K. P. Burnham and D. Anderson. Model Selection and Multi-Model Inference. SpringerVerlag, London, 2002.
[7] Chad Clark, Attila Egri-Nagy, A. R. Francis, and Volker Gebhardt. Bacterial phylogeny
in the Cayley graph. preprint arXiv:1601.04398, 2016.
[8] Aaron E Darling, István Miklós, and Mark A Ragan. Dynamics of genome rearrangement
in bacterial populations. PLoS Genet, 4(7):e1000128, 2008.
[9] AJ Drummond and A Rambaut. Beast: Bayesian evolutionary analysis by sampling trees.
BMC Evol. Biol., 7:214, 2007.
[10] Attila Egri-Nagy, Volker Gebhardt, Mark M Tanaka, and Andrew R Francis. Grouptheoretic models of the inversion process in bacterial genomes. J. Math. Biol., 69(1):243–
265, 2014.
[11] Niklas Eriksen and Axel Hultman. Estimating the expected reversal distance after a fixed
number of reversals. Advances in Applied Mathematics, 32(3):439–453, 2004.
[12] James S Farris. Methods for computing wagner trees. Systematic Biology, 19(1):83–92,
1970.
[13] J. Felsenstein. Cases in which parsimony or compatibility methods will be positively
misleading. Syst. Zool., 27:401–410, 1978.
[14] J. Felsenstein. Evolutionary trees from DNA sequences: a maximum likelihood approach.
J. Mol. Evol., 17:368–376, 1981.
[15] J. Felsenstein. Inferring Phylogenies. Sinauer Associates, Sunderland, 2004.
[16] Walter M Fitch. Toward defining the course of evolution: minimum change for a specific
tree topology. Systematic Biology, 20(4):406–416, 1971.
[17] Andrew R. Francis. An algebraic view of bacterial genome evolution. Journal of Mathematical Biology, 69(6):1693–1718, 2014.
[18] O. Gascuel, editor. Mathematics of Evolution and Phylogenetics. Oxford University Press,
2005.
[19] Nigel D F Grindley, Katrine L Whiteson, and Phoebe A Rice. Mechanisms of site-specific
recombination. Annual review of biochemistry, 75(1):567–605, 2006.
[20] Stéphane Guindon, Jean-François Dufayard, Vincent Lefort, Maria Anisimova, Wim
Hordijk, and Olivier Gascuel. New algorithms and methods to estimate maximumlikelihood phylogenies: assessing the performance of PhyML 3.0. Systematic Biology,
59(3):307–321, 2010.
[21] Sridhar Hannenhalli and Pavel A Pevzner. Transforming cabbage into turnip: polynomial
algorithm for sorting signed permutations by reversals. Journal of the ACM (JACM),
46(1):1–27, 1999.
[22] D. H. Huson and D. Bryant. Application of phylogenetic networks in evolutionary studies.
Mol. Biol. Evol., 23:254–267, 2005.
[23] T. H. Jukes and C. R. Cantor. Mammalian protein metabolism., chapter Evolution of
Protein Molecules, pages 21–132. New York: Academic Press, 1969.
12

[24] M. Kimura. A simple method for estimating evolutionary rates of base substitutions
through comparative studies of nucleotide sequences. J. Mol. Evol., 17:111–120, 1980.
[25] A Klippel, R Kanaar, R Kahmann, and N R Cozzarelli. Analysis of strand exchange and
DNA binding of enhancer-independent Gin recombinase mutants. The EMBO journal,
12(3):1047–1057, March 1993.
[26] P. J. Lockhart, M. A. Steel, M. D. Hendy, and D. Penny. Recovering evolutionary trees
under a more realistic model of sequence evolution. Mol. Biol. Evol., 11:605–612, 1994.
[27] D. Penny. Relativity for molecular clocks. Nature, 436:183–184, 2005.
[28] B. E. Sagan. The Symmetric Group: Representations, Combinatorial Algorithms, and
Symmetric Functions. Second Edition. Graduate Texts in Mathematics. Springer, New
York, 2001.
[29] N. Saitou and M. Nei. The neighbor-joining method for reconstructing phylogenetic trees.
Mol. Biol. Evol., 4(4):406–425, 1987.
[30] David Sankoff and Mathieu Blanchette. The median problem for breakpoints in comparative genomics. In International Computing and Combinatorics Conference, pages 251–263.
Springer, 1997.
[31] C. Semple and M. Steel. Phylogenetics. Oxford Press, 2003.
[32] Stuart Serdoz, Attila Egri-Nagy, J. G. Sumner, Barbara Holland, P. D. Jarvis, Mark M.
Tanaka, and A. R. Francis. Maximum likelihood estimates of rearrangement distances.
eprint arXiv:1602.03962, 2016.
[33] Peter HA Sneath, Robert R Sokal, et al. Numerical taxonomy. The principles and practice
of numerical classification. 1973.
[34] J. G. Sumner, Jesús Fernández-Sánchez, and P. D. Jarvis. Lie Markov models. J. Theor.
Biol., 298:16–31, 2012.
[35] L.N. Trefethen and D. Bau. Numerical Linear Algebra. Society for Industrial and Applied
Mathematics, 1997.
[36] Li-San Wang, Tandy Warnow, Bernard ME Moret, Robert K Jansen, and Linda A Raubeson. Distance-based genome rearrangement phylogeny. Journal of Molecular Evolution,
63(4):473–483, 2006.
[37] RA Welch, V Burland, GIII Plunkett, P Redford, P Roesch, D Rasko, EL Buckles, SR Liou,
A Boutin, J Hackett, et al. Extensive mosaic structure revealed by the complete genome
sequence of uropathogenic Escherichia coli. Proceedings of the National Academy of Sciences, 99(26):17020–17024, 2002.
[38] Wei Yang. Topoisomerases and site-specific recombinases: similarities in structure and
mechanism. Critical Reviews in Biochemistry and Molecular Biology, 45(6):520–534, July
2010.
[39] Ziheng Yang and Bruce Rannala. Bayesian phylogenetic inference using dna sequences: a
markov chain monte carlo method. Molecular Biology and Evolution, 14(7):717–724, 1997.

13

