|Home | About | Journals | Submit | Contact Us | Français|
A reparameterization of the torsional parameters for the glycosidic dihedral angle, χ, for the AMBER99 force field in RNA nucleosides is used to provide a modified force field, AMBER99χ. Molecular dynamics simulations of cytidine, uridine, adenosine, and guanosine in aqueous solution using the AMBER99 and AMBER99χ force fields are compared with NMR results. For each nucleoside and force field, 10 individual molecular dynamics simulations of 30 ns each were run. For cytidine with AMBER99χ force field, each molecular dynamics simulation time was extended to 120 ns for convergence purposes. Nuclear magnetic resonance (NMR) spectroscopy, including one-dimensional (1D) 1H, steady-state 1D 1H nuclear Overhauser effect (NOE), and transient 1D 1H NOE, was used to determine the sugar puckering and preferred base orientation with respect to the ribose of cytidine and uridine. The AMBER99 force field overestimates the population of syn conformations of the base orientation and of C2′-endo sugar puckering of the pyrimidines, while the AMBER99χ force field’s predictions are more consistent with NMR results. Moreover, the AMBER99 force field prefers high anti conformations with glycosidic dihedral angles around 310° for the base orientation of purines. The AMBER99χ force field prefers anti conformations around 185°, which is more consistent with the quantum mechanical calculations and known 3D structures of folded ribonucleic acids (RNAs). Evidently, the AMBER99χ force field predicts the structural characteristics of ribonucleosides better than the AMBER99 force field and should improve structural and thermodynamic predictions of RNA structures.
Understanding the physical interactions governing the structure and dynamics of ribonucleosides should improve the accuracy of simulations of ribonucleic acid (RNA) molecules. Methods for simulating biological systems include residue-centered force fields (coarse-grained),(1) atom-centered force fields (AMBER,(2) CHARMM,3,4 GROMOS),5,6 approximate quantum mechanics,7,8 and mixed quantum mechanics/molecular mechanics methods (QM/MM).9−18 With advances in computer power, it is possible to run simulations at least as long as milliseconds and microseconds with coarse-grained and atom-centered potentials, respectively.19−23 The AMBER force fields are particularly widely used for simulations of RNA. They have provided satisfactory descriptions of structural and thermodynamic properties for some RNA and deoxyribonucleic acid (DNA) systems,24−29 while some challenging systems still provide difficulty.30−32 Predictions for the individual ribonucleosides have not been extensively used as benchmarks for AMBER force fields. A fundamental understanding of nucleosides is crucial to simulate the behavior of residues in single strands, noncanonical base pairs, and hairpins. Mimicking the real behavior of ribonucleosides in simulations should improve predictions of RNA properties.
Due to limitations of computer power, small model systems were used to parametrize the glycosidic dihedral angle in the AMBER94 force field.(2) In this article, the glycosidic dihedral angle, χ, of ribonucleic acids is reparameterized by extending the quantum mechanical (QM) fitting protocol, and new parameters are used in a revised force field, AMBER99χ. Structural and thermodynamic results are extracted from molecular dynamics (MD) simulations using AMBER99(33) and AMBER99χ force fields.
Previous experimental work on nucleosides and nucleotides has classified the behavior of individual torsion angles.34−41 Structures of modified and unmodified nucleosides/nucleotides have been interrogated by one-dimensional (1D) 1H nuclear magnetic resonance (NMR) and steady-state 1D 1H nuclear Overhauser effect (NOE) difference spectroscopy (SSNOE).42−50 In this work, transient 1D 1H NOE spectroscopy(51) and sugar proton coupling constants extracted from 1D 1H NMR spectra for cytidine (C) and uridine (U) are used to quantitatively deduce the preferred conformations of the glycosidic dihedral angle and the sugar pucker, respectively. These results are compared to computational predictions. The AMBER99 force field overestimates the fraction of syn conformations for the base orientation and of C2′-endo sugar puckering of the pyrimidines, while the results of AMBER99χ are more consistent with that of the experimental NMR data. Simulations on adenosine (A) and guanosine (G) show that AMBER99 prefers high anti conformations around 310°, while AMBER99χ prefers anti conformations around 185°. The latter is more consistent with QM energy profiles and is the typical anti region seen in crystal structures of nucleic acids.
C, U, A, and G were purchased from Sigma Aldrich. Solutions of 0.2, 1, and 5 mM nucleosides were made in H2O with an NMR buffer consisting of 80 mM NaCl, 10 mM sodium phosphate, and 0.5 mM disodium EDTA at a pH of 7.0. Two lyophilizations were performed on each sample, reconstituting each time with 99.9% D2O (Cambridge Isotopes Laboratories). One final lyophilization was performed, and each sample was reconstituted with 99.990% D2O (Sigma Aldrich).
NMR experiments were performed with Varian Unity Inova 500 and 600 MHz spectrometers. Chemical shift data were extracted from 1D 1H NMR (see Supporting Information). For A, the chemical shifts of H8, H2, H1′, and H2′ protons vary with concentration, implying that there is base stacking and/or base pairing interactions (see Supporting Information).(52) The 5′-guanosine monophosphate is known to form quadruplex structures and other kinds of aggregates in solution,53−55 and presumably, guanosine does the same. Aggregation and even precipitation was seen in 5 mM G solutions. Thus, the NMR spectra for nucleosides of A and G were not interpreted, except that 3J spin−spin couplings of 0.2 mM samples were measured as a function of temperature (see Supporting Information).
For C and U, transient 1D NOE measurements were performed with a selective inversion−recovery experiment in which the frequency of the selective inversion pulse was alternated between on resonance with the H6 proton and 2000 Hz downfield, where no resonances are present. The on/off resonance spectra were subtracted, and the integral of the resulting NOE peaks was divided by peak integrals in a 1D spectrum to obtain percent enhancement. Steady-state 1D NOE spectra were acquired in a similar manner with the inversion−recovery replaced by low-power irradiation for 10 s that was on/off the H6 resonance.
Initial geometries were chosen to represent experimental conformations. The γ dihedral angle (O5′−C5′−C4′−C3′) was set to 54°, which is the observed γ value for A-form RNA. The δ dihedral angle (C5′−C4′−C3′−O3′) was set to either 140° or 81°, which is C2′- or C3′-endo sugar pucker, respectively. The O4′−C1′−C2′−C3′ dihedral was set to either 32° or −24° to force the sugar pucker to stay in C2′- or C3′-endo conformations, respectively. In ribonucleosides, there are three OH groups (5′, 3′, and 2′) that are free to rotate in solution. The 3′ OH group will not interact with the base as much as the 5′ and 2′ OH groups. Thus, different conformations of 5′ and 2′ OH groups were included in the fitting.
For each nucleoside (Figure (Figure1),1), four different sugar conformations (Table (Table1)1) were chosen for QM calculations with Gaussian03.(56) For each sugar conformation, a potential energy surface (PES) scan was done around the glycosidic dihedral angle with increments of 5°, yielding 4 × 72 = 288 conformations for each nucleoside. For each conformation in the PES scan, the structures were first optimized with HF/6-31G* level of theory. During the optimization, most dihedrals were frozen in order to have a smooth energy profile with respect to the χ torsion angle (see Supporting Information). Then, QM energies, EQM, were calculated with MP2/6-31G* level of theory.
Atom notations of nucleosides: (a) cytidine, (b) uridine, (c) adenosine, and (d) guanosine. For C and U, χ is the dihedral angle defined by O4′−C1′−N1−C2, and for A and G, χ is defined by O4′−C1′−N9−C4. ...
Dihedral Angles Used to Create the Four Sugar Conformations (sc) for Each Nucleoside
The molecular mechanics (MM) energies, EMM(noCHI), of each conformation were calculated by restraining the dihedral angles to the values of the optimized QM geometries with a force constant of 1500 kcal/mol·A2 using the AMBER99(33) force field parameters, except χ torsion parameters were set to zero (see Supporting Information). AMBER9(57) was used to calculate the MM energies, which use the default 1−4 vdW and electrostatic screening factors of 2.0 and 1.2, respectively.
The energy difference, EQM − EMM(noCHI), represents the potential energy due to χ torsion:
Here, ϕ1 and ϕ2 are the dihedral angles of O4′−C1′−N1−C6 (O4′−C1′−N9−C8) and C2′−C1′−N1−C6 (C2′−C1′−N9−C8), respectively. Vn1 and Vn2 are the potential energy barriers of O4′−C1′−N1−C6 (O4′−C1′−N9−C8) and C2′−C1′−N1−C6 (C2′−C1′−N9−C8) torsions. For each nucleoside, a separate fitting was done to calculate the χ torsion energy barriers, Vn1 and Vn2. The new χ torsion parameters are listed in Table Table22.
New χ Torsion Parameters for Adenosine, Guanosine, Cytidine, and Uridine
Each structure was created with the xleap module of AMBER9.(57) Two conformations were used as initial structures: C3′-endo sugar puckering with base orientations of anti or syn. C, U, A, and G were solvated with TIP3P water molecules(58) in a truncated octahedral box, having 458, 451, 427, and 430 water molecules, respectively.
The structures were minimized in two steps: (i) With the nucleoside held fixed with a restraint force of 500 kcal/mol·Å2, steepest descent minimization of 500 steps was followed by a conjugate gradient minimization of 500 steps. (ii) 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 8 Å.
After minimization, two steps of pressure equilibration were done with the SANDER module in AMBER9: (i) Nucleosides were held fixed with a restraint force of 10 kcal/mol·Å2. Constant volume dynamics with a long-range cutoff of 8 Å was used. SHAKE(59) was turned on for bonds involving hydrogen atoms. The 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. (ii) The above conditions were chosen, except the constant pressure dynamics with isotropic position scaling was turned on. The 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. The particle mesh Ewald (PME) method was used for all simulations.
The production run was similar to the second step of the pressure equilibration described above. Constant pressure dynamics was chosen with a long-range cutoff of 8 Å. SHAKE was turned on for bonds involving hydrogen atoms. For each nucleoside, a total of 30 ns of MD were run with a 1 fs time step. For cytidine with AMBER99χ force field, the simulation time was 120 ns for convergence purposes. In production runs, simulations were carried out with the PMEMD module in AMBER9.(57) Trajectory files were written at each 250 fs time step.
Simulations were performed for systems prepared with the AMBER99 and AMBER99χ force fields. For C, U, A, and G, and each force field, 10 separate simulations of 30 ns each were run at 300 K yielding a total of 300 ns of explicit solvent MD simulation (see Supporting Information). Five of the 10 MD simulations had a starting structure of anti type, while the other five had a starting structure of syn type (see Supporting Information). For C with AMBER99χ force field, the simulations were extended to 11 separate simulations with 120 ns each (see Supporting Information). The fractions of anti and syn conformations observed were essentially independent of the starting structure as were values obtained for C when the time for each of the 11 simulations was extended from 30 to 60 ns and 120 ns (see Supporting Information).
Ultrasonic relaxation studies in aqueous solution revealed a relaxation time of 3 ns for A and no relaxation signal for pyrimidines.60,61 The relaxation signal is attributed to the syn→anti transformation of the χ torsion. Evidently, 300 ns of MD simulations of the nucleosides is sufficient to sample adequately the syn→anti transformation.
In solution, nucleosides have two important regions that describe their structures: (i) the glycosidic dihedral angle, and (ii) the sugar pucker. NMR NOE experiments were done to analyze the structures of C and U.
The magnitudes of NOEs are proportional to 1/(rij)6, where rij is the distance between the protons of i and j. When the base of a pyrimidine is oriented in an anti conformation, the H6 proton is about 3.5 Å from the H1′ proton, essentially independent of sugar pucker.(62) Thus, irradiation of H6 yields a moderate NOE to H1′. When the base of a pyrimidine is oriented in a syn conformation, however, the H6 proton is about 2.1 Å from H1′, yielding a strong NOE to H1′ when H6 is irradiated.(62) In pyrimidines, the distance between the H5 and H6 protons is constant at 2.48 Å, which can be used as a reference for calculating interproton distances from NOESY or transient NOE experiments according to eq 3:(63)
Here, NOEij is the NOE between protons i and j, NOEH5H6 is the NOE between H5 and H6 protons, and rH5H6 is the distance between the H5 and H6 protons, i.e. 2.48 Å.
Transient NOE spectroscopy(51) with different mixing times was used to quantitatively analyze the preferences for anti/syn populations, and the results are presented in Table Table33 (also see Supporting Information). Transient NOE is similar to NOESY NMR except that it is 1D. To minimize spin diffusion effects and maximize signal-to-noise ratio, mixing times in the linear region of intensity vs mixing time plots were used to estimate distances between protons (see Supporting Information). A two-state model described by the following equation, which assumes that the structure is in either syn or anti conformations, was used to determine the proportions of each conformation:
Experimentally Deduced and Force Field Predicted Base Orientation and Sugar Puckering for C, U, A, and G, and ΔG° (in kcal/mol) of Syn→Anti and C2′-endo→C3′-endo Transformations for C and Ua
Here, NOEH1′H6 is the NOE between the protons of H1′ and H6, Fanti and Fsyn are the fractions of anti and syn conformations satisfying Fanti + Fsyn = 1, rH1′H6,anti and rH1′H6,syn are the distances between the protons of H1′ and H6 when the structures are in anti and syn conformations, respectively, which are 3.48 Å and 2.12 Å, corresponding to the distances extracted from the minimum energy structures of the PES scans for C and U (see Methods Section). As can be seen from Table Table3,3, the anti orientation is favored over syn. Comparison of NMR results for C at 2 and 10 °C show that the fraction of anti base orientation is essentially independent of temperature (Supporting Information). Higher temperature could not be used because of the overlap of the H1′ and H5 peaks (see Supporting Information). SSNOE spectroscopy confirms that anti is favored over syn base orientation (see Supporting Information).
Sugar proton coupling constants extracted from 1D 1H NMR spectra were used to determine the sugar puckering on the basis of the following equation:(64)
where 3J1′2′ and 3J3′4′ are 3J spin−spin couplings between H1′ and H2′ and between H3′ and H4′ protons, respectively. The proportion of C2′-endo sugar puckering is equal to (1 − fraction of C3′-endo). Sugar pucker (±2%) is independent of temperature from 5 to 40 °C (Supporting Information), and results at 30 °C are presented in Table Table33.
Figures Figures22−5 show the QM, MMAMBER99, MMAMBER99χ, and MMOde(65) energy profiles with respect to the glycosidic dihedral angle of all the structures used in the fitting protocol for the nucleosides, where AMBER99, AMBER99χ, and Ode(65) force fields were used to calculate MMAMBER99, MMAMBER99χ, and MMOde energies, respectively. In all the plots, energy profiles of the AMBER99χ force field describe the QM energy profiles best, although the Ode force field’s energy profile is also similar to the QM energy profiles. The differences between the predictions of the AMBER99χ and Ode force fields is likely due to the Ode force field using CH3, H2C−CH3 and H2C−O−CH3 as model systems to represent the sugar, while the AMBER99χ force field used the entire ribose with four different sugar conformations to calculate the χ torsional parameters. The Ode force field also uses more parameters. Yet, both Ode and AMBER99χ force fields should provide similar predictions for structural and/or dynamical properties of RNA. Comparisons of the force fields to QM calculations on eight sugar conformations not included in the fitting showed that AMBER99χ also describes those QM energy profiles better than AMBER99 and Ode force fields (see Supporting Information).
For comparison with NMR results, predictions of population distributions of χ dihedral angle and sugar pucker were analyzed for C, U, A, and G using the combined trajectories of the 10 individual MD simulations with AMBER99 and AMBER99χ force fields (see Methods). Population distribution plots in 2D of χ dihedral and pseudorotation angles for each nucleoside are shown in Figures Figures66 and and7.7. Table Table33 shows the force field predictions of base orientation and sugar pucker for each nucleoside (also see Table Table44 and Supporting Information). Analyses of the individual MD simulations show at least seven syn↔anti transformations for each (see Supporting Information).
Population Analysis Results for C, U, A and G of the AMBER99 and AMBER99χ Force Fieldsa
Table Table33 shows the experimental results for C and U as well as the predictions of AMBER99 and AMBER99χ force fields of the base orientation and the sugar pucker for C, U, A, and G. For the syn→anti equilibrium of C and U, NMR indicates 87% and 93% anti conformation, respectively, corresponding to ΔG°syn→anti of −1.07 and −1.45 kcal/mol. The AMBER99 force field predicts 30% and 28% anti conformation, respectively, corresponding to ΔG°syn→anti of 0.49 and 0.55 kcal/mol. In comparison, the AMBER99χ force field predicts 66% and 83% anti conformation, respectively, corresponding to ΔG°syn→anti of −0.45 and −0.95 kcal/mol, closer to the NMR results. Evidently, AMBER99 overestimates the syn conformations of C and U (see Figure Figure66).
For the C2′-endo→C3′-endo equilibrium of C and U, NMR indicates 60% and 56% C3′-endo sugar puckering at 30 °C, respectively, corresponding to free energy differences, ΔG°C2′→C3′, of −0.24 and −0.15 kcal/mol (Table (Table3).3). The percentages are essentially independent of temperature from 5 to 40 °C (see Supporting Information). The AMBER99 force field predicts 27% and 35% C3′-endo sugar pucker at 27 °C, respectively, corresponding to ΔG°C2′→C3′ of 0.58 and 0.36 kcal/mol. In comparison, the AMBER99χ force field predicts 54% and 55% C3′-endo sugar pucker at 27 °C, respectively, corresponding to ΔG°C2′→C3′ of −0.11 and −0.13 kcal/mol, which is close to the experimental values. Evidently, AMBER99 underestimates C3′-endo sugar puckering of C and U (see Figure Figure66).
The AMBER99χ force field predicts A and G to have 13% and 24% anti conformation (Table (Table3),3), respectively, with a χ dihedral angle around 185°, which is consistent with QM calculations and typical of the anti region seen in crystal structures of RNA.(66) The AMBER99 force field predicts 15% and 11% anti conformation (Table (Table3),3), respectively, but with a χ dihedral angle around 310° (Figure (Figure7),7), which is the high anti region. QM PES scans did not find any minimum around 310° but rather between 180−250° for three different sugar puckers for A and G (Figures (Figures44−5 and Supporting Information, where the x-axis, however, is χ + 180°).
The concentration dependence of chemical shifts for A and G indicated aggregation at concentrations required to determine NOEs with enough signal-to-noise to determine the base orientation quantitatively. Pioneering studies of 2′- and 3′- AMP and GMP at high concentrations, however, indicated syn populations well over 50%.67,68
The AMBER99 force field predicts A and G to have 24% and 35% C3′-endo sugar puckering, respectively, while AMBER99χ predicts 32% and 54%. Chemical shift data of 0.2, 1.0, and 5.0 mM A implies base stacking that differs with concentration. The differences of the chemical shifts between 0.2 and 1.0 mM samples are small, however. Therefore, the 0.2 mM samples of A and G were used to calculate 3J spin−spin couplings to estimate the sugar puckering (see Supporting Information). At room temperature, the C3′-endo sugar puckering of A and G is about 40% (Table (Table33 and Supporting Information). For A, both force fields’ predictions are similar to the experimental results. For G, the AMBER99 force field apparently predicts better than AMBER99χ does. It is known, however, that guanosine monophosphate forms quadruplex structures and other aggregates in solution.53−55 Aggregation and precipitation were seen by eye in the 5 mM G NMR samples. Thus, it is not conclusive whether 0.2 mM G can be used to reveal the sugar puckering of monomer G.
There are several reasons why the AMBER99χ force field improves predictions for nucleosides. When the χ torsions were parametrized for AMBER99, model systems for adenosine and thymidine were used, and the results were generalized for all DNA/RNA residues.(2) Moreover, the model systems mimicked deoxyribose C2′-endo sugar puckering. At that time, QM calculations were limited by computer power and only 8−9 data points were used in the QM fitting. Also, in the AMBER99 force field, the original Cornell force field parameters for χ torsions were changed without doing any fitting. The V2 term of χ torsion parameters was zeroed to improve the C2′-endo sugar puckering phase angle for DNA residues.(69) This effect, however, changes the whole predicted potential energy surface of the nucleosides, which, therefore, does not represent the QM energy surface well.
For the AMBER99χ force field, the χ torsions of C, U, A, and G were reparameterized individually. A multiconformational fitting that included the entire nucleoside with different sugar puckering was done to provide the χ torsion parameters. In the PES scan, a total of 4 × 72 = 288 data points were used in the fitting protocol for each nucleoside. The new parameter set was tested on 12 different sugar conformations (four separate conformations for each of C2′-endo, C3′-endo, and O4′-endo sugar puckering) for each nucleoside and shown to predict well the QM energy surface for these conformations (see Figures Figures22−5 and Supporting Information). The shape of the QM energy surfaces of these conformations is also predicted well by the Ode force field,(65) although not quite as well as by AMBER99χ (see Figures Figures22−5 and Supporting Information). As a result, there should not be any big difference between AMBER99χ and Ode force field(65) predictions for structural and thermodynamic properties of nucleosides.
Many reasonable combinations of parameters were tested for approximating the QM PES representing the four major conformations of each nucleoside. For instance, we tried fitting to two dihedrals with three cosine terms, four dihedrals with two cosine terms, and four dihedrals with three cosine terms. Two dihedrals with four cosine terms provided excellent fits, and more terms gave minimal improvement. As a comparison, the Ode force field(65) uses 3 dihedrals (a total of 13 Vi parameters) to represent the χ torsions, while we use 2 dihedrals (a total of 8 Vi parameters), but comparisons of the force fields to the QM potential energy surfaces shown in Figures Figures22−5 and Supporting Information reveal that AMBER99χ provides a better fit. This may be because the calculations for AMBER99χ included the entire ribose group.
It is crucial to use a force field that appropriately models the true behavior of RNA systems. Otherwise, during MD simulations, sampling space will include unphysical regions, which will cause errors in predictions. With the AMBER99χ modification, significant improvements are seen in the structural and thermodynamic predictions for cytidine and uridine in solution (Table (Table3).3). This modification should be particularly important for non-Watson−Crick regions and terminal base pairs because sampling will not include unrealistic populations of syn conformations or of C2′-endo sugar puckering. In Watson−Crick regions, the χ torsion is restricted by hydrogen bonding in base pairs, so little effect should be seen there. Thus, the AMBER99χ force field should improve structural and thermodynamic predictions for RNA.
We are grateful to the Center for Research Computing at the University of Rochester for providing the necessary computing systems (BlueHive cluster) and to the personnel to enable the research presented in this manuscript. This work was supported by NIH grant GM22939 (D.H.T.).
National Institutes of Health, United States
Frozen dihedrals during QM optimization in PES scan; restrained dihedral angles and sample restraint file used in MM minimization; chemical shift data of C, U, A and G; sugar pucker conformations used to test AMBER99, AMBER99χ, and Ode force fields; NOE data from transient NOE and SSNOE for C and U; details of population analysis results for C, U, A, and G; details of ΔG° values of C2′→C3′ and syn→anti transformations for C and U; 3J spin−spin couplings and experimentally deduced sugar puckering of C, U, A, and G; detailed analysis of Table Table3;3; analysis of individual simulations of C, U, A, and G; comparison of AMBER99, AMBER99χ, and Ode force fields to QM energy surface predictions for C, U, A, and G; intensity vs mixing time plots of C and U; rmsd vs time plots of MD simulations of C, U, A, and G; convergence analysis of the population distributions of C, U, A, and G. This material is available free of charge via the Internet at http://pubs.acs.org.