11footnotetext:  Department of Statistics, University of Oxford
(UK).22footnotetext:  Bioinformatics Centre, Section for Computational and RNA
Biology, Department of Biology, University of Copenhagen
(Denmark).33footnotetext:  Department of Mathematical Sciences, University of
Copenhagen (Denmark).44footnotetext:  Department of Mathematics, University of
Leeds (UK).55footnotetext:  Image Section, Department of Computer Science,
University of Copenhagen (Denmark).66footnotetext: Corresponding author.
e-mail<EMAIL_ADDRESS>
# A generative angular model of protein structure evolution

Michael Golden1,6 Eduardo García-Portugués2,3 Michael Sørensen3 Kanti V.
Mardia4 Thomas Hamelryck2,5 Jotun Hein1

###### 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
modelled using a novel angular diffusion process on the two-dimensional torus.
Coupling sequence and structure evolution in our model allows for modelling
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.

Keywords: Evolution; protein structure; probabilistic model; directional
statistics.

##  1 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, they may rely on a single fixed alignment of
proteins, which is highly unlikely to be the true underlying alignment. Such a
fixed alignment may bias the inference of the 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. This
extends the models presented in Challis and Schmidler (2012) and Herman et al.
(2014), which emphasise estimation of evolutionary parameters such as the
evolutionary time between species, tree topologies and alignment.

ETDBN is motivated by the non-evolutionary TorusDBN model (Boomsma et al.,
2008). TorusDBN is a first-order Hidden Markov Model (HMM) that represents a
single protein structure as a sequence of $\phi,\psi$ dihedral angle pairs,
which are modelled using continuous bivariate angular distributions (Frellsen
et al., 2012). Likewise, ETDBN treats protein structure evolution as a random
walk in space, again making use of the $\phi$ and $\psi$ dihedral angles (top
of Figure 1).

Figure 1: Above: dihedral angle representation. A small section of a single
protein backbone (three amino acids) with $\phi$ and $\psi$ dihedral angles
shown, together with $C\alpha$ 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 $H$
along the horizontal axis represents an evolutionary hidden node. The
horizontal edges between evolutionary hidden nodes encode neighbouring
dependencies between aligned sites. The arrows between the evolutionary hidden
nodes and evolutionary regime pair nodes encode the conditional independence
between the observation pair variables $A_{a}^{x_{i}},A_{b}^{y_{i}}$ (amino
acid site pair),
$X_{a}^{x_{i}},X_{b}^{y_{i}}=\langle\phi_{a}^{x_{i}},\psi_{a}^{x_{i}}\rangle,\langle\phi_{b}^{y_{i}},\psi_{b}^{y_{i}}\rangle$
(dihedral angle site pair) and $S_{a}^{x_{i}},S_{b}^{y_{i}}$ (secondary
structure class site pair) The circles represent continuous variables and the
rectangles represent discrete variables.

Each amino acid in a peptide chain is covalently bonded to the next via a
chemical bond referred to as 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 $\phi$ and $\psi$ dihedral angles that
defines the relationship between the planes in three-dimensional space. A
benefit of this representation is that it bypasses the need for structural
alignment, unlike in Challis and Schmidler (2012) and Herman et al. (2014),
which is required to simultaneously superimpose structures in three-
dimensional space for comparison purposes. Furthermore, dihedral angles
provide a natural character state for aligning positions in the absence of
amino acid characters.

The evolution of dihedral angles in ETDBN is modelled using a novel stochastic
diffusion process developed in García-Portugué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 us to
capture changes in amino acid that are 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
modelled in order to account for alignment uncertainty (Thorne et al., 1992).

Given that naturally occurring protein sequences are evolutionarily related,
it is only natural to consider models that take such evolution into account.
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 main factors: the time since the common ancestor and the rate of
evolution. If the rate of evolution is assumed constant, recently separated
proteins will typically be similar in both sequence and structure, whereas
proteins which share a common ancestor in the distant past will plausibly
appear independent of one another.

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 modelling. An example of a popular homology modelling software
environment is SWISS-MODEL (Arnold et al., 2006). Stochastic models such as
ETDBN are not expected to compete with packages such as SWISS-MODEL. However,
they allow for estimation of evolutionary parameters and statements about
uncertainty to be made in a statistically rigorous manner.

Parameters of ETDBN, particularly those governing the amino acid and dihedral
angle evolutionary processes, were learnt during a training phase using 1200
homologous protein pairs from the HOMSTRAD database (Mizuguchi et al., 1998).
This resulted in a realistic prior distribution over proteins compared to
previous stochastic models, enabling 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.

##  2 Evolutionary model

###  2.1 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 Figure 1). Each hidden node
of the HMM, 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 neighbouring dependencies between adjacent evolutionary states. For
example, it is expected that evolutionary hidden states specifying
predominantly $\alpha$-helical secondary structure configurations will occur
adjacent to one another. Finally, each hidden state specifies a distribution
over a pair of evolutionary regimes at each aligned position. This introduces
a stronger coupling between changes in amino acid that are informative of
changes in dihedral angle or secondary structure between proteins, whilst
simultaneously enabling the possibility of a jump – a large change in dihedral
angle or secondary structure, (e.g. helix to sheet) coupled with specific
amino acid changes at an aligned position.

###  2.2 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 $i$th site observation pair,
$O^{i}=(O_{a}^{x(i)},O_{b}^{y(i)})$, is associated with every aligned site $i$
in an alignment $M_{ab}$ of $p_{a}$ and $p_{b}$, where
$M_{ab}^{i}\in\left\\{\left(\begin{subarray}{c}x\\\
y\end{subarray}\right),\left(\begin{subarray}{c}x\\\
-\end{subarray}\right),\left(\begin{subarray}{c}-\\\
y\end{subarray}\right)\right\\}$ specifies the homology relationship at
position $i$ of the alignment (homologous, deletion with respect to $p_{a}$
and insertion with respect to $p_{a}$, respectively.), $i$ is taken to run
from 1 to $m$, $m$ is the length of the alignment $M_{ab}$, and
$x\in\\{1,\ldots,|p_{a}|\\}$ and $y\in\\{1,\ldots,|p_{b}|\\}$ specify the
indices of the positions in $p_{a}$ and $p_{b}$, respectively. $|p_{a}|$ and
$|p_{a}|$ specify the number of sites in $p_{a}$ and $p_{b}$, respectively.

Each site observation, $O_{a}^{x(i)}$ and $O_{b}^{y(i)}$, contains amino acid
and structural information corresponding to the two $C\alpha$ atoms at aligned
site $i$ belonging to each of the two proteins. A site observation
corresponding to a particular protein at aligned site $i$, $O_{a}^{x(i)}$, is
comprised of three different data types associated with the $C\alpha$ atom: an
amino acid ($A_{a}^{x(i)}$; discrete, one of twenty canonical amino acids),
$\phi$ and $\psi$ dihedral angles
($X_{a}^{x(i)}=\langle\phi_{a}^{x(i)},\psi_{a}^{x(i)}\rangle$; continuous,
bivariate), and a secondary structure classification ($S_{a}^{x(i)}$;
discrete, one of three classes: helix (H), sheet (S) or coil (C)). Therefore,
$O_{a}^{x(i)}=(A_{a}^{x(i)},X_{a}^{x(i)},S_{a}^{x(i)})$ and
$O_{b}^{y(i)}=(A_{b}^{y(i)},X_{b}^{y(i)},S_{b}^{y(i)})$.

###  2.3 Model structure

The sequence of hidden nodes in the HMM is written as
$H=(H^{1},H^{2},\ldots,H^{m}$). Each hidden node $H^{i}$ in the HMM
corresponds to a site observation pair, $O_{a}^{x(i)}$ and $O_{b}^{y(i)}$, at
an aligned site $i$ in the alignment $M_{ab}$. Initially we treat the
alignment $M_{ab}$ as given a priori, but later modify the HMM to marginalise
out an unobserved alignment.

The model is parametrised by $h$ hidden states. Every hidden node $H^{i}$
corresponding to an aligned site $i$ takes an integer value from $1$ to $h$
for the hidden state at node $H^{i}$. In turn, each hidden state specifies a
distribution over an evolutionary regime pair: $(r_{a}^{i},r_{b}^{i})$ as a
function of evolutionary time. An evolutionary regime pair consists of two
evolutionary regimes: $r_{a}^{i}$ and $r_{b}^{i}$. Each of the two
evolutionary regimes takes an integer value 1 or 2, i.e.
$(r_{a}^{i},r_{b}^{i})\in\\{(1,1),(1,2),(2,1),(2,2)\\}$. We return to the
specific role of the evolutionary regimes pairs in the next section.

