Improving Computational Predictions of Single-Stranded RNA Tetramers with Revised α/γ Torsional Parameters for the Amber Force Field

Excited-Ground-State Transition of the RNA Strand Slippage Mechanism Captured by the Base-Specific Force Field

Excited-ground-state transition and strand slippage of RNA play key roles in transcription and translation of central dogma. Due to limitation of current experimental techniques, the dynamic structure ensembles of RNA remain inadequately understood. Molecular dynamics simulations offer a promising complementary approach, whose accuracy depends on the force field. Here, we develop the new version of RNA base-specific force field (BSFF2) to address underestimation of base pairing stability and artificial backbone conformations. Extensive evaluations on typical RNA systems have comprehensively confirmed the accuracy of BSFF2. Furthermore, BSFF2 demonstrates exceptional efficiency in de novo folding of tetraloops and reproducing base pair reshuffling transition between RNA excited and ground states. Then, we explored the RNA strand slippage mechanism with BSFF2. We conducted a comprehensive three-dimensional structural investigation into the strand slippage of the most complex r(G4C2)9 repeat element and presented the molecular details in the dynamic transition along with the underlying mechanism. Our results of capturing the strand slippage, excited-ground transition, de novo folding, and simulations for various typical RNA motifs indicate that BSFF2 should be one of valuable tools for dynamic conformation research and structure prediction of RNA, and a future contribution to RNA-targeted drug design as well as RNA therapy development.

Prediction of 3D RNA Structures from Sequence Using Energy Landscapes of RNA Dimers: Application to RNA Tetraloops

Access to the three-dimensional structure of RNA enables an ability to gain a more profound understanding of its biological mechanisms, as well as the ability to design RNA-targeting drugs, which can take advantage of the unique chemical environment imposed by a folded RNA structure. Due to the dynamic and structurally complex properties of RNA, both experimental and traditional computational methods have difficulty in determining RNA's 3D structure. Herein, we introduce TAPERSS (Theoretical Analyses, Prediction, and Evaluation of RNA Structures from Sequence), a physics-based fragment assembly method for predicting 3D RNA structures from sequence. Using a fragment library created using discrete path sampling calculations of RNA dinucleoside monophosphates, TAPERSS can sample the physics-based energy landscapes of any RNA sequence with relatively low computational complexity. We have benchmarked TAPERSS on 21 RNA tetraloops, using a combinatorial algorithm as a proof-of-concept. We show that TAPERSS was successfully able to predict the apo-state structures of all 21 RNA hairpins, with 16 of those structures also having low predicted energies as well. We demonstrate that TAPERSS performs most accurately on GNRA-like tetraloops with mostly stacked loop-nucleotides, while having limited success with more dynamic UNCG and CUYG tetraloops, most likely due to the influence of the RNA force field used to create the fragment library. Moreover, we show that TAPERSS can successfully predict the majority of the experimental non-apo states, highlighting its potential in anticipating biologically significant yet unobserved states. This holds great promise for future applications in drug design and related studies. With discussed improvements and implementation of more efficient sampling algorithms, we believe TAPERSS may serve as a useful tool for a physics-based conformational sampling of large RNA structures.

Conformational dynamics of RNA G4C2 and G2C4 repeat expansions causing ALS/FTD using NMR and molecular dynamics studies

G4C2 and G2C4 repeat expansions in chromosome 9 open reading frame 72 (C9orf72) are the most common cause of genetically defined amyotrophic lateral sclerosis (ALS) and frontotemporal dementia (FTD), or c9ALS/FTD. The gene is bidirectionally transcribed, producing G4C2 repeats [r(G4C2)exp] and G2C4 repeats [r(G2C4)exp]. The c9ALS/FTD repeat expansions are highly structured, and structural studies showed that r(G4C2)exp predominantly folds into a hairpin with a periodic array of 1 × 1 G/G internal loops and a G-quadruplex. A small molecule probe revealed that r(G4C2)exp also adopts a hairpin structure with 2 × 2 GG/GG internal loops. We studied the conformational dynamics adopted by 2 × 2 GG/GG loops using temperature replica exchange molecular dynamics (T-REMD) and further characterized the structure and underlying dynamics using traditional 2D NMR techniques. These studies showed that the loop's closing base pairs influence both structure and dynamics, particularly the configuration adopted around the glycosidic bond. Interestingly, r(G2C4) repeats, which fold into an array of 2 × 2 CC/CC internal loops, are not as dynamic. Collectively, these studies emphasize the unique sensitivity of r(G4C2)exp to small changes in stacking interactions, which is not observed in r(G2C4)exp, providing important considerations for further principles in structure-based drug design.

