rna is being increasingly recognized as important for many functions 
 in the cell and as a potential target for therapeutics .
molecular dynamics ( md ) simulations can be used to 
 predict structure , mechanism , and dynamics which can then provide insight into function .
rna 
 md has been done both atomistically and with coarse - graining and has a wide spectrum of applications .
md uses computationally inexpensive classical mechanics to 
 approximate 
 molecular movements and equilibria . in 
 md simulations ,
potential energy is expressed in terms of bond lengths , 
 angles between bonds , torsion angles , lennard - jones , and coulomb pairwise interactions .
the lack of consistent agreement between 
 md simulations of rna and experimental results , however , has revealed 
 weaknesses in force field parameters .
additionally , the lennard - jones term between bases may not be optimal and not well 
 balanced with solution . here
, md simulations with amber parm99 , parmbsc0 , parm99_yil , ff10 , and parmtor force fields are tested against nmr spectra 
 for rna tetramers r(aaaa ) , r(caau ) , r(gacc ) , and r(uuuu ) , 
 similarly to previous work done on proteins .
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 j 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 ) and r(gacc ) do not exhibit cooperative 
 transitions and thus are single - stranded in solution .
while intramolecular 
 base  base hydrogen bonding is theoretically possible in a tetramer , there is no nmr evidence for it in the tetramers 
 studied here .
poly(a ) stacks strongly , while poly(u ) 
 is unstacked in aqueous solution as shown by light scattering , viscometry , 
 sedimentation velocity , and saxs and smfret , so each force field s predictions of 
 stacking were tested as well on r(aaaa ) and r(uuuu ) .
crick base - paired regions , but single stranded regions 
 can be variable . 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 .
the results presented here expand 
 the sequence dependence of comparisons and provide working hypotheses 
 for potential improvements in force fields .
the force fields use different nucleic 
 acid backbone torsion energy functions , 
 as briefly described below 
 and in figures 1 and 2 and table 1.i.ff99 . in vacuo mp2/6 - 31g(d ) 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 c2endo 
 sugar pucker , rather than rna s preferred c3endo 
 sugar pucker .
amber ff99 has been used 
 extensively in diverse applications.ii.parmbsc0 . 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.iii.parm99_yil .
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 , and r(gacc).iv.ff10 , also known as  parmol  .
this force 
 field 
 used mp2 theory with the complete basis 
 set extrapolation and a cosmo solvation model 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.v.parmtor .
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 
isoc 
 base pairs , where isog and isoc have amino and carbonyl groups transposed 
 relative to g and c. ff99 . in vacuo mp2/6 - 31g(d ) 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 c2endo 
 sugar pucker , rather than rna s preferred c3endo 
 sugar pucker .
parmbsc0 . 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 .
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 , and r(gacc ) .
this force 
 field 
 used mp2 theory with the complete basis 
 set extrapolation and a cosmo solvation model 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 .
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 
isoc 
 base pairs , where isog and isoc have amino and carbonyl groups transposed 
 relative to g and c. rna atom names and torsions are defined in panel 
 a. note that  = c5c4c3o3 
 while 3 = c2c3c4o4. 
 panel b illustrates c3endo ( n ) and c2endo 
 ( s ) conformations .
parm99_yil and ff10 use different quantum mechanical methods 
 to derive . all force fields are based on the ff94 partial 
 charge and van der waals model .
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 ) .
a r(gacc ) h  c hmqc spectrum was collected on a varian inova 500 mhz spectrometer , as c chemical shifts have been shown to differentiate  and .all 
 tetramers spectra were collected at 2 c to maximize 
 noe volumes .
p hetcor , and h  c hmqc data 
 for r(caau ) were collected on a varian inova 600 mhz spectrometer 
 and processed with nmrpipe .