The state of $H^{i}$ together with the evolutionary regime pair,
$(r_{a}^{i},r_{b}^{i})$, 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: $A^{i}=(A_{a}^{x(i)},A_{b}^{y(i)})$,
$X^{i}=(X_{a}^{x(i)},X_{b}^{y(i)})$ and $S^{i}=(S_{a}^{x(i)},S_{b}^{y(i)})$.
This conditional independence structure allows the likelihood of a site
observation pair at an aligned site $i$ to be written as follows:

$\displaystyle p(O^{i}|H^{i},r_{a}^{i},r_{b}^{i},t_{ab})=\overset{\text{amino
acid
evolution}}{\overbrace{p(A^{i}|H^{i},r_{a}^{i},r_{b}^{i},t_{ab})}}\times\overset{\text{dihedral
angle
evolution}}{\overbrace{p(X^{i}|H^{i},r_{a}^{i},r_{b}^{i},t_{ab})}}\times\overset{\text{secondary
structure evolution}}{\overbrace{p(S^{i}|H^{i},r_{a}^{i},r_{b}^{i},t_{ab})}}.$
(1)

The assumption of conditional independence provides computational
tractability, allowing us to avoid costly marginalisation when certain
combinations of data are missing (e.g. amino acid sequences present, but
secondary structures and dihedral angles missing).

###  2.4 Stochastic processes: modelling evolutionary dependencies

Each evolutionary regime couples together three time-reversible stochastic
processes that separately describe the evolution of the three pairs of
observation types, as in equation (1). Each evolutionary regime is intended to
capture different features of sequence and structural evolution. Parameters
that correspond to a particular evolutionary regime are termed regime-
specific, whereas parameters that are shared across all evolutionary regimes
are termed global.

#### Amino acid evolution

As is typical with models of sequence evolution, amino acid evolution,
$p(A_{a}^{x(i)},A_{b}^{y(i)}|H^{i},r_{a}^{i},r_{b}^{i},t_{ab})$, is described
by a Continuous-Time Markov Chain (CTMC). Each amino acid CTMC is parametrised
in the following way: the exchangeability of amino acids is described by a
$20\times 20$ symmetric global exchangeability matrix $S$ (190 free
parameters; Whelan and Goldman (2001)), a regime-specific set of 20 amino acid
equilibrium frequencies
$\Pi_{r}^{h}=\text{diag}\\{\pi_{1},\pi_{2},\ldots,\pi_{20}\\}$ (19 free
parameters per evolutionary regime) and a regime-specific scaling factor
$\Lambda_{r}^{h}$ (1 free parameter per evolutionary regime). Together these
parameters define a regime-specific time-reversible amino acid rate matrix
$Q_{r}^{h}=\Lambda_{r}^{h}S\Pi_{r}^{h}$. The stationary distribution of
$Q_{r}^{h}$ is given by the amino acid equilibrium frequencies: $\Pi_{r}^{h}$.

#### Secondary structure evolution

Secondary structure evolution,
$p(S_{a}^{x(i)},S_{b}^{y(i)}|H^{i},r_{a}^{i},r_{b}^{i},t_{ab})$, is also
described by a CTMC. In our model we use 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\times 3$ symmetric global exchangeability matrix $V$ and a regime-
specific set of 3 secondary structure equilibrium frequencies
$\Omega_{r}^{h}=\text{diag}\\{\pi_{1},\pi_{2},\pi_{3}\\}$. Together they
define a regime-specific time-reversible secondary structure rate matrix
$R_{r}^{h}=V\Omega_{r}^{h}$, with stationary distribution: $\Omega_{r}^{h}$.

#### Dihedral angle evolution

Central to our model is evolutionary dependence between dihedral angles,
$p(X_{a}^{x(i)},X_{b}^{y(i)}|H^{i},r_{a}^{i},r_{b}^{i},t_{ab})$. Typically,
the continuous-time evolution of the continuous-state random variables is
modelled by a diffusive process such as the Ornstein-Uhlenbeck (OU) process,
as in Challis and Schmidler (2012). However, an OU 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-Portugués et al. (2017).

Figure 2: Drift vector field for the WN diffusion with $A=(1,0.5;0.5,0.5)$,
$\mu=(0,0)$ and $\Sigma=(1.5)^{2}I$. The colour 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
colour 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 $\alpha_{3}$) between the $\phi$ and $\psi$ dihedral
angles.

Topologically, the WN diffusion (see Figure 2 for a pictorial example) can be
thought of as the analogue of the OU process on the torus
$\mathbb{T}^{2}=[-\pi,\pi)\times[-\pi,\pi)$. The WN diffusion arises as the
wrapping on $\mathbb{T}^{2}$ of the following Euclidean diffusion:

$\displaystyle dX_{t}=\overbrace{A\underset{k\in\mathbb{Z}^{2}}{\sum}(\mu-
X_{t}-2k\pi)w_{k}(X_{t})}^{\begin{subarray}{c}\text{drift}\\\
\text{coefficient}\end{subarray}}dt\;+\overbrace{\Sigma^{\frac{1}{2}}}^{\begin{subarray}{c}\text{diffusion}\\\
\text{coefficient}\end{subarray}}dW_{t},$ (2)

where $W_{t}$ is the two-dimensional Wiener process, $A$ is the drift matrix,
$\mu\in\mathbb{T}^{2}$ is the stationary mean, $\Sigma$ is the infinitesimal
covariance matrix and

$w_{k}(\theta)=\frac{\phi_{\frac{1}{2}A^{-1}\Sigma}(\theta-\mu+2k\pi)}{\underset{m\in\mathbb{Z}^{2}}{\sum}\phi_{\frac{1}{2}A^{-1}\Sigma}(\theta-\mu+2m\pi)},$
(3)

$\theta\in\mathbb{T}^{2}$, is a probability density function (pdf) for
$k\in\mathbb{Z}^{2}$. $\phi_{\Sigma}$ stands for the pdf of a bivariate
Gaussian $\mathcal{N}(0,\Sigma)$. The pdf (3) weights the linear drifts of (2)
such that they become smooth and periodic.

It is shown in García-Portugués et al. (2017) that the stationary distribution
of the WN diffusion is a $\mathrm{WN}(\mu,\Sigma)$, which has pdf:

$p_{\mathrm{WN}}(\theta|\mu,\Sigma)=\sum_{k\in\mathbb{Z}^{2}}\phi_{\Sigma}(\theta-\mu+2k\pi).$
(4)

Despite involving an infinite sum over $\mathbb{Z}^{2}$, 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 (2), yielding the following multimodal pseudo-tpd:

$\displaystyle\tilde{p}(\theta_{2}|\theta_{1},A,\mu,\Sigma,t)=\sum_{m\in\mathbb{Z}^{2}}p_{\mathrm{WN}}(\theta_{2}|\mu_{t}^{m},\Gamma_{t})w_{m}(\theta_{2}),$
(5)

with $\theta_{1},\theta_{2}\in\mathbb{T}^{2}$,
$\mu^{m}_{t}=\mu+e^{-tA}(\theta_{1}-\mu+2\pi m)$ and
$\Gamma_{t}=\int_{0}^{t}e^{-sA}\Sigma e^{-sA^{T}}ds$. The pseudo-tpd provides
a good approximation to the true tpd in key circumstances: i) $t\to 0$, since
it collapses in the Dirac’s delta; ii) $t\to\infty$, since it converges to the
stationary distribution; iii) high concentration, since the WN diffusion
becomes an OU process. Furthermore, it is shown in García-Portugués et al.
(2017) that the pseudo-tpd 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 $(X_{a}^{x(i)},X_{b}^{y(i)})$ is

$\displaystyle\tilde{p}(X_{b}^{y(i)}|X_{a}^{x(i)},A,\mu,\Sigma,t_{ab})p_{\mathrm{WN}}(X_{a}^{x(i)}|\mu,\tfrac{1}{2}A^{-1}\Sigma)=\tilde{p}(X_{a}^{x(i)}|X_{b}^{y(i)};A,\mu,\Sigma,t_{ab})p_{\mathrm{WN}}(X_{b}^{y(i)}|\mu,\tfrac{1}{2}A^{-1}\Sigma).$

