ARTICLE

TOP ABSTRACT ARTICLE REFERENCES

The biological functions of proteins require internal motions which can be characterized by the atomic fluctuations and their cross-correlations (Brooks et al., 1988). The determination of the atomic fluctuations from molecular dynamics simulations is complicated by global rotation, which is difficult to remove since the moments of inertia change with time for a non-rigid object like a protein. The most common method for removing rotational motion is to minimize the root-mean-squared deviation (RMSD) through finite rotations of all configurations with respect to the first configuration (Kabsch, 1976; Brooks et al., 1983). This method, although exact for a rigid body, is only approximate for a flexible system. Furthermore, a question has been raised as to whether the RMSD should be minimized for the most rigid part of the system (Hünengeberger et al., 1995), or for the entire system (Ichiye and Karplus, 1991; Karplus and Ichiye, 1996). Different implementations have been found to yield significantly different results for long-range cross-correlations of the internal motions (Hünengeberger et al., 1995; Ichiye and Karplus, 1991; Karplus and Ichiye, 1996). Thus, it is of interest to examine the internal motions of proteins under the constraint of a zero-total angular momentum. One way of achieving this is by doing normal mode dynamics to determine the internal motions (Brooks and Karplus, 1983; Ichiye and Karplus, 1991). Here, we perform molecular dynamics at zero-angular momentum for a protein model and show that the results can be used to obtain a well-defined separation of the internal motions and the overall rotation.

The total angular and linear momentum can be constrained to zero by using a constant-energy molecular dynamics simulation in which the total angular and linear momentum are conserved (Allen and Tildesley, 1987). Consequently, if total angular and linear momentum are set to zero at the beginning of a simulation, a trajectory with zero total angular and linear momentum will be obtained if numerical errors do not accumulate. In this paper, we obtain such a trajectory for a simple model three-helix bundle protein (Zhou and Karplus, 1997; Zhou and Karplus, 1999a,b).

The model protein (Fig. 1) consists of 46 freely jointed beads each of which represents an amino acid residue that can interact with other residues via a square-well potential. The square-well depth is ( > 0) if the interaction pair are in contact in the global-minimum structure and is 0.7, otherwise. Details of the model and the discrete molecular dynamics method used for the simulations can be found in Zhou and Karplus (1999b). Despite the simplicity of the model, it possesses many essential features of proteins, such as the existence of both disordered and ordered globule states, a two-state folding transition and a low-temperature frozen state similar in structure to the native state (Zhou and Karplus, 1997). The folding behavior of the model has been used to study general aspects of the folding of -helical proteins (Zhou and Karplus, 1999a,b).

View larger version (93K): [in this window] [in a new window]   FIGURE 1   The global minimum structure of the model with helix I to III colored light to dark.

The method for obtaining a constant-energy trajectory is as follows. First, the value for the total energy of the system is calculated from constant-temperature simulations. Based on results from the constant-temperature simulations (Zhou and Karplus, 1997), the model is known to be in its native state at a reduced temperature of T* = 0.25 (T* = k_BT/, k_B, the Boltzmann constant and T, the temperature). The corresponding average total energy TE from a simulation at that temperature is found to be 212.37 from the equation

(1)

where KE and PE are the average kinetic and potential energies, respectively; N = 46 is the total number of residues. Note that six global rotational and translational degrees of freedom are explicitly removed from the kinetic energy in Eq. 1 and they do not contribute to the potential energy in this isolated system. The instantaneous total angular momentum L of the system about the center of mass is (Goldstein, 1980)

(2)

where r_i and v_i are the position and velocity of residue i, respectively and M is the mass of a residue; the same mass is assigned to all residues. The instantaneous angular velocity  of the system satisfies the equation

(3)

where I is the moment of inertia tensor; I_xx = M(r_i²  x_i²) and I_xy = Mx_iy_i; etc.; principal moments of inertia correspond to the eigenvalues of the inertia tensor I (Goldstein, 1980). In this paper, all quantities are reported in reduced units; for the total angular momentum and the inertia tensor, they are L* = L/ and I* = I/M², respectively. For time, t* = t. Here,  is the hard core diameter which is equal to 4.27 Å (Zhou and Karplus, 1997).

The initial configuration and velocities are obtained from a phase point on the equilibrium constant-temperature simulation at T* = 0.25. The total linear and angular momenta are removed by removing the linear (v_i/N) and angular ( × r_i) components of the velocities (v_i^new = v_i^old  _j v_j^old/N   × r_i) (Brooks et al., 1983). The resulting velocities are then scaled by a constant multiplication factor so that the total energy is equal to 212.37, the average total energy for the system (see above). The first 100 million steps (equilibration collisions) are discarded. The equilibrium trajectory is sampled every 10 reduced time units (~10,000 collisions). Each reduced time unit corresponds approximately to 1 ns (Zhou and Karplus, 1999a) based on the experimental time of collapse in protein folding (~1 µs) (Ballew et al., 1996; Muñoz et al., 1997).