r(aaaa ) , 
 r(caau ) , and r(uuuu ) noesy spectra were collected at 200 , 400 , 600 , 
 and 800 ms with p 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 p decoupling on a varian inova 600 mhz 
 spectrometer to determine jh5p5 , jh5p5 , and jh3p3. noesy volumes , v , were integrated with the box 
 method in sparky 3.113 .
nmr distances 
 between any h atoms i and j were calculated from eq 1,1errors were quantified with eqs 1a and 1b.1a1b the scaling constant c and standard deviation csd were determined from the calculated scaling 
 constants of known h 
the x term was set to a minimum of 1.0 and increased until all 
 of the known h 
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
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 .
the derivation is shown in supporting 
 information section 2.1.2 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 , and  anti  or  syn  refers to the  conformation ( figure 1 ) .
a - form starting structures were generated with 
 nucgen and the others were generated 
 ( figure 3 ) with a simulated annealing procedure .
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 .
starting structures 
 were solvated in an 8.65  tip3p 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
details on nmr agreement of starting 
 structures are presented in supporting information section 7 .
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 .
for 
 an ensemble of conformations in rapid exchange , noesy spectra report average properties .
md predictions 
 of j 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 ) , l(caau ) , where h  h md 
 distances were computed via eq 3:3where n is the number of 
 md trajectory points and ri 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 .
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.i.measured noes have measured 
 upper and lower distances , cf .
supporting information tables 12 and 13 , which are calculated via eqs 1a and 1b.ii.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 .iii.predicted but not 
 observed 
 noes , i.e. , absent noes , are h 
h pairs predicted to have an noe distance 5.0  but 
 are not visible in noesy spectra .
all potential nonsequential internucleotide 
 noes , between nucleotides 13 , 14 , and 24 , 
 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 314 
 intercalation , but no noes are observed between nucleotide 1 and nucleotides 
 3 and 4 .
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
h pairs predicted to have an noe distance 5.0  but 
 are not visible in noesy spectra .
all potential nonsequential internucleotide 
 noes , between nucleotides 13 , 14 , and 24 , 
 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 314 
 intercalation , but no noes are observed between nucleotide 1 and nucleotides 
 3 and 4 .
sugar pucker 
 can be quantified with three j scalar 
 couplings following marino s nomenclature : 1 = o1c1c2c3 , 
 2 = c1c2c3c4 , 
 and 3 = c2c3c4o1. 
 haasnoot parameters were replaced with parameters 
 fit to empirical h 
h torsion angles 
 from the uucg tetraloop determined via 
 dipolar cross - coupled relaxation rates ( supporting information section 4 ) .
torsion 
 angles and j were fit to a karplus equation 
 to yield eq 4 .
these functions are compared 
 to davies and haasnoot predictions in supporting information section 4 . here ,
all md predictions for all ribose torsion angles 
 are expressed as scalar couplings via eq 4:4where n is 1 , 2 , or 3 , i is the md trajectory point , and n is 
 the total number of simulation data points . similarly , both j values associated 
 with  can be predicted:5a5bboth values 
 must be within error bounds of 
 observed j to be counted correct .
the 
 reported error is the least of the two possible errors in hz . similarly 
 for ,6a6band for 7note 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.8without considering predicted 
 but not observed 
 noes , intercalated structures show unreasonably high nmr agreement .
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 . another stacking definition uses cartesian coordinates 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 .
i.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 and adenine dimers , we consider bases to be unstacked if d0 > 5.0  . if d0
if 3.5  < d0 < 
 5.0  , then the score is decreased as r from 1 to 0 , where r is the distance 
 between the centers of mass . if d0
