Benchmarking AMBER Force Fields for RNA: Comparisons to NMR Spectra for Single-Stranded r(GACC) Are Improved by Revised χ Torsions

2.1. NMR

GACC was purchased from Dharmacon Inc. The NMR sample had 3.2 mM GACC in 80 mM NaCl, 10 mM sodium phosphate, and 0.5 mM disodium EDTA at pH 7.0. Two deuterium exchanges were performed on the sample with 99.9% D₂O (Cambridge Isotopes Laboratories). One final deuterium exchange was performed with 99.990% D₂O (Sigma Aldrich).

All spectra were acquired on Varian Inova 500 or 600 MHz NMR spectrometers. The nonexchangeable protons of GACC were assigned from NOESY, ¹H–³¹P HETCOR, DQF-COSY, and TOCSY at 2 and 20 °C (Table S1). NOESY spectra were recorded with mixing times (τ_m) of 200 and 800 ms. Relevant acquisition parameters are shown in Table S2.

2.2. Simulated Annealing Protocol for NMR Structure Modeling

³J_1′2′ spin–spin coupling constants were used to deduce the sugar puckers for each residue (Table S3). Proton–proton distance restraints were calculated from magnitudes of NOEs in 200 ms mixing time NOESY spectra at 2 °C (Table S4) using the two-spin approximation.(57) Restraint limits were calculated by allowing for a ± 4-fold range in each NOE volume to account for relaxation effects, spin diffusion, and/or baseline distortions. The magnitudes of NOEs are proportional to 1/(r_ij)⁶, where r_ij is the distance between protons i and j. H1′ to H2′ NOEs were used as reference NOEs in the distance calculations with H1′ to H2′ distances taken as 2.75 and 3.00 Å, respectively, for C3′-endo and C2′-endo conformations (see Supporting Information for detail).

Cross-peaks in NOESY spectra reveal protons separated by <5 Å. The upper limit of some of the distance restraints was set to 10.0 Å in order to search a broader conformational space even though the magnitudes of these NOE peaks are medium in size indicating a distance <5 Å. Similarly, the lower and upper limits of some distance restraints were set to 4.5 and 99.0 Å, respectively, when no NOE was seen (Table S4).

Simulated annealing MD simulations with NMR restraints were done with sander(58) using the AMBER99 force field. The first MD simulation was initiated from a starting structure of GACC in A-form. Another 199 simulations were performed with the starting structure for each taken from the final structure of the previous simulation. Different random seed numbers were used in each simulation.

The generalized Born implicit solvent model⁵⁹⁻⁶¹ with a 1 M salt concentration was used in the simulated annealing protocol. For each simulation, a total of 50 ps of MD were run with a 1 fs time step. A long-range cutoff of 20.0 Å was used for nonbonded interactions. Berendsen temperature control(62) with a coupling time of 1 ps was used in all simulated annealing MD simulations. Distance and sugar pucker restraints extracted from NOE data were applied using square bottom wells with parabolic sides with force constants of 50.0 kcal/mol Ų and 50.0 kcal/mol rad², respectively (Table S5). Chirality restraints with 50.0 kcal/mol rad² were used to keep chiral atoms in the correct orientations so that the covalent bonds do not flip at high temperatures in the simulated annealing protocol (Table S5). No additional torsional restraints were used. Temperature was increased to 3000 K in 5 ps and gradually cooled to 0 K in 45 ps. During the first 5 ps, relative weights of the NMR restraints were set to 0.1 to make sure that the structure was independent of the starting conformation (see Supporting Information for structural parameters as a function of time during simulation). The NMR restraint weight was set to 1.0 during the rest of the simulation (Table S6).

2.3. Molecular Dynamics Simulations of GACC

A-form single-stranded GACC RNA was created with the nucgen module of AMBER9.(58) GACC was solvated with TIP3P water molecules(63) in a truncated octahedral box with 1076 water molecules. Three Na⁺ ions^58,64 were used to neutralize the system.

The system was minimized in two steps: (1) With GACC held fixed with a restraint force of 500 kcal/mol Ų, steepest descent minimization of 500 steps was followed by a conjugate gradient minimization of 500 steps in order to minimize the potential energy due to water molecules. (2) With all restraints removed, steepest descent minimization of 1000 steps was followed by a conjugate gradient minimization of 1500 steps. The long-range cutoff for nonbonded interactions during the minimization was 10.0 Å.

