Embeddability of Kimura 3ST Markov matrices

TL;DR

This paper characterizes when Kimura 3ST Markov matrices can be embedded, computes their volume within all Markov matrices, and discusses mutation rate identifiability and symmetry issues.

Contribution

It provides a eigenvalue-based characterization of embeddability and analyzes mutation rate identifiability in Kimura 3ST models.

Findings

Eigenvalue conditions determine embeddability.

Computed volume of embeddable matrices within all Markov matrices.

Mutation rates are not always identifiable from substitution probabilities.

Abstract

In this note, we characterize the embeddability of generic Kimura 3ST Markov matrices in terms of their eigenvalues. As a consequence, we are able to compute the volume of such matrices relative to the volume of all Markov matrices within the model. We also provide examples showing that, in general, mutation rates are not identifiable from substitution probabilities. These examples also illustrate that symmetries between mutation probabilities do not necessarily arise from symmetries between the corresponding mutation rates.

Tables

Table 1. Table 1: Volumes and relative volumes of the embeddable K3 matrices. The relative volumes of embeddable matrices within each spaces are shown in the second row of the table and are obtained as the quotients V ( Δ ε ∩ ⋅ ) / V ( Δ ⋅ ) V(\Delta_{\varepsilon}\cap\cdot)/V(\Delta_{\cdot}) .

	${\Delta }_{\epsilon }$	$\Delta$	${\Delta }_{\ast }$	${\Delta }_{+}$	${\Delta }_d$
$V(\cdot )$	1/4	8/3	2/3	1/2	1/3
relative vol. of embeddable	1	3/32	3/8	1/2	0.61008

Equations

$$e^X=\sum _{n\ge 0}\frac{X^n}{n!}.$$

$$\underset{n\to \infty }{\lim }{e}^{Q_n}=e^{{\lim }_{n\to \infty }{Q}_n}=e^{(0)}=Id,$$

$$M=\left(\begin{array}{cccc}a& b& c& d\\ b& a& d& c\\ c& d& a& b\\ d& c& b& a\end{array}\right).$$

$$\displaystyle S:=\left(\begin{array}[]{rrrr}1&1&1&1\\ 1&1&-1&-1\\ 1&-1&1&-1\\ 1&-1&-1&1\end{array}\right).$$

$$x:=a+b-c-d,y:=a-b+c-d,z:=a-b-c+d.$$

$$K(a,b,c,d)\cdot K(a^{\prime },b^{\prime },c^{\prime },d^{\prime })=(SD{S}^{-1})\cdot (S{D}^{\prime }{S}^{-1})=SD{D}^{\prime }{S}^{-1},$$

$$e^{Q_{1}}e^{Q_{2}}=exp\Big{(}Q_{1}+Q_{2}+\frac{1}{2}[Q_{1},Q_{2}]+\dots\Big{)}=exp\big{(}Q_{1}+Q_{2}\big{)},$$

$$x\ge yz,y\ge xz,z\ge xy.$$

$$\beta \ge 0\iff \frac{log(x)-log(y)-log(z)}{4}\ge 0\iff log\left(\frac{x}{yz}\right)\ge 0\iff x\ge yz.$$

$$x\ge yz,y\ge xz,z\ge xy.$$

$$\displaystyle M=\frac{1}{4}\left(\begin{array}[]{cccc}1+e^{-\lambda}-2e^{-\mu}&1+e^{-\lambda}+2e^{-\mu}&1-e^{-\lambda}&1-e^{-\lambda}\\ 1+e^{-\lambda}+2e^{-\mu}&1+e^{-\lambda}-2e^{-\mu}&1-e^{-\lambda}&1-e^{-\lambda}\\ 1-e^{-\lambda}&1-e^{-\lambda}&1+e^{-\lambda}-2e^{-\mu}&1+e^{-\lambda}+2e^{-\mu}\\ 1-e^{-\lambda}&1-e^{-\lambda}&1+e^{-\lambda}+2e^{-\mu}&1+e^{-\lambda}-2e^{-\mu}\end{array}\right)$$

$$\displaystyle Q_{k}=\overline{S}\;\mathrm{diag}\Big{(}0,-\lambda,-\mu+(2k+1)\pi i,-\mu-(2k+1)\pi i\Big{)}\;\overline{S}^{-1},\qquad k\in\mathbb{Z}$$

$$\overline{S}=\left(\begin{array}{cccc}1& 1& 1+i& 1-i\\ 1& 1& -1-i& -1+i\\ 1& -1& 1-i& 1+i\\ 1& -1& -1+i& -1-i\end{array}\right)$$

$$Q_k=\frac{1}{4}\left(\begin{array}{cccc}-\lambda -2\mu & -\lambda +2\mu & \lambda -2\pi (2k+1)& \lambda +2\pi (2k+1)\\ -\lambda +2\mu & -\lambda -2\mu & \lambda +2\pi (2k+1)& \lambda -2\pi (2k+1)\\ \lambda +2\pi (2k+1)& \lambda -2\pi (2k+1)& -\lambda -2\mu & -\lambda +2\mu \\ \lambda -2\pi (2k+1)& \lambda +2\pi (2k+1)& -\lambda +2\mu & -\lambda -2\mu \end{array}\right).$$

$$2\pi \le \lambda \le 2\mu .$$

$$\displaystyle\left(\begin{array}[]{cccc}0.2410701&0.2593858&0.2497720&0.2497720\\ 0.2593858&0.2410701&0.2497720&0.2497720\\ 0.2497720&0.2497720&0.2410701&0.2593858\\ 0.2497720&0.2497720&0.2593858&0.2410701\end{array}\right).$$