these distances were 
 chosen based on x - ray statistics ( supporting information table 2).ii.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 of(a)the distance between the centers of 
 mass ( d0),(b)the length of the 5 base s 
 normal vector ( d1 ) , and(c)the distance between the closest normal 
 vector of the 5 base to the 3 base s center 
 of mass ( d2).the 
  angle is then computed by the law of cosines:9if   25 , 
 the stack score 
 is incremented + 1 . if 25 < 
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.iii. 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: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 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 and adenine dimers , we consider bases to be unstacked if d0 > 5.0  . if d0
if 3.5  < d0 < 
 5.0  , then the score is decreased as r from 1 to 0 , where r is the distance 
 between the centers of mass . if d0
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 of(a)the distance between the centers of 
 mass ( d0),(b)the length of the 5 base s 
 normal vector ( d1 ) , and(c)the distance between the closest normal 
 vector of the 5 base to the 3 base s center 
 of mass ( d2).the 
  angle is then computed by the law of cosines:9if   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: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 .
pdb structure 157d(112 ) residues 
 13c and 14 g 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. the stack score by this method ranges from 2 to 2 . 

two bases 
 are considered stacked if the stacking percentage is ( > 50% ) .
  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 ) .
all h , c , and p 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 .
distances with upper and 
 lower bounds from measured noes are shown in supporting 
 information tables 12 and 13 .
noesy spectra for tetramers r(aaaa ) , r(caau ) , 
 and r(uuuu ) at 800 
 ms mixing time with p decoupling , drawn with sparky 3.113 .
r(gacc ) spectra 
 are from yildirim et al . and are not 
 repeated here .
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 
 nonoverlapped jh5/5p5 couplings showed 90100% fractions in the  conformer , together with all three 
  torsions being at approximately 210 ( table 2 ) .
it is likely that h exchange with solvent and 
 the lack of a phosphate group distort 1 to a more highly populated 
  state .
most jh1h2 are 2 hz , implying that 
 sugar puckers are mostly c3endo in solution ( 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 . nmr j 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 ) .
j scalar coupling 
 values for the 3 terminal nucleotides and for all r(uuuu ) 
 nucleotides imply they have significant populations in the c2endo 
 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 .
there is no 
 significance to the order of j with h5/h5 
 for deducing  and  populations .
backbone  ,  , and  
 torsions observed via nmr 
 fall within standard crystallographic a - form ranges ( table 2).j scalar couplings ( table 2 ) and 
 noesy distances ( supporting information table 13 and supporting information figure 
 10 ) indicate that r(aaaa ) , r(gacc ) , and 
 r(caau ) are majority a - form in solution ( see methods ) .
p chemical shifts are consistent with a - form for 
 r(aaaa ) , r(caau ) , and r(gacc ) . several noes in r(gacc ) , 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 . 

few were seen with r(uuuu ) , 
 consistent with r(uuuu ) likely being unstacked and/or disordered .
h pairs are visible in 
 the spectrum ( cf . supporting information tables 12 and 13 ) .
the noes ( cf . section 2.3 ) and j 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 .
simulations for 
 all four tetramers with four starting structures each were completed 
 with parmbsc0 , parm99_yil , ff10 , and parmtor .
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 .
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 .
mean % is an average 
 of 16 simulations 
 for a given force field .
< 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 69 ) .
this over representation 
 is particularly large for r(uuuu ) and more generally for parmtor .
r(aaaa ) 
 simulation data shows percentage 
 of time the simulation showed stacks , intercalations , random coil , 
 as well as rmsd of predicted j couplings 
 relative to nmr .
thus , a 12 stack represents a structure 
 with the bases of nucleotides 1 and 2 stacked on each other , and a 
 132 intercalation represents stacking of base 3 on 
 bases 1 and 2 .
similarly , 
 a 24 stack includes 142 , 243 , 
 and 324 intercalations .
sugar pucker rmsd is the sum 
 of all errors in  : jh1h2 + jh2h3 , + jh3h4. 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 j couplings .
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 .
the order is r(caau ) > r(gacc )  r(aaaa ) 
 > r(uuuu ) .
tables 69 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 , 23 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 , 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 13 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 ) .
representative structures 
 from r(uuuu ) simulation starting a - form 
 with ff10 . structure ( b ) illustrates a 13 stack .
