Reparameterization of RNA chi Torsion Parameters for the AMBER Force Field and Comparison to NMR Spectra for Cytidine and Uridine.

A reparameterization of the torsional parameters for the glycosidic dihedral angle, chi, for the AMBER99 force field in RNA nucleosides is used to provide a modified force field, AMBER99chi. Molecular dynamics simulations of cytidine, uridine, adenosine, and guanosine in aqueous solution using the AMBER99 and AMBER99chi 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 AMBER99chi force field, each molecular dynamics simulation time was extended to 120 ns for convergence purposes. Nuclear magnetic resonance (NMR) spectroscopy, including one-dimensional (1D) (1)H, steady-state 1D (1)H nuclear Overhauser effect (NOE), and transient 1D (1)H 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 AMBER99chi 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 degrees for the base orientation of purines. The AMBER99chi force field prefers anti conformations around 185 degrees , which is more consistent with the quantum mechanical calculations and known 3D structures of folded ribonucleic acids (RNAs). Evidently, the AMBER99chi force field predicts the structural characteristics of ribonucleosides better than the AMBER99 force field and should improve structural and thermodynamic predictions of RNA structures.

Introduction

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.

Methods

NMR

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.

Ab Initio Potential Energy Surface Scan of χ

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 1), four different sugar conformations (Table 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.

Figure 1

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. These particular structures in a−d have anti χ angles and C2′-endo sugar conformations.

Table 1

Dihedral Angles Used to Create the Four Sugar Conformations (sc) for Each Nucleoside

	C2′-endo		C3′-endo
dihedral	sc 1	sc 2	sc 3	sc 4
H5T−O5′−C5′−C4′	60	60	174	174
O5′−C5′−C4′−C3′	54	54	54	54
C5′−C4′−C3′−O3′	140	140	81	81
C4′−C3′−O3′−H3T	−148	−148	−148	−148
O4′−C1′−C2′−C3′	32	32	−24	−24
C1′−C2′−O2′−HO′2	−61	21	−153	93

Force Field Fitting of χ Torsions

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: For each nucleoside, the 4 × 72 = 288 data points from eq 1 were fitted by linear least-squares to the Fourier series shown in eq 2.

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. V and V 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, V and V. The new χ torsion parameters are listed in Table 2.

Table 2

New χ Torsion Parameters for Adenosine, Guanosine, Cytidine, and Uridine

nucleoside	torsion	n	Vn	nucleoside	torsion	n	Vn
adenosine	O4′−C1′−N9−C8	1	1.355570	cytidine	O4′−C1′−N1−C6	1	0.331762
		2	0.504875			2	0.592225
		3	−1.699430			3	−3.108180
		4	0.152425			4	−0.116806
	C2′−C1′−N9−C8	1	1.603540		C2′−C1′−N1−C6	1	1.724800
		2	−0.278197			2	−0.62684
		3	1.267980			3	2.287890
		4	0.228818			4	0.0664267
guanosine	O4′−C1′−N9−C8	1	0.835436	uridine	O4′−C1′−N1−C6	1	0.0409516
		2	0.789849			2	0.604617
		3	0.351892			3	−2.686990
		4	0.183535			4	−0.0104774
	C2′−C1′−N9−C8	1	1.047920		C2′−C1′−N1−C6	1	1.235900
		2	−0.0516452			2	−0.683638
		3	−0.905523			3	2.277010
		4	0.131907			4	0.147500

MD Simulations of Cytidine, Uridine, Adenosine, and Guanosine

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.

Results

NMR Results for Cytidine and Uridine

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/(r)6, where r 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, NOE 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 3 (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:

Table 3

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

	base orientation, % anti			sugar pucker, % C3′-endo
	(ΔG°syn→anti, kcal/mol)			(ΔG°C2′→C3′, kcal/mol)
	AMBER99	AMBER99χ	NMRb	AMBER99	AMBER99χ	NMRc
C	30	66	87	27	54	60
	(0.49)	(−0.45)	(−1.07)	(0.58)	(−0.11)	(−0.24)
U	28	83	93	35	55	56
	(0.55)	(−0.95)	(−1.45)	(0.36)	(−0.13)	(−0.15)
A	15e	13f	−	24	32	37d
G	11e	24f	−	35	54	41d

For a transformation of A→B, ΔG°A→B = −RTln(K), where R = 1.987 cal K−1 mol−1, T is the temperature in kelvins, and K is the ratio of the concentrations of each species, [B]/[A] (see Supporting Information).

Measurements of the syn/anti proportions of pyrimidines were extracted from transient NOE experiments at 10 °C, while the simulations were done at 300 K (27 °C). NMR spectra for C at 2 and 10 °C indicate essentially no temperature dependence for the syn→anti equilibrium (see Supporting Information), so all ΔG°’s were calculated at 300 K.

These values are for 30 °C (see Supporting Information).

These values are for 0.2 mM samples of A and G at 30 °C where there may be some association (see Supporting Information).

These values represent populations of high anti conformations with χ ≈ 310° (see Supporting Information and Figure 7).

These values represent populations of anti conformations with χ ≈ 185° (see Supporting Information and Figure 7).

Figure 7

Population distribution of adenosine and guanosine using AMBER99 (a and b, respectively), and AMBER99χ (c and d, respectively) force fields. PSE (y-axis) and CHI (x-axis) stand for the pseudorotation and χ dihedral angles. Table 4 shows the predicted populations of (i), (ii), (iii), (iv), (v), and (vi). PSE angles of 18° and 162° correspond to C3′-endo and C2′-endo sugar pucker, respectively. χ angles of 200°, 300°, and 60° correspond to anti, high-anti, and syn conformations, respectively.(70)

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 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 3.

Comparison of Force Field to QM Energies

Figures 2−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).

Figure 2

Total energy (in kcal/mol) vs O4′−C1′−N1−C6 of cytidine with AMBER99 (black), AMBER99χ (red), QM (green), and Ode force field (blue) for: (a) sc 1, (b) sc 2, (c) sc 3, and (d) sc 4 (see Table 1). For visualization purposes, minimum energies of each curve are set to zero. Anti, high anti, and syn base orientations correspond to x-axis ranges of 0−70°, 100−180°, and 200−300°, respectively, because the x-axis is χ + 180° to be consistent with the AMBER94 force field.(2)

Figure 5

Total energy (in kcal/mol) vs O4′−C1′−N9−C8 of guanosine with AMBER99 (black), AMBER99χ (red), QM (green), and Ode force field (blue) for: (a) sc 1, (b) sc 2, (c) sc 3, and (d) sc 4 (see Table 1). For visualization purposes, minimum energies of each curve are set to zero. Anti, high anti, and syn base orientations correspond to x-axis ranges of 0−70°, 100−180°, and 200−300°, respectively, because the x-axis is χ + 180° to be consistent with the AMBER94 force field.(2)

Total energy (in kcal/mol) vs O4′−C1′−N1−C6 of uridine with AMBER99 (black), AMBER99χ (red), QM (green), and Ode force field (blue) for: (a) sc 1, (b) sc 2, (c) sc 3, and (d) sc 4 (see Table 1). For visualization purposes, minimum energies of each curve are set to zero. Anti, high anti, and syn base orientations correspond to x-axis ranges of 0−70°, 100−180°, and 200−300°, respectively, because the x-axis is χ + 180° to be consistent with the AMBER94 force field.(2)

MD Simulations of Cytidine, Uridine, Adenosine, and Guanosine with AMBER99 and AMBER99χ

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 6 and 7. Table 3 shows the force field predictions of base orientation and sugar pucker for each nucleoside (also see Table 4 and Supporting Information). Analyses of the individual MD simulations show at least seven syn↔anti transformations for each (see Supporting Information).

Figure 6

Population distribution of cytidine and uridine using AMBER99 (a and b, respectively) and AMBER99χ (c and d, respectively) force fields. PSE (y-axis) and CHI (x-axis) stand for the pseudorotation and χ dihedral angles. Table 4 shows the predicted populations of (i−iv). PSE angles of 18° and 162° correspond to C3′-endo and C2′-endo sugar pucker, respectively. χ angles of 200°, 300°, and 60° correspond to anti, high-anti, and syn conformations, respectively.(70)

Table 4

Population Analysis Results for C, U, A and G of the AMBER99 and AMBER99χ Force Fieldsa

	i (%)	ii (%)	iii (%)	iv (%)	v (%)	vi (%)
AMBER99
cytidine	52	16	19	11	−	−
uridine	47	24	17	11	−	−
adenosine	57	21	−	−	12	3
guanosine	54	31	−	−	7	4
AMBER99χ
cytidine	20	11	23	43	−	−
uridine	9	8	36	47	−	−
adenosine	59	24	5	8	−	−
guanosine	33	39	9	15	2	−

Regions of (i) syn/C2′-endo, (ii) syn/C3′-endo, (iii) anti/C2′-endo, (iv) anti/C3′-endo, (v) high anti/C2′-endo, and (vi) high-anti/C3′-endo (Figures 6 and 7).

Discussion

Table 3 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 6).

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 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 6).

The AMBER99χ force field predicts A and G to have 13% and 24% anti conformation (Table 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 3), respectively, but with a χ dihedral angle around 310° (Figure 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 4−5 and Supporting Information, where the x-axis, however, is χ + 180°).

Figure 4

Total energy (in kcal/mol) vs O4′−C1′−N9−C8 of adenosine with AMBER99 (black), AMBER99χ (red), QM (green), and Ode force field (blue) for: (a) sc 1, (b) sc 2, (c) sc 3, and (d) sc 4 (see Table 1). For visualization purposes, minimum energies of each curve are set to zero. Anti, high anti, and syn base orientations correspond to x-axis ranges of 0−70°, 100−180°, and 200−300°, respectively, because the x-axis is χ + 180° to be consistent with the AMBER94 force field.(2)

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 3 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 2−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 2−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 2−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 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.