After minimization, two steps of pressure equilibration were done with the sander module in AMBER9: (1) GACC was held fixed with a restraint force of 10 kcal/mol-Ų. Constant volume dynamics with an atom-based cutoff of 10.0 Å was used for the direct space sum. SHAKE(65) was turned on for bonds involving hydrogen atoms. Temperature was raised from 0 to 300 K in 20 ps. Langevin dynamics with a collision frequency of 1 ps^–1 was used. A total of 20 ps of MD were run with a 2 fs time step. (2) The above conditions were chosen, except that constant pressure dynamics with isotropic position scaling was turned on. Reference pressure was 1 atm with a pressure relaxation time of 2 ps. A total of 100 ps of MD were run with a 2 fs time step. Particle mesh Ewald (PME)^66,67 was used for all calculations.

Each production run was similar to the second step of the pressure equilibration described above. Constant pressure dynamics was performed with a long-range cutoff of 10.0 Å. SHAKE was turned on for bonds involving hydrogen atoms. A total of 50 ns of MD was run at 300 K with a 1 fs time step. In production runs, simulations were carried out with the PMEMD module in AMBER9.(58) Trajectory files were written at each 500 fs time step. The first 1 ns of each MD simulation was omitted from calculations.

Simulations were done on systems using AMBER99(9) and AMBER99χ(44) force fields. For each system and force field, 15 separate 50 ns simulations were run. Different random seed numbers were used for each simulation.

In order to examine the time evolution of GACC with the AMBER99χ force field, a single simulation of 1.9 μs was also run at 300 K. This simulation took over six months on a Dell Precision T5400 workstation using the eight 2.33 MHz Intel Xeon dual quad cores, and was stopped after no new stable structures appeared during the last 800 ns. Root-mean-square deviations of whole structure and backbone from A-form GACC single strand as well as χ and δ dihedral angles of each residue were analyzed. Measured rate constants for stacking and unstacking of nucleotides in AC, CA, CC, A₇, poly A, and poly C range from 5 × 10⁶ to 2 × 10⁷ s^–1,⁵³⁻⁵⁶ indicating that structural transitions occur in roughly 100 ns. Thus, a single 1.9 μs simulation should allow several transitions.

3.1. NMR

³J_1′2′ spin–spin coupling constants from 1D ¹H NMR spectra of GACC revealed that from 2 to 55 °C all the residues prefer C3′-endo sugar pucker, especially the middle nucleotides, A2 and C3 (Tables 1 and S3 and Figure S1). Figure 2 shows the NOESY walk for GACC in an 800 ms mixing time NOESY spectrum. In the 200 ms NOESY spectrum, the NOEs of G1H1′-H8, A2H1′-H8, C3H1′-H6, and C4H1′-H6 are all weak, implying anti base orientation for the residues (Table 1 and Figure S2). NOEs of G1H2′/H3′-A2H8, A2H2′/H3′-C3H6, and C3H2′/H3′-C4H6 are consistent with base stacking in A-form structures (Figure 3). Weak C3H6-C4H2′ and C3H3′-C4H2′ cross-peaks are seen in 200 and 800 ms NOESY spectra (Figures 3 and S3), however, which suggests a population of conformations with a flipped, upside down C4 sugar. Such reversal of ribose has been observed previously⁶⁸⁻⁷¹ but is rare. Thus, no C3H6-C4H2′ and C3H3′-C4H2′ distance restraints were used in the simulated annealing protocol. The chemical shifts of all phosphorus resonances were within 0.2 ppm indicating no unusual backbone conformations (Table S1).

3.2. Dynamics of A-Form GACC in MD Simulations with the AMBER99 and AMBER99χ Force Fields

For each of the 50 ns unrestrained MD simulations, snapshots were extracted from the trajectories at intervals of 5 ns and analyzed in order to cluster them into structural groups. No preferred conformation was found for the fifteen 50 ns AMBER99 MD simulations, and G1 of GACC went to a syn base orientation within 0.5 ns on average. Coupled with the force field’s preference for C2′-endo sugar pucker,(44) this creates distortions in the backbone and the other residues. As a result, random-coil conformations are observed within less than 10 ns with the AMBER99 force field.