no nmr evidence is seen for 14 stacks , but parm99_yil , 
 ff10 , and parmtor all predict that > 9% of r(gacc ) structures will 
 have 14 stacking .
structures of r(gacc ) generated with an 
 m - remd approach using ff10 suggest roughly 25% will have 14 
 stacks .
there is also no nmr evidence 
 for intercalated structures , where 
  intercalations  are defined as a nucleotide stack ni 
this 
 also 
 contrasts with md predictions , as seen previously for r(gacc ) .
tables 6 , 8 , and 9 show that significant 132 and 314 
 intercalated structures are often predicted for r(aaaa ) , r(caau ) , 
 and r(gacc ) .
, with ff10 starting with an a - form structure for r(caau ) , 
 a 314 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 314 intercalation at 540 ns and 
 did not become free as of 9.6 s .
one example of a 314 
 intercalation : r(caau ) starting 
 from an a - form structure with force field ff10 intercalates at 490 
 ns and remains intercalated .
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 . )
the intercalations occurred simultaneously with sharp decreases 
 in total solvent - accessible surface area ( s.a.s.a . ) for both ff10 ( figure 8) and parmtor . 

this possibly indicates that base  base interactions within 
 the rna are too strong and/or that nucleotide 
. 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 .
the hydrogen bonds observed 
 for r(caau ) were also seen with r(aaaa ) 
 ( cf .
unlike r(cccc ) , however , in r(caau ) 
 there were no carbonyl - amino group electrostatic attractions favoring 
 intercalation .
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 . noted 
 that  molecular mechanics simulations fail to adequately describe 
 ... ribose ring puckers .
 nevertheless , parm99_yil , ff10 , and parmtor 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 so that quantum mechanical reparameterization 
 of  also lowered the rmsd between measured and predicted ribose j scalar coupling .
even with this improvement , 
 however , only a plurality of c3endo structures are 
 predicted instead of the majority c3endo 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 c3endo 
 conformation and simplifies computations .
the agreement with nmr for 
 r(caau ) simulations can be compared to that for simulations of the 
 lna , l(caau ) . 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 . 

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 ) .
tables 69 compare predicted 
 and measured j scalar 
 couplings for  ,  , and . the parmbsc0 , parm99_yil , and ff10 force fields use the same 
 parameters for these torsions .
parmtor revised  ,  , and  via mp2/6 - 31g(d ) quantum calculations 
 ( cf .
table 1 ) and the revision of  resulted 
 in poor agreement between measured jh5p5 and jh5p5 couplings and those predicted 
 with eqs 6 .
in general , the range in differences 
 between measured and predicted j couplings 
 is too small to choose between the parmbsc0 , parm99_yil , and 
 ff10 force fields , especially because of uncertainties in the karplus 
 equations ( eqs 47 ) as discussed in supporting information section 4.5 .
accurate determination of rna 
 secondary structure is relatively rapid , 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 .
accurate force fields for rna would facilitate modeling loop structures , 
 dynamics , and docking with other molecules .
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 .
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 .
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 ) . 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 . because parameters for different components 
 of a force field are interdependent , 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 69 ) . 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
the  reparametrization in ff10 improved 
 predictions of uucg and gnra hairpins relative to parmbsc0 .
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 132 and 314 
 intercalations .
both have been seen in previous simulations of r(gacc ) , and 314 was seen in r(cccc ) simulations .
therefore , 
 r(uuuu ) is the only sequence thus far which is predicted not to have 
 very favorable intercalation .
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 23 
 stacking than 12 or 34 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 14 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 ) .
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 ) , or r(u40 ) .
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 .
the amber force field may have base  base van der 
 waals interactions that are too strong .
alternatively , or in addition , nucleotide  water and water  water 
 interactions may not be balanced well .
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 .
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 .
 molecular 
 dynamics ( md ) simulations for rna tetramers r(aaaa ) , 
 r(caau ) , r(gacc ) , and r(uuuu ) are benchmarked against 1h1h 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
Methods
Results
Discussion