Stacking in RNA: NMR of Four Tetramers Benchmark Molecular Dynamics.

Molecular dynamics (MD) simulations for RNA tetramers r(AAAA), r(CAAU), r(GACC), and r(UUUU) are benchmarked against 1H-1H NOESY distances and 3J scalar couplings to test effects of RNA torsion parametrizations. Four different starting structures were used for r(AAAA), r(CAAU), and r(GACC), while five starting structures were used for r(UUUU). On the basis of X-ray structures, criteria are reported for quantifying stacking. The force fields, AMBER ff99, parmbsc0, parm99χ_Yil, ff10, and parmTor, all predict experimentally unobserved stacks and intercalations, e.g., base 1 stacked between bases 3 and 4, and incorrect χ, ϵ, and sugar pucker populations. The intercalated structures are particularly stable, often lasting several microseconds. Parmbsc0, parm99χ_Yil, and ff10 give similar agreement with NMR, but the best agreement is only 46%. Experimentally unobserved intercalations typically are associated with reduced solvent accessible surface area along with amino and hydroxyl hydrogen bonds to phosphate nonbridging oxygens. Results from an extensive set of MD simulations suggest that recent force field parametrizations improve predictions, but further improvements are necessary to provide reasonable agreement with NMR. In particular, intramolecular stacking and hydrogen bonding interactions may not be well balanced with the TIP3P water model. NMR data and the scoring method presented here provide rigorous benchmarks for future changes in force fields and MD methods.

Introduction

RNA is being increasingly recognized as important for many functions in the cell[1−6] and as a potential target for therapeutics.[7−11] Molecular dynamics (MD) simulations can be used to predict structure,[12,13] mechanism,[14] and dynamics[15−18] which can then provide insight into function. RNA MD has been done both atomistically[19−21] and with coarse-graining[22] and has a wide spectrum of applications.[23−34]

MD uses computationally inexpensive classical mechanics to approximate molecular movements and equilibria.[35] In MD simulations, potential energy is expressed in terms of bond lengths, angles between bonds, torsion angles, Lennard-Jones,[36] and Coulomb pairwise interactions. RNA parametrizations have been developed for CHARMM[37] and AMBER[35,38−40] MD packages. The lack of consistent agreement between MD simulations of RNA and experimental results, however, has revealed weaknesses in force field parameters. For example, reparametrizations of torsion terms have improved agreement with experiment.[41−44] Additionally, the Lennard-Jones term between bases[12,45,46] may not be optimal and not well balanced with solution.

Here, MD simulations with AMBER parm99,[39] parmbsc0,[41] parm99χ_Yil,[42] ff10,[43] and parmTor[44] force fields are tested against NMR spectra for RNA tetramers r(AAAA), r(CAAU), r(GACC),[47−49] and r(UUUU), similarly to previous work done on proteins.[50,51] Tetramers were chosen to minimize NMR overlap, ensuring that NOE cross-peaks between sugars and bases could be identified relatively easily and a maximum number of 3J coupling constants could be measured. Solvated tetramers are relatively small systems compared to biologically relevant RNAs, enabling more rapid computations than possible with larger systems. Sequences were chosen to preclude duplex formation in NMR experiments because single-stranded systems allow conformational freedom and thus rigorous testing of particular aspects of the force field. The temperature dependence of 1D NMR spectra showed that r(CAAU)[52] and r(GACC)[47] do not exhibit cooperative transitions and thus are single-stranded in solution. While intramolecular base–base hydrogen bonding is theoretically possible in a tetramer,[49] there is no NMR evidence for it in the tetramers studied here. Thus, the comparisons focus on other aspects of the force field.

Poly(A) stacks strongly,[53] while poly(U) is unstacked in aqueous solution as shown by light scattering, viscometry, sedimentation velocity,[54,55] and SAXS and smFRET,[56] so each force field’s predictions of stacking were tested as well on r(AAAA) and r(UUUU). Three-dimensional RNA structures tend to be A-form in Watson–Crick[57] base-paired regions, but single stranded regions can be variable.[58−62] Thus, accurate modeling of single-stranded regions requires accurate force field parametrization. Comparisons between MD simulations and experiments for a few systems indicate this has yet to be achieved.[47,49,63] The results presented here expand the sequence dependence of comparisons and provide working hypotheses for potential improvements in force fields.

RNA Force Field Parametrizations

The force fields use different nucleic acid backbone torsion[64] energy functions, as briefly described below and in Figures 1 and 2 and Table 1.

Figure 1

RNA atom names and torsions[64,65] are defined in panel A. Note that δ = C5′–C4′–C3′–O3′ while ν3 = C2′–C3′–C4′–O4′. Panel B illustrates C3′–endo (N) and C2′–endo (S) conformations. Hydrogens have been omitted for clarity.[66]

Figure 2

IUPAC torsion angle regions.[93−95]

Table 1

Summary of Force Field Torsion Reparametrizations Relative to AMBER ff99a

parametrization	altered torsions	reference
parmbsc0	α, γ	(41)
parm99χ_Yil	χYil	(42)
ff10	α, γ, χ ff10	(43)
parmTor	α, β, γ, ϵ, ζ, χ Yil	(44)

All α/γ changes are from parmbsc0 (cf. Figure 1 for torsion definitions). Parm99χ_Yil and ff10 use different quantum mechanical methods to derive χ. All force fields are based on the ff94 partial charge[67] and van der Waals model.

ff99.[39] In vacuo MP2/6-31G(d)[68,69] quantum mechanics (QM) was used to estimate the χ torsion energy profiles for a small model system designed to mimic DNA rather than RNA and did not include the entire sugar. The atoms were arranged to model a C2′–endo sugar pucker, rather than RNA’s preferred C3′–endo sugar pucker.[65] AMBER ff99 has been used extensively in diverse applications.[70−75]

