Markov-modulated continuous-time Markov chains to identify site- and branch-specific evolutionary variation

TL;DR

This paper introduces Markov-modulated models (MMMs) for phylogenetic inference that account for time-varying evolutionary processes at individual sites, improving model fit and tree estimation accuracy.

Contribution

It develops a flexible MMM framework in BEAST, enabling modeling of site-specific evolutionary variation over time, which was previously overlooked.

Findings

MMMs improve phylogenetic tree estimation accuracy.

Model fit is substantially better with MMMs compared to standard models.

Simulations confirm accurate model identification without bias toward complex models.

Abstract

Markov models of character substitution on phylogenies form the foundation of phylogenetic inference frameworks. Early models made the simplifying assumption that the substitution process is homogeneous over time and across sites in the molecular sequence alignment. While standard practice adopts extensions that accommodate heterogeneity of substitution rates across sites, heterogeneity in the process over time in a site-specific manner remains frequently overlooked. This is problematic, as evolutionary processes that act at the molecular level are highly variable, subjecting different sites to different selective constraints over time, impacting their substitution behaviour. We propose incorporating time variability through Markov-modulated models (MMMs) that allow the substitution process (including relative character exchange rates as well as the overall substitution rate) that models the evolution at an individual site to vary across lineages. We implement a general MMM framework in BEAST, a popular Bayesian phylogenetic inference software package, allowing researchers to compose a wide range of MMMs through flexible XML specification. Using examples from bacterial, viral and plastid genome evolution, we show that MMMs impact phylogenetic tree estimation and can substantially improve model fit compared to standard substitution models. Through simulations, we show that marginal likelihood estimation accurately identifies the generative model and does not systematically prefer the more parameter-rich MMMs. In order to mitigate the increased computational demands associated with MMMs, our implementation exploits recently developed updates to BEAGLE, a high-performance computational library for phylogenetic inference.

Equations

$$P(X(t+b)=j\mid X(b)=i,\{X(u):0\le u<b\})=P(X(t+b)=j\mid X(b)=i).$$

$$P(X(t+b)=j\mid X(b)=i)=P(X(t)=j\mid X(0)=i),$$

$$\mathbf{P}(t)=\exp (\mathbf{Q}t),$$

$$\mathbf{\Phi}_{GTR}=\delta\left(\begin{array}[]{ccccc}-&\eta_{2}R_{1\leftrightarrow 2}&\cdots&\eta_{K}R_{1\leftrightarrow K}\\ \eta_{1}R_{1\leftrightarrow 2}&-&&\vdots\\ \vdots&&\ddots&\eta_{K-1}R_{K-1\leftrightarrow K}\\ \eta_{1}R_{1\leftrightarrow K}&\cdots&\eta_{K-1}R_{K-1\leftrightarrow K}&-\end{array}\right),$$

$$\mathbf{\Phi}=\left(\begin{array}[]{ccccc}-\sum_{k\neq 1}\phi_{1k}&\phi_{12}&\cdots&\phi_{1K}\\ \phi_{21}&-\sum_{k\neq 2}\phi_{2k}&&\vdots\\ \vdots&&\ddots&\phi_{K-1,K}\\ \phi_{K1}&\cdots&\phi_{K,K-1}&-\sum_{k\neq K}\phi_{Kk}\end{array}\right),$$

$$\mathbf{\Lambda }$$

$$=\text{diag}\left({\rho }_1{\mathbf{Q}}_1,\dots ,{\rho }_K{\mathbf{Q}}_K\right)+\mathbf{\Phi }\otimes \mathbf{I},$$

$${\rho }_1=\cdots ={\rho }_K=1.$$

$${\mathbf{Q}}_k=\mathbf{Q}\text{ for all }k.$$

$${\mathbf{Q}}_k$$

$${\varphi }_{k{k}^{\prime }}$$

$${\varphi }_{k{k}^{\prime }}=0\text{ for all }|k-k^{\prime }|>1,$$