While simulations with AMBER99 generated rapid and large structural fluctuations of GACC within 50 ns, simulations with AMBER99χ did not. For the fifteen 50 ns AMBER99χ MD simulations, 107 and 11 out of 150 structures cluster into A-form-like major and minor conformations, respectively. To further explore the time dependence of dynamics with AMBER99χ, a single 1.9 μs MD simulation was run. Time evolutions of rms deviations of all atoms and backbone heavy atoms with respect to a pure A-form GACC conformation and of the χ and δ torsion angles of each residue are shown in Figure 4. This decomposition provides insight into the individual conformational transformations of the structural dynamics of GACC. Cluster analysis was done on this single 1.9 μs MD simulation and revealed that 63 and 18 out of 380 structures represent major and minor conformations, respectively, while 77 structures represent a stable non-A-form conformation (see Figures 4 and S4), which is not consistent with NMR spectra. The latter structure appears after 700 ns. In the first 700 ns, 68% of the structures are A-form-like (Figure 4). The rest of the structures in the 1.9 μs MD simulation do not cluster to any conformation and are in non-A-form conformations. The results show that for 700 ns AMBER99χ generates structures consistent with NMR but that additional revisions are necessary. In particular, GACC did not return to an A-form-like structure during the 800 ns between 1.1 and 1.9 μs when the predicted structures were random coil.

3.3. Predictions of Sugar Pucker and Base Orientation of Individual Residues of GACC

Structures extracted every 0.5 ps from the unrestrained MD trajectories were analyzed to quantify the sugar pucker and base orientation of GACC residues predicted by each force field. Population distribution plots of pseudorotation phase angles and χ torsions for each residue were created (Figures S5 and S6) to obtain the percentages of C3′-endo sugar pucker and anti type base orientation (Table 1). The sugar pucker of a five-membered ring can be described by the pseudorotation phase angle.^72,73 In RNA and DNA, pseudorotation phase angles of 18° and 162° that correspond to C3′-endo and C2′-endo sugar pucker, respectively, are mainly found.

These percentages can be quantitatively compared to NMR results (Table 1). At 25 °C, NMR indicates around 65%, 85%, 85%, and 75% C3′-endo sugar pucker for G1, A2, C3, and C4, respectively. For a series of fifteen 50 ns simulations that start in A-form, the AMBER99 force field predicts 49%, 33%, 64%, and 41% C3′-endo sugar pucker for G1, A2, C3, and C4, respectively (Table 1 and Figure S5). Evidently, AMBER99 underestimates the favorability of the C3′-endo sugar pucker for GACC. The 1.9 μs AMBER99χ simulation predicts 76%, 81%, 83%, and 67% C3′-endo sugar pucker for G1, A2, C3, and C4, respectively, which agrees reasonably well with the NMR data (Table 1).

The predicted percentages of C3′-endo sugar pucker, however, depend on the time domain within the 1.9 μs simulation (Table 1). For example, they are 90%, 85%, 92%, and 74% for G1, A2, C3, and C4, respectively, in the first 700 ns, where GACC samples A-form-like and random coil conformations (Figure 4). In contrast, they are 52%, 69%, 67%, and 65% in the 1101–1900 ns region where MD predicts random coil.

Quantitative comparisons can also be made between predicted and NMR determined orientation of bases relative to sugars. NMR reveals anti conformations for all the residues of GACC. For a series of fifteen 50 ns simulations that start in A-form, the AMBER99 force field predicts GACC to have 37%, 97%, 90%, and 92% anti/high-anti conformations for G1, A2, C3, and C4, respectively. The anti range of the χ torsions observed in AMBER99 simulations was 180–360° and 180–300° for G1/A2 and C3/C4, respectively (see Figure S6). The broad range for AMBER99 is due to generation of both anti and high-anti conformations. The 1.9 μs AMBER99χ simulation predicts 53%, 57%, 93%, and 90%, anti conformations for G1, A2, C3, and C4, respectively, with only 2–6% of high-anti conformations for these residues. Most of the syn conformations, however, occur after 1.1 μs when GACC is not in a stable structure (Table 1 and Figure 4). Prior to 1.1 μs, G1 and A2 are 80% anti on average, with the rest either syn or high-anti.