$$Q_k$$

$$4\pi \le \lambda \le 2\mu ,$$

$$\displaystyle\footnotesize M=a\;\begin{pmatrix}1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1\\ \end{pmatrix}+b\;\begin{pmatrix}0&1&0&0\\ 1&0&0&0\\ 0&0&0&1\\ 0&0&1&0\\ \end{pmatrix}+c\;\begin{pmatrix}0&0&1&0\\ 0&0&0&1\\ 1&0&0&0\\ 0&1&0&0\\ \end{pmatrix}+d\;\begin{pmatrix}0&0&0&1\\ 0&0&1&0\\ 0&1&0&0\\ 1&0&0&0\\ \end{pmatrix},\quad\begin{array}[]{c}a,b,c,d\geq 0,\\ a+b+c+d=1\end{array}$$

$$M=\frac{1}{4}\left(\begin{array}{cccc}1& 1& 1& 1\\ 1& 1& 1& 1\\ 1& 1& 1& 1\\ 1& 1& 1& 1\end{array}\right).$$

$$\Delta$$

$$V({\Delta }_d)=\frac{1}{6}|det(\overrightarrow{p_1{q}_2},\overrightarrow{p_1{q}_3},\overrightarrow{p_1{q}_4})|=1/3.$$

$$\frac{1}{6}|det(\overrightarrow{O{q}_i},\overrightarrow{O{q}_j},\overrightarrow{O{p}_{1ij}})|=1/18.$$

$$V({\Delta }_{\mathcal{E}})={\int }_0^1\left({\int }_0^z{\int }_{zx}^{x/z}dydx+{\int }_z^1{\int }_{xz}^{z/x}dydx\right)dz$$

$${\Delta }^{JC}$$

$${\Delta }_{\ast }^{JC}$$

$${\Delta }_{+}^{JC}={\Delta }_{\epsilon }^{JC}$$

$${\Delta }_d^{JC}$$

Full text

Keywords: Markov matrix; Markov generator; eigenvalues; evolutionary model; embeddability

1 Introduction

Genomic data expressed by means of sequence alignments is widely used to infer phylogenetic relationships between species. Substitution models are used to describe the evolutionary process that leads from one DNA sequence to another. These models are usually given in terms of a family of Markov matrices with a prescribed structure. The entries of these matrices represent the conditional probabilities of nucleotide substitution between one sequence and the other, and can be obtained either by counting the relative frequencies of these substitutions or fitting the parameters of the model using maximum likelihood. Usually, the structure imposed by the model is motivated by some biological / biochemical properties observed (e.g. the Kimura 3ST model [Kim81]) or some computational / mathematical convenient assumptions to deal with the model (e.g the GTR model [Tav86] or Lie Markov models [SFSJ12]). Moreover, evolution is usually modelled by means of Markov chains, together with the additional assumption that all sites in the sequences evolve independently and according to the same probabilities.

A general approach in modelling evolution corresponds to regarding time as a continuous variable where substitution events always happen at the same rate, which remains constant throughout the whole evolutionary process. This leads to the homogeneous continuous-time substitution models, where only Markov matrices that are the exponential of a rate matrix are considered. Clearly this is used as an approximation to biological reality where it is well known that transition rates vary over time [HPCD05, HSP*+*07] and also among the different branches of the phylogenetic tree [LSB*+*98]. However, given the bias / variance compensation of the statistical analysis [BA02], modelling phylogenetic evolution as a non-homogeneous process is not statistically feasible in practice (cf. [SFSJ12]).

A different approach appears when one regards the evolutionary process as a whole and only takes into account the conditional probabilities between the original and the final sequences, without caring about rates of mutation. When these probabilities are taken as the parameters of the model, we deal with the so-called algebraic models111Here, “algebraic” refers to the fact that the probabilities of pattern observation at the leaves of a phylogenetic tree evolving under these models are given by algebraic expressions (only sums and products) in terms of the parameters of the model.. Algebraic models have been used in a number of theoretical papers, including [AR08, SS05, DK08, CFS10].

If one attempts to connect both approaches, a natural question is to decide whether a given Markov matrix is the exponential of some rate matrix, whose entries would be some kind of average of the rates involved throughout the evolutionary process. In this case, we say the matrix is embeddable and this question is known in the literature as the embedding problem for Markov matrices. An easier version of this problem is to decide whether the rate matrices associated to the embeddable matrices of a particular (algebraic) model $\mathcal{M}$ should keep the same symmetries as the model ($\mathcal{M}$-embeddability, see definition in Section 2.2). The embedding problem is relevant even if restricted to continuous-time models since it is not true in general that the product of embeddable matrices is necessarily embeddable (indeed, the Baker-Campbell-Hausdorff formula [Cam97] leads to ask whether some series of matrices is convergent or not, which is not always true [BC04]). These questions are closely related to the problem of the multiplicative closure of continuous-time models, namely whether the product of matrices $e^{Q_{1}}e^{Q_{2}}$ where $Q_{1}$ and $Q_{2}$ are rate matrices in one particular (continuous-time) model can be obtained as some $e^{Q}$ for some rate matrix $Q$ in the same model. After [SFSJ12, Sum17], it is known that there are popular models which are not multiplicatively closed, notably including the GTR model and the HKY model.

