help button home button Biophys. J.
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
Originally published as Biophys J. BioFAST on August 12, 2005.
doi:10.1529/biophysj.105.065664
This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow All Versions of this Article:
biophysj.105.065664v1
89/5/2939    most recent
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Van Wynsberghe, A. W.
Right arrow Articles by Cui, Q.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Van Wynsberghe, A. W.
Right arrow Articles by Cui, Q.
Biophysical Journal 89:2939-2949 (2005)
© 2005 The Biophysical Society

Comparison of Mode Analyses at Different Resolutions Applied to Nucleic Acid Systems

Adam W. Van Wynsberghe * and Qiang Cui * {dagger}

* Graduate Program in Biophysics and {dagger} Department of Chemistry and Theoretical Chemistry Institute, University of Wisconsin, Madison, Wisconsin

Correspondence: Address reprint requests to Q. Cui, Tel.: 608-262-9801; E-mail: cui{at}chem.wisc.edu.


   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 METHODS
 RESULTS AND DISCUSSION
 CONCLUDING REMARKS
 ACKNOWLEDGEMENTS
 REFERENCES
 
More than two decades of different types of mode analyses has shown that these techniques can be useful in describing large-scale motions in protein systems. A number of mode analyses are available and include quasiharmonics, classical normal mode, block normal mode, and the elastic network model. Each of these methods has been validated for protein systems and this variety allows researchers to choose the technique that gives the best compromise between computational cost and the level of detail in the calculation. These same techniques have not been systematically tested for nucleic acid systems, however. Given the differences in interactions and structural features between nucleic acid and protein systems, the validity of these techniques in the protein regime cannot be directly translated into validity in the nucleic acid realm. In this work, we investigate the usefulness of the above mode analyses as applied to two RNA systems, i.e., the hammerhead ribozyme and a guanine riboswitch. We show that classical normal-mode analysis can match the magnitude and direction of residue fluctuations from the more detailed, anharmonic technique, quasiharmonic analysis of a molecular dynamics trajectory. The block normal-mode approximation is shown to hold in the nucleic acid systems studied. Only the mode analysis at the lowest level of detail, the elastic network model, produced mixed results in our calculations. We present data that suggest that the elastic network model, with the popular parameterization, is not best suited for systems that do not have a close packed structure; this observation also hints at why the elastic network model has been found to be valid for many globular protein systems. The different behaviors of block normal-mode analysis and the elastic network model, which invoke similar degrees of coarse-graining to the dynamics but use different potentials, suggest the importance of applying a heterogeneous potential function in a robust analysis of the dynamics of biomolecules, especially those that are not closely packed. In addition to these comparisons, we briefly discuss insights into the conformational space available to the hammerhead ribozyme.


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
METHODS
 RESULTS AND DISCUSSION
 CONCLUDING REMARKS
 ACKNOWLEDGEMENTS
 REFERENCES
 
Studying biologically relevant problems via computational methods is always a compromise between the level of detail and accuracy needed to answer a specific question and the availability of computational resources. This is due to both the large size of biological polymers and the long timescales over which they exhibit interesting behavior. Processes such as protein folding or conformational change often can occur on the millisecond timescale, whereas conventional molecular dynamics (MD) simulations need to propagate in femtosecond time-steps to simulate atomic motion accurately. This necessitates 1012 evaluations of the energy and its gradient, which is an impractical if not impossible requirement given commonly available computing power. Consequently, considerable effort has been devoted to developing techniques to gain access to specific quantities of interest while treating other degrees of freedom that do not affect those quantities in less detail (1Go–5Go). The most commonly coarse-grained aspect of molecular simulation is solvent degrees of freedom (1Go,3Go), but other methods also approximate properties of the macromolecule. For example, to study enzyme mechanisms, hybrid QM/MM potential energy functions have been developed (6Go) that treat the chemically important residues quantum mechanically, whereas the remaining residues' electronic structure is approximated with a classical force field. Biased dynamics such as targeted MD (2Go) have been designed to force specific transitions to take place during nanosecond simulations; attempts have also been made (7Go–9Go) to recover the equilibrium thermodynamics using Jarzynski's equality (10Go) from such nonequilibrium simulations. The successes of these few examples vary, but in each case, simplifying approximations are made to access the desired properties. These approximations can only be useful if they do not affect the observable of interest. Obviously, this must be checked through careful test calculations using representative systems.

The computational studies in this work test the validity of the approximations made in several mode analyses, techniques that are generally used to identify mobile structural motifs in large macromolecular systems. Since this observable has major contributions from motions that occur at long timescales, only highly approximate methods are able to access this type of information for very large systems. One such method that has been successfully applied to protein systems for more than three decades is classical normal-mode analysis (CNMA) (11Go–14Go). Traditionally, this method represents the potential energy surface of a system by a harmonic approximation around a single minimum in an all-atom molecular mechanics force field. Construction and subsequent diagonalization of the mass-weighted second derivative matrix (the Hessian matrix) allows the equations of motion of this simplified system to be solved analytically (see Methods for details). Although this level of approximation is not appropriate for the analysis of detailed side-chain motions, useful information about large-scale motion and mobile structural motifs can be readily obtained; a number of studies (15Go–20Go) have illustrated that CNMA often gives rather reliable results in this context.

Unfortunately, even for an approximate method such as CNMA, the larger biological complexes such as the ribosome, ATPases, and RNA polymerases require such significant memory storage and time for diagonalization as to be intractable on conventional computers. To treat such systems, the blocked Hessian method was developed (21Go) and subsequently improved upon (17Go). This method transforms the Hessian into a reduced space spanned by the rotational and translational degrees of freedom of user-defined blocks, usually taken as a single residue or nucleotide. This reduced space results in a much smaller Hessian that can be more easily stored and diagonalized. The block approximation was found to be valid in calculations of residue fluctuations in protein systems where only the low-frequency modes are of interest. Currently known as the block normal-mode analysis, this technique has made it possible to study the thermal fluctuations of large macromolecular systems while still using an all-atom molecular mechanics force field. In addition, other methods have also been devised to deal with the problem of storage and diagonalization of large matrices (22Go).

For even larger systems or when computational power is limited and calculational speed is of the utmost importance, a further simplified method of mode analysis, the elastic network model (ENM), has been proposed (23Go–29Go). It uses a simplified potential to create a network of harmonic springs that connect atoms that are within a given cutoff distance. This method has been used successfully to study multiple RNA polymerases from different organisms (30Go) as well as the extremely large ribosomal multisubunit complex (31Go,32Go). As it is the least detailed method described here, it also has the most severe limitations. The arbitrary nature of the utilized spring constant removes any ability to predict the magnitude of thermal fluctuations, and the results also have the lowest resolution, usually at the residue level.