$$\phi_{kk^{\prime}}=\left\{\begin{array}[]{l}\phi\text{ if }|k-k^{\prime}|=1,\\ 0\text{ otherwise. }\end{array}\right.$$

$$L(\mathbf{\theta },\mathbf{\varphi },\mathcal{T},\mathcal{M}\mid \mathbf{Y})=\prod _i\left(\sum _{(k,s)}{\pi }_k{\pi }_s{L}_i^R((k,s),\mathbf{\theta },\mathbf{\varphi },\mathcal{T},\mathcal{M}\mid {\mathbf{Y}}_i)\right),$$

$$Q_{ij}^{(k)}={\pi }_{kj}{R}_{i\leftrightarrow j}^{(k)}.$$

$${\mathbf{Q}}_k=\mathbf{Q},k=1,\dots ,4.$$

Full text

Markov-modulated continuous-time Markov chains

to identify site- and branch-specific

evolutionary variation

Guy Baele1, Mandev S. Gill1, Philippe Lemey1,

and Marc A. Suchard2,3,4

1*Department of Microbiology, Immunology and Transplantation, Rega Institute, KU Leuven, Herestraat 49, 3000 Leuven, Belgium

2Department of Biostatistics, Jonathan and Karin Fielding School of Public Health, University of California, Los Angeles, United States*

3*Department of Biomathematics, David Geffen School of Medicine at UCLA, University of California, Los Angeles, United States

4Department of Human Genetics, David Geffen School of Medicine at UCLA, Universtiy of California, Los Angeles, United States*

Corresponding author: Guy Baele, Department of Microbiology, Immunology and Transplantation, Rega Institute, KU Leuven, Herestraat 49, 3000 Leuven, Belgium; E-mail: guy.baele@kuleuven.be

Abstract

Markov models of character substitution on phylogenies form the foundation of phylogenetic inference frameworks. Early models made the simplifying assumption that the substitution process is homogeneous over time and across sites in the molecular sequence alignment. While standard practice adopts extensions that accommodate heterogeneity of substitution rates across sites, heterogeneity in the process over time in a site-specific manner remains frequently overlooked. This is problematic, as evolutionary processes that act at the molecular level are highly variable, subjecting different sites to different selective constraints over time, impacting their substitution behaviour. We propose incorporating time variability through Markov-modulated models (MMMs) that allow the substitution process (including relative character exchange rates as well as the overall substitution rate) that models the evolution at an individual site to vary across lineages. We implement a general MMM framework in BEAST, a popular Bayesian phylogenetic inference software package, allowing researchers to compose a wide range of MMMs through flexible XML specification. Using examples from bacterial, viral and plastid genome evolution, we show that MMMs impact phylogenetic tree estimation and can substantially improve model fit compared to standard substitution models. Through simulations, we show that marginal likelihood estimation accurately identifies the generative model and does not systematically prefer the more parameter-rich MMMs. In order to mitigate the increased computational demands associated with MMMs, our implementation exploits recently developed updates to BEAGLE, a high-performance computational library for phylogenetic inference.

1 Introduction

Molecular sequence evolution is typically modeled by Markov models of character substitution acting along the branches of a phylogenetic tree. These models are phenomenological descriptions of the evolution of DNA as a string of a number of discrete character states, with models of nucleotide substitution among four states being the most widely used in statistical phylogenetics. The Markovian property within such a model reflects the common assumption that evolution has no memory. Further, it is standard to assume that the Markov model is time-homogeneous, so that it can be characterized by a generator or instantaneous rate matrix $\mathbf{Q}$ that remains constant during evolution (Gascuel and Guindon 2007). Starting with the original nucleotide substitution model of Jukes and Cantor (Jukes and Cantor 1969), which does not require any parameter to be estimated from the data, a number of different Markov models of DNA sequence evolution have been proposed over the years, of which the HKY model (Hasegawa et al. 1985) and the general time-reversible model (Lanave et al. 1984; Tavaré 1986) are among the most widely used.

The stochastic process $\{X(t):t\geq 0\}$ over discrete state space $\mathcal{S}$ is called a continuous-time Markov chain (CTMC) if for all $t\geq 0,b\geq 0$ and $i,j\in\mathcal{S}$,

[TABLE]

In other words, a CTMC is a stochastic process having the Markovian property that the conditional distribution of the future $X(t+b)$ given the present $X(b)$ and the past $X(u),0\leq u<b$, depends only on the present and is independent of the past (Tolver 2016). If additionally, $P(X(t+b)=j\mid X(b)=i)$ is independent of $b$, then the CTMC has stationary or time-homogeneous transition probabilities. In this case, we have

[TABLE]

and we denote this transition probability by $P_{ij}(t)$. The transition probabilities for all $i,j\in\mathcal{S}$ can be obtained through computing the matrix exponential

[TABLE]

which is generally accomplished via diagonalization of $\mathbf{Q}$ (Felsenstein 2004) and leads to an $S\times S$ matrix called the transition probability matrix where $S=|\mathcal{S}|$.

Early probabilistic phylogenetic reconstruction methods assumed a single substitution model that acted independently across all sites and lineages. However, the evolutionary processes that act at the molecular level are highly variable, as some sequence positions evolve rapidly while others barely undergo any changes. Different sites are subject to different selective constraints, as exemplified for different codon positions in protein coding genes where third codon positions usually evolve faster than first positions, which, in turn, evolve faster than second positions. Accommodating different rates at different sites has been one of the most influential modeling innovations in molecular phylogenetics (Yang 1996). However, sites may also evolve in a qualitatively different manner, in which case simply accommodating among-site rate variation is insufficient. Pagel and Meade (2004) incorporated such “pattern heterogeneity” through a flexible mixture model that posits that, among a collection of prespecified substitution models, each substitution model applies to any given alignment site with a specific probability, or “weight.” The likelihood for each site is then the weighted sum over these different substitution models. Notably, the Pagel and Meade (2004) mixture model for pattern heterogeneity can be combined with the popular among-site rate variation model of Yang (1994) that assumes site-specific substitution rates drawn from a discretized gamma distribution.

Markov-modulated models (MMMs) naturally extend mixture models to account for time variability, and were introduced into phylogenetics (Tuffley and Steel 1998; Lockhart et al. 1998) to accommodate observations that the evolutionary rate of a particular site in coding sequences can vary across the phylogeny because sites that are critical for the function of a macromolecule may change within the nucleotide sequence over time (Fitch and Markowitz 1970). MMMs allow the process that governs the evolution of an individual site to change over time along a phylogeny by augmenting the usual state space of the substitution model (consisting of sequence characters such as nucleotides or codons) to include categories corresponding to different evolutionary models. These categories also change following a homogeneous, stationary and time-reversible Markovian process, similar to standard CTMCs of character substitution. Gascuel and Guindon (2007) presented a general time-reversible model (GTR) for $K$ categories:

[TABLE]

where each row sums to 0, $\boldsymbol{\eta}=\left(\eta_{1},\ldots,\eta_{K}\right)$ is the stationary distribution of the categories and $\delta$ is an additional parameter that expresses the global rate of changes between all categories. Typically, the $R$ coefficients are normalized such that $\delta$ is the expected number of category changes during one time unit (Gascuel and Guindon 2007). From Equation 4, it follows that an MMM reduces to a mixture model with weights $\boldsymbol{\eta}$ when the rate $\delta$ at which changes between categories occur tends to zero (Guindon et al. 2004; Gascuel and Guindon 2007).

Nearly fifty years ago, Fitch and Markowitz (1970) proposed their “concomitantly variable codons” (or “covarion”) hypothesis, which states that at any given time some sites are invariable due to functional or structural constraints, but that as mutations are fixed elsewhere in the sequence these constraints may change, so that sites that were previously invariable may become variable and vice versa. The covarion hypothesis therefore implies that the pool of variable sites is changing with time. Inspired by these observations, Tuffley and Steel (1998) put forth a simple covarion-style model where a two-state Markov process acts as a switch that turns sites “on” (variable) and “off” (invariable), but does not trigger a concomitant reaction in any other sites. Under their model, substitutions occur at a fixed rate for sites that are turned “on”, while no substitutions occur for sites that are turned “off.” Huelsenbeck (2002) generalized this model by allowing each site to have a fixed rate drawn from a discrete gamma distribution, with sites assuming the fixed rate when turned “on” and assuming a zero rate when turned “off”. Zhou et al. (2010) further extended the Huelsenbeck (2002) model by relaxing the assumption that the two category transition rates (one from “on” to “off”, and the other from “off” to “on”) are constant across sites, employing a Dirichlet process mixture to achieve rate heterogeneity. Galtier (2001) presented an alternative generalization of the Tuffley and Steel (1998) model by assuming that a fixed subset of sites evolves according to a covarion process, with rates drawn from a discrete gamma distribution and switching allowed between different rate classes at a single fixed rate. The remaining sites that do not evolve under a covarion process are posited to have site-specific rates also drawn from a discrete gamma distribution. Wang et al. (2007, 2009) introduced a general covarion model that combines the features of the Tuffley and Steel (1998), Galtier (2001) and Huelsenbeck (2002) models.

These covarion-style models can capture variations in lineage-specific evolutionary rates over time and as such constitute special instances of the Markov-modulated class of Markov processes (Fischer and Meier-Hellstern 1993), whose generator matrices are well described by Kronecker product formalism (Galtier and Jean-Marie 2004). These models suffer from potentially high dimensionality, as the transition matrix has a dimension equal to the product of the number of character states and the number of rate categories, making matrix exponentiation and likelihood calculation computationally challenging. To address these difficulties, Galtier and Jean-Marie (2004) proposed a fast algorithm for diagonalizing the generator matrix. The authors show their algorithm to be tractable even for a large number of rate classes and a large number of states, employing up to 20 rate classes combined with full $61\times 61$ codon models, leading to a $1220\times 1220$ generator matrix of their MMM process.

The development of covarion-style MMMs has inspired the introduction of more general MMMs that permit switching between evolutionary models that differ in relative character exchange rates as well in overall substitution rate. Guindon et al. (2004) introduced a codon-based MMM that allows substitutions between codons to occur under either a purifying, neutral, or positive selection regime and allowing changes (or switches) between selection classes to occur according to a three-state continuous-time Markov process. The authors showed in an analysis of eight HIV-1 env sequence data sets that their model provides a significantly better fit to the data than a model that does not accommodate switches between selection patterns. Gascuel and Guindon (2007) proposed an extension to a general MMM framework for the analysis of nucleotide, amino acid or codon data. Notably, they allow for switching between any evolutionary models, so long as all models share the same stationary distribution of characters. Whelan (2008) proposed a similar general MMM framework for nucleotide substitution that also permits the different substitution models to have different stationary distributions.

In spite of all of these developments, MMMs remain relatively under-utilised due to their prohibitive computational cost and a lack of readily accessible implementations in widely used software packages for phylogenetic inference. There may also exist a general unfamiliarity with how MMMs may improve phylogenetic tree reconstruction. We introduce a Bayesian inference framework for MMMs, implemented in BEAST (Suchard et al. 2018), a software package for Bayesian evolutionary analysis. We strive for optimal generality by allowing switching between evolutionary models that have different substitution rates, relative character exchange rates and stationary distributions. Further, our implementation accommodates phylogenetic uncertainty and is equipped for analysis of any data type. Researchers can easily compose and customize MMMs through flexible XML specification, as well as use an MMM as part of a larger Bayesian phylodynamic analysis by combining it with any of the wide range of models for demographic reconstruction, phylogeographic inference, and phenotypic trait evolution available in BEAST. To cope with the higher computational burden imposed by MMMs, our implementation leverages the power of the BEAGLE library for high-performance likelihood computation (Ayres et al. 2019), and we show considerable performance increases when employing graphics cards aimed at scientific computing. We illustrate the flexibility and utility of our MMM implementation by employing a variety of different model compositions in examples of different complexity from bacterial, plastid, and viral evolution. To the best of our knowledge, no comparably general MMM framework is currently available in a phylogenetic inference software package.

2 Methods

2.1 Markov-modulated model structure

Consider an MMM composed of $K$ evolutionary models (irrespective of those models being nucleotide, amino acid or codon models). Each evolutionary model is defined by a relative substitution rate multiplier $\rho_{k}$ and a substitution model characterized by an instantaneous rate matrix $\mathbf{Q}_{k}=\left\{Q^{(k)}_{ss^{\prime}}\right\}$, of dimension $S$, and stationary distribution $\boldsymbol{\Pi}_{k}=\left(\pi_{k1},\ldots,\pi_{kS}\right)$. We also adopt the usual constraint -$\sum_{s=1}^{S}Q^{(k)}_{ss}\pi_{ks}=1$. The switching process between the $K$ models is defined by a $K$-state continuous-time Markov process with rate matrix

[TABLE]

where the element $\phi_{ij}$ corresponds to the rate of switching from substitution model $i$ to substitution model $j$, and the diagonal elements are fixed such that the rows sum to [math]. We denote the stationary distribution of this switching process by $\boldsymbol{\Psi}=\left(\psi_{1},\ldots,\psi_{K}\right)$. These model switches follow a homogeneous, stationary - but not necessarily time-reversible - Markovian process. In Equation 5, we do not make use of an additional parameter $\delta$ that expresses the global rate of change between the evolutionary models (as in equation 4) because this is a deterministic parameter obtained by normalizing the model-switching process (Guindon et al. 2004; Gascuel and Guindon 2007).

The MMM is characterized by a $KS\times KS$ rate matrix $\mathbf{\Lambda}$(Fischer and Meier-Hellstern 1993):

[TABLE]

where $\mathbf{I}$ is an $S\times S$ identity matrix and $\otimes$ denotes the Kronecker product. The MMM can therefore be considered a single Markov process with a state space equal to the Cartesian product of the state space of the switching process (the evolutionary models) and the state space of the evolutionary models, with cardinality $KS$ and stationary distribution $\boldsymbol{\Pi}_{\mathbf{\Lambda}}=\left(\psi_{1}\pi_{11},\ldots,\psi_{1}\pi_{1S},\ldots,\psi_{K}\pi_{K1},\ldots,\psi_{K}\pi_{KS}\right)$ (Guindon et al. 2004). As noted by Gascuel and Guindon (2007), the MMM in Equation 6 allows for every compound state $(k,s)$ to either: (1) stay in model $k$ and transition to $(k,s^{\prime})$ with rate defined by $\rho_{k}\mathbf{Q}_{k}$, or (2) change evolutionary models and transition to $(k^{\prime},s)$ with rate $\phi_{kk^{\prime}}$. All rows in $\mathbf{\Lambda}$ sum to [math], and because $\boldsymbol{\Psi}\mathbf{\Phi}=0$ and $\boldsymbol{\Pi}_{k}\mathbf{Q}_{k}=0$, it follows that $\boldsymbol{\Pi}_{\mathbf{\Lambda}}\mathbf{\Lambda}=0$.

Standard nucleotide and codon substitution models are often superimposed with a model that employs a discretized gamma distribution to accommodate among-site rate variation (Yang 1994). We will refer to this specific model as the ASRV model. Combining the ASRV model with MMMs can result in a model that is not identifiable. For example, consider a rooted tree of two taxa, $K=2$ with $\mathbf{Q}_{1}=\mathbf{Q}_{2}$ and an ASRV model with two rate categories with rate multipliers $\rho_{1}$ and $\rho_{2}$. Further, assume a switching process with instantaneous rates $\phi_{12}=\phi_{21}=\phi$ and hence equilibrium frequencies $\psi_{1}=\psi_{2}=1/2$. Given this setup, it becomes impossible to distinguish between sites that evolve according to model $\mathbf{Q}_{1}$ (i.e. when the switching process puts them in the first category) with rate $\rho_{2}$ and sites that evolve according to model $\mathbf{Q}_{2}$ with rate $\rho_{1}$.

However, we can still superimpose the ASRV model on an MMM while maintaining identifiability by restricting the MMM relative rate multipliers:

[TABLE]

The resulting model is a generalization of the model proposed by Guindon et al. (2004), but with a more general switching process. Note that the covarion model of Tuffley and Steel (1998) is also nested within the general MMM formulation by setting $K=2$, $\rho_{1}=0$ and $\rho_{2}=1$. The covarion model proposed by Galtier and Jean-Marie (2004) is also nested within our framework by positing one possible substitution model, i.e.

[TABLE]

We note that our proposed MMM can be parameterized as the ubiquitously-used ASRV model, which employs a discretized gamma distribution to determine its relative rates (Yang 1994), solely through setting its parameters as follows:

[TABLE]

and prescribing that $\rho_{k}$ for all $k$ equal the $(k-\frac{1}{2})/K\times 100^{\text{\tiny\,th}}$ percentile of a gamma distribution with fixed expectation 1 and unknown shape $\alpha$.

Other common restrictions on the switching process, as extensively used in telecommunication and teletraffic networks (Fischer and Meier-Hellstern 1993), include an ordered model random walk:

[TABLE]

and an ordered model random walk with single rate:

[TABLE]

2.2 Likelihood

In this section, we adopt a similar notation to Gascuel and Guindon (2007) to describe the data likelihood under an MMM. Likelihood calculations for MMMs employ a standard pruning approach (Felsenstein 1981), with integration over the compound states (i.e. the evolutionary model and character state) at the internal nodes of the tree, and integration over the unobserved categories at the tips. Let $\mathbf{Y}=(\mathbf{Y}_{1},\ldots,\mathbf{Y}_{L})$ where $\mathbf{Y}_{\ell}$ are the extant characters observed at aligned site $\ell$ for $\ell=1,\ldots,L$, and let ${\cal T}$ denote the phylogenetic tree with its branch lengths. Let ${\cal M(\boldsymbol{\theta},\boldsymbol{\phi})}$ denote the MMM that models the evolutionary process for all sites, where $\boldsymbol{\theta}=\{\boldsymbol{\theta}_{1},\ldots,\boldsymbol{\theta}_{K}\}$ and $\boldsymbol{\theta}_{k}$ represents parameters for the $k$th evolutionary model, and $\boldsymbol{\phi}$ parameters of the switching process. The observed data likelihood is:

[TABLE]

where the product is taken over every site $i$ in the alignment, with each site assumed to evolve independently. The sum over the compound states $(k,s)$ replaces the sum over the nucleotide characters $s$ that is performed for standard nucleotide substitution models (Gascuel and Guindon 2007). Here, $L_{i}^{R}((k,s),\boldsymbol{\theta},\boldsymbol{\phi},{\cal T},{\cal M}\mid\mathbf{Y})$ is the likelihood of the data at site $i$ under category $k$ and given that state $s$ is observed at site $i$ of the root node $R$. We can generalize this notation as $L_{i}^{v}((k,s),\boldsymbol{\theta},\boldsymbol{\phi},{\cal T},{\cal M}\mid\mathbf{Y})$ for node $v$ to express the partial likelihood of observing the characters at site $i$ in the extant sequences descending from $v$. This notation can be shortened to $L_{i}^{v}(k,s)$ because $\boldsymbol{\theta}$, $\boldsymbol{\phi}$, ${\cal T}$, ${\cal M}$ and $\mathbf{Y}$ are the same for all sites and nodes. Let $l$ and $r$ be the left and right descendants of $v$ and $t_{v}$ the length of the branch connecting $v$ to its parent. Each partial likelihood is then defined as follows (taking into account that the evolutionary categories are unobservable; Gascuel and Guindon (2007)):

L_{i}^{v}(k,s)=\left\{\begin{tabular}[]{l}1 if v$is a leaf with nucleotide character$s$\\ 0 if$v$is a leaf with nucleotide character different from$s$\\$\left(\sum_{(k^{\prime},s^{\prime})}P_{(k,s)(k^{\prime},s^{\prime})}(t_{l})L_{i}^{l}(k^{\prime},s^{\prime})\right)\left(\sum_{(k^{\prime},s^{\prime})}P_{(k,s)(k^{\prime},s^{\prime})}(t_{r})L_{i}^{r}(k^{\prime},s^{\prime})\right) otherwise.\end{tabular}\right.

The substitution probabilities $P_{(k,s)(k^{\prime},s^{\prime})}(t)$ are computed using Equation 3, which requires a diagonalisation of $\mathbf{\Lambda}$ that has computational complexity ${\cal O}\left(K^{3}S^{3}\right)$ (Pan and Chen 1999), although lower complexity may be achieved depending on the Kronecker structure for $\boldsymbol{\Lambda}$ (but see the Discussion section). Computing these probabilities for all ${\cal O}\left(N\right)$ branches in the phylogeny therefore has a complexity of ${\cal O}\left(NK^{3}S^{3}\right)$. Evaluating the $L$ site likelihoods through the tree-pruning (or peeling) algorithm (Felsenstein 1981) amounts to a complexity of ${\cal O}\left(NLK^{2}S^{2}\right)$. Taken together, with a relatively small cost ${\cal O}\left(L\right)$ for taking logarithm of site likelihoods and summing over sites results in a computational complexity of ${\cal O}\left(NK^{3}S^{3}+NLK^{2}S^{2}\right)$ for the log-likelihood of the observed data. We confront these computationally demanding tasks by heavily parallelizing their computation using powerful multi- and many-core hardware through the BEAGLE library (Ayres et al. 2019).

2.3 Models and priors

We use the HKY (Hasegawa et al. 1985) and GTR (Lanave et al. 1984; Tavaré 1986) nucleotide substitution models in our examples. We adopt a Yule pure birth process (Yule 1924) tree prior on the speciation process for the isochronous data sets, with a normally distributed prior on the log-transformed birth rate with mean 1.0 and standard deviation 1.5. We assume a constant population size coalescent tree prior for the heterochronous influenza data set, with a normal prior on the log-transformed population size with mean 1.0 and standard deviation 1.5. For the HKY model, we assume a normally distributed prior on the log-transformed transition-transversion ratio with mean 1.0 and standard deviation 1.25. For the GTR model, we assume diffuse gamma distributed priors on its relative rate parameters. In particular, we employ a gamma distribution with shape 0.05 and scale 20.0 for $r_{AG}$ and gamma distribution with shape 0.05 and scale 10.0 for $r_{AC}$, $r_{AT}$, $r_{CG}$ and $r_{GT}$. We assume an exponential prior with mean 0.5 on the rate heterogeneity parameter of the ASRV model (Yang 1994) as well as on the parameter describing the discretized gamma distribution to generate relative rate multipliers $\rho_{i}$. Additionally, we assume a CTMC reference prior on the strict molecular clock rate in the heterochronous data set (Ferreira and Suchard 2008) and an exponential prior with mean 1/3 on the standard deviation of the relaxed clock model with an underlying lognormal distribution for the isochronous data sets. Finally, we assume an exponential prior with mean 1.0 on each of the switching rate parameters $\phi_{ij}$ and a Dirichlet prior with concentration parameters equal to 1.0 on each set of nucleotide base frequencies. For all analyses performed in the results section, the effective sample size (ESS) values for all estimated parameters were well above 200, as identified using Tracer (Rambaut et al. 2018).

2.4 Bayesian model comparison

Bayesian phylogenetic inference focuses on the joint posterior probability distribution of the tree and all evolutionary model parameters conditional on the observed sequence data. Employing Markov chain Monte Carlo (MCMC) methods to approximate this posterior distribution handily avoids computing the marginal likelihood $p(\mathbf{Y}\mid{\cal M})$ that is notoriously difficult to estimate. In order to compare the relative performance of our proposed MMMs to standard nucleotide substitution models for phylogenetic inference, however, estimating the marginal likelihood remains essential. In all of our examples, we employ generalized stepping-stone sampling (GSS; Fan et al. (2011); Baele et al. (2016)), an estimation procedure that accommodates phylogenetic uncertainty and provides accurate marginal likelihood estimates in a feasible amount of time. GSS requires working distributions (Fan et al. 2011; Baele et al. 2016) for all parameters based on the posterior samples collected when exploring the preceding Markov chain. BEAST (Suchard et al. 2018) constructs these distributions using kernel density estimation (KDE) based on Gaussian kernels and Silverman’s rule to determine the bandwidth (Silverman 1986). Using MMMs will often require estimating a large collection of parameters, which may lead to poor ESS values for those parameters and also a more challenging KDE. In order to avoid manually inspecting every MCMC analysis prior to starting the marginal likelihood estimation process, we suggest using the prior (distribution) as the working distribution for the MMM parameters as a general approach in conjunction with the specific tree working priors proposed in Baele et al. (2016).

3 Results

We here consider substitution models that are time-reversible and therefore substitution model $k$ will have instantaneous rates $Q^{(k)}_{ij}$ that can be expressed in terms of base frequencies $\pi_{kj}$ and symmetric rate parameters $R^{(k)}_{i\leftrightarrow j}=R^{(k)}_{j\leftrightarrow i}$ as follows:

[TABLE]

Thus a substitution model can be specified in terms of its base frequencies and symmetric rate parameters $\textbf{R}_{k}=\{R^{(k)}_{i\leftrightarrow j}|i\neq j,(i,j)\in\mathcal{S}\}$.

We adopt the following notation: MMM($M$)ijk, where $M$ denotes the type of substitution model and $i$, $j$, and $k$ denote the numbers of distinct sets of symmetric rate parameters, sets of base frequencies, and relative rate multipliers, respectively. For example, an MMM(HKY)222 refers to an MMM featuring two different HKY substitution models, each with its own set of symmetric rate parameters and set of base frequencies, and two different relative rate multipliers. An MMM(HKY)122 refers to an MMM featuring two different HKY substitution models that share the same set of symmetric rate parameters but have different sets of base frequencies, along with two different relative rate multipliers. When the relative rate multipliers are all fixed to 1 to superimpose an ASRV model (Yang 1994, 1996), the $k$ subscript is omitted (for example, MMM(HKY)22).

3.1 Bacterial 16S ribosomal RNA

Differences in base composition throughout the genome can bias phylogenetic inference when not properly taken into account. Often, the proportion of A+T in a genome differs from that of G+C, and different organisms exhibit different patterns of base composition. At the level of the entire genome, GC content varies greatly within and among major groups of organisms, which can skew phylogenetic reconstruction if not properly unaccounted for (Mooers and Holmes 2000). Two different evolutionary processes have been singled out as possible explanations for varying patterns of base composition: biases in the underlying process of mutation, as similar levels of GC content are often found in regions with different functional constraints, and natural selection, with increased global GC content in bacteria possibly being selected for by UV exposure (Singer and Ames 1970).

Environmental variation shaping nucleotide composition may cause unrelated taxa to share similar base composition and therefore be grouped together within a clade. In order to accurately reconstruct evolutionary histories through phylogenetic inference, these potentially differing base compositions need to be accommodated in an explicit manner by the nucleotide substitution model. To address this, Blanquart and Lartillot (2006) developed a non-stationary and non-homogeneous model accounting for compositional biases, allowing the composition to change at random points in the tree, with the total number of change points across the tree being inferred from the data. Through a Bayesian analysis of eubacterial 16S rRNA and BAS1 gene yeast data sets, the authors show that in most cases, the stationarity assumption was rejected in favor of their non-stationary model.

We evaluate our MMM framework on 16S ribosomal RNA of five bacterial sequences: Deinococcus radiodurans, Thermus thermophilus, Thermotoga maritima, Aquifex pyrophilus and Bacillus subtilis (GenBank accession numbers: Y11332.1, AJ251939.1, NR_029163.1, M83548.2 and CP009796.1). We use standard nucleotide substitution models as well as MMMs to infer their evolutionary history while fixing the Aquifex pyrophilus sequence as an outgroup. The true tree topology of this eubacterial data set is believed to group D. radiodurans and T. thermophilus together to the exclusion of B. subtilis, T. maritima, and A. pyrophilus, given that D. radiodurans and T. thermophilus share the same peptidoglycan and menaquinone type (Murray 1992). However, phylogenetic reconstruction under stationary models has a tendency to erroneously group D. radiodurans and B. subtilis together, because these mesophiles have similar, relatively low GC content.

Figure 1 shows the results of these phylogenetic reconstructions, with the HKY and GTR models – both featuring an ASRV model and a relaxed molecular clock with an underlying lognormal distribution – yielding similar (log) marginal likelihoods. Note that, because we will include an ASRV model in all of these MMMs, we set all $\rho_{i}$ in equation 6 to 1 to ensure identifiability. Both the HKY and GTR models express strong support in favor of a clustering of D. radiodurans and B. subtilis (see Figure 1), with the GTR model yielding a small increase in model fit to the data over the HKY model (log BF $<$ 1). As such, both models yield an incorrect clustering, which appears to be primarily based on both sequences being mesophilic (low GC content), whereas the three other sequences are considered thermophilic (high GC content). While an MMM of the type introduced by Tuffley and Steel (1998) offers no improvement over these models when $\mathbf{Q}_{1}$ is parameterised as an HKY model, a significant improvement in model fit can be obtained when $\mathbf{Q}_{1}$ is parameterised as a GTR model (log BF = 19). However, any MMM with two sets of base frequencies and with either a single set of symmetric rate parameters (an MMM($*$)12) or with two different sets of symmetric rate parameters (an MMM($*$)22) offers a further improvement in model fit compared to the standard nucleotide substitution models tested (25 $<$ log BF $<$ 35). This can be attributed to the fact that MMMs are able to accommodate differing base compositions throughout the tree topology, and consequently yield an accurate phylogenetic reconstruction of the bacterial relationships, with the D. radiodurans and T. thermophilus clustering together (see Figure 1) (Embley et al. 1993; Mooers and Holmes 2000).

The base frequency estimates for the CTMC models within the MMM reflect the presence of mesophilic sequences (low GC content; orange in Figure 1) and thermophilic sequences (high GC content; blue in Figure 1) in our data, illustrating the capabilities of these MMMs. Despite the fact that only eight branches connect the observed sequences, alignment sites switch up to four times between CTMC models across the phylogeny, indicating evolutionary patterns that cannot possibly be accommodated using standard nucleotide substitution models. Over forty percent of the alignment sites undergoing at least one switch between CTMC models in a highly asymmetric manner (see Figure 1), with the models also exhibiting drastically differing transition/transversion ratios.

In conclusion, we show that appropriately modeling compositional heterogeneity for these eubacterial sequences enables inference of the correct phylogeny as well as base frequency compositions that reflect the presence of both mesophilic and thermophilic sequences in the data set. We note that similar observations regarding phylogenetic inference under compositional heterogeneity have been made in previous studies. For example, Foster (2004) used a heterogeneous model that accommodates compositional differences over the tree, without requiring that each branch is equipped with its own independent base composition, allowing modeling of compositional heterogeneity with few additional parameters. This reflects the assumption that compositional differences do not vary continuously over the tree. The compositional heterogeneity model of Foster (2004) yields the correct topology, as does the non-stationary and non-homogeneous model of Blanquart and Lartillot (2006), while standard nucleotide substitution models fail to do so.

3.2 Plant plastid genes

We consider nucleotide sequence data from the protein-coding genes of 23 completely sequenced plant plastid genomes, previously analysed by Ané et al. (2005) to measure the independence of the substitution process between two groups of taxa as a means of detecting covarion evolution. Assuming a fixed underlying reference tree that represents the likely relationships of plant taxa for which complete chloroplast sequences were available at the time, the covarion test of Ané et al. (2005) detected significant covarion evolution ($P<0.0005$) in 14 of 57 genes analysed across all positions. We here analyse the psaB and ndhD genes with standard nucleotide substitution models and MMMs and compare the inferred phylogenies and model fit.

A comparison of standard nucleotide substitution models reveals that the combination of a GTR model and an ASRV model, along with a relaxed clock assuming an underlying lognormal distribution, yields the highest (log) marginal likelihood for both data sets. We conduct analyses with MMMs that feature an HKY or GTR substitution model with a single set of symmetric rate parameters along with two or three different sets of base frequencies (i.e. MMM($*$)12 and MMM($*$)13 models), as well as generalizations of these MMMs that feature as many different sets of symmetric rate parameters as sets of base frequencies (i.e. MMM($*$)22 and MMM($*$)33 models). For all of these models, we set all $\rho_{i}$ in equation 6 to 1 to ensure identifiability when using an ASRV model in combination with MMMs. We also analyse the data with a nucleotide covarion model (Tuffley and Steel 1998), which we can easily compose within our MMM framework through XML specification.

For the psaB data set, the covarion-style model is strongly preferred over a standard GTR+ASRV substitution model by a log Bayes factor of 208. The MMM(GTR)12 and MMM(GTR)22 yield log Bayes factors of 63 and 316, respectively, over the standard GTR+ASRV model. MMM(GTR)13 and MMM(GTR)33 parameterisations yield further increases in model fit of 327 and 356, respectively, over the GTR+ASRV model. Because additional categories within the MMM offer diminishing returns in terms of model fit at the expense of additional computation time, we did not explore MMMs with even higher dimensions. Figure 2 shows the maximum clade credibility (MCC) trees obtained under the standard GTR+ASRV model and the MMM(GTR)33 that generated the highest (log) marginal likelihood. While the clustering within the seed plants is identical under both models, substantial differences in posterior support can be observed for specific clades. In the remaining part of the tree, these models result in completely different clustering patterns with strong support for many clades under the MMM(GTR)33 model.

In Figure 3, we illustrate the complex substitution patterns across all sites on the maximum a posteriori (MAP) *psaB *phylogeny, using the most probable hidden state realisations of the MMM(GTR)33. We use a simple counting procedure to reconstruct which sites evolve according to which CTMC within the MMM(GTR)33, and we observe a relatively small amount of CTMC switching throughout the phylogeny. The reconstructed patterns go beyond mere codon position partitioning, as we observe different substitution dynamics per codon position. In particular, the third codon position is the only position that evolves according to a particular CTMC a majority of the time, and it also exhibits the greatest degree of switching between CTMC realisations. We also depict the mean posterior instantaneous substitution rates of the various MMM components in Figure 3, showing a clearly asymmetric CTMC switching process and three distinct GTR model realisations within the MMM. This complex interplay of model components is consistent with the strong Bayes factor support of the MMM(GTR)33 over all other models tested.

For the ndhD data set, the covarion-style model is strongly preferred over a standard GTR+ASRV substitution model, yielding a log Bayes factor of 84 in favor of the former. However, the MMM(GTR)12 and MMM(GTR)22 yield a log Bayes factor of respectively 138 and 155 over the standard GTR+G model. Again, further increases in model fit can again be obtained using the MMM(GTR)13 and MMM(GTR)33 parameterisations, of respectively 171 and 184 over the GTR+ASRV model. We deemed the performance increase of the MMM(GTR)33 over the MMM(GTR)22 insufficient to warrant testing even higher-dimensional MMMs. Figure 4 shows the inferred phylogenies for the standard GTR+ASRV model and the MMM(GTR)33 that generated the highest (log) marginal likelihood. The clustering within the angiosperms is identical under both models, but substantial differences can again be observed in the posterior support of particular clades. Basal to the angiosperm clade, the clustering of Physcomitrella patens with Marchantia polymorpha is retrieved under both models but other differences can clearly be observed, often with strong posterior support under the MMM(GTR)33.

In Figure 5, we illustrate the complex substitution patterns across all sites on the maximum a posteriori (MAP) *ndhD *phylogeny, using the most probable hidden state realisations of the MMM(GTR)33. We again use a simple counting procedure to reconstruct per-site evolutionary patterns within the MMM(GTR)33, and we see that a more substantial amount of CTMC switching occurs throughout the phylogeny for all positions than for the *psaB *data set. The third codon position exhibits the greatest degree of switching between CTMC realisations, as with the *psaB *data set. On the other hand, the first codon position in this *ndhD *data set evolves according to a predominant CTMC, rather than the third codon position, as in the *psaB *data set. Given the substantial amount of CTMC switching across the *ndhD *phylogeny, models that assume a fixed CTMC per position cannot accurately describe the complex substitution dynamics we observe here, explaining the large Bayes factors in favor of MMMs over traditional substitution models. The mean posterior instantaneous substitution rates of the various MMM components show a less asymmetric CTMC switching process than for the *psaB *data set, with again three clearly distinct GTR model realisations within the MMM.

In conclusion, we observe that MMMs infer substantially different phylogenies than standard nucleotide substitution models, and the differences yielded by MMMs are supported by significant increases in model fit. In order to show that these large differences are not artifacts of using such high-dimensional models, we perform a small simulation study to assess the ability of MMMs to retrieve the generative models of simulated sequence alignments and to quantify their increase in model fit when the MMM was the generative model (Appendix A). Our simulation study shows similar differences in model fit to the ones obtained in this section for the psaB and ndhD genes, as well as the ability of state-of-the-art Bayesian model selection to select the generative substitution model even when compared with similar model parameterisations. Importantly, when simulating data under a standard GTR model, MMMs exhibit a worse model fit than the GTR model that generated the data. These analyses of simulated data show that MMMs can easily be used in combination with recent developments in Bayesian model selection, and provide additional support for our conclusions that these models can yield substantial increases in model fit over standard nucleotide substitution models.

We note that each additional CTMC within an MMM leads to (much) increased computational demands, and that a search for the optimal MMM may hence prove time-consuming for complex large data sets. In order to make such computations manageable, BEAST can however exploit the BEAGLE library (Ayres et al. 2019) to offload the large matrix exponentiations and multiplications onto powerful multi-core hardware solutions. In particular, the use of graphics cards for scientific computing yields significant performance gains over standard multi-core processors, rendering phylogenetic inference under these MMMs feasible despite their complexity (Appendix B).

3.3 Assessing heterotachy using Markov-modulated models

Phylogenetic inference using molecular clock models has proven to be crucial for accurate reconstruction of the origin and spread influenza A virus (IAV) within and between hosts (Rambaut et al. 2008). Recently, Worobey et al. (2014) showed that estimates under the popular uncorrelated relaxed clock model can be flawed when discrete rate variation exists among clades. In a similar fashion, traditional substitution models may not be able to fully capture complex substitution dynamics over time that lead to non-stationary base distributions. We here analyse a sequence alignment comprising the external antigenic haemagglutinin (HA) gene (Worobey et al. 2014), encompassing 454 taxa and spanning 1815 base pairs, providing for a computationally demanding analysis using MMMs, to test if the evolutionary rate of sites varies over time, a process known as heterotachy (Lopez et al. 2002).

We consider an MMM(HKY)114 for this analysis, as we specifically focus on the patterns of evolutionary rate variation across the phylogeny, i.e. we impose

[TABLE]

We draw the relative rate multipliers ($\rho_{k},\text{with }k=1,\ldots,4$) from a discretized gamma distribution $\Gamma(\alpha)$, with its shape $\alpha$ estimated from the data and its scale $\beta=\alpha$. Note that this discretization is identical to how among-site rate variation is typically modelled (Yang 1994, 1996). We compare the performance of this MMM(HKY)114 to that of a standard HKY substitution model (with and without an ASRV model). We also assess the performance of partitioning according to codon position for each of the models tested.

Figure 6 shows the relative performance of the different models tested on this data set along with their clock rate and root height estimates (95% credibility intervals are shown in brackets). As expected, accommodating ASRV substantially improves model fit for a standard nucleotide substitution model (with and without codon position partitioning), and also has a significant impact on root height estimation. Partitioning according to codon position consistently improves model fit as well, as this has clear biological motivation (Shapiro et al. 2006). Modeling heterotachy using an MMM(HKY)114 yields further substantial improvements in model fit over ASRV models, both with and without codon position partitioning (log Bayes factors of 1841 and 2610, respectively), and results in a root height that is much older than the one estimated using a standard nucleotide substitution model without an ASRV model. We did not, however, uncover systematic differences in the clock rate estimation with MMMs compared to standard substitution models (with and without an ASRV model).

We observe heterogeneous rate switching behaviour across codon positions for the top-performing MMM+CP model (Figure 6). This is consistent with the Bayes factor support of MMMs over ASRV models and illustrates that the assumption of constant rates for sites across the underlying phylogeny is an oversimplification. We specifically observe strong asymmetric rate switching at the third codon position, indicative of reduced functional constraints at these positions. These results shows the flexibility of our MMM implementation as well as the potential for MMMs to capture complex processes, such as heterotachy, that cannot be fully captured by popular models that are in frequent use today. Finally, these analyses also show that despite the large number of sequences in this HA data set, Bayesian inference with MMMs – with accompanying (log) marginal likelihood estimation – is feasible on large data sets when employing a high-performance computing library on GPU (Ayres et al. 2019).

4 Discussion

Markov-modulated rate processes were originally developed for information-handling purposes in communication and computer systems, where the arrival rate of information fluctuates over time (Stern and Elwalid 1991). Consisting of a superposition of a number of independent finite state reversible Markov processes, the Markov-modulated Poisson process has seen extensive use in many applications (see e.g. Fischer and Meier-Hellstern (1993)), before finding its way into phylogenetics (see the introduction for an extensive overview). Given the high dimensionality of MMMs, initial applications in phylogenetics focused on efficient approaches for MMM matrix diagonalization and exponentiation (Galtier 2001; Galtier and Jean-Marie 2004) and targeted MMMs with only a few free parameters (Guindon et al. 2004).

Specifically, when an MMM comprises evolutionary models that share the same instantaneous rate matrix (i.e. $\mathbf{Q}_{k}=\textbf{Q}$ for all $k$), the computational complexity of instantaneous rate matrix diagonalization (required for calculation of the data likelihood) is reduced via an efficient Kronecker formulation from ${\cal O}\left(S^{3}K^{3}\right)$ to ${\cal O}\left(S^{3}\right)+{\cal O}\left(SK^{3}\right)$, with $S$ ranging from 4 for nucleotides, over 20 for amino acids, to 61 for codons (Galtier and Jean-Marie 2004). However, this gain in computational efficiency is not applicable in the case of more than one unique $\mathbf{Q}_{k}$, which was the case in all of our examples, with the exception of the heterotachy test for the antigenic haemagglutinin example. For such cases, more general matrix exponentiation via eigenvalue decomposition is required, for which the computational demands increase tremendously with the MMM’s dimension. However, substantial performance improvements in computational hardware now enable the use of MMMs featuring dozens of estimable parameters. In particular, high-performance computational libraries can alleviate the computational burden through offloading the most demanding routines to powerful graphics cards (Ayres et al. 2019). Efficient scaling of higher dimensional MMMs may also be aided by innovations in MCMC sampling. For example, adaptive MCMC transition kernels may allow for more efficient posterior exploration for parameter-rich MMMs, as such kernels have been shown to substantially increase ESS values for standard nucleotide substitution model parameters (Baele et al. 2017).

We have shown that MMMs can offer great improvements in model fit for a wide range of scenarios, from bacterial and viral to plant plastid data sets. Notably, we demonstrated that MMMs can outperform the ASRV model, which has become ubiquitous in phylogenetic reconstruction (Yang 1994, 1996). An MMM is in fact a generalization of the ASRV model, made possible via constructing $\mathbf{\Lambda}$ directly from scaled versions of $\mathbf{Q}_{k}$ (as we have shown in this manuscript) while also allowing for branch rate variation. We have employed the same discretization of an underlying gamma distribution for determining the relative rate multipliers of the $\mathbf{Q}_{k}$, although other approaches may also prove useful towards this end.

A major strength of our framework is its generality. Researchers can compose a wide variety of MMMs through a flexible XML specification for analyses in BEAST. An MMM can comprise an arbitrary number of evolutionary models that can differ in relative substitution rates, relative character exchange rates, and stationary distributions. However, determining the ideal composition of an MMM can be challenging. It may not be obvious how many and which evolutionary models to use as its building blocks. To mitigate the need for such difficult a priori decisions and to avoid the pitfalls of adopting unsuitable models, we plan on employing Bayesian nonparametrics to infer the number of evolutionary models (Beal et al. 2002) and Bayesian model averaging techniques to account for evolutionary model uncertainty (Wu et al. 2013).

Most of the development of MMMs in evolutionary inference has focused on their utility in phylogenetic reconstruction from molecular sequence data. However, there is also a need for MMMs to analyse phenotypic traits. For example, Beaulieu et al. (2013) investigated the evolution of growth habit in the angiosperm clade Campanulidae by employing a covarion-like model that modulates between “fast” and “slow” hidden regimes. We anticipate that the more general MMM formulation we consider here will also prove useful in the context of phenotypic trait evolution, and we plan to extend our implementation to accommodate (discrete) trait data.

Given the potential impact of MMMs on model fit and phylogenetic reconstruction, it is important be able to examine the realisations of such MMMs throughout evolutionary histories, which can be achieved through stochastic mapping. Stochastic mapping of mutations or ancestral trait changes on phylogenies is a widely used method to aid evolutionary hypothesis testing (Nielsen 2002). Beaulieu et al. (2013) developed a stochastic mapping method for MMMs for binary traits that feature two evolutionary rate classes. Irvahn (2016) extended this approach to allow for an arbitrary number of evolutionary rate classes, exploiting efficient sampling techniques that avoid costly matrix exponentiations (Irvahn and Minin 2014). We will build on this work to develop stochastic mapping methods for more general MMMs that can accommodate traits that may assume an arbitrary number of states.

5 Supplementary Material

We have made available an online tutorial on how to construct XML files to perform phylogenetic inference using Markov-modulated models in BEAST: http://beast.community/markov_modulated.html

6 Acknowledgments

We are grateful to Cécile Ané for kindly providing the plant plastid genome data sets. GB acknowledges support from the Interne Fondsen KU Leuven / Internal Funds KU Leuven under grant agreement C14/18/094. The research leading to these results has received funding from the European Research Council under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 725422-ReservoirDOCS). The Artic Network receives funding from the Wellcome Trust through project 206298/Z/17/Z. PL acknowledges support by the Research Foundation – Flanders (‘Fonds voor Wetenschappelijk Onderzoek – Vlaanderen’, G066215N, G0D5117N and G0B9317N). MAS is partially supported through National Science Foundation grant DMS 1264153 and National Institutes of Health grants R01 AI107034 and U19 AI135995. We gratefully acknowledge support from NVIDIA Corporation with the donation of parallel computing resources used for this research.

7 Appendix A: simulation study

Given the large differences in (log) marginal likelihoods obtained in favour of MMMs in the data sets we tested, for example for the psaB and ndhD analyses, we here assess the suitability of Bayesian model testing through (log) marginal likelihood estimation to compare MMMs to standard nucleotide substitution models. To this end, we have simulated three data sets – similar in complexity to the psaB plant plastid data set – using an updated version (1.4) of the $\pi$BUSS sequence simulation package (Bielejec et al. 2014); we provide these XML files in the Supplementary Materials. Each simulation was performed on a randomly drawn tree from the posterior distribution of the psaB plant plastid data set (see section 3.2), as generated under the corresponding substitution model. For each generated data set, we assess the performance of a range of substitution models: an independent general time-reversible model with empirical base frequencies (GTR) and estimated base frequencies (GTR+F), MMMs with two sets of base frequencies and one or two sets of symmetric rate parameters (MMM(GTR)12 and MMM(GTR)22), and MMMs with three sets of base frequencies and either one or three sets of symmetric rate parameters (MMM(GTR)13 and MMM(GTR)33). The model fit to the data for each data set was evaluated by estimating (log) marginal likelihoods using generalized stepping-stone sampling (GSS; Baele et al. (2016)), a recently introduced powerful and flexible approach to perform Bayesian model selection. Computational demands for GSS were increased in a stepwise manner until convergence of the (log) marginal likelihoods was attained, resulting in an initial Markov chain of 20 million iterations followed by 100 power posteriors of 1 million iterations each.

Figure 7 summarizes the outcome of (log) marginal likelihood estimation for a range of possible substitution models under three generative models: GTR, MMM(GTR)22 and MMM(GTR)33. When the nucleotide data was simulated under a GTR substitution model, (log) marginal likelihood estimation using GSS (Baele et al. 2016) selected the GTR model with empirical (or fixed) base frequencies as the model that best fits the data. All other models are penalised for their abundance in parameters, illustrating that models with more parameters are not systematically preferred as the standard GTR model is a nested model within those other models. For the data simulated under an MMM(GTR)22, only the MMM(GTR)33 is able to estimate similar substitution patterns, but yields a slightly lower (log) marginal likelihood, as the substitution process can be adequately captured with fewer parameters by the MMM(GTR)22. When the data is simulated under an MMM(GTR)33, the generative model is retrieved with very strong (log) Bayes factor support as no simpler model is able to capture its complexity. In conclusion, our simulations show that when a Markov-modulated process correctly captures the evolutionary patterns within a data set, (log) Bayes factor support of MMMs over standard nucleotide models will be very strong and differences in (log) marginal likelihood can amount up to hundreds of log units.

8 Appendix B: computational performance

Integral to computing the observed data likelihood are the finite-time transition probabilities that characterize how states transition across a branch in the phylogeny. As discussed earlier, computation of the finite-time transition probability matrix (see Equation 3) through numerical eigendecomposition of $\boldsymbol{\Lambda}$ can be computationally demanding for MMMs owing to their high dimensionality and the lack of closed-form expressions. While the required diagonalization of $\boldsymbol{\Lambda}$ occurs infrequently, the increased dimensionality of MMMs constitutes a computational bottleneck because of the large matrix multiplications associated with the numerical eigendecomposition. Both diagonalization and matrix multiplications can be sped up significantly by exploiting parallelization on state-of-the-art multi-core hardware devices in order to perform matrix multiplications in parallel (Suchard and Rambaut 2009). Recent developments in BEAST (version 1.10; Suchard et al. (2018)) and its accompanying high-performance computational library BEAGLE (version 3; Ayres et al. (2019)) have brought about further improvements in dealing with such demanding use cases, particularly in the presence of multiple partitions.

We benchmark MMMs with varying dimensions against standard nucleotide substitution models, on a range of multi-core CPU and GPU hardware configurations. In terms of raw computing power, highly parallel computing technologies such as GPUs have overtaken traditional CPUs, offering a very cost-effective option to achieve peak performance, with an increasing amount of on-board memory and ever-faster memory bandwidth. We include two multi-core CPU configurations: a 20-core (i.e. $2\times 10$) Intel Xeon (R) Xeon E5-2680 v2 system (Ivy Bridge microarchitecture) and a 24-core (i.e. $2\times 12$) Intel Xeon (R) Xeon E5-2680 v3 system (Haswell microarchitecture), the former with 60 GB/s memory bandwidth and the latter with 68GB/s memory bandwidth. We stack these CPU architectures up against three GPU devices: a Tesla K20X with 2688 CUDA cores and memory bandwidth of 250 GB/s, a Tesla K40 with 2880 CUDA cores and memory bandwidth of 288 GB/s; and a Tesla P100 with 3584 CUDA cores and memory bandwidth of 720 GB/s.

Figure 8 shows a performance comparison of these various architectures on the psaB and ndhD plant plastid data sets analysed in this paper. Across the board, the Tesla P100 GPU outperforms all other hardware configurations, both for standard nucleotide substitution models and MMMs. CPU architectures offer poor performance when evaluating MMMs and are greatly outperformed by GPU cards aimed at the scientific computing market, with graphics cards that are several years old outperforming multi-core CPUs nearly tenfold. Comparing the performance among different GPUs shows that recent developments in the GPU market continue to offer staggering performance improvements, with the Tesla P100 GPU substantially outperforming its predecessors and drastically limiting the loss of performance associated with MMMs. More recently developed GPUs, such as the Tesla V100, promise to deliver even further performance enhancements, but could not yet be included into this study due to lack of availability.