The reader is referred to [Dav10] for a nice overview of the embedding problem from a mathematical point of view. In a more biological and applied setting, the paper by Verbyla et al. [VYP*+*13] deals with the possible consequences for phylogenetic inference. Also, the paper [SJFS*+*12] and the more recent paper [WSL*+*17] deal with the incidental question of how the lack of (multiplicative) closure in substitution models have consequences for the phylogenetic analysis of data.

In this paper, we deal with the embedding problem from a theoretical perspective. The main goal is to obtain a characterization for the embeddability of generic matrices of the Kimura 3ST model [Kim81]. From our results, we will be able to compute the whole volume of embeddable Kimura 3ST matrices and compare it with the volume of the whole space of Kimura 3ST Markov matrices. At the same time, we provide a number of examples showing matrices that are embeddable but for which the mutation rates are not identifiable or do not keep the same structure of the model. The recent paper [KK17] deals with the similar question of characterizing embeddable matrices of symmetric group-based phylogenetic models, but focusing on the existence of rate matrices strictly in the model.

The organization of the paper is as follows. In section 2, we recall some definitions and basic facts concerning the embedding problem and the Kimura 3 parameter model. Here, we also show that any embeddable matrix is biologically relevant since it can be seen as the transition matrix of a concatenation of realistic evolutionary processes (“realistic” here means a process whose transition matrix is close to the identity matrix, see Theorem 2.2). In section 3, we prove the main theorem which characterizes under the (generic) assumption of having different eigenvalues the Kimura 3ST embeddable matrices in terms of inequalities to be satisfied by the eigenvalues. We devote as well some attention to the case of matrices with repeated eigenvalues as they present certain situations that may be interesting from a theoretical and applied point of view. Namely, these matrices show that the identifiability of the mutation rates is not a generic property for the Kimura 2ST model or the Jukes-Cantor model, as well as that there are embeddable matrices with rate matrices that do not keep the same symmetries of the model (see Theorem 3.9). As a consequence of the characterization mentioned above, in section 4 we are able to compute the volume of embeddable matrices and compare it to the volume of all Kimura 3ST Markov matrices. Finally, Section 5 discusses implications and possibilities for future work.

2 Preliminaries

2.1 Embedding problem of Markov matrices

We denote by $M_{k}(\mathbb{K})$ the space of all square $k$-matrices with entries in a field $\mathbb{K}$, where $\mathbb{K}$ is $\mathbb{R}$ or $\mathbb{C}$. Given a matrix $A\in M_{k}(\mathbb{K})$, we say that $B\in M_{k}(\mathbb{K})$ is a logarithm of $A$ if $e^{B}=A$, where the exponential of a matrix is defined as

[TABLE]

A classical result states that $det(e^{X})=e^{tr(X)}$, so the determinant of any matrix of the form $e^{X}$ is never 0. Given a non-negative complex number $x\in\mathbb{C}\setminus\mathbb{R}^{-}$, we will denote by $log(x)$ its principal logarithm, that is, the only logarithm of $x$ that lies in the strip $\{z\mid-\pi<Im(z)<\pi\}$. Although the exponential map of matrices is not injective, it is known that if $A$ is a matrix with no negative eigenvalues, there is a unique logarithm $X$ of $A$ all of whose eigenvalues are given by the principal logarithm of the eigenvalues of $A$ (Theorem 1.31 of [Hig08]). We will refer to this as the principal logarithm of $A$ and we will denote it by $Log(A)$. In the particular case where the matrix $A$ is diagonalizable, $A=S\,D\,S^{-1}$ then $Log(A)=SLog(D)S^{-1}$, where $Log(D)$ is the diagonal matrix with diagonal entries equal to the principal logarithm of the eigenvalues of $A$.

Definition 2.1**.**

A matrix $M\in M_{k}(\mathbb{R})$ is said to be a Markov matrix if all the entries are non-negative and the rows sum to one. A matrix $Q\in M_{k}(\mathbb{R})$ is said to be a rate matrix if all the non-diagonal entries are non-negative and the rows sum to zero.

If $Q$ is a rate matrix, it is well-known that $e^{tQ}=\sum_{n\geq 0}\frac{t^{n}Q^{n}}{n!}$ is a Markov matrix for all $t\geq 0$. That is why rate matrices are also referred as Markov generators [Dav10]. However, not every Markov matrix can be obtained in this way. A Markov matrix $M$ is said to be embeddable if $M=e^{Q}$ for some rate matrix $Q$. The embedding problem attempts to decide which (Markov) matrices are embeddable, that is, which matrices can be written as $M=e^{Q}$, where $Q$ is a rate matrix. We would like to point out that every embeddable matrix can be obtained as the substituion matrix of a long-running biologically realistic Markov process. Namely,

Theorem 2.2**.**

Every embeddable matrix is the product of embeddable matrices close to the identity matrix.

Proof.

Assume that $M$ is an embeddable Markov matrix: $M=e^{Q}$. Clearly, $Q_{n}:=\frac{1}{n}Q$ is still a rate matrix for any $n\geq 1$, so $M=(e^{Q_{n}})^{n}$ appears as the $n$-th power of a Markov matrix. Moreover, since

[TABLE]

we can take $n$ big enough so that $e^{Q_{n}}$ is as close to $Id$ as wanted. ∎

2.2 Kimura models

In this work we deal with the substitution model introduced by Kimura in [Kim81]. The Kimura 3ST model assigns three parameters to different type of substitutions: one parameter for transitions, i.e. substitutions between purines ($\tt A\leftrightarrow\tt G$) or pyrimidines ($\tt C\leftrightarrow\tt T$), and two parameters for transversions, i.e. substitutions that change the type of nucleotide: from purine to pyrimidine or vice versa. Ordering the set of nucleotides as ${\tt A},{\tt G},{\tt C},{\tt T}$, the Markov matrices within the model are described by the following structure:

