short communications Acta Crystallographica Section A Foundations of

A common measure of conformational similarity in structural bioinformatics is the minimum root mean square deviation (RMSD) between the coordinates of two macromolecules. In many applications, the rotations relating the structures are not needed. Several common algorithms for calculating RMSDs require the computationally costly procedures of determining either the eigen decomposition or matrix inversion of a 3x3 or 4x4 matrix. Using a quaternion-based method, here a simple algorithm is developed that rapidly and stably determines RMSDs by circumventing the decomposition and inversion problems.


Introduction
Orthogonal rotations are commonly used for comparing macromolecular structures, and the root mean square deviation (RMSD) is a natural metric for quantitating the similarity of two optimally rotated structures (Flower, 1999). Several least-squares algorithms have been proposed for finding the rotation that minimizes the squared distances between corresponding atoms in the structures. The most efficient of these methods requires the eigen decomposition (or inversion) of a key matrix constructed from the sums and squares of the atomic coordinates in the structures being compared. The popular method of Kabsch (1978) involves a 3 Â 3 key matrix, whereas the methods of Diamond (1988), Horn (1987) and Kearsley (1989) use quaternion representations of the rotations with 4 Â 4 key matrices. In each of these methods, the minimizing rotation is formed from the eigenvectors of the key matrices and the minimizing RMSD is found from the eigenvalues of the key matrices.
In Kabsch's method, improper rotation matrices (rotoinversions) may arise and one must verify that the determinant of the rotation matrix is positive (Kabsch, 1978). Quaternion-based methods lack this complication since quaternion rotations are always proper. Thus, with quaternion-based methods the minimum RMSD can be calculated using the eigenvalues of the key matrix alone, without knowing the eigenvectors. If only RMSDs are desired, the quaternion-based methods can save a considerable amount of computation by avoiding a complete eigen decomposition of the key matrix. Diamond (1988) proposed the fastest method known for determining minimum RMSDs using an iterative algorithm. Diamond's method requires one 3 Â 3 matrix inversion per cycle in order to find the largest eigenvalue of a quaternion-based 4 Â 4 matrix. While matrix inversions are generally expensive, the inversion of a 3 Â 3 matrix can be calculated analytically and relatively quickly. Diamond's method, however, is unstable when the minimizing rotation is near 180 . This problem is especially grave when the coordinate differences are of magnitude similar to the coordinates themselves, such as is common when superimposing relatively small fragments of proteins or nucleic acids.
Based on Horn's quaternion superposition method (Horn, 1987), which is relatively unknown in structural biology, an extremely fast, numerically stable and accurate algorithm is presented here for determining minimum RMSDs. This algorithm does not require separate consideration of special problematic cases and it bypasses the computationally costly diagonalization and inversions commonly used to find the needed eigenvalues. Rather, a quick and simple Newton-Raphson root-finding method is used to determine the appropriate eigenvalue from the characteristic polynomial of the key matrix. This quaternion-based characteristic polynomial (QCP) algorithm requires significantly less computation than the alternate methods, including Diamond's.