The equilibrium simulation of the model protein lasts for 10⁵ reduced time units, which corresponds to roughly 102 million collisions. The average reduced temperature is found to be 0.255 and the average potential energy, PE, is 229.2. This is in good agreement with PE = 227.5 at the same average temperature obtained from a constant-temperature simulation. Independent constant-energy simulations (with different initial configurations and velocities but with the same total energy) yield essentially the same equilibrium results.

The total angular momentum is found to increase systematically from 10¹⁶ reduced units (the limit of numerical precisions) to 10¹² reduced units. To avoid possible error accumulations, the angular and linear momentum are reset to zero after every collision in the production simulation. Since the resetting changes only global rigid-body motions, the equilibrium averages of the thermodynamic properties (i.e., energy, heat capacity) are essentially unchanged, as expected. Due to numerical precision, the instantaneous total angular momentum is not exactly zero even though it is reset to zero at every collision. Fig. 2 shows a nearly symmetric essentially Gaussian distribution around zero with a width at half-height equal to 3.5 × 10¹⁶ reduced units for the instantaneous total angular momentum. This indicates that the numerical errors in the value of the angular momentum are random and of the order expected from the numerical precision. Thus, there should be no cumulative error that could lead to global rotations.

View larger version (14K): [in this window] [in a new window]   FIGURE 2   An example of the distribution of the reduced total angular momentum L*x, L*y, and L*z found in a constant-energy simulation during which the linear and angular momentum are set to zero after every collision.

Fig. 3 shows six instantaneous structures obtained from the constant-energy simulation under the constraint of zero-total angular and linear momentum. As time increases, the structure of the model protein appears to gradually "rotate" in a counter-clockwise direction away from the starting structure; no translational motion is observed. The various structures differ from the t* = 0 structure by only a very small (0.15  0.21) root-mean-square deviation after optimal superposition. Thus, the observed rotation is approximately "rigid-body" rotation, even though the system has essentially zero total angular momentum (see below). The rotational matrix elements used to minimize the RMSD with respect to the starting structure are shown as a function of time in Fig. 4. The overall rotational matrix elements are small for the first 1000 reduced time units and increase in a nonuniform manner.

View larger version (48K): [in this window] [in a new window]   FIGURE 3   The instantaneous structures of the constant-energy simulation at (a) t* = 0, (b) 2 × 104, (c) 4 × 104, (d) 6 × 104, (e) 8 × 104, and (f) 105, respectively. The three helices (I-III) are colored light to dark.

View larger version (26K): [in this window] [in a new window]   FIGURE 4   The nine rotational matrix elements, uij, that are obtained when the structure is rotated to minimize the RMSD with respect to the starting structure. The results of t* from 0 to 5000 reduced units are shown at the top. The results for the entire simulation are shown on the bottom.

Where does the long-time-scale rotation come from? The principal moments of inertia (I^*_xx, I^*_yy, and I^*_zz) for the starting configuration are 91.77, 96.12, and 60.38, respectively. To estimate the maximum amount of rotation caused by error accumulation, we consider a rigid system that rotates with the reduced angular momentum L^*_x = 0, L^*_y = 0, and L^*_z = 10¹³ for t* = 10⁵. The value 10¹³, which is two orders of magnitude larger than the maximum value of angular momentum shown in Fig. 2, is used to estimate an upper bound for the accumulated error. Such a rigid system would rotate over L^*_zt*/I^*_zz ~ 10⁸ degrees during the entire simulation. This is to be compared with the actual rotation of 42° (the Euler angle ) at t* = 10⁵. Thus, small numerical errors in angular momentum cannot explain the large rotation found in the simulation, even if they were to accumulate.

Another possibility is that the rotation arises from the presence of systematic errors in the program. To ensure that this is not the case, we simulated a rigid three-helix model using the same program. The rigid three-helix model corresponds to t* = 0 structure used in the equilibrium constant-energy simulation. The interaction between the pair of residues i and j is given by an infinitely deep square-well potential,

(4)

where R_ij is the distance between residues i and j in the t* = 0 structure and  is the bond flexibility parameter. All residues are connected via a nearly rigid bond. For  = 0, we have a rigid three-helix bundle with fixed distances between all pairs of residues. Over the time period used for the flexible protein model (10⁵ reduced time units), the rigid model with  = 0.01 rotates very little (Fig. 5 a); the rotation is only 2° in terms of the Euler angle . This is reflected by the fact that the elements of the rotational matrix for optimal superposition is either very close to 1 (u_ii) or 0 (u_ij) during the entire 10⁵ reduced time unit period. As the flexibility of the model increases ( increases from 0.01 to 0.1), relative rotation between different instantaneous structures starts to appear (Fig. 5 b). Since  = 0.1 is used for bond constraint in the actual model for the protein, this points to protein flexibility as the source of the time-dependent rotation shown in Figs. 3 and 4.

