Improving Computational Predictions of Single-Stranded RNA

Subscriber access provided by University of Newcastle, Australia

Article

Improving Computational Predictions of Single-stranded RNA Tetramers with Revised #/# Torsional Parameters for the Amber Force Field David J. Wales, and Ilyas Yildirim J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b00819 • Publication Date (Web): 20 Mar 2017 Downloaded from http://pubs.acs.org on March 29, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Wales and Yildirim, p1

Improving Computational Predictions of Single-stranded RNA Tetramers with Revised α/γ Torsional Parameters for the Amber Force Field David J. Wales† and Ilyas Yildirim*‡§ †

Department of Chemistry, University of Cambridge, Cambridge, Cambridgeshire CB2 1EW, United Kingdom. ‡

Department of Chemistry and Biochemistry, Florida Atlantic University, Jupiter, FL 33458, USA. §

Scripps Research Institute, Jupiter, FL 33458, USA.

*

To whom correspondence should be addressed. Phone: (561) 799-8325. E-mail: [email protected]. ABSTRACT With current advancements in RNA based therapeutics, it is becoming crucial to utilize theoretical and computational methods to describe properly the physical properties of RNA molecules. NMR and X-ray crystallography are two powerful techniques for investigating structural properties. However, if the RNA molecules are complex or dynamic, these methods might not be adequate. For computational approaches, the quality of the force field will determine accuracy of our predictions. In this contribution, we revise the α/γ torsional parameters of RNA for amber force field using a model system representing an RNA dimer backbone. Combined with revised ߯ torsional parameters, previously shown to improve computational predictions, we benchmarked the revised force field on five single-stranded RNA (ssRNA) tetramers, three RNA dodecamer duplexes, and an RNA hairpin. A total of 60 µs of molecular dynamics (MD) simulations were run. We also employ the Discrete Path Sampling (DPS) approach to compare the predictions for the revised amber force field with those for amber10. Our results indicate that the unphysical states observed with amber10 in ssRNA MD simulations do not prevail with the ACS Paragon Plus Environment

`


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 32


revised amber force field. In line with NMR experimental observations, incorporation of the revised α/γ and ߯ torsional parameters leads to A-form-like conformational states as the most favorable ssRNA tetramer conformations. Furthermore, the revised force field maintains A-form geometry in regular RNA duplexes. Our revised amber force field for RNA should therefore improve structural and thermodynamic predictions for challenging RNA systems. INTRODUCTION The central dogma of biology states that genetic information stored in DNA is first transcribed into messenger RNA, and then translated into proteins. This workflow, however, does not include other roles of RNA in the cell such as retroviruses having RNA genomes,1-2 guide RNA (gRNA) directing the insertion and deletion of uridine bases,3 riboswitches regulating the mRNA level,4-5 retrotransposon amplifying specific genes via RNA,6 RNA interference inhibiting gene expression,7-8 ribozyme catalyzing specific biological activities,9-10 and RNA aptamers targeting different types of molecules.11 Furthermore, recent studies on anti-microRNA12 and CRISPR13 technologies open the door to next-gen therapeutics, which have the potential to treat cancer and genetic diseases. Hence, computational studies of RNA are of upmost importance to comprehend atomistic details of structure and dynamics of challenging RNA systems, which are not accessible with current experimental methods. However, this approach requires accurate computational models, capable of making physically meaningful predictions. Molecular mechanics (MM) methods are extensively used for simulations of RNA molecules because with the improvements of computational technology they can provide tools to run MD simulations comparable to biological time scales. MM methods describe the potential energy of a biomolecular system using a set of parameters that together constitute a force field, such as the amber model.14 Recently it has been shown that revision of ߯ torsional parameters for RNA is essential to reconcile computational predictions with optical melting and NMR results for RNA mononucleosides, single-stranded RNA (ssRNA) tetramers, and RNA molecules with modified residues.15-21 Even though revision of ߯ torsional parameters of RNA improved the predictions, MD simulations of ssRNA ACS Paragon Plus Environment

`

Page 3 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



tetramers still displayed some unphysical conformations, which could not be verified by NMR experiments.17, 22-24 Other issues in RNA tetraloops have been identified by others even after inclusion of revised ߯ torsional parameters.25-26 A recent study shows that revising the backbone torsions improve predictions for single-stranded RNA systems.21 Figure 1A shows an unphysical yet stable conformational state observed in MD simulations of ssRNA, GACC, with the amber10 force field. This intercalated conformation has been observed in other ssRNA tetramers as well, even after applying the revised ߯ torsional parameters. Thus, further improvements in the amber force field are required for a correct description of ssRNA tetramers.

Figure 1. (A) An unphysical intercalated state previously observed in MD simulations of the ssRNA GACC tetramer with the amber force field. (B) Model system used to revise the α/γ torsional parameters (see Figure S1 and Table S1 for details). To revise the α/γ torsional parameters of the amber force field for RNA, we use the model system shown in Figure 1B. Combined with the ߯ torsional parameters,16 the revised force field was benchmarked on five ssRNA tetramers, three RNA dodecamer duplexes, and an RNA hairpin. The revised force field maintains the A-form geometry in RNA dodecamer duplexes, and the Watson-Crick region of the RNA hairpin, while predicting the correct structural properties of ssRNA tetramers. The unphysical intercalated state (Figure 1A) no longer dominates for ssRNA tetramers with the revised

ACS Paragon Plus Environment

`


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 32


amber force field. These results are verified by discreet path sampling (DPS), where we analyze the complex energy landscape of the ssRNA GACC. Furthermore, MD simulations of the ssRNA tetramers with the revised force field display structures with approximately linear geometries. Nevertheless, in the MD simulations of ssRNA tetramers with adenosine residues (such as AAAA and CAAU), there are stable states, which cannot be verified by NMR. Our results therefore suggest that further force field revisions will be needed to resolve these issues. Overall, the revised α/γ description combined with revised ߯ torsional parameters represent RNA molecules better than the amber10 force field. METHODS Parameterization. The model system shown in Figure 1B was used to revise α and γ torsions for RNA. To create this model, an A-form dimer RNA was initially constructed. Later, the bases were replaced with –CH3 groups (Figure S1). RESP charges for the model system were calculated as previously described (Table S1).27 In the RESP protocol, only the charges for –CH3 groups were fit, and the original resp charges for the rest of the atoms were kept (see Table S1). Ab Initio Potential Energy Surface (PES) Scan of α/γ. The α and γ torsions were revised on the basis of 2D potential energy surface (PES) scans on two conformational states of the model system: i) Both sugars in C3ʹ-endo (C3ʹ-C3ʹ), and ii) 5ʹ- sugar in C3ʹ-endo and 3ʹ- sugar in C2ʹ-endo (C3ʹ-C2ʹ). In the model, the other torsions were kept in A-form dihedral values using torsional restraints. For each conformation, α and γ were rotated with increments of 10°, yielding 2×(36×36)=2592 data points for the α/γ revision. The γ torsional parameters for the 5ʹ-terminal end were also revised, where the 5ʹ-OH group was rotated in increments of 10°, yielding 36 data points. For each conformation in the PES scan, the structures were first optimized at the HF/6-31G* level of theory followed by a single point energy calculation at the MP2/6-31G* level.28-30 The final energy of each conformation, EQM , was then used in force field fitting (Eq 1). All the QM calculations were done with Gaussian 03.31 Force Field Fitting of α/γ Torsions. A similar protocol to that previously described was followed.16 (noAG ) The molecular mechanics (MM) energies, EMM where AG stands for α/γ, of each conformation were


`