The superposition problem
A molecular structure consisting of N atoms can be represented as an N Â 3 matrix, where the three columns correspond to the x, y and z coordinates. The superposition objective is to find the orthogonal rotation and translation that minimizes the squared Euclidean distance between the rows of two matrices corresponding to the two structures being compared. The translational component of the problem can be removed from the outset by translating each structure so that its respective centroid is at the origin. When superimposing structure B onto structure A, in matrix notation the problem is to minimize the sum of squared errors E with respect to the orthogonal rotation R: where kXk F denotes the Frobenius (or Euclidean) norm of the matrix X, and X 0 denotes the transpose of X and tr X denotes the matrix trace of X. Upon expansion of equation (1), it can be shown that where G A is the inner product of structure A, and M is the inner product of the matrices A and B. The inner product matrix M is given by where 3. A quaternion-based characteristic quartic polynomial for calculating RMSDs Horn (1987) has shown that the sum of squared errors E can be minimized by solving for the largest positive eigenvalue max of a 4 Â 4 symmetric key matrix K: ðS xx þ S yy þ S zz Þ S yz À S zy S zx À S xz S xy À S yx S yz À S zy ðS xx À S yy À S zz Þ S xy þ S yx S zx þ S xz S zx À S xz S xy þ S yx ðÀS xx þ S yy À S zz Þ S yz þ S zy S xy À S yx S zx þ S xz S yz þ S zy ðÀS xx À S yy þ S zz Þ 2 6 6 4 3 7 7 5 : The eigenvector corresponding to the largest eigenvalue of K is a quaternion equivalent to the optimal rotation, where Eigenvalues are usually determined by diagonalization of the matrix. Alternatively, one can find the eigenvalues of a symmetric matrix by locating the roots of the characteristic polynomial of the matrix. With jXj denoting the determinant of matrix X, the characteristic equation for the key matrix K is a quartic given by yz S zy þ S yy S zx S xz þ S zz S xy S yx Þ À 8ðS xx S yy S zz þ S yz S zx S xy þ S zy S yx S xz Þ zy þ 2ðS yy S zz À S yz S zy Þ F ¼ ½ÀðS xz þ S zx ÞðS yz À S zy Þ þ ðS xy À S yx ÞðS xx À S yy À S zz Þ Â ½ÀðS xz À S zx ÞðS yz þ S zy Þ þ ðS xy À S yx ÞðS xx À S yy þ S zz Þ G ¼ ½ÀðS xz þ S zx ÞðS yz þ S zy Þ À ðS xy þ S yx ÞðS xx þ S yy À S zz Þ Â ½ÀðS xz À S zx ÞðS yz À S zy Þ À ðS xy þ S yx ÞðS xx þ S yy þ S zz Þ H ¼ ½ðS xy þ S yx ÞðS yz þ S zy Þ þ ðS xz þ S zx ÞðS xx À S yy þ S zz Þ Â ½ÀðS xy À S yx ÞðS yz À S zy Þ þ ðS xz þ S zx ÞðS xx þ S yy þ S zz Þ I ¼ ½ðS xy þ S yx ÞðS yz À S zy Þ þ ðS xz À S zx ÞðS xx À S yy À S zz Þ Â ½ÀðS xy À S yx ÞðS yz þ S zy Þ þ ðS xz À S zx ÞðS xx þ S yy À S zz Þ: Using these coefficients, the characteristic equation can be reduced to the following simple form: Although the coefficient equations may appear somewhat horrendous, coding them is straightforward. With careful programming, by avoiding duplicate evaluations of several terms (such as S yz À S zy which appears four times), calculating the coefficients from the sums found in equation (4) involves at most only 66 floating point operations (FLOPs). Horn (1987) proposed solving this characteristic equation using Ferrari's method (Abramowitz & Stegun, 1972) for quartics. While analytically exact, Ferrari's method has seen little use as it is computationally slow and notoriously numerically unstable. Additionally, Horn (1987) omits expressions for C 0 and an equation given for C 2 contains an important typographical error.
Fortunately in this case, the largest positive root of the characteristic polynomial given in equation (8) can be found very quickly and with excellent numerical stability using the Newton-Raphson algorithm. Successful use of the Newton-Raphson method requires a good initial guess for the desired root. From equation (3) and noting that the largest positive eigenvalue of the key matrix K is equal to trðMRÞ with the optimal rotation, the largest positive root of the characteristic polynomial attains its maximum possible value when the two structures are in identical conformations (i.e. when RMSD ¼ 0). Thus, the maximum possible value of the largest positive root is constrained to be less than or equal to the average of the inner products of the the two structures, ðG A þ G B Þ=2. Using this value as an initial guess in the Newton-Raphson algorithm is essentially foolproof, since the algorithm will converge first on the largest root of the polynomial. The Newton-Raphson QCP algorithm can be described by the following simple pseudocode: whileðj À old j < precisionÞ old ¼ ; where the first derivative of the polynomial is given by dPðÞ=d ¼ 4 3 þ 2C 2 þ C 1 . The Newton-Raphson QCP algorithm given above is much faster than other common methods that require determination of eigenvectors, and it is also faster than methods using eigen decomposition for only eigenvalues. The efficiency and speed of the QCP Newton-Raphson algorithm were compared to the popular method using Householder reduction to tridiagonal form followed by QL decomposition with implicit shifts (H-QL method; Golub & Van Loan, 1996). Each iteration of the QCP algorithm requires merely 11 FLOPS, of which none are square roots. In contrast, the H-QL algorithm requires a comparable number of total FLOPs but over two dozen costly square roots. Diamond's (1988) method is likewise much more efficient than the eigen decomposition methods and also requires no square roots. Each cycle of the Diamond iteration uses % 70 FLOPS (accounting for symmetry considerations), plus 15 FLOPS for set-up. Since both the QCP and Diamond's methods take an average of about five iterations to converge to a relative precision of 10 À6 , the QCP method uses only about one third of the operations of Diamond's method for eigenvalue determination. Thus, essentially all of the computational effort in the QCP method is expended in calculating the summations of the key matrix M given in equation (4), a step that is common to all methods and dependent on the number of atoms compared.
The QCP algorithm has been benchmarked on several computing platforms (OSX, Linux and Irix) for >10 8 independent superposition calculations with protein fragments ranging from 5 to 500 residues. In all cases, the same RMSD solution was found, within error, as the Kabsch and Kearsley methods. In contrast, especially with smaller fragments, Diamond's method converged on an incorrect solution (corresponding to the second largest eigenvalue) in a significant fraction of cases. In the author's experience, if the construction of the key matrices is ignored, the QCP algorithm is over four times as fast as Diamond's method and 30 to 70 times faster than the eigen decomposition methods. An ANSI C library implementing the QCP method is available at http://monkshood.colorado.edu/QCP/ or by request from the author. I thank B. K. P. Horn for helpful discussion, an anonymous referee for alerting me to Diamond's scalar method and the American Cancer Society for funding.