View larger version (24K): [in this window] [in a new window]   FIGURE 5   As in Fig. 4 but for the rigid model with the flexibility parameter (a)  = 0.01 and (b)  = 0.1, respectively.

That a flexible system can rotate with zero total angular momentum is well established and has been used to explain such everyday observations as that divers, gymnasts, and cats can rotate under overall zero-angular momentum conditions (Frohlich, 1980). For proteins, the movement of this flexible system at each time step can be decomposed conceptually into a rigid rotational motion in which the distance between all pairs of residues are unchanged and the individual translational motion of each residue that moves the residues from the rigidly rotated positions to their final positions. A rigid rotation under zero-total angular momentum is possible when the angular momentum associated with the rigid rotation is cancelled by the angular momentum associated with the individual atomic motions.

The zero-angular momentum trajectory from the protein model permits us to analyze how overall (rigid) rotations affect the cross-correlation, c_ij, of the internal motions. The cross-correlation of the internal motions for residues (atoms) i and j (Ichiye and Karplus, 1991) is defined by the equation

(5)

where   denotes configurational averages. Fig. 6 shows c_ij as a function of the average distance between residues i, j, r_ij, obtained via two different methods for removing rotational motions from the trajectory. When the rotational motion is removed by rotating each coordinate to minimize the RMSD for the entire molecule from the initial coordinate set (Fig. 6 a), the c_ij values show positive correlations at small and large separations and negative correlations at intermediate separations. This is in agreement with all-atom molecular dynamics simulations of bovine pancreatic trypsin inhibitor (BPTI) as well as normal mode results (Ichiye and Karplus, 1991). The use of rigidly packed residues (core residues here) as the reference frame for minimizing the RMSD changes the long range correlation to nearly zero (Fig. 6 b), a result similar to that obtained in all-atom simulations of BPTI and lysozyme by Hünengeberger et al. (1995).

View larger version (20K): [in this window] [in a new window]   FIGURE 6   The cross-correlation of residue motions as a function of the average distance between a pair of residues average over entire equilibrium simulation (t* = 0 to 105). Overall rotations are removed (a) via minimizing RMSD for all 46 residues and (b) via removing RMSD for the 9 residues that are closest to the center of mass on average.

In Fig. 7, the cross-correlation matrix c_ij is calculated directly from the zero-total angular momentum trajectory. The trajectory is analyzed for eight different time intervals: 0-1000, 0-2000, 0-3000, 0-4000, 0-5000, 0-10000, 0-50000 and 0-100000. The matrix elements for each time interval are calculated from the configurations obtained in that length of time. On the short time scale (up to 4000 reduced time units, ~4,000,000 collisions), the overall shape of the c_ij versus distance curve is essentially the same as that obtained by minimizing the RMSD for all 46 residues. As time increases, the correlation of c_ij at long distance changes from positive to negative and the distribution of c_ij at a specific distance becomes broader. A "pure" rigid rotational motion has large negative correlation between the parts that move in opposite directions. This leads to significant increases of negative correlations as the protein rotates away from its original position (Fig. 3). This overall rotation is "random" and, thus, it is expected to average to zero, if the system is simulated long enough to explore fully its rotational degrees of freedom. The similarity between c_ij in Fig. 6 a and the "short" time behavior of c_ij in Fig. 7 suggests that removing rotational motions by minimizing the RMSD for the entire system is more appropriate than doing it by minimizing the RMSD only for the more rigid part of the system, in accord with the conclusions of Karplus and Ichiye (1996) and Abseher and Nilges (1998).

View larger version (30K): [in this window] [in a new window]   FIGURE 7   The cross-correlation of residue motions as a function of the average distance between a pair of residues average over entire equilibrium simulation. The zero-angular momentum trajectory was used directly. The trajectory was analyzed for eight different lengths of time intervals: 0-1000, 0-2000, 0-3000, 0-4000, 0-5000, 0-10000, 0-50000 and 0-100000, as labeled.

We have shown that it is impossible to eliminate overall global rotation exactly in molecular dynamics simulations of a flexible protein-like model and by inference for a protein. Such behavior has been observed in a zero angular momentum simulation of BPTI (B. Brooks, unpublished). This is a consequence of the fact that even under the constraint of zero total angular momentum, global rotation still occurs over a long time period due to the flexibility of the protein. The short-time behavior of the zero total momentum trajectory, which is essentially rotation-free, suggests that removing rotations by minimizing the RMSD for the entire system is more appropriate than minimizing the RMSD for the more rigid part of the system.