Definition 2.3**.**

A matrix $M\in M_{4}(\mathbb{C})$ is Kimura 3ST (is K3 or has K3 form, for short) if it has the following structure:

[TABLE]

For ease of reading we will use the notation $M=K(a,b,c,d)$ to denote a matrix with the structure in (1).

When $M$ is a Markov matrix, the structure above describes the symmetry between the substitution probabilities of the Kimura 3ST model. Keeping the order of the nucleotides, if $i,j={\tt A,G,C,T}$, the $(i,j)$-entry corresponds to the probability of nucleotide $i$ being replaced by nucleotide $j$. The submodels of Kimura 3ST model, namely Kimura 2ST [Kim80] and Jukes-Cantor [JC69], appear when more symmetries are considered: if $c=d$, we say that the matrix $M$ is Kimura 2ST (K2, for short); if $b=c=d$, we say that $M$ is Jukes-Cantor (JC, for short).

If one restricts the embedding problem to one particular model $\mathcal{M}$ given by some equalities between the entries of the Markov matrices (such as in Kimura 3ST, Kimura 2ST, Jukes-Cantor), it is natural to ask whether these matrices have Markov generators fulfilling the same equalities. In this case, we say that these matrices are $\cal{M}$-embeddable. Example 3.6 of the next section shows that this is not true in general.

The following lemma is fundamental for our study. It essentially claims that all K3 matrices are diagonalized through the following Hadamard matrix:

[TABLE]

Note that $S^{2}=4\cdot Id$, thus $S^{-1}=\frac{1}{4}S$.

Lemma 2.4**.**

A matrix is K3 if and only if it can be diagonalized through $S$. In this case, $K(a,b,c,d)=S\,D\,S^{-1}$, where $D=\operatorname{diag}(a+b+c+d,\ a+b-c-d,a-b+c-d,\ a-b-c+d)$. In particular, a K3 matrix is real if and only if its eigenvalues are all real numbers.

Proof.

The proof is straightforward and follows by direct computation. ∎

Since the rows of K3 Markov matrices sum to $a+b+c+d=1$, the first eigenvalue must be equal to one. Hence, we derive that every K3 Markov matrix is determined by the set of the other eigenvalues:

[TABLE]

Moreover, it is immediate to check that a matrix $M$ as in (1) is K2 (resp. JC) if and only if $y=z$ (resp. $x=y=z$).

It also follows that

Theorem 2.5**.**

The product of $K3$-embeddable matrices is $K3$-embeddable.

Proof.

By virtue of Lemma 2.4, we have

[TABLE]

which is a K3 matrix. Since $D$ and $D^{\prime}$ are diagonal matrices, we have that $D\,D^{\prime}=D^{\prime}\,D$ and from this, the product of K3 matrices is commutative. Therefore, if $Q_{1}$, $Q_{2}$ are K3 rate matrices, the Baker-Campbell-Hausdorff formula [Cam97] gives

[TABLE]

where $[A,B]=AB-BA$ is the Lie bracket. Since $Q_{1}+Q_{2}$ is also a K3 rate matrix, the claim follows. ∎

3 Embeddability of Kimura Markov matrices

The first result of this section describes how to compute the principal logarithm of a K3 Markov matrix and characterizes when it is a Markov generator.

First, we need a lemma which follows from the definition of the exponential matrix.

Lemma 3.1**.**

Let $Q$ be a logarithm of $M$. If $v$ is an eigenvector of $Q$ with eigenvalue $\mu$, then $v$ is an eigenvector of $M$ with eigenvalue $e^{\mu}$.

Remark 3.2**.**

Note that the converse is not true in general. For instance, any vector of $\mathbb{C}^{2}$ is an eigenvector of M:=\left(\begin{array}[]{cc}-1&0\\ 0&-1\end{array}\right) with eigenvalue $-1$. However, the matrix Q=\left(\begin{array}[]{cc}0&-1\\ 1&0\end{array}\right) is a logarithm of $M$, with eigenvalues $i$ and $-i$ and corresponding eigenspaces $[(i,1)]$ and $[(-i,1)]$. In particular, any vector not in these subspaces cannot be an eigenvector of $Q$.

The following result solves the $K3$-embeddability.

Theorem 3.3**.**

Let $M$ be a K3 Markov matrix with eigenvalues $1,x,y,z$. Then,

1. i)

$M$* has a real logarithm $Q$ with K3 form if and only if $x,y,z>0$. In this case, $Q$ is necessarily the principal logarithm $Log(M)=S\;\mathrm{diag}(0,log(x),log(y),log(z))\;S^{-1}$.* 2. ii)

$M$* is $K3$-embeddable (i.e. $Log(M)$ is a Markov generator) if and only if*

[TABLE]

Proof.

(i) First of all, if $Q$ is a real logarithm of $M$ with K3 form, then the eigenvalues of $Q$ are real (Lemma 2.4) and, by Lemma 3.1, the eigenvalues of $M$ have to be positive. Conversely, if $x,y,z>0$ then we can take principal logarithms of $x,y,z$ and define $Q$ as the principal logarithm of $M$, i.e. $Q:=Log(M)=S\;\mathrm{diag}(0,log(x),log(y),log(z))\;S^{-1}$. Then, $Q$ is a real matrix, $e^{Q}=M$ and it has K3 form because of Lemma 2.4.

