BACKGROUND
conformational heterogeneity is a common theme in protein structures and relevant in a wide range of different contexts.
proteins are flexible macromolecules whose function is often accompanied by a structural transition  <cit> .
allostery, for example, is a ubiquitous mechanism in signal transduction  <cit>  and continues to be a controversial field of research  <cit> .
structural heterogeneity may also stem from a lack of data.
nmr structures are usually represented as ensembles of conformers that fit the data equally well  <cit> .
here structural heterogeneity mainly reflects a scarcity of restraints and not necessarily true conformational flexibility.
when comparing protein structures in different conformational states, one is mainly interested in internal structural changes rather than differences that can be accounted for by a rigid-body movement.
the separation of external from internal movements directly relates to the problem of how to compare and superimpose protein structures in different conformational states.
the hallmark of protein structure comparison is the root mean square deviation  between equivalent atom positions after the rigid modes of structural change have been removed.
the rmsd defines an optimality criterion to determine the rotation and translation that best separate rigid-body from internal movements.
how to minimize the rmsd over all possible translations and rotations is a classical problem in structural bioinformatics and has been treated by many authors  <cit> .
after early accounts by diamond and mclachlan  <cit> , kabsch has given a closed analytical expression for the optimal translation and rotation in terms of a singular value decomposition  <cit> .
kearsley and others have provided an alternative solution based on quaternions that improves the kabsch algorithm in terms of speed and stability  <cit> .
a physical justification for superimposing protein structures by rmsd minimization originates in the theory of dynamics in semi-rigid molecules.
eckart has derived  <cit>  conditions for the separation of external  from internal modes of movement, if the molecule is subject to small-amplitude vibrational motions.
recently, it has been pointed out that structure comparison by rmsd minimization is equivalent to searching for the frame of reference that satisfies eckart's conditions  <cit> .
therefore, if one considers a set of heterogenous structures an ensemble of fluctuating states, rmsd minimization is the physically correct method for removing rigid-body displacements.
however, in many situations one is interested in a frame of reference that is different from the eckart frame.
such a situation may occur if we want to compare proteins that undergo structural transitions upon interaction with other molecules.
this is a non-equilibrium situation in which the protein is driven to a different energy basin.
a classical example is adenylate kinase comprising three rigid domains that undergo an opening-closing conformational transition upon the binding of substrate  <cit> .
here, rmsd fitting fails to highlight the relative rigid-body movements leading to domain closure.
because the assumption of vibrational conformational changes is not fulfilled over the entire polypeptide chain, atoms that belong to the mobile domains would appear as "outliers" that cannot be described by vibrational dynamics.
rmsd minimization is a least-squares technique and therefore suffers from the same problems that least-squares methods have in other data analysis applications, namely sensitivity to outliers.
problems with outliers in protein structure comparison have been treated in a number of ways.
one simple fix is to extend the kabsch formula to weighted rmsd superposition.
a weight is assigned to every atom and applied when summing over the distances between equivalent positions.
by choosing small weights for outliers their dominance can be alleviated.
other approaches build on different metrics than euclidean distances.
lesk  <cit> , for example, uses the chebyshev distance to identify common substructures in proteins.
the lms fit algorithm  <cit>  minimizes the median rather the mean squared deviation between equivalent atoms.
another class of methods seeks to find a flexible, rather than a rigid alignment between protein structures  <cit> .
the possibility to weight atoms individually is already mentioned in kabsch's paper but no rule for setting the weights is given.
as a consequence, different schemes for weighting atoms have been proposed.
an important class of algorithms iteratively filters out atoms whose deviation exceeds a predefined threshold  <cit> .
this strategy tries to identify those atoms that make up an invariant structural core.
more recent applications assign continuous weights to every atom, such that all atoms contribute to the weighted rmsd.
the gaussian-weighted rmsd  method  <cit>  updates the weights iteratively by plugging distances between equivalent positions into a gaussian distribution.
wu and wu  <cit>  set the weight to a theoretical temperature factor predicted with a gaussian network model.
here we elaborate on the idea of using a weighted rmsd to find a superposition that identifies geometrically similar substructures in heterogeneous protein structures.
the method is intended for the comparison of structures that exhibit significant structural disparity or undergo large-scale conformational change.
we introduce probabilistic models that describe very generic properties of the large-amplitude structural changes that we want to account for.
we show that learning these models solves the superposition problem also in the presence of gross structural transitions and is equivalent to a weighted rmsd whose weights are updated iteratively.
a problem with existing methods minimizing a weighted rmsd is that the choice of the weights is heuristic and often depends on user-defined parameters such as thresholds  <cit>  and decay constants  <cit> .
our approach is principled and provides objective rules for setting the weights.
by averaging analytically over the weights we show that our models for structural displacements are heavy-tailed distributions, which are often used to describe extreme events, for example, in economics.
model comparison techniques allow us to choose among the different models and objectively select the one that is most supported by the structures.
we illustrate the applicability of the approach in various contexts ranging from structural changes in proteins to nmr ensembles.
from a practical point of view, the new method provides an objective and robust basis for protein structure comparison and superposition.
RESULTS
frames of reference in protein structure comparison
rather than adopting a coordinate system prescribed by physical principles such as the eckart frame, we are interested in finding a frame of reference that overlays heterogenous protein structures such that they share a tight, maximally large structural core at the cost of large-scale outliers.
to illustrate the difference consider a normal-mode simulation of elongation factor g , which contains five compact domains.
we used the elastic network model implemented in mmtk  <cit>  to sample conformations about 1fnm.
the elastic network model imposes purely vibrational interactions between restrained atoms.
snapshots along the first principal mode are, by construction, related only via internal displacements.
this is indeed verified by calculating the rmsd between, for example, the first and last conformation: the optimal rotation is the identity matrix, the translation vector zero.
in the weighted rmsd frame that maximizes the overlay of spatially invariant positions, the differences between the two structures are interpreted as a finite rigid-body movement .
to describe this situation quantitatively, we assume that alternative positions of the ith atom are encoded by three-dimensional vectors xi and yi.
the most general relation between equivalent positions is given by the generative model:
where the rigid-body transformation involves a rotation matrix r and translation vector t and the vectors di are the non-rigid displacements.
in the eckart frame, the structures are displaced around the center of mass, whereas in a superposition maximizing the common structural core only domain iv seems to move .
consequently, we observe different distributions for the displacements depending on whether we superimpose by unweighted or core-weighted rmsd.
in the weighted fit, the distribution of displacements exhibits a narrow central peak corresponding to the well-fitting core and broad tails accounting for large-amplitude movements .
in the eckart frame, the displacements are distributed more homogeneously.
we observe similar large-amplitude displacements in a comparison of groel in atp-free and atp-bound state  demonstrating that the shape of the distribution is universal irrespective of the specific context.
if we knew the displacements exactly, we could obtain the rotation and translation by solving the above system of equations.
in general, however, the internal displacements are unknown and need to be estimated simultaneously with the superposition.
this is a chicken-and-egg problem: estimation of the rigid body transformation requires knowledge of the internal changes, which themselves can be calculated only if the superposition is known.
in general, the decomposition into rigid-body and internal changes is not unique.
we need additional principles to estimate the decomposition from a set of heterogenous structures.
modeling non-rigid structural change in proteins
to separate internal from external structural changes we need to make assumptions about their properties.
our assumptions will be of statistical nature, because we do not know the correct displacements and can only infer them from the given structures.
we encode our assumptions in a probability distribution over the displacement vectors f. given this distribution, a statistical approach to protein structure comparison proceeds by plugging the displacements calculated according to generative model into the distribution.
this results in the total likelihood for the separation into rigid and non-rigid structural changes under the assumed model f:
and depends on the choice of the rigid-body transformation.
the optimal separation is obtained by maximizing l over all rotations and translations to implicitly obtain displacements whose distribution best matches f.
what are reasonable assumptions about non-rigid structural changes in proteins that are realistic and, at the same time, simple enough to allow for efficient computation?
one straightforward property of f is that its mean is zero, because if it were not, the mean could not be distinguished from the overall translation.
second, it should not matter which of the two structures we superimpose onto the other.
the displacements resulting from either superposition should follow the same distribution.
the reverse of the generative model is xi = rtyi - t - rtdi.
that is,  = -rtdi are the displacements according to the reverse superposition expressed through the parameters of the original superposition.
we demand f() = f from which follows that f is isotropic, i.e.
it depends only on the norm of the displacements, not their direction.
an intuitive quantity to characterize the displacements are their expected amplitudes ai which are the second moments f. because we assume f to be isotropic, we only need to consider the average squared norm of the displacements  rather than the full covariance matrix.
for a gaussian model, the amplitude is directly related to the isotropic variance σ2: ai = 3σ <dig>  the assumption underlying rmsd fitting is that the unknown displacements occur in a homogenous fashion: all amplitudes vary on the same scale σ, which is a statistical analog of the eckart assumption.
however, as discussed before, we are interested in modeling displacements that are most of the time small but occasionally huge.
a natural extension is to associate a separate variance  with each displacement.
for mathematical convenience, we work with precisions si =  rather than variances; si can be viewed as a measure of local stiffiness.
by allowing the si to adopt values of diverse orders of magnitude, we can model displacements that vary on different scales.
according to this model, highly restrained atoms will be assigned large si as is the case, for example, for core atoms that are restricted in their mobility by nearest neighbor interactions.
on the opposite end of the spectrum, large-scale displacements as, for example, in t <dig> rna polymerase  <cit>  would be described with si's that are close to zero.
under a multi-scale model of structural change, the optimal rotation and translation is obtained by maximizing the corresponding total likelihood or rather minimizing its negative logarithm:
where the inverse amplitudes si =  are atom specific weights.
amplitude spectra of large-scale conformational changes
if the amplitudes of the displacements were known, the optimal rigid body transformation could be determined by a weighted superposition.
however, the scales si are unknown, and it seems that we just shifted the problem.
we need to estimate the scales and therefore must make assumptions about their distribution g. because the scales are non-zero, their distribution should only have support on the positive axis.
if we interpret the si as force constants, the spectrum should span a wide range of variability owing to the fact that protein structures are subject to internal forces of different strengths.
we use a gamma and an inverse gamma distribution for g. these distributions are flexible enough to gradually switch between different situations in which we have more or less homogeneously distributed amplitudes.
the functional form of gamma and inverse gamma distributions is governed by two parameters, a shape parameter α and a scale β, which is an overall scale of the displacement amplitudes.
what is the distribution of the displacements f implied by our choice of g?
to answer this question we need to average over all scales.
figure  <dig> illustrates this averaging process.
the effective distribution of non-rigid displacements exhibits a dominant central peak but also allows outliers to occur occasionally, which is reflected in the elevated "fatness" of the tails when compared with a gaussian distribution.
in the statistical literature, distributions of this form are called heavy-tailed.
heavy-tailed distributions are robust against outliers in data  <cit>  and can often be represented as averages over gaussian distributions with zero mean and increasing width, so-called scale mixtures of gaussian distributions  <cit> .
from a pragmatic point of view, the average distribution of internal displacements has the convenient property that it is heavy-tailed and thereby accommodates large-amplitude structural changes.
for a gamma and inverse gamma-distribution, the effective distribution f can be calculated analytically.
the former corresponds to a student t distribution, the latter is a member of the k distribution family  <cit>  comprising the laplace distribution as a special case.
in their superposition algorithm theseus  <cit> , theobald and wuttke use a related model.
they introduce a multivariate gaussian with a full-covariance structure for inter-positional dependencies and estimate the covariance matrix during the superposition.
the eigenvalues of the covariance matrix are assumed to be distributed according to an inverse gamma distribution.
in case of a diagonal covariance matrix, the model is equivalent to a gamma prior on the scales  and thus to the student t model.
algorithms for disentangling rigid and non-rigid structural changes
the generalized optimization problem that we face when maximizing l is to minimize:
that is, f implies the generalized metric - log f for comparing equivalent atom positions.
for gaussian f, this metric is the squared euclidean distance.
for laplace-distributed displacements, it is the euclidean distance .
analytical expressions for the best rigid transformation under such metrics are not known.
one could optimize the negative log-likelihood numerically.
however, we will pursue an alternative approach based on the scale mixture representations that we discussed in the previous section.
the basic idea is illustrated in a flowchart .
instead of averaging analytically over the scales, we estimate them simultaneously with the rigid transformation and the parameters α and β that determine the shape of g. the algorithm proceeds iteratively by updating the different groups of parameters, {si},  and , separately.
updates of the rigid transformation involve a weighted rmsd fit.
how to update the weights si, the shape α and the scale β is detailed in the methods section and additional file  <dig>
we have developed two version of the updates: a deterministic and a stochastic algorithm.
the latter, a gibbs sampler, allows us to estimate not only the optimal parameters, but also their uncertainty and enables comparison of different models, as will be shown later.
both algorithms converge very rapidly within the first  <dig> iterations, after which the likelihood does no longer improve significantly.
on a modern computer the best em fit is obtained within a few fractions of a second.
the gibbs sampler is approximately ten times slower.
running times for both algorithms on different protein structure pairs are provided in additional file  <dig>
robust fitting of proteins subject to large conformational changes
to demonstrate the validity of our approach, we tested the framework on a number of proteins undergoing large conformational changes.
figure  <dig> shows the superposition of two different conformations of groel  <cit>  based on a gaussian, student t, and k distribution.
the transition of groel from an unbound to a bound state involves rigid body movements of the intermediate and apical domains.
the superposition based on a gaussian model fails to reveal the relative movements of the domains, whereas the heavy-tailed models converge on a tight fit of the equatorial domain.
this observation is supported by an analysis of the individual domains.
although the overall rmsd increases from  <dig>  Å for the least-squares superposition to  <dig>  Å for both non-gaussian methods, the rmsd of the equatorial domain drops from  <dig>  Å down to  <dig>  Å and  <dig>  Å for the student t and k distribution, respectively.
both values are close to the optimal rmsd of  <dig>  Å when fitting the equatorial domain alone.
the reliability of the superposition can be assessed through the structure ensembles obtained by applying the random transformations generated during gibbs sampling.
the ensemble generated by gaussian superposition is broad reflecting a high degree of uncertainty.
in contrast, the ensembles based on heavy-tailed distributions are narrow, which indicates that the superposition is very well defined .
the findings are confirmed by looking at the histograms of non-rigid displacements and their parametric fits .
only the heavy-tailed distributions fit the displacments reasonably well, whereas a gaussian fails to describe the simultaneous occurance of many well-fitting positions and a few large-scale outliers.
these fits are obtained automatically during the superposition - once a model has been chosen, all unknown parameters are estimated self-consistently .
in contrast, other superposition algorithms based on a weighted rmsd involve adjustable parameters that are set heuristically.
this may lead to problems with identifying the optimal structural core .
the domain architecture of groel is highlighted by a trace plot of the local scales .
large weights are assigned to the well fitting equatorial domain and almost all weights for the intermediate and apical domain are close to zero.
heavy-tailed models are well-suited to describe even large conformational changes.
this is exemplified for pneumolysin  <cit>  which, upon membrane insertion, refolds two of its four domains leading to an invariant core of about 30% of all residues only.
figure  <dig> compares the least-squares superposition with superpositions based on a student t and k distribution.
again only a non-gaussian superposition is able to locate and fit the invariant region.
superposition of nmr ensembles
nmr structures are usually represented as ensembles that reflect the quality and completeness of the data as well as the local precision of the structure  <cit> .
often termini and loops show high variability either due to protein dynamics or missing data.
if one wants to assess the precision of an nmr structure, superposition by rmsd minimization often fails to reflect local differences due to variations in restraint density  <cit> .
as a consequence, no generally accepted way to fit ensembles exist.
the superposition is often determined based on secondary structure elements or subjective criteria such as a small number of manually defined positions.
our framework provides a more objective, robust and model-driven alternative to such practice.
a particularly suited example to demonstrate the violation of the least-squares assumption is the nmr ensemble of calmodulin .
calmodulin is involved in cellular regulation and consists of two identical domains connected by a flexible hinge region  <cit> .
the flexibility is also reflected in the nmr ensemble.
if one domain is superimposed, the other undergoes a large relative motion, which makes a joint superposition impossible.
figure  <dig> highlights differences between a gaussian and an outlier-tolerant superposition for this example.
the heavy-tailed models favor a superposition onto the n-terminal domain with an internal rmsd of  <dig>  Å, whereas the least-squares algorithm fails to find a tight superposition.
by tweaking the initial weights, it is possible to find the superposition onto the c-terminal domain with an rmsd of  <dig>  Å, albeit this superposition has a lower likelihood.
the higher likelihood of a superposition onto the n-terminal domain is consistent with the observation by kuboniwa et al.
<cit>  that the n-terminal is better defined than c-terminal domain and shows a lower internal rmsd.
bayesian model comparison
so far, we a given model of the displacements to a set of structures.
bayesian inference allows us to go further and infer which model is the most appropriate given a set of structures.
this information is provided by the evidence or marginal likelihood p, the probability of model m  given data d   <cit> .
the evidence is computed by integration over all possible parameter values.
two models m <dig> and m <dig> are ranked relative to each other through the bayes factor p/p   <cit> .
if the bayes factor is significantly greater than one, the data favor model m <dig> over model m <dig> and vice versa.
here we seek to assess whether to choose the student t or k distribution over a gaussian model for structure superposition and comparison.
because the evidence is not amenable to analytical evaluation, we use an estimator calculated from the posterior samples obtained with gibbs sampling  <cit> .
logarithm of the marginal likelihood p of the different displacement models obtained for nine structure pairs undergoing domain movements.
highlighted in boldface are the maximum log-marginal likelihood values.
the last two rows at the bottom report the log-marginal likelihoods for synthetic data generated according to a student t and a gauss distribution.
CONCLUSIONS
we present a robust probabilistic approach to protein structure superposition and comparison.
the approach builds on heavy-tailed distributions to model non-rigid displacements between protein structures.
to estimate these distributions and an optimal superposition we employ a scale mixture representation of the heavy-tailed models.
practically, this amounts to introducing weights for each atom position and to estimate the weights iteratively during structure superposition.
in contrast to other weight-based superposition methods, the scale mixture framework provides a firm statistical basis for setting the weights.
moreover, the link to the closed form helps to interpret the weighting scheme in terms of heavy-tailed models for structural displacements.

 conformational heterogeneity is a common theme in protein structures and relevant in a wide range of different contexts. proteins are flexible macromolecules whose function is often accompanied by a structural transition. allostery, for example, is a ubiquitous mechanism in signal transduction and continues to be a controversial field of research. structural heterogeneity may also stem from a lack of data. nmr structures are usually represented as ensembles of conformers that fit the data equally well. here structural heterogeneity mainly reflects a scarcity of restraints and not necessarily true conformational flexibility. when comparing protein structures in different conformational states, one is mainly interested in internal structural changes rather than differences that can be accounted for by a rigid-body movement. the separation of external from internal movements directly relates to the problem of how to compare and superimpose protein structures in different conformational states. the hallmark of protein structure comparison is the root mean square deviation between equivalent atom positions after the rigid modes of structural change have been removed. the rmsd defines an optimality criterion to determine the rotation and translation that best separate rigid-body from internal movements. how to minimize the rmsd over all possible translations and rotations is a classical problem in structural bioinformatics and has been treated by many authors. after early accounts by diamond and mclachlan, kabsch has given a closed analytical expression for the optimal translation and rotation in terms of a singular value decomposition. kearsley and others have provided an alternative solution based on quaternions that improves the kabsch algorithm in terms of speed and stability. a physical justification for superimposing protein structures by rmsd minimization originates in the theory of dynamics in semi-rigid molecules. eckart has derived conditions for the separation of external from internal modes of movement, if the molecule is subject to small-amplitude vibrational motions. recently, it has been pointed out that structure comparison by rmsd minimization is equivalent to searching for the frame of reference that satisfies eckart's conditions. therefore, if one considers a set of heterogenous structures an ensemble of fluctuating states, rmsd minimization is the physically correct method for removing rigid-body displacements. however, in many situations one is interested in a frame of reference that is different from the eckart frame. such a situation may occur if we want to compare proteins that undergo structural transitions upon interaction with other molecules. this is a non-equilibrium situation in which the protein is driven to a different energy basin. a classical example is adenylate kinase comprising three rigid domains that undergo an opening-closing conformational transition upon the binding of substrate. here, rmsd fitting fails to highlight the relative rigid-body movements leading to domain closure. because the assumption of vibrational conformational changes is not fulfilled over the entire polypeptide chain, atoms that belong to the mobile domains would appear as "outliers" that cannot be described by vibrational dynamics. rmsd minimization is a least-squares technique and therefore suffers from the same problems that least-squares methods have in other data analysis applications, namely sensitivity to outliers. problems with outliers in protein structure comparison have been treated in a number of ways. one simple fix is to extend the kabsch formula to weighted rmsd superposition. a weight is assigned to every atom and applied when summing over the distances between equivalent positions. by choosing small weights for outliers their dominance can be alleviated. other approaches build on different metrics than euclidean distances. lesk, for example, uses the chebyshev distance to identify common substructures in proteins. the lms fit algorithm minimizes the median rather the mean squared deviation between equivalent atoms. another class of methods seeks to find a flexible, rather than a rigid alignment between protein structures. the possibility to weight atoms individually is already mentioned in kabsch's paper but no rule for setting the weights is given. as a consequence, different schemes for weighting atoms have been proposed. an important class of algorithms iteratively filters out atoms whose deviation exceeds a predefined threshold. this strategy tries to identify those atoms that make up an invariant structural core. more recent applications assign continuous weights to every atom, such that all atoms contribute to the weighted rmsd. the gaussian-weighted rmsd method updates the weights iteratively by plugging distances between equivalent positions into a gaussian distribution. wu and wu set the weight to a theoretical temperature factor predicted with a gaussian network model. here we elaborate on the idea of using a weighted rmsd to find a superposition that identifies geometrically similar substructures in heterogeneous protein structures. the method is intended for the comparison of structures that exhibit significant structural disparity or undergo large-scale conformational change. we introduce probabilistic models that describe very generic properties of the large-amplitude structural changes that we want to account for. we show that learning these models solves the superposition problem also in the presence of gross structural transitions and is equivalent to a weighted rmsd whose weights are updated iteratively. a problem with existing methods minimizing a weighted rmsd is that the choice of the weights is heuristic and often depends on user-defined parameters such as thresholds and decay constants. our approach is principled and provides objective rules for setting the weights. by averaging analytically over the weights we show that our models for structural displacements are heavy-tailed distributions, which are often used to describe extreme events, for example, in economics. model comparison techniques allow us to choose among the different models and objectively select the one that is most supported by the structures. we illustrate the applicability of the approach in various contexts ranging from structural changes in proteins to nmr ensembles. from a practical point of view, the new method provides an objective and robust basis for protein structure comparison and superposition.  frames of reference in protein structure comparison rather than adopting a coordinate system prescribed by physical principles such as the eckart frame, we are interested in finding a frame of reference that overlays heterogenous protein structures such that they share a tight, maximally large structural core at the cost of large-scale outliers. to illustrate the difference consider a normal-mode simulation of elongation factor g, which contains five compact domains. we used the elastic network model implemented in mmtk to sample conformations about 1fnm. the elastic network model imposes purely vibrational interactions between restrained atoms. snapshots along the first principal mode are, by construction, related only via internal displacements. this is indeed verified by calculating the rmsd between, for example, the first and last conformation: the optimal rotation is the identity matrix, the translation vector zero. in the weighted rmsd frame that maximizes the overlay of spatially invariant positions, the differences between the two structures are interpreted as a finite rigid-body movement. to describe this situation quantitatively, we assume that alternative positions of the ith atom are encoded by three-dimensional vectors xi and yi. the most general relation between equivalent positions is given by the generative model: where the rigid-body transformation involves a rotation matrix r and translation vector t and the vectors di are the non-rigid displacements. in the eckart frame, the structures are displaced around the center of mass, whereas in a superposition maximizing the common structural core only domain iv seems to move. consequently, we observe different distributions for the displacements depending on whether we superimpose by unweighted or core-weighted rmsd. in the weighted fit, the distribution of displacements exhibits a narrow central peak corresponding to the well-fitting core and broad tails accounting for large-amplitude movements. in the eckart frame, the displacements are distributed more homogeneously. we observe similar large-amplitude displacements in a comparison of groel in atp-free and atp-bound state demonstrating that the shape of the distribution is universal irrespective of the specific context. if we knew the displacements exactly, we could obtain the rotation and translation by solving the above system of equations. in general, however, the internal displacements are unknown and need to be estimated simultaneously with the superposition. this is a chicken-and-egg problem: estimation of the rigid body transformation requires knowledge of the internal changes, which themselves can be calculated only if the superposition is known. in general, the decomposition into rigid-body and internal changes is not unique. we need additional principles to estimate the decomposition from a set of heterogenous structures. modeling non-rigid structural change in proteins to separate internal from external structural changes we need to make assumptions about their properties. our assumptions will be of statistical nature, because we do not know the correct displacements and can only infer them from the given structures. we encode our assumptions in a probability distribution over the displacement vectors f. given this distribution, a statistical approach to protein structure comparison proceeds by plugging the displacements calculated according to generative model into the distribution. this results in the total likelihood for the separation into rigid and non-rigid structural changes under the assumed model f: and depends on the choice of the rigid-body transformation. the optimal separation is obtained by maximizing l over all rotations and translations to implicitly obtain displacements whose distribution best matches f. what are reasonable assumptions about non-rigid structural changes in proteins that are realistic and, at the same time, simple enough to allow for efficient computation? one straightforward property of f is that its mean is zero, because if it were not, the mean could not be distinguished from the overall translation. second, it should not matter which of the two structures we superimpose onto the other. the displacements resulting from either superposition should follow the same distribution. the reverse of the generative model is xi = rtyi - t - rtdi. that is, = -rtdi are the displacements according to the reverse superposition expressed through the parameters of the original superposition. we demand f() = f from which follows that f is isotropic, i.e. it depends only on the norm of the displacements, not their direction. an intuitive quantity to characterize the displacements are their expected amplitudes ai which are the second moments f. because we assume f to be isotropic, we only need to consider the average squared norm of the displacements rather than the full covariance matrix. for a gaussian model, the amplitude is directly related to the isotropic variance σ2: ai = 3σ the assumption underlying rmsd fitting is that the unknown displacements occur in a homogenous fashion: all amplitudes vary on the same scale σ, which is a statistical analog of the eckart assumption. however, as discussed before, we are interested in modeling displacements that are most of the time small but occasionally huge. a natural extension is to associate a separate variance with each displacement. for mathematical convenience, we work with precisions si = rather than variances; si can be viewed as a measure of local stiffiness. by allowing the si to adopt values of diverse orders of magnitude, we can model displacements that vary on different scales. according to this model, highly restrained atoms will be assigned large si as is the case, for example, for core atoms that are restricted in their mobility by nearest neighbor interactions. on the opposite end of the spectrum, large-scale displacements as, for example, in t rna polymerase would be described with si's that are close to zero. under a multi-scale model of structural change, the optimal rotation and translation is obtained by maximizing the corresponding total likelihood or rather minimizing its negative logarithm: where the inverse amplitudes si = are atom specific weights. amplitude spectra of large-scale conformational changes if the amplitudes of the displacements were known, the optimal rigid body transformation could be determined by a weighted superposition. however, the scales si are unknown, and it seems that we just shifted the problem. we need to estimate the scales and therefore must make assumptions about their distribution g. because the scales are non-zero, their distribution should only have support on the positive axis. if we interpret the si as force constants, the spectrum should span a wide range of variability owing to the fact that protein structures are subject to internal forces of different strengths. we use a gamma and an inverse gamma distribution for g. these distributions are flexible enough to gradually switch between different situations in which we have more or less homogeneously distributed amplitudes. the functional form of gamma and inverse gamma distributions is governed by two parameters, a shape parameter α and a scale β, which is an overall scale of the displacement amplitudes. what is the distribution of the displacements f implied by our choice of g? to answer this question we need to average over all scales. figure illustrates this averaging process. the effective distribution of non-rigid displacements exhibits a dominant central peak but also allows outliers to occur occasionally, which is reflected in the elevated "fatness" of the tails when compared with a gaussian distribution. in the statistical literature, distributions of this form are called heavy-tailed. heavy-tailed distributions are robust against outliers in data and can often be represented as averages over gaussian distributions with zero mean and increasing width, so-called scale mixtures of gaussian distributions. from a pragmatic point of view, the average distribution of internal displacements has the convenient property that it is heavy-tailed and thereby accommodates large-amplitude structural changes. for a gamma and inverse gamma-distribution, the effective distribution f can be calculated analytically. the former corresponds to a student t distribution, the latter is a member of the k distribution family comprising the laplace distribution as a special case. in their superposition algorithm theseus, theobald and wuttke use a related model. they introduce a multivariate gaussian with a full-covariance structure for inter-positional dependencies and estimate the covariance matrix during the superposition. the eigenvalues of the covariance matrix are assumed to be distributed according to an inverse gamma distribution. in case of a diagonal covariance matrix, the model is equivalent to a gamma prior on the scales and thus to the student t model. algorithms for disentangling rigid and non-rigid structural changes the generalized optimization problem that we face when maximizing l is to minimize: that is, f implies the generalized metric - log f for comparing equivalent atom positions. for gaussian f, this metric is the squared euclidean distance. for laplace-distributed displacements, it is the euclidean distance. analytical expressions for the best rigid transformation under such metrics are not known. one could optimize the negative log-likelihood numerically. however, we will pursue an alternative approach based on the scale mixture representations that we discussed in the previous section. the basic idea is illustrated in a flowchart. instead of averaging analytically over the scales, we estimate them simultaneously with the rigid transformation and the parameters α and β that determine the shape of g. the algorithm proceeds iteratively by updating the different groups of parameters, {si}, and, separately. updates of the rigid transformation involve a weighted rmsd fit. how to update the weights si, the shape α and the scale β is detailed in the methods section and additional file we have developed two version of the updates: a deterministic and a stochastic algorithm. the latter, a gibbs sampler, allows us to estimate not only the optimal parameters, but also their uncertainty and enables comparison of different models, as will be shown later. both algorithms converge very rapidly within the first iterations, after which the likelihood does no longer improve significantly. on a modern computer the best em fit is obtained within a few fractions of a second. the gibbs sampler is approximately ten times slower. running times for both algorithms on different protein structure pairs are provided in additional file robust fitting of proteins subject to large conformational changes to demonstrate the validity of our approach, we tested the framework on a number of proteins undergoing large conformational changes. figure shows the superposition of two different conformations of groel based on a gaussian, student t, and k distribution. the transition of groel from an unbound to a bound state involves rigid body movements of the intermediate and apical domains. the superposition based on a gaussian model fails to reveal the relative movements of the domains, whereas the heavy-tailed models converge on a tight fit of the equatorial domain. this observation is supported by an analysis of the individual domains. although the overall rmsd increases from Å for the least-squares superposition to Å for both non-gaussian methods, the rmsd of the equatorial domain drops from Å down to Å and Å for the student t and k distribution, respectively. both values are close to the optimal rmsd of Å when fitting the equatorial domain alone. the reliability of the superposition can be assessed through the structure ensembles obtained by applying the random transformations generated during gibbs sampling. the ensemble generated by gaussian superposition is broad reflecting a high degree of uncertainty. in contrast, the ensembles based on heavy-tailed distributions are narrow, which indicates that the superposition is very well defined. the findings are confirmed by looking at the histograms of non-rigid displacements and their parametric fits. only the heavy-tailed distributions fit the displacments reasonably well, whereas a gaussian fails to describe the simultaneous occurance of many well-fitting positions and a few large-scale outliers. these fits are obtained automatically during the superposition - once a model has been chosen, all unknown parameters are estimated self-consistently. in contrast, other superposition algorithms based on a weighted rmsd involve adjustable parameters that are set heuristically. this may lead to problems with identifying the optimal structural core. the domain architecture of groel is highlighted by a trace plot of the local scales. large weights are assigned to the well fitting equatorial domain and almost all weights for the intermediate and apical domain are close to zero. heavy-tailed models are well-suited to describe even large conformational changes. this is exemplified for pneumolysin which, upon membrane insertion, refolds two of its four domains leading to an invariant core of about 30% of all residues only. figure compares the least-squares superposition with superpositions based on a student t and k distribution. again only a non-gaussian superposition is able to locate and fit the invariant region. superposition of nmr ensembles nmr structures are usually represented as ensembles that reflect the quality and completeness of the data as well as the local precision of the structure. often termini and loops show high variability either due to protein dynamics or missing data. if one wants to assess the precision of an nmr structure, superposition by rmsd minimization often fails to reflect local differences due to variations in restraint density. as a consequence, no generally accepted way to fit ensembles exist. the superposition is often determined based on secondary structure elements or subjective criteria such as a small number of manually defined positions. our framework provides a more objective, robust and model-driven alternative to such practice. a particularly suited example to demonstrate the violation of the least-squares assumption is the nmr ensemble of calmodulin. calmodulin is involved in cellular regulation and consists of two identical domains connected by a flexible hinge region. the flexibility is also reflected in the nmr ensemble. if one domain is superimposed, the other undergoes a large relative motion, which makes a joint superposition impossible. figure highlights differences between a gaussian and an outlier-tolerant superposition for this example. the heavy-tailed models favor a superposition onto the n-terminal domain with an internal rmsd of Å, whereas the least-squares algorithm fails to find a tight superposition. by tweaking the initial weights, it is possible to find the superposition onto the c-terminal domain with an rmsd of Å, albeit this superposition has a lower likelihood. the higher likelihood of a superposition onto the n-terminal domain is consistent with the observation by kuboniwa et al. that the n-terminal is better defined than c-terminal domain and shows a lower internal rmsd. bayesian model comparison so far, we a given model of the displacements to a set of structures. bayesian inference allows us to go further and infer which model is the most appropriate given a set of structures. this information is provided by the evidence or marginal likelihood p, the probability of model m given data d. the evidence is computed by integration over all possible parameter values. two models m and m are ranked relative to each other through the bayes factor p/p. if the bayes factor is significantly greater than one, the data favor model m over model m and vice versa. here we seek to assess whether to choose the student t or k distribution over a gaussian model for structure superposition and comparison. because the evidence is not amenable to analytical evaluation, we use an estimator calculated from the posterior samples obtained with gibbs sampling. logarithm of the marginal likelihood p of the different displacement models obtained for nine structure pairs undergoing domain movements. highlighted in boldface are the maximum log-marginal likelihood values. the last two rows at the bottom report the log-marginal likelihoods for synthetic data generated according to a student t and a gauss distribution.  we present a robust probabilistic approach to protein structure superposition and comparison. the approach builds on heavy-tailed distributions to model non-rigid displacements between protein structures. to estimate these distributions and an optimal superposition we employ a scale mixture representation of the heavy-tailed models. practically, this amounts to introducing weights for each atom position and to estimate the weights iteratively during structure superposition. in contrast to other weight-based superposition methods, the scale mixture framework provides a firm statistical basis for setting the weights. moreover, the link to the closed form helps to interpret the weighting scheme in terms of heavy-tailed models for structural displacements. 