Indeed, the WN diffusion is the unique time-reversible diffusion with the
stationary pdf (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 time-reversibility 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 marginalisation of
an unobserved ancestor.

The likelihood of a dihedral angle observation pair
$(X_{a}^{x(i)},X_{b}^{y(i)})$, assuming that $X_{a}^{x(i)}$ is drawn from the
stationary distribution, is given by:

$\displaystyle p(X_{a}^{x(i)},X_{b}^{y(i)}|H^{i},r_{a}^{i},r_{b}^{i},t_{ab})$
$\displaystyle=p(X_{a}^{x(i)},X_{b}^{y(i)}|A,\mu,\Sigma,t_{ab})$
$\displaystyle\approx\tilde{p}(X_{b}^{y(i)}|X_{a}^{x(i)},A,\mu,\Sigma,t_{ab})p_{\mathrm{WN}}(X_{a}^{x(i)}|\mu,\tfrac{1}{2}A^{-1}\Sigma),$
(6)

$A$ and $\Sigma$ are constrained to yield a covariance matrix $A^{-1}\Sigma$.
A parametrization that achieves this is
$\Sigma=\mathrm{diag}(\sigma_{1}^{2},\sigma_{2}^{2})$ and
$A=\big{(}\alpha_{1},\tfrac{\sigma_{1}}{\sigma_{2}}\alpha_{3};\tfrac{\sigma_{2}}{\sigma_{1}}\alpha_{3},\alpha_{2}\big{)}$,
$\alpha_{1}\alpha_{2}>\alpha_{3}^{2}$. $\alpha_{1}$ and $\alpha_{2}$ are the
drift components for the $\phi$ and $\psi$ dihedral angles, respectively.
Dependence (correlation) between the dihedral angles is captured by
$\alpha_{3}$. A depiction of a WN diffusion with given drift and diffusion
parameters is depicted in Figure 2.

The computationally tractable nature of the stationary density and tpd
required in (6) is key to enabling efficient training of and sampling under
the model. Specifically, the computation of (5) involves evaluating $e^{-tA}$
and $\Gamma_{t}$, which we can work out explicitly. First,
$e^{-tA}=a(t)I-b(t)A$ with $a(t)=e^{-rt}(\cosh(qt)+r\frac{\sinh(qt)}{q})$,
$b(t)=e^{-rt}\frac{\sinh(qt)}{q}$, $r=\tfrac{\mathrm{tr}(A)}{2}$ and
$q=\sqrt{|\det(A-rI)|}$. Second, since $A^{-1}\Sigma$ is symmetric,

$\Gamma_{t}=s(t)\tfrac{1}{2}A^{-1}\Sigma+i(t)\Sigma,$

with $s(t)=1-a(2t)$ and $i(t)=\tfrac{b(2t)}{2}$. This gives a neat
interpolation of the stationary and infinitesimal covariance matrices,
particularly convenient for efficiently evaluating (5) at different $t$’s.

###  2.5 Evolutionary regimes: modelling shift and drift

We now turn to the meaning of the evolutionary regime pairs. Two modes of
evolution are modelled: constant evolution and jump evolution. Constant
evolution occurs when the evolutionary regime starting in protein $p_{a}$ at
aligned site $i$, $r_{a}^{i}$, is the same as the evolutionary regime ending
in protein $p_{b}$ at aligned site $i$, $r_{b}^{i}$, i.e.
$r_{a}^{i}=r_{b}^{i}$. Conversely, jump evolution occurs when $r_{a}^{i}\neq
r_{b}^{i}$. Constant evolution is intended to capture angular drift (changes
in dihedral angles localised to a region of the Ramachandran plot), whereas
jump evolution is intended 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
the evolutionary regime pairs:

$p(r_{a}^{i},r_{b}^{i}|H^{i},t_{ab})=p(r_{a}^{i}|H^{i},r_{b}^{i},t_{ab})p(r_{b}^{i}|H^{i}),$
(7)

where

$\displaystyle
p(r_{a}^{i}|H^{i},r_{b}^{i},t_{ab})=\begin{cases}e^{-\gamma_{H^{i}}t_{ab}}+\pi_{H^{i},r_{b}^{i}}(1-e^{-\gamma_{H^{i}}t_{ab}}),&\text{if
}r_{a}^{i}=r_{b}^{i},\\\
\pi_{H^{i},r_{b}^{i}}(1-e^{-\gamma_{H^{i}}t_{ab}}),&\text{if }r_{a}^{i}\neq
r_{b}^{i},\\\ \end{cases}$

and $p(r_{a}^{i}|H^{i})=\pi_{H^{i},r_{a}^{i}}$ and
$p(r_{b}^{i}|H^{i})=\pi_{H^{i},r_{b}^{i}}$. $\pi_{H^{i},r_{a}^{i}}$ and
$\pi_{H^{i},r_{b}^{i}}$ are model parameters specifying the probability of
starting in regime $r_{a}^{i}$ or $r_{b}^{i}$, respectively, corresponding to
the hidden state specified by node $H^{i}$. $\gamma_{H^{i}}>0$ is a model
parameter specifying the jump rate corresponding to the hidden state specified
by node $H^{i}$.

The regime pair jump probabilities have been chosen so that time-reversibility
holds, in other words:

$\displaystyle
p(r_{a}^{i}|H^{i},r_{b}^{i},t_{ab})p(r_{b}^{i}|H^{i})=p(r_{b}^{i}|H^{i},r_{a}^{i},t_{ab})p(r_{a}^{i}|H^{i}).$

The hidden state at node $H^{i}$, together with a regime pair
$(r_{a}^{i},r_{b}^{i})$ and the evolutionary time $t_{ab}$, specifies the
joint likelihood over site observation pairs:

$\displaystyle p$
$\displaystyle(O_{a}^{x(i)},O_{b}^{y(i)}|H^{i},r_{a}^{i},r_{b}^{i},t_{ab})=\begin{cases}p(O_{a}^{x(i)},O_{b}^{y(i)}|H^{i},r_{c}^{i},t_{ab}),&\text{if
}r_{a}^{i}=r_{b}^{i}=r_{c}^{i},\\\
p(O_{a}^{x(i)}|H^{i},r_{a}^{i})p(O_{b}^{y(i)}|H^{i},r_{b}^{i}),&\text{if
}r_{a}^{i}\neq r_{b}^{i}.\\\ \end{cases}$ (8)

In the case of constant evolution, evolution at aligned $i$ is described in
terms of the same evolutionary regime $r_{c}^{i}$. 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 jump evolution, the evolutionary processes are, after the
evolutionary jump, restarted independently in the stationary distribution of
the new evolutionary regime. Thus the site observations $O_{a}^{x(i)}$ and
$O_{b}^{y(i)}$ are assumed to be drawn from the stationary distributions of
two separate stochastic processes corresponding to evolutionary regimes
$r_{a}^{i}$ and $r_{b}^{i}$, respectively. This implies that, conditional on a
jump, the likelihood of the observations is no longer dependent on $t_{ab}$.
Furthermore, there is no longer an evolutionary trajectory linking the two
site observations, hence there is no need to perform a computationally
expensive marginalisation over all possible trajectories, as might be
necessary in a model with continuous-time Markovian switching between
evolution regimes. The likelihood of an observation pair is now simply a sum
over the four possible regime pairs:

$\displaystyle
p(O^{x(i)},O^{y(i)}|H^{i},t_{ab})=\sum_{(r_{a}^{i},r_{b}^{i})\in{R}}p(O^{x(i)},O^{y(i)}|H^{i},r_{a}^{i},r_{b}^{i},t_{ab})p(r_{a}^{i},r_{b}^{i}|H^{i},t_{ab}),$

where $R=\\{(1,1),(1,2),(2,1),(2,2)\\}$ is the set of four regime pairs,
$p(O^{x(i)},O^{y(i)}|H^{i},r_{a}^{i},r_{b}^{i})$ is given by (8) and
$p(r_{a}^{i},r_{b}^{i}|H^{i},t_{ab})$ is given by (7).

###  2.6 Identification of evolutionary motifs encoding jump evolution

In order to identify aligned sites having potential evolutionary motifs
encoding jump evolution, 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,
$(A_{a},A_{b},X_{a},X_{b})$. Homologous sites corresponding to a single hidden
state and with evidence of jump evolution ($r_{a}^{i}\neq r_{b}^{i}$) at
posterior probability $>0.90$ were identified, that is, the $i$’s such that
$p(H^{i},r_{a}^{i}\neq r_{b}^{i}|A_{a},A_{b},X_{a},X_{b})>0.90$.

In a second filtering step, amino acid sequences and a single set of dihedral
angles corresponding to one of the proteins were used ($A_{a},A_{b},X_{a}$ or
$A_{a},A_{b},X_{b}$) to infer the posterior probability, this time at a lower
threshold: $p(H^{i},r_{a}^{i}\neq r_{b}^{i}|A_{a},A_{b},X_{a})>0.50$ or
$p(H^{i},r_{a}^{i}\neq r_{b}^{i}|A_{a},A_{b},X_{b})>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 further
downstream analysis.

###  2.7 Statistical alignment: modelling 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
for indels, a modified pairwise TKF92 alignment HMM based on Miklós et al.
(2004) was implemented. The TKF92 alignment HMM was augmented with additional
evolutionary hidden states intended to capture local sequence and structure
evolutionary dependencies. Furthermore, it was modified such that neighbouring
dependencies amongst hidden states at adjacent alignment sites were modelled.

Whilst it is possible to fix the alignment in advance by pre-aligning the
sequences using one of the many available alignment methods (Katoh et al.
(2002); Edgar (2004)) or using a curated alignment (such as from the HOMSTRAD
database), doing so ignores alignment uncertainty. For more detail we refer to
the supplementary material.

###  2.8 Training and test dataset

A training dataset of 1200 protein pairs (2400 proteins; 417,870 site
observation pairs) and a test dataset of 38 protein pairs (76 proteins; 14125
site observation pairs) were assembled from 1032 protein families in the
HOMSTRAD database.

For each protein family in HOMSTRAD (ranging in size from 2 to 22 homologous
proteins each), a phylogenetic tree was inferred from the HOMSTRAD protein
family sequence alignment using FastTree (Price et al., 2010). Each protein
family tree was taken and protein pairs selected such that the sum of the
branches between pairs was maximised, whilst ensuring that no pair of proteins
in the set shared an overlapping evolutionary history. This was done in order
to maximise the amount of information, whilst minimising dependencies due to
shared evolutionary history.

Dihedral angles were computed from the PDB coordinates of each protein
structure using the BioPython.PDB package (Hamelryck and Manderick, 2003).
Furthermore, each protein was taken and the corresponding full length protein
structure obtained from the PDB database and the secondary structure annotated
at each amino acid position using DSSP (Touw et al., 2015). For some proteins
we were not able to obtain the full length PDB structures, which may result in
misclassification of secondary structure interactions due to truncation. For
these structures the secondary structure annotations were treated as missing.
All 38 proteins in the test dataset had secondary structure annotations and
were from distinct protein families.

###  2.9 Model training

Maximum likelihood estimation of the model parameters, $\hat{\Psi}$, was done
using Stochastic Expectation Maximisation (StEM). Forward Filtering Backward
Sampling (FFBS) was used to jointly sample alignment configurations
($M_{ab}$), hidden node states ($H$) and evolutionary regimes,
($r_{a},r_{b}$). The Metropolis-Hastings algorithm was used to sample the four
pair-specific continuous parameters
$\theta_{ab}=\\{t_{ab},\lambda_{ab},\mu_{ab},r_{ab}\\}$: evolutionary time
($t_{ab}$), insertion ($\lambda_{ab}$), deletion ($\mu_{ab}$) and geometric-
extension ($r_{ab}$) rates. In other words, at iteration $k$ for each pair of
unaligned observation sequences $O_{a}$ and $O_{b}$ we draw samples, from the
following joint-distribution:

$\displaystyle Z_{ab}^{(k)}\sim
p(M_{ab},H,r_{a},r_{b},\theta_{ab}|O_{a},O_{b},\Psi^{(k)}).$

In the M-step the samples from the previous E-step, were used to update the
hidden node parameters ($\hat{\Psi}$) using efficient sufficient statistics
and COBYLA optimization algorithm (Powell, 1994) in the NLOpt library
(Johnson, 2014) was used.

###  2.10 Angular distances

For benchmarking purposes, the angular cosine distance was used to measure
distances between pairs of dihedral angles, $\langle\phi_{a},\psi_{a}\rangle$
and $\langle\phi_{b},\psi_{b}\rangle$. It is defined as follows (Downs and
Mardia, 2002):

$\displaystyle d$
$\displaystyle(\langle\phi_{a},\psi_{a}\rangle,\langle\phi_{b},\psi_{b}\rangle)=\sqrt{4-2\cos(\phi_{a}-\phi_{b})-2\cos(\psi_{a}-\psi_{b})}.$
(9)

The maximum possible distance is $\sqrt{8}\approx 2.828$. It has the property
that when $\phi_{a}-\phi_{b}\approx 0$ and $\psi_{a}-\psi_{b}\approx 0$ are
near zero it may be approximated by the Euclidean distance – using the small
angle approximation for cosine ($\cos\theta\approx 1-\frac{\theta^{2}}{2}$
when $\theta$ is near zero):

$\displaystyle
d(\langle\phi_{a},\psi_{a}\rangle,\langle\phi_{b},\psi_{b}\rangle)\approx\sqrt{4-2(1-(\phi_{a}-\phi_{b})^{2}/2)-2(1-(\psi_{a}-\psi_{b})^{2}/2)}=\sqrt{(\phi_{a}-\phi_{b})^{2}+(\psi_{a}-\psi_{b})^{2}}.$

##  3 Results and discussion

###  3.1 Selecting the number of hidden states

Models with 8, 16, 32, 48, 64, 80, 96 and 112 hidden states were trained until
convergence for three different repetitions using different initial random
number seeds. The highest log-likelihood model of the three repetitions for
each number of hidden states was selected for downstream analysis.

Following that, marginal likelihoods $p(D|\text{model})$, i.e. model evidence
for each of the 38 protein pairs in the test dataset, were computed under each
model by fixing the alignments to the respective HOMSTRAD alignments. The
alignments were fixed a priori in order to make computation of the marginal
likelihoods computationally tractable. Additionally, predictive accuracies
under a homology modelling scenario,
$p(X_{b}|A_{a},A_{b},X_{a},\text{model})$, were calculated for the same 38
protein pairs.

The 112 hidden state model had the highest total marginal log-likelihood for
protein pairs in the test dataset and predictive accuracy under the homology
modelling scenario comparable to the predictive accuracies of the 16, 32, 48,
64, 80 and 96 hidden state models (all were with standard error of the mean,
see supplementary data). Only the model with 8 hidden states had significantly
poorer predictive accuracy on the test dataset.

Figure 3: 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.

Although the model with 112 hidden states had the highest total marginal log-
likelihoods, the model with 64 hidden states was selected as the final model.
A 64 hidden state model was chosen as a trade-off between total marginal log-
likelihood on the test dataset, computational time and the number of model
parameters. Inference under the model scales $\mathcal{O}(h^{2})$ in the
number of hidden states and the higher total marginal log-likelihood of the
112 hidden state model did not justify the increase in computational time when
sampling from the model.

###  3.2 Stationary distributions over dihedral angles capture the empirical
distribution

Figure 3 illustrates the sampled and empirical dihedral angle distributions.
There is a good correspondence between dihedral angles sampled under the model
(Figure 3, left) and the empirical distribution of dihedral angles in our
training dataset (Figure 3, 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, however, it is a good indication that there are a lack of obvious
errors in the model implementation.

###  3.3 Estimates of evolutionary time from dihedral angles are consistent
with estimates from sequence

Figure 4 compares evolutionary times estimated using only pairs of homologous
amino acid sequences only versus pairs of homologous dihedral angles only. 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.

Figure 4: Scatterplot comparing evolutionary times estimated using pairs of
homologous amino acid sequences only versus pairs of homologous sets of
dihedral angles only for $N=38$ proteins pairs in the test dataset. The
$x$-coordinate of each point gives the estimated evolutionary time based only
on the amino acid sequence, whereas the $y$-coordinate gives the estimated
evolutionary time based only on the dihedral angles. The diagonal line
represents $y=x$.

A paired $t$-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 vs. 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 under-estimate or over-estimate the evolutionary times when either
sequence or dihedral angles are used.

Interestingly, the variance in the sampled evolutionary times is higher when
dihedral angles only are used, as compared to sequence only (see Figure 11).

###  3.4 The relationship between evolutionary time and angular distance is
adequately modelled

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 (Figure
5).

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 $\sqrt{8}\approx 2.828$.

Figure 5: Evolutionary time vs. angular distances between real and
corresponding sampled proteins pairs in the training dataset at 50
representative evolutionary times. Mean angular distances between the real
dihedral angles in a protein pair (red, $X_{a}$ and $X_{b}$) and sampled
dihedral angles in a sampled protein pair (blue, $X_{a}$ and $\hat{X}_{b}$)
were compared to test how well the sampled dihedral angles reproduced the real
angular distances. The dihedral angles ($\hat{X}_{b}$) of each sampled protein
pair were sampled by conditioning on both amino acid sequences and the
homologous dihedral angles ($A_{a}$, $A_{b}$, $X_{a}$), and the estimated
evolutionary time ($\hat{t}_{ab}$) for the real protein pair. The regression
curves were obtained by a quadratic LOcally-weighted regrESSion (LOESS), with
smoothing parameter chosen by leave-one-out cross-validation. The 95%
confidence intervals for the mean assume error normality.

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
Figure 5), 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.

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 under-estimate angular distances. This likely
reflects the fact that the tpd of the WN diffusion specified is localised
around its mean when $t_{ab}\to\infty$ and 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 is still somewhat limited in its flexibility, as it only
allows jumps between two evolutionary regimes.

###  3.5 Evaluation of the model

The conditional independence structure in (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. The average angular distance (9) between the sampled and
known dihedral angles was used as the measure of predictive accuracy.

Figure 6: Cartoon structure representations of E.coli
glyceraldehyde-3-phosphate dehydrogenase structure (PDB 1gad) are depicted in
each panel, overlaid with predictive accuracy when using different
combinations of observed data to predict missing dihedral angles in 1gad.
Thermus aquaticus glyceraldehyde-3-phosphate dehydrogenase (PDB 1cer) was used
as a homolog for the purposes of prediction. Predictive accuracy is indicated
using a colour gradient depicting the mean angular distance between the true
dihedral angle ($X_{\text{1gad}}^{i}$) and the predicted (sampled) dihedral
angles ($\hat{X}_{\text{1gad}}^{i}$) at each amino acid position. The label at
the bottom of each panel indicates the data combination used. In A, no data
was used for prediction. In B, only the amino acid sequence corresponding to
1gad ($A_{\text{1gad}}$) was used. In C, the amino acid sequence of 1gad
($A_{\text{1gad}}$) and the amino acid sequence of the homologous protein
($A_{\text{1cer}}$) were used. In D, both amino acids sequences
($A_{\text{1cer}}$ and $A_{\text{1gad}}$) and the secondary structure of the
homologous protein ($S_{\text{1cer}}$) were used. In E, both the amino acid
sequences ($A_{\text{1cer}}$ and $A_{\text{1gad}}$) and the dihedral angles of
the homologous protein ($X_{\text{1cer}}$) were used. Finally, in panel F the
same combination of observations was used as in E, but the alignment was
treated as known a priori.

Figure 6 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 7 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 the key patterns
identified in Figures 6 and 7 here.

Figure 7: Benchmarks of predictive accuracy (measured using angular distance,
lower is better) on a random subset of ten protein pairs in the test dataset,
giving a representative view of predictive accuracy under six different
combinations of observations. The dihedral angles $X_{b}$ of protein $b$ were
treated as missing and were sampled under the model, whereas $p_{a}$ was a
homologous protein used for the purposes of prediction. See the legend of
Figure 6 for a description of each combination (A-F). The final set of bars,
denoted ‘Mean (N=38)’, are the mean values for the entire test dataset of N=38
protein pairs. The error bars are the standard errors.

Combination 1 refers to random sampling from the model, implying no data
observations were conditioned on besides the respective lengths of proteins
$p_{a}$ and $p_{b}$. The average angular distance between the true and
predicted dihedral angles was 1.6. Random sampling acts as a baseline for
predictive accuracy. It is apparent from Figure 6 that the model has a
propensity to predict right-handed $\alpha$-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 6 and 7 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 modelling 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) 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
marginalisation 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 correctness
of the alignment HMM.

The alignment HMM is valuable accounts for alignment uncertainty in a
principled manner, this is particularly useful when an appropriate alignment
is unavailable. However, it should be noted that inference scales
$\mathcal{O}(|p_{a}||p_{b}|h^{2})$ when treating the alignment as unobserved,
where $|p_{a}|$ and $|p_{b}|$ are the lengths of $p_{a}$ and $p_{b}$,
respectively. Inference scales $\mathcal{O}(mh^{2})$ when the alignment is
fixed a priori, where $m$ is the length of alignment $M_{ab}$ and is typically
much smaller than $|p_{a}||p_{b}|$.

It should be emphasised that we do not expect ETDBN to compete with structure
prediction packages such as Rosetta (Rohl et al., 2004) or homology modelling
software such as Arnold et al. (2006) in terms of predictive accuracy. Our
current model is a local model of structure evolution – it is not even
expected capture constraints such as the radius of gyration of a protein or
other global features typical of proteins.

###  3.6 Evolutionary hidden states reveal a common evolutionary motif

One benefit of ETDBN is that the 64 evolutionary hidden states learned during
the training phase are interpretable. We give an example of a hidden state
encoding jump evolution that was subsequently found to represent an
evolutionary motif present in a large number of protein pairs in our test and
training datasets.

Evolutionary hidden state 3 (Figure 8) was selected from the 64 hidden states
as an example of a hidden state encoding jump evolution and capturing angular
shift (a large change in dihedral angle). A notable feature of this hidden
state is that the change in dihedral angles between evolutionary regimes
$r_{1}$ and $r_{2}$ is associated with specific amino acid changes. In regime
$r_{1}$ the amino acid frequencies are relatively spread out amongst a number
of amino acids, whereas in regime $r_{2}$ the frequencies are particularly
concentrated in favour of glycine (Gly) and asparagine (Asp), with glycine
being significantly more probable in regime $r_{2}$ than $r_{1}$. This
suggests that conditioned on hidden state 3, an exchange between a glycine to
another amino acid is likely indicative of a jump and hence a corresponding
change in dihedral angle. This is consistent with what we find in a subsequent
analysis of evolutionary motifs. This particular jump occurs in coil regions.

Having selected hidden state 3, positions in 238 protein pairs were analysed
for evidence of the corresponding evolutionary motif. 38 protein pairs in the
test dataset and a further 200 from the training dataset were analysed using
the criteria described in the Methods section. Using the first criterion, 84
protein sites in 59 protein pairs corresponding to $H^{i}=3$ (evolutionary
hidden state 3) were identified. Of the 84 protein sites, 34 protein sites met
the second criterion.

We give an example of a homologous protein pair illustrating the identified
evolutionary motif. Two histidine-containing phosphocarriers, 1pch (M.
capricolum) and 1poh (E. coli), were identified as having the evolutionary
motif (Figure 9) at homologous site E39/G39.

Figure 8: Depiction of evolutionary hidden state 3. This hidden state was
sampled at $0.69\%$ of sites (the average was $1.56\%$). The equilibrium
frequencies of $r_{1}$ and $r_{2}$ were ${\pi_{1}}=0.691$ and
${\pi_{2}}=0.309$, respectively. The jump rate was $\gamma=31.76$. The
corresponding regime pair probabilities are depicted to the right as a
function of evolutionary time. In the main figure, the two rows depict the
parameters encoded by the two evolutionary regimes, respectively. Columns 1
and 3 depict the parameters governing the the amino acid and secondary
structure stochastic processes, respectively. The secondary structure classes
correspond to H=helix, S=sheet and C=coil. Column 2 depicts the WN diffusions.
The stationary distributions of the WN diffusions are shown using black
contour lines, the direction of the drifts are indicated by the arrows and the
magnitude of the drifts at each position indicated using the colour gradient.

Most positions in the homologous pair have very low jump probabilities
($\approx 0.0$), with the exception of positions N38/N38 and E39/G39, which
both have high posterior jump probabilities ($\approx 1.0$). The exchange
between a glutamate (at position 39 in 1poh) and a glycine (at position 39 in
1pch) appears to be responsible for the shift in dihedral angle. This exchange
corresponds to a significant jump in dihedral angle
($\langle\phi_{\text{1poh,E39}}=-1.63,\psi_{\text{1pch,E39}}=-0.06\rangle\rightarrow\langle\phi_{\text{1poh,G39}}=1.40,\psi_{\text{1pch,G39}}=0.22\rangle$).
The angular distance between the two dihedral angles is 2.01. This is
consistent with the amino acid frequency parameters specified by the two
regimes for hidden state 3 (Figure 8).

Regime $r_{1}$ indicates that a number of amino acids (alanine, aspartic acid,
glycine, histidine, lysine, asparagine, proline, glutamine, arginine, serine
and theorine) other than glutamate plausibly coincide with the particular
dihedral angle conformation specified by regime $r_{1}$. The involvement of
glycine in a jump is not surprising as it is a small and flexible amino acid,
whereas the role of asparagine is less clear. In our analysis of 238 protein
pairs we found that of the seven positions meeting the criteria for hidden
state 3 and involving an exchange with asparagine (Asn), four were an exchange
between an asparagine and a glycine, whereas the remaining three were between
asparagine and one of lysine, histidine or serine.

Figure 9: Depiction of two histidine-containing phosphocarriers, PDB 1pch and
1poh, superimposed. On the left is a cartoon representation of the two
proteins corresponding to regions F29-K49 and F29-A42, respectively, with
posterior jump probabilities at each position overlaid. On the right is a
ball-and-stick representation giving atomic detail for a smaller region
(I36-G42 and T36-A42, respectively). The exchange between a glutamate (E39 in
1poh) and a glycine (G39 in 1pch) is associated with a large change in
dihedral angle as indicated by the curved arrows.

###  3.7 Using dihedral angles for alignment

A valuable feature of our model is its ability to account for alignment
uncertainty by summing over possible pairwise alignments using the TKF92 model
as a prior distribution over indel histories, whilst simultaneously taking
into account neighbouring dependencies amongst aligned sites. Doing so results
in a sample of alignments rather than a single alignment. Nevertheless, a
single Maximum A Posteriori (MAP) pairwise alignment may be obtained from the
alignment samples and used for downstream analysis.

Figure 10: Alignment benchmarks. A: averages of alignment similarity across
different methods and combinations of data observations, where the
observations were simulated from the ETDBN model conditioned on a single
alignment, hence the underlying alignment was known. B: averages of alignment
similarity, where the 38 HOMSTRAD alignments in the test dataset were taken as
the true alignments. C: averages of precision in predicting homologous site
pairs in 38 sequence pairs from the test dataset, where homologous site pairs
in the HOMSTRAD alignments were taken to be the true homologous site pairs.
Note: combination 11 in panels B and C denotes that the StatAlign structural
alignment method (Herman et al., 2014) was used to align the two proteins
based on the $C\alpha$ coordinates.

ETDBN and several other alignment methods (namely StatAlign, BAli-Phy, MUSCLE
and MAFFT) were used to infer pairwise alignments from simulated and real data
under various combinations of data observations, for example: an amino acid
sequence pair ($A_{a},A_{b}$), a secondary structure sequence pair
($S_{a},S_{b}$), a dihedral angle sequence pair ($X_{a},X_{b}$) and
combinations thereof.

In the first set of benchmarks (Figure 10A), pairs of proteins were simulated
from the ETDBN model conditioned on 38 different pairwise alignments and
corresponding evolutionary times. This resulted in a set of 38 simulated
pairwise alignments together with corresponding observations, implying that
the true underlying alignments were known for each of the simulated protein
pairs. ETDBN and a number other alignment methods were used to infer pairwise
alignments for each. The alignment similarity metric (Schwartz et al., 2005)
was used to measure the similarity between the inferred alignments and the
true alignments, where higher similarity indicates better predictions. It was
found that, when using the simulated amino acid sequences alone, ETDBN (10A.5)
outperformed all four other methods tested (10A.1 MUSCLE, 10A.2 MAFFT, 10A.3
StatAlign, 10A.4 BAli-Phy). The greater performance of ETDBN compared to other
methods can not be considered a fair comparison, as the data were simulated
under the ETDBN model.

More revealing in Figure 10A was the alignment similarity under ETDBN when
using different combinations of simulated data observations. It was found that
secondary structure alone (10A.6) performed the worst, which is unsurprising
given that only three states were available to align the proteins. The second
worst in terms of alignment similarity was amino acid sequences alone (10A.5),
followed by amino acid sequences and secondary structures (10A.7).
Interestingly, using dihedral angles only (10A.8) outperformed both 10A.5
(sequences only) and 10A.6 (secondary structures only). Finally, using amino
acid sequence together with dihedral angles (10A.9) or all three data types
combined (10A.10) outperformed all other combinations. This illustrates that,
at least under simulation conditions, increasing the number of data
observations results in better alignment accuracy.

Following that, the various alignment methods were benchmarked against 38
pairwise alignments consisting of real sequence and structure observations in
the test dataset. These pairwise alignments were obtained from the HOMSTRAD
alignments. The sequence identity of these pairwise alignments ranged from
$10\%$ to $93\%$, with an average sequence identity of $39\%$. When
benchmarking the MAP estimated alignments against the HOMSTRAD alignments
(Figure 10B), using real sequences alone for inference $(A_{a},A_{b})$, ETDBN
(10B.5) had a similar degree of accuracy when compared to several other
sequence-based methods (10B.1 StatAlign, 10B.2 BaliPhy, 10B.3 MUSCLE and 10B.4
MAFFT). This demonstrates that ETDBN has performance comparable to that of
other commonly-used sequence alignment methods.

Using $(S_{a},S_{b})$ alone, ETDBN (10B.6) had substantially lower alignment
similarity compared to sequence only, which was expected given that a similar
result was obtained for the simulated data (10A.6). However, when including
the real sequences (10B.7) the predictive accuracy was once again comparable
to sequence only inferences (10B.1–10B.5).

When using ($X_{a}X_{b}$) alone (10B.8), the alignment similarity was found to
be somewhat worse than the sequence only cases. Furthermore, when introducing
the sequences (10.9) and secondary structures (10B.10) in addition to the
dihedral angles, the similarity still remained worse than the sequence only
methods (10B.1–10B.5), despite the additional information. These results are
in contrast to the results we obtained for simulated data (10A.8–10A.10).

Upon further investigation, the trend was found to reverse (Figure 10C) when
the precision of predicting homologous sites was calculated (the fraction of
sites which were predicted as homologous and were correctly predicted as
such). Therefore when only dihedral angle observations are used, ETDBN
underpredicts the number of homologous sites, however, when a homologous site
is predicted, it is correctly predicted more often than when using only amino
acid sequences. In particular, ETDBN predicted fewer homologous sites with
coiled secondary structure compared to homologous sites with helical or sheet
secondary structure. This pattern of results may be in part due to the WN
diffusion used to model evolution of dihedral angles. The WN diffusion is
suitable for modelling angular drift (small changes in angle localised around
a region of the Ramachandran plot) but does not sufficiently capture angular
shift (large changes in angle between regions of the Ramachandran plot, which
are more likely in coiled regions). The jump model to some extent captures
angular shift, but it can only do so for a single evolutionary jump between
two regions of the Ramachandran plot. In other words, it does not consider
multiple different homologous angular shift events caused by amino acid
mutation.

When interpreting these results it important to note the HOMSTRAD alignments
should not be considered the true underlying alignments and may even be
strongly biased (e.g. they may favour the most parsimonious alignments, with
the fewest number of indels). In practice, it is extremely difficult to obtain
the true underlying alignment, because it would require an experiment where
every indel event since the common ancestor is observed, a seeming impossible
task outside of simulation or laboratory conditions.

##  4 Concluding remarks

The main achievement of this work is a computationally tractable, generative
and interpretable probabilistic model of protein sequence and structure
evolution on a local scale.

Previous stochastic models of protein sequence and structure evolution
emphasised estimation of evolutionary parameters (Challis and Schmidler, 2012;
Herman et al., 2014). ETDBN is somewhat of a departure from these previous
models, but is likewise capable of estimating evolutionary parameters. We show
that estimates of evolutionary times inferred under ETDBN are consistent
regardless of whether amino acid sequence or dihedral angle observations are
used. In addition, the relationship between evolutionary time and angular
distance in real proteins is adequately recapitulated in protein pairs sampled
under the model, albeit the angular distance is under-estimated for larger
evolutionary times, which might be explained by the limited flexibility of the
jump model. Like previous models, ETDBN is capable of dealing with alignment
uncertainty by marginalising over indel histories; it predicts pairwise MAP
consensus alignments with accuracy similar to that of score-based and
statistical alignment methods.

The generative nature of ETDBN allows us to demonstrate that the underlying
empirical distributions over dihedral angles (depicted using Ramachandran
plots) are captured and that the model is capable of predicting missing
observations, such as dihedral angles, from a variety of different data types.
For example, an amino acid sequence, a homologous amino acid sequence, a
homologous secondary structure, a homologous set of dihedral angles or any
combination thereof.

Due to its local nature, ETDBN does not constitute a homology modelling method
in itself. Rather, it can be used as a building block, much like fragment
libraries model local structure in protein structure prediction methods. ETDBN
places the homology modelling problem on a statistical footing, enabling a
number of approaches to later be used, such as multi-level modelling, i.e.
combining fine-grained distributions (for example, distributions over dihedral
angles, such as ETDBN) and coarse-grained distributions (for example,
distributions describing the global properties of proteins, such as
compactness).

In addition to multi-level modelling, probabilistic models such as ETDBN allow
one to account for and to make statements about uncertainty (e.g. with respect
to evolutionary time, alignment, etc.) in a rigorous manner. One immediate use
of ETDBN from a structure prediction or homology modelling perspective is as
an efficient proposal distribution. ETDBN could be used to sample protein
structures (possibly conditioned on various data observations) in a highly
computationally efficient manner, such that the resulting samples are expected
to be located in regions of high probability density with respect to the true
underlying distribution.

A final key feature of our evolutionary model is its interpretable nature.
This interpretability enables the identification of potential evolutionary
motifs – common patterns of sequence-structure evolution. We identify one such
evolutionary motif in 34 different homologous protein pairs. A major direction
for future research is the further identification of such evolutionary motifs.
Understanding these evolutionary motifs, may i) improve homology modelling
predictions; ii) provide more accurate estimates of evolutionary parameters;
and iii) produce better models of protein evolution that more realistically
capture evolutionary trajectories through sequence and structure space, which
may help identify functionally relevant positions that may be possible drug
targets.

