The Theory Behind Root Mean Squared Deviation (RMSD)

(note: this page was adapted, with permission, from Bosco Ho's discussion of RMSD, originally available at http://bosco.infogami.com/Root_Mean_Square_Deviation)

If it's your day job to push proteins in silico then you will one day have to calculate the RMSD of a protein. For example, you've just simulated the protein GinormousA for a whole micro-second, but you don't even know if GinormousA is stable. Sure, you could load up a protein viewer and eyeball the stability of the protein throughout the simulation. But eye-balling is just not going to cut it with your boss -- you need to slap a number on it. This is where RMSD comes in. RMSD is a quantitative measurement of difference between two atomic structures.

A molecule's conformation is basically a set of 3-dimensional coordinates. This means that a conformation is a set of $N$ vectors, $x_n$, where each $x_n$ has three components. As the molecular conformation changes, therefore, the result is a new set of vectors $y_n$. We can capture the the differences between conformations using these two sets of vectors.

To make life easier, we first shift the center of mass of $y_n$ and $x_n$ to the origin of the coordinate system (we can always shift them back afterwards). Once "centered," the differences between the vector sets can be captured by the total squared distance between them:

\[ E = \sum_n{|x_n - y_n|^2} \]

The root mean-square deviation (RMSD) is then:

\[ RMSD = \sqrt{\frac{E}{N}} \]

But hola, you say...this mean-square measure doesn't measure similarity very well! After all, any rotation of the set of $y_n$ (which doesn't change the internal arrangement of $y_n$), would distort the RMSD! What we really need, then, is to find the best rotation of $y_n$ with respect to $x_n$ before taking the RMSD. To rotate the vectors $y_n$, we apply a rotation matrix $U$ to get $y'_n = U \cdot y_n$. We can then express $E$ as a function of the rotation $U$:

\[ E = \sum_n{|x_n - U \cdot y_n|^2} \]

Thus, to find the RMSD, we must first find the rotation $U_{min}$ that minimizes $E$ -- the optimal superposition of the two structures.

Matching Sets of Vectors

The classic papers of Wolfgang Kabsch showed how to minimize $E$ to get this rotation. The algorithm described in those papers have since become part of the staple diet of protein analysis. As the original proof for the Kabsch equations (using 6 different Lagrange multipliers) is rather tedious, we'll describe a much simpler proof using standard linear algebra.

First, we expand $E$ (which is now a function of $U$):

\begin{eqnarray*} E &=& \sum_n{ |x_n|^2 + |y_n|^2 } - 2 \sum_n{ x_n \cdot U y_n }\\ &=& E_0 - 2 \sum_n{x_n \cdot U y_n} \end{eqnarray*}

The first part, $E_0$, is invariant with respect to a rotation $U$. The variation resides in the last term, $L = \sum_n{x_n \cdot U y_n}$. The goal, therefore, is to find the matrix $U_{min}$ that gives the largest possible value of $L$, from which we obtain:

\[ RMSD = \sqrt{\frac{E_{min}}{N}} = \sqrt{\frac{E_0 - 2 \cdot L_{max}}{N}} \]

Maximizing L

The trick here is to realize that we can promote the set of vectors $x_n$ and $y_n$ into two $N \times 3$ matrices $X$ and $Y$. Then we can write $L$ as a matrix trace:

\[ L = Tr( X^T U Y ) \]

where $X^T$ is a $3 \times N$ matrix and $Y$ is an $N \times 3$ matrix. Juggling the matrices inside the trace, we get:

\begin{eqnarray*} L &=& Tr( X^T U Y )\\ &=& Tr( U Y X^T )\\ &=& Tr( U R ) \end{eqnarray*}

Now, we still don't know what $U$ is, but we can study the matrix $R = YX^T$, which is the $3 \times 3$ correlation matrix between $X$ and $Y$. Furthermore, by invoking the powerful Singular Value Decomposition theorem, we know that we can always write R as:

\[ R = V S W \]

and obtain the $3 \times 3$ matrices $V$ and $W$ (which are orthogonal), and $S$ (which is positive diagonal, with elements $\sigma_i$). We can then substitute these matrices into the expression for $L$:

\begin{eqnarray*} L &=& Tr( U V S W^T )\\ &=& Tr( S W^T U V )\\ &=& Tr( S T ) \end{eqnarray*}

Because $S$ is diagonal, we can get a simple expression for $L$ in terms of the $\sigma_i$:

\[ L = Tr( S T ) = \sigma_1 \cdot T_{11} + \sigma_2 \cdot T_{22} + \sigma_3 \cdot T_{33} \]

Determining the Rotation

Writing L in this form, we can figure out the maximum value of $L$. The secret sauce is in $T$. If we look at the definition of $T$, we can see that it is a product of orthogonal matrices, namely, the matrices $U$, $W^T$ and $V$. Thus, $T$ must also be orthogonal.