3.4. Predictions of Structures Satisfying H–H NOE Distance Restraints

For each force field, the structures used in Part 3.3 were also used to calculate the percentages of structures consistent with the NMR distance restraints for some of the H–H NOEs (Table 1). These predictions were extracted from the population distribution plots shown in Figures S7–S9 and provide another measure of prediction accuracy.

In 50 ns simulations starting with A-form and using the AMBER99 force field, the percentages of structures consistent with NMR distance restraints for G1H2′-A2H8, A2H2′-C3H6, A2H3′-C3H6, A2H2-C3H1′, C3H2′-C4H6, and C3H3′-C4H6 are 19%, 24%, 49%, 9%, 41%, and 65%, respectively. With the 1.9 μs AMBER99χ simulation, the fractions are 31%, 54%, 60%, 63%, 35%, and 72%, respectively (Table 1). The NOE cross-peak A2H2-C3H1′ represents stacking of the middle bases A2 and C3, and this cross-peak is underestimated much more by AMBER99 than by the AMBER99χ force field, even when the AMBER99χ simulation is run to 1.9 μs. As expected, the percentages of structures in the 1.9 μs simulation that are consistent with NMR distance restraints depend on the time domain (Table 1). For example, between 701 and 1100 ns, the G1H2′-A2H8 and C3H2′-C4H6 distances never fall in the range consistent with NMR data. As discussed below, the stable structure generated between 701 and 1100 ns is inconsistent with NMR spectra.

3.5. NOE “Distances” Measured by NMR and Predicted by AMBER99/AMBER99χ Force Fields

Comparison of NMR measured and MD predicted NOE “distances” (R) is a more direct test of the force fields. These “distances” are not true distances if GACC populates more than one conformation. The structures used in Part 3.3 were used to calculate the NOEs (1/r_ij⁶) for each snapshot, which then yielded average NOE and, therefore, average “distance”, R. That is, the “distance” is an ensemble average with each exact distance weighted by its inverse sixth power. For the combined fifteen 50 ns MD simulations, rms deviations of six NOE “distances” between NMR and, respectively, AMBER99 and AMBER99χ are 0.84 and 0.39 Å. For the single 1.9 μs MD simulation of AMBER99χ, this rms deviation is 0.24 Å (Tables 2 and S7). Thus, AMBER99χ performs better with this benchmark.

Again as expected, the distance comparison for the 1.9 μs simulation depends on the time domain (Table 2). The distances predicted from the first 700 ns agree better with the NMR than those predicted from the 701 to 1100 ns and 1101 to 1900 ns domains.

3.6. Conformational Fluctuations in the Single 1.9 μs MD Simulation with AMBER99χ Force Field

Figure 4 shows the time evolution for the 1.9 μs MD simulation. In Figure 4, regions highlighted with red and green represent sampling of major and minor A-form-like conformations, respectively, whereas regions highlighted with blue represent a stable non-A-form conformation (Figure S4 and Table S8). The rest of the regions highlighted with black represent sampling of random-coil non-A-form conformations.

During the first 700 ns MD, there are multiple major to minor and A-form-like to non-A-form conformational transitions (Figure 4). The switch from major to minor conformation is due to the conformational change of the backbone (Figure 4). During these periods, χ dihedrals of G1, A2, C3, and C4 are mostly around 180°, which is the anti region. G1 goes to the syn conformation a few times during these periods, but returns to anti conformation almost instantly (Figure 4). Between 200 and 375 ns, 550–600 ns, and around 520 ns, however, the structure resembles neither major nor minor conformations. At the beginning of these times, base orientation of G1 goes to either high-anti or syn, which distorts the backbone. This allows A2 to change its base orientation to either high-anti or syn (Figure 4). During these times, C3 and C4 mostly stay in anti base orientation, but G1, A2, and C3 start sampling C2′-endo sugar puckering (Figure 4). During the whole simulation time, C4 rapidly fluctuates between C3′-endo and C2′-endo sugar puckers (Figure 4).