##  5 Future challenges

For reasons of computational tractability the implemented model is pairwise,
but it is theoretically possible to generalise it to a phylogeny, such as in
Herman et al. (2014). In practice, for three or more sequences on a phylogeny
it is necessary to marginalise out the unobserved ancestral protein states in
order to compute likelihoods. Felsenstein’s algorithm can be used to
marginalise over discrete ancestral states, such as amino acids in a
computationally efficient manner. However, we have not yet established whether
a similar efficient algorithm exists for marginalising the continuous
ancestral dihedral angle states under the WN diffusion, thereby necessitating
a more expensive MCMC algorithm. A possibly greater computational hindrance to
considering a phylogeny is the alignment problem, which scales
$\mathcal{O}(l_{1}\times{l_{2}}\times\ldots\times{l_{N}})$, where $l_{i}$ is
the length of sequence $i$ and $N$ is the number of sequences, although MCMC
approaches are also possible (Herman et al., 2014).

Although we believe our model provides a substantial improvement over current
stochastic models of sequence and structural evolution, there is still scope
for improvement. The WN diffusions used to model dihedral angle evolution
adequately capture angular drift (small local changes in dihedral angle), but
are less capable of capturing angular shift (large changes in dihedral angle).
This problem is somewhat mitigated by our jump model, but even this only
allows jumps between two possible evolutionary regimes conditioned on a
particular hidden state. A more realistic model would allow an arbitrary
number of switches between evolutionary regimes together with neighbouring
dependencies amongst adjacent sites along the evolutionary trajectory.
Continuous-Time Bayesian Networks (Nodelman et al., 2002) provide such a
framework for large state space models with certain sparsity conditions, but
will likely incur a large penalty in terms of computational efficiency.