Given this, the largest possible value for any element of $T$ is 1, or $T_{ij} \le 1$. And since $L$ is a product of terms linear in the diagonal elements of $T$, the maximum value of $L$ occurs when the $T_{ii} = 1$:

\[ L_{max} = Tr( S T_{max} ) = \sigma_1 + \sigma_2 + \sigma_3 = Tr( S ) \]

Plugging this into the equation for RMSD, we get

\[ RMSD = \sqrt{ \frac{E_0 - 2 (\sigma_1 + \sigma_2 + \sigma_3 )}{N} } \]

From this analysis, we can also determine the rotation that yields the optimal superposition of the vectors. When $T_{ii} = 1$, then $T = I$, the identity matrix. This immediately gives us an equation for the optimal rotation $U_{min}$ in terms of the left and right orthogonal vectors of the matrix R:

\[ T_{max} = W^T U_{min} V = I \]

\[ U_{min} = W V^T \]

How to Solve for R

In the previous section we've shown that RMSD can be expressed in terms of the singular values of $R$, but we still don't know how to find them! Never fear. Kabsch, who first derived the solution, showed that to solve for $\sigma_i$, we should solve the matrix $R^TR$, which is a nice symmetric real matrix, and is consequently much easier to solve. Expanding on the singular value decomposition from the previous section

\[ R^TR = (W^T S^T V^T)(V S W) = W^T S^2 W \]

we can see that the singular value term of $R^TR$ is $S^2$. In other words, the singular values of $R^TR$ ($\lambda_i$) are related to the singular values of $R$ ($\sigma_i$) by the relationship $\lambda_i = \sigma_i^2$. And also, the right singular vectors $W$ are the same for both $R$ and $R^TR$.

Furthermore, since $R^TR$ is a symmetric real matrix, its eigenvalues are identical to its singular values $\lambda_i$. Thus, we can use any standard method (e.g. the Jacobi transformation) to find the eigenvalues of $R^TR$, and by substituting $\lambda_i = \sigma_i^2$ into the RMSD formula:

\[ RMSD = \sqrt { \frac{ E_0 - 2 ( \sqrt{\lambda_1} + \sqrt{\lambda_2} + \sqrt{\lambda_3} )}{N} } \]

The Curse of Chirality

Unfortunately, there's one more ambiguity that we will have to consider before we get to the Kabsch equation for RMSD. This ambiguity arises from a degeneracy in the resulting $U_{min}$ from the SVD theorem. You see, $U_{min}$ is guaranteed to be orthogonal, but for our RMSD calculation, we don't just want orthogonal matrices, we want rotation matrices. And while orthogonal matrices can have determinants of +1 or -1, "proper" rotation matrices have determinants of +1 (unless they're also reflected, in which case they'll have determinants of -1, and are termed "improper" rotations).

When does this happen? To answer this question, take a look at your hands. Assuming that you're a primate, you have two of them, both alike, but reflected across the plane of symmetry that runs down the center of your body. Now, try to superimpose your hands. Go ahead. Try. If you're observant, you'll notice that you simply can't superimpose your hands. Even if you cut them off of your arms, the best you can hope to do is to align them so that their outlines are coincident -- but the palm of one hand will still be touching the other. You can never perfectly superimpose your hands, because your hands are chiral.

The same thing happens in molecular systems -- molecules are frequently chiral, and when two molecules are related by chirality, it sometimes happens that their optimal superposition involves a rotation, followed by a reflection. Thus, the "improper" rotation. Unfortunately, a reflection in $U_{min}$ will result in an incorrect RMSD value using the stated algorithm, so we need to detect this situation before it becomes a problem.

To restrict the solutions of $U_{min} = WV^T$ to proper rotations, we check the determinants of $W$ and $V$:

\[ \omega = det( V ) \cdot det( W ) \]

If $\omega = -1$ then there is a reflection in $U_{min}$, and we must eliminate the reflection. The easiest way to do this is to reflect the smallest principal axis in $U_{min}$, by reversing the sign of the third (smallest) eignevalue, which yields the final Kabsch formula for RMSD in terms of $\lambda_i$ and $\omega$:

\[ RMSD = \sqrt { \frac{ E_0 - 2 ( \sqrt{\lambda_1} + \sqrt{\lambda_2} + \omega \sqrt{\lambda_3} )}{N} } \]

Implementation in the BTK

The BTK implements an optimized version of the Kabsch algorithm in its leastsquares_superposition() methods. These methods solve a system of cubic equations to determine the eigenvalues of $R^TR$, rather than using an iterative approach (such as Jacobi). If only RMSD is desired (i.e. the calculation of a rotation matrix is unecessary), the get_rmsd() method is more efficient.

In addition to these methods, a get_trivial_rmsd() function is provided to calculate the RMSD value between two structures without attempting to superimpose their atoms.

Generated on Sun Jul 15 20:46:28 2007 for BTK Core by  doxygen 1.5.1