Embeddability of centrosymmetric matrices capturing the double-helix structure in natural and synthetic DNA

Journal of Mathematical Biology (2023) 86:69 https://doi.org/10.1007/s00285-023-01895-8 Mathematical Biology Embeddability of centrosymmetric matrices capturing the double-helix structure in natural and synthetic DNA Muhammad Ardiyansyah1·Dimitra Kosta2·Jordi Roca-Lacostena3 Received: 24 February 2022 / Revised: 3 March 2023 / Accepted: 6 March 2023 / Published online: 5 April 2023 © The Author(s) 2023 Abstract In this paper, we discuss the embedding problem for centrosymmetric matrices, which are higher order generalizations of the matrices occurring in strand symmetric models. These models capture the substitution symmetries arising from the double helix structure of the DNA. Deciding whether a transition matrix is embeddable or not enables us to know if the observed substitution probabilities are consistent with a homogeneous continuous time substitution model, such as the Kimura models, the Jukes-Cantor model or the general time-reversible model. On the other hand, the generalization to higher order matrices is motivated by the setting of synthetic biology, which works with different sizes of genetic alphabets. Keywords Evolutionary model ·Embedding problem ·Markov matrix · Centrosymmetric matrix Mathematics Subject Classification 60J10 ·60J27 ·15B51 ·15A16 ·92D15 BMuhammad Ardiyansyah [email protected] Dimitra Kosta D.K[email protected] Jordi Roca-Lacostena [email protected] 1Department of Mathematics and Systems Analysis, Aalto University, Espoo, Finland 2School of Mathematics, University of Edinburgh and Maxwell Institute for Mathematical Sciences, Edinburgh, UK 3Universitat Politécnica de Catalunya, Barcelona, Catalunya, Spain 123 69 Page 2 of 37 M. Ardiyansyah et al. 1 Introduction Phylogenetics is the study of evolutionary relationships among biological entities, also known as taxa, that aims to infer the evolutionary history among them. In order to model evolution, we consider a directed acyclic graph, called a phylogenetic tree, depicting the evolutionary relationships amongst a selected set of taxa. Phylogenetic trees consist of vertices and edges. Vertices represent taxa, while edges between vertices represent the evolutionary processes between the taxa. In order to describe the real evolutionary process along an edge of a phylogenetic tree, one often assumes that the evolutionary data occurred following a Markov process. A Markov process is a random process in which the future is independent of the past, given the present. Under this Markov process, transitions between nstates given by conditional probabilities are presented in a n×nMarkov matrix M, that is a square matrix whose entries are nonnegative and rows sum to one. A well-known problem in probability theory is the so-called embedding problem which was initially posed by Elfving (Elfving 1937). The embedding problem asks whether given a Markov matrix M, one can find a real square matrix Qwith rows summing to zero and nonnegative off-diagonal entries, such that M=exp(Q). The matrix Qis called a Markov generator. In the complex setting, the embedding problem is completely solved by Higham (2008); a complex matrix Ais embeddable if and only if Ais invertible. However, as our motivation arises from molecular models of evolution we are interested in the embedding problem over the real numbers, so from now on we will denote by Ma real Markov matrix. It was shown by Kingman (1962) that if an n×nreal Markov matrix Mis embeddable, then the matrix Mhas det M>0. Moreover, in the same work by Kingman it was shown that det M>0 is a necessary and sufficient condition for a 2×2 Markov matrix Mto be embeddable. For 3 ×3 Markov matrices a complete solution of the embedding problem is provided in a series of papers (James 1973; Johansen 1974; Carette 1995; Chen and Chen 2011), where the characterisation of embeddable matrices depends on the Jordan decomposition of the Markov matrix. For 4×4Markov matrices the embedding problem is completely settled in a series of papers (Casanellas et al. 2020a,2023; Roca-Lacostena and Fernández-Sánchez 2018a), where similarly to the 3 ×3 case the full characterisation of embeddable matrices is distinguished into cases depending on the Jordan form of the Markov matrices. For the general case of n×nMarkov matrices, there are several results; some presenting necessary conditions (Elfving 1937; Kingman 1962; Runnenberg 1962), while others sufficient conditions (James 1973; Fuglede 1988; Goodman 1970;Davies et al. 2010) for embeddability of Markov matrices. Moreover, the embedding problem has been solved for special n×nmatrices with a biological interest such as equal-input and circulant matrices (Baake and Sumner 2020), group-based models (Ardiyansyah et al. 2021) and time-reversible models (Jia 2016). Despite the fact that there is no theoretical explicit solution for the embeddability of general n×nMarkov matrices, there are results (Casanellas et al. 2023) that enable us to decide whether a n×n Markov matrix with distinct eigenvalues is embeddable or not. This is achieved by providing an algorithm that outputs all Markov generators of such a Markov matrix (Casanellas et al. 2023; Roca-Lacostena 2021). 123 Embeddability of centrosymmetric matrices capturing... Page 3 of 37 69 In this paper, we focus on the embedding problem for n×nmatrices that are symmetric about their center and are called centrosymmetric matrices (see Definition 2). We also study a variation of the famous embedding problem called model embeddability, where apart from the requirement that the Markov matrix is the matrix exponential of a rate matrix, we additionally ask that the rate matrix follows the model structure. For instance, for centrosymmetric matrices, model embeddability means that the rate matrix is also centrosymmetric. The motivation for studying centrosymmetric matrices comes from evolutionary biology, as the most general nucleotide substitution model when considering both DNA strands admits any n×ncentrosymmetric Markov matrix as a transition matrix, where nis the even number of nucleotides. For instance, by considering the four natural nucleotides A-T, C-G we arrive at the strand symmetric Markov model, a well-known phylogenetic model whose substitution probabilities reflect the symmetry arising from the complementarity between the two strands that the DNA is composed of (see (Casanellas and Sullivant 2005)). In particular, a strand symmetric model for DNA must have the following equalities of probabilities in the root distribution: πA=πTand πC=πG(1.1) and the following equalities of probabilities in the transition matrices (θij) θAA =θTT,θ AC =θTG,θ AG =θTC,θ AT =θTA, θCA =θGT,θ CC =θGG,θ CG =θGC,θ CT =θGA. Therefore, the corresponding transition matrices of this model are 4×4 centrosymmetric matrices, usually called strand symmetric Markov matrices in this context. In this article, we will use the terminology 4 ×4 centrosymmetric Markov matrix and strand symmetric Markov matrix interchangeably. In the strand symmetric model there are less restrictions on the way genes mutate from ancestor to child compared to other widely known molecular models of evolution. In fact, special cases of the strand symmetric model are the group-based phylogenetic models such as the Jukes-Cantor (JC) model, the Kimura 2-parameter (K2P) and Kimura 3-parameter (K3P) models. The algebraic structure of strand symmetric models was initially studied in (Casanellas and Sullivant 2005), where it was argued that strand symmetric models capture more biologically meaningful features of real DNA sequences than the commonly used group-based models, as for instance, in any group-based model, the stationary distribution of bases for a single species is always the uniform distribution, while computational evidence in (Yap and Pachter 2004) suggests that the stationary distribution of bases for a single species is rarely uniform, but must always satisfy the symmetries (1.1) arising from nucleotide complementarity, as assumed by the strand symmetric model. In this article, we also explore higher order centrosymmetric matrices for which n>4, which is justified by the use of synthetic nucleotides. One of main goals of synthetic biology is to expand the genetic alphabet to include an unnatural or synthetic base pair. The more letters in a genetic system could possibly lead to an increased potential for retrievable information storage and bar-coding and combinatorial tagging 123 69 Page 4 of 37 M. Ardiyansyah et al. (Benner and Sismour 2005). Naturally the four-letter genetic alphabet consists of just two pairs, A-T and G-C. In 2012, a genetic system comprising of three base pairs was introduced in (Malyshev et al. 2012). In addition to the natural base pairs, the third, unnatural or synthetic base pair 5SICS-MMO2 was proven to be functionally equivalent to a natural base pair. Moreover, when it is combined with the natural base pairs, 5SICS-MMO2 provides a fully functional six-letter genetic alphabet. Namely, six-letter genetic alphabets can be copied (Yang et al. 2007), polymerase chain reaction (PCR)-amplified and sequenced (Sismour et al. 2004; Yang et al. 2011), transcribed to six-letter RNA and back to six-letter DNA (Leal et al. 2015), and used to encode proteins with added amino acids (Bain et al. 1992). This biological importance and relevance of the above six-letter genetic alphabets motivates us to particularly study the 6 ×6 Markov matrices describing the probabilities of changing base pairs in the six-letter genetic system in Sect. 6. When considering both DNA strands, each substitution is observed twice due to the complementarity between both strands, and hence the resulting transition matrix is centrosymmetric. Moreover there are other synthetic analogs to natural DNA which justify studying centrosymmetric matrices for n>6. For instance, hachimoji DNA is a synthetic DNA that uses four synthetic nucleotides B,Z,P,S in addition to the four natural ones A,C,G,T. With the additional four synthetic ones, hachimoji DNA forms four types of base pairs, two of which are unnatural: Pbinds with Zand Bbinds with S.The complementarity between both strands of the DNA implies that the transition matrix is centrosymmetric. Moreover, the research group responsible for the hachimoji DNA system had also studied a synthetic DNA analog system that used twelve different nucleotides, including the four found in DNA (see Yang et al. 2006). Although the biological models which motivate the study of centrosymmetric matrices in this paper require nto be an even number due to the double-helix structure of DNA, in Sect. 5, we include the case of nbeing odd for completeness. Apart from embeddability, that is existence of Markov generators, it is also natural to ask about uniqueness of a Markov generator which is called the rate identifiability problem. Identifiability is a property which a model must satisfy in order for precise statistical inference to be possible. A class of phylogenetic models is identifiable if any two models in the class produce different data distributions. In this article, we further develop the results on rate identifiability of the Kimura two parameter model (Casanellas et al. 2020a) to study rate identifiability for strand symmetric models. We also show that there are embeddable strand symmetric Markov matrices with non identifiable rates, namely the Markov generator is not unique. Moreover, we show that strand symmetric Markov matrices are not generically identifiable, that is, there exists a positive measure subset of strand symmetric Markov matrices containing embeddable matrices whose rates are not identifiable. This paper is organised as follows. In Sect. 2, we introduce the basic definitions and results on embeddability. In Sect. 3, we give a characterisation for a 4 ×4 centrosymmetric Markov matrix Mwith four distinct real nonnegative eigenvalues to be embeddable providing necessary and sufficient conditions in Theorem 2, while we also discuss their rate identifiability property in Proposition 3. Moreover in Sect. 4,using the conditions of our main result Theorem 2, we compute the relative volume of all 4×4 centrosymmetric Markov matrices relative to the 4×4 centrosymmetric Markov 123 Embeddability of centrosymmetric matrices capturing... Page 5 of 37 69 matrices with positive eigenvalues and >0, as well as the relative volume of all 4×4 centrosymmetric Markov matrices relative to the 4×4 centrosymmetric Markov matrices with four distinct eigenvalues and >0. We also compare the results on relative volumes obtained using our method with the algorithm suggested in Casanellas et al. (2023) to showcase the advantages of our method. In Sect. 5, we study higher order centrosymmetric matrices and motivate their use in Sect. 6by exploring the case of synthetic nucleotides where the phylogenetic models admit 6 ×6 centrosymmetric mutation matrices. Finally, Sect. 7discusses implications and possibilities for future work. 2 Preliminaries In this section we will introduce the definitions and results that will be required throughout the paper. We will denote by Mn(K)the set of n×nsquare matrices with entries in the field K=Ror C. The subset of non-singular matrices in Mn(K)will be denoted by GLn(K). Definition 1 AMarkov (or transition) matrix is a non-negative real square matrix with rows summing to one. A rate matrix is a real square matrix with rows summing to zero and non-negative off-diagonal entries. In this paper, we are focusing on a subset of Markov matrices called centrosymmetric Markov matrices. Definition 2 A real n×nmatrix A=(ai,j)is said to be centrosymmetric (CS) if ai,j=an+1−i,n+1−j for every 1 ≤i,j≤n. Definition 2reveals that a CS matrix is nothing more than a square matrix which is symmetric about its center. This class of matrices has been previously studied, for instance, in (Aitken 2017, page 124) and Weaver (1985). Examples of CS matrices for n=5 and n=6, are the following two matrices respectively: ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ a11 a12 a13 a14 a15 a21 a22 a23 a24 a25 a31 a32 a33 a32 a31 a25 a24 a23 a22 a21 a15 a14 a13 a12 a11 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ and ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ a11 a12 a13 a14 a15 a16 a21 a22 a23 a24 a25 a26 a31 a32 a33 a34 a35 a36 a36 a35 a34 a33 a32 a31 a26 a25 a24 a23 a22 a21 a16 a15 a14 a13 a12 a11 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . The class of CS matrices plays an important role in the study of Markov processes since they are indeed transition matrices for some processes in evolutionary biology. For instance, in Kimura (1957), centrosymmetric matrices are used to study the random assortment phenomena of subunits in chromosome division. Furthermore, in Schensted (1958), the same centrosymmetric matrices appear as the transition matrices in 123 69 Page 6 of 37 M. Ardiyansyah et al. the model of subnuclear segregation in the macronucleus of ciliates. Finally, the work (Iosifescu 2014) examines a special case of the random genetic drift phenomenon, which consists of a population of individuals that are able to produce a single type of gamete. In this case, the transition matrices of the associated Markov chain are given by centrosymmetric matrices. The embedding problem is directly related to the notions of matrix exponential and logarithm which we introduce for completeness below. Definition 3 We define the exponential exp(A)of a matrix A, using the Taylor power series of the function f(x)=ex,as exp(A)=∞  k=0 Ak k!, where A0=Inand Indenotes the n×nidentity matrix. If A=Pdiag(λ1,...,λ n) P−1is an eigendecomposition of A, then exp(A)=Pdiag(eλ1,...,eλn)P−1. Given a matrix A∈Mn(K), a matrix B∈Mn(K)is said to be a logarithm of A if exp(B)=A.Ifvis an eigenvector corresponding to the eigenvalue λof A, then vis an eigenvector corresponding to the eigenvalue eλof exp(A). A Markov matrix Mis called embeddable if it can be written as the exponential of aratematrixQ, namely M=exp(Q). Then any rate matrix Qsatisfying the equation M=exp(Q)is called a Markov generator of M. Remark 1 Embeddable Markov matrices occur when we assume a continuous time Markov chain, in which case the Markov matrices have the form M=exp(tQ), where t≥0 represents time and Qis a rate matrix. However, in the rest of the paper, we assume that tis incorporated in the rate matrix Q. The existence of multiple logarithms is a direct consequence of the distinct branches of the logarithmic function in the complex field. Definition 4 Given z∈C\R≤0and k∈Z,thek-th branch of the logarithm of zis logk(z):= log |z|+(Arg(z)+2πk)i, where log is the logarithmic function on the real field and Arg(z)∈(−π,π) denotes the principal argument of z. The logarithmic function arising from the branch log0(z)is called the principal logarithm of zand is denoted as log(z). It is known that if Ais a matrix with no negative eigenvalues, then there is a unique logarithm of Aall of whose eigenvalues are given by the principal logarithm of the eigenvalues of A(Higham 2008, Theorem 1.31). We refer to this unique logarithm as the principal logarithm of A, denoted by Log(A). By definition, the Markov generators of a Markov matrix Mare those logarithms of Mthat are rate matrices. In particular they are real logarithms of M.Thefollowing result enumerates all the real logarithms with rows summing to zero of any 123 Embeddability of centrosymmetric matrices capturing... Page 7 of 37 69 given Markov matrix with positive determinant and distinct eigenvalues. Therefore, all Markov generators of such a matrix are necessarily of this form. Proposition 1 (Casanellas et al. 2023, Proposition 4.3). Let M =Pdiag 1,λ 1,..., λt,μ 1, μ1,...,μ s, μsP−1be an n ×n Markov matrix with P ∈GLn(C) and distinct eigenvalues λi∈R>0for i =1,...,t and μj∈{z∈C: Im(z)>0}for j =1,...,s, all of them pairwise distinct. Then, a matrix Q is a real logarithm of M with rows summing to zero if and only if Q = Pdiag 0,log(λ1),...,log(λt), logk1(μ1), logk1(μ1), . . . , logks(μs), logks(μs) P−1for some k1,...,kj∈Z. Remark 2 In particular, the principal logarithm of Mcan be computed as Log(M)=Pdiag 0,log(λ1),...,log(λt), log(μ1), log(μ1),...,log(μs), log(μs)P−1. In this paper, we focus on the embedding problem for the class of centrosymmetric matrices. In Sect. 3, we will first study the embeddability of 4 ×4 centrosymmetric Markov matrices, which include the K3P, K2P and JC Markov matrices. In Sect. 5 and Sect.6, we will further study the embeddability of higher order centrosymmetric Markov matrices. 3 Embeddability of 4 ×4 centrosymmetric matrices In this section, we begin our study by analyzing the embeddability of 4×4 centrosymmetric matrices also known as strand symmetric matrices. We will provide necessary and sufficient conditions for 4×4 centrosymmetric matrices to be embeddable. Moreover, we will discuss their rate identifiability problem as well. The transition matrices of 4 ×4 centrosymmetric matrices are assumed to have the form M=⎛ ⎜ ⎜ ⎝ m11 m12 m13 m14 m21 m22 m23 m24 m24 m23 m22 m21 m14 m13 m12 m11 ⎞ ⎟ ⎟ ⎠ , where m11 +m12 +m13 +m14 =1=m21 +m22 +m23 +m24 and mij ≥0. Recall that the K3P matrices are assumed to have the form M=⎛ ⎜ ⎜ ⎝ m11 m12 m13 m14 m12 m11 m14 m13 m13 m14 m11 m12 m14 m13 m12 m11 ⎞ ⎟ ⎟ ⎠ . 123 69 Page 8 of 37 M. Ardiyansyah et al. In the case of the K2P matrices, we additionally have m12 =m13, while in the case of JC matrices, m12 =m13 =m14. It can be easily seen that K3P, K2P, and JC Markov (rate) matrices are centrosymmetric. Let us define the following matrix S=⎛ ⎜ ⎜ ⎝ 10 0 1 01 1 0 01−10 10 0 −1 ⎞ ⎟ ⎟ ⎠;(3.1) compare (Casanellas and Kedzierska 2013, Section 6). For a 4 ×4CSMarkovmatrix M, we define F(M):= S−1MS. By direct computation, it can be checked that F(M) is a block diagonal matrix F(M)=⎛ ⎜ ⎜ ⎝ λ1−λ00 1−μμ 00 00αα  00ββ ⎞ ⎟ ⎟ ⎠ ,(3.2) where λ=m11 +m14,μ=m22 +m23, α=m22 −m23,α =m21 −m24, β=m11 −m14,β =m12 −m13. (3.3) Define two matrices, M1:= λ1−λ 1−μμ and M2:= αα  ββ, which are the upper and lower block matrices in (3.2), respectively. Similarly, the rate matrices in strand symmetric models are assumed to have the 4×4 centrosymmetric form Q=⎛ ⎜ ⎜ ⎝ q11 q12 q13 q14 q21 q22 q23 q24 q24 q23 q22 q21 q14 q13 q12 q11 ⎞ ⎟ ⎟ ⎠ , where q11 +q12 +q13 +q14 =0=q21 +q22 +q23 +q24 and qij ≥0fori= j. So, for a 4 ×4CSratematrixQ, we can also define F(Q):= S−1QS. By direct computation, it can be checked that F(Q)=⎛ ⎜ ⎜ ⎝ −ρρ 00 σ−σ00 00δδ  00γγ ⎞ ⎟ ⎟ ⎠ ,(3.4) 123 Embeddability of centrosymmetric matrices capturing... Page 9 of 37 69 where ρ=q12 +q13,σ=q21 +q24, δ=q22 −q23,δ =q21 −q24, γ=q11 −q14,γ =q12 −q13. Define two matrices, Q1:= −ρρ σ−σand Q2:= δδ  γγ, which are the upper and lower block matrices in (3.4), respectively. The following results provide necessary conditions for a 4 ×4CSMarkovmatrix to be embeddable. Lemma 1 Let M =(mij)be a 4×4CS Markov matrix and M =exp(Q)for some CS rate matrix Q. Then 1. m11 +m14 +m22 +m23 >1and 2. (m22 −m23)(m11 −m14)>(m24 −m21)(m13 −m12). Proof We have that F(M)=S−1MS =S−1exp (Q)S=exp(S−1QS)=exp (F(Q)). Then M10 0M2=exp(F(Q)) =exp(Q1)0 0exp(Q2). Thus, M1is an embeddable 2 ×2 Markov matrix. Using the embeddability criteria of 2 ×2 Markov matrices in Kingman (1962), we have that 1 <tr(M1)=λ+μ, which is the desired inequality. Additionally, since M2=exp(Q2),det(M2)>0as desired.  Lemma 2 Let M =(mij)be a 4×4CS Markov matrix and M =exp(Q)for some CS rate matrix Q =(qij).Ifλ+μ= 2, then q12 +q13 =−λ+1 λ+μ−2ln(λ +μ−1) and q21 +q24 =−μ+1 λ+μ−2ln(λ +μ−1). Proof By direct computations and the proof of Lemma 1, M1=exp(Q1)=1 ρ+σe−ρ−σρ+σ−e−ρ−σρ+ρ −e−ρ−σσ+σe−ρ−σσ+ρ. 123 69 Page 16 of 37 M. Ardiyansyah et al. Theorem 1). Moreover, they are also equivalent to the restriction to centrosymmetricmatrices of the embeddability criteria for 4 ×4 Markov matrices with different eigenvalues given in (Casanellas et al. 2023, Theorem 1.1) In the last part of this section, we discuss the rate identifiability problem for 4×4 centrosymmetric matrices. If a centrosymmetric Markov matrix arises from a continuous-time model, then we want to determine its corresponding substitution rates. In other words, given an embeddable 4 ×4 CS matrix, we want to know if we can uniquely identify its Markov generator. It is worth noting that Markov matrices with repeated real eigenvalues may admit more than one Markov generator (e.g. examples 4.2 and 4.3 in (Casanellas et al. 2020a) show embeddable K2P matrices with more than one Markov generator). Nonetheless, this is not possible if the Markov matrix has distinct eigenvalues, because in this case its only possible real logarithm would be the principal logarithm (Culver 1966). As one considers less restrictions in a model, the measure of the set of matrices with repeated real eigenvalues decreases, eventually becoming a measure zero set. For example, this is the case within the K3P model, where both its submodels (the K2P model and the JC model) consist of matrices with repeated eigenvalues and have positive measure subsets of embeddable matrices with non-identifiable rates. However, when considering the whole set of K3P Markov matrices, the subset of embeddable matrices with more than one Markov generator has measure zero (see Chapter 4 in (Roca-Lacostena 2021)). Nevertheless, this behaviour only holds if the Markov matrices within the model have real eigenvalues. Proposition 3 There is a positive measure subset of 4×4CS Markov matrices that are embeddable and whose rates are not identifiable. Moreover, all the Markov generators of the matrices in this set are also CS matrices. Proof Given P=⎛ ⎜ ⎜ ⎝ 1−51−i1+i 12 −ii 12 i−i 1−5−1+i−1−i ⎞ ⎟ ⎟ ⎠ , let us consider the following matrices M=Pdiag(1,e−7π,e−4πi,−e−4πi)P−1,Q=Pdiag(0,−7π, −4π−3π 2i,−4π+3π 2i)P−1. A straightforward computation shows that Mis a CS Markov matrix and Qis a CS rate matrix. Moreover they both have non-zero entries. By applying the exponential series to Q, we get that exp(Q)=M. This means that Mis embeddable and Qis a Markov generator of M. Since Qis a rate matrix, so is Qt for any t∈R≥0. Therefore, exp(Qt)is an embeddable Markov matrix, because the exponential of any rate matrix is necessarily a Markov matrix. See (Pachter and Sturmfels 2005, Theorem 4.19) for more details. 123 Embeddability of centrosymmetric matrices capturing... Page 17 of 37 69 Moreover, we have that S−1P=⎛ ⎜ ⎜ ⎝ 1−50 0 12 0 0 00 −ii 001−i1+i ⎞ ⎟ ⎟ ⎠ , so S−1exp(Qt)Sis a 2-block diagonal matrix. Hence, by Proposition 2we have that exp(Qt)is an embeddable strand symmetric Markov matrix for all t∈R>0. Now, let us define V=Pdiag(0,0,2πi,−2πi)P−1. Note that Qand Vdiagonalize simultaneously via Pand hence they commute. Therefore, exp(Q+V)=exp(Q)exp(V)=MI 4=M by the Baker-Campbell-Haussdorff formula. Moreover, exp(Qt +kV)=exp(Qt)exp(kV)=exp(Qt)I4=exp(Qt) for all k∈Z. Note that kV is a bounded matrix for any given kand hence, given t large enough, it holds that Qt +mV is a rate matrix for any mbetween 0 and k. This shows that, for tlarge enough, exp(Qt)is an embeddable CS Markov matrix with at least k+1 different CS Markov generators. Moreover, exp(Qt)and all its generators have no null entries by construction and they can therefore be perturbed as in Theorem 3.3 in Casanellas et al. (2020b) to obtain a positive measure subset of embeddable CS Markov matrices that have k+1 CS Markov generators.  Remark 4 The perturbation presented in Theorem 3.3 in Casanellas et al. (2020b) consists of small changes on the real and complex parts of the eigenvalues and eigenvectors of Mother than the eigenvector (1,...,1)and its corresponding eigenvalue 1. If those changes are small enough then the resulting transition matrix and all its generators still satisfy the stochastic constraints of Markov/rate matrices. Remark 5 Using the same notation as in the proposition above and given C∈GL2(C), let us define Q(C)=Pdiag(I2,C)diag 1,−7π, −4π−3π 2i,−4π+3π 2idiag(I2,C−1)P−1. Since Q(I2)=Qis a CS rate matrix with no null entries, so is Q(C)for C∈ GL2(C)close enough to I2. Moreover, by construction we have that exp(2tQ(C)) = exp(2tQ)for all t∈N. Therefore, for t∈Nwe have that exp(2tQ)has uncountably many Markov generators (i.e. 2tQ(C)with Cclose to I2) and all of them are CS matrices (Culver 1966, Corollary 1). It is worth noting that according to (Culver 1966, Corollary 1), if a matrix has uncountably many logarithms, then it necessarily has repeated real eigenvalues. Therefore, the subset of embeddable CS Markov matrices with uncountably many generators has measure zero within the set of all matrices. 123 69 Page 18 of 37 M. Ardiyansyah et al. 4 Volumes of 4 ×4 CS Markov matrices In this section, we compute the relative volumes of embeddable 4 ×4CSMarkov matrices within some meaningful subsets of Markov matrices. The aim of this section is to describe how large the different sets of matrices are compared to each other. Let VMarkov 4be the set of all 4 ×4 CS Markov matrices. We use the following description VMarkov 4={(b,c,d,e,g,h)T∈R6:b,c,d,e,g,h≥0,1−b−c−d≥0, 1−e−g−h≥0}. More explicitly, we identify the 4 ×4 CS Markov matrix ⎛ ⎜ ⎜ ⎝ 1−b−c−db c d e1−e−g−hg h hg1−e−g−he dcb1−b−c−d ⎞ ⎟ ⎟ ⎠ with a point (b,c,d,e,g,h)∈VMarkov 4.Let V+be the set of all CS Markov matrices having real positive eigenvalues, where =((1−e−2g−h)−(1−b−c−2d))2+4(e−h)(b−c), is the discriminant of the matrix M2as stated in Sect. 3.WehaveV+⊆VMarkov 4. More explicitly, V+={(b,c,d,e,g,h)∈R6:b,c,d,e,g,h≥0, 1−b−c−d≥0,1−e−g−h≥0,1−b−c−e−h>0, (2−b−c−2d−e−2g−h)+>0,(2−b−c−2d−e−2g−h)−>0,>0}. Let Vem+be the set of all embeddable 4 ×4 CS Markov matrices with four distinct real positive eigenvalues. We have Vem+⊆V+. Therefore, by Theorem 2, Vem+={(b,c,d,e,g,h)∈R6:b,c,d,e,g,h≥0,1−b−c−d≥0,1−e−g−h≥0, 1−b−c−e−h>0, (2−b−c−2d−e−2g−h)+>0,(2−b−c−2d−e−2g−h)−>0,>0, |φ(0)|≤−α1,|ε(0)|≤−β1,δ(0)≤β1,γ(0)≤α1}. Finally, we consider the following two biologically relevant subsets of VMarkov 4. Let VDLC be the set of diagonally largest in column (DLC) Markov matrices, which is the subset of VMarkov 4containing all CS Markov matrices such that the diagonal element is the largest element in each column. These matrices are related to matrix parameter identifiability in phylogenetics (Chang 1996). Secondly, we let VDD be the set of diagonally dominant (DD) Markov matrices, which is the subset of VMarkov 4 matrices containing all CS Markov matrices such that in each row the diagonal element is at least the sum of all the other elements in the row. Biologically, the subspace VDD 123 Embeddability of centrosymmetric matrices capturing... Page 19 of 37 69 consists of matrices with probability of not mutating at least as large as the probability of mutating. If a diagonally dominant matrix is embeddable, it has an identifiable rate matrix (Cuthbert 1972; James 1973). By the definition of each set, we have the inclusion VDD ⊆VDLC. Remark 6 The sets V+,Vem+,VDLC,VDD that we consider in this section are all subsets of the set VMarkov 4of all 4×4 CS Markov matrices, but we can use the same definition to refer to the equivalent subsets of n×nCS Markov matrices. Therefore, we will use the same notation V+,Vem+,VDLC,VDD to refer to the equivalent subsets of the set VMarkov nof n×nCS Markov matrices without confusion in the following sections. In the rest of this section, the number v(A)denotes the Euclidean volume of the set A. By definition, VMarkov 4,VDLC and VDD are polytopes, since they are defined by the linear inequalities in R6. Hence, we can use Polymake (Gawrilow and Joswig 2000) to compute their exact volumes and obtain that v(VMarkov 4)=1 36,v(VDLC)=1 576 and v(VDD)=1 2304. Hence, we see that VDLC and VDD constitute roughly only 6.25% and 1.56% of VMarkov 4, respectively. On the other hand, we will estimate the volume of the sets V+,Vem+,VDLC ∩ V+,VDLC ∩Vem+,VDD ∩V+,and VDD ∩Vem+using the hit-and-miss Monte Carlo integration method (Hammersley 2013) with sufficiently many sample points in Mathematica (Inc. 2022). Theoretically, Theorem 2enables us to compute the exact volume of these relevant sets. For example in the case of K3P matrices, such exact computation of volumes has been feasible in Roca-Lacostena and FernándezSánchez (2018b). However, while for the K3P matrices, the embeddability criterion is given by three quadratic polynomial inequalities, in the case of CS matrices the presence of nonlinear and nonpolynomial constraints imposed on each set, makes the exact computation of the volume of these sets intractable. Therefore, we need to approximate the volume of these sets. Given a subset A⊆VMarkov 4, the volume estimate of v(A)computed using the hit-and-miss Monte Carlo integration method with nsample points is given by the number of points belonging to Aout of nsample points. For computational purposes, in the formula of φ(0)and ε(0), we use the fact that y−z=log (2−b−c−2d−e−2g−h)+√ (2−b−c−2d−e−2g−h)−√. =log ((2−b−c−2d−e−2g−h)+√)2 (2−b−c−2d−e−2g−h)2−. =log ((2−b−c−2d−e−2g−h)+(b+c+2d−e−2g−h)2+4(e−h)(b−c))2 (2−b−c−2d−e−2g−h)2−((b+c+2d−e−2g−h)2+4(e−h)(b−c))  All codes for the computations implemented Mathematica and Polymake can be found at the following address: https://github.com/ardiyam1/Embeddability-andrate-identifiability-of-centrosymmetric-matrices. 123 69 Page 20 of 37 M. Ardiyansyah et al. Table 1 Number of samples in the sets n104105106107 Samples in VMarkov 4280 2767 27829 277628 Samples V+23 192 1999 20601 Samples in Vem+3 34 359 3511 Samples in VDLC ∩V+19 154 1541 15830 Samples in VDD ∩V+3 31 262 2889 Samples in VDLC ∩Vem+3 34 357 3503 Samples in VDD ∩Vem+1 15 105 1011 V+,Vem+,VDLC ∩V+,VDLC ∩Vem+,VDD ∩V+and VDD ∩Vem+ using hit-and-miss methods and Theorem 2. Table 2 Relative volumes ratio between the relevant subsets obtained using hit-and-miss method and Theorem 2. The volumes were estimated as the quotient of the sample sizes in Table 1 n104105106107 v(Vem+) v(V+)0.130435 0.177083 0.17959 0.170429 v(VDLC∩Vem+) v(VDLC∩V+)0.157895 0.220779 0.231668 0.221289 v(VDLC∩Vem+) v(V+)0.130435 0.177083 0.178589 0.17004 v(VDLC∩Vem+) v(Vem+)1 1 0.994429 0.997721 v(VDD∩Vem+) v(VDD∩V+)0.333333 0.483871 0.400763 0.349948 v(VDD∩Vem+) v(V+)0.0434783 0.078125 0.052563 0.0490753 v(VDD∩Vem+) v(Vem+)0.333333 0.441176 0.292479 0.287952 The results of these estimations using the hit-and-miss Monte Carlo integration implemented in Mathematica with nsample points are presented in Table 1, while Table 2provides an estimated volume ratio between relevant subsets of centrosymmetric Markov matrices using again the hit-and-miss Monte Carlo integration with n sample points. In Table 1, we firstly generate ncentrosymmetric matrices whose offdiagonal entries were sampled uniformly in [0,1]and forced the rows of the matrix to sum to one. Out of these nmatrices, we test how many of them are actually Markov matrices (i.e. the diagonal entries are non-negative) and then out of these how many have positive eigenvalues. In particular, for n=107sample points containing 277628 centrosymmetric Markov matrices, Table 2suggests that there are approximately 1.7% of centrosymmetric Markov matrices with distinct positive eigenvalues that are embeddable. Moreover, we can see that for n=107, out of all embeddable centrosymmetric Markov matrices with distinct positive eigenvalues, almost all are diagonally largest in column, while only 28% are diagonally dominant. An alternative approach for approximating the number of embeddable matrices within the model is to use Algorithm 5.8 in Casanellas et al. (2023) to test the embeddability of the sample points. Tables 4and 5below are analogous to Tables 1and 2,but Table 4was obtained using the sampling method in (Roca-Lacostena 2021, Appendix 123 Embeddability of centrosymmetric matrices capturing... Page 21 of 37 69 Table 3 Number of samples in V+,VDLC ∩V+,and VDD ∩V+obtained by using the sampling method in (Roca-Lacostena 2021, Appendix A) Samples in V+104105106107 Samples in VDLC ∩V+8531 85446 854709 8549100 Samples in VDD ∩V+1464 14538 144546 1448720 Table 4 Number of samples in Vem+,VDLC ∩Vem+and VDD ∩Vem+obtained by applying either Theorem 2or the results in Casanellas et al. (2023) on the sample set in Table 3 Samples in Vem+1877 18663 185357 1862413 Samples in VDLC ∩Vem+1869 18586 184555 1854592 Samples in VDD ∩Vem+516 5164 50058 504304 Table 5 Relative volumes ratio between the relevant subsets obtained using hit-and-miss method and either Algorithm 5.8 in Casanellas et al. (2023)or Theorem 2. The volumes were estimated as the quotient of the sample sizes in Tables 3and 4 n104105106107 v(Vem+) v(V+)0.1877 0.18663 0.185357 0.1862413 v(VDLC∩Vem+) v(VDLC∩V+)0.2191 0.2175 0.2159 0.2169 v(VDLC∩Vem+) v(V+)0.1869 0.18586 0.184555 0.1854592 v(VDLC∩Vem+) v(Vem+)0.9957 0.9959 0.99567 0.99580 v(VDD∩Vem+) v(VDD∩V+)0.3524 0.3552 0.3463 0.3481 v(VDD∩Vem+) v(V+)0.0516 0.05164 0.050058 0.0504 v(VDD∩Vem+) v(Vem+)0.2749 0.2767 0.2701 0.2708 A), while using either Algorithm 5.8 in Casanellas et al. (2023) or the inequalities in Theorem 2yields identical results which are provided in Tables 4and 5. We used the python implementation of Algorithm 5.8 in Casanellas et al. (2023) provided in (Roca-Lacostena 2021, Appendix A) and modified it to sample on the set of 4×4 CS Markov matrices with positive eigenvalues. The original sampling method used in (Roca-Lacostena 2021, Appendix A) consisted of sampling uniformly on the set of 4 ×4 centrosymmetric-Markov matrices until we obtained nsamples (or as many samples as we require) with positive eigenvalues. Despite the fact that Theorem 2and Algorithm 5.8 in Casanellas et al. (2023)were originally implemented using different programming languages (Wolfram Mathematica and Python respectively) and were tested with different sample sets, the results obtained are quite similar as illustrated by Tables 2and 5. In fact, when we apply both Algorithm 5.8 in Casanellas et al. (2023) and Theorem 2on the same sample set in Table 3, we obtain identical results which are displayed in Tables 4and 5. It is worth noting that the embeddability criteria given in Theorem 2use inequalities depending on the entries of the matrix, whereas Algorithm 5.8 in Casanellas et al. (2023) relies on the computation of its principal logarithm and its eigenvalues 123 69 Page 22 of 37 M. Ardiyansyah et al. Table 6 Running times for the Python implementation of the embeddability criterion arising from Theorem 2and from Algorithm 5.8 in Casanellas et al. (2023). The simulations were run using a computer with 8GB of memory 104105106107 Sampling time 12.5s 121.5s (2 min) 1222s (20min) 12141.8s (3h 22min) Embedding criteria (Theorem 2) 28.3s 273.2s (4min 30s) 2703s (45min) 27413s (7h 37min) Embedding criteria (Algorithm 5.8) 84.2s 840.5s (15 min) 8358 (2h 19min) 83786s (23h 16min) Table 7 Embeddable matrices within 4 ×4 CS Markov matrices and its intersection with DLC matrices and DD matrices Samples Embeddable samples Proportion of embeddable VMarkov 4107173455 0.0173455 VDLC 1021195 172380 0.1688022 VDD 156637 49471 0.3158321 and eigenvector, which may cause numerical issues when working with matrices with determinant close to 0. Moreover, the computation of logarithms can be computationally expensive. As a consequence, the algorithm implementing the criterion for embeddability arising from Theorem 2is faster. Table 6shows the running times for the implementation of both embeddability criteria used to obtain Table 5. The Python implementation of Algorithm 5.8 in Casanellas et al. (2023) provided in (Roca-Lacostena 2021, Appendix A) can also be used to test the embeddability of any 4 ×4 CS Markov matrix (including those with non-real eigenvalues) without modifying the embeddability criteria. To do so, it is enough to apply the algorithm to a set of Markov matrices with different eigenvalues sampled uniformly from the set of all 4 ×4 CS Markov matrix. As hinted in Remark 3, this would also be possible using the embedability criterion in Theorem 2together with the boundaries for kprovided in (Casanellas et al. 2023, Theorem 5.5). Table 7shows the results obtained when applying Algorithm 5.8 in Casanellas et al. (2023) to a set of 1074×4CSMarkov matrices sampled uniformly. As most DLC and DD matrices have positive eigenvalues, the proportion of embeddable matrices within these subsets is almost the same when admitting matrices with non-positive eigenvalues (as in Table 7instead of only considering matrices with positive eigenvalues as we did in Tables 2and 5. On the other hand, the proportion of 4×4 embeddable CS matrices is much smaller in this case. 5 Centrosymmetric matrices and generalized Fourier transformation In Sects. 3and 4we have seen the embeddability criteria for 4 ×4 centrosymmetric Markov matrices and the volume of their relevant subsets. In this section, we are extending this framework to larger matrices. The importance of this extension is rel123 Embeddability of centrosymmetric matrices capturing... Page 23 of 37 69 evant to the goal of synthetic biology which aims to expand the genetic alphabet. For several decades, scientists have been cultivating ways to create novel forms of life with basic biochemical components and properties far removed from anything found in nature. In particular, they are working to expand the number of amino acids which is only possible if they are able to expand the genetic alphabet (see for example (Hoshika et al. 2019)). 5.1 Properties of centrosymmetric matrices For a fixed n∈N,letVndenote the set of all centrosymmetric matrices of order n.Moreover, let VMarkov nand Vrate ndenote the set of all centrosymmetric Markov and rate matrices of order n, respectively. As a subspace of the set of all n×nreal matrices, for neven, dim(Vn)=n2 2while for nodd, dim(Vn)=n 2(n+1)+1. We will now mention some geometric properties of the sets VMarkov nand Vrate n.Furthermore, for any real number x,xand xdenote the floor and the ceiling function of x, respectively. Proposition 4 1. For n even, V Markov n⊆R n(n−1) 2 ≥0is a Cartesian product of n 2standard (n−1)-simplices and its volume is 1 (n−1)!n 2 .For n odd, V Markov n⊆Rn 2n ≥0is a Cartesian product of n 2standard (n−1)-simplices and the n 2-simplex with vertices {0,ei 2}1≤i≤n 2∪{en 2+1}, where eiis the i-th standard unit vector in Rn. Hence, the volume of V Markov nis 1 2n 2(n 2)!(n−1)!n 2. 2. For n even, V rate n=R n(n−1) 2 ≥0and for n odd, V rate n=Rn 2n ≥0. Proof Here we consider the following identification for an n×ncentrosymmetric matrix M.Forneven, Mcan be thought as a point (M1,...,Mn 2)∈(Rn ≥0)n 2where the point Mi∈Rn ≥0corresponds to the i-th row of M. Similarly, for nodd, we identify Mas a point in (Rn ≥0)n 2×Rn 2+1 ≥0. Since Mis a Markov matrix, under this identification, each point Milies in some simplices. Therefore, VMarkov nis a Cartesian product of some simplices. For neven, these simplices are the standard (n−1)-dimensional simplex: x1+···+xn=1, xi≥0,1≤i≤n⇔x1+···+xn−1≤1, xi≥0,1≤i≤n−1(5.1) For nodd and 1 ≤i≤n 2, the point Mibelongs to standard (n−1)-simplex above and the point Mn 2+1belongs to the simplex 2x1+···+2xn 2+xn 2+1=1 xi≥0,1≤i≤n 2+1⇔x1+···+xn 2≤1 2, xi≥0,1≤i≤n 2(5.2) We now compute the volume of VMarkov n. Let us recall the fact that the volume of the Cartesian product of spaces is equal to the product of volumes of each factor 123 69 Page 24 of 37 M. Ardiyansyah et al. space if the volume of each factor space is bounded. Moreover, the (n−1)-dimensional volume of the standard simplex in Eq. (5.1)inRn−1is 1 (n−1)!.For neven, the statement follows immediately. For nodd, we use the fact that the n 2-dimensional volume of the simplex in Eq. (5.2)is 1 2n 2(n 2)!.We refer the reader to Stein (1966) for an introductory text on the volume of simplices. For the second statement, we use the fact that if Qis a rate matrix, then qii = −j=iqij where qij ≥0fori= j. In the rest of this section, let Jnbe the n×nanti-diagonal matrix, i.e. the (i,j)- entries are one if i+j=n+1 and zero otherwise. The following proposition provides some properties of the matrix Jnthat can be checked easily. Proposition 5 Let A =(aij)∈Mn(R). Then 1. (AJn)ij =ai,n+1−jand (JnA)ij =an+1−i,j. 2. A is a centrosymmetric matrix if only if JnAJn=A. In Sect. 3, we have seen that 4 ×4 CS matrices can be block-diagonalized through the matrix S. Now we will present a construction of generalized Fourier matrices to block-diagonalize any centrosymmetric matrices. Let us consider the following recursive construction of the n×nmatrix Sn: S1=1,S2=11 1−1and Sn:= ⎛ ⎝ 10 1 0Sn−20 10−1⎞ ⎠,for n≥3.(5.3) Proposition 6 For each natural number n ≥3,S nis invertible and its inverse is given by S−1 n=⎛ ⎝ 1 201 2 0S−1 n−20 1 20−1 2 ⎞ ⎠. Proof The proposition easily follows from the definition of Sn. Indeed, we have ⎛ ⎝ 1 201 2 0S−1 n−20 1 20−1 2 ⎞ ⎠Sn=⎛ ⎝ 1 201 2 0S−1 n−20 1 20−1 2 ⎞ ⎠⎛ ⎝ 10 1 0Sn−20 10−1⎞ ⎠=In.  The following proposition provides another block decomposition of the matrix Sn and its inverse. Proposition 7 Let n ≥2. 1. For n even, Sn=In 2Jn 2 Jn 2−In 2, while for n odd, Sn=⎛ ⎝ In 20Jn 2 01 0 Jn 20−In 2⎞ ⎠. 123 Embeddability of centrosymmetric matrices capturing... Page 25 of 37 69 2. Using these block partitions, S−1 n=1 2Snfor n even, while S−1 n=⎛ ⎜ ⎝ 1 2In 201 2Jn 2 01 0 1 2Jn 20−1 2In 2 ⎞ ⎟ ⎠ for n odd. Proof The proof follows from induction on nand the fact that J2 n=In. We will call a vector v∈Rnsymmetric if vi=vn+1−ifor every 1 ≤i≤n, i.e. Jnv=v. Moreover, we call a vector w∈Rnanti-symmetric if vi=−vn+1−ifor every 1 ≤i≤n, i.e. Jnv=−v. The following technical proposition will be used in what follows in order to simplify a centrosymmetric matrix. Proposition 8 Let n ≥2. Let v∈Rnbe a symmetric vector and w∈Rnbe an anti-symmetric vector. 1. The last n 2entries of Snvand vTSnare zero. Similarly, the last n 2entries of S−1 nvand vTS−1 nare zero. 2. The first n 2entries of Snwand wTSnare zero. Similarly, the first n 2entries of S−1 nwand wTS−1 nare zero. 3. Then the sum of the entries of Snvand vTSnis the sum of the entries of v. 4. Then the sum of the entries of S−1 nvand vTS−1 nis the sum of the first n 2entries of v. Proof We will only prove the first part of item (1) in the proposition using mathematical induction on n.The base case for n=2 can be easily obtained. Suppose now that the proposition holds for all k<n.Let v=⎛ ⎝ v1 v v1⎞ ⎠∈Rnbe a symmetric element. Then v∈Rn−2is also symmetric. By direct computation we obtain Snv=⎛ ⎝ 10 1 0Sn−20 10−1⎞ ⎠⎛ ⎝ v1 v v1⎞ ⎠=⎛ ⎝ 2v1 Sn−2v 0⎞ ⎠. The last n−2 2entries of Sn−2vare zero. Thus, the last n−2 2+1=n 2entries of Snvare zero as well. The proof of the other statements can be obtained analogously using induction. In particular, let us note that the proof given for item (1) directly implies item (3).  For a fixed number n, let us define the following map: Fn:Mn(R)→Mn(R) A→ Fn(A):= S−1 nASn. For n=4,we have seen that if Ais a CS matrix, then F4(A)is a block-diagonal matrix where each block is of size 2 ×2 and is given by A1and A2. Moreover, the upper block is a Markov matrix. The following lemma provides a generalization to these results. 123 69 Page 32 of 37 M. Ardiyansyah et al. Case 2 In this case Ahas exactly one conjugate pair of complex eigenvalues and we obtain the following criterion by adapting Corollary 5.6 in Casanellas et al. (2023) to our framework: Proposition 13 Given the matrix V := Pdiag(0,0,0,0,2πi,−2πi)P−1define: L:= max (i,j):i=j,Vi,j>0−Log(A)i,j Vi,j,U:= min (i,j):i=j,Vi,j<0−Log(A)i,j Vi,j and set N:= {(i,j):i= j,Vi,j=0and Log(A)i,j<0}.Then, 1. A is embeddable if and only if N=∅and L≤U. 2. the set of Markov generators for A is Q=Log(A)+kV :k∈Zsuch that L≤k≤U. Proof The proof of this theorem is analogous to the proof of Theorem 5.5 in Casanellas et al. (2020a) but considering the matrix Vas defined here. According to Proposition 1, any Markov generator of Ais of the form Logk(A)=Pdiag(0,log(λ1), log(λ2), log(μ), logk(γ1), logk(γ1))P−1 =Pdiag(0,log(λ1), log(λ2), log(μ), logk(γ1)+2πki,logk(γ1)−2πki)P−1. Such a logarithm can be rewritten as Log(A)+kV. Using this, we will prove that Logk(A)=Log(A)+kV is a rate matrix if and only if N=∅and L≤k≤U. Suppose that there exists k∈Zsuch that Logk(A)is a rate matrix. Hence, Log(A)i,j+kVi,j≥0 for all i= j.Fori= j,wehave: (a) Log(A)i,j≥0 for all i= jsuch that Vi,j=0. This means that N=∅. (b) −Log(A)i,j Vi,j≤kfor all i= jsuch that Vi,j>0. This means that L≤k. (c) −Log(A)i,j Vi,j≥kfor all i= jsuch that Vi,j<0. This means that k≤U. Conversely, suppose that N=∅and and that there is k∈Zsuch that L≤k≤U.We want to check that Logk(A)is a rate matrix. According to Proposition 1, each row of Logk(A)sums to 0. Moreover, for i= j,wehave: (a) if Vi,j=0, then Logk(A)i,j=Log(A)i,j. Since N=∅,Logk(A)i,j=Log(A)i,j ≥0. (b) if Vi,j>0, then Logk(A)i,j=Log(A)i,j+kVi,j≥Log(A)i,j+LVi,j≥Log(A)i,j +(−Log(A)i,j Vi,j)Vi,j=0. (c) if Vi,j<0, then −Logk(A)i,j=−Log(A)i,j−kVi,j≤−Log(A)i,j−UVi,j ≤−Log(A)i,j−(−Log(A)i,j Vi,j)Vi,j=0. The proof is now complete.  123 Embeddability of centrosymmetric matrices capturing... Page 33 of 37 69 Case 3 As in Case 2, Ahas exactly one conjugate pair of eigenvalues and hence its embeddability (and all its generators) can be determined by using Proposition 13 but defining the matrix Vas V=Pdiag(0,0,0,0,2πi,−2πi)P−1. However in Case 3 the conjugate pair of eigenvalues lie in A1which is a Markov matrix. This allows us to use the results regarding the embeddability of 3 ×3 Markov matrices to obtain an alternative criterion to test the embeddability of A. To this end we define Log−1(A):= Pdiag(0,z,z,log(μ), log(γ1)log(γ2)) P−1(6.2) where z:= log−1(λ1). Proposition 14 The matrix A is embeddable if and only if Log(A)or Log−1(A)are rate matrices. Proof Note that exp(Log(A)) =exp(Log−1(A)) =Aso one of the implications is immediate to prove. To prove the other implication, we assume that Ais embeddable and let Qbe a Markov generator for it. Proposition 1yields that Q=Pdiag(0,logk1(λ1), logk2(λ2), logk3(μ), logk4(γ1), logk5(γ2)) P−1, for some integers k1,...,k5∈Z. Therefore, F(Q)=Q10 0Q2where Q1and Q2 are real logarithms of A1and A2respectively. Since A2is a real matrix with distinct positive eigenvalues, its only real logarithm is its principal logarithm. This implies that k3=k4=k5=0 (so that Q2=Log(A2)). Now, recall that A1is a Markov matrix (see Lemma 5). Using Proposition 1again, we obtain that Q1is a rate matrix, thus A1is embeddable. To conclude the proof it is enough to recall Theorem 4 in James (1973), which yields that A1is embeddable if and only if Log(A1)or P1diag(0,z,z)P−1 1is a rate matrix.  Case 4 In this case, the solution to the embedding problem can be obtained as a byproduct of the results for the previous cases: Proposition 15 Let Log0,0(A)denote the principal logarithm of A and Log−1,0(A) denote the matrix in (6.2). Given the matrix V := Pdiag(0,0,0,0,2πi,−2πi)P−1 and k ∈{0,−1}define: Lk:= max (i,j):i=j,Vi,j>0−Logk,0(A)i,j Vi,j,Uk:= min (i,j):i=j,Vi,j<0−Logk,0(A)i,j Vi,j and set Nk:= {(i,j):i= j,Vi,j=0and Logk,0(A)i,j<0}.Then, 1. A is embeddable if and only if Nk=∅and Lk≤Ukfor k =0or k =−1. 123 69 Page 34 of 37 M. Ardiyansyah et al. 2. If A is embeddable, then at least one of its Markov generator can be written as Logk,k2(A):= Pdiag(0,logk(λ1), logk(λ1), log(μ), logk2(γ1), logk2(γ1)) P−1 with k ∈{0,−1}and k2∈Zsuch that Lk≤k2≤Uk. Proof The matrix Ais embeddable if and only if it admits a Markov generator. According to Proposition 1, if such a generator Qexists then it can be written as Logk1,k2(A) for some k1,k2∈Z. Therefore, Lemma 5implies that F(A)=A10 0A2for some matrices A1and A2. Moreover, F(Q)=Q10 0Q2where Q1and Q2are real logarithms of A1and A2respectively. As shown in the proof of Proposition 14,A1is actually a Markov matrix and Q1 is a Markov generator for it (see also Lemma 5). Moreover, by Theorem 4 in James (1973), A1is embeddable if and only if Log(A1)or Log−1(A1)are rate matrices. This implies that Logk1,k2(A)is a rate matrix if and only if Log0,k2(A)or Log−1,k2are rate matrices. To conclude the proof we proceed as in the proof of Proposition 13. Indeed, note that for k∈{0,−1},Logk,k2(A)=Logk,0(A)+k2V. Using this, it is immediate to check that Logk,k2(A)is a rate matrix if and only if Nk=∅and Lk≤k2≤Uk. 7 Discussion The central symmetry is motivated by the complementarity between both strands of the DNA. When a nucleotide substitution occurs in one strand, there is also a substitution between the corresponding complementary nucleotides on the other strand. Therefore, working with centrosymmetric Markov matrices is the most general approach when considering both DNA strands. In this paper, we have discussed the embedding problem for centrosymmetric Markov matrices. In Theorem 2, we have obtained a characterization of the embeddabilty of 4 ×4 centrosymmetric Markov matrices which are exactly the strand symmetric Markov matrices. In particular, we have also shown that if a 4 ×4CS Markov matrix is embeddable, then any of its Markov generators is also a CS matrix. Furthermore, In Sect. 6, we have discussed the embeddability criteria for larger centrosymmetric matrices. As a consequence of the characterization of Theorem 2, we have been able to compute and compare the volume of the embeddable 4 ×4 CS Markov matrices within some subspaces of 4 ×4 CS Markov matrices. These volume comparisons can be seen in Table 2and Table 7. For larger matrices, using the results in Sect. 6, we have estimated the proportion of embeddable matrices within the set of all 6 ×6 centrosymmetric Markov matrices and within the subsets of DLC and DD matrices. This is summarized in Table 9below. The computations were repeated several times obtaining results with small differences in the values but the same order of magnitude and starting digits. 123 Embeddability of centrosymmetric matrices capturing... Page 35 of 37 69 Table 9 Relative volume of embeddable matrices within relevant subsets of 6×6 centrosymmetric Markov matrices. The results were obtained using the hit-and-miss Monte Carlo integration with 107sample points Set Sample points Embeddable sample points Rel. vol. of embeddable matrices VMarkov 61081370 0.0000137 VDLC 1034607 1362 0.0013164 VDD 3048 84 0.0275590 As we have seen in Sect. 3and 6, we have only considered in detail the embeddability of CS Markov matrices of size n=4 and n=6. We expect that the proportion of the embeddable CS Markov matrices within the subset of Markov matrices in larger dimension tends to zero as ngrows larger as indicated by Tables 2,7,8, and 9. These results together with the results obtained for the strand symmetric model (see Table 7) indicate that restricting to homogeneous Markov processes in continuoustime is a very strong restriction because non-embeddable matrices are discarded and their proportion is much larger than that of embeddable matrices. For instance, in the 2 ×2 case exactly 50% of the matrices are discarded (Ardiyansyah et al. 2021, Table 5), while in the case of 4 ×4 matrices up to 98.26545% of the matrices are discarded (see Table 7) and in the case of 6 ×6 matrices the amount of discarded matrices is about 99.99863% as indicated in Table 9. However, when restricting to subsets of Markov matrices which are mathematically more meaningful in biological terms, such as DD or DLC matrices, the proportion of embeddable matrices is much higher so that we are discarding less matrices (e.g. for DD we discard 68.41679% of 4×4 matrices and 97.2441% of 6 ×6 matrices). This is not to say that it makes no sense to use continuous-time models but to highlight that one should take the above restrictions into consideration when working with these models. Conversely, when working with the whole set of Markov matrices one has to be aware that they might end up considering lots of non-meaningful matrices. Acknowledgements Dimitra Kosta was partially supported by a Royal Society Dorothy Hodgkin Research Fellowship DHF\R1\201246. Jordi Roca-Lacostena was partially funded by Secretaria d’Universitats i Recerca de la Generalitat de Catalunya (AGAUR 2018FI_B_0094). Muhammad Ardiyansyah is partially supported by the Academy of Finland Grant No. 323416. Funding Open Access funding provided by Aalto University. Declarations Conflict of interest The authors declare no conflict of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 123 69 Page 36 of 37 M. Ardiyansyah et al. References Aitken AC (2017) Determinants and matrices. Read Books Ltd Ardiyansyah M, Kosta D, Kubjas K (2021) The model-specific Markov embedding problem for symmetric group-based models. J Math Biol 83(3):1–26 Baake M, Sumner J (2020) Notes on markov embedding. Linear Algebra Appl 594:262–299 Bain J, Switzer C, Chamberlin R, Benner SA (1992) Ribosome-mediated incorporation of a non-standard amino acid into a peptide through expansion of the genetic code. Nature 356(6369):537–539 Benner SA, Sismour AM (2005) Synthetic Biol. Nat Rev Genet 6(7):533–543 Cantoni A, Butler P (1976) Eigenvalues and eigenvectors of symmetric centrosymmetric matrices. Linear Algebra Appl 13(3):275–288 Carette P (1995) Characterizations of embeddable 3×3 stochastic matrices with a negative eigenvalue. New York J Math 1:120–129 Casanellas M, Kedzierska AM (2013) Generating Markov evolutionary matrices for a given branch length. Linear Algebra Appl 438(5):2484–2499 Casanellas M, Sullivant S (2005) The strand symmetric model. In: Pachter L, Sturmfels B (eds) Algebraic statistics for computational biology. Cambridge University Press, New York Casanellas M, Fernández-Sánchez J, Roca-Lacostena J (2020) Embeddability and rate identifiability of Kimura 2-parameter matrices. J Math Biol 80(4):995–1019 Casanellas M, Fernández-Sánchez J, Roca-Lacostena J (2020b) An open set of 4×4 embeddable matrices whose principal logarithm is not a markov generator. Linear and Multilinear Algebra pp 1–12 Casanellas M, Fernández-Sánchez J, Roca-Lacostena J (2023) The embedding problem for Markov matrices. Publicacions Matemàtiques 67:411–445 Chang JT (1996) Full reconstruction of Markov models on evolutionary trees: identifiability and consistency. Math Biosci 137(1):51–73 Chen Y, Chen J (2011) On the imbedding problem for three-state time homogeneous Markov chains with coinciding negative eigenvalues. J Theor Probab 24:928–938 Culver WJ (1966) On the existence and uniqueness of the real logarithm of a matrix. Proceed Am Math Soc 17:1146–1151 Cuthbert JR (1972) On uniqueness of the logarithm for Markov semi-groups. J Lond Math Soc 2(4):623–630 Davies E et al. (2010) Embeddable Markov matrices. Electron J Probab 15:1474–1486 Elfving G (1937) Zur theorie der Markoffschen ketten. Acta Societatis Scientiarum FennicæNova Series A 2(8):17 pages Fuglede B (1988) On the imbedding problem for stochastic and doubly stochastic matrices. Probab Theory Relat Fields 80:241–260 Gawrilow E, Joswig M (2000) Polymake: a framework for analyzing convex polytopes. In: Polytopescombinatorics and computation, Springer, pp 43–73 Goodman GS (1970) An intrinsic time for non-stationary finite Markov chains. Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete 16:165–180 Hammersley J (2013) Monte carlo methods. Springer Science & Business Media Higham NJ (2008) Functions of matrices: Theory and computation, vol 104. SIAM Hoshika S, Leal NA, Kim MJ, Kim MS, Karalkar NB, Kim HJ, Bates AM, Watkins NE, SantaLucia HA, Meyer AJ et al. (2019) Hachimoji DNA and RNA: a genetic system with eight building blocks. Science 363(6429):884–887 Inc WR (2022) Mathematica, Version 13.1 https://www.wolfram.com/mathematica Iosifescu M (2014) Finite Markov processes and their applications. Courier Corporation James CR (1973) The logarithm function for finite-state Markov semi-groups. J Lond Math Soc 2(3):524– 532 Jia C (2016) A solution to the reversible embedding problem for finite markov chains. Statist Probab Lett 116:122–130 Johansen S (1974) Some results on the imbedding problem for finite Markov chains. J London Math Soc s2-8(2):345–351 Kimura M (1957) Some problems of stochastic processes in genetics. Ann Math Statistics pp 882–901 Kingman JFC (1962) The imbedding problem for finite Markov chains. Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete 1(1):14–24 Leal NA, Kim HJ, Hoshika S, Kim MJ, Carrigan MA, Benner SA (2015) Transcription, reverse transcription, and analysis of RNA containing artificial genetic components. ACS Synth Biol 4(4):407–413 123 Embeddability of centrosymmetric matrices capturing... Page 37 of 37 69 Malyshev DA, Dhami K, Quach HT, Lavergne T, Ordoukhanian P, Torkamani A, Romesberg FE (2012) Efficient and sequence-independent replication of DNA containing a third base pair establishes a functional six-letter genetic alphabet. Proceedings of the National Academy of Sciences 109(30):12,005–12,010 Pachter L, Sturmfels B (2005) Algebraic statistics for computational biology, vol 13. Cambridge University Press, UK Roca-Lacostena J (2021) The embedding problem for Markov matrices. PhD thesis, Universitat Politècnica de Catalunya Roca-Lacostena J, Fernández-Sánchez J (2018) Embeddability of Kimura 3ST Markov matrices. J Theor Biol 445:128–135 Roca-Lacostena J, Fernández-Sánchez J (2018) Embeddability of Kimura 3st markov matrices. J Theor Biol 445:128–135 Runnenberg JT (1962) On Elfving’s problem of imbedding a time-discrete markov chain in a timecontinuous one for finitely many states. Proceedings of the KNAW - Series A, Mathematical Sciences 65:536–541 Schensted IV (1958) Appendix model of subnuclear segregation in the macronucleus of ciliates. Am Nat 92(864):161–170 Sismour AM, Lutz S, Park JH, Lutz MJ, Boyer PL, Hughes SH, Benner SA (2004) Pcr amplification of DNA containing non-standard base pairs by variants of reverse transcriptase from human immunodeficiency virus-1. Nucleic Acids Res 32(2):728–735 Stein P (1966) A note on the volume of a simplex. Am Math Mon 73(3):299–301 Weaver JR (1985) Centrosymmetric (cross-symmetric) matrices, their basic properties, eigenvalues, and eigenvectors. Am Math Mon 92(10):711–717 Yang Z, Hutter D, Sheng P, Sismour A, Benner S (2006) Artificially expanded genetic information system: a new base pair with an alternative hydrogen bonding pattern. Nucleic Acids Res 34(21):6095–101 Yang Z, Sismour AM, Sheng P, Puskar NL, Benner SA (2007) Enzymatic incorporation of a third nucleobase pair. Nucleic Acids Res 35(13):4238–4249 Yang Z, Chen F, Alvarado JB, Benner SA (2011) Amplification, mutation, and sequencing of a six-letter synthetic genetic system. J Am Chem Soc 133(38):15,105–15,112 Yap VB, Pachter L (2004) Identification of evolutionary hotspots in the rodent genomes. Genome Res 14(4):574–579 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123