## Software availability

Julia code (tested on both Windows and Linux platforms) is available at:

http://www.computingforbiology.org/software/etdbn

## References

  * Arnold et al. (2006) Arnold, K., Bordoli, L., Kopp, J., and Schwede, T. 2006.  The SWISS-MODEL workspace: a web-based environment for protein structure homology modelling.  Bioinformatics, 22(2): 195–201.
  * Boomsma et al. (2008) Boomsma, W., Mardia, K. V., Taylor, C. C., Ferkinghoff-Borg, J., Krogh, A., and Hamelryck, T. 2008.  A generative, probabilistic model of local protein structure.  Proceedings of the National Academy of Sciences, 105(26): 8932–8937.
  * Challis and Schmidler (2012) Challis, C. J. and Schmidler, S. C. 2012.  A stochastic evolutionary model for protein structure alignment and phylogeny.  Molecular biology and evolution, 29(11): 3575–3587.
  * Downs and Mardia (2002) Downs, T. D. and Mardia, K. 2002.  Circular regression.  Biometrika, 89(3): 683–698.
  * Felsenstein (1985) Felsenstein, J. 1985.  Phylogenies and the comparative method.  American Naturalist, 125(1): 1–15.
  * Frellsen et al. (2012) Frellsen, J., Mardia, K. V., Borg, M., Ferkinghoff-Borg, J., and Hamelryck, T. 2012\.  Towards a general probabilistic model of protein structure: the reference ratio method.  In Bayesian methods in structural bioinformatics, pages 125–134. Springer.
  * García-Portugués et al. (2017) García-Portugués, E., Sørensen, M., Mardia, K. V., and Hamelryck, T. 2017\.  Langevin diffusions on the torus: estimation and applications.  arXiv:1705.00296.
  * Grishin (2001) Grishin, N. V. 2001.  Fold change in evolution of protein structures.  Journal of structural biology, 134(2): 167–185.
  * Hamelryck and Manderick (2003) Hamelryck, T. and Manderick, B. 2003.  Pdb file parser and structure class implemented in python.  Bioinformatics, 19(17): 2308–2310.
  * Herman et al. (2014) Herman, J. L., Challis, C. J., Novák, Á., Hein, J., and Schmidler, S. C. 2014.  Simultaneous bayesian estimation of alignment and phylogeny under a joint model of protein sequence and structure.  Molecular biology and evolution, 31(9): 2251–2266.
  * Johnson (2014) Johnson, S. G. 2014.  The NLopt nonlinear-optimization package.
  * Miklós et al. (2004) Miklós, I., Lunter, G., and Holmes, I. 2004.  A “long indel” model for evolutionary sequence alignment.  Molecular Biology and Evolution, 21(3): 529–540.
  * Mizuguchi et al. (1998) Mizuguchi, K., Deane, C. M., Blundell, T. L., and Overington, J. P. 1998.  Homstrad: a database of protein structure alignments for homologous families.  Protein science, 7(11): 2469–2471.
  * Nodelman et al. (2002) Nodelman, U., Shelton, C. R., and Koller, D. 2002.  Continuous time bayesian networks.  In Proceedings of the Eighteenth conference on Uncertainty in artificial intelligence, pages 378–387. Morgan Kaufmann Publishers Inc.
  * Powell (1994) Powell, M. J. D. 1994.  A Direct Search Optimization Method That Models the Objective and Constraint Functions by Linear Interpolation, pages 51–67.  Springer Netherlands.
  * Price et al. (2010) Price, M. N., Dehal, P. S., and Arkin, A. P. 2010.  FastTree 2–approximately maximum-likelihood trees for large alignments.  PloS one, 5(3): e9490.
  * Rohl et al. (2004) Rohl, C. A., Strauss, C. E., Misura, K. M., and Baker, D. 2004.  Protein structure prediction using Rosetta.  Methods in enzymology, 383: 66–93.
  * Schwartz et al. (2005) Schwartz, A. S., Myers, E. W., and Pachter, L. 2005.  Alignment metric accuracy.  arXiv:q-bio/0510052.
  * Thorne et al. (1992) Thorne, J. L., Kishino, H., and Felsenstein, J. 1992.  Inching toward reality: an improved likelihood model of sequence evolution.  Journal of molecular evolution, 34(1): 3–16.
  * Touw et al. (2015) Touw, W. G., Baakman, C., Black, J., te Beek, T. A., Krieger, E., Joosten, R. P., and Vriend, G. 2015.  A series of PDB-related databanks for everyday needs.  Nucleic acids research, 43(D1): D364–D368.
  * Whelan and Goldman (2001) Whelan, S. and Goldman, N. 2001.  A general empirical model of protein evolution derived from multiple protein families using a maximum-likelihood approach.  Molecular biology and evolution, 18(5): 691–699.