We proceed to show that $Log(M)$ is the only possible real logarithm with K3 form. By Lemma 2.4, if $Q^{\prime}$ is any other real logarithm with K3 form, then it can be written as $Q^{\prime}=S\;\mathrm{diag}(r,s,u,v)\;S^{-1}$, for some $r,s,u,v\in\mathbb{R}$. Because of Lemma 3.1, these values are real logarithms of $1,x,y,z$. Necessarily, $Q^{\prime}$ must be the principal logarithm of $M$.

(ii) Using Lemma 2.4 we have $Log(M)=K(\alpha,\beta,\gamma,\delta)$, where \alpha=\frac{1}{4}\big{(}log(x)+log(y)+log(z)\big{)}, \beta=\frac{1}{4}\big{(}log(x)-log(y)-log(z)\big{)}, \gamma=\frac{1}{4}\big{(}log(y)-log(x)-log(z)\big{)} and \delta=\frac{1}{4}\big{(}log(z)-log(x)-log(y)\big{)}. It is immediate that $\alpha+\beta+\gamma+\delta=0$. Therefore, we only need to check that the non-diagonal entries $\beta$, $\gamma$ and $\delta$ of $Log(M)$ are non-negative if and only if the inequalities (4) are satisfied. This is straightforward: for instance,

[TABLE]

The other inequalities are proved similarly. ∎

Lemma 3.4**.**

If $M$ is a K3 matrix with no repeated eigenvalues and $Q$ is a real logarithm of $M$, then $Q$ is the principal logarithm of $M$. In particular, $Q$ has $K3$ form.

Proof.

Because of Lemma 3.1, if the matrix $M$ has no repeated eigenvalues, then so does the matrix $Q$. It follows that $Q$ diagonalizes, and that $M$ and $Q$ have the same eigenvectors. In particular they both diagonalize through the matrix $S$ and hence $Q$ must be K3 (Lemma 2.4). Now, it is enogh to apply Theorem 3.3. ∎

Corollary 3.5**.**

Let $M$ be a K3 Markov matrix with no repeated eigenvalues. Then, the following are equivalent:

1. (i)

$M$* is embeddable;* 2. (ii)

$M$* is $K3$-embeddable;* 3. (iii)

the eigenvalues $x,y,z$ of $M$ are strictly positive, and satisfy

[TABLE]

Proof.

It follows directly from theorem 3.3 and Lemma 3.4. ∎

The case of repeated eigenvalues

After the previous result, it is natural to ask what can be said in the case of repeated eigenvalues. In this case, there are some theoretically interesting examples showing that:

There are embeddable K3 matrices with repeated negative eigenvalues. Remarkably, these matrices do not admit K3 Markov generators and their rates are not identifiable. See forthcoming Example 3.6. 2. 2.

Restricted to the case of repeated positive eigenvalues, there are K3 Markov matrices for which rates are not identifiable.

Although the matrices presented in the forthcoming examples are close to saturation and have a big mutation rate, they have still biological interest by virtue of Theorem 2.2.

Repeated negative eigenvalues

It is well known that a matrix with negative eigenvalues has a real logarithm if and only if the negative eigenvalues have even multiplicity [Cul66]. Consider a K3 Markov matrix (actually, it is a $K2$ matrix) with eigenvalues $1$, $e^{-\lambda}$ and $-e^{-\mu}$ with multiplicity $2$, where $\lambda,\mu\geq 0$:

[TABLE]

By virtue of Theorem 3.3, $M$ has no real logarithm with K3 or K2 form.

According to Theorem 1.28 in [Hig08], distinct logarithms for the matrix $M$ are obtained as $\overline{S}D\overline{S}^{-1}$, where $D$ is a diagonal matrix with distinct determinations of the logarithm of the eigenvalues of $M$ and the columns of $\overline{S}$ form a basis of eigenvectors of $M$ (see also Chap. VIII§8 in [Gan59]). Hence, if we choose a pair of conjugated eigenvectors $v,\overline{v}$ of the eigenvalue $-e^{-\mu}$ as the third and fourth columns of the matrix $\overline{S}$, then all the matrices of the form

[TABLE]

are real logarithms of $M$. For example, we can take

[TABLE]

to obtain

[TABLE]

It is straightforward to check that this is a Markov generator if and only if $2\pi|2k+1|\leq\lambda$ and $\lambda\leq 2\mu$. In particular, we see that $M$ is embeddable if

[TABLE]

Moreover, in this case, it is enough to take $k=0$ and $k=-1$ to obtain a pair of Markov generators for $M$.

We illustrate this construction with a numerical example.

Example 3.6**.**

Let us take $\lambda=7$ and $\mu=4$, so that $2\pi\leq\lambda\leq 2\mu$. Then the matrix $M$ is (rounding off to 7 decimals)

[TABLE]

As mentioned above, $Log(M)$ is not a real matrix and hence is not a rate matrix either. In spite of that we are still able to find a pair of Markov generators for $M$ by taking $k=0$ and $k=-1$ respectively:

$Q_{0}=\frac{1}{4}\begin{pmatrix}-15&1&7-2\pi&7+2\pi\\ 1&-15&7+2\pi&7-2\pi\\ 7+2\pi&7-2\pi&-15&1\\ 7-2\pi&7+2\pi&1&-15\\ \end{pmatrix}\qquad Q_{-1}=\frac{1}{4}\begin{pmatrix}-15&1&7+2\pi&7-2\pi\\ 1&-15&7-2\pi&7+2\pi\\ 7-2\pi&7+2\pi&-15&1\\ 7+2\pi&7-2\pi&1&-15\\ \end{pmatrix}$