parmbsc0.[41] Using a larger model system than the original parm99, a similar MP2/6-31G(d) potential energy surface scan was done to model DNA α/γ backbone torsions, resulting in more accurate and longer time scale simulations for many DNA and RNA sequences.[76−79]

parm99χ_Yil.[42,47] Complete RNA nucleoside model systems were used to gain more accurate χ energy profiles with MP2/6-31G(d) QM. The new parameters result in increased barrier heights around χ and alter the equilibria to favor T, or anti χ, orientations. This greatly improved agreement with NMR spectra for cytosine, uridine,[42] and r(GACC).[47]

ff10,[43] also known as “parmOL”. This force field used MP2 theory[68] with the complete basis set extrapolation[80,81] and a COSMO solvation model[82] to reparametrize χ at a higher QM level than parm99χ_Yil. FF10 includes parmbsc0’s revised α/γ. This improved simulation accuracy for several tetraloop hairpins, especially by preventing ladder-like structures.[83]

parmTor.[44] New β, ϵ, and ζ torsion profiles were generated with methods similar to parm99χ_Yil. ParmTor also incorporates parmbsc0’s α and γ parameters and was able to improve prediction of the differences in thermodynamic stability of tetramer duplexes formed by G–C or isoG–isoC base pairs, where isoG and isoC have amino and carbonyl groups transposed relative to G and C.

Methods

NMR

Tetramers, r(AAAA) and r(UUUU), were synthesized in Poznan and 1.3 μmol of r(CAAU) was purchased from Dharmacon. Samples were dissolved in 300 μL of 30 mM phosphate buffer, 150 mM NaCl, and 0.5 mM Na2EDTA, at pH 7.4 in water. The buffer was vacuum centrifuged twice with 600 μL of 99% D2O (Cambridge Isotope Laboratories) and finally dissolved in 99.99% D2O (Cambridge Isotope Laboratories). r(GACC) 1H chemical shift assignments and NOESY volumes[84] are from Yildirim et al.[47] A r(GACC) 1H–13C HMQC[85,86] spectrum was collected on a Varian Inova 500 MHz spectrometer, as 13C chemical shifts have been shown to differentiate χG+ and χT.[87]All tetramers’ spectra were collected at 2 °C to maximize NOE volumes. 1H–1H NOE, 1H–31P HETCOR, and 1H–13C HMQC data for r(CAAU) were collected on a Varian Inova 600 MHz spectrometer and processed with NMRPipe.[88] r(AAAA), r(CAAU), and r(UUUU) NOESY spectra were collected at 200, 400, 600, and 800 ms with 31P decoupling to minimize overlap and maximize signal-to-noise of H3′, H5′, and H5″. NOESY spectra at 2 °C and an 800 ms mixing time were obtained with and without 31P decoupling on a Varian Inova 600 MHz spectrometer to determine 3JH5′–P5′, 3JH5″–P5′, and 3JH3′–P3′. NOESY volumes, V, were integrated with the box method in Sparky 3.113.[89]

NMR distances between any 1H atoms i and j were calculated from eq 1,Errors were quantified with eqs 1a and 1b.[90,91]

The scaling constant c and standard deviation cSD were determined from the calculated scaling constants of known 1H–1H distances (Supporting Information Table 3). The x term was set to a minimum of 1.0 and increased until all of the known 1H–1H distances are within the error bounds (Supporting Information Table 4). Verr is the standard deviation of 20 randomly selected areas of the spectrum, each area approximately the same as the NOE cross-peak. To minimize the influence of spin-diffusion, 2D NOESYs with a mixing time of 200 ms were used for NOE volume measurements. Volume errors due to partial saturation effects are estimated to be within the error bounds.

The proportion of C2′/C3′-endo conformation of a sugar was estimated with eq 2.[92] The derivation is shown in Supporting Information Section 2.1.

Molecular Dynamics

Multiple starting structures can test a force field’s ability to model an RNA accurately. For example, if a starting structure is in a conformation not present in a solution ensemble, then the force field should be able to correct that given enough time. In addition, for a converged simulation, the ensemble should not depend on starting structure. Therefore, multiple starting structures not representative of the experimental ensemble were used in assessing the force fields.

For all four tetramers, four starting structures were used, A-form (N-anti), N-syn, S-anti, and S-syn. Here “N” refers to a ribose C3′-endo pseudorotation phase angle,[65] and “anti” or “syn” refers to the χ conformation (Figure 1). A-form starting structures were generated with nucgen[97] and the others were generated (Figure 3) with a simulated annealing procedure.[90,91] An additional starting structure called “>5 Å” was generated with >5 Å restraints on nH1′–(n + 1)H6 to provide a nonstacked structure expected for r(UUUU). All 17 structures were visually inspected to ensure no chiral inversions had occurred during simulated annealing. FF99 was run only with r(CAAU) and r(GACC), where its predictions were worse than parm99χ_Yil, as also seen with previous tests.[42,44,47,90] Starting structures were solvated in an 8.65 Å TIP3P[98] truncated octahedral box and run at 275 K to match NMR spectra. Each system was neutralized with three Na+ ions in AMBER’s LEAP program. Running times for all the 76 simulations range from 7.7 to 11.9 μs and are listed in Supporting Information Table 34. Details on NMR agreement of starting structures are presented in Supporting Information Section 7. This work tests four force fields (cf. Table 1) on four tetramers with each in four starting structures and a fifth for r(UUUU). Minimization, equilibration, and production runs of all tetramers followed a published protocol.[91] AMBER input files are included in Supporting Information Section 6.

Figure 3

MD starting structures for r(CAAU) and r(GACC). Image was generated with PyMOL 1.5.0.1.[96]

Comparison between MD and NMR