the drawback is that under a heavy-tailed model analytical expressions for the optimal superposition no longer exist. we develop two iterative procedures, an expectation maximization algorithm and a gibbs sampler, to estimate the local weights, the optimal superposition, and the parameters of the heavy-tailed distribution. the standard dissimilarity measure for protein structures is the root mean square deviation of representative atom positions such as α-carbons. to circumvent this problem we work with a scale mixture representation, which implies a weighted rmsd. a scale mixture representation facilitates the fitting of these distributions and enables outlier-tolerant superposition. to do so, we model non-rigid displacements in protein structures with outlier-tolerant probability distributions. protein structure comparison is a central issue in structural bioinformatics. these distributions exhibit heavier tails than the gaussian distribution underlying standard rmsd minimization and thus accommodate highly divergent structural regions. we present a probabilistic method for robust superposition and comparison of protein structures. heavy-tailed distributions are well-suited to describe large-scale conformational differences in protein structures. applications demonstrate that heavy-tailed models capture differences between structures undergoing substantial conformational changes and can be used to assess the precision of nmr structures. our method aims to identify the largest structurally invariant core. how to evaluate optimal fits becomes a matter of debate, if the structures contain regions which differ largely - a situation encountered in nmr ensembles and proteins undergoing large-scale conformational transitions.
 protein structure comparison is a central issue in structural bioinformatics. the standard dissimilarity measure for protein structures is the root mean square deviation of representative atom positions such as α-carbons. to evaluate the rmsd the structures under comparison must be superimposed optimally so as to minimize the rmsd. how to evaluate optimal fits becomes a matter of debate, if the structures contain regions which differ largely - a situation encountered in nmr ensembles and proteins undergoing large-scale conformational transitions.  we present a probabilistic method for robust superposition and comparison of protein structures. our method aims to identify the largest structurally invariant core. to do so, we model non-rigid displacements in protein structures with outlier-tolerant probability distributions. these distributions exhibit heavier tails than the gaussian distribution underlying standard rmsd minimization and thus accommodate highly divergent structural regions. the drawback is that under a heavy-tailed model analytical expressions for the optimal superposition no longer exist. to circumvent this problem we work with a scale mixture representation, which implies a weighted rmsd. we develop two iterative procedures, an expectation maximization algorithm and a gibbs sampler, to estimate the local weights, the optimal superposition, and the parameters of the heavy-tailed distribution. applications demonstrate that heavy-tailed models capture differences between structures undergoing substantial conformational changes and can be used to assess the precision of nmr structures. by comparing bayes factors we can automatically choose the most adequate model. therefore our method is parameter-free.  heavy-tailed distributions are well-suited to describe large-scale conformational differences in protein structures. a scale mixture representation facilitates the fitting of these distributions and enables outlier-tolerant superposition. 