Between 701 and 1100 ns, a stable and unique non-A-form conformation is sampled by AMBER99χ (Figures 4 and S4). At around 700 ns, the base orientation and sugar puckering of all four residues are anti and C3′-endo, respectively, but the backbone changes dramatically and creates a conformation with G1 between C3 and C4 (Figure S4). This stable conformation has unique potential NOE cross-peaks such as G1H8-C3/C4H2′/H3′, G1H8-C3/C4H5/H6, and C3H6-G1H2′/H3′ due to proton–proton distances <5 Ǻ. These potential cross-peaks would be easy to detect because they would appear in unique and uncrowded NOESY regions but no such cross-peaks have been detected.

After 1100 ns, the structure transforms to random coil conformations. Again, first G1 goes to a syn base orientation that gives enough flexibility to the backbone to search other possible configurations. This allows A2 to go to a syn base orientation. C3 and C4 rarely change their base orientations but sample syn a couple of times (Figure 4). Similar to the random-coil conformations seen within 700 ns, fluctuations of the δ torsions are high after 1100 ns where all the residues start to sample C2′-endo sugar puckering significantly (Figure 4).

The terminal G1 and C4 nucleotides are flexible enough to change their conformations from anti to syn and/or C3′-endo to C2′-endo. These changes are seen in Figure 4, but are not the main reason that the structure samples random-coil conformations; this broader sampling is due primarily to conformational changes of A2. C3 and C4, however, are found primarily in anti conformations as expected for cytidine,(44) even when the structure is in a random-coil conformation. One aspect of Figure 4 is that in order for the single-stranded GACC to have a stable conformation, it has to have all the bases in anti orientations (Figure 4). Otherwise, high fluctuations on the χ and δ torsions yield random-coil conformations.

3.7. Simulated Annealing of GACC using NMR Restraints

GACC presumably has an ensemble of structures in rapid exchange that contribute to the NMR spectrum. Simulated annealing with NMR restraints was used to see if any particular structures are consistent with the NMR data. The protocol described in the Methods section generated 182 out of 200 structures with constraint energies less than 0.1 kcal/mol. Of these structures, two had energies greater than the mean plus the standard deviation. As a result, 180 out of 200 structures were analyzed. Cluster analyses showed that 116 and 32 represent A-form-like major and first minor conformations, respectively (Figures 5 and S4 and Tables S9 and S10). Another 18 structures represent a second A-form-like minor conformation with an average energy less favorable than both the major and first minor conformations (Figure S4). Both minor conformations have a flipped, upside down C4 sugar, but with different backbone orientations (Figure S4). The flipped C4 sugar conformation is consistent with weak C3H6–C4H2′ and C3H3′-C4H2′ NOEs that were observed, but not used as restraints. The major and first minor conformations have stacking between G1/A2, A2/C3, and C3/C4, while the C3/C4 stacking is missing in the second minor conformation. With the exception of the flipped C4, the structures generated by simulated annealing are largely A-form-like.

AMBER99χ Torsion Parameters Improve Agreement between NMR and Simulations

The NMR results for GACC can be qualitatively compared to predictions of the AMBER99 and AMBER99χ force fields. In the AMBER99 force field unrestrained simulations, there are large fluctuations within 50 ns in rms deviations of individual residues from anti/C3′-endo conformation and in backbone from A-form (Table 1). Evidently, anti/C3′-endo conformations and A-form are not preferred by the AMBER99 force field and transition barriers allow rapid conformational changes (see Figures S10 and S11). In contrast, all residues of GACC in the fifteen 50 ns AMBER99χ force field simulations show stable conformations similar to anti/C3′-endo. Because few transitions are seen within 50 ns with AMBER99χ, the 1.9 μs simulation provides a better comparison. In that simulation, a stable backbone similar to A-form is retained for almost 68% of the first 700 ns (Figure 4). After 700 ns, various random-coil and a non-A-form stable structure inconsistent with NMR spectra are seen.