Dynamic docking of small molecules targeting RNA CUG repeats causing myotonic dystrophy type 1

Expansion of RNA CUG repeats causes myotonic dystrophy type 1 (DM1). Once transcribed, the expanded CUG repeats strongly attract muscleblind-like 1 (MBNL1) proteins and disturb their functions in cells. Because of its unique structural form, expanded RNA CUG repeats are prospective drug targets, where small molecules can be utilized to target RNA CUG repeats to inhibit MBNL1 binding and ameliorate DM1-associated defects. In this contribution, we developed two physics-based dynamic docking approaches (DynaD and DynaD/Auto) and applied them to nine small molecules known to specifically target RNA CUG repeats. While DynaD uses a distance-based reaction coordinate to study the binding phenomenon, DynaD/Auto combines results of umbrella sampling calculations performed on 1 × 1 UU internal loops and AutoDock calculations to efficiently sample the energy landscape of binding. Predictions are compared with experimental data, displaying a positive correlation with correlation coefficient (R) values of 0.70 and 0.81 for DynaD and DynaD/Auto, respectively. Furthermore, we found that the best correlation was achieved with MM/3D-RISM calculations, highlighting the importance of solvation in binding calculations. Moreover, we detected that DynaD/Auto performed better than DynaD because of the use of prior knowledge about the binding site arising from umbrella sampling calculations. Finally, we developed dendrograms to present how bound states are connected to each other in a binding process. Results are exciting, as DynaD and DynaD/Auto will allow researchers to utilize two novel physics-based and computer-aided drug-design methodologies to perform in silico calculations on drug-like molecules aiming to target complex RNA loops.

Evaluating Geometric Definitions of Stacking for RNA Dinucleoside Monophosphates Using Molecular Mechanics Calculations

RNA modulation via small molecules is a novel approach in pharmacotherapies, where the determination of the structural properties of RNA motifs is considered a promising way to develop drugs capable of targeting RNA structures to control diseases. However, due to the complexity and dynamic nature of RNA molecules, the determination of RNA structures using experimental approaches is not always feasible, and computational models employing force fields can provide important insight. The quality of the force field will determine how well the predictions are compared to experimental observables. Stacking in nucleic acids is one such structural property, originating mainly from London dispersion forces, which are quantum mechanical and are included in molecular mechanics force fields through nonbonded interactions. Geometric descriptions are utilized to decide if two residues are stacked and hence to calculate the stacking free energies for RNA dinucleoside monophosphates (DNMPs) through statistical mechanics for comparison with experimental thermodynamics data. Here, we benchmark four different stacking definitions using molecular dynamics (MD) trajectories for 16 RNA DNMPs produced by two different force fields (RNA-IL and ff99OL3) and show that our stacking definition better correlates with the experimental thermodynamics data. While predictions within an accuracy of 0.2 kcal/mol at 300 K were observed in RNA CC, CU, UC, AG, GA, and GG, stacked states of purine-pyrimidine and pyrimidine-purine DNMPs, respectively, were typically underpredicted and overpredicted. Additionally, population distributions of RNA UU DNMPs were poorly predicted by both force fields, implying a requirement for further force field revisions. We further discuss the differences predicted by each RNA force field. Finally, we show that discrete path sampling (DPS) calculations can provide valuable information and complement the MD simulations. We propose the use of experimental thermodynamics data for RNA DNMPs as benchmarks for testing RNA force fields.

Computational Investigation of Bending Properties of RNA AUUCU, CCUG, CAG, and CUG Repeat Expansions Associated With Neuromuscular Disorders