Page 5 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



calculated by restraining the dihedral angles with a force constant of 100,000 kcal/(mol-rad2) to the values observed in the optimized QM geometries. In the calculations, the α/γ torsional parameters were (noAG ) (noAG ) set to zero yielding EMM (Eq 1). The energy difference, EQM − EMM , represents the potential energy

due to the α/γ torsions. ( noAG ) EQM − EMM = EAG

(1)

The data points of Eq 1 were fit by linear least-squares32-33 to the Fourier series shown in Eq 2 by utilizing trigonometric identities: 4

[

(

)]

[

(

fit (α , γ ) = ∑Vnα 1 + cos nα − Pnα + Vnγ 1 + cos nγ − Pnγ E AG

)]

(2)

n =1

Here, α and γ are the dihedral angles defined as O3A-P-O5B-C5B and O5B-C5B-C4B-C3B, respectively, (see Figure S1), n is the periodicity, Vnα and Vnγ are the potential energy barriers for the α and γ torsions, and Pnα and Pnγ are the phase angles for periodicity n of α and γ torsions. A similar approach was followed to revise the γ torsional parameters of the 5ʹ-terminal end. The revised α/γ parameters are collected in Table 1. Comparison of α and γ profiles between amber10 and the revised α/γ force fields is shown in Figure 2. Comparisons between the 2D QM and the MM profiles with the α/γ revision are shown in Figure 3. The 2D MM profiles with the amber10 force field are in Figure S2. The results show that with the revised α/γ torsional parameter set the RMSD of 2D energy profiles with respect to QM for C3ʹ-C2ʹ and C3ʹ-C3ʹ are 1.02 and 0.85 kcal/mol, respectively, while they are 1.74 and 2.10 kcal/mol, respectively, with the amber10 force field (Figures 3 and S2).


`


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 32


Table 1. Revised α/γ torsional parameters for RNA. Torsion α

γ

γ-terminalb

a

Torsional

potential

energy

na Via (kcal/mol) 1 1.75635 2 1.71456 3 0.45866 4 0.35268 1 0.45773 2 1.05380 3 1.40341 4 0.02753 1 2.28202 2 0.989337 3 1.10808 4 0.0200738 in the AMBER

Pia (°) 42.5639 35.9304 352.3032 317.2117 204.0770 243.0240 329.1598 90.6059 161.139 234.231 346.035 294.774 force field

is

calculated

as

4

EMM ,tor (φ ) = ∑Vi [1 + cos(nφ − Pi )] where Vi is the potential energy barrier, φ is the dihedral angle, Pi is n =1

the phase, and n is the periodicity. b For the 5ʹ-terminal, separate fitting was performed for the γ-torsion (see Methods for details).

Figure 2. Energy profiles for the α and γ torsions with the revised (red) and amber10 (blue) force fields. Dashed red profile represents the γ torsion of the 5ʹ-terminal end.


`

Page 7 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Preparation of the systems. Five ssRNA tetramers (GACC, CCCC, UUUU, CAAU, AAAA), three RNA dodecamer duplexes, and one RNA hairpin were investigated (see Table 2 for details). The initial conformations of the ssRNA structures were either A-form or intercalated, which were created with either the nucgen module of AMBER12 or by homology modeling by replacing the RNA residues. The NMR and X-ray structures of two of the dodecamer duplexes (PDB ID 1AL5,34 and 1SDR35) were used to study these RNA duplexes. To investigate the B-form → A-form transformation, we utilized the Xray structure of a dodecamer (PDB ID 1BNA36). Homology modeling was used to build the RNA version of 1BNA. Finally, the NMR structure of PDB ID 2KOC37 was used to investigate the properties of a hairpin. The leap module of the AMBER MD package14 was used to neutralize all the systems with Na+ ions.38 After neutralization, five Na+ and Cl- ions were added to each system. Structures were solvated with TIP3P39 water molecules in a truncated octahedral box (see Table 2 for details). Previously, the χ torsional revision has been shown to improve predictions.15-17, 19, 23-24, 40-42 Thus, the amber force field43 with revised χ16 and α/γ torsional parameters (Table 1) was used in the MD simulations. For testing purposes, ssRNA GACC was also analyzed with the amber10 force field (see Table 2 for details).


`


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 32


Figure 3. 2D PES scans of α (x-axis) and γ (y-axis). (A) and (C) represent the QM profiles while (B) and (D) represent the MM profiles calculated with the revised α/γ torsional parameters. (A) and (B) are calculated for the model system (Figure 1B) with sugars in C3ʹ-endo and C2ʹ-endo conformations, while (C) and (D) are calculated for the model system with both sugars in C3ʹ-endo conformations. Energy scales are in kcal/mol. The RMSD of (B) and (D) with respect to (A) and (C) is 1.02 and 0.85 kcal/mol, respectively. See Figure S2 for a comparison with the amber10 predictions.


`

Page 9 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Table 2. RNA structures investigated. System GACC GACC GACC CCCC UUUU CAAU AAAA 1AL5e 1SDRf 1BNAg 2KOCh a

Sequence 5ʹ-GACC-3ʹ 5ʹ-GACC-3ʹ 5ʹ-GACC-3ʹ 5ʹ-CCCC-3ʹ 5ʹ-UUUU-3ʹ 5ʹ-CAAU-3ʹ 5ʹ-AAAA-3ʹ 5ʹ-CGCAAAUUUGCG-3ʹ 3ʹ-GCGUUUAAACGC-5ʹ 5ʹ-UAAGGAGGUGAU-3ʹ 3ʹ-AUUCCUCCACUA-5ʹ 5ʹ-CGCGAAUUCGCG-3ʹ 3ʹ-GCGCUUAAGCGC-5ʹ 5ʹ-GGCACUUCGGUGCC-3ʹ

Volumed (×103 Å3) 63.0 63.0 63.0 63.0 63.0 63.0 63.0 214.2

Initial Structure Intercalated Intercalated A-form Intercalated Intercalated Intercalated Intercalated A-form

Force Field amber99+߯+α/γ amber10 amber10 amber99+߯+α/γ amber99+߯+α/γ amber99+߯+α/γ amber99+߯+α/γ amber99+߯+α/γ

