Received

Efficiency of Rotational Operators for Geometric Manipulation of Chain MoleculesChaok Seok* and Evangelos A. Coutsias†,*Department of Chemistry, College of Natural Sciences, Seoul National University, Seoul 151-747, Korea*E-mail: chaok@snu.ac.kr†Department of Mathematics and Statistics, University of New Mexico, Albuquerque, New Mexico 87131, USA*E-mail: vageli@Received July 23, 2007Geometric manipulation of molecules is an essential elementary component in computational modelingprograms for molecular structure, stability, dynamics, and design. The computational complexity oftransformation of internal coordinates to Cartesian coordinates was discussed before.1The use of rotationmatrices was found to be slightly more efficient than that of quaternion although quaternion operators havebeen widely advertised for rotational operations, especially in molecular dynamics simulations of liquids wherethe orientation is a dynamical variable.2 The discussion on computational efficiency is extended here to a moregeneral case in which bond angles and sidechain torsion angles are allowed to vary. The algorithm ofThompson3 is derived again in terms of quaternions as well as rotation matrices, and an algorithm with optimalefficiency is described. The algorithm based on rotation matrices is again found to be slightly more efficientthan that based on quaternions.Key Words : Rotation operator, Rotation matrix, Quaternion, Internal coordinatesIntroductionCartesian coordinate modeling is popular in molecular studies,for example, in molecular dynamics simulations of biomolecules.4,5In many other situations, internal coordi-nates are a natural choice when constraints on internal coordinates such as bond lengths and bond angles are to be implemented to reduce the number of degrees of freedom. Additional dihedral angle constraints from experiments or structures of homologous protein sequences can be easily incorporated in internal coordinates.If the dihedral angles are taken from the actual structure, and the bond lengths and bond angles are fixed at ideal values, the degree of deviation from the actual structure increases with the chain length. For example, a protein structure with 800 amino acids can show RMSD (Root-Mean-Square Deviation) as large as 20 Å if all the bond lengths and bond angles are replaced with ideal values. However, changes in the bond lengths and angles can be absorbed into small variations in dihedral angles to accurately represent the structure.In Ref. 6, it is shown that protein structures can be represented very accurately (within 0.3 Å) with ideal bond lengths and bond angles if dihedral angles are slightly modified to minimize RMSD.For this reason, strategies to take a reduced number of degrees of freedom in internal coordinates are popular in protein structure prediction and design studies.When internal coordinates are employed for geometry representation, conversion of the internal coordinates to Cartesian coordinates is necessary because most of the realistic energy functions or scoring functions involve terms that are conveniently evaluated in Cartesian coordinates. The transformation requires rotation operators: quatern-ions,7,8 rotation matrices,9 and other transformations10 have been used. (See Ref. 1 for more detailed discussions.) Alvarado and Kazerounian1showed that operation in rotation matrix is the most efficient for typical protein chains, and operation using quaternion in the matrix form is almost equivalent. A more general system than that in Ref. 1 is considered here in that bond angles are not required to be kept constant and side chains are not rigid. An efficient algorithm, which is equivalent to that of Thompson,3is derived by employing a simple reference state and perform-ing rotations in the reverse order. It is shown again that the method using rotation matrix is slightly more efficient than quaternion.A method that requires minimal number of rotation operations is described below. The method can be described in abstract terms, considering rotations as building blocks. The method is first described in the context of quaternion operators for convenience, and the same method is stated in terms of rotation matrices later. Efficiencies of the two methods using quaternions and rotation matrices are then compared.Method of Constructing a Chain Moleculewith Minimal Number of OperationsA linear chain involving N atoms is considered. Branches can be easily added once the backbone is constructed. This method can be combined with an analytical loop closure algorithm11to manipulate chains with internal loops, for example, polypeptides with disulfide bonds or cyclic chains.A linear chain can be fully described by N sets of bond lengths b i, bond angles θi, and dihedral angles φi as follows:, (1) P N = b i, θi, φi(), i = 1, …, N{}where the backbone atoms of the chain are numbered from 1to N , b i is the bond length between the nodes at i −1 and i , θi is the bond angle formed by the three nodes i −2, i −1, and i ,and φi is the dihedral angle defined by the four nodes i −3,i −2, i −1, and i . Six degrees of freedom, b 1, θ1, θ2, φ1, φ2, and φ3, do not affect the internal geometry of the chain, but define the absolute location and orientation of the chain in space with respect to ‘virtual’ nodes i =−2, −1, and 0.Therefore, the internal geometry of the chain can be fully described by N −1 bond lengths, N −2 bond angles, and N −3dihedral angles.To construct the molecule, simple rotations are applied starting from a simple initial arrangement. A straight, linear conformation which is laid along the positive x -axis is taken as the initial geometry. The initial arrangement is defined by.(2)Next, Cartesian coordinates for the nodes are determinedsuccessively. Placing successive atoms requires a rotation at each step, to the desirable θi and φi values. There is a “best”way for performing these rotations.Recall that 2,12 rotating a vector v about axis u by angle w is affected through conjugation by the quaternion q,(3)with q given by.(4)Here, ‘*’ is the operator of quaternion multiplication,defined below (see Eq. 20). The i -th link in the chain is placed by performing, in succession, all indicated rotations along the chain about the local axes u k by w k , where k = 1, 2,…, i .The net quaternion defining this rotation, q i , is constructed by the product of individual quaternions as,(5)where p 1 and p 2 can be set to unit quaternions (p *p =1) if the node 1 and 2 are to be placed on the x -axis for simplicity. It is important to note that in this construction, the axes of each successive rotation, u k are determined by the effect of allprevious rotations on the initial vector, as.(6)However, by carrying out the rotations in the reverse order, the same net rotation can be affected by utilizing the initial axes as.(7)The reason for this, obvious geometrically, can be also seen as follows: (8).is the multiplication of , rotation about the x -axisby φi , and p z (θi ), rotation about the z -axis by π−θi :. (9)This construction is understood as follows: start morphingthe chain from its end, node 1; rotate the bond vector (b i , 0,0) about the positive z -axis by angle π−θi (this enforces correct bond angle at node i ); rotate about the positive x -axis by angle φi to create the proper dihedral angle; then apply the rotation accumulated by placement of all the previous atoms.The rotations are simple, either around the x - or z -axis, and always applied to a vector on the x -axis. The algorithm can thus be implemented very efficiently. The calculation of the Cartesian coordinate R i for the i th node is now collected as(10), (11)where i =1, 2, …, N . R 0 and first few quaternions can be chosen to give simple arrangement for the first few atoms when only internal arrangement is of interest, or be chosen to place the chain in a precise location and orientation depending on the purposes.Comparison of Rotation Matrix and Quaternion Methods Rotation matrix method. One can construct a 3×3rotation matrix U i for each (θi ,φi ) pair, and multiply by a vector to get,(12)where r 0 is an arbitrary unit vector. The i -th rotation matrix U i depends on all the angles indexed 1 through i because the axis of rotation changes with previous angles. This can be avoided if rotations are performed in the reverse order:, (13)where the superscript 0 implies rotation about fixed axes.The rotation matrix can be expressed as a multiplication of two rotation matrices and , where U z is matrix for rotation by π−θi around the z -axis, and U x for rotation by φi around the x -axis. This simplification is madepossible by setting r 0 equal to a unit vector along the x -axis. The resulting isP N o= b i , π, 0(), i = 1, …, N {}v v ′→ = q *v *q cq = cos w 2---, u sin w2---⎝⎠⎛⎞q i = p i *…*p 2*p 1, i = 1, 2, …, N u k 0u k = q k 1–*u k 0*q k 1–cq i = p i *…*p 2*p 1 = p 10*p 20…*p i 0q i = p i *q i 1–= cos w i 2----, u i sin wi 2----⎝⎠⎛⎞*q i 1–= cos w i 2----, q i 1–*u i *q i 1–c[]sin w i 2----⎝⎠⎛⎞*q i 1–= q i 1–*cos w i 2----, u i 0 sin w i 2----⎝⎠⎛⎞*q i 1–c*q i 1–= q i 1–*p i 0p i 0p x φi ()p i 0= p x φi ()*p z θi ()= cos φi 2---, eˆx sin φi 2---⎝⎠⎛⎞*cos πθi –2-------------, e ˆz sin πθi –2-------------⎝⎠⎛⎞q i = q i 1–*cos φi 2---, eˆx sin φi 2---⎝⎠⎛⎞*cos πθi –2-------------, e ˆz sin πθi –2-------------⎝⎠⎛⎞R i = R i 1– + b i q i eˆx q i cR i = R i 1– + b i U i θi ,φi ,…,θ1,φ1()U i 1–θi 1–,φi 1–,…,θ1,φ1()…U 1θi ,φi ()r 0R i = R i 1– + b i U 10θ1,φ1()U 20θ2,φ2()…U i 0θi ,φi ()e ˆx U i 0U z θi ()U x φi ()eˆx U i 0.(14)Each node i is added to a growing chain with i −1 nodes by updating the accumulated rotation,,(15)applying the rotation on he vector (b i ,0,0), and translating the resulting vector by R i -1.The computational cost of this procedure is counted as follows:1. Calculation of cosines and sines, cos φi , sin φi , cos θi , and sin θi .2. Rotation update, (12 multiplications and 6 additions).3. Rotation update, (12 multiplications and 6 additions).4. First column of (no cost).5. Scaling by b i and shift by R i −1 (3 multiplications and 3additions).In fact, the above algorithm is equivalent to the algorithm described by Thompson 3 in 1967 using rotation matrices in a local coordinate system. A physical procedure of chain construction is provided here, and an algorithm using quaternions can be easily derived using this procedure, as shown above.The computational cost of adding an atom at a branching point is the same as the above if the rotation matrix V i has to be updated for later atoms connected to it. If an atom is terminal, only the first column of that rotation matrix is required, which is. (16)The cost of adding a terminal node is therefore 11 multi-plications and 6 additions for Eq. (16) and 3 multiplications and 3 additions for scaling by b i and shift by R i −1, in addition to the cost of calculating cosines and sines.Quaternion method. Quaternions, instead of rotation matrices, can be used to accumulate successive rotations as follows:, (17)where is given by Eq. (9). The four components of thequaternions and are,(18),(19)Noting 12 that multiplication of two quaternions a = (a 0, a )and b =(b 0, b ) is given by,(20)and the rotation matrix corresponding to quaternion q = (q 0,q 1, q 2, q 3) is, (21)the cost of computation for adding atom i is counted as follows:1. Division of the angles θi and φi by 2 (2 divisions).2. Calculation of cosines and sines, and .3. Quaternion update, (8 multiplications,4 additions).4. Quaternion update, (8 multiplications, 4additions).5. First column of U (q i )/2 (6 multiplications, 4 additions).6. Scaling by 2b i and shift by R i −1 (4 multiplications, 3additions).The computational cost of adding a terminal atom is the same as above because q i needs to be calculated to get the first column of U (q i )/2.Comparison of computational efficiency. The minimum number of mathematical operations required for calculation of the position of one atom to be added to a growing chain is summarized in Table 1.The total number of operations for constructing the whole chain is approximately the number of atoms times the above numbers. Note that this linear dependence on the number of atoms comes from the reverse rotation. If the operations cos and sin are assumed to be 20 times more expensive than +,−, ×, and ÷, which are assumed to be equal in computational time, the relative computational costs are 122and 123 for rotation matrix and quaternion methods,respectively. If, alternatively, the computation of cos and sin is performed using the familiar half-tangent formulasU i 0θi ,φi () = U x φi ()U z θi ()= 1000cos φi sin φi –0sin φi cos φi⎝⎠⎜⎟⎜⎟⎜⎟⎛⎞cos θi –sin θi –0sin θi cos θi –0001⎝⎠⎜⎟⎜⎟⎜⎟⎛⎞V i = U 10U 20…U i 1–0U i 0= V i 1–U i 0V i ′ = V i 1–U x φi ()V i = V i ′U z θi ()U i 0V i 1–cos θi –sin θi cos φisin θi sin φi ⎝⎠⎜⎟⎜⎟⎜⎟⎛⎞R i = R i 1– + b i p 10p 20…p i 0()eˆx b i p 10p 20…p i 0()cp i 0p x φi ()p z φi ()p x φi () = cos φi 2---, sin φi 2---, 0, 0⎝⎠⎛⎞p z θi () = cos πθi –2-------------, 0, 0, sin πθi –2-------------⎝⎠⎛⎞a *b = a 0b 0−a b ⋅, a 0b + b 0a + a b ×()U q () = 2q 02q 121/2–+q 1q 2q 0q 3–q 1q 3q 0q 2–q 1q 2+q 0q 3q 02q 221/2–+q 2q 3q 0q 1–q 1q 3q 0q 2–q 2q 3+q 0q 1q 02q 321/2–+⎝⎠⎜⎟⎜⎟⎜⎟⎜⎟⎛⎞cos φi 2---, sin φi 2---, cos θi 2---,sin θi2---q i ′ = q i 1–p x φi ()q i = q i ′p z θi ()Table 1. Comparison of the number of operations for the two methodsMethod + or −×÷cos sin Rotation matrix 1527022Quaternion1526222Table 2. Comparison of the number of operations for the twomethods when applied to two proteins of different sizes (cos or sin operations were not counted here.)Method Kinase Keratin Rotation matrix 3039176931Quaternion40549103759,(22)the costs of the two methods then become 96 and 95 FLOPs,respectively (again assuming the cost of the tangent as 20operations).Using rotation matrices or quaternions for accumulating rotations for non-terminal atoms are therefore almost equivalent in terms of computational cost, but quaternions require more operations for terminal atoms. Overall, rotation matrix calculation is more efficient, the exact cost depending on the number of terminal atoms. Table 2 shows a com-parison of the number of FLOPs when the two methods are applied to two proteins, Abl tyrosine kinase and heparin-binding growth factor, following Ref. 1. The kinase has 62amino acids: MNDPNLFV ALYDFV ASGDNTLSITKGEK-LRVLGYNHNGEWCEAQTKNGQGWVPSNYITPVNS.Keratin has 154 amino acids: MAAGSITTLPALPEDGG-SGAFPPGHFKDPKRLYCKNGGFFLRIHPDGRVDGVR EKSDPHIKLQLQAEERGVVSIKGVSANRYLAMKEDG RLLASKSVTDECFFFERLESNNYNTYRSRKYTSWV A LKRTGQYKLGSKTGPGQKAILFLPMSAKS. According to Table 2, the rotation matrix method is about 25% more efficient, because of the differences in FLOPs required to invert terminal atoms.DiscussionAn algorithm with minimal rotation operations for inter-conversions between internal and Cartesian coordinate systems for a molecular chain is described, starting from a simple, straight reference state and applying rotations on the chain backwards. The rotation matrix formula is equivalent to that of Thompson.3 The number of algebraic operations is reduced compared to the methods described in Ref. 1 due to the use of a simpler reference. In Ref. 1, the cost is 70operations (+, −, ×, and ÷) for calculation of operator and 15 for placement of atom when rotation matrices are used,and 74 and 15, respectively, when quaternions are used. In the algorithms described here, the cost is 36 and 6 with rotation matrices, and 37 and 6 with quaternions. (All methods require equal numbers of sin/cos evaluations). The conclusion that the use of rotation matrices results in slightly more efficient computation than the use of quaternions is maintained. Efficient transformation of coordinate systems can facilitate studies on protein structure especially when local constraints are incorporated in natural manners.Acknowledgement. EAC and CS would like to acknow-ledge the hospitality of the Dill group at UCSF during their stay. CS acknowledges grants from the Seoul R&BD pro-gram, and a grant from MarineBio21, Ministry of Maritime Affairs and Fisheries, Korea.References1.Alvarado, C.; Kazerounian, K. Prot. Eng. 2003, 16, 717.2.Rappaport, D. C. The Art of Molecular Dynamics Simulation ;Cambridge Univ. Press: New York, 1995.3.Thompson, H. B. J. Chem. Phys. 1967, 47, 3407.4.Ko, H.; Shim, G .; Kim, Y . Bull. Korean Chem. Soc. 2005, 26, 2001.5.Choi, Y .; Lee, J.; Cho, K. W.; Hwang, S.; Jeong, K.; Jung, S. Bull.Korean Chem. Soc. 2005, 26, 1203.6.Coutsias, E. A.; Seok, C.; Dill, K. A. J. Comput. Chem . 2004, 25,1849.7.Jain, A.; V aidehi, N.; Rodriguesz, G . J. Comput. Phys. 1993, 106, 258. 8.Kneller, G .; Hinsen, K. Phys. Rev. E 1994, 50, 1559.9.Dullweber, A.; Leimkuhler, B.; McLachlan, R. J. Chem. Phys .1997, 107, 5840.10.Manocha, D.; Zhu, Y .; Wright, W. Comput. Appl. Biol. Sci. 1995,11, 71.11.Coutsias, E. A.; Seok, C.; Jacobson, M. P.; Dill, K. A. J. Comput.Chem . 2004, 25, 510.12.Coutsias, E. A.; Romero, L. The Quaternions with Applications toRigid Body Dynamics ; Sandia National Laboratories, Technical Report SAND2004-0153, 2004.cos θ = 1u 2–1u 2+--------------, sin θ = 2u 1u 2+--------------, where u = tan θ2--⎝⎠⎛⎞。