Expansions of RNA AUUCU, CCUG, CAG, and CUG repeats cause spinocerebellar ataxia type 10, myotonic dystrophy type 2, Huntington’s disease, and myotonic dystrophy type 1, respectively. By performing extensive molecular dynamic simulations, we investigated the bending propensities and conformational landscapes adopted by 3×3, 2×2, and 1×1 internal loops observed in RNA AUUCU, CCUG, CAG, and CUG repeat expansions using model systems having biologically relevant repeat sizes. We show that the conformational variability experienced by these loops is more complex than previous reports where a variety of unconventional hydrogen bonds are formed. At the global scale, strong bending propensity was observed in r(AUUCU)₁₀, r(CCUG)₁₅, r(CAG)_20, and r(CUG)_20, and, to a lesser extent, in r(AUUCU)₄, r(CCUG)₁₀, r(CAG)₁₀, and r(CUG)₁₀. Furthermore, RNA CAG repeats exhibit a tendency toward bent states with more than 50% of observed conformations having bending angles greater than 50°, while RNA CUG repeats display relatively linear-like conformations with extremely bent conformations accounting for less than 25% of the observed structures. Conformations experienced by RNA AUUCU repeats are a combination of strongly bent and kinked structures. The bent states in RNA CCUG repeats mostly fall into the moderately bent category with a marginal ensemble experiencing extreme bending. The general pattern observed in all the bent structures indicates the collapse of the major groove width as the mechanical trigger for bending, which is caused by alteration of base pair step parameters at multiple locations along the RNA due to local distortions at the loop sites. Overextension is also observed in all the RNA repeats that is attributed to widening of the major groove width as well as undertwisting phenomenon. This information and the rich structural repository could be applied for structure based small molecule design targeting disease-causing RNAs. The bending propensities of these constructs, at the global level, could also have implications on how expanded RNA repeats interact with proteins.

Toward Convergence in Folding Simulations of RNA Tetraloops: Comparison of Enhanced Sampling Techniques and Effects of Force Field Modifications

Atomistic molecular dynamics simulations represent an established technique for investigation of RNA structural dynamics. Despite continuous development, contemporary RNA simulations still suffer from suboptimal accuracy of empirical potentials (force fields, ffs) and sampling limitations. Development of efficient enhanced sampling techniques is important for two reasons. First, they allow us to overcome the sampling limitations, and second, they can be used to quantify ff imbalances provided they reach a sufficient convergence. Here, we study two RNA tetraloops (TLs), namely the GAGA and UUCG motifs. We perform extensive folding simulations and calculate folding free energies (ΔG_fold^°) with the aim to compare different enhanced sampling techniques and to test several modifications of the nonbonded terms extending the AMBER OL3 RNA ff. We demonstrate that replica-exchange solute tempering (REST2) simulations with 12-16 replicas do not show any sign of convergence even when extended to a timescale of 120 μs per replica. However, the combination of REST2 with well-tempered metadynamics (ST-MetaD) achieves good convergence on a timescale of 5-10 μs per replica, improving the sampling efficiency by at least 2 orders of magnitude. Effects of ff modifications on ΔG_fold^° energies were initially explored by the reweighting approach and then validated by new simulations. We tested several manually prepared variants of the gHBfix potential which improve stability of the native state of both TLs by ∼2 kcal/mol. This is sufficient to conveniently stabilize the folded GAGA TL while the UUCG TL still remains under-stabilized. Appropriate adjustment of van der Waals parameters for C-H···O5' base-phosphate interaction may further stabilize the native states of both TLs by ∼0.6 kcal/mol.

Nuclear Magnetic Resonance Spectra and AMBER OL3 and ROC-RNA Simulations of UCUCGU Reveal Force Field Strengths and Weaknesses for Single-Stranded RNA

Single-stranded regions of RNA are important for folding of sequences into 3D structures and for design of therapeutics targeting RNA. Prediction of ensembles of 3D structures for single-stranded regions often involves classical mechanical approximations of interactions defined by quantum mechanical calculations on small model systems. Nuclear magnetic resonance (NMR) spectra and molecular dynamics (MD) simulations of short single strands provide tests for how well the approximations model many of the interactions. Here, the NMR spectra for UCUCGU at 2, 15, and 30 °C are compared to simulations with the AMBER force fields, OL3 and ROC-RNA. This is the first such comparison to an oligoribonucleotide containing an internal guanosine nucleotide (G). G is particularly interesting because of its many H-bonding groups, large dipole moment, and proclivity for both syn and anti conformations. Results reveal formation of a G amino to phosphate non-bridging oxygen H-bond. The results also demonstrate dramatic differences in details of the predicted structures. The variations emphasize the dependence of predictions on individual parameters and their balance with the rest of the force field. The NMR data can serve as a benchmark for future force fields.