For an ensemble of conformations in rapid exchange, NOESY[99] spectra report average properties. MD predictions are benchmarked against 1H–1H NOE volumes and 3J scalar couplings. MD predictions of 3J scalar couplings were given an error of ±1.5 Hz to account for measurement and Karplus function errors. Each structure starts at varying accuracy with regard to NMR spectra, and NMR agreement was scored after the first 500 ns of a simulation.

Each of the NMR observables was averaged over an MD simulation and the score graded as a percentage. This was done similarly to the Locked Nucleic Acid (LNA),[100] L(CAAU),[91] where 1H–1H MD distances were computed via eq 3:where N is the number of MD trajectory points and r is the distance at time i. The negative sixth power is necessary to weight the ensemble, as an arithmetic mean would would not model this effect. L(CAAU) was evaluated as the mean of distances calculated for each 20 ns interval, while here N represents the entire simulation. All reported NMR scores are for the entire simulation after 500 ns. To provide insight into the degree of predicted variation around an “average” structure as was done for L(CAAU), Supporting Information Section 8 contains plots of percentage NMR agreement for each 20 ns window vs simulation time. Here, eq 3 is used to evaluate all NOEs: measured, overlapped, and predicted but not observed.

Measured NOEs have measured upper and lower distances, cf. Supporting Information Tables 12 and 13, which are calculated via eqs 1a and 1b.

Overlapped NOEs are visible in NOESY spectra, but do not have reliable volume information. Overlapped NOEs are counted as accurate if the MD distance is <7.0 Å.

Predicted but not observed NOEs, i.e., absent NOEs, are 1H–1H pairs predicted to have an NOE distance ≤5.0 Å but are not visible in NOESY spectra. All potential nonsequential internucleotide NOEs, between nucleotides 1–3, 1–4, and 2–4, are listed as possible NOEs. If an MD simulation predicts an NOE with average distance < 5.0 Å calculated via eq 3, but the predicted NOE is not observed in NOESY spectra, then the predicted NOE is scored as incorrect. This term is important because a large number of intercalated structures are predicted by MD but not observed by NMR. For example, many predicted structures have base 1 stacked between bases 3 and 4, forming a 3–1–4 intercalation, but no NOEs are observed between nucleotide 1 and nucleotides 3 and 4.

Sugar pucker can be quantified with three 3J scalar couplings following Marino’s nomenclature:[101] ν1 = O1′–C1′–C2′–C3′, ν2 = C1′–C2′–C3′–C4′, and ν3 = C2′–C3′–C4′–O1′. Haasnoot parameters[102,103] were replaced with parameters fit to empirical 1H–1H torsion angles from the UUCG tetraloop[104] determined via dipolar cross-coupled relaxation rates[105] (Supporting Information Section 4). Torsion angles and 3J were fit to a Karplus equation to yield eq 4. These functions are compared to Davies[106] and Haasnoot[102] predictions in Supporting Information Section 4. Here, all MD predictions for all ribose torsion angles are expressed as scalar couplings via eq 4:where n is 1, 2, or 3, i is the MD trajectory point, and N is the total number of simulation data points.

Similarly, both 3J values associated with γ can be predicted:[107]Both values must be within error bounds of observed 3J to be counted correct. The reported error is the least of the two possible errors in Hz. Similarly for β,[108,109]and for ϵ[108,109]Note that all scalar coupling equations have numerous possible degeneracies.

Agreement between MD simulations and NMR spectra was quantified via eq 8, which is similar to positive predictive value.[110]Without considering predicted but not observed NOEs, intercalated structures show unreasonably high NMR agreement.

Criteria for Stacking

A qualitative way to rapidly consider the structures generated during simulations is to score base stacking. While the word “stacking” is frequently used in RNA literature, there is no single definition. One definition states a “stacked” nucleotide has ≤4.0 Å between non-hydrogen atoms, dihedral angle between the planes of ≤30°, and overlap between the bases.[111] Another stacking definition uses Cartesian coordinates[63] but lacks a way to quantify the amount of stacking.

Here, we define stacking in terms of three criteria based on A-form X-ray structures (Figure 4 and Supporting Information Table 2). Each base plane is defined by vectors, a⃗ and b⃗, whose cross products a⃗ × b⃗ and b⃗ × a⃗ define each base’s normal vectors (Figure 4). Vectors a⃗ and b⃗ are far apart from one another to minimize out-of-plane distortions. For adenine, a⃗ is defined from the Center of Mass (CoM) to C8 and b⃗ is defined from the CoM to N6. Similarly, for guanosine a⃗ = (CoM → C8) and b⃗ = (CoM → O6). For cytosine, a⃗ = (CoM → O2) and b⃗ = (CoM → N6), and for uracil a⃗ = (CoM → O2) and b⃗ = (CoM → O6). As described below, the distance between CoMs and the angles, ω and χ, are used to provide a rough measure of stacking. Because Ξ can be positive or negative, depending on whether the base–base alignment is either relatively parallel or perpendicular (Figure 5), the stacking score ranges from 100% to −100%, respectively. Only positive percentages score as stacking.

Figure 4

PDB structure 157D(112) residues 13C and 14G shown to illustrate the distance, d0, as well as ω and Ξ angles used to quantify stacking. The 5′ and 3′ subscripts specify the 5′ and 3′ bases, respectively. For this case, d0 = 4.5 Å, ω = 40.7°, and Ξ = 17.3°.

Figure 5

Ξ and ω compared.

The first variable to consider whether two bases are “stacked” is the distance, d0, between their centers of mass. Based on CCSD(T) studies of uracil[63] and adenine[45] dimers, we consider bases to be unstacked if d0 > 5.0 Å. If d0 ≤ 3.5 Å, the stacking score is incremented +1 for d0. If 3.5 Å < d0 < 5.0 Å, then the score is decreased as r–3 from 1 to 0, where r is the distance between the centers of mass. If d0 > 5.0 Å, ω and Ξ are not computed. These distances were chosen based on X-ray statistics (Supporting Information Table 2).

