A Generative Angular Model of Protein Structure Evolution

A Generative Angular Model of Protein Structure Evolution Michael Golden,* ,1 Eduardo Garc ıa-Portugue´s, 2,3,5 Michael Sørensen, 3 Kanti V. Mardia, 1,4 Thomas Hamelryck, 5,6 and Jotun Hein 1 1 Department of Statistics, University of Oxford, Oxford, United Kingdom 2 Department of Statistics, Carlos III University of Madrid, Madrid, Spain 3 Department of Mathematical Sciences, University of Copenhagen, Copenhagen, Denmark 4 Department of Mathematics, University of Leeds, Leeds, United Kingdom 5 Bioinformatics Centre, Section for Computational and RNA Biology, Department of Biology, University of Copenhagen, Copenhagen, Denmark 6 Image Section, Department of Computer Science, University of Copenhagen, Copenhagen, Denmark *Corresponding author: E-mail: [email protected]. Associate editor: Jeffrey Thorne Protein families were obtained from the HOMSTRAD database. Abstract Recently described stochastic models of protein evolution have demonstrated that the inclusion of structural information in addition to amino acid sequences leads to a more reliable estimation of evolutionary parameters. We present a generative, evolutionary model of protein structure and sequence that is valid on a local length scale. The model concerns the local dependencies between sequence and structure evolution in a pair of homologous proteins. The evolutionary trajectory between the two structures in the protein pair is treated as a random walk in dihedral angle space, which is modeled using a novel angular diffusion process on the two-dimensional torus. Coupling sequence and structure evolution in our model allows for modeling both “smooth” conformational changes and “catastrophic” conformational jumps, conditioned on the amino acid changes. The model has interpretable parameters and is comparatively more realistic than previous stochastic models, providing new insights into the relationship between sequence and structure evolution. For example, using the trained model we were able to identify an apparent sequence–structure evolutionary motif present in a large number of homologous protein pairs. The generative nature of our model enables us to evaluate its validity and its ability to simulate aspects of protein evolution conditioned on an amino acid sequence, a related amino acid sequence, a related structure or any combination thereof. Key words: evolution, protein structure, probabilistic model, directional statistics. Introduction Recently, several studies (Challis and Schmidler 2012;Herman et al. 2014) have proposed joint stochastic models of evolution which take into account simultaneous alignment of protein sequence and structure. These studies point out the limitations of earlier non-probabilistic methods, which often rely on heuristic procedures to infer parameters of interest. A major disadvantage of using heuristic procedures is that they typically fail to account for sources of uncertainty. For example, relying on a single fixed alignment, which is highly unlikely to be the true underlying alignment, may bias the inference of the posterior distribution over evolutionary trees. We present a generative evolutionary model, ETDBN (Evolutionary Torus Dynamic Bayesian Network) for pairs of homologous proteins. ETDBN captures dependencies between sequence and structure evolution, accounts for alignment uncertainty, and models the local dependencies between aligned sites. A key step in modeling protein structure evolution is selecting a suitable structural representation and corresponding evolutionary model. Early works by Gutin and Badretdinov (1994) and Grishin (1997) represented protein structure using three-dimensional Cartesian coordinates of protein backbone atoms and used diffusions processes to model the relationship between structural distance (measured using RMSD) and sequence similarity. More recent publications by Challis and Schmidler (2012) and Herman et al. (2014) likewise used the three-dimensional Cartesian coordinates of amino acid C a atoms to represent protein structure and additionally used Ornstein-Uhlenbeck (OU) processes to construct Bayesian probabilistic models of protein structure evolution. These models emphasize estimation of evolutionary parameters such as the evolutionary time between species, tree topologies and alignment, and attempt to fully account for sources of uncertainty. For the sake of computational tractability, the aforementioned approaches treat the Cartesian coordinates associated with atoms as evolving independently of another. A non-probabilistic approach by Echave (2008) and Echave and Fern andez (2010) referred to as the Linearly Forced Elastic Network Model (LFENM) treats protein structures as a collection of C a atoms connected by spring forces. The major Article ßThe Author 2017. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons. org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Open Access Mol. Biol. Evol. 34(8):2085–2100 doi:10.1093/molbev/msx137 Advance Access publication April 27, 2017 2085 Downloaded from https://academic.oup.com/mbe/article-abstract/34/8/2085/3772139 by guest on 16 March 2019 benefit of LFENMs is that they do not assume independence of atomic coordinates and take into account non-local dependencies due to physical interactions. In their current formulation LFENMs do not distinguish between the differing chemical nature of different amino acids and therefore do not account for the variable effect of sequence mutation on protein structure evolution. Rather than using a Cartesian coordinate representation, our model, ETDBN, uses a dihedral angle representation motivated by the non-evolutionary TorusDBN model (Boomsma et al. 2008,2014). TorusDBN represents a single protein structure as a sequence of ð/;wÞdihedral angle pairs, which are modeled using continuous bivariate angular distributions (Frellsen et al. 2012). Likewise, ETDBN treats protein structure as a random walk in space, again making use of the /and w dihedral angles (top of fig. 1). The dihedral angle representation is informed by the chemical nature of peptide bonds. Each amino acid in a protein peptide chain is covalently bonded to the next via a peptide bond. Peptide bonds have a partial double bond nature that results in a planar configuration of atoms in space. This configuration allows the protein backbone structure to be largely described in terms of a series of /and wdihedral angles that defines the relationship between the planes in threedimensional space. A benefit of this representation is that it bypasses the need for structural alignment, unlike in models on Cartesian coordinates which typically need to additionally superimpose the structures for comparison purposes (Herman et al. 2014). Accordingly, having to account for superimposition introduces an additional source of uncertainty. A further advantage of the dihedral angle representation is that there are fewer degrees of freedom per amino acid and therefore typically fewer parameters required in order to model their evolution. The evolution of dihedral angles in ETDBN is modeled using a novel stochastic diffusion process developed in Garc ıa-Portugue´s et al. (2017). In addition to this, a coupling is introduced such that an amino acid change can lead to a jump in dihedral angles and a change in diffusion process, allowing the model to capture changes in amino acid that are directionally coupled with changes in dihedral angle or secondary structure. As in Challis and Schmidler (2012) and Herman et al. (2014), the insertion and deletion (indel) evolutionary process is also modeled in order to account for alignment uncertainty (Thorne et al. 1992). The OU processes used in Challis and Schmidler (2012) and Herman et al. (2014) ignore bond lengths and treat C a atoms as evolving independently for the sake of computationally tractability. Furthermore, the OU process makes Gaussian assumptions. From a generative perspective these properties will lead to evolved proteins with C a atoms that are unnaturally dispersed in space. Bond lengths are also ignored in ETDBN, but can be plausibly fixed or modeled. As a result, it is expected that the use of angular diffusions will much more naturally capture the underlying protein structure manifold. Two or more homologous proteins will share a common ancestor, which leads to underlying tree-like dependencies. These dependencies manifest themselves most noticeably in the degree of amino acid sequence similarity between two homologous proteins. The strength of these dependencies is assumed to be a result of two major factors: the time since the common ancestor and the rate of evolution. Failing to account for evolutionary dependencies can lead to false conclusions (Felsenstein 1985), whereas accounting for evolutionary dependencies allows information from homologous proteins to be incorporated in a principled manner. This can lead to more accurate inferences, such as the prediction of a protein structure from a homologous protein sequence and structure, known as homology modeling (Arnold et al. 2006). Stochastic models such as ETDBN are not expected to compete with homology modeling software such as SWISSMODEL (Arnold et al. 2006). However, they allow for estimation of evolutionary parameters and statements about uncertainty to be made in a statistically rigorous manner. Most models of structural evolution ignore dependencies amongst sites because of the increased computational demand and complexity associated with such models. These dependencies are expected to influence patterns of evolution, specifically patterns of amino acid substitution. The current model deals with local dependencies only—dependencies that are expected to arise due to interactions between neighboring amino acids, for example, between amino acids in an a-helix. ETDBN does not account for global dependencies—dependencies that result in the globular nature of proteins (Boomsma et al. 2008). In ETDBN, we attempt to model local dependencies only by using a Hidden Markov Model (HMM) to capture dependencies amongst neighboring aligned positions. HMMs such as PASSML (Li o et al. 1998) have been successfully used to predict protein secondary structure from aligned sequences, however, these models typically havethedisadvantagethattheyassumeacanonicalsecondary structure shared amongst all the sequences being analyzed. This restricts analysis to closely related sequences where conservation of secondary structure is a reasonable assumption. ETDBN does not assume a canonical secondary structure, but instead uses a phylogenetic HMM approach, similar to Siepel and Haussler (2004), that assumes dependencies between evolutionary processes at neighboring aligned positions. Parameters of ETDBN were estimated using 1,200 homologous protein pairs from the HOMSTRAD database (Mizuguchi et al. 1998). The resulting model provides a realistic prior distribution over proteins and protein structure evolution in comparison to previous stochastic models. Doing so enables biological insights into the relationship between sequence and structure evolution, such as patterns of amino acid change that are informative of patterns of structural change (Grishin 2001). It was with these features in mind that ETDBN was developed. Evolutionary Model Overview ETDBN is a dynamic Bayesian network model of local protein sequence and structure evolution along a pair of aligned homologous proteins p a and p b . ETDBN can be can be viewed as an HMM (see fig. 1). Each hidden node of the HMM, Golden et al. .doi:10.1093/molbev/msx137 MBE 2086 Downloaded from https://academic.oup.com/mbe/article-abstract/34/8/2085/3772139 by guest on 16 March 2019 corresponding to an aligned position, adopts an evolutionary hidden state specifying a distribution over three different observations pairs: a pair of amino acid characters, a pair of dihedral angles and a pair of secondary structures classifications. A transition probability matrix specifies neighboring dependencies between adjacent evolutionary states. For example, transitions along the alignment between hidden states encoding predominantly a-helix evolution would be expected to occur more frequently than transitions between an evolutionary hidden state encoding predominantly a-helix evolution and another encoding predominantly b-sheet evolution. Ideally, the underlying hidden states would not just vary across the length of the alignment as captured by the HMM in the current model, but also evolve along the branches of the phylogenetic tree. This remains computationally intractable at present. Allowing the hidden states to evolve along the tree would allow capturing large structural changes, even induced by a single mutation. For now we model such events using a jump model (see below). Partially in order to mitigate this, each hidden state specifies a distribution over a pair of site-classes at each aligned position. This gives rise to the possibility of a ‘jump event’. A jump event allows a large change in dihedral angle or secondary structure (e.g. helix to sheet) to occur at a given aligned position and also introduces a directional coupling between changes in amino acid that are informative of changes in dihedral angle or secondary structure conformation. Observation Types The two proteins, p a and p b , in a homologous pair are associated with a pair of observation sequences O a and O b obtained from experimental data, respectively. An ith site observation pair, Oi¼ðOxðiÞ a;OyðiÞ bÞ, is associated with every aligned site iin an alignment M ab of p a and p b ,whereMi ab 2fy x ; x  ;y  ðÞ gspecifies the homology relationship at position iof the alignment (homologous, deletion with respect to p a and insertion with respect to p a , respectively), iis taken to run from 1 to m,mis the length of the alignment M ab ,andx2f1;...;jpajg and y2f1;...;jpbjg specify the indices of the positions in p a and p b , respectively. jp a j and jp a jgive the number of sites in p a and p b , respectively. Each site observation, OxðiÞ aand OyðiÞ b, contains amino acid and structural information corresponding to the two C a atoms at aligned site ibelonging to each of the two proteins. A site observation corresponding to a particular protein at aligned site i,OxðiÞ a, is comprised of three different data types FIG.1.Above: dihedral angle representation. A small section of a single protein backbone (three amino acids) with /and wdihedral angles shown, together with C a atoms which attach to the amino acid side-chains. Each amino acid side-chain determines the characteristic nature of each amino acid. Every amino acid position corresponds to a hidden node in the HMM below. Note that we only show a single protein, whereas the model considers a pair. Below: depiction of HMM architecture of ETDBN where each Halong the horizontal axis represents an evolutionary hidden node. The horizontal edges between evolutionary hidden nodes encode neighboring dependencies between aligned sites. The arrows between the evolutionary hidden nodes and site-class pair nodes encode the conditional independence between the observation pair variables Axi a;Ayi b(amino acid site pair), Xxi a¼h/xi a;wxi ai;Xyi b¼h/yi b;wyi bi(dihedral angle site pair) and Sxi a;Syi b(secondary structure class site pair). The circles represent continuous variables and the rectangles represent discrete variables. Evolutionary TorusDBN .doi:10.1093/molbev/msx137 MBE 2087 Downloaded from https://academic.oup.com/mbe/article-abstract/34/8/2085/3772139 by guest on 16 March 2019 associated with the C a atom: an amino acid (AxðiÞ a;discrete, one of twenty canonical amino acids), /and wdihedral angles (XxðiÞ a¼h/xðiÞ a;wxðiÞ ai; continuous, bivariate), and a secondary structure classification (SxðiÞ a; discrete, one of three classes: helix (H), sheet (S), or coil (C)). Therefore, OxðiÞ a¼ ðAxðiÞ a;XxðiÞ a;SxðiÞ aÞand OyðiÞ b¼ðAyðiÞ b;XyðiÞ b;SyðiÞ bÞ. Model Structure The sequence of hidden nodes in the HMM is written as H¼ðH1;H2;...;Hm). Each hidden node H i in the HMM corresponds to a site observation pair, OxðiÞ aand OyðiÞ b,atan aligned site iin the alignment M ab . Initially we treat the alignment M ab as given a priori, but later modify the HMM to marginalize out an unobserved alignment. The model is parameterized by hhidden states. Every hidden node H i corresponding to an aligned site itakes an integer value from 1 to qfor the hidden state at node H i .Inturn,each hidden state specifies a distribution over a site-class pair: ðri a;ri bÞas a function of evolutionary time. A site-class pair consists of two site-classes: ri aand ri b.Eachofthetwo site-classes takes an integer value 1 or 2, that is, ðri a;ri bÞ2fð1;1Þ;ð1;2Þ;ð2;1Þ;ð2;2Þg. We return to the specific role of the site-classes pairs in the next section. The state of H i together with the site-class pair, ðri a;ri bÞ, and the evolutionary time separating proteins p a and p b ,t ab , specify a distribution over three conditionally independent stochastic processes describing each of the three types of site observation pairs: Ai¼ðAxðiÞ a;AyðiÞ bÞ;Xi¼ðXxðiÞ a;XyðiÞ bÞand Si¼ðSxðiÞ a;SyðiÞ bÞ. This conditional independence structure allows the likelihood of a site observation pair at an aligned site i to be written as follows: pðOijHi;ri a;ri b;tabÞ¼pðAijHi;ri a;ri b;tabÞ z}|{ amino acid evolution pðXijHi;ri a;ri b;tabÞ z}|{ dihedral angle evolution pðSijHi;ri a;ri b;tabÞ z}|{ secondary structure evolution :(1) The assumption of conditional independence provides computational tractability, allowing us to avoid costly marginalization when certain combinations of data are missing (e.g. amino acid sequences present, but secondary structures and dihedral angles missing). Stochastic Processes: Modeling Evolutionary Dependencies Each site-class couples together three time-reversible stochastic processes that separately describe the evolution of the three pairs of observation types, as in equation (1).Each site-class is intended to capture both physical and evolutionary features pertaining to sequence and structure. Parameters that correspond to a particular site-class are termed site-class specific, whereas parameters that are shared across all siteclasses are termed global. The use of site-class specific parameters, such as site-class specific amino acid frequencies and dihedral angle diffusion parameters, as described in the next section, is intended to model site-specific physical–chemical properties (Halpern and Bruno 1998;Koshi and Goldstein 1998;Lartillot and Philippe 2004). Amino Acid Evolution As is typical with models of sequence evolution, amino acid evolution, pðAxðiÞ a;AyðiÞ bjHi;ri a;ri b;tabÞ,isdescribedbya Continuous-Time Markov Chain (CTMC). Each amino acid CTMC is parameterized in the following way: the exchangeability of amino acids is described by a 20 20 symmetric global exchangeability matrix S(190 free parameters; Whelan and Goldman 2001), a site-class specific set of 20 amino acid equilibrium frequencies Ph r¼diagfp1;p2;...;p20g(19 free parameters per site-class) and a site-class specific scaling factor Kh r(one free parameter per site-class). Together these parameters define a site-class specific time-reversible amino acid rate matrix Qh r¼Kh rSPh r. The stationary distribution of Qh ris given by the amino acid equilibrium frequencies: Ph r. Secondary Structure Evolution Secondary structure evolution, pðSxðiÞ a;SyðiÞ bjHi;ri a;ri b;tabÞ,is also described by a CTMC. For the sake of simplicity we use only three discrete classes to describe secondary structure at each position: helix (H), sheet (S), and random coil (C). The exchangeability of secondary structure classes at a position is described by a 3 3 symmetric global exchangeability matrix Vand a site-class specific set of three secondary structure equilibrium frequencies Xh r¼diagfp1;p2;p3g.Together they define a site-class specific time-reversible secondary structure rate matrix Rh r¼VXh r, with stationary distribution: Xh r. Dihedral Angle Evolution Central to our model is evolutionary dependence between dihedral angles, pðXxðiÞ a;XyðiÞ bjHi;ri a;ri b;tabÞ. Typically, the continuous-time evolution of the continuous-state random variablesismodeledbyadiffusiveprocesssuchastheOU process, as in Challis and Schmidler (2012).However,anOU process is not appropriate for dihedral angles as they have a natural periodicity. For this reason, a bivariate diffusion that captures the periodic nature of dihedral angles, the Wrapped Normal (WN) diffusion, was specifically developed for this paper in Garc  ıa-Portugue´s et al. (2017). Topologically, the WN diffusion (see fig. 2 for a pictorial example)canbethoughtofastheanalogueoftheOUprocess on the torus T2¼½p;pÞ½p;pÞ. The WN diffusion arises as the wrapping on T2of the following Euclidean diffusion: dXt¼AX k2Z2ðlXt2kpÞwkðXtÞ z}|{ coefficient drift dt þR1 2 z}|{ coefficient diffusion dWt; (2) where W t is the two-dimensional Wiener process, Ais the drift matrix, l2T2is the stationary mean, Ris the infinitesimal covariance matrix and Golden et al. .doi:10.1093/molbev/msx137 MBE 2088 Downloaded from https://academic.oup.com/mbe/article-abstract/34/8/2085/3772139 by guest on 16 March 2019 wkðhÞ¼ /1 2A1Rðhlþ2kpÞ P m2Z2 /1 2A1Rðhlþ2mpÞ;h2T2;(3) is a probability density function (pdf) for k2Z2./Rstands for the pdf of a bivariate Gaussian Nð0;RÞ.Thepdf(eq. 3) weights the linear drifts of equation (2) such that they become smooth and periodic. It is shown in Garc  ıa-Portugue´s et al. (2017) that the stationary distribution of the WN diffusion is a WNðl;RÞ,which has pdf: pWNðhjl;RÞ¼X k2Z2 /Rðhlþ2kpÞ:(4) Despite involving an infinite sum over Z2, taking just the first few terms of this sum provides a tractable and accurate approximation to the stationary density for most of the realistic parameter values. Maximum Likelihood Estimation (MLE) for diffusions is based on the transition probability density (tpd), which only has a tractable analytical form for very few specific processes. A highly tractable and accurate approximation to the tpd is given for the WN diffusion. This approximation results from weighting the tpd of the OU process in the same fashion as the linear drifts are weighted in equation (2), yielding the following multimodal pseudo-tpd: ~ pðh2jh1;A;l;R;tÞ¼X m2Z2 pWNðh2jlm t;CtÞwmðh1Þ;(5) with h1;h22T2;lm t¼lþetAðh1lþ2pmÞand Ct¼Ðt 0esAResATds. The pseudo-tpd provides a good approximation to the true tpd in key circumstances: 1) t!0, since it collapses in the Dirac delta; 2) t!1, since it converges to the stationary distribution; 3) high concentration, since the WN diffusion becomes an OU process. Furthermore, it is shown in Garc  ıa-Portugue´s et al. (2017) that the pseudotpd has a lower Kullback–Leibler divergence with respect to the true tpd than the Euler and Shoji-Ozaki pseudo-tpds, for most typical scenarios and discretization times in the diffusion trajectory. A further desirable property of the pseudo-tpd is that it obeys the time-reversibility equation, which in terms of ðXxðiÞ a; XyðiÞ bÞis ~ pðXyðiÞ bjXxðiÞ a;A;l;R;tabÞpWNðXxðiÞ ajl;1 2A1RÞ ¼~ pðXxðiÞ ajXyðiÞ b;A;l;R;tabÞpWNðXyðiÞ bjl;1 2A1RÞ: Indeed, the WN diffusion is the unique time-reversible diffusion with constant diffusion coefficient and stationary pdf (eq. 4), in the same way the OU is with respect to a Gaussian. Time-reversibility is an assumption of the overall model and many other models of sequence evolution. A benefit of timereversibility in a pairwise model such as ETDBN is that one of the proteins in a pair may be arbitrarily chosen as the ancestor, thus avoiding a computationally expensive marginalization of an unobserved ancestor. The likelihood of a dihedral angle observation pair ðXxðiÞ a;XyðiÞ bÞ, assuming that XxðiÞ ais drawn from the stationary distribution, is given by: pðXxðiÞ a;XyðiÞ bjHi;ri a;ri b;tabÞ ¼pðXxðiÞ a;XyðiÞ bjA;l;R;tabÞ ~ pðXyðiÞ bjXxðiÞ a;A;l;R;tabÞpWNðXxðiÞ ajl;1 2A1RÞ; (6) Aand Rare constrained to yield a covariance matrix A1R. A parameterization that achieves this is R¼diagðr2 1 ;r2 2Þand A¼ða1;r1 r2a3;r2 r1a3;a2Þ;a1a2>a2 3.a 1 and a 2 are the drift components for the /and wdihedral angles, respectively. Dependence (correlation) between the dihedral angles is captured by a 3 . A depiction of a WN diffusion with given drift and diffusion parameters is shown in figure 2. Site-Classes: Constant Evolution and Jump Events We now turn to the meaning of the site-class pairs. Two modes of evolution are modeled: constant evolution and jump events. Constant evolution occurs when the site-class starting in protein p a at aligned site i,ri a,isthesameasthe site-class ending in protein p b at aligned site i,ri b,thatis, ri a¼ri b. Thus the distribution over observation pairs at a site is specified by a single site-class. As already stated, a site-class specifies the parameters of the three conditionally independent stochastic processes describing amino acid, dihedral angle, and secondary structure evolution. A limitation of “constant evolution” is that the coupling between the three stochastic processes is somewhat weak. This in part stems from the time-reversibility of the stochastic processes—swapping the order of one of the three observation pairs at a homologous site, e.g. (Glycine, Proline) instead of (Proline, Glycine), does not alter the likelihood in equation (1). Alternatively restated from a generative perspective: a ‘directional coupling’ of an amino acid interchange does not inform the direction of change in dihedral angle or secondary structure. For example, replacing a glycine in an ahelix in one protein with a proline at the homologous position in a second protein would be expected to break the ahelix in the second protein and to strongly inform the plausible dihedral angle conformations in the second protein. Ideally, we would consider a model in which the underlying site-classes were not fixed over the evolutionary trajectory separating the two proteins, as in the case of constant evolution as described above, but instead were able to ‘evolve’ in time. This would allow occasional switches in the underlying site-class at a particular homologous site, which would create a stronger dependency between amino acid, dihedral angle and secondary structure evolution, that furthermore captures the directional coupling we desire. Such an approach is considered computationally intractable due to the introduction of context-dependence when having to consider neighboring dependencies amongst evolutionary trajectories at adjacent sites. In order to approximate this ‘ideal’ model in a computationally efficient manner we introduce the notion of a jump Evolutionary TorusDBN .doi:10.1093/molbev/msx137 MBE 2089 Downloaded from https://academic.oup.com/mbe/article-abstract/34/8/2085/3772139 by guest on 16 March 2019 event. A jump event occurs when ri a6¼ ri b.Whereasconstant evolution is intended to capture angular drift (changes in dihedral angles localized to a region of the Ramachandran plot), a jump event is intended to create a directional coupling between amino acid and structure evolution, and is also expected to capture angular shift (large changes in dihedral angles, possibly between distant regions of the Ramachandran plot). The hidden state at node H i , together with the evolutionary time t ab separating proteins p a and p b , specifies a joint distribution over a site-class pair: pðri a;ri bjHi;tabÞ¼pðri ajHi;ri b;tabÞpðri bjHiÞ;(7) where pðri ajHi;ri b;tabÞ ¼ ecHitab þpHi;ri bð1ecHitab Þ;if ri a¼ri b; pHi;ri bð1ecHitab Þ;if ri a6¼ ri b; 8 < : and pðri ajHiÞ¼pHi;ri aand pðri bjHiÞ¼pHi;ri b.pHi;ri aand pHi;ri b are model parameters specifying the probability of starting in site-class ri aor ri b, respectively, corresponding to the hidden state specified by node H i .cHi>0 is a model parameter giving the jump rate corresponding to the hidden state given by node H i . The site-class jump probabilities have been chosen so that time-reversibility holds, in other words: pðri ajHi;ri b;tabÞpðri bjHiÞ¼pðri bjHi;ri a;tabÞpðri ajHiÞ: The hidden state at node H i , together with a site-class pair ðri a;ri bÞand the evolutionary time t ab , specifies the joint likelihood over site observation pairs: pðOxðiÞ a;OyðiÞ bjHi;ri a;ri b;tabÞ ¼ pðOxðiÞ a;OyðiÞ bjHi;ri c;tabÞ;if ri a¼ri b¼ri c; pðOxðiÞ ajHi;ri aÞpðOyðiÞ bjHi;ri bÞ;if ri a6¼ ri b: 8 > < > : (8) In the case of constant evolution, evolution at aligned iis described in terms of the same site-class ri c. Evolution is considered constant because each observation type is drawn from a single stochastic process specified by H i and r c .Note that the strength of the evolutionary dependency within an observation pair is a function of the evolutionary time t ab . In the case of a jump event, the evolutionary processes are, after the evolutionary jump, restarted independently in the stationary distribution of the new site-class. Thus the site observations OxðiÞ aand OyðiÞ bareassumedtobedrawnfrom the stationary distributions of two separate stochastic processes corresponding to site-classes ri aand ri b, respectively. This implies that, conditional on a jump, the likelihood of the observations is no longer dependent on t ab .Ajumpevent is an abstraction that captures the end-points of the evolutionary process, but ignores the potential evolutionary trajectory linking the two site observations. The advantage of abstracting the evolutionary trajectory is that there is no need to perform a computationally expensive marginalization over all possible trajectories, as might be necessary in a model where the hidden states evolve along a tree. The likelihood of an observation pair is now simply a sum over the four possible site-class pairs: pðOxðiÞ;OyðiÞjHi;tabÞ ¼P ðri a;ri bÞ2R pðOxðiÞ;OyðiÞjHi;ri a;ri b;tabÞpðri a;ri bjHi;tabÞ; where R¼fð1;1Þ;ð1;2Þ;ð2;1Þ;ð2;2Þg is the set of four site-class pairs, pðOxðiÞ;OyðiÞjHi;ri a;ri bÞis given by equation (8) and pðri a;ri bjHi;tabÞis given by equation (7). Identification of Evolutionary Motifs Encoding Jump Events In order to identify aligned sites having potential evolutionary motifs encoding jump events, a specific criterion was developed. For a particular protein pair, inference was performed under the model conditioned on the amino acid sequence and dihedral angles for both proteins ðAa;Ab;Xa;XbÞ. Homologous sites corresponding to a single hidden state andwithevidenceofajumpevent(ri a6¼ ri b)atposterior probability >0.90 were identified, that is, the i’s such that pðHi;ri a6¼ ri bjAa;Ab;Xa;XbÞ>0:90. In a second filtering step, amino acid sequences and a single set of dihedral angles corresponding to one of the FIG.2.Drift vector field for the WN diffusion with A¼ð1;0:5;0:5; 0:5Þ;l¼ð0;0Þand R¼ð1:5Þ2I. The color gradient represents the Euclidean norm of the drift. The contour lines represent the stationary distribution. An example trajectory starting at x 0 ¼(0,0) and ending at x 2 , running in the time interval [0, 2] is depicted using a white to red color gradient indicating the progression of time. The periodic nature of the diffusion can be seen by the wrapping of both the stationary diffusion and the trajectory at the boundaries of the square plane. The fact that stationary distribution is not aligned with the horizontal and vertical axes illustrates the dependence (given by a 3 ) between the uand wdihedral angles. Golden et al. .doi:10.1093/molbev/msx137 MBE 2090 Downloaded from https://academic.oup.com/mbe/article-abstract/34/8/2085/3772139 by guest on 16 March 2019 proteins were used (Aa;Ab;Xaor Aa;Ab;Xb)toinferthe posterior probability, this time at a lower threshold: pðHi;ri a6¼ ri bjAa;Ab;XaÞ>0:50 or pðHi;ri a6¼ ri bjAa;Ab;XbÞ>0:50. This second criterion ensured that the evolutionary motif was identifiable under typical conditions where one has limited access to structural information (in this case a single protein structure in a pair). Only those aligned sites meeting both criteria were selected for downstream analysis. Statistical Alignment: Modeling Insertions and Deletions Protein sequences can not only undergo amino transitions due to underlying nucleotide mutations in the coding sequence, but also indel events. To account this, a modified pairwise TKF92 alignment HMM based on Mikl os et al. (2004) was implemented. The TKF92 alignment HMM was augmented with the ETDBN evolutionary hidden states in order to capture local sequence and structure evolutionary dependencies. Furthermore, it was modified such that neighboring dependencies amongst hidden states at adjacent alignment sites were modeled. For more details see ‘Statistical alignment’ in Supplementary Material online. Whilst it is possible to fix the alignment in advance by prealigning the sequences using one of the many available alignment methods (Katoh et al. 2002;Edgar 2004)orusinga curated alignment (such as from the HOMSTRAD database), doing so ignores alignment uncertainty. Training and Test Datasets A training dataset of 1,200 protein pairs (2,400 proteins; 417,870 site observation pairs) and a test dataset of 38 protein pairs (76 proteins; 14,125 site observation pairs) were assembled from 1,032 protein families in the HOMSTRAD database. For further details see ‘Construction of test and training datasets’ in the supplementary Material, Supplementary Material online. Model Training and Selection Maximum Likelihood Estimation of the model parameters, ^ W, was done using Stochastic Expectation Maximization (StEM, Gilks et al. 1995). For further details of the Eand M-steps of the StEM algorithm and for details about model selection please refer to ‘Model training and selection’ in the supplementary Material, Supplementary Material online. Results and Discussion Selecting the Number of Hidden States Fourteen models were trained (8, 16, 32, 48, 52, 56, 60, 64, 68, 72, 76, 80, 96, and 112 hidden state models). The 64 hidden state model was chosen as the best model, as it had the lowest Bayesian Information Criterion (BIC, fig. 3). Stationary Distributions over Dihedral Angles Capture the Empirical Distribution Figure 4 illustrates the sampled and empirical dihedral angle distributions. There is a good correspondence between dihedral angles sampled under the model (fig. 4,left)andthe empirical distribution of dihedral angles in our training dataset (fig. 4, right) for all three cases illustrated (all amino acids, glycine only, and proline only). The correspondence is not surprising given that ETDBN is effectively a mixture model with a large number of mixture components. Estimates of Evolutionary Time from Dihedral Angles Are Consistent with Estimates from Sequence Whilst ETDBN has a general scope with respect to applications (including acting as proposal distribution or as a building block in a homology modeling application), we envision the primary application being inference of evolutionary parameters. Figure 5 compares evolutionary times estimated using only pairs of homologous amino acid sequences versus only pairs of homologous dihedral angles. As desired, the two estimates of evolutionary time for each protein pair are similar, as can be seen by the proximity of the points to the identity line. Apairedt-test gave a p-value of 0.578, thus failing to reject the null hypothesis that there is no difference between branch lengths estimated using sequence only versus angles only. This indicates that there is sufficient evolutionary information in the dihedral angles to estimate the evolutionary times and that the model is consistent in its estimates, lacking a significant tendency to underestimate or overestimate the evolutionary times when either sequence or dihedral angles are used. Interestingly, the variance in the sampled evolutionary timesishigherwhendihedralanglesonlyareused,ascompared with sequence only (see fig. S1, Supplementary Material online). The Relationship between Evolutionary Time and Angular Distance Is Adequately Modeled We investigated the relationship between evolutionary time and angular distance between real protein pairs and protein pairs where the dihedral angles of p b (X b ) were treated as missing and hence sampled (fig. 6). As expected, for both real and sampled pairs, angular distance tends to increase as a function evolutionary time. For larger evolutionary times a plateau begins to emerge, which is expected as the maximum possible theoretical angular distance is ffiffiffi8 p2:828. When the evolutionary time is exactly zero (t ab ¼0) under our model, the angular distance between sampled dihedral angles is exactly zero (not shown in fig. 6). However, this is not expected to be the case for real protein pairs when the two sequences are identical (due to the inherently flexible nature of proteins, different experimental conditions, experimental noise, etc.). It is therefore not surprising that the regression curve for the real protein pairs does not pass through zero. Evolutionary TorusDBN .doi:10.1093/molbev/msx137 MBE 2091 Downloaded from https://academic.oup.com/mbe/article-abstract/34/8/2085/3772139 by guest on 16 March 2019 For small evolutionary times (<0.2) the curves for the real and sampled protein pairs show a good correspondence, however, for larger evolutionary times the model tends to underestimate angular distances. This may reflect the fact that the tpd of the WN diffusion specified is localized around its mean, even when the evolutionary time is large; therefore, dihedral angles distant from this mean are unlikely to be sampled. To a certain extent this is mitigated by the jump model, which occasionally allows for large changes in dihedral angle, but may still be somewhat limited in its flexibility, as jumps can only occur between two site-classes. The majority of protein pairs in our training dataset represent smaller evolutionary times (81.7% of evolutionary times are smaller than 0.4) and therefore protein pairs with larger evolutionary times and their associated jumps are underrepresented in our dataset, which may also explain the underestimation. An additional possibility is that ETDBN does not attempt to model global dependencies. Echave and Fern andez (2010) use a LFENM model (which does take into account global dependencies) and provide evidence showing that the majority of structural changes is due to collective global deformations rather than local deformations. A local model such as ETDBN, by definition, does not take into account global dependencies and therefore does not fully account for their contribution to structural divergence. Evaluation of the Model The conditional independence structure in (eq. 1)enables computationally efficient sampling from the model under different combinations of observed or missing data. For example, ETDBN can be used to sample (i.e., predict) the dihedral angles of a protein from its corresponding amino acid sequence, a homologous amino acid sequence, a homologous set of dihedral angles, the corresponding secondary structure, a homologous secondary structure, or any combination of them. Predictive accuracy was measured using 38 homologous protein pairs in the test dataset. For every protein pair (p a ,p b ), the dihedral angles of p b in each pair were treated as missing, and these missing dihedral angles were sampled under the model given a particular combination of observation types. Theaverageangulardistancebetweenthesampledand known dihedral angles was used as the measure of predictive accuracy. Figure 7 gives an example of predictive accuracy under different combinations of observations types overlaid on a cartoon structure of the protein structure being predicted, whereas figure 8 provides a representative view of predictive accuracy across 10 different protein pairs in the test dataset for different combinations of observations types. We highlight some of the key patterns identified in figures 7 and 8as follows. Combination 1 refers to random sampling from the model, implying no data observations were conditioned on besides therespectivelengthsofproteinsp a and p b .Theaverage angular distance between the true and predicted dihedral angleswas1.6.Randomsamplingactsasabaselineforpredictive accuracy. It is apparent from figure 7 that the model has a propensity to predict right-handed a-helices, which is the most populated region in the Ramachandran plot. Under combination 2, only the amino acid sequence corresponding to p b is observed. As expected in figures 7 and 8 there is an increase in predictive accuracy with the addition of the amino acid sequence relative to combination 1. Under combination 3, we add in the amino acid sequence of a homologous protein (p a ). In all ten cases there is an improvement in predictive accuracy. The improvement in predictive accuracy is reasonable, as knowledge of the sequence evolutionary trajectory is expected to encode information about structure evolution and hence will inform the dihedral angle conformational possibilities. Under combination 4, in addition to the two amino acid sequences we treat the homologous secondary structure as observed. This results in a substantial improvement in predictive accuracy as one would expect. Knowledge of the amino acid sequence and a homologous secondary structure strongly informs regions of the Ramachandran plot that are likely to be occupied. Under combination 5 (which we consider the canonical combination—the standard homology modeling scenario), we treat both amino acid sequences as observed, as well as the dihedral angles of the homologous protein (p a )—in all cases the predictive accuracy improves over combination 4. This is anticipated as the homologous dihedral angles are expected to be the best proxy for missing dihedral angles and are therefore expected to be more informative than secondary structure alone. Note that the availability of a homologous amino acid sequence pair here and in combination 4 is consequential as it informs the evolutionary time t ab parameter, which will typically constrain the distribution over dihedral angles and reduce the associated uncertainty. Finally, in combination 6, the same data observations as in combination 5 are used, except the alignment is treated as given a priori (by the HOMSTRAD alignment) FIG.3.BIC scores (points) and number of free parameters (curve) as a function of the number of hidden states across 14 models (indicated above the dotted vertical lines) trained using the 1,200 protein pairs in the training dataset. A 64 hidden state model had the lowest BIC score. Each model represents the best of several attempts. Golden et al. .doi:10.1093/molbev/msx137 MBE 2092 Downloaded from https://academic.oup.com/mbe/article-abstract/34/8/2085/3772139 by guest on 16 March 2019 rather than as unobserved. The HOMSTRAD alignment is based on a structural and sequence alignment of p a and p b and therefore is expected to encode a higher degree of homology and structural information than combination 5 (where the alignment is treated as unobserved and therefore a marginalization over alignments is performed). On average, there is a slight improvement in predictive accuracy when fixing the alignment, albeit the magnitude of improvement is not substantial. This demonstrates the accuracy of the alignment HMM. FIG.4.Ramachandran plots depicting sampled and empirical dihedral angle distributions. The top row depicts the distributions for all amino acids, the middle for glycine only and the bottom for proline only. The leftmost plots show dihedral angles sampled under the jump model, whereas the rightmost plots show the empirical distributions of dihedral angles in the training dataset. Evolutionary TorusDBN .doi:10.1093/molbev/msx137 MBE 2093 Downloaded from https://academic.oup.com/mbe/article-abstract/34/8/2085/3772139 by guest on 16 March 2019 Supplementary Material Supplementary material are available at Molecular Biology and Evolution online. Acknowledgments The authors acknowledge funding from the University of Copenhagen 2016 Excellence Programme for Interdisciplinary Research (UCPH2016-DSIN). The second author acknowledges support from project MTM2016-76969-P from the Spanish Ministry of Economy and Competitiveness and ERDF. Authors acknowledge valuable comments from three referees that led to substantial improvements in the manuscript, as well as initial discussions with Christian Havn, Ian Lim and Mathias Cronj€ ager. References Arnold K, Bordoli L, Kopp J, Schwede T. 2006. The SWISS-MODEL workspace: a web-based environment for protein structure homology modelling. Bioinformatics 22(2):195–201. Boomsma W, Mardia KV, Taylor CC, Ferkinghoff-Borg J, Krogh A, Hamelryck T. 2008. A generative, probabilistic model of local protein structure. Proc Natl Acad Sci U S A. 105(26):8932–8937. Boomsma W, Tian P, Frellsen J, Ferkinghoff-Borg J, Hamelryck T, LindorffLarsen K, Vendruscolo M. 2014. Equilibrium simulations of proteins using molecular fragment replacement and NMR chemical shifts. Proc Natl Acad Sci U S A. 111(38):13852–13857. Challis CJ, Schmidler SC. 2012. A stochastic evolutionary model for protein structure alignment and phylogeny. MolBiolEvol. 29(11): 3575–3587. Echave J. 2008. Evolutionary divergence of protein structure: the linearly forced elastic network model. Chem Phys Lett. 457(4):413–416. Echave J, Fern andez FM. 2010. A perturbative view of protein structural variation. Proteins: Struct Funct Bioinform. 78(1):173–180. Edgar RC. 2004. MUSCLE: multiple sequence alignment with high accuracyandhighthroughput.Nucleic Acids Res. 32(5):1792–1797. Felsenstein J. 1985. Phylogenies and the comparative method. Am Nat. 125(1):1–15. Frellsen J, Mardia KV, Borg M, Ferkinghoff-Borg J, Hamelryck T. 2012. Towards a general probabilistic model of protein structure: the reference ratio method. In: Bayesian methods in structural bioinformatics. New York: Springer, p. 125–134. Garcıa-Portugue´s E, Sørensen M, Mardia KV, Hamelryck T, 2017. Langevin diffusions on the torus: estimation and applications. arXiv:1705.00296. Gilks WR, Richardson S, Spiegelhalter D. 1995. Markov chain Monte Carlo in practice. London: CRC Press. Grishin NV. 1997. Estimation of evolutionary distances from protein spatial structures. J Mol Evol. 45(4):359–369. Grishin NV. 2001. Fold change in evolution of protein structures. JStruct Biol. 134(2):167–185. Gutin AM, Badretdinov AY. 1994. Evolution of protein 3D structures as diffusion in multidimensional conformational space. JMolEvol. 39(2):206–209. Halpern AL, Bruno WJ. 1998. Evolutionary distances for protein-coding sequences: modeling site-specific residue frequencies. MolBiolEvol. 15(7):910–917. Hamelryck T, Borg M, Paluszewski M, Paulsen J, Frellsen J, Andreetta C, Boomsma W, Bottaro S, Ferkinghoff-Borg J. 2010. Potentials of mean force for protein structure prediction vindicated, formalized and generalized. PLoS ONE. 5(11):e13714. HermanJL,ChallisCJ,Nov  ak  A, Hein J, Schmidler SC. 2014. Simultaneous Bayesian estimation of alignment and phylogeny under a joint model of protein sequence and structure. MolBiolEvol. 31(9):2251–2266. Holm L, Rosenstro¨m P. 2010. Dali server: conservation mapping in 3D. Nucleic Acids Res. 38(Web Server issue):W545–W549. KatohK.etal.2002.MAFFT:anovelmethod for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 30(14):3059–3066. Koshi JM, Goldstein RA. 1998. Models of natural mutations including site heterogeneity. Proteins 32(3):289–295. Lartillot N, Philippe H. 2004. A Bayesian mixture model for across-site heterogeneities in the amino-acid replacement process. Mol Biol Evol. 21(6):1095–1109. Li o P, Goldman N, Thorne JL, Jones DT. 1998. PASSML: combining evolutionary inference and protein secondary structure prediction. Bioinformatics 14(8):726–733. Mikl os I, Lunter G, Holmes I. 2004. A long indel model for evolutionary sequence alignment. MolBiolEvol. 21(3):529–540. Mizuguchi K, Deane CM, Blundell TL, Overington JP. 1998. HOMSTRAD: a database of protein structure alignments for homologous families. Protein Sci. 7(11):2469–2471. Ortiz AR, Strauss CE, Olmea O. 2002. Mammoth (matching molecular models obtained from theory): an automated method for model comparison. Protein Sci. 11(11):2606–2621. Robinson DM, Jones DT, Kishino H, Goldman N, Thorne JL. 2003. Protein evolution with dependence among codons due to tertiary structure. Mol Biol Evol. 20(10):1692–1704. Rohl CA, Strauss CE, Misura KM, Baker D. 2004. Protein structure prediction using Rosetta. Methods Enzymol. 383:66–93. Schwartz AS, Myers EW, Pachter L. 2005. Alignment metric accuracy. arXiv preprint q-bio/0510052. Shindyalov IN, Bourne PE. 1998. Protein structure alignment by incremental combinatorial extension (CE) of the optimal path. Protein Eng. 11(9):739–747. Siepel A, Haussler D. 2004. Combining phylogenetic and hidden Markov models in biosequence analysis. JComputBiol. 11(2–3):413–428. Thorne JL, Kishino H, Felsenstein J. 1992. Inching toward reality: an improved likelihood model of sequence evolution. JMolEvol. 34(1):3–16. Whelan S, Goldman N. 2001. A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach. MolBiolEvol. 18(5):691–699. Yu J, Thorne JL. 2006. Dependence among sites in RNA evolution. Mol Biol Evol. 23(8):1525–1537. ZhangY,SkolnickJ.2005.TM-align:aproteinstructurealignment algorithm based on the TM-score. Nucleic Acids Res. 33(7):2302–2309. Golden et al. .doi:10.1093/molbev/msx137 MBE 2100 Downloaded from https://academic.oup.com/mbe/article-abstract/34/8/2085/3772139 by guest on 16 March 2019