The previous example shows that a K3 Markov matrix $M$ can be embeddable, even if it does not have any Markov generator with K3 form (see Theorem 3.3 (i)). Thus, embeddability of $K3$ matrices does not imply $K3$-embeddability, and is not determined by the principal logarithm (cf. [VYP*+*13, KK17]). This fact exhibits that the structure of the K3 model, which imposes certain symmetries between transitions and between transversions, is not always captured by the same symmetries between the mutation rates (cf. [Kim80, Kim81]). The reader may note that the expected number of substitutions for a process ruled by a matrix $M$ as in (9) is $-\frac{1}{4}tr(Q_{k})=\frac{1}{4}(\lambda+2\mu)\geq\pi.$

Remark 3.7**.**

After the previous example, we derive easily that the product of embeddable matrices within the Kimura 3ST model is not necessarily embeddable (cf. Theorem 2.5). Indeed, it is enough to consider the embeddable matrix $M$ shown in (15) and any K3 matrix $N$ with positive eigenvalues $1,x,y,z$ satisfying the inequalities (4) so that $N$ is embeddable. The product of $M$ and $N$ is clearly a K3 Markov matrix, whose eigenvalues are the product of the eigenvalues of $M$ and $N$. Thus, $MN$ has two different negative eigenvalues and another positive eigenvalue. By virtue of [Cul66], $MN$ has no real logarithm so, in particular, it cannot be embeddable.

Repeated positive eigenvalues

A similar computation to that for negative eigenvalues can be done by assuming positive and repeated eigenvalues: $1$, $e^{-\lambda}$ and $e^{-\mu}$ with multiplicity 2. In this case, we can produce real logarithms by taking

[TABLE]

Such a matrix is a Markov generator if and only if $4\pi|k|\leq\lambda$ and $\lambda\leq 2\mu$. In particular, we see that if $0\leq\lambda\leq 2\mu$, then $M$ is $K2$-embeddable since $Q_{0}$ is a Markov generator. If in addition,

[TABLE]

then we have (at least) three Markov generators, which correspond to $k=0,\pm 1$.

Note that the Markov generator $Q_{0}$ is a K3 matrix (corresponding to the principal logarithm) while $Q_{1}$ and $Q_{-1}$ are not. The expected number of substitutions for a process ruled by these matrices is $-\frac{1}{4}tr(Q_{k})=\frac{1}{4}(\lambda+2\mu)\geq 2\pi$.

Remark 3.8**.**

Following the construction of Theorem 2.2, for any $n\geq 1$ the matrix $e^{(1/n)Q_{0}}$ is a K2 Markov matrix that approaches to $Id$ when $n$ grows. If a substitution process is ruled by such a matrix for some long time $t>0$, the rates of the resulting Markov matrix $e^{t(Q_{0}/n)}$ become unidentifiable at some point. This situation cannot occur for K3 embeddable matrices with different eigenvalues (see Theorem 3.3).

The existence of two or more Markov generators for the same Markov matrix (both in the case of negative and positive eigenvalues) exhibit that, in general, mutation rates are not identifiable from the mutation probabilities. Even more, if we restrict to the Kimura 3ST submodels, such matrices do not appear as marginal cases. Indeed, we have seen that identifiability of rates of a K2 embeddable matrix $M$ with eigenvalues $1,e^{-\lambda},e^{-\mu}$ (with multiplicity 2) does not hold in the subspace defined by the inequalities (16), which has positive measure within the space of all K2 Markov matrices. Furthermore, we have also seen that those Markov matrices with a negative repeated eigenvalue $1,e^{-\lambda},-e^{-\mu}$ (with multiplicity 2) satisfying (10) are embeddable but not K2-embeddable (nor K3-embeddable). The space of such matrices has also positive measure within space of K2 Markov matrices.

Similarly, for a Jukes-Cantor matrix with eigenvalues $1$ and $e^{-\lambda}$ (with multiplicity 3), identifiability of rates does not hold in the subspace defined by $4\pi\leq\lambda$, which has positive measure within the space of all JC matrices. On the other hand, every embeddable matrix is $JC$-embeddable since by [Cul66], a necessary condition for a Markov Jukes-Cantor matrix to be embeddable is that its eigenvalues $1,x$ are positive, and then Theorem 3.3 ensures that the principal logarithm (which is a JC matrix) is a Markov generator (it is enough to check that $x\geq x^{2}$, which is true since $x\in[0,1]$).

To conclude this section, we state the following theorem which summarizes the consequences of the previous examples:

Theorem 3.9**.**

For the Kimura 3ST model, $K3$-embeddability and the identifiability of rates are generic properties of embeddable matrices.222A generic property is a property that holds almost everywhere, that is, except for a set of measure zero. 2. 2.

For the Kimura 2ST model, $K2$-embeddability and the identifiability of rates are not generic properties of embeddable matrices. 3. 3.

For the Jukes-Cantor model, every embeddable matrix is JC-embeddable but identifiability of rates is not a generic property of embeddable matrices.

4 The volume of embeddable K3 matrices