Z-DNA as a Touchstone for Additive Empirical Force Fields and a Refinement of the Alpha/Gamma DNA Torsions for AMBER

Although current AMBER force fields are relatively accurate for canonical B-DNA, many noncanonical structures are still described incorrectly. As noncanonical motifs are attracting increasing attention due to the role they play in living organisms, further improvement is desirable. Here, we have chosen the Z-DNA molecule, which can be considered a touchstone of the universality of empirical force fields, since the noncanonical α and γ backbone conformations native to Z-DNA are also found in protein-DNA complexes, i-motif DNA, and other noncanonical DNAs. We show that spurious α/γ conformations occurring in simulations with current AMBER force fields, OL15 and bsc1, are largely due to inaccurate α/γ parametrization. Moreover, stabilization of native Z-DNA substates involving γ = trans conformations appears to be in conflict with the correct description of the canonical B-DNA structure. Because the balance of the native and spurious conformations is influenced by nonadditive effects, this is a difficult case for an additive dihedral energy scheme such as AMBER. We propose new α/γ parameters, denoted OL21, and show that they improve the stability of native α/γ Z-DNA substates while keeping the canonical DNA description virtually unchanged, thus representing a reasonable compromise within the additive force field framework. Although further extensive testing is needed, the new modification appears to be a promising step toward a more reliable description of noncanonical DNA motifs and provides the best performance for Z-DNA molecules among current AMBER force fields.

Role and Perspective of Molecular Simulation-Based Investigation of RNA-Ligand Interaction: From Small Molecules and Peptides to Photoswitchable RNA Binding

Aberrant RNA–protein complexes are formed in a variety of diseases. Identifying the ligands that interfere with their formation is a valuable therapeutic strategy. Molecular simulation, validated against experimental data, has recently emerged as a powerful tool to predict both the pose and energetics of such ligands. Thus, the use of molecular simulation may provide insight into aberrant molecular interactions in diseases and, from a drug design perspective, may allow for the employment of less wet lab resources than traditional in vitro compound screening approaches. With regard to basic research questions, molecular simulation can support the understanding of the exact molecular interaction and binding mode. Here, we focus on examples targeting RNA–protein complexes in neurodegenerative diseases and viral infections. These examples illustrate that the strategy is rather general and could be applied to different pharmacologically relevant approaches. We close this study by outlining one of these approaches, namely the light-controllable association of small molecules with RNA, as an emerging approach in RNA-targeting therapy.

A general RNA force field: comprehensive analysis of energy minima of molecular fragments of RNA

Force fields are actively used to study RNA. Development of accurate force fields relies on a knowledge of how the variation of properties of molecules depends on their structure. Detailed scrutiny of RNA's conformational preferences is needed to guide such development. Towards this end, minimum energy structures for each of a set of 16 small RNA-derived molecules were obtained by geometry optimization at the HF/6-31G(d,p), B3LYP/apc-1, and MP2/cc-pVDZ levels of theory. The number of minima computed for a given fragment was found to be related to both its size and flexibility. Atomic electrostatic multipole moments of atoms occurring in the [HO-P(O₃)-CH₂-] fragment of 30 sugar-phosphate-sugar geometries were calculated at the HF/6-31G(d,p) and B3LYP/apc-1 levels of theory, and the transferability of these properties between different conformations was investigated. The atomic multipole moments were found to be highly transferable between different conformations with small standard deviations. These results indicate necessary elements of the development of accurate RNA force fields.

Molecular Dynamics Simulation of the Conformational Preferences of Pseudouridine Derivatives: Improving the Distribution in the Glycosidic Torsion Space