Na /Cl /H2O 8 / 5 / 2037 8 / 5 / 2037 8 / 5 / 2037 8 / 5 / 2037 8 / 5 / 2037 8 / 5 / 2037 8 / 5 / 2037 27 / 5 / 6848

MD Timec 10 µs 1 µs 1 µs 10 µs 10 µs 10 µs 10 µs 1 µs

A-form

amber99+߯+α/γ

27 / 5 / 6367

1 µs

199.6

B-form

amber99+߯+α/γ

27 / 5 / 6417

1.5 µs

201.1

Hairpin

amber99+߯+α/γ

18 / 5 / 4212

3.9 µs

131.6

a

amber99+߯+α/γ represents the amber99 force field with revised

+

߯16

-

b

and α/γ torsional parameters.

b

Total number of Na+, Cl-, and TIP3P water molecules in each system. c A 2 fs time step was used in all the MD simulations. d Volume of each unit cell after equilibration. e NMR structure with PDB accession code 1AL534 was used to build the initial conformation. f Crystal structure with PDB accession code 1SDR35 was used to build the initial conformation. g Crystal structure with PDB accession code 1BNA36 was used to homology model the initial B-form RNA conformation.

h

NMR structure with PDB

accession code 2KOC37 was used to build the initial conformation, which is a UUCG hairpin RNA. Molecular dynamics simulations. The structures were minimized in two steps as described before.27, 40 In the first step, heavy atoms were restrained by harmonic restraints with a force constant of 10 kcal/(mol-Å2), while in the second step no restraints were applied. Both steepest-descent and conjugategradient minimization were employed. The nonbonded cutoff during the minimization was 8.0 Å (see Table S2 for sample input files). After minimization, each system was equilibrated in two steps (see Table S2 for sample input files).15, 17, 40

In the first step, the temperature was gradually increased from 0 to 300 K within 20 ps under

constant volume dynamics (NVT) while restraining all the heavy atoms. In the second step, another 2 ns of MD at 300 K was run under constant pressure dynamics (NPT), while still imposing restraints on the heavy atoms (see Table S2 for details).


`


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 32


After equilibration, NPT dynamics with isotropic positional scaling was used in all the production runs. The reference pressure was set to 1 atm with a pressure relaxation time of 2 ps. Bonds involving hydrogen atoms were constrained with SHAKE.44 In all the MD simulations an atom-based long-range hard cutoff of 8.0 Å was used. MD simulation times for ssRNA systems were 10 µs. The MD simulations of ssRNA GACC using the amber10 force field were run for 1 µs, and those for the RNA dodecamer duplexes were run for 1 or 1.5 µs. The MD simulation of the RNA hairpin was close to 4 µs long (see Table 2 for details). The GPU version of pmemd (pmemd.cuda)14 was used to run all these MD simulations. Discrete path sampling. The discrete path sampling (DPS)45-46 method provides an approach based on geometry optimization to scan configuration space efficiently, yielding global and local minimum conformational states of RNA systems as well as transformation pathways that might not be accessible with regular MD simulations. DPS has been applied successfully to 1×1 UU internal loops in RNA CUG repeat expansions.47 In the present work, DPS was applied to GACC to compare the predictions of the revised ߯ and α/γ parameter set with the predictions of amber10. For the DPS calculations initial conformations utilized frames extracted from the MD trajectories of GACC that included both A-form and intercalated states. An initial database was built by minimizing the starting conformations with a modified version of the LBFGS algorithm.48 The convergence criterion for the root mean square gradient was 10−6 kcal mol−1. Connection attempts were made between the minima in a pairwise fashion. In the end, a disconnectivity graph49-51 analysis was performed to display the stable conformational states. Further refinement of the database was performed with the UNTRAP scheme.52 For all these calculations, the OPTIM and PATHSAMPLE software was used (see Table S3 for sample input files). The final stationary point database for the revised α/γ parameter set contains ~95K minima and ~151K transition states, while the database for the amber10 force field contains ~140K minima and


`

Page 11 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



~194K transition states for GACC. Free energies were estimated from the database by employing the harmonic superposition approximation.53 Analysis. The ptraj module of AMBER1214 was used for dihedral, root-mean-square deviation (rmsd), and cluster analyses. X3DNA54 was used to extract global structural properties of RNA duplexes from the MD trajectories. DSSR55 was used to extract base stacking properties in ssRNA structures from the MD trajectories. Disconnectivity graphs49-50 were plotted to visualize the free energy landscapes for GACC.51, 56 The disconnectivity graphs depicted in this work were constructed from the stationary point databases using the disconnectionDPS code.57 RESULTS and DISCUSSIONS Comparison of energy profiles for the revised α/γ parameters and amber10. The α/γ torsional parameters in amber10 were previously revised with a model system of one deoxyribose with a phosphate group at the end.58 Here, the model system shown in Figure 1B represents the RNA backbone more realistically. Figure 2 compares the energy profiles for the revised α and γ parameters with those in amber10. There are several differences between the parameter sets that can be summarized as follows: 1) In amber10, the energy profile of α has two minima (at ~70° and ~270°) with almost identical energies. Furthermore, the energy barriers separating these two minima are between 2 and 3 kcal/mol (at ~0° and ~180°). In contrast, the energy profile of the revised α parameter set has a global minimum at ~300° and a local minimum at ~120°. Furthermore, there are two pathways in the revised α torsional profile with maxima at ~0° and ~180° with energy barriers of ~6.5 and ~3 kcal/mol, respectively. This result basically implies that the revised α parameter set will prefer the pathway at ~180° while amber10 supports both pathways at ~0° and ~180° for transformations. 2) Both amber10 and the revised γ torsional profiles exhibit two local minima at ~180° and ~300° and the same global minimum at ~50°. However, the energy barriers are different (see Figure 2). The differences seen in the α and γ torsional profiles will have distinct effects on how the RNA backbone behaves in MD simulations, as we have observed in the current work. ACS Paragon Plus Environment

`


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 32


DPS analysis shows that the revised (amber99+߯+α/γ) and amber10 force fields prefer the (physical) A-form and the (unphysical) intercalated states, respectively, as the global minimum structures for the GACC. Figure 4 displays the free energy disconnectivity graphs of GACC using the revised and amber10 force fields. According to the DPS analysis, the revised ߯ and α/γ torsional parameter set predicts the A-form conformation to be the global minimum structure (structure # 1, Figure 5). On the other hand, amber10 predicts the intercalated state to be the global minimum (structure # 3*, Figure S3) as previously seen by others.22 With amber10, the intercalated state is ~2 kcal/mol more stable than Aform, while in the revised force field A-form is ~ 1 kcal/mol more stable than the intercalated state in line with experimental results (Figure 4). According to NMR analysis, GACC prefers A-form-like conformations.17 Thus, the preference for the unphysical intercalated state with amber10 is eliminated by the revised α/γ torsional parameters.

Figure 4. Free energy disconnectivity graphs of single-stranded GACC RNA at 300 K using (A) revised ߯ + α/γ torsional parameters, and (B) amber10 force field. No regrouping threshold was used in the analysis. For definitions of the structures highlighted in the figure, see Figure 5. Structure 3* is shown in ACS Paragon Plus Environment

`