The non-A-form conformation is likely an artifact of the force field, and suggests that reparameterization beyond χ is necessary. The high-dimensional character of the conformational space, however, does not allow definitive identification of the torsions causing this non-A-form conformation to be stable enough to be sampled significantly by AMBER99χ. Histogram analyses of RNA backbones in crystal structures show preferred regions for α, β, γ, δ, ε, and ζ,^68,69 and most dihedrals in the 1.9 μs MD occur in these regions (Figure S12). All the residues of the stable non-A-form conformation have ζ between 7° and 71° (Table S8), whereas an A-form RNA has ζ around −70° and in crystal structures, ζ is rarely between −40° and +40°.(68) It is possible that a revision of the ζ torsion for RNA might improve predictions. Alternatively, the stable non-A-form and the random-coil conformations may actually exist in solution, but are not prevalent enough to influence the NMR spectra. Testing this hypothesis, however, would require running much longer MD simulations to allow more switches in conformation.

Quantitative comparisons show that the sugar pucker, base orientation, and proton–proton NOE distances derived from NMR results are predicted better by AMBER99χ than AMBER99 (Tables 1 and 2). Even weak C3H6- C4H2′ and C3H3′- C4H2′ NOEs are consistent with the flipped C4 generated by unrestrained AMBER99χ simulations.

In order to test the effects of reparameterizing only χ, the parmbsc0 force field,(46) which has revised parameters for α and γ, was not used in MD simulations. Preliminary results on a different benchmark, however, indicate that combining parmbsc0 with AMBER99χ further improves agreement with experiment (I.Y., H.A.S., D.H.T., unpublished). A different reparameterization of χ produced the Ode force field.(43) Although the Ode force field has more parameters than AMBER99χ, it predicts similar potential energy surfaces for nucleosides.(44) Thus, it should provide predictions similar to those of AMBER99χ, although tests on UUCG and GNRA tetraloop hairpins indicate that the Ode force field overpredicts high-anti conformations while AMBER99χ does not.(45) Additional tests are required to see if this difference is general.

AMBER99χ Torsion Parameters Improve Agreement between Predicted and Expected Dynamics

The structural fluctuations observed with AMBER99 occur in less than 10 ns, which is much faster than the ∼100 ns expected on the basis of experimental measurements on single-stranded oligo and polynucleotides.⁵³⁻⁵⁶ In contrast, AMBER99χ generated structures consistent with NMR spectra are typically stable for 50–100 ns (Figure 4). Evidently, AMBER99χ reasonably approximates barriers to conformational change.

Major and Minor Conformations Generated in the Initial 700 ns of Unrestrained AMBER99χ MD Are Similar to Structures Generated by Simulated Annealing with AMBER99 Restrained with NMR Data

Because the dynamics of unpaired oligonucleotides occur in less than 1 μs, NMR reports a population weighted average of spectral properties for the ensemble of structures. Somewhat surprisingly, the NMR spectra for GACC are consistent with expectations for a high population of A-form-like conformations. Simulated annealing with AMBER99 and NMR restraints generates two A-form-like structures with the minor structure having an upside down terminal ribose that is consistent with minor NOESY peaks. This suggests that the ensemble has largely A-form conformations, rather than a high population of random conformations. The first 700 ns of AMBER99χ unrestrained MD is largely consistent with this interpretation (Figure 4). After 700 ns, however, a structure inconsistent with NMR spectra is generated, suggesting that further refinement of the force field will be important. Experimental NMR restraints, however, appear able to compensate for limitations of the current AMBER force field.

Summary

In order to sample physical regions of conformational space in MD simulations, the force field must accurately describe the interactions that give rise to conformational preferences. Such sampling is important for predicting the structure and dynamics of RNA and the thermodynamics of folding an RNA from an unpaired single strand. Modification of χ parameters in the AMBER99 force field improved agreement between simulations and NMR spectra of single-stranded GACC RNA, as summarized in Figure 6. The modified χ torsion parameters are therefore expected to improve sampling of non-Watson–Crick regions and terminal base pairs in RNA simulations. In Watson–Crick regions of an RNA duplex, the χ torsions are restricted by cross-strand hydrogen bonding, so the AMBER99χ torsion parameters should have less effect. A stable non-A-form and random-coil conformations are observed significantly during the 1.9 μs MD simulation, suggesting refinement of other parameters will be important. In principle, such refinement can rely on QM calculations and the parmbsc0 revisions(46) will likely help. Many subtle interactions, however, define RNA.(42) Thus, some combination of knowledge based force fields(74) with QM may prove useful. Finally, NMR spectra for single-stranded regions of RNA provide useful benchmarks for testing any revisions.