Roughly speaking, the goal of this section is to measure how many K3 Markov matrices we are considering when the continuous-time approach is taken and compare this value with the corresponding value of K3 Markov matrices with no further restriction. These values are expressed in terms of the volume of the corresponding subspaces. At the same time, it is direct to obtain the volume of subspaces of K3 Markov matrices with some constraints to make them biologically realistic. To this aim, we proceed to represent K3 Markov matrices in a geometrical way as follows (cf. [CFS08]): keeping the notation introduced in (3), Lemma 2.4 allows us to identify the K3 Markov matrices with the coordinates $(x,y,z)$ of a 3-dimensional space. Moreover, since every K3 matrix is a convex combination of the identity matrix and permutation matrices:

[TABLE]

the space of all K3 Markov matrices describes the 3-dimensional simplex (a regular tetrahedron) with vertices given by the corresponding eigenvalues: $p_{1}=(1,1,1)$, $p_{2}=(1,-1,-1)$, $p_{3}=(-1,1,-1)$ and $p_{4}=(-1,-1,1)$ (see Figure 1). The centroid of this simplex has coordinates (eigenvalues) $O=(0,0,0)$ and corresponds to the matrix

[TABLE]

According to this representation, the Jukes-Cantor matrices [JC69] ($b=c=d$) correspond to the line determined by the identity vertex and the centroid of the simplex ($x=y=z)$, while Kimura 2ST matrices [Kim80] ($c=d$) correspond to a plane section of the simplex ($y=z$).

We proceed to compute the volume of a number of subspaces of K3 Markov matrices that can be interesting from a biological perspective. First of all, we introduce some notation:

[TABLE]

By Lemma 2.4, it is straighforward to check that $\Delta_{d}\subset\Delta_{+}\subset\Delta_{*}\subset\Delta$. Note that the examples of the preceding section show that $\Delta_{\varepsilon}\not\subset\Delta_{*}$ (see the matrix in (15)). However, according to Corollary 3.5, all embeddable matrices that are not in $\Delta_{+}$ correspond to the case of repeated eigenvalues, so they are a marginal case with measure (i.e. volume) 0 within the whole space of K3 matrices. This is because these matrices are constrained by nontrivial algebraic constraints, which make the dimension of the corresponding subspace necessarily smaller. Nevertheless, as shown in (e) and (f) of the next result, there are lots of embeddable matrices with no repeated eigenvalues that are not diagonal dominant.

Theorem 4.1**.**

*We have the following:

(a) $V(\Delta_{d})=1/3$; (b) $V(\Delta_{+})=1/2$; (c) $V(\Delta_{*})=2/3$; (d) $V(\Delta)=8/3$; (e) $V(\Delta_{\mathcal{E}})=1/4$;

(f) $V(\Delta_{\varepsilon}\cap\Delta_{d})\simeq 0.20336$ *

Proof.

(a) The space $\Delta_{d}$ of diagonal dominant matrices is defined by the inequality $a\geq b+c+d$. Since $a+b+c+d=1$, this is equivalent to $a\geq 1/2$. Thus, $\Delta_{d}$ is the regular simplex with vertices $p_{1}$, $q_{2}$, $q_{3}$ and $q_{4}$ (see Figure 1), and its volume is given by the well known formula

[TABLE]

(b) The space $\Delta_{+}$ of matrices with positive eigenvalues is composed of $\Delta_{d}$ together with the simplex with vertices $q_{1}$, $q_{2}$, $q_{3}$ and the centroid $O$. The formula above gives that this last simplex has volume $1/6$. Therefore, $V(\Delta_{+})=V(\Delta_{d})+1/6=1/2$.

(c) The three inequalities defining $\Delta_{*}$ are equivalent to $x+y\geq 0$, $x+z\geq 0$ and $y+z\geq 0$. If we denote by $p_{ijk}$ the centroid of the triangle defined by $p_{i}$, $p_{j}$ and $p_{k}$, it is straighforward to see that $\Delta_{*}$ is composed of $\Delta_{+}$ together with the three simplices defined by $\{O,q_{2},q_{3},p_{123}\}$, $\{O,q_{3},q_{4},p_{134}\}$ and $\{O,q_{4},q_{2},p_{124}\}$. These three simplices have the same volume:

[TABLE]

It follows that $V(\Delta_{*})=1/2+3\,(1/18)=2/3$.

(d) follows similarly to the computation of (a).

(e) The computation of $\Delta_{\varepsilon}$ is more involved. First of all, as noted above, we can restrict the computation to matrices with no repeated eigenvalues. Therefore, it is enough to compute the volume of the space $\mathcal{E}$ defined by the inequalities (4). By cutting the space $\mathcal{E}$ with planes of the form $z=a$ with $a\in[0,1]$, we obtain a shape like Figure 2. The computation follows by integrating the corresponding area for all values of $a$. Given $a\in[0,1]$ and $x\in[0,a]$, the range of values for $y$ is between $a\,x$ and $x/a$, while for $x\in[a,1]$, the value of $y$ lies between $a\,x$ and $a/x$. We are led to compute the following integral

[TABLE]

which can be easily shown to be equal to $1/4$.

(f) The computation of the volume of the intersection $\Delta_{d}\cap\Delta_{\varepsilon}$ is similar to the computation of $V(\Delta_{\varepsilon})$ but for each plane section $z=a$, $a\in[0,1]$ we have to remove the area of the space below the line $x+y=1-a$ (corresponding to non-embeddable matrices). This leads to two different situations: for $a\in[0,3-2\sqrt{2}]$ the line cuts the hyperbola $xy=a$; while for $a\in[3-2\sqrt{2},1]$ line and hyperbola does not meet. The computation of the corresponding integrals is tedious and we do not include it here. The final value has been obtained using the mathematical software SAGE [Dev]. ∎