Page 13 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



Figure S3. Note that amber10 overstabilizes 3*, which is known to be an unphysical state. On the other hand, the revised ߯ + α/γ torsional parameter set improves predictions, with physical A-form states stabilized over 3*.

Figure 5. Selected conformations observed in the ssRNA MD simulations (see Table 3 for descriptions). For illustration purposes, only GACC structures are shown here (for details see Figures S5-S9). MD simulations of GACC with the amber10 force field verify the DPS results. The quality of the amber10 force field was tested on GACC by running two independent 1 µs long MD simulations. The initial structures were the A-form and intercalated (see Figure 1A, and structure # 3 in Figure 5). The RMSD with respect to the A-form and intercalated states were calculated and are shown in Figure S4. Results indicate that the intercalated state is the most dominant for GACC with the amber10 force field (Figure S4). Within 100 ns, the A-form structure transforms into structure 3* (Figure S3), stays there for over 600 ns, and transforms into the intercalated state shown in Figure 1A (Figure S4A, and Movie S1). In the MD simulation starting with the intercalated conformation (Figure 1A), GACC transforms into the A-form briefly (for 100 ns), and then transforms back into the intercalated state (Figure S4B). Similar results were previously observed for single-stranded RNA.17,


`

22, 42

The MD results verify the


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 32


DPS results discussed above for amber10, where the unphysical intercalated states dominate the MD simulations. MD simulations of GACC with the revised ߯+α/γ torsional parameters predict A-form structures to be the dominant states. While the ߯ torsional parameters for RNA used in the current work improved predictions of NMR and optical melting results for different RNA systems,15-17 other problems remained, including intercalated states being a stable conformation in MD simulations.17,20-22 Thus, predicting the favorable conformations of single-stranded RNA correctly will increase trust in RNA force fields, and further their application to challenging systems. To further test the revised ߯ and α/γ torsional parameters, a 10 μs MD simulation was run for GACC using intercalated state as starting structure (Figure 6A). The unphysical intercalated state transformed into A-form conformations within 250 ns (Figure 6A). The MD trajectory exhibits other conformational states, but the A-form conformations are the dominant structures (see Figure S5 for the structures). Cluster analysis showed that GACC prefers the major and minor A-form conformations 74.4% and 6.9%, respectively. These structures are consistent with the experimental NMR spectra. In contrast, the intercalated state is present for only 1.6% of the simulation (Table 3). Unstacked G1 and C4 states are observed for 3.8% and 2.5%, respectively, which is reasonable because of their flexibility compared to the internal residues A2 and C3. The bulged out A2 and C3 states (structures 4 and 7 in Figure 5) account for 2.1% and 4.3%, respectively. Similar results were observed after doing stacking analyses on the MD trajectory (Table 4). These results basically indicate that with the revised force field, the simulation predicts over 85% “extended” A-form conformations for GACC, while bulged out and intercalated states account for less than 10% (Tables 3 and 4).


`

Page 15 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48



Figure 6. RMSD analysis of single-stranded (A) GACC, (B) CCCC, (C) UUUU, (D) CAAU, and (E) AAAA RNA. A color code is used to highlight the conformations observed in the MD trajectories (see Table 3 and Figures 5, S5-S9).


`


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 32


Tabel 3. Observed conformations in single-stranded RNA simulations with the revised RNA force field.a A color code is used in Figure 6 to highlight the conformations seen in the MD trajectories. Clustering is performed with RMSD analysis. Percentage (%) Structure Description

Color

GACC

CCCC

UUUU

CAAU

AAAA

1

Major A-form Conformation

Black

74.4

78.1

39.5

37.4

31.8

2

Minor A-form Conformation

Red

6.9

5.2

36.0

25.1

12.5

3

Intercalated structure

Blue

1.6

2.9

1.4

18.0

11.5

4

1 unstacked

Green

3.8

2.6

2.8

3.3

9.7

5

1-2 and 3-4 stacked

Yellow

0.9

1.2

4.4

3.5

1.5

6

2 bulged out

Orange

2.1

2.8

-

3.2

0.5

7

4 unstacked

Grey

2.5

3.5

8.4

1.6

1.1

8

3 bulged out

Brown

4.3

-

0.9

1.1

3.8

9

Linear conformations (seen in

Magenta

-

-

-

3.0

23.6

Cyan

-

-

-

-

0.2

Lime

-

-

3.0

-

-

AAAA and CAAU) 10

Only 2-3 stacked (seen in AAAA)

11

Other conformations seen in UUUU

a

Structures 1-8 are illustrated in Figure 5. Structures 9-11 are illustrated in S7-S9.


`

Page 17 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48


Wales and Yildirim, p17 a

Table 4. Percentage (%) of observed base-base stackings in single-stranded RNA simulations. ssRNA

1-2-3-4

1-2-3

2-3-4

2-3-1-4

1-2 & 3-4

1-2-4

1-3-4

1-2

2-3

3-4

1-3

2-4

1-4

1-3-2

GACC

74.4

7.7

5.9

1.0

2.6

4.1

1.6

0.5

0.5

0.2

0.5

0.2

0.0

0.4

CCCC

76.0

6.8

4.7

2.1

5.1

0.0

2.3

0.6

0.5

0.5

0.5

0.0

0.1

0.2

UUUU

25.2

17.1

8.2

0.0

19.1

0.4

0.4

11.4

6.5

6.3

0.5

0.5

0.2

0.3

CAAU

48.1

7.5

10.2

13.8

7.0

0.8

2.5

1.2

2.1

1.5

0.5

0.2

0.1

3.6

AAAA

50.6

5.6

14.0

10.8

6.1

3.4

0.5

0.9

1.5

2.4

0.3

0.2

0.0

2.9

Comparable Structures according to Table 3b Structure a

1, 2

7

4

3

5

8

6

10

Base-base stacking calculations are performed with DSSR.55 Numbers displayed in the stacking analysis correspond to residue numbers, e.g.

1-2-3-4 correspond to R1-R2-R3-R4 where R stands for the residue name. Note also that 1-2-3-4 represent all the residues stacked in a linear or A-form conformations. b The numbers displayed here correspond to the structures shown in Table 3.

`



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 32


Similar to GACC, MD simulations of CCCC with the revised ߯+α/γ torsional parameters predict A-form structures to be the dominant states. Figure 6B displays the RMSD analysis of ssRNA CCCC. Similar to GACC, the dominant structures observed in the MD trajectory of CCCC are A-form-like. The initial intercalated state transforms into an A-form structure within 250 ns (Figure 6B). Cluster analysis showed that CCCC prefers major and minor A-form conformational states 78.1% and 5.2% of the time, respectively, while the intercalated state is sampled for only 2.9% of the MD trajectory (see Figures 5 and 6B). Unstacked C1 and C4 states are observed for 2.6% and 3.5%, respectively, while bulged out C2 is observed for 2.8% (Table 4 and Figure 5). No bulged out C3 state was seen in the MD trajectory of CCCC. Similar results were obtained in stacking analysis (Table 4). The computational predictions are consistent with the NMR results of CCCC where major and minor A-form conformations are observed.23 Similar to the GACC results, the revised ߯+α/γ torsional parameters predict over 90% Aform-like conformations for ssRNA CCCC in a 10 µs long MD simulation (see Tables 3 and 4). Compared to GACC and CCCC, UUUU displays flexional behavior, even though the dominant states are still A-form-like. Figure 6C shows the RMSD analysis of UUUU. The initial intercalated structure immediately transformed into an “extended” conformation and was almost never sampled again (Figure 6C). Although the MD trajectory displays flexional behavior for UUUU, cluster analysis shows that it sampled major and minor A-form states for 39.5% and 36.0% of the time, respectively (Figure 6C and Table 3). The unstacked U1 and U4 states are sampled for a total of 11.2% of the time while only bulged out U3 was observed less than 1% (Table 3). The 1-2&3-4 stacked states are observed for 4.4%; the highest value within the ssRNA systems studied (Table 3). This result implies that bending in the middle of UUUU is predicted to be much easier compared to the other ssRNA tetramers we have considered. The stacking analysis is particularly important because it verifies the flexional behavior seen in UUUU. Analysis of the MD trajectory reveals UUUU to have 25.2% 1-2-3-4 stacking (Table 4). Furthermore, analysis indicates that 1-2-3, 2-3-4, and 1-2&3-4 stacking is observed for 17.1%, 8.2%, and 19.1% of the time, respectively (Table 4). Moreover, 1-2 and 2-3 stacking are observed for 11.4% and 6.5%, `


Page 19 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



respectively (Table 4). All results combined for UUUU indicate that it prefers “extended” conformations; yet, the flexional behavior happens in such a way that it transforms rapidly from one Aform-like state to another (Figure 6C). The predictions are particularly interesting because similar results were observed in the NMR analysis of UUUU.42 Only one NOE was observed for UUUU, indicating solvent exposure for the residues.42 NMR results for UUUU also indicate that the structure is either unstacked or disordered. Furthermore, SAXS and smFRET experiments indicate poly(U) having less backbone flexibility compared to poly(T).59 Thus, the predictions with the revised ߯+α/γ torsional parameters are in line with the NMR results, in the sense that they predict a combination of “extended” and disordered-linear RNA states. The dominant structures in CAAU are A-form-like, but the MD trajectory has the highest percentage of intercalated states within the set of ssRNA structures studied here. Figure 6D shows the RMSD analysis of ssRNA CAAU, which indicates that the major and minor conformational states observed in the MD trajectory account for 37.4% and 25.1% of the time, respectively (Table 3). The initial intercalated state transformed into the A-form conformation within 50 ns; however, the intercalated states were observed for 18% of the MD trajectory (Table 3). The other structures observed in the trajectory each appeared less than 4% of the time (see Table 3). Similar results were observed in stacking analyses (Table 4). According to the NMR results, CAAU mostly prefers A-form-like conformations (both major and minor A-form structures). Previous MD simulations of CAAU with the amber10 force field displayed a high percentage of intercalated state.42 However, with the revised ߯+α/γ torsional parameters, the predictions are more in line with the experimental results. Nevertheless, there might be other issues with the RNA force field causing CAAU to sample the intercalated state significantly. As for CAAU, AAAA has a relatively high population of intercalated states in the MD trajectory. In the MD simulation of AAAA, the initial intercalated structure transforms into the A-form within 400 ns (Figure 6E). According to the cluster analyses, AAAA samples the major and minor A-form conformers for 31.8% and 12.5% of the time, respectively, while the intercalated state is sampled for 11.5% (Table `



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 32


3). Furthermore, the MD trajectory exhibits 23.6% of conformations in “extended” but non-A-form conformations (Figure 6E and Table 3). Even though the MD trajectory of AAAA displays 55.1% combined A-form-like structures (major, minor, and residues 1 and 4 unstacked structures, see Table 3) our results indicate further issues with the RNA force field, as apparent from CAAU results (see above). The “extended” non-A-form conformations seen in the MD trajectory of AAAA have adenosine residues in syn conformational states. According to the NMR results, AAAA has A-form-like conformations (structures 1 and 2 shown in Figure 5) where all the adenosines are in anticonformational states.42 Seeing unphysical states in ssRNA tetramers of CAAU and AAAA implies that the force field parameters of adenosine might need improvement. Recently it was proposed that the Lennard-Jones parameters of the adenosine base atoms required revisions to explain some experimental observations.60 Furthermore, the water model used in the MD studies might also affect the predictions. For example, Cheatham and co-workers recently benchmarked a revised force field for the OPC water model on GACC and CCCC.22 Even though the revised α/γ torsional parameter set lowers the percentage of intercalated structures seen in MD simulations, further revisions are required, especially on adenosine residues. Revised ߯ and α/γ parameters preserve the A-form structure of 1AL5 dodecamer duplex (Table 2). To investigate the effects of revised ߯ and α/γ parameters on RNA duplexes, the revised force field was benchmarked on several RNA systems (see Table 2). 1AL5, r(CGCAAAUUUGCG)2, is an NMR structure. The RMSD of the MD trajectory with respect to A-form is shown in Figure S10A. Overlap of the A-form, NMR, and average structure calculated from the MD trajectory are shown in Figure 7A. Results indicate that the RNA duplex is A-form and stable for at least 1 µs. with the revised ߯ and α/γ parameter set. Properties of base-pair steps extracted from the MD trajectory were also analyzed (Table S4). While the MD calculation predicts a minor groove width (mgw) between 15.3 and 15.7 Å, the major groove width (Mgw) is between 15.8 and 17.1 Å (Table S4). These results are in line with the structural properties of the NMR structure, which has a mgw and Mgw of 14.5 and 17.3 Å, respectively `


Page 21 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



(Table S4). The pure A-form version of 1AL5 was also created and compared to both the NMR and predicted average structures (Figure 7A). Pure A-form structure has a mgw and Mgw of 16.9 and 12.3 Å, respectively, which differ from the NMR structure (Table S4). In solution, biological systems are dynamic, and water will interact with the functional groups of RNA residues. This interaction will, therefore, change the groove widths, especially the Mgw, as observed in the analyses. Thus, one has to be careful when comparing X-ray structures to average structures extracted from MD simulations.

Figure 7. Overlap of MD average (Blue), A-form (Red), and experimental X-RAY/NMR (Orange) structures of (A) 1AL5, (B) 1SDR, and (C) 1BNA. The average structures are predicted from MD trajectories, and overlap well with A-form conformations (see Tables S4-S6 for details). In (C), the initial structure of 1BNA is B-form, and transforms to A-form in the MD simulation (see Figure S10). As for the 1AL5 results, revised ߯ and α/γ parameters preserve the A-form structure of the 1SDR dodecamer duplex (Table 2). We also benchmarked the X-ray structure of the RNA duplex 1SDR, r(UAAGGAGGUGAU)2, which has a resolution of 2.6 Å. The difference between 1AL5 and 1SDR is the order of the base-pairs in the duplex. In 1AL5, six AU are flanked by three GC base-pairs. On the other hand, 1SDR has six GC base-pairs flanked by three AU pairs (Table 2). Figure S10B shows the RMSD for the MD trajectory with respect to the A-form. Furthermore, the overlaps of the A-form, X-

`



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 32


ray, and average structure calculated from the MD trajectory are shown in Figure 7B. As for 1AL5, the structure of 1SDR stays in the A-form throughout. Base-pair step, mgw, and Mgw analyses also display A-form properties (Table S5). Similar to the 1AL5 results, MD simulation of 1SDR displays mgw between 15.3 and 15.6 Å, which is comparable to the X-ray value of 15.7 Å (Table S5). On the other hand, the simulation produces Mgw between 17.0 and 17.8 Å while the X-ray structure has Mgw around 12.3 Å (Table S5). As discussed above, this observation is expected because X-ray represent a frozen packed structure, while MD simulations will include dynamic effects. The MD calculated Mgw values are comparable to the results discussed above for 1AL5. B-form dodecamer RNA transforms into A-form in a 1.5 µs long MD simulation. As a different test, a Bform dodecamer duplex, 1BNA (Table 2), was transformed into RNA by replacing deoxyribose with ribose sugars and T with U. The B-form RNA was transformed into an A-form-like state immediately in the MD simulation. However, the complete transformation does not happen until 1µs (Figures 7C and S10C). The B→A transformation was directed with the transformation of the sugar puckers from C2ʹendo to C3ʹ-endo. At the end, all the sugars were C3ʹ-endo and the base-pair step parameters agreed with A-form properties (Table S6). The mgw and Mgw values predicted for the 1BNA structure were ~15.5 and 17.5 Å, respectively, which are again consistent with typical A-form (Table S6). MD simulation of UUCG hairpin RNA with the revised ߯ and α/γ parameters displays A-form properties. The revised amber force field was also tested against the NMR hairpin structure of r(GGCACUUCGGUGCC), which contains a UUCG loop and a stem region with five Watson-Crick base pairs (PDB ID 2KOC, Figure 8A). A 3.9 µs MD simulation was run to investigate the dynamic properties of the stem and loop regions. RMSD analysis shows that the stem stays in A-form (Figure S11A). The RMSD of the average structure of the stem region calculated from the MD trajectory with respect to the NMR geometry is 1.0 Å (Figure 8). The hairpin residues of U6, U7, and C8 were not as dynamic as G9, and stayed in their initial conformations on the MD time scale (Figure S11B,C,D). G9, however, unstacked from the initial state around 0.75 µs, and stayed in this unstacked state for the rest of `


Page 23 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



the simulation (Figure 8B and S11E). This observation might imply other issues with the RNA force field, such as the water model, the van der Waals parameters for RNA residues, and parameters describing ions, or the MD simulation time is not enough to fully investigate the hairpin region.

Figure 8. (A) Initial, and (B) final structures of 2KOC observed in MD simulation. Note that G9 is unstacked from the hairpin region after around 0.75 µs (Figure S11). The stem part and the other residues in the hairpin region maintain their initial conformation (see also Figure S11). SUMMARY AND CONCLUSIONS We have revised the α/γ parameters of the amber force field for RNA, and benchmarked the force field against five single-stranded RNA tetramers, three RNA dodecamer duplexes, and one RNA hairpin. Combined with revised ߯ parameters for RNA described previously,16 the new α/γ parameters improve MD simulations of ssRNA by greatly reducing the unphysical intercalated state and predicting the Aform-like conformational states to be dominant. The simulations also predict the ssRNA structures to be dynamic, sampling conformational states such as stacked, unstacked, and bulged out RNA residues. `



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 32


Such states are not excluded by the NMR experiments published thus far. Moreover, application of the DPS approach for GACC make predictions similar to the MD simulations, and suggests the potential application of DPS to study challenging RNA systems, which have complex energy landscapes. Benchmarks against RNA dodecamer duplexes and hairpin reveal that the revised ߯ and α/γ parameters keep RNA duplexes in the A-form, which is the native conformational state of RNA. Furthermore, starting with a B-form RNA dodecamer duplex, the combination of revised ߯ and α/γ parameters transforms the duplex into an A-form conformation through transformations of the sugar puckers from C2ʹ-endo to C3ʹ-endo. Finally, even though inclusion of revised ߯ and α/γ parameters in the RNA force fields improves predictions there are details suggesting other issues with the force field for RNA. The possible improvements still required can be summarized as i) revision of van der Waals parameters for adenosine residue, ii) optimization of ion and water models, and iii) revision of other backbone angles. Overall, our results show that the combination of revised ߯ and α/γ parameters significantly improves our structural and thermodynamic predictions for ssRNA. ASSOCIATED CONTENT Supporting Information. Movie description; atom types/names/charges of the model system; 2D PES scans, MD trajectories, and conformations observed with amber10; structures observed in ssRNA simulations; RMSD analyses; sample AMBER/OPTIM/PATHSAMPLE/disconnectionDPS input files; helical properties of RNA duplexes. The Supporting Information is available free of charge on the ACS Publication website. ACKNOWLEDGMENT We thank Prof. Douglas H. Turner for helpful discussions and editing the manuscript. Computations were done in the Theory Group Computing Clusters at the University of Cambridge. This work was supported by the Department of Chemistry and Biochemistry, Florida Atlantic University, (IY) and the EPSRC (grant reference N035003) (DJW). `


Page 25 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



REFERENCES 1. Bishop, J. M. Cellular Oncogenes and Retroviruses. Annu. Rev. Biochem. 1983, 52, 301-354. 2. Gifford, R.; Tristem, M. The Evolution, Distribution and Diversity of Endogenous Retroviruses. Virus Genes 2003, 26, 291-316. 3. Kable, M. L.; Seiwert, S. D.; Heidmann, S.; Stuart, K. RNA Editing: A Mechanism for gRNASpecified Uridylate Insertion into Precursor mRNA. Science 1996, 273, 1189-1195. 4. Haller, A.; Souliere, M. F.; Micura, R. The Dynamic Nature of RNA as Key to Understanding Riboswitch Mechanisms. Acc. Chem. Res. 2011, 44, 1339-1348. 5. Smith, A. M.; Fuchs, R. T.; Grundy, F. J.; Henkin, T. M. Riboswitch RNAs Regulation of Gene Expression by Direct Monitoring of a Physiological Signal. RNA Biol. 2010, 7, 104-110. 6. Sabot, F.; Schulman, A. H. Parasitism and the Retrotransposon Life Cycle in Plants: A Hitchhiker's Guide to the Genome. Heredity 2006, 97, 381-388. 7. Hannon, G. J. RNA Interference. Nature 2002, 418, 244-251. 8. Kole, R.; Krainer, A. R.; Altman, S. RNA Therapeutics: Beyond RNA Interference and Antisense Oligonucleotides. Nat. Rev. Drug Discovery 2012, 11, 125-140. 9. Birikh, K. R.; Heaton, P. A.; Eckstein, F. The Structure, Function and Application of the Hammerhead Ribozyme. Eur. J. Biochem. 1997, 245, 1-16. 10. Frank, D. N.; Pace, N. R. Ribonuclease P: Unity and Diversity in a tRNA Processing Ribozyme. Annu. Rev. Biochem. 1998, 67, 153-180. 11. Patel, D. J.; Suri, A. K.; Jiang, F.; Jiang, L. C.; Fan, P.; Kumar, R. A.; Nonin, S. Structure, Recognition and Adaptive Binding in RNA Aptamer Complexes. J. Mol. Biol. 1997, 272, 645-664. 12. Lu, Y. J.; Xiao, J. N.; Lin, H. X.; Bai, Y. L.; Luo, X. B.; Wang, Z. G.; Yang, B. F. A Single AntimicroRNA Antisense Oligodeoxyribo-Nucleotide (AMO) Targeting Multiple microRNAs Offers an Improved Approach for microRNA Interference. Nucleic Acids Res. 2009, 37, e24.

`



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 32


13. Cong, L.; Ran, F. A.; Cox, D.; Lin, S. L.; Barretto, R.; Habib, N.; Hsu, P. D.; Wu, X. B.; Jiang, W. Y.; Marraffini, L. A., et al. Multiplex Genome Engineering Using Crispr/Cas Systems. Science 2013, 339, 819-823. 14. Case, D. A.; Darden, T. A.; Cheatham, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M., et al. Amber 12, University of California, San Francisco: San Francisco, CA, 2012. 15. Yildirim, I.; Kennedy, S. D.; Stern, H. A.; Hart, J. M.; Kierzek, R.; Turner, D. H. Revision of Amber Torsional Parameters for RNA Improves Free Energy Predictions for Tetramer Duplexes with GC and iGiC Base Pairs. J. Chem. Theory Comput. 2012, 8, 172-181. 16. Yildirim, I.; Stern, H. A.; Kennedy, S. D.; Tubbs, J. D.; Turner, D. H. Reparameterization of RNA Χ Torsion Parameters for the Amber Force Field and Comparison to NMR Spectra for Cytidine and Uridine. J. Chem. Theory Comput. 2010, 6, 1520-1531. 17. Yildirim, I.; Stern, H. A.; Tubbs, J. D.; Kennedy, S. D.; Turner, D. H. Benchmarking Amber Force Fields for RNA: Comparisons to NMR Spectra for Single-Stranded R(GACC) Are Improved by Revised Χ Torsions. J. Phys. Chem. B 2011, 115, 9261–9270. 18. Zgarbova, M.; Otyepka, M.; Sponer, J.; Mladek, A.; Banas, P.; Cheatham, T. E.; Jurecka, P. Refinement of the Cornell Et Al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles. J. Chem. Theory Comput. 2011, 7, 2886-2902. 19. Deb, I.; Sarzynska, J.; Nilsson, L.; Lahiri, A. Conformational Preferences of Modified Uridines: Comparison of Amber Derived Force Fields. J. Chem. Inf. Model. 2014, 54, 1129-1142. 20. Bergonzo, C.; Henriksen, N. M.; Roe, D. R.; Cheatham, T. E. Highly Sampled Tetranucleotide and Tetraloop Motifs Enable Evaluation of Common RNA Force Fields. RNA 2015, 21, 1578-1590. 21. Aytenfisu, A. H.; Spasic, A.; Grossfield, A.; Stern, H. A.; Mathews, D. H. Revised RNA Dihedral Parameters for the Amber Force Field Improve RNA Molecular Dynamics. J. Chem. Theory Comput. 2017, 13, 900-915. `


Page 27 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



22. Bergonzo, C.; Cheatham, T. E. Improved Force Field Parameters Lead to a Better Description of RNA Structure. J. Chem. Theory Comput. 2015, 11, 3969-3972. 23. Tubbs, J. D.; Condon, D. E.; Kennedy, S. D.; Hauser, M.; Bevilacqua, P. C.; Turner, D. H. The Nuclear Magnetic Resonance of CCCC RNA Reveals a Right-Handed Helix, and Revised Parameters for Amber Force Field Torsions Improve Structural Predictions from Molecular Dynamics. Biochemistry 2013, 52, 996-1010. 24. Condon, D. E.; Yildirim, I.; Kennedy, S. D.; Mort, B. C.; Kierzek, R.; Turner, D. H. Optimization of an Amber Force Field for the Artificial Nucleic Acid, LNA, and Benchmarking with NMR of L(CAAU). J. Phys. Chem. B 2014, 118, 1216-1228. 25. Bottaro, S.; Banas, P.; Sponer, J.; Bussi, G. Free Energy Landscape of Gaga and UUCG RNA Tetraloops. Journal of Physical Chemistry Letters 2016, 7, 4032-4038. 26. Haldar, S.; Kuhrova, P.; Banas, P.; Spiwok, V.; Sponer, J.; Hobza, P.; Otyepka, M. Insights into Stability and Folding of GNRA and Uncg Tetra Loops Revealed by Microsecond Molecular Dynamics and Well-Tempered Metadynamics. J. Chem. Theory Comput. 2015, 11, 3866-3877. 27. Yildirim, I.; Stern, H. A.; Sponer, J.; Spackova, N.; Turner, D. H. Effects of Restrained Sampling Space and Non-Planar Amino Groups on Free Energy Predictions for RNA with Imino and Sheared Tandem GA Base Pairs Flanked by GC, CG, iGiC or iCiG Base Pairs. J. Chem. Theory Comput. 2009, 5, 2088-2100. 28. Francl, M. M.; Pietro, W. J.; Hehre, W. J.; Binkley, J. S.; Gordon, M. S.; Defrees, D. J.; Pople, J. A. Self-Consistent Molecular-Orbital Methods .23. A Polarization-Type Basis Set for 2nd-Row Elements. J. Chem. Phys. 1982, 77, 3654-3665. 29. Fock, V. Approximation Method for the Solution of the Quantum Mechanical Multibody Problems. Z. Phys. 1930, 61, 126-148. 30. Moller, C.; Plesset, M. S. Note on an Approximation Treatment for Many-Electron Systems. Phys. Rev. 1934, 46, 0618-0622. `



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 32


31. Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J., J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C., et al. Gaussian 03, Revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. 32. Hopkins, C. W.; Roitberg, A. E. Fitting of Dihedral Terms in Classical Force Fields as an Analytic Linear Least-Squares Problem. J. Chem. Inf. Model. 2014, 54, 1978-1986. 33. Vanommeslaeghe, K.; Yang, M. J.; MacKerell, A. D. Robustness in the Fitting of Molecular Mechanics Parameters. J. Comput. Chem. 2015, 36, 1083-1101. 34. Conte, M. R.; Conn, G. L.; Brown, T.; Lane, A. N. Conformational Properties and Thermodynamics of the RNA Duplex R(CGCAAAUUUGCG)(2): Comparison with the DNA Analogue D(CGCAAATTTGCG)(2). Nucleic Acids Res. 1997, 25, 2627-2634. 35. Schindelin, H.; Zhang, M.; Bald, R.; Furste, J. P.; Erdmann, V. A.; Heinemann, U. Crystal-Structure of an RNA Dodecamer Containing the Escherichia-Coli Shine-Dalgarno Sequence. J. Mol. Biol. 1995, 249, 595-603. 36. Drew, H. R.; Wing, R. M.; Takano, T.; Broka, C.; Tanaka, S.; Itakura, K.; Dickerson, R. E. Structure of a B-DNA Dodecamer - Conformation and Dynamics .1. Proc. Natl. Acad. Sci. U. S. A. 1981, 78, 2179-2183. 37. Nozinovic, S.; Furtig, B.; Jonker, H. R. A.; Richter, C.; Schwalbe, H. High-Resolution NMR Structure of an RNA Model System: The 14-Mer cUUCGg Tetraloop Hairpin RNA. Nucleic Acids Res. 2010, 38, 683-694. 38. Joung, I. S.; Cheatham, T. E. Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations. J. Phys. Chem. B 2008, 112, 9020-9041. 39. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926-935.

`


Page 29 of 32

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



40. Yildirim, I.; Park, H.; Disney, M. D.; Schatz, G. C. A Dynamic Structural Model of Expanded RNA CAG Repeats: A Refined X-Ray Structure and Computational Investigations Using Molecular Dynamics and Umbrella Sampling Simulations. J. Am. Chem. Soc. 2013, 135, 3528-3538. 41. Banas, P.; Hollas, D.; Zgarbova, M.; Jurecka, P.; Orozco, M.; Cheatham, T. E.; Sponer, J.; Otyepka, M. Performance of Molecular Mechanics Force Fields for RNA Simulations: Stability of UUCG and GNRA Hairpins. J. Chem. Theory Comput. 2010, 6, 3836-3849. 42. Condon, D. E.; Kennedy, S. D.; Mort, B. C.; Kierzek, R.; Yildirim, I.; Turner, D. H. Stacking in RNA: NMR of Four Tetramers Benchmark Molecular Dynamics. J. Chem. Theory Comput. 2015, 11, 2729-2742. 43. Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179-5197. 44. Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical-Integration of Cartesian Equations of Motion of a System with Constraints: Molecular-Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327-341. 45. Wales, D. J. Discrete Path Sampling. Mol. Phys. 2002, 100, 3285-3305. 46. Wales, D. J. Some Further Applications of Discrete Path Sampling to Cluster Isomerization. Mol. Phys. 2004, 102, 891-908. 47. Yildirim, I.; Chakraborty, D.; Disney, M. D.; Wales, D. J.; Schatz, G. C. Computational Investigation of RNA CUG Repeats Responsible for Myotonic Dystrophy 1. J. Chem. Theory Comput. 2015, 11, 4943-4958. 48. Liu, D. C.; Nocedal, J. On the Limited Memory Bfgs Method for Large-Scale Optimization. Math. Program. 1989, 45, 503-528. 49. Wales, D. J.; Miller, M. A.; Walsh, T. R. Archetypal Energy Landscapes. Nature 1998, 394, 758760. `



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 32


50. Becker, O. M.; Karplus, M. The Topology of Multidimensional Potential Energy Surfaces: Theory and Application to Peptide Structure and Kinetics. J. Chem. Phys. 1997, 106, 1495-1517. 51. Krivov, S. V.; Karplus, M. Free Energy Disconnectivity Graphs: Application to Peptide Models. J. Chem. Phys. 2002, 117, 10894-10903. 52. Strodel, B.; Whittleston, C. S.; Wales, D. J. Thermodynamics and Kinetics of Aggregation for the GNNQQNY Peptide. J. Am. Chem. Soc. 2007, 129, 16005-16014. 53. Strodel, B.; Wales, D. J. Free Energy Surfaces from an Extended Harmonic Superposition Approach and Kinetics for Alanine Dipeptide. Chem. Phys. Lett. 2008, 466, 105-115. 54. Lu, X. J.; Olson, W. K. 3DNA: A Versatile, Integrated Software System for the Analysis, Rebuilding and Visualization of Three-Dimensional Nucleic-Acid Structures. Nat. Protoc. 2008, 3, 1213-1227. 55. Lu, X. J.; Bussemaker, H. J.; Olson, W. K. DSSR: An Integrated Software Tool for Dissecting the Spatial Structure of RNA. Nucleic Acids Res. 2015, 43, e142. 56. Evans, D. A.; Wales, D. J. Free Energy Landscapes of Model Peptides and Proteins. J. Chem. Phys. 2003, 118, 3891-3897. 57. Miller,

M.;

Wales,

D.

J.;

de

Souza,

V.

Disconnectiondps,

http://www-

wales.ch.cam.ac.uk/software.html (accessed March 31, 2015). 58. Perez, A.; Marchan, I.; Svozil, D.; Sponer, J.; Cheatham, T. E.; Laughton, C. A.; Orozco, M. Refinement of the Amber Force Field for Nucleic Acids: Improving the Description of Alpha/Gamma Conformers. Biophys. J. 2007, 92, 3817-3829. 59. Chen, H.; Meisburger, S. P.; Pabit, S. A.; Sutton, J. L.; Webb, W. W.; Pollack, L. Ionic StrengthDependent Persistence Lengths of Single-Stranded RNA and DNA. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 799-804. 60. Chen, A. A.; Garcia, A. E. High-Resolution Reversible Folding of Hyperstable RNA Tetraloops Using Molecular Dynamics Simulations. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 16820-16825. `


Page 31 of 32



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

TOC Graphics

`



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

414x580mm (72 x 72 DPI)


Page 32 of 32

Improving Computational Predictions of Single-Stranded RNA

Recommend Documents