There are only four derivatives of pseudouridine (Ψ) that are known to occur naturally in RNA as post-transcriptional modifications. We have studied the conformational consequences of pseudouridylation and further modifications using replica exchange molecular dynamics simulations at the nucleoside level, and the simulated conformational preferences were compared with the available experimental (NMR) data. We found that the existing AMBER FF99-derived parameters for these nucleosides did not reproduce the observed experimental features and while the recommended bsc0 correction could be combined with these parameters leading to an improvement in the description of sugar pucker distributions, the χ_OL3 correction could not be applied to these nucleosides as such because of base isomerization. On the other hand, the revised χ torsion parameters (χ_IDRP) for Ψ developed earlier by us (Deb, I., J. Comput. Chem., 2016, 37, 1576-1588) in combination with the AMBER provided parameters and the revised γ torsion parameters generated conformational distributions, which generally were in better agreement with the experimental data. A significant shift of the distribution of base orientation toward the syn conformation was observed with our revised parameter sets compared to the large excess of anti conformation predicted by the FF99 parameters. Overall, our observations indicated that our revised set of parameters (χ_IDRP) for Ψ were also able to generate conformational distributions for all of the derivatives of Ψ in better agreement with the experimental data.

Free-Energy Calculation of Ribonucleic Inosines and Its Application to Nearest-Neighbor Parameters

Can current simulations quantitatively predict the stability of ribonucleic acids (RNAs)? In this research, we apply a free-energy perturbation simulation of RNAs containing inosine, a modified ribonucleic base, to the derivation of RNA nearest-neighbor parameters. A parameter set derived solely from 30 simulations was used to predict the free-energy difference of the RNA duplex with a mean unbiased error of 0.70 kcal/mol, which is a level of accuracy comparable to that obtained with parameters derived from 25 experiments. We further show that the error can be lowered to 0.60 kcal/mol by combining the simulation-derived free-energy differences with experimentally measured differences. This protocol can be used as a versatile method for deriving nearest-neighbor parameters of RNAs with various modified bases.

Development and Testing of the OPLS-AA/M Force Field for RNA

Significant improvements have been made to the OPLS-AA force field for modeling RNA. New torsional potentials were optimized based on density functional theory (DFT) scans at the ωB97X-D/6-311++G(d,p) level for potential energy surfaces of the backbone α and γ dihedral angles. In combination with previously reported improvements for the sugar puckering and glycosidic torsion terms, the new force field was validated through diverse molecular dynamics simulations for RNAs in aqueous solution. Results for dinucleotides and tetranucleotides revealed both accurate reproduction of 3 J couplings from NMR and the avoidance of several unphysical states observed with other force fields. Simulations of larger systems with noncanonical motifs showed significant structural improvements over the previous OPLS-AA parameters. The new force field, OPLS-AA/M, is expected to perform competitively with other recent RNA force fields and to be compatible with OPLS-AA models for proteins and small molecules.

Computational Investigation of RNA A-Bulges Related to the Microtubule-Associated Protein Tau Causing Frontotemporal Dementia and Parkinsonism

Mutations in the human tau gene result in alternative splicing of the tau protein, which causes frontotemporal dementia and Parkinsonism. One disease mechanism is linked to the stability of a hairpin within the microtubule-associated protein tau (MAPT) mRNA, which contains an A-bulge. Here we employ computational methods to investigate the structural and thermodynamic properties of several A-bulge RNAs with different closing base-pairs. We find that the current amber RNA force field has a preference to overstabilize base-triple over stacked states, even though some of the A-bulges are known to prefer stacked states according to NMR studies. We further determined that if the neighboring base-pairs of A-bulges are AU, this situation can lead to base slippage. However, when the 3'-side of the A-bulge has an UA base-pair, the stacked state is stabilized by an extra interaction that is not observed in the other sequences. We suggest that these A-bulge RNA systems could be used as benchmarks to improve the current RNA force fields.

Improvement of RNA Simulations with Torsional Revisions of the AMBER Force Field

Our current knowledge on the unique roles of RNA in cells makes it vital to investigate the properties of RNA systems using computational methods because of the potential pharmaceutical applications. With the continuous advancement of computer technology, it is now possible to study RNA folding. Molecular mechanics calculations are useful in discovering the structural and thermodynamic properties of RNA systems. Yet, the predictions depend on the quality of the RNA force field, which is a set of parameters describing the potential energy of the system. Torsional parameters are one of the terms in a force field that can be revised using physics-based approaches. This chapter focuses on improvements provided by revisions of torsional parameters of the AMBER (Assisted Model Building with Energy Refinement) RNA force field. The theory behind torsional revisions and re-parameterization of several RNA torsions is briefly described. Applications of the revised torsional parameters to study RNA nucleosides, single-stranded RNA tetramers, and RNA repeat expansions are described in detail. It is concluded that RNA force fields require constant revisions and should be benchmarked against diverse RNA systems such as single strands and internal loops in order to test their qualities.