Supplement to “A generative angular model of protein structure evolution”

Michael Golden1,6 Eduardo García-Portugués2,3 Michael Sørensen3 Kanti V.
Mardia4 Thomas Hamelryck2,5 Jotun Hein1

11footnotetext:  Department of Statistics, University of Oxford
(UK).22footnotetext:  Bioinformatics Centre, Section for Computational and RNA
Biology, Department of Biology, University of Copenhagen
(Denmark).33footnotetext:  Department of Mathematical Sciences, University of
Copenhagen (Denmark).44footnotetext:  Department of Mathematics, University of
Leeds (UK).55footnotetext:  Image Section, Department of Computer Science,
University of Copenhagen (Denmark).66footnotetext: Corresponding author.
e-mail<EMAIL_ADDRESS>Figure 11: Diagrammatic representation of
modified TKF92 alignment with neighbour-dependent evolutionary hidden states.
The small white nodes represent non-emitting nodes. The three large oval nodes
represent emitting nodes, emitting a deletion with respect to the first
sequence $\binom{x}{-}$, a insertion with respect to the first sequence
$\binom{-}{y}$ or a homologous pair of amino acids $\binom{x}{y}$. Edges with
filled black circles indicate that the evolutionary hidden state is permitted
to transition to a potentially different evolutionary hidden state, such
transitions are permissible on all edges which exit one of the three emitting
states. Parameters determining the transition probabilities are as follows:
the insertion rate ($\lambda>0$), the deletion rate ($\mu>0$, the extension
rate ($0<r<1$) and the evolutionary time $t_{ab}$. It is required that
$\lambda<\mu$, such that the distribution ($t\rightarrow\infty$) of sequence
lengths is finite at stationarity. The evolutionary hidden state transition
probabilities, $P(H^{i}|H^{i-1})$, are given by a $h\times h$ transition
probability matrix and the initial probabilities at the first alignment site,
$p(H^{i})$, is given by a length $h$ probability vector.

##  Appendix A Supplementary methods

###  A.1 Statistical alignment: modelling insertions and deletions with
neighbouring dependencies

Protein sequences can not only undergo point mutation events, but also also
insertion and deletion (indel) events. We describe a modified pairwise TKF92
alignment HMM that models both local sequence/structure evolution and sequence
alignment based on (Miklós et al., 2004).

Whilst it is possible to fix the alignment in advance by pre-aligning the
sequences using one of the many available optimisation-based alignment methods
(Katoh et al., 2002; Edgar, 2004) or using a curated alignment (such as from
the HOMSTRAD database), doing so ignores alignment uncertainty.

An alignment can be thought of as a statement about homology, such that when
amino acid positions are aligned in order to indicate homology they are
considered to have evolved solely via mutation along the evolutionary
trajectory linking them and therefore not arising via an indel. As the
evolutionary trajectory of indels is rarely observed in practice, it is
difficult to make statements about the true underlying alignment (homology
relationships), especially when the sequences being compared are distantly
related and/or the rate indel of evolution is high. The TKF92 model (Thorne et
al., 1992) gives a suitable distribution describing indel evolution.

For the pairwise case, the TKF92 model can be represented as an HMM using the
formulation described in Miklós et al. (2008). This HMM formulation allows one
to sum over all possible pairwise alignments in $\mathcal{O}(nm)$ time using
the HMM forward-backward algorithm, where $n$ and $m$ are the respective
lengths of the two sequences. Thereby accounting for alignment uncertainty due
to insertions or deletions.

We implemented a modified version of this HMM implementation (Figure 11) such
that each emitted pair of characters is drawn from one of $h$ evolutionary
hidden states (Figure 11). Additionally, we encode neighbouring dependencies
amongst evolutionary hidden states along the alignment, by specifying a
probability transition matrix $p(H^{i},H^{i+1})$, that allows the hidden
states to transition at “Insertion”, “Deletion” or “Match” nodes. These hidden
states are intended to encode local sequence and structure evolution. The
introduction of evolutionary hidden states with neighbouring dependencies
increases the computational complexity from $\mathcal{O}(nm)$ in a model
without evolutionary hidden states to $\mathcal{O}(nmh^{2})$ in a model with
evolutionary hidden states, where $h$ is the number of evolutionary hidden
states. The likelihood of observation pair under this model depends on the
homology relationship, $M_{ab}^{i}$, at a given position is as follows:

$\displaystyle
p(O^{i}|M_{ab}^{i},H^{i},t_{ab})=\begin{cases}p(O_{a}^{x(i)},O_{b}^{y(i)}|H^{i},t_{ab}),&\text{if
observations at positions}\\\ &\text{$x(i)$ and $y(i)$ in proteins $a$ and
$b$,}\\\ &\text{respectively, are homologous;}\\\
p(O_{a}^{x(i)}|H^{i}),&\text{if the observation at}\\\ &\text{position x(i) in
protein $a$}\\\ &\text{is the result of an indel;}\\\
p(O_{b}^{y(i)}|H^{i}),&\text{if the observation at}\\\ &\text{position y(i) in
protein $b$}\\\ &\text{is the result of an indel.}\end{cases}$

$M_{ab}^{i}\in\left\\{\left(\begin{subarray}{c}x\\\
y\end{subarray}\right),\left(\begin{subarray}{c}x\\\
-\end{subarray}\right),\left(\begin{subarray}{c}-\\\
y\end{subarray}\right)\right\\}$ specifies one of three possible homology
relationships at position $i$ in the alignment (homologous amino acids,
deletion with respect to protein $a$, and insertion with respect to protein
$a$, respectively.). $x\in\\{1,\ldots,|a|\\}$ and $y\in\\{1,\ldots,|b|\\}$
specifies the indices of the positions in proteins $a$ and $b$, respectively.

Figure 12: Dihedral angle representation. A very small section of the protein
backbone is displayed in atomic detail. The dihedral angles, $\phi$ and
$\psi$, are shown. The partial double bond nature of the peptide bond, gives
rise to a planar configuration of O, C$\beta$, N and H atoms. This allows the
protein backbone structure to be largely described in terms of a series of
$\phi$ and $\psi$ dihedral angles that relate the planes in three-dimensional
space.

Note that a disadvantage of our approach is that neighbouring amino acid
positions in the presence of a deletion or series of deletions are no longer
treated as directly adjacent by our HMM with respect to that particular
protein, whereas in physical reality the amino acids positions would be
directly adjacent to one another. However, the original HMM formulation where
we treat the alignment as given a priori has the same shortcoming.

###  A.2 Parallelisation of model training

The StEM algorithm is easily parallelised. In the E-step, the parameters,
alignment and hidden states of protein pair proteins can be independently
sampled in parallel, conditioned on the parameters, $\Psi^{(r)}$, from the
M-step. Whereas, in the M-step the parameters corresponding to each hidden
state can be independently updated in parallel conditioned on the samples from
the E-step.

###  A.3 Time-reversibility

The three stochastic processes are assumed to be time-reversible. This,
together with assumption of time-reversibility in jumping between evolutionary
regimes in equation ensures overall time-reversibility. This allows to treat
the phylogenetic tree relating proteins $p_{a}$ and $p_{b}$ as unrooted,
implying we can arbitrarily pick $p_{a}$ or $p_{b}$ as a root of the
phylogenetic tree (Felsenstein, 1981). This avoids the need to marginalise
over the common ancestor protein of proteins $p_{a}$ and $p_{b}$. Note that
whilst time-reversibility of the evolutionary processes at each site holds,
this is different from reversibility of the HMM. The transition probability
matrix of the HMM is not restricted to be reversible and therefore detailed-
balance does not necessarily hold.

Figure 13: Time-reversibility allows an arbitrary rooting of the phylogeny.
The stochastic processes describing the observations have all been chosen such
that time-reversibility holds. This is particularly useful in the case of a
pairwise phylogeny, because it permits us to arbitrarily pick one of the
extant proteins (protein $a$ or protein $b$) as the root of the phylogeny,
without changing the likelihood of the data. This avoids computationally
costly marginalisation of the unobserved common ancestor protein that in
reality is shared by both of the extant proteins. Furthermore, only a single
parameter, $t_{ab}$, the evolutionary time, need be marginalised when training
the model or performing inference.

##  Appendix B Supplementary results

###  B.1 Estimates of evolutionary time from dihedral angles are consistent
with estimates from sequence

Figure 14 compares evolutionary times estimated using pairs of homologous
amino acid sequences only versus pairs of homologous dihedral angles only. The
mean 90% confidence interval (CI) widths for the sampled evolutionary times
were smaller for amino acid sequences (mean 90% CI width $=0.139$) than
dihedral angles (mean 90% CI width $=0.191$).

## References

  * Edgar (2004) Edgar, R. C. 2004.  MUSCLE: multiple sequence alignment with high accuracy and high throughput.  Nucleic acids research, 32(5): 1792–1797.
  * Felsenstein (1981) Felsenstein, J. 1981.  Evolutionary trees from DNA sequences: a maximum likelihood approach.  Journal of molecular evolution, 17(6): 368–376.
  * Katoh et al. (2002) Katoh, K., Misawa, K., Kuma, K.-i., and Miyata, T. 2002.  Mafft: a novel method for rapid multiple sequence alignment based on fast fourier transform.  Nucleic acids research, 30(14): 3059–3066.
  * Miklós et al. (2004) Miklós, I., Lunter, G., and Holmes, I. 2004.  A “long indel” model for evolutionary sequence alignment.  Molecular Biology and Evolution, 21(3): 529–540.
  * Thorne et al. (1992) Thorne, J. L., Kishino, H., and Felsenstein, J. 1992.  Inching toward reality: an improved likelihood model of sequence evolution.  Journal of molecular evolution, 34(1): 3–16.

Figure 14: Evolutionary times were estimated for 39 protein pairs under a
model with jump evolution under two different conditions. Under the first
condition where only the amino acid sequences were treat as observed, whereas
under the second only the dihedral angles were treated as observed.