The values of these volumes illustrate the relative size between the spaces of Markov matrices considered above. Table 1 shows these volumes and the relative volume of embeddable matrices in each of the above subspaces of K3 Markov matrices. These relative volumes are a measure of how many matrices are rejected when taking the continuous-time approach instead of considering other subspaces of matrices. The figures in the second row of the table show, in particular, that there is a big difference between these volumes. For example, embeddable matrices suppose half of the matrices with positive eigenvalues. Similarly, only three out of eight K3 Markov matrices satisfying $a\geq b,c,d$ are embeddable, while we observe that they represent more than the 60% of the diagonal dominant matrices. However, if we only consider diagonal dominant matrices we are rejecting a non-negligible number of embeddable matrices (the difference of volumes between embeddable and diagonal dominant embeddable is $1/4-0.20336=0.046641$, see Theorem 4.1).

Remark 4.2**.**

From Theorem 3.9 (i) it follows that all the values of Table 1 remain the same if we only consider $K3$-embeddable matrices (instead of embeddable matrices).

Since we have no characterization for embeddable matrices with repeated eigenvalues, and this set has positive measure within the Kimura 2ST model (Theorem 3.9), we are not able to compute volumes within this model.

Nevertheless, for the Jukes-Cantor model, Theorem 3.9 states that embeddability is equivalent to JC-embeddability, which has been characterized in Theorem 3.3. Adapting the notation used for the K3 matrices to the JC model, we get that

[TABLE]

A straightforward computation shows that the space of embeddable Jukes-Cantor matrices has volume (length) $\sqrt{3}$ while the space of all Markov Jukes-Cantor matrices has volume $\frac{4}{3}\sqrt{3}$. That is, three out of four Jukes-Cantor matrices are embeddable.

5 Discussion

As suggested in the introduction, there are four connected problems relative to the algebraic and continuous-time evolutionary models. Namely,

whether a given Markov matrix is $e^{Q}$ for some rate matrix $Q$ (embedding problem); 2. 2.

whether a given Markov matrix $M$ within an algebraic model $\mathcal{M}$ is $e^{Q}$ for some rate matrix $Q$ with the same symmetries as $M$ ($\mathcal{M}$-embedding problem); 3. 3.

whether the product of embeddable matrices is embeddable, i.e. if $Q_{1}$ and $Q_{2}$ are rate matrices, then is it true that $e^{Q_{1}}e^{Q_{2}}=e^{Q_{3}}$ for some rate matrix $Q_{3}$? (multiplicative closure of embeddable matrices); 4. 4.

same question as in 3. but with $Q_{1}$,$Q_{2}$ and $Q_{3}$ within a particular continuous-time model (multiplicative closure of $\mathcal{M}$-embeddable matrices).

The first two questions are motivated by the connection between the algebraic and the continuous-time models. The last two are more intrinsic of the continuous-time approach. In this paper, we have discussed the embedding problem for Markov matrices within the Kimura 3-parameter model (1.). Under a generic assumption (that is, eigenvalues should be different), we have obtained a characterization of the embeddability of the matrices for this model in terms of some inequalities relative to their eigenvalues (Corollary 3.5). Moreover, Theorem 3.9 shows that in this case, rates are identifiable and keep the K3 symmetries, so that embeddability holds if and only if so does $K3$-embeddability, which has been characterized in Theorem 3.3 (2.). As for the problem (3.), Remark 3.7 shows that the product of embeddable matrices within the Kimura 3ST is not embeddable in general. However, under the assumption of different eigenvalues again, problems (3.) and (4.) become equivalent and Theorem 2.5 gives an affirmative answer to both questions.

As a consequence of the characterization of Corollary 3.5, we have been able to compute and compare the volume of embeddable and K3-embeddable matrices within some subspaces of K3 Markov matrices that may be regarded as a compromise solution between the continuous-time approach and the whole space of K3 Markov matrices (see Theorem 4.1 and Table 1).

Despite the fact that the evolutionary submodels of Kimura 3ST (Kimura 2ST and Jukes-Cantor) have repeated eigenvalues, the results obtained here also give information about them. While a characterization of the the embeddability of the Jukes-Cantor model follows easily from the general case of K3 matrices, the embeddability of Kimura 2ST model is more involved and remains an open question. Under this model, we have provided examples of matrices with repeated eigenvalues showing that there are embeddable matrices within the K2ST model that have no Markov generator with K3 form (see Theorem 3.9) and although these matrices are close to saturation, Theorem 2.2 shows they are still biologically relevant.

Another open last issue is the identifiability of rates for K3 embeddable matrices with repeated eigenvalues. In this case, the inequalities (10) for negative eigenvalues and (16) for positive eigenvalues are sufficient for a K3 Markov matrix to have different Markov generators. We believe that the these inequalities are actually necessary, but a much more technical analysis is required. Moreover, we conjecture that under the Kimura 3ST model, the rates of an embeddable matrix that is not $K3$-embeddable are not identifiable. An affirmative answer to these questions would allows us to characterize the embeddability of any K3 Markov matrix (cf. Theorem 3.4), and show that for Markov matrices with positive eigenvalues, embeddability is equivalent to $\mathcal{M}$-embeddability for K3ST and any of its submodels. We defer such a study for a future publication.

6 Acknowledgement

The authors want to thank Marta Casanellas and Jeremy Sumner for conversations held on this topic. They alse wish to thank the reviewers for useful comments and suggestions.

JRL and JFS are partially supported by Spanish government MTM2015-69135-P. Second author is also supported by Generalitat de Catalunya 2014SGR634.