The Hairpin Form of r(G4C2)exp in c9ALS/FTD Is Repeat-Associated Non-ATG Translated and a Target for Bioactive Small Molecules

The most common genetic cause of amyotrophic lateral sclerosis (ALS) and frontotemporal dementia (FTD) is an expanded G₄C₂ repeat [(G₄C₂)^exp] in C9ORF72. ALS/FTD-associated toxicity has been traced to the RNA transcribed from the repeat expansion [r(G₄C₂)^exp], which sequesters RNA-binding proteins (RBPs) and undergoes repeat-associated non-ATG (RAN) translation to generate toxic dipeptide repeats. Using in vitro and cell-based assays, we identified a small molecule (4) that selectively bound r(G₄C₂)^exp, prevented sequestration of an RBP, and inhibited RAN translation. Indeed, biophysical characterization showed that 4 selectively bound the hairpin form of r(G₄C₂)^exp, and nuclear magnetic resonance spectroscopy studies and molecular dynamics simulations defined this molecular recognition event. Cellular imaging revealed that 4 localized to r(G₄C₂)^exp cytoplasmic foci, the putative sites of RAN translation. Collectively, these studies highlight that the hairpin structure of r(G₄C₂)^exp is a therapeutically relevant target and small molecules that bind it can ameliorate c9ALS/FTD-associated toxicity.

RNA Structural Dynamics As Captured by Molecular Simulations: A Comprehensive Overview

With both catalytic and genetic functions, ribonucleic acid (RNA) is perhaps the most pluripotent chemical species in molecular biology, and its functions are intimately linked to its structure and dynamics. Computer simulations, and in particular atomistic molecular dynamics (MD), allow structural dynamics of biomolecular systems to be investigated with unprecedented temporal and spatial resolution. We here provide a comprehensive overview of the fast-developing field of MD simulations of RNA molecules. We begin with an in-depth, evaluatory coverage of the most fundamental methodological challenges that set the basis for the future development of the field, in particular, the current developments and inherent physical limitations of the atomistic force fields and the recent advances in a broad spectrum of enhanced sampling methods. We also survey the closely related field of coarse-grained modeling of RNA systems. After dealing with the methodological aspects, we provide an exhaustive overview of the available RNA simulation literature, ranging from studies of the smallest RNA oligonucleotides to investigations of the entire ribosome. Our review encompasses tetranucleotides, tetraloops, a number of small RNA motifs, A-helix RNA, kissing-loop complexes, the TAR RNA element, the decoding center and other important regions of the ribosome, as well as assorted others systems. Extended sections are devoted to RNA-ion interactions, ribozymes, riboswitches, and protein/RNA complexes. Our overview is written for as broad of an audience as possible, aiming to provide a much-needed interdisciplinary bridge between computation and experiment, together with a perspective on the future of the field.

Mapping the Universe of RNA Tetraloop Folds

We report a map of RNA tetraloop conformations constructed by calculating pairwise distances among all experimentally determined four-nucleotide hairpin loops. Tetraloops with similar structures are clustered together and, as expected, the two largest clusters are the canonical GNRA and UNCG folds. We identify clusters corresponding to known tetraloop folds such as GGUG, RNYA, AGNN, and CUUG. These clusters are represented in a simple two-dimensional projection that recapitulates the relationship among the different folds. The cluster analysis also identifies 20 novel tetraloop folds that are peculiar to specific positions in ribosomal RNAs and that are stabilized by tertiary interactions. In our RNA tetraloop database we find a significant number of non-GNRA and non-UNCG sequences adopting the canonical GNRA and UNCG folds. Conversely, we find a significant number of GNRA and UNCG sequences adopting non-GNRA and non-UNCG folds. Our analysis demonstrates that there is not a simple one-to-one, but rather a many-to-many mapping between tetraloop sequence and tetraloop fold.