The angle ω (“oh-mega” for “overlap”), is a measure of overlap between bases, analogous to the angle between steps on a staircase (Figures 4 and 5). The ω angle is defined in terms ofThe ω angle is then computed by the law of cosines:If ω ≤ 25°, the stack score is incremented +1. If 25° < ω ≤ 50°, the score is linearly decreased from 1 to 0. These angles were chosen based on X-ray statistics (Supporting Information Table 2). If ω > 50°, the base is not considered “stacked” and Ξ is not computed.

the distance between the centers of mass (d0),

the length of the 5′ base’s normal vector (d1), and

the distance between the closest normal vector of the 5′ base to the 3′ base’s center of mass (d2).

Ξ angles are the minimum angle between each base’s normal vectors (eq 10 and Figure 5). Ξ = 0° indicates vectors normal to each base plane are parallel, and the bases stack similarly to the shape of the Greek letter Ξ. The cross-product identity |a⃗ × b⃗| = |a⃗||b⃗| sin θ, where θ is the angle between vectors a⃗ and b⃗, implies Ξ can be computed via eq 10:Ξ quantifies whether a given base pair is a parallel stack or T-shape. If 45.0° < χ < 135.0°, then the stack score is multiplied by −1, which indicates a T-shape.

The stack score by this method ranges from −2 to 2. It is reported as a percentage from −100% to 100%. Here, two bases are considered stacked if the stacking percentage is (>50%).

Results

Comparison of MD simulations with NMR spectra requires assignment of NMR resonances to particular atoms. Typical NOESY walk spectra for this purpose are shown in Figure 6 for r(AAAA), r(CAAU), and r(UUUU). Additional spectra are shown in Supporting Information Section 5. All 1H, 13C, and 31P measured chemical shifts are reported in Supporting Information Table 6. C atoms not bonded to H are not detected by HMQC and were not assigned. There was insufficient concentration to measure natural abundance 15N chemical shifts. Distances with upper and lower bounds from measured NOEs are shown in Supporting Information Tables 12 and 13.

Figure 6

NOESY spectra for tetramers r(AAAA), r(CAAU), and r(UUUU) at 800 ms mixing time with 31P decoupling, drawn with Sparky 3.113.[89] NOESY walks are outlined in red. r(GACC) spectra are from Yildirim et al.[47] and are not repeated here. More images of NMR spectra are available in Supporting Information Section 5. Note that the dispersion of chemical shifts for r(UUUU) is much less than for r(AAAA) and r(CAAU), consistent with r(UUUU) residues having much more similar chemical environments.

NMR scalar couplings provide information on torsion angles, which can be compared to MD simulations. All 3J scalar couplings are given in Table 2. All nonoverlapped 3JH5′/5″–P5′ couplings showed 90–100% fractions in the βT conformer,[103] together with all three ϵ torsions being at approximately 210° (Table 2). All γ are majority γG+, as implied from Wijmenga’s Table 4.[103] It is likely that 2H exchange with solvent and the lack of a phosphate group distort γ1 to a more highly populated γT state. Most 3JH1′–H2′ are ≤2 Hz, implying that sugar puckers are mostly C3′–endo in solution[92] (Table 2). The exceptions are r(UUUU) and each 3′ terminal ribose, which show higher populations of C2′-endo. This is likely due to solvent exposure. For each tetramer, H1′, H2′, C1′, A and G H8 and C8, and C and U H6 chemical shifts (Supporting Information Table 6) imply that all four nucleotides are χT.[87]

Table 2

Torsional Parameters and/or 3J for A-Form RNA and Measured for Tetramersa

torsion	A-form[101,104,113]	r(AAAA)	r(CAAU)	r(GACC)	r(UUUU)
β2	T	3.8, 1.0; T	3.7, 2.2, ∼96% T	3.7, 0.9; ∼100% T	2.0, 2.0; ∼100% T
β3	T	3.0, 1.0; T	3.5, <1, ∼100% T	4.0, 2.0; ∼100% T	2.0, 2.0; ∼100% T
β4	T	3.2, 1.0; T	3.8–4.3, 3.3, ∼92% T	4.4, 2.0; ∼95% T	2.0, 2.0; ∼100% T
γ1	G+	3.8, 2.0; ∼77% G+	1.6, ∼3.75, ∼82% G+	4.6, 2.3; ∼66% G+	3.7, 2.8; ∼70% G+
γ2	G+	∼2, ∼1; ∼100% G+	∼2.8,∼2.1; ∼87% G+	∼2, ∼2; ∼96% G+	≈2.5, ≈2.5; ∼86% G+
γ3	G+	2.0, 2.0; ∼96% G+	2.0, 1.3; ∼93% G+	1.5, <1; ∼100% G+	<≈2.5, ≈2.5; ∼86% G+
γ4	G+	2.0, 2.0; ∼96% G+	1.0, <1; ∼70% G+	1.8, 1.3; ∼100% G+	overlap
ϵ1	185–280°	8.0–8.9; 207–233°	9.4–9.8; 218–230°	9.3; 213–230°	8.2; 210–218°
ϵ2	185–280°	8.6–8.8; 210–225°	8.4; 211–219°	9.1; 213–229°	7.8; 207–215°
ϵ3	185–280°	8.3–8.4; 209–221°	8.3; 210–219°	9.0; 212–227°	7.7; 207–215°
ν11	0.2–1.3	2.0	≈1	1.8	4.5
ν21	3.3–6.0	3.4	4.5	4.4	4.7
ν31	8.8–11.5	5.4–6.0	8.4	7.9	5.3
ν12	0.2–1.3	1.0	≈1	<1	5.2
ν22	3.3–6.0	3.6	4.6	4.4	4.6
ν32	8.8–11.5	7.3	8.3	8.9	5.5
ν13	0.2–1.3	1.0	1.4	<1	5.2
ν23	3.3–6.0	4.0	4.3	4.0	4.8
ν33	8.8–11.5	8.0	8.2–8.4	8.8	5.4
ν14	0.2–1.3	2.75	3.7	2.5	3.7
ν24	3.3–6.0	overlap	3.7	4.6	overlap
ν34	8.8–11.5	6.0	overlap	7.2	overlap