All of the above methods have been extensively tested for and used in protein systems. They have not been rigorously tested for nucleic acids, however. Since the early 1980s, it has been clear that nucleic acids, particularly RNA, can play cellular roles other than just genomic storage (33Go–36Go). Crystal and NMR structures of RNA molecules that exhibit enzymatic and regulatory functions have been numerous in the last 10 years (37Go–45Go). Discussions of the mechanisms and the mode of action of these systems have led to questions about the inherent flexibilities in these structures (46Go–54Go). Single-molecule experiments have suggested that the dynamics of even small RNA systems can be strikingly complex and intimately coupled to function (55Go). Computational efforts that could assist in understanding the flexibility present in theses systems have the potential to be very useful. It is not clear, however, which, if any, of the mode analyses can be reliably applied to nucleic acids. The stabilizing interactions in nucleic acids are mainly hydrogen-bonding and base-stacking, whereas in proteins hydrophobic effects dominate. Protein structures are generally globular in shape and densely packed, whereas nucleic acids tend to be nonspherical with elongated and loosely packed global structures. With these differing structural features and interactions, it is not known a priori if these model techniques, which have worked surprisingly well in proteins, are reliable for nucleic acids. To our knowledge, the only tests of mode analysis applied to nucleic acid systems are a work from Viduna et al. (56Go) and a work from Zacharias (57Go). The latter work compared a quasiharmonic analysis of a molecular dynamics trajectory of a simple double-stranded DNA octamer with a classical normal-mode calculation. In this article, we hope to extend this work to nucleic acid systems with more complex tertiary structures as well as to other types of mode analyses.


   METHODS
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
RESULTS AND DISCUSSION
 CONCLUDING REMARKS
 ACKNOWLEDGEMENTS
 REFERENCES
 
All calculations were performed with the CHARMM (58Go) suite of programs compiled to use the Cornell force field (59Go), normally associated with AMBER (60Go). This choice of force field was necessitated by our inability to acquire stable MD trajectories with the CHARMM27 nucleic acid force field (61Go) (data not shown). The CHARMM-formatted Cornell parameters were obtained from Dr. Tom Cheatham III. Nonbonded parameters for Mg2+ were adapted from Aqvist (62Go). Both equilibration and productive simulations were performed using periodic boundary conditions in the microcanonical ensemble. Electrostatic interactions were treated using the particle-mesh Ewald summation (63Go), and the van der Waals potential was smoothly shifted to zero at 12.0 Å (64Go).

Two different systems were studied in this work: the hammerhead ribozyme (PDB code #301D) (37Go) and the guanine riboswitch (PDB code #1U8D) (40Go) and are shown in Fig. 1. All four calculation types—quasiharmonic analysis of a molecular dynamics trajectory, classical normal-mode analysis, block normal-mode analysis, and the elastic network model—were applied to the hammerhead ribozyme. For the guanine riboswitch, only the latter three calculations were performed.



View larger version (48K):
[in this window]
[in a new window]
  		FIGURE 1  The hammerhead ribozyme (A) and the guanine riboswitch (B) are shown in surface representations. Atoms closer to the center of mass of each system are colored red, whereas those far away are colored blue with a linear gradient for intermediate distances. Note that the guanine riboswitch is more densely packed than the hammerhead ribozyme. This figure was created with VMD (88Go).

 
Molecular dynamics simulations and quasiharmonic analysis (QHA)
Nucleic acid systems, especially those with marginal stability such as the hammerhead ribozyme, must be carefully prepared to ensure a stable molecular dynamics (MD) simulation. In particular, because these systems are so highly charged, accurate treatment of the electrostatics is essential. To prepare the system, the hammerhead ribozyme along with its five associated Mg2+ ions was immersed in a rhombic dodecahedron (RHDO) of thermalized TIP3P waters (65Go) with its characteristic dimension set at 78.28 Å. Twenty-nine Na+ ions were added to neutralize the total charge via Solvate 1.0 (66Go).

A careful equilibration scheme was adapted from the work of Hermann et al. (67Go). The equilibration in this work was accomplished in four stages. Initially, the ribozyme and the counterion positions were held fixed, while the water coordinates as well as the lattice dimension of the RHDO unit cell were optimized with 2000 steepest-descent (SD) steps. The optimized RHDO dimension was 78.03 Å and remained fixed throughout the rest of equilibration and productive simulations. Keeping the constraints as before, the system was then subjected to 75 ps of dynamics starting at 50 K with the temperature increasing in 50 K increments every 5 ps until 300 K was reached. Temperature was checked and equilibrated (if necessary) every 1 ps. The second stage of equilibration removed the constraint from the counterions to allow their positions to adjust to the ribozyme and the water. Both counterion and water positions were optimized with 2000 SD steps followed by another 75 ps of heating dynamics as before. In the third stage of equilibration, the ribozyme was allowed to move but was restrained by a harmonic potential to allow unfavorable contacts to dissipate without disturbing the structure significantly. Initially a mass-weighted force constant of 25 kcal/mol/amu/Å2 was applied to all atoms in the ribozyme and 2000 SD steps were performed. Using the same heating procedure as before, the system was heated to 300 K. The harmonic constraints were then gradually relaxed with 75 ps of simulation at each distinct constraint value. Velocities were reassigned at 300 K between each constraint. The values for the force constant were (all in kcal/mol/amu/Å2): 20, 15, 10, 5, 1, 0.1, 0.01, and 0.001. The last stage of equilibration was 75 ps of simulation, with the entire system evolving without any constraints along with velocity rescaling every 1 ps if the temperature had veered >10° from 300 K. Productive simulation differed from this last stage of equilibration only by the absence of this velocity rescaling. Productive simulations for the hammerhead system were carried out for a total of 13.5 ns with coordinates saved every 100 fs. For the majority of this time, the system exhibited stable dynamics with an all-atom RMSD of ~2.75 Å, as shown in Fig. 2. This value is somewhat larger than is typically seen in protein simulations because these small RNA systems are fairly flexible. It should also be noted that the Stem I conformation observed in the crystal structure used in these simulations, PDB #301D, is not the only low-energy conformer of the hammerhead. A curved conformation of Stems I and II was observed in a crystal structure of a RNA-DNA hammerhead ribozyme-inhibitor complex (68Go). The evolution of the structure away from the conformation seen in 301D to this curved conformation is reflected in the increase in RMSD in the first nanosecond of simulation. At ~12.5 ns, the system developed a local instability in the unpaired region between Stems I and III, which then quickly caused alterations in the structure and a sharp increase in the RMSD. Removing both the beginning and the end of this trajectory still results in >10 ns of stable simulation in which the system evolves as a solution-phase ribozyme. Thus, we chose to perform the quasiharmonic analysis on only this center part of the trajectory.



View larger version (21K):
[in this window]
[in a new window]
  		FIGURE 2  All-atom RMSD calculated every 100 fs with respect to the crystal structure (301D) of the hammerhead ribozyme. The outlined box denotes the portion of the trajectory used for the quasiharmonic analysis.

 
Principal component analysis is a general statistical technique that allows one to analyze a distribution to extract information about its variances and the directions of those variances (69Go,70Go). When applied to molecular dynamics it is usually known as essential dynamics or quasiharmonic analysis (QHA). In this specific case, the distribution of interest is that of the atomic positions, with the variance relating to the thermal fluctuations of the system. QHA provides data like any of the other mode analyses discussed in this article: directions that correspond to certain modes of motion and associated variances that describe how energetically costly it is to move in those directions. QHA is considered the more exact of the techniques discussed in this article, since it is derived from a MD trajectory sampling an atomic and anharmonic potential energy function. Because of this, results from QHA are a harmonic approximation of the free energy surface, whereas the other methods discussed are harmonic approximations to a potential energy surface. Nevertheless, it is understood that the reliability of QHA depends on the quality of the force field as well as on the length of the MD simulation (71Go).

The directions of motion and associated variances from QHA are the eigenvectors and eigenvalues of the covariance matrix from the MD trajectory, C (72Go),

(1)
where qi is the ith mass-weighted displacement coordinate from the minimized structure and the brackets average over all time-steps. Fluctuations can be calculated directly from the MD trajectory.

Classical normal-mode analysis (CNMA)
Classical normal-mode analysis is the second most detailed calculation in this work, but it is far more approximate than QHA. It does not require a MD trajectory, so the calculations are far easier to accomplish and can be finished in less than a day for small systems. The only preliminary step for CNMA is a local energy minimization. This minimization was completed using cycles of the adapted-basis Newton-Raphson method with gradually decreasing harmonic constraints to remove local steric clashes without perturbing the structure significantly. Minimizations were considered complete when only six modes (corresponding to rotational and translational degrees of freedom) had frequencies less than or equal to zero. For the hammerhead ribozyme, this required minimizing the RMS energy gradient to 0.0001 kcal/mol per Å, whereas the guanine riboswitch was minimized to 0.001 kcal/mol per Å. This minimization step is necessary for CNMA so that the linear term in the Taylor expansion of the potential is zero (13Go). Truncating at the second-order leaves a system whose equations of motion can be analytically solved by diagonalizing the Hessian matrix, H,

(2)
where V is the potential and qi is a mass-weighted displacement coordinate. The diagonalization of this matrix reduces the problem from a system of 3N coupled harmonic oscillators, a rather difficult problem, to a system of 3N uncoupled harmonic oscillators, a trivial problem to solve. This procedure yields eigenvectors and eigenvalues that can be used to calculate atomic fluctuations using

(3)
where kB is Boltzmann's constant, T is the temperature, mk is the atomic mass, and Wki is the kth component of the eigenvector with associated eigenvalue {omega}i.

Block normal-mode analysis (BNMA)
Block normal-mode analysis was originally suggested by Tama et al. (21Go) and improved by Li and Cui (17Go), and is a technique that extends CNMA to larger systems. BNMA constructs the same Hessian, Eq. 2, as in CNMA, but then the Hessian is projected onto a space spanned by the rotational and vibrational degrees of freedom of predefined blocks, reducing the Hessian storage requirement by approximately a factor of 200 and the diagonalization time by a factor of 3000 if one nucleotide equals one block. The resulting eigenvectors are then projected back to the all-atom space and have the same form as the eigenvectors from CNMA. This procedure perturbs the magnitudes of the eigenvalues, but in a linear fashion for the low-frequency modes (17Go,21Go). A plot of CNMA eigenvalues versus BNMA eigenvalues returns a straight line whose slope can be used to scale the BNMA results; the scaling factor has been observed to be insensitive to the system being studied and correlates with the size of the defined block (17Go,21Go). Fluctuations are again calculated using Eq. 3, but with these scaled frequencies.

Elastic network model (ENM)
The elastic network model was originally suggested by Tirion (23Go) and is the most coarse-grained, and subsequently, is the least computationally expensive of all the methods discussed in this article. It requires no minimization and uses a simplified representation of the macromolecule and the intramolecular potential. In this work, each base was represented by a single sphere located at the position of its phosphorus atom. Any pair of spheres within a 20 Å cutoff radius, rcut, of each other was connected with a Hookean spring at its equilibrium position with a force constant of 0.95 kcal/mol/Å2 (24Go). Previous work has shown that slight modifications to the coarse-graining scheme, e.g., using a sphere for both the phosphate and the N1 or N9 base nitrogen or changing rcut, does not affect the results significantly (73Go). Fluctuations are calculated with Eq. 3. Because the spring constant used is entirely arbitrary, ENM cannot predict the absolute magnitudes of the fluctuations. To compare the trends in the ENM root mean-square fluctuations (RMSFs) to the other calculations, the ENM RMSFs must be scaled. A simple, uniform scale factor was chosen so that the sum of the ENM RMSFs equaled the sum of the CNMA RMSFs.

Methods of comparison
In addition to comparing RMSFs from different mode analyses, the directional nature of the modes can also be examined. One can compare the subspaces each set of low-frequency modes spans by examining the spanning coefficients as defined by (17Go)

(4)
where the Wi values correspond to individual modes from the appropriate technique. Obviously, if the sum in Eq. 4 went over all modes, then each SPANi would be equal to 1 and this metric would be useless. The useful comparison is for the low-frequency modes only, because these modes significantly contribute to the thermal fluctuations of the system. Therefore, the sum is limited to only a few modes, and in this work, only the 41 lowest-frequency modes are used.

One can make a more detailed comparison of the sets of modes by looking at the individual terms that are summed over in Eq. 4, the expansion coefficients. This quantity, given by the absolute value of the expression within the parentheses in Eq. 4, is simply the absolute value of the dot product between two eigenvectors and will be referred to as the overlap. Calculating all possible overlaps between two sets of modes gives rise to an overlap-matrix (74Go). If the two sets of modes were identical, the matrix would have ones on the minor diagonal and zeros everywhere else. Spanning coefficients can be generated from an overlap matrix by summing the squares of an individual column.

Note that the ENM eigenvectors have a different dimensionality than the other, all-atom calculations. To compare the eigenvectors of ENM with the other calculations directly, the components corresponding to the phosphorus atoms in the all-atom modes are abstracted and then renormalized.


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS AND DISCUSSION
CONCLUDING REMARKS
 ACKNOWLEDGEMENTS
 REFERENCES
 
Hammerhead ribozyme
The hammerhead ribozyme is a particularly useful and interesting system for the purpose of this work because it is relatively small and therefore more easily studied computationally with multiple techniques. It is also an attractive system because questions regarding its mechanism can be explored using molecular dynamics. Although this system has been studied for more than 15 years, the detailed mechanism is still not completely elucidated (75Go). In this work, we attempt to investigate one of the remaining questions: Is the P9/G10.1 Mg2+ able to play a role in the enzymatic mechanism?

It was observed that the abolishment of the binding site of P9/G10.1 Mg2+ via a phosphorthioate substitution results in a decrease in catalytic activity of nearly 1000-fold (51Go). However, various crystal structures have consistently placed the binding site of this Mg2+ nearly 20 Å away from the cleavable phosphodiester bond (37Go,68Go,76Go). Attempting to explain this discrepancy, Peracchi et al. proposed that the hammerhead ribozyme undergoes a conformational change that brings this Mg2+ close to the active site (51Go). To test this hypothesis, we ran 12.5 ns of molecular dynamics and analyzed the trajectory with quasiharmonic analysis. This type of analysis allows one to visualize the soft motions of a macromolecular system and thus determine which parts of the system are mobile (i.e., exhibit large thermal fluctuations) and which are rigid (i.e., small thermal fluctuations). In Fig. 3, the root mean-square fluctuations (RMSFs) as computed from the MD trajectory are mapped onto the phosphate backbone. As expected, the most mobile parts of the ribozyme are the three stem regions at the extremities away from the core. The region known as Domain I contains the active site and seems to be the most rigid as it has the lowest calculated fluctuations. Domain II has larger thermal fluctuations, but it is still not as mobile as the three stems. Using QHA, it is possible to not only calculate the magnitude of the nanosecond timescale fluctuation but also obtain its direction. Fig. 4, A and B, show end-structures of the two quasiharmonic modes with the highest variances (lowest associated frequencies). These two modes show that Stems I and II primarily flex into and out of the cleft between them with a motion much like a pair of pliers. Stem III's major direction of fluctuation is in and out of the plane formed by the three helices.



View larger version (24K):
[in this window]
[in a new window]
  		FIGURE 3  RMSFs from the MD trajectory for each phosphate atom is plotted onto a tube structure of the hammerhead ribozyme. The conformation corresponds to the average structure from the portion of the MD trajectory shaded in Fig. 2. The black sphere highlights the binding site of the P9/G10.1 Mg2+. This figure was created with MOLSCRIPT 2.1.2 (89Go).

 


View larger version (32K):
[in this window]
[in a new window]
  		FIGURE 4  The extremes of the quasiharmonic modes with the lowest ({omega} = 0.26 cm–1) (A) and the next-to-lowest ({omega} = 1.38 cm–1) (B) frequencies. In A, the mode is shown with a temperature of 300 K, whereas in B, the mode is shown with a temperature of 2500 K. Since mode 1 has a much lower frequency, the amplitude of motion allowed at a given temperature is much greater. B is shown at the elevated temperature to make the direction of motion easier to visualize. This figure was created with MOLSCRIPT 2.1.2 (89Go).

 
What can these simulations tell us about the possible large-scale motion in the hammerhead ribozyme? The major conclusion from these simulations is that Stems I and II are highly mobile as reflected by their large magnitude of thermal fluctuations. As has been suggested by earlier experiments, we directly observe the conformational freedom of these two helices (77Go). This motion is characterized by a hinge near the active site and allowed the distance between Stems I and II to vary on the order of 5 Å during the MD trajectory. Perhaps surprisingly, we see little motion in the active site or Domain II. The claim that, during catalysis, Domain II approaches Domain I remains controversial, however (78Go). In this simulation, we see no evidence that such a conformational transition takes place. Unfortunately, however, our simulation result is limited by its timescale. With only 12.5 ns of simulation, any conformation change or rare event that takes place on longer timescales cannot be observed. It is possible that, in order to bring Domains I and II into close proximity, a large reorganization of the hammerhead structure must take place. Recent cross-linking studies have suggested that this indeed may be the case (79Go).

Mode analysis comparison: hammerhead ribozyme
The primary goal of this work is to study how well the approximations made in various types of mode analysis perform for nucleic acid systems. Because BNMA is a direct approximation to CNMA, the appropriate test for BNMA is a comparison to CNMA. For consistency, all other techniques are compared to CNMA as well, although in principle, QHA is the most detailed and physically motivated technique presented in this work. QHA, CNMA, and BNMA are all-atom methods and their results depend on the quality of the molecule force field with QHA being most sensitive; ENM is not an all-atom technique, but is computationally efficient and involves the least amount of parameterization.

Since these techniques are primarily used to identify mobile regions in a macromolecular system, the first and most relevant comparison to make is an examination of the RMSFs. As shown in Fig. 5, the CNMA results match the MD RMSFs in a qualitative way, as the trends match up very well. The magnitude predicted by CNMA is in fairly good agreement but seem to underestimate the MD fluctuations. This is reasonable, given that the CNMA results come from a harmonic approximation to the potential energy surface. The results also show that BNMA matches both the trends and the magnitudes of the RMSFs calculated from the MD and CNMA very well, and suggest that the BNMA approximation at least holds for the hammerhead ribozyme. The last method, the ENM, with the popular parameterization, performs rather poorly. The overall magnitude of the ENM fluctuations is correct, but they have been scaled to allow easier comparison. They would not be predictive if the CNMA results were not already known (see Methods for details). Even the qualitative trends in the RMSFs were not reproduced. In particular, it misses the minimum in RMSF at residue index 8 and poorly miscalculates the trend and magnitude of Stem I of the enzyme strand.



View larger version (23K):
[in this window]
[in a new window]
  		FIGURE 5  RMSFs for the hammerhead ribozyme from the four mode analyses. The Residue Index is as follows: 1–15 corresponds to the phosphorus atom from nucleotides within the enzyme strand, going from 5' to 3' and 16–39 are from the substrate strand, also going from 5' to 3'. Note that the 5' end of each strand is terminated with a hydroxyl, not a phosphate, so that no fluctuation is plotted for this residue. (See Fig. 1 from Scott et al. (37Go) for accepted nucleotide numbering.)

 
As to the eigenvectors, Fig. 6 shows that the BNMA modes span the CNMA modes very well, with the first 11 modes having a spanning coefficient >0.9 with a slow decrease as the frequency of the mode increases. The QHA modes are nearly equally as impressive, illustrating that CNMA, and by extension BNMA, modes well represent the conformational space accessed by the MD trajectory. This again gives credence that both CNMA and BNMA are valid approximations for this system. The ENM results are not quite as impressive, however. The spanning coefficients for the lowest-frequency modes, the ones important in calculating the thermal fluctuations, are much lower than observed in the QHA or BNMA. The 10 lowest-frequency modes' spanning coefficients fluctuate at ~0.75 and then fall off. This raises serious concerns that ENM may not be reproducing the structural fluctuations observed in the MD trajectory.



View larger version (20K):
[in this window]
[in a new window]
  		FIGURE 6  Spanning coefficients (Eq. 4) for the hammerhead ribozyme. Only the 41 lowest-frequency modes are included in the sum for this quantity.

 
As seen in Fig. 7 A, the individual CNMA modes match up reasonably well with those from the QHA. Because BNMA is a direct approximation of CNMA, the correlation between the two sets is also very good, as the overlap matrix in Fig. 7 B shows high overlap only near the diagonal, especially at the low-frequency modes. Fig. 7 C, however, shows the very poor overlap between the ENM modes and the CNMA modes. There is a weak correlation between the two sets at the very lowest frequencies, but very little elsewhere. It should be noted that these comparisons are only meaningful because the eigenvalues associated with the discussed eigenvectors are well-correlated (data not shown).



View larger version (103K):
[in this window]
[in a new window]
  		FIGURE 7  Overlap matrices for the hammerhead ribozyme. Perfect correspondence between two techniques would result in a black line along the minor diagonal. The scale along each axis is 41 modes.

 
The most immediate conclusion from these data is that CNMA seems to be a reliable way to identify mobile regions (i.e., exhibit large thermal fluctuations) in nucleic acids. Even with its harmonic approximation and the lack of explicit solvent, the RMSFs from CNMA predict which parts of the hammerhead ribozyme are mobile and which are rigid, and do so in a semiquantitative manner. The directional motion is also reproduced, as is illustrated in both the spanning coefficients and the overlap matrices. The data suggest, perhaps unsurprisingly, that the BNMA approximation holds up well, since BNMA results match those from CNMA. Unfortunately, the same success does not hold for the ENM. The RMSFs have some qualitative inaccuracies and the modes generated from ENM are not very similar to the other three methods. It is perhaps not too surprising that CNMA, and therefore BNMA, are valid for nucleic acids as well as for proteins systems. The changes in the interactions are coded into the parameters used in the actual force field. Obviously the changes in relative weights of certain types of interactions do not perturb the validity of the harmonic approximation used in CNMA. In the ENM, however, there is only one parameter that actually matters: the cutoff distance, rcut. Physically, the ENM can be viewed as a density detector; any atoms within rcut are included in a measure of the density of adjacent atoms. The denser the packing, the more springs are added, and the more rigidly the given atom is held. This correlates well to the trend in thermal fluctuations in proteins where, for the most part, the systems have a dense core and the mobile regions are protrusions. Previous work has shown that ENM has difficulties dealing with protrusions that are pathologically extended from the majority of the protein (80Go). In the case of the hammerhead ribozyme, nearly the entire system is a loosely packed, extended system. This is generally true for most small nucleic acid systems, where the burial of hydrophobic surface does not dominate the folding, as is often the case with proteins. We hypothesize that it is this irregularly packed structure, on the scale of rcut, that causes the breakdown of the ENM in the hammerhead system. To test this hypothesis, we performed the same calculations (excepting the MD simulation) on a more densely packed system, the guanine riboswitch.

Mode analysis comparison: guanine riboswitch
The guanine riboswitch has a more compact structure as compared to the hammerhead ribozyme, as shown in Fig. 1. Instead of three elongated double-helical regions made from two separate RNA strands as in the hammerhead, the riboswitch is made from a single self-complementary RNA strand that forms a tightly packed pancake-shaped structure. For the purposes of this work, this system is a small RNA like the hammerhead but is somewhat more compact. Using the same analysis tools as discussed in the above section on the hammerhead ribozyme, we show that the ENM is able to predict atomic fluctuations for this densely packed system more accurately.

As illustrated in Fig. 8, the BNMA calculations match both the magnitude and the trends of the CNMA RMSFs very well for the guanine riboswitch. The quality of this match is similar to that seen in Fig. 5 for the hammerhead ribozyme. However, unlike in the hammerhead ribozyme, the ENM results have no major inaccuracies and for the most part capture the trends from the CNMA well. Although there are some minor discrepancies, in this measure, the ENM performs much better for the guanine riboswitch than in the hammerhead ribozyme.



View larger version (20K):
[in this window]
[in a new window]
  		FIGURE 8  RMSFs for the guanine riboswitch. Only three curves are shown, since no MD trajectory was run for this system. The Residue Index corresponds to the single RNA strand going from 5' to 3'. This structure also lacked a 5' phosphate group, and no fluctuation from the 5' terminal nucleotide is plotted.

 
As seen in the hammerhead system, Fig. 9 shows that BNMA approximates the CNMA riboswitch results very well with high spanning coefficients, especially at low frequencies. The ENM results still do not match those from BNMA, and are similar in quality as for the hammerhead ribozyme. The BNMA/CNMA overlap matrix shown in Fig. 10 A illustrates that the two sets of modes are, as expected, highly correlated. Fig. 10 B shows that the ENM results are still not very good, but at low frequencies, the ENM results do correlate with CNMA very well. This correlation quickly dies off as the frequency of the modes increases, however. Overall, the ENM results from the guanine riboswitch are better than those from the hammerhead ribozyme. RMSFs for the riboswitch match up qualitatively with CNMA, whereas those for the hammerhead have some discrepancies. The ENM for the riboswitch still does not do a great job at predicting the direction of the structural mobility, but ENM applied to the hammerhead is worse. ENM does not work particularly well with either system, though, when compared to BNMA. From these data and previous analysis using a range of different parameters in the ENM model (73Go), we suggest that the ENM does not work well for systems with loosely packed structures, and therefore is not an appropriate choice for small RNA systems.



View larger version (17K):
[in this window]
[in a new window]
  		FIGURE 9  Spanning coefficients (Eq. 4) for the guanine riboswitch. Only the 41 lowest-frequency modes are used to calculate this quantity.

 


View larger version (65K):
[in this window]
[in a new window]
  		FIGURE 10  Overlap matrices for the guanine riboswitch. The scale along each axis is 41 modes.

 

   CONCLUDING REMARKS
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS AND DISCUSSION
 CONCLUDING REMARKS
ACKNOWLEDGEMENTS
 REFERENCES
 
A stable, 12.5-ns MD simulation of the hammerhead ribozyme was performed, the longest reported hammerhead simulation by >10 ns. Analysis of this trajectory allows us to identify mobile structural motifs in this small RNA system. Our results suggest that the most mobile regions are the three stems, while the active site is fairly rigid. We can also suggest that, in order for Domain II to approach the active site, a fairly extensive reorganization of the hammerhead must take place.

We have also used this MD trajectory as a benchmark to test the validity of several types of mode analyses in a complex nucleic acid system, within the limit of 10-ns-scale simulation time. The widely used CNMA provided qualitative and semiquantitative prediction of the magnitude of thermal fluctuations and the directionality of the soft motions. The BNMA approximation also was shown to hold for the two small RNAs studied here, and we expect that this observation will transfer to other RNA systems. The ENM model, however, had difficulties in matching the trends in fluctuations from the MD trajectory and the other two mode analyses based on physical potentials; this result was not found to be sensitive to the variation in the cutoff or identity of nodes in the network (73Go), although it is possible that more extensive reparameterization of the model will improve the results (I. Bahar, private communication). Another set of mode analyses performed on the guanine riboswitch suggests that the reliability of the ENM model may be improved for larger, more compact systems; indeed, the agreement between ENM and BNMA for 30S and 50S ribosome was found to be quite impressive (73Go).

It is interesting to discuss these observations briefly, in the context of constructing appropriate coarse-grained models for biological molecules (81Go). This is a topic of tremendous interest considering the rapid development of structural biology, which has made structures of biological systems at various resolutions and of increasing size and complexity readily available. In an increasing number of research articles, it has been demonstrated that a simple model like the elastic network model works extremely well for many large macromolecules, even when only low-resolution structural data were available (82Go,83Go). An interesting question, then, is "why does ENM work?" On physical grounds, one may argue that global flexibility is, by definition, a collective property and therefore is not expected to be sensitive to fine structural details but depends mainly on the overall shape of the macromolecule. This argument explains why models like BNMA, which neglect motions at small scales but maintain physical interactions, work well for many protein and nucleic acid systems. It remains intriguing why ENM, which assumes a homogeneous interaction strength (i.e., uniform spring constant in different parts of the macromolecule), works well for many protein complexes. The observation in this work that the quality of ENM depends on the compactness of the RNA system offers an interesting clue: in compact systems, such as most globular proteins (84Go), interatomic interactions, on average, are fairly homogeneous at the residue-level. That BNMA—which coarse-grained the dynamics of the system to a similar degree to ENM but kept physical interactions—produced reliable results, further corroborates the importance of heterogeneous interactions, even at the residue level, in systems that are not compact. A more quantitative connection between the compactness of the system and appropriate level of coarse-graining is under study. Finally, the fact that small changes in sequence (due to either mutation or organism origin) at hinge residues can often have profound impact on the displacement along functionally important direction of biomolecules (e.g., point mutations can abolish the motility of molecular motors (85Go,86Go)) also suggests the importance of heterogeneous interactions, although it is possible that such mutations mainly modify the free energy cost for moving along a particular direction instead of significantly modifying the character of low-frequency modes. Recent work has shown that these hinges can be identified by introducing nonhomogeneous interactions into the ENM (87Go), which actually relies on the variation of low-frequency modes upon changing the elastic force constants associated with specific residues (87Go). In other words, although homogeneous ENM is indeed incredibly useful in many studies, it is also of great interest to develop coarse-grained models that go beyond the limit of the homogeneous elastic model for robustly describing the flexibility of large biomolecular systems, for which the resolution of experimental characterization is typically low and the cost of calculations is an important factor.


   ACKNOWLEDGEMENTS
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS AND DISCUSSION
 CONCLUDING REMARKS
 ACKNOWLEDGEMENTS
REFERENCES
 
The authors thank Tom E. Cheatham 3rd for helpful discussions regarding nucleic acid simulations as well as coordinates used for initial testing of simulation protocols. The authors also thank Sam Butcher for suggesting the guanine riboswitch as a particularly compact small RNA structure for study. Guohui Li's assistance with the block normal-mode and elastic network mode code was also greatly appreciated. Discussions with Profs. I. Bahar and C. Brooks are appreciated.

A.W.V.W. was supported by a National Science Foundation predoctoral fellowship. Q.C. is an Alfred P. Sloan Research Fellow. Computational resources from the National Center for Supercomputing Applications at the University of Illinois are greatly appreciated.

Submitted on May 4, 2005; accepted for publication August 1, 2005.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 METHODS
 RESULTS AND DISCUSSION
 CONCLUDING REMARKS
 ACKNOWLEDGEMENTS
 REFERENCES
 
1. Feig, M., and C. L. Brooks. 2004. Recent advances in the development and application of implicit solvent models in biomolecule simulations. Curr. Opin. Struct. Biol. 14:217–224.[CrossRef][Medline]

2. Schlitter, J., M. Engels, and P. Kruger. 1994. Targeted molecular-dynamics—a new approach for searching pathways of conformational transitions. J. Mol. Graph. 12:84–89.[CrossRef][Medline]

3. Im, W., S. Berneche, and B. Roux. 2001. Generalized solvent boundary potential for computer simulations. J. Chem. Phys. 114:2924–2937.[CrossRef]

4. Shurki, A., and A. Warshel. 2003. Structure/function correlations of proteins using MM, QM/MM, and related approaches: methods, concepts, pitfalls, and current progress. Adv. Protein Chem. 66:249–313.[Medline]

5. Gao, J. L. 1996. Hybrid quantum and molecular mechanical simulations: an alternative avenue to solvent effects in organic chemistry. Acct. Chem. Res. 29:298–305.[CrossRef]

6. Field, M. J., P. A. Bash, and M. Karplus. 1990. A combined quantum-mechanical and molecular mechanical potential for molecular-dynamics simulations. J. Comput. Chem. 11:700–733.[CrossRef]

7. Rodriguez-Gomez, D., E. Darve, and A. Pohorille. 2004. Assessing the efficiency of free energy calculation methods. J. Chem. Phys. 120:3563–3578.[CrossRef][Medline]

8. Park, S., and K. Schulten. 2004. Calculating potentials of mean force from steered molecular dynamics simulations. J. Chem. Phys. 120:5946–5961.[CrossRef][Medline]

9. Jensen, M. O., S. Park, E. Tajkhorshid, and K. Schulten. 2002. Energetics of glycerol conduction through aquaglyceroporin GlpF. Proc. Natl. Acad. Sci. USA. 99:6731–6736.[Abstract/Free Full Text]

10. Jarzynski, C. 1997. Nonequilibrium equality for free energy differences. Phys. Rev. Lett. 78:2690–2693.[CrossRef]

11. Brooks, B., and M. Karplus. 1983. Harmonic dynamics of proteins—normal modes and fluctuations in bovine pancreatic trypsin-inhibitor. Proc. Natl. Acad. Sci. USA. 80:6571–6575.[Abstract/Free Full Text]

12. Go, N., T. Noguti, and T. Nishikawa. 1983. Dynamics of a small globular protein in terms of low-frequency vibrational modes. Proc. Natl. Acad. Sci. USA. 80:3696–3700.[Abstract/Free Full Text]

13. Hayward, S. 2001. Normal mode analysis of biological molecules. In Computational Biochemistry and Biophysics. O.M. Becker, A.D. MacKerell Jr., B. Roux, and M. Watanabe, editors. Marcel Dekker, New York.

14. Levitt, M., C. Sander, and P. S. Stern. 1983. The normal modes of a protein—native bovine pancreatic trypsin inhibitor. Int. J. Quantum Chem. 10:181–199.[CrossRef]

15. Ma, J. P. 2005. Usefulness and limitations of normal mode analysis in modeling dynamics of biomolecular complexes. Structure. 13:373–380.[Medline]

16. Cui, Q., G. H. Li, J. P. Ma, and M. Karplus. 2004. A normal mode analysis of structural plasticity in the biomolecular motor F-1-ATPase. J. Mol. Biol. 340:345–372.[CrossRef][Medline]

17. Li, G. H., and Q. Cui. 2002. A coarse-grained normal mode approach for macromolecules: an efficient implementation and application to Ca2+-ATPase. Biophys. J. 83:2457–2474.[Abstract/Free Full Text]

18. Krebs, W. G., V. Alexandrov, C. A. Wilson, N. Echols, H. Y. Yu, and M. Gerstein. 2002. Normal-mode analysis of macromolecular motions in a database framework: developing mode concentration as a useful classifying statistic. Proteins. 48:682–695.[CrossRef][Medline]

19. Thomas, A., K. Hinsen, M. J. Field, and D. Perahia. 1999. Tertiary and quaternary conformational changes in aspartate transcarbamylase: a normal mode study. Proteins. 34:96–112.[CrossRef][Medline]

20. Seno, Y., and N. Go. 1990. Deoxymyoglobin studied by the conformational normal mode analysis. I. Dynamics of globin and the hemoglobin interaction. J. Mol. Biol. 216:95–109.[CrossRef][Medline]

21. Tama, F., F. X. Gadea, O. Marques, and Y. H. Sanejouand. 2000. Building-block approach for determining low-frequency normal modes of macromolecules. Proteins. 41:1–7.[Medline]

22. Mouawad, L., and D. Perahia. 1993. Diagonalization in a mixed basis—a method to compute low-frequency normal modes for large macromolecules. Biopolymers. 33:599–611.[CrossRef]

23. Tirion, M. M. 1996. Large amplitude elastic motions in proteins from a single-parameter, atomic analysis. Phys. Rev. Lett. 77:1905–1908.[CrossRef][Medline]

24. Doruker, P., R. L. Jernigan, and I. Bahar. 2002. Dynamics of large proteins through hierarchical levels of coarse-grained structures. J. Comput. Chem. 23:119–127.[CrossRef][Medline]

25. Haliloglu, T., I. Bahar, and B. Erman. 1997. Gaussian dynamics of folded proteins. Phys. Rev. Lett. 79:3090–3093.[CrossRef]

26. Bahar, I., A. R. Atilgan, and B. Erman. 1997. Direct evaluation of thermal fluctuations in proteins using a single-parameter harmonic potential. Fold. Des. 2:173–181.[CrossRef][Medline]

27. Tama, F., and Y. H. Sanejouand. 2001. Conformational change of proteins arising from normal mode calculations. Protein Eng. 14:1–6.[Abstract/Free Full Text]

28. Atilgan, A. R., S. R. Durell, R. L. Jernigan, M. C. Demirel, O. Keskin, and I. Bahar. 2001. Anisotropy of fluctuation dynamics of proteins with an elastic network model. Biophys. J. 80:505–515.[Abstract/Free Full Text]

29. Doruker, P., A. R. Atilgan, and I. Bahar. 2000. Dynamics of proteins predicted by molecular dynamics simulations and analytical approaches: application to {alpha}-amylase inhibitor. Proteins. 40:512–524.[CrossRef][Medline]

30. Delarue, M., and Y. H. Sanejouand. 2002. Simplified normal mode analysis of conformational transitions in DNA-dependent polymerases: the elastic network model. J. Mol. Biol. 320:1011–1024.[CrossRef][Medline]

31. Tama, F., M. Valle, J. Frank, and C. L. Brooks. 2003. Dynamic reorganization of the functionally active ribosome explored by normal mode analysis and cryo-electron microscopy. Proc. Natl. Acad. Sci. USA. 100:9319–9323.[Abstract/Free Full Text]

32. Wang, Y. M., A. J. Rader, I. Bahar, and R. L. Jernigan. 2004. Global ribosome motions revealed with elastic network model. J. Struct. Biol. 147:302–314.[CrossRef][Medline]

33. Cech, T. R. 2002. Ribozymes, the first 20 years. Biochem. Soc. Trans. 30:1162–1166.[CrossRef][Medline]

34. Guerriertakada, C., K. Gardiner, T. Marsh, N. Pace, and S. Altman. 1983. The RNA moiety of ribonuclease-P is the catalytic subunit of the enzyme. Cell. 35:849–857.[CrossRef][Medline]

35. Cech, T. R., A. J. Zaug, and P. J. Grabowski. 1981. In vitro splicing of the ribosomal-RNA precursor of Tetrahymena—involvement of a guanosine nucleotide in the excision of the intervening sequence. Cell. 27:487–496.[CrossRef][Medline]

36. Winkler, W., A. Nahvi, and R. R. Breaker. 2002. Thiamine derivatives bind messenger RNAs directly to regulate bacterial gene expression. Nature. 419:952–956.[CrossRef][Medline]

37. Scott, W. G., J. B. Murray, J. R. P. Arnold, B. L. Stoddard, and A. Klug. 1996. Capturing the structure of a catalytic RNA intermediate: the hammerhead ribozyme. Science. 274:2065–2069.[Abstract/Free Full Text]

38. Serganov, A., Y. R. Yuan, O. Pikovskaya, A. Polonskaia, L. Malinina, A. T. Phan, C. Hobartner, R. Micura, R. R. Breaker, and D. J. Patel. 2004. Structural basis for discriminative regulation of gene expression by adenine- and guanine-sensing mRNAs. Chem. Biol. 11:1729–1741.[CrossRef][Medline]

39. Egli, M. 2004. Nucleic acid crystallography: current progress. Curr. Opin. Chem. Biol. 8:580–591.[CrossRef][Medline]

40. Batey, R. T., S. D. Gilbert, and R. K. Montange. 2004. Structure of a natural guanine-responsive riboswitch complexed with the metabolite hypoxanthine. Nature. 432:411–415.[CrossRef][Medline]

41. Adams, P. L., M. R. Stahley, A. B. Kosek, J. M. Wang, and S. A. Strobel. 2004. Crystal structure of a self-splicing group I intron with both exons. Nature. 430:45–50.[CrossRef][Medline]

42. Ke, A. L., K. H. Zhou, F. Ding, J. H. D. Cate, and J. A. Doudna. 2004. A conformational switch controls hepatitis {delta}-virus ribozyme catalysis. Nature. 429:201–205.[CrossRef][Medline]

43. Krasilnikov, A. S., X. J. Yang, T. Pan, and A. Mondragon. 2003. Crystal structure of the specificity domain of ribonuclease P. Nature. 421:760–764.[CrossRef][Medline]

44. Rupert, P. B., A. P. Massey, S. T. Sigurdsson, and A. R. Ferre-D'Amare. 2002. Transition state stabilization by a catalytic RNA. Science. 298:1421–1424.[Abstract/Free Full Text]

45. Zhang, L., and J. A. Doudna. 2002. Structural insights into group II intron catalysis and branch-site selection. Science. 295:2084–2088.[Abstract/Free Full Text]

46. Ando, T., T. Tanaka, and Y. Kikuchi. 2003. Substrate shape specificity of E. coli RNase P ribozyme is dependent on the concentration of magnesium ion. J. Biochem. (Tokyo). 133:445–451.[Abstract/Free Full Text]

47. Reiter, N. J., H. Blad, F. Abildgaard, and S. E. Butcher. 2004. Dynamics in the U6 RNA intramolecular stem-loop: a base-flipping conformational change. Biochemistry. 43:13739–13747.[CrossRef][Medline]

48. Hampel, K. J., and J. M. Burke. 2003. Solvent protection of the hammerhead ribozyme in the ground state: evidence for a cation-assisted conformational change leading to catalysis. Biochemistry. 42:4421–4429.[CrossRef][Medline]

49. Shan, S. O., and D. Herschlag. 2002. Dissection of a metal-ion-mediated conformational change in Tetrahymena ribozyme catalysis. RNA. 8:861–872.[Abstract]

50. Pereira, M. J. B., D. A. Harris, D. Rueda, and N. G. Walter. 2002. Reaction pathway of the trans-acting hepatitis {delta}-virus ribozyme: a conformational change accompanies catalysis. Biochemistry. 41:730–740.[CrossRef][Medline]

51. Peracchi, A., L. Beigelman, E. C. Scott, O. C. Uhlenbeck, and D. Herschlag. 1997. Involvement of a specific metal ion in the transition of the hammerhead ribozyme to its catalytic conformation. J. Biol. Chem. 272:26822–26826.[Abstract/Free Full Text]

52. Kim, H. D., G. U. Nienhaus, T. Ha, J. W. Orr, J. R. Williamson, and S. Chu. 2002. Mg2+-dependent conformational change of RNA studied by fluorescence correlation and FRET on immobilized single molecules. Proc. Natl. Acad. Sci. USA. 99:4284–4289.[Abstract/Free Full Text]

53. Pinard, R., K. J. Hampel, J. E. Heckman, D. Lambert, P. A. Chan, F. Major, and J. M. Burke. 2001. Functional involvement of G8 in the hairpin ribozyme cleavage mechanism. EMBO J. 20:6434–6442.[CrossRef][Medline]

54. Peracchi, A., A. Karpeisky, L. Maloney, L. Beigelman, and D. Herschlag. 1998. A core-folding model for catalysis by the hammerhead ribozyme accounts for its extraordinary sensitivity to basic mutations. Biochemistry. 37:14765–14775.[CrossRef][Medline]

55. Zhuang, X. W., H. Kim, M. J. B. Pereira, H. P. Babcock, N. G. Walter, and S. Chu. 2002. Correlating structural dynamics and function in single ribozyme molecules. Science. 296:1473–1476.[Abstract/Free Full Text]

56. Viduna, D., K. Hinsen, and G. Kneller. 2000. Influence of molecular flexibility on DNA radiosensitivity: a simulation study. Phys. Rev. E. 62:3986–3990.[CrossRef]

57. Zacharias, M. 2000. Comparison of molecular dynamics and harmonic mode calculations on RNA. Biopolymers. 54:547–560.[CrossRef][Medline]

58. Brooks, B. R., R. E. Bruccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus. 1983. CHARMM—a program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 4:187–217.[CrossRef]

59. Cornell, W. D., P. Cieplak, C. I. Bayly, I. R. Gould, K. M. Merz, D. M. Ferguson, D. C. Spellmeyer, T. Fox, J. W. Caldwell, and P. A. Kollman. 1996. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. J. Am. Chem. Soc. 117:5179–5197.[CrossRef]

60. Pearlman, D. A., D. A. Case, J. W. Caldwell, W. S. Ross, T. E. Cheatham, S. Debolt, D. Ferguson, G. Seibel, and P. Kollman. 1995. AMBER, a package of computer programs for applying molecular mechanics, normal-mode analysis, molecular dynamics and free-energy calculations to simulate the structural and energetic properties of molecules. Comput. Phys. Commun. 91:1–41.[CrossRef]

61. Foloppe, N., and A. D. MacKerell. 2000. All-atom empirical force field for nucleic acids. I. Parameter optimization based on small molecule and condensed phase macromolecular target data. J. Comput. Chem. 21:86–104.[CrossRef]

62. Aqvist, J. 1990. Ion-water interaction potentials derived from free-energy perturbation simulations. J. Phys. Chem. 94:8021–8024.[CrossRef]

63. Darden, T., D. York, and L. Pedersen. 1993. Particle mesh Ewald—an N·Log(N) method for Ewald sums in large systems. J. Chem. Phys. 98:10089–10092.[CrossRef]

64. Steinbach, P. J., and B. R. Brooks. 1994. New spherical cutoff methods for long-range forces in macromolecular simulation. J. Comput. Chem. 15:667–683.[CrossRef]

65. Jorgensen, W. L., J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein. 1983. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79:926–935.[CrossRef]

66. Grubmüller, H. 1996. SOLVATE 1.0: a program to create atomic solvent models. http://www.mpibpc.mpg.de/groups/grubmueller/solvate.html.

67. Hermann, T., P. Auffinger, and E. Westhof. 1998. Molecular dynamics investigations of hammerhead ribozyme RNA. Eur. Biophys. J. Biophys. Lett. 27:153–165.

68. Pley, H. W., K. M. Flaherty, and D. B. McKay. 1994. Three-dimensional structure of a hammerhead ribozyme. Nature. 372:68–74.[CrossRef][Medline]

69. Jolliffe, I. T. 2002. Principal Component Analysis. Springer-Verlag, New York.

70. Hayward, S., A. Kitao, and N. Go. 1994. Harmonic and anharmonic aspects in the dynamics of BPTI—a normal-mode analysis and principal component analysis. Protein Sci. 3:936–943.[Abstract]

71. Balsera, M. A., W. Wriggers, Y. Oono, and K. Schulten. 1996. Principal component analysis and long-time protein dynamics. J. Phys. Chem. 100:2567–2572.[CrossRef]

72. Becker, O. 2001. Conformational analysis. In Computational Biochemistry and Biophysics. O.M. Becker, A.D. MacKerell Jr., B. Roux, and M. Watanabe, editors. Marcel Dekker, New York.

73. Li, G. H., A. Van Wynsberghe, O. N. A. Demerdash, and Q. Cui. 2005. Normal Mode Analysis: Theory and Applications to Biological and Chemical Systems. Q. Cui and I. Bahar, editors. CRC Press, Boca Raton, FL.

74. Hinsen, K. 1998. Analysis of domain motions by approximate normal mode calculations. Proteins. 33:417–429.[CrossRef][Medline]

75. Blount, K. F., and O. C. Uhlenbeck. 2005. The structure-function dilemma of the hammerhead ribozyme. Annu. Rev. Biophys. Biomol. Struct. 34:415–440.[CrossRef][Medline]

76. Scott, W. G., J. T. Finch, and A. Klug. 1995. The crystal-structure of an all-RNA hammerhead ribozyme—a proposed mechanism for RNA catalytic cleavage. Cell. 81:991–1002.[CrossRef][Medline]

77. Rueda, D., K. Wick, S. E. McDowell, and N. G. Walter. 2003. Diffusely bound Mg2+ ions slightly reorient stems I and II of the hammerhead ribozyme to increase the probability of formation of the catalytic core. Biochemistry. 42:9924–9936.[CrossRef][Medline]

78. Suzumura, K., Y. Takagi, M. Orita, and K. Taira. 2004. NMR-based reappraisal of the coordination of a metal ion at the pro-Rp oxygen of the A9/G10.1 site in a hammerhead ribozyme. J. Am. Chem. Soc. 126:15504–15511.[Medline]

79. Heckman, J. E., D. Lambert, and J. M. Burke. 2005. Photocrosslinking detects a compact, active structure of the hammerhead ribozyme. Biochemistry. 44:4148–4156.[CrossRef][Medline]

80. Van Wynsberghe, A., G. H. Li, and Q. Cui. 2004. Normal-mode analysis suggests protein flexibility modulation throughout RNA polymerase's functional cycle. Biochemistry. 43:13083–13096.[CrossRef][Medline]

81. Tozzini, V. 2005. Coarse-grained models for proteins. Curr. Opin. Struct. Biol. 15:144–150.[CrossRef][Medline]