NMR 3J scalar coupling values imply A-form structure for r(AAAA), r(CAAU), and r(GACC), where nucleotides are numbered 1 to 4 starting at the 5′ nucleotide (cf. Figure 3). All values are in Hz unless specified. For ribose torsions, 3JH1′–H2′ = ν1, 3JH2′–H3′ = ν2, and 3JH3′–H4′ = ν3. 3J scalar coupling values for the 3′ terminal nucleotides and for all r(UUUU) nucleotides imply they have significant populations in the C2′–endo state (eq 2). The % of β in the T conformation can be estimated by eq 11a, and the % of γ in the G+ conformation can be estimated with eq 11b.[103] There is no significance to the order of 3J with H5′/H5″ for deducing β and γ populations.

Backbone β, γ, and ϵ torsions observed via NMR fall within standard crystallographic A-form ranges (Table 2).[113]3J scalar couplings (Table 2) and NOESY distances (Supporting Information Table 13 and Supporting Information Figure 10) indicate that r(AAAA), r(GACC),[47] and r(CAAU) are majority A-form in solution (see Methods). 31P chemical shifts are consistent with A-form for r(AAAA), r(CAAU), and r(GACC). Several NOEs in r(GACC),[47] r(AAAA), and r(CAAU) suggest their 3′ nucleotide is sometimes inverted. Table 3 summarizes the number of NMR observables that can be compared to MD predictions. The most NOEs were seen with r(CAAU). Few were seen with r(UUUU), consistent with r(UUUU) likely being unstacked and/or disordered.

Table 3

Number of NMR Observablesa

tetramer	measured NOEs	overlap NOEs	β	γ	ribose	ϵ	total
r(AAAA)	23	3	3	4	11	3	47
r(CAAU)	42	6	3	4	11	3	69
r(GACC)	23	9	3	4	12	3	54
r(UUUU)	1	5	3	3	10	3	25

Number of measured NOEs depends on how many 1H–1H pairs are visible in the spectrum (cf. Supporting Information Tables 12 and 13). The NOEs (cf. Section 2.3) and 3J scalar coupling allow a direct comparison between NMR results and computational predictions. Overlap NOEs could not be definitively measured in NOESY spectra (cf. Figure 6). Only one NOE was measured with r(UUUU), which is consistent with r(UUUU) being far from an A-form structure. NOEs with H5′ or H5″ were omitted.

Molecular Dynamics

Simulations for all four tetramers with four starting structures each were completed with parmbsc0,[41] parm99χ_Yil,[42] ff10,[43] and parmTor.[44] Tetramers r(CAAU) and r(GACC) were also simulated with ff99 and a fifth starting structure was tested for r(UUUU), for a total of 76 simulations and 737 μs of simulation time (cf. Tables 4 and 5). Table 4 and Supporting Information Table 32 show for each simulation the percentage agreement with NMR, including in Table 4 the number of NOEs predicted but not observed.

Table 4

Agreement of Simulations with NMRa

NMR Mean Scores (%) and Number of NOEs Predicted but Not Observed
	force field	A-form (43)	N-syn (30)	S-anti (24)	S-syn (21)	mean
r(AAAA)	parmbsc0	40 (5)	34 (6)	37 (2)	33 (10)	36 (5.8)
	parm99χ_Yil	32 (9)	26 (31)	32 (16)	27 (15)	29 (17.8)
	ff10	37 (5)	38 (18)	37 (21)	38 (13)	38 (14.2)
	parmTor	27 (16)	14 (40)	20 (52)	21 (58)	20 (41.5)

Each starting structure’s NMR agreement is given in parentheses in the column headings. NMR agreement via eq 8 is shown in black, while the number of NOEs predicted to be below 5 Å but absent in spectra is shown in italics in parentheses. Explicit details in starting structure NMR accuracy are presented in Supporting Information section 7.

Table 5

Total Simulation Time for Each Force Field and Summary of Results from Table 4a

force field	min. score (%)	mean (%)	max. score (%)	total time (μs)
parmbsc0	18	31	43	165
parm99χ_Yil	19	31	46	166
ff10	21	33	42	162
parmTor	14	23	36	166

Mean % is an average of 16 simulations for a given force field. Simulation for r(UUUU) with “<5 Å” starting structure is also included.

The number of predicted but not observed NOEs can be a significant or even, for r(UUUU), a majority part of the denominator of eq 8. This is due to a large fraction of simulation structures that contain stacks between bases that are not nearest neighbors in the sequence (Tables 6–9). This over representation is particularly large for r(UUUU) and more generally for parmTor.

Table 6

r(AAAA) Simulation Dataa

	NMR	parmbsc0	parm99χ_Yil	ff10	parmTor
1–2 stack	present	36.1%	19.0%	24.5%	20.3%
2–3 stack	present	36.5%	43.3%	50.3%	36.9%
3–4 stack	present	30.2%	33.4%	22.6%	42.9%
1–3 stack	none	1.6%	4.7%	7.5%	24.3%
1–4 stack	none	9.0%	3.2%	4.3%	11.2%
2–4 stack	none	1.4%	7.9%	5.4%	9.5%
1–3–2 intercalation	none	0.0%	2.7%	5.6%	0.8%
1–4–2 intercalation	none	0.0%	0.6%	0.0%	1.8%
2–1–3 intercalation	none	0.0%	0.0%	0.0%	5.2%
2–4–3 intercalation	none	0.0%	2.0%	0.6%	2.3%
3–1–4 intercalation	none	0.1%	13.0%	22.2%	3.2%
3–2–4 intercalation	none	0.0%	0.7%	0.0%	2.0%
random coil	none	21.8%	15.1%	11.8%	8.2%
β RMSD (Hz)	T	0.5	0.3	0.4	6.0
ϵ RMSD (Hz)	cf. Table 2	2.5	2.7	3.0	3.1
γ RMSD (Hz)	cf. Table 2	0.8	1.4	0.7	1.2
suger pucker RMSD (Hz)	cf. Table 2	2.4	1.8	1.4	2.2

r(AAAA) simulation data shows percentage of time the simulation showed stacks, intercalations, random coil, as well as RMSD of predicted 3J couplings relative to NMR. Nucleotides are numbered 1 to 4 starting with the 5′ nucleotide. Thus, a 1–2 stack represents a structure with the bases of nucleotides 1 and 2 stacked on each other, and a 1–3–2 intercalation represents stacking of base 3 on bases 1 and 2. Intercalations with maximum <1% for all force fields are not listed. “Random coil” means conformations where all nucleotides failed the stacking definition. Note that some intercalations can occur simultaneously, viz., 1–3–2/3–2–4 and 1–3–2/3–1–4. Non-sequential stacks are only counted if there is no intercalation present. Similarly, a 2–4 stack includes 1–4–2, 2–4–3, and 3–2–4 intercalations. Sugar pucker RMSD is the sum of all errors in ν: 3JH1′–H2′ + 3JH2′–H3′, + 3JH3′–H4′.

Table 9

r(GACC)a

	NMR	parm99	parmbsc0	parm99χ_Yil	ff10	parmTor
1–2 stack	present	42.6%	55.4%	51.4%	35.4%	17.4%
2–3 stack	present	33.1%	27.4%	68.0%	68.3%	76.6%
3–4 stack	present	11.9%	11.7%	47.4%	31.1%	13.3%
1–3 stack	none	19.5%	12.1%	5.0%	12.2%	3.1%
1–4 stack	none	4.6%	4.3%	9.7%	21.9%	17.9%
2–4 stack	none	17.7%	13.2%	7.3%	8.5%	3.3%
1–3–2 intercalation	none	0.1%	0.1%	7.7%	10.7%	30.1%
3–1–4 intercalation	none	0.1%	0.0%	4.2%	4.9%	13.8%
random coil	none	14.9%	14.0%	3.3%	3.2%	3.7%
β RMSD (Hz)	T	1.0	1.2	0.4	0.6	3.3
ϵ RMSD (Hz)	cf. Table 2	3.5	3.4	2.3	2.6	5.1
γ RMSD (Hz)	cf. Table 2	1.4	0.9	1.1	0.8	1.3
suger pucker RMSD (Hz)	cf. Table 2	3.2	2.8	0.4	1.3	1.7

Cf. Table 6.

Table 5 summarizes average agreement with NMR for four force fields and lists the total simulation time for each. Supporting Information Table 32 gives RMSD values between predicted and measured distances and 3J couplings. The parmTor force field shows the least average agreement with NMR. Parmbsc0, parm99χ_Yil, and ff10 performed similarly on average when compared to NMR spectra for the four single stranded tetramers studied here (Table 5).

When parmbsc0, parm99χ_Yil, and ff10 are considered, the mean agreement with NMR results depends on sequence (cf. Table 4). The order is r(CAAU) > r(GACC) ≈ r(AAAA) > r(UUUU). None, however, averages >40% NMR agreement.

Tables 6–9 provide qualitative insight into the structures seen in the MD simulations. For r(AAAA), r(CAAU), and r(GACC), with some exceptions for parmbsc0 and parmTor, 2–3 stacks have the highest percentage of stacking between two bases. Presumably, this reflects the fact that bases 1 and 4 are more exposed to water. As expected from studies of poly U,[54,55] r(UUUU) has the least predicted stacking (as defined in Figures 4 and 5) between nearest neighbor bases and the largest predicted percentage of random coil (cf. Figure 7), where “random-coil” is defined as having no stacking above 50%. Unexpectedly, however, r(UUUU) is predicted to have a large fraction of 1–3 stacks with structures similar to that shown in Figure 7B. Additionally, r(UUUU) has a large number of predicted but not observed NOEs, which implies that predicted structures are more collapsed than expected from small-angle X-ray scattering measurements on r(U40).[56]

Figure 7

Representative structures from r(UUUU) simulation starting A-form with ff10. Structure (B) illustrates a 1–3 stack. Hydrogen bonds are shown in green.

No NMR evidence is seen for 1–4 stacks, but parm99χ_Yil, ff10, and parmTor all predict that >9% of r(GACC) structures will have 1–4 stacking. Structures of r(GACC) generated with an M-REMD approach using ff10 suggest roughly 25% will have 1–4 stacks.[49]

There is also no NMR evidence for intercalated structures, where “intercalations” are defined as a nucleotide stack N–N–N, where j < i or j > i + 1. This also contrasts with MD predictions, as seen previously for r(GACC).[47−49] Tables 6, 8, and 9 show that significant 1–3–2 and 3–1–4 intercalated structures are often predicted for r(AAAA), r(CAAU), and r(GACC). Intercalations can be very stable in simulations. For example, with ff10 starting with an A-form structure for r(CAAU), a 3–1–4 intercalation started at ≈490 ns and did not break free as of 9.2 μs (Figure 8). Similarly, parmTor’s A-form starting structure of r(CAAU) started a 3–1–4 intercalation at ≈540 ns and did not become free as of 9.6 μs.

Table 8

r(CAAU)a

	NMR	parm99	parmbsc0	parm99χ_Yil	ff10	parmTor
1–2 stack	present	29.6%	25.2%	6.9%	6.2%	8.8%
2–3 stack	present	60.1%	40.5%	75.1%	83.0%	47.4%
3–4 stack	present	32.8%	12.2%	10.6%	1.8%	19.5%
1–3 stack	none	2.7%	6.4%	2.3%	3.1%	18.3%
1–4 stack	none	2.4%	1.4%	0.5%	1.8%	1.4%
2–4 stack	none	2.8%	6.4%	2.1%	0.6%	22.7%
1–3–2 intercalation	none	0.0%	0.5%	15.6%	35.0%	17.4%
2–1–3 intercalation	none	0.0%	0.1%	0.0%	0.0%	1.7%
3–1–4 intercalation	none	0.1%	2.0%	44.6%	49.5%	17.1%
3–2–4 intercalation	none	0.1%	0.0%	1.7%	0.0%	0.0%
random coil	none	18.4%	29.7%	10.0%	2.3%	8.0%
β RMSD (Hz)	T	0.9	0.2	1.3	1.2	8.0
ϵ RMSD (Hz)	cf. Table 2	3.2	3.2	3.8	4.3	3.7
γ RMSD (Hz)	cf. Table 2	1.5	1.2	1.1	1.1	1.2
suger pucker RMSD (Hz)	cf. Table 2	3.3	2.7	0.8	0.5	2.4

Cf. Table 6.

Figure 8

One example of a 3–1–4 intercalation: r(CAAU) starting from an A-form structure with force field ff10 intercalates at ≈490 ns and remains intercalated. Structure shown is visualized at 5996.0 ns. Hydrogen bonds are shown in green, which lock the intercalation in place. Both ζ1 and α2 transition from G– to G+ as the intercalation occurs and Solvent-Accessible Surface Area (S.A.S.A.) decreases. Figure was drawn with Grace 5.1.23 and Avogadro.[66]

The intercalations occurred simultaneously with sharp decreases in total Solvent-Accessible Surface Area[114] (S.A.S.A.) for both ff10 (Figure 8) and parmTor. This possibly indicates that base–base interactions within the RNA are too strong[12,45] and/or that nucleotide–solvent interactions are not strong enough compared to water–water interactions.[115]

Inspection of the intercalated r(CAAU) ff10 structure at 5996.0 ns (Figure 8) shows two principal probable causes for this structure’s longevity: (i) a hydrogen bond between C1’s amino group and a U4 phosphate nonbridging oxygen and (ii) two hydrogen bonds between A3 phosphate nonbridging oxygens and C1’s HO5′ and U4’s HO3′. The defect in hydrogen bonding may be especially significant when modeling splicing reactions, which depend on the positioning of the phosphate backbone.[116−118] Such base–phosphate interactions are present in biological RNAs.[119]

The hydrogen bonds observed for r(CAAU) were also seen with r(AAAA) (cf. Figure 9). Simulations of r(CCCC) also generated amino to phosphate hydrogen bonds.[90] Unlike r(CCCC),[90] however, in r(CAAU) there were no carbonyl-amino group electrostatic attractions favoring intercalation.

Figure 9

r(AAAA) with ff10 and N-syn starting structure has a very similar intercalation as r(CAAU) in Figure 8. Unlike r(CAAU), the intercalation was present at the beginning of the production run and broke free at 7.5 μs.

Comparisons between predicted and observed sugar puckers and various torsion angles provide additional insight into strengths and weaknesses of force fields. Huang et al.[120] noted that “molecular mechanics simulations fail to adequately describe ... ribose ring puckers.” Nevertheless, parm99χ_Yil,[42] ff10,[43] and parmTor[44] all improved sugar puckers relative to parm99 and parmbsc0 simulations of r(CAAU) and r(GACC) (cf. Tables 8 and 9). Presumably, this improvement reflects the close coupling between χ and sugar pucker[65,121−124] so that quantum mechanical reparameterization of χ also lowered the RMSD between measured and predicted ribose 3J scalar coupling. Even with this improvement, however, only a plurality of C3′–endo structures are predicted instead of the majority C3′–endo seen in NMR (cf. Supporting Information Section 8).

LNAs have a CH2 bridge connecting ribose O2′ and C4′ atoms, which restricts the sugar to a C3′–endo conformation and simplifies computations. The agreement with NMR for r(CAAU) simulations can be compared to that for simulations of the LNA, L(CAAU).[91] In a 3 μs simulation starting with the A-form structure of L(CAAU), no intercalations were observed, and NMR agreement was 66% with the parm99_LNA force field. The latter had partial charges and χ reparametrized for LNA. This can be roughly compared to the 36% NMR agreement for r(CAAU) with A-form starting structure and parm99χ_Yil force field (Table 4). The comparison implies AMBER sugar parameters can be improved.[120]

Tables 6–9 compare predicted and measured 3J scalar couplings for β, ϵ, and γ. The parmbsc0,[41] parm99χ_Yil,[42] and ff10[43] force fields use the same parameters for these torsions. ParmTor[44] revised β, ϵ, and ζ via MP2/6-31G(d) quantum calculations (cf. Table 1) and the revision of β resulted in poor agreement between measured 3JH5′–P5′ and 3JH5″–P5′ couplings and those predicted with eqs 6. In general, the range in differences between measured and predicted 3J couplings is too small to choose between the parmbsc0, parm99χ_Yil, and ff10 force fields, especially because of uncertainties in the Karplus equations (eqs 4–7) as discussed in Supporting Information Section 4.5.

Discussion

RNA has many functions[6] and increasingly is being targeted with therapeutics.[125−128] Accurate determination of RNA secondary structure is relatively rapid,[129−134] but determination of 3D structure and dynamics is slow. While double helical regions can be reasonably modeled as A-form RNA, accurately modeling loop regions remains a challenge.[135] Accurate force fields for RNA would facilitate modeling loop structures, dynamics, and docking with other molecules.

Practical force fields must approximate well the many interactions driving RNA folding.[136] These include stacking, hydrogen bonding, solvation, counterion interactions, and torsion potentials. The tetramers presented here are too short to have base–base hydrogen bonds and counterion condensation.[137,138] Thus, the benchmarks presented here focus on other interactions.

Parametrization of torsions is the main difference between the force fields tested (Table 1). ParmTor had the most modifications (Table 1) but least agreement with NMR (Tables 4 and 5). The ParmTor β reparametrization gave the most dramatic change in energy vs torsion angle.[44] In particular, the β energy landscape is much flatter for parmTor than for the other force fields. While parmTor gives a worse performance here, it gave essentially the same agreement as parm99χ_Yil and better agreement than parm99 when less extensive comparisons were made between simulations and NMR spectra of r(CCCC).[90] In a comparison of force fields to predict differences in ΔG27° for forming tetramer duplexes with all GC or isoG–isoC base pairs, parmTor gave better agreement with experiment than parm99 and parm99χ_Yil.[44] Because parameters for different components of a force field are interdependent,[12] it is not clear that the new β parametrization is fundamentally worse than the parm99 one, also employed in parmbsc0, parm99χ_Yil, and ff10.

The parmbsc0, parm99χ_Yil, and ff10 simulations had similar agreement with NMR observables (Tables 4 and 5), although none had means > 40%. The structural ensembles, however, are somewhat different (Tables 6–9). For example, parmbsc0 generates the fewest intercalated structures for r(AAAA), r(CAAU), and r(GACC). Because there is no NMR evidence for intercalated structures, parmbsc0 therefore has the fewest absent NOEs lowering NMR agreement with MD via eq 8. Parm99χ_Yil and ff10 both differ from parmbsc0 by having χ reparametrized, which favors an anti conformation over a syn conformation. Predicted intercalated bases have anti conformations. In a different benchmark, the χ reparametrization in ff10 improved predictions of UUCG and GNRA hairpins relative to parmbsc0.[83] Evidently, the force fields must be tested against multiple benchmarks.

For the parm99χ_Yil and ff10 simulations of r(AAAA), r(CAAU), and r(GACC), the most common structures not in agreement with NMR (Tables 6, 8, and 9) are 1–3–2 and 3–1–4 intercalations. Both have been seen in previous simulations of r(GACC),[47−49] and 3–1–4 was seen in r(CCCC) simulations.[90] Figures 8 and 9 show examples for r(CAAU) and r(AAAA). Therefore, r(UUUU) is the only sequence thus far which is predicted not to have very favorable intercalation. It is also the only sequence lacking amino groups. Evidently, amino to phosphate interactions are one key for stabilizing predicted intercalations (Figures 8 and 9).

In addition to hydrogen bonds from the intercalated base amino group, there are hydrogen bonds from 5′-terminal HO and 3′-terminal HO groups to phosphate nonbridging oxygens that stabilize the intercalated structures generated by MD (Figures 8 and 9). Moreover, for r(AAAA), r(CAAU), and r(GACC), intercalation is also accompanied by a decrease in Solvent-Accessible Surface Area (Figure 8).

Parm99χ_Yil and ff10 predict that there is more 2–3 stacking than 1–2 or 3–4 stacking in r(AAAA), r(CAAU), and r(GACC). This is reasonable because the middle two nucleotides can be held in place by two stacking interactions while the terminal nucleotides have water on one side. r(GACC) is special in having more predicted 1–4 stacking than other sequences, which suggests extra favorable stacking between G and C. Total stacking and intercalation are least for r(UUUU) as expected from studies of poly(U).[54,55] Nevertheless, r(UUUU) has the largest number of absent NOEs, suggesting the simulations predict a relatively disordered, but collapsed, structure not consistent with NMR of r(UUUU), experiments on poly(U),[54,55] or r(U40).[56]

The parm99χ_Yil and ff10 simulations of r(AAAA), r(CAAU), and r(GACC) reveal predicted stacking and hydrogen-bonding interactions not consistent with NMR spectra. There are several possible reasons for this. The AMBER force field may have base–base van der Waals interactions that are too strong.[12] Alternatively, or in addition, nucleotide–water and water–water interactions may not be balanced well. The same reasons may explain the collapsed structures observed for r(UUUU). Approximating all the forces driving RNA conformations is difficult.[136]

RNA force fields are often used to refine 3D structures determined by NMR and crystallography or to predict 3D structure and/or dynamics when only secondary structure is known. Often, the RNA is complexed with protein. In all these cases, MD simulations are likely to perform better than reported here for unrestricted tetramers because of experimental and/or covalent restraints and additional volume exclusion. Thus, the tetramers provide a particularly rigorous benchmark for testing the approximations inherent in a classical force field.

Table 7

r(UUUU)a

	NMR	parmbsc0	parm99χ_Yil	ff10	parmTor
1–2 stack	none	5.5%	15.7%	13.8%	21.0%
2–3 stack	none	6.3%	22.1%	12.3%	20.6%
3–4 stack	none	4.5%	9.1%	8.3%	21.4%
1–3 stack	none	50.6%	36.1%	33.2%	15.4%
1–4 stack	none	1.8%	2.3%	2.9%	7.2%
2–4 stack	none	30.8%	11.5%	12.8%	7.2%
1–3–2 intercalation	none	0.0%	3.1%	0.3%	3.8%
3–1–4 intercalation	none	0.1%	4.2%	3.3%	1.5%
random coil	possible	27.6%	21.6%	30.2%	28.1%
β RMSD (Hz)	T	0.1	0.1	0.1	6.4
ϵ RMSD (Hz)	cf. Table 2	1.3	1.2	1.5	2.2
γ RMSD (Hz)	cf. Table 2	0.9	1.4	0.7	0.5
suger pucker RMSD (Hz)	cf. Table 2	2.3	0.6	1.2	0.7

Cf. Table 6.