Quantum chemical studies of nucleic acids: can we construct a bridge to the RNA structural biology and bioinformatics communities?

New classifications for quantum bioinformatics: Q-bioinformatics, QCt-bioinformatics, QCg-bioinformatics, and QCr-bioinformatics

Bioinformatics has revolutionized biology and medicine by using computational methods to analyze and interpret biological data. Quantum mechanics has recently emerged as a promising tool for the analysis of biological systems, leading to the development of quantum bioinformatics. This new field employs the principles of quantum mechanics, quantum algorithms, and quantum computing to solve complex problems in molecular biology, drug design, and protein folding. However, the intersection of bioinformatics, biology, and quantum mechanics presents unique challenges. One significant challenge is the possibility of confusion among scientists between quantum bioinformatics and quantum biology, which have similar goals and concepts. Additionally, the diverse calculations in each field make it difficult to establish boundaries and identify purely quantum effects from other factors that may affect biological processes. This review provides an overview of the concepts of quantum biology and quantum mechanics and their intersection in quantum bioinformatics. We examine the challenges and unique features of this field and propose a classification of quantum bioinformatics to promote interdisciplinary collaboration and accelerate progress. By unlocking the full potential of quantum bioinformatics, this review aims to contribute to our understanding of quantum mechanics in biological systems.

Polarizable Molecular Block Model: Toward the Development of an Induced Dipole Force Field for DNA

For flexible and highly ionized macromolecules such as DNA, it is important to correctly evaluate the intramolecular polarization in an induced dipole force field. In a proposed polarizable molecular block (PMB) model, a large molecule is divided into several molecular blocks. The atomic charges of the blocks are optimized by using the respective electrostatic potentials (ESPs) on the molecular surface. By using the capped hydrogen removal operation, the total charge of the blocks is controlled exactly to have an integer charge. The atomic polarizabilities of the blocks are optimized by using the respective polarized one-electron potentials that are the differences between ESPs with and without an external test charge. Induced dipole-charge interactions between the blocks are all included, but those interactions within the blocks are strictly excluded. All dipole-dipole interactions are included, but the damping functions are applied to the close dipole-dipole pairs. Several types of damping (simple scaling, exponential, linear, and Gaussian) are evaluated. The validity of the PMB model was verified by using trinucleotide duplexes which have A-, B-, and Z-DNA forms. The reference energies of trinucleotide duplexes including counterions (GGT3Na-ACC3Na, GAC3Na-GTC3Na, and GCG3Na-CGC3Na) are calculated using ωB97XD/aug-cc-pVDZ. All damping types reproduced well the reference interaction energies, dipole moments, and ESPs. Among them, the simple scaling with strong attenuation to 1-2 atomic pairs showed the highest stability against the polarization catastrophe. This study shows that it is possible to develop a high-quality polarizable force field by treating the intramolecular polarization on a block-by-block basis.

Short-Range Imbalances in the AMBER Lennard-Jones Potential for (Deoxy)Ribose···Nucleobase Lone-Pair···π Contacts in Nucleic Acids

The lone-pair···π (lp···π) (deoxy)ribose···nucleobase stacking is a recurring interaction in Z-DNA and RNAs that is characterized by sub-van der Waals lp···π contacts (<3.0 Å). It is a part of the structural signature of CpG Z-step motifs in Z-DNA and r(UNCG) tetraloops that are known to behave poorly in molecular dynamics (MD) simulations. Although the exact origin of the MD simulation issues remains unclear, a significant part of the problem might be due to an imbalanced description of nonbonded interactions, including the characteristic lp···π stacking. To gain insights into the links between lp···π stacking and MD, we present an in-depth comparison between accurate large-basis-set double-hybrid Kohn-Sham density functional theory calculations DSD-BLYP-D3/ma-def2-QZVPP (DHDF-D3) and data obtained with the nonbonded potential of the AMBER force field (AFF) for NpN Z-steps (N = G, A, C, and U). Among other differences, we found that the AFF overestimates the DHDF-D3 lp···π distances by ∼0.1-0.2 Å, while the deviation between the DHDF-D3 and AFF descriptions sharply increases in the short-range region of the interaction. Based on atom-in-molecule polarizabilities and symmetry-adapted perturbation theory analysis, we inferred that the DHDF-D3 versus AFF differences partly originate in identical nucleobase carbon atom Lennard-Jones (LJ) parameters despite the presence/absence of connected electron-withdrawing groups that lead to different effective volumes or vdW radii. Thus, to precisely model the very short CpG lp···π contact distances, we recommend revision of the nucleobase atom LJ parameters. Additionally, we suggest that the large discrepancy between DHDF-D3 and AFF short-range repulsive part of the interaction energy potential may significantly contribute to the poor performances of MD simulations of nucleic acid systems containing Z-steps. Understanding where, and if possible why, the point-charge-type effective potentials reach their limits is vital for developing next-generation FFs and for addressing specific issues in contemporary MD simulations.

Chemical reactivity and binding interactions in ribonucleic acid-peptide complexes

The covalent and noncovalent backbone binding interactions in RNA-peptide complexes were studied by DFT methods. Four RNA structures R1(GGCUAGCC), R2(AAUCGAUU), R3(GGGAUCCC), and R4(AAAGCUUU) has been selected for eight protonated peptides (DR, ER, GR, KR, NGR, RR, tmeGnd (tme), and VR) interactions based on an experimental study (Anal Chem. 2019; 91:1659-1664). Chemical reactivity theory is used to study the reactivity of eight peptides with global descriptors. Lower hardness values reflected low stability and high reactivity for the protonated peptides. DR, ER, GR, KR, NGR, RR, and VR show lower value of ω, μ while tme has high value of ω, μ. Larger highest occupied molecular orbital (HOMO)-lowest unoccupied molecular orbital (LUMO) gap for ER, GR, and KR showed greater structural stability for peptides. AutoDock and PatchDock results indicated that R1, R2, and R4 retain hairpin structures while interacting with peptide complexes. The calculated binding energies of (R1-R4)-peptide complexes from AutoDock tools are (1.49-11.12) kcal/mol. Results showed that the noncovalent interactions are stronger than the covalent interactions for R1-peptide complexes. The reason might be the transfer of proton from protonated ligand to deprotonated RNA, which has initiated the loss of the ligand. Also it has been observed that proton transfer has become energetically unfavorable in presence of additional hydrogen bonds which is predicted in the experimental results.

Revisiting the Potential Energy Surface of the Stacked Cytosine Dimer: FNO-CCSD(T) Interaction Energies, SAPT Decompositions, and Benchmarking

Nucleobase stacking interactions are crucial for the stability of nucleic acids. This study investigates base stacking energies of the cytosine homodimer in different configurations, including intermolecular separation plots, detailed twist dependence, and displaced structures. Highly accurate ab initio quantum chemical single point energies using an energy function based on MP2 complete basis set extrapolation ([6 → 7]ZaPa-NR) and a CCSD(T)/cc-pVTZ-F12 high-level correction are presented as new reference data, providing the most accurate stacking energies of nucleobase dimers currently available. Accurate SAPT2+(3)δMP2 energy decomposition is used to obtain detailed insights into the nature of base stacking interactions at varying vertical distances and twist values. The ab initio symmetry adapted perturbation theory (SAPT) energy decomposition suggests that the base stacking originates from an intricate interplay between dispersion attraction, short-range exchange-repulsion, and Coulomb interaction. The interpretation of the SAPT data is a complex issue as key energy terms vary substantially in the region of optimal (low energy) base stacking geometries. Thus, attempts to highlight one leading stabilizing SAPT base stacking term may be misleading and the outcome strongly depends on the used geometries within the range of geometries sampled in nucleic acids upon thermal fluctuations. Modern dispersion-corrected density functional theory (among them DSD-BLYP-D3, ωB97M-V, and ωB97M-D3BJ) is benchmarked and often reaches up to spectroscopic accuracy (below 1 kJ/mol). The classical AMBER force field is benchmarked with multiple different sets of point-charges (e.g. HF, DFT, and MP2-based) and is found to produce reasonable agreement with the benchmark data.

Gene editing in the context of an increasingly complex genome

The reporting of the first draft of the human genome in 2000 brought with it much hope for the future in what was felt as a paradigm shift toward improved health outcomes. Indeed, we have now mapped the majority of variation across human populations with landmark projects such as 1000 Genomes; in cancer, we have catalogued mutations across the primary carcinomas; whilst, for other diseases, we have identified the genetic variants with strongest association. Despite this, we are still awaiting the genetic revolution in healthcare to materialise and translate itself into the health benefits for which we had hoped. A major problem we face relates to our underestimation of the complexity of the genome, and that of biological mechanisms, generally. Fixation on DNA sequence alone and a 'rigid' mode of thinking about the genome has meant that the folding and structure of the DNA molecule -and how these relate to regulation- have been underappreciated. Projects like ENCODE have additionally taught us that regulation at the level of RNA is just as important as that at the spatiotemporal level of chromatin.In this review, we chart the course of the major advances in the biomedical sciences in the era pre- and post the release of the first draft sequence of the human genome, taking a focus on technology and how its development has influenced these. We additionally focus on gene editing via CRISPR/Cas9 as a key technique, in particular its use in the context of complex biological mechanisms. Our aim is to shift the mode of thinking about the genome to that which encompasses a greater appreciation of the folding of the DNA molecule, DNA- RNA/protein interactions, and how these regulate expression and elaborate disease mechanisms.Through the composition of our work, we recognise that technological improvement is conducive to a greater understanding of biological processes and life within the cell. We believe we now have the technology at our disposal that permits a better understanding of disease mechanisms, achievable through integrative data analyses. Finally, only with greater understanding of disease mechanisms can techniques such as gene editing be faithfully conducted.

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.

Physics-based all-atom modeling of RNA energetics and structure

The database of RNA sequences is exploding, but knowledge of energetics, structures, and dynamics lags behind. All-atom computational methods, such as molecular dynamics, hold promise for closing this gap. New algorithms and faster computers have accelerated progress in improving the reliability and accuracy of predictions. Currently, the methods can facilitate refinement of experimentally determined nuclear magnetic resonance and x-ray structures, but are 'unreliable' for predictions based only on sequence. Much remains to be discovered, however, about the many molecular interactions driving RNA folding and the best way to approximate them quantitatively. The large number of parameters required means that a wide variety of experimental results will be required to benchmark force fields and different approaches. As computational methods become more reliable and accessible, they will be used by an increasing number of biologists, much as x-ray crystallography has expanded. Thus, many fundamental physical principles underlying the computational methods are described. This review presents a summary of the current state of molecular dynamics as applied to RNA. It is designed to be helpful to students, postdoctoral fellows, and faculty who are considering or starting computational studies of RNA. WIREs RNA 2017, 8:e1422. doi: 10.1002/wrna.1422.

Ensemble simulations: folding, unfolding and misfolding of a high-efficiency frameshifting RNA pseudoknot

Massive all-atom molecular dynamics simulations were conducted across a distributed computing network to study the folding, unfolding, misfolding and conformational plasticity of the high-efficiency frameshifting double mutant of the 26 nt potato leaf roll virus RNA pseudoknot. Our robust sampling, which included over 40 starting structures spanning the spectrum from the extended unfolded state to the native fold, yielded nearly 120 μs of cumulative sampling time. Conformational microstate transitions on the 1.0 ns to 10.0 μs timescales were observed, with post-equilibration sampling providing detailed representations of the conformational free energy landscape and the complex folding mechanism inherent to the pseudoknot motif. Herein, we identify and characterize two alternative native structures, three intermediate states, and numerous misfolded states, the latter of which have not previously been characterized via atomistic simulation techniques. While in line with previous thermodynamics-based models of a general RNA folding mechanism, our observations indicate that stem-strand-sequence-separation may serve as an alternative predictor of the order of stem formation during pseudoknot folding. Our results contradict a model of frameshifting based on structural rigidity and resistance to mechanical unfolding, and instead strongly support more recent studies in which conformational plasticity is identified as a determining factor in frameshifting efficiency.

Structural landscape of base pairs containing post-transcriptional modifications in RNA

Base pairs involving post-transcriptionally modified nucleobases are believed to play important roles in a wide variety of functional RNAs. Here we present our attempts toward understanding the structural and functional role of naturally occurring modified base pairs using a combination of X-ray crystal structure database analysis, sequence analysis, and advanced quantum chemical methods. Our bioinformatics analysis reveals that despite their presence in all major secondary structural elements, modified base pairs are most prevalent in tRNA crystal structures and most commonly involve guanine or uridine modifications. Further, analysis of tRNA sequences reveals additional examples of modified base pairs at structurally conserved tRNA regions and highlights the conservation patterns of these base pairs in three domains of life. Comparison of structures and binding energies of modified base pairs with their unmodified counterparts, using quantum chemical methods, allowed us to classify the base modifications in terms of the nature of their electronic structure effects on base-pairing. Analysis of specific structural contexts of modified base pairs in RNA crystal structures revealed several interesting scenarios, including those at the tRNA:rRNA interface, antibiotic-binding sites on the ribosome, and the three-way junctions within tRNA. These scenarios, when analyzed in the context of available experimental data, allowed us to correlate the occurrence and strength of modified base pairs with their specific functional roles. Overall, our study highlights the structural importance of modified base pairs in RNA and points toward the need for greater appreciation of the role of modified bases and their interactions, in the context of many biological processes involving RNA.

How Does Mg2+ Modulate the RNA Folding Mechanism: A Case Study of the G:C W:W Trans Basepair

Reverse Watson-Crick G:C basepairs (G:C W:W Trans) occur frequently in different functional RNAs. This is one of the few basepairs whose gas-phase-optimized isolated geometry is inconsistent with the corresponding experimental geometry. Several earlier studies indicate that through post-transcriptional modification, direct protonation, or coordination with Mg²⁺, accumulation of positive charge near N7 of guanine can stabilize the experimental geometry. Interestingly, recent studies reveal significant variation in the position of putatively bound Mg²⁺. This, in conjunction with recently raised doubts regarding some of the Mg²⁺ assignments near the imino nitrogen of guanine, is suggestive of the existence of multiple Mg²⁺ binding modes for this basepair. Our detailed investigation of Mg²⁺-bound G:C W:W Trans pairs occurring in high-resolution RNA crystal structures shows that they are found in 14 different contexts, eight of which display Mg²⁺ binding at the Hoogsteen edge of guanine. Further examination of occurrences in these eight contexts led to the characterization of three different Mg²⁺ binding modes: 1) direct binding via N7 coordination, 2) direct binding via O6 coordination, and 3) binding via hydrogen-bonding interaction with the first-shell water molecules. In the crystal structures, the latter two modes are associated with a buckled and propeller-twisted geometry of the basepair. Interestingly, respective optimized geometries of these different Mg²⁺ binding modes (optimized using six different DFT functionals) are consistent with their corresponding experimental geometries. Subsequent interaction energy calculations at the MP2 level, and decomposition of its components, suggest that for G:C W:W Trans , Mg²⁺ binding can fine tune the basepair geometries without compromising with their stability. Our results, therefore, underline the importance of the mode of binding of Mg²⁺ ions in shaping RNA structure, folding and function.

Polarizable Multipole-Based Force Field for Aromatic Molecules and Nucleobases

Aromatic molecules with π electrons are commonly involved in chemical and biological recognitions. For example, nucleobases play central roles in DNA/RNA structure and their interactions with proteins. The delocalization of the π electrons is responsible for the high polarizability of aromatic molecules. In this work, the AMOEBA force field has been developed and applied to 5 regular nucleobases and 12 aromatic molecules. The permanent electrostatic energy is expressed as atomic multipole interactions between atom pairs, and many-body polarization is accounted for by mutually induced atomic dipoles. We have systematically investigated aromatic ring stacking and aromatic-water interactions for nucleobases and aromatic molecules, as well as base–base hydrogen-bonding pair interactions, all at various distances and orientations. van der Waals parameters were determined by comparison to the quantum mechanical interaction energy of these dimers and fine-tuned using condensed phase simulation. By comparing to quantum mechanical calculations, we show that the resulting classical potential is able to accurately describe molecular polarizability, molecular vibrational frequency, and dimer interaction energy of these aromatic systems. Condensed phase properties, including hydration free energy, liquid density, and heat of vaporization, are also in good overall agreement with experimental values. The structures of benzene liquid phase and benzene-water solution were also investigated by simulation and compared with experimental and PDB structure derived statistical results.

Higher order structures involving post transcriptionally modified nucleobases in RNA

Quantum chemical studies are carried out to understand the structures and stabilities of higher order structures involving post-transcriptionally modified nucleobases in RNA.

Comparative Assessment of Different RNA Tetranucleotides from the DFT-D3 and Force Field Perspective

Classical force field (FF) molecular dynamics (MD) simulations of RNA tetranucleotides have substantial problems in reproducing conformer populations indicated by NMR experiments. To provide more information about the possible sources of errors, we performed quantum mechanical (QM, TPSS-D3/def2-TZVP) and molecular mechanics (MM, AMBER parm99bsc0+χ_OL3) calculations of different r(CCCC), r(GACC), and r(UUUU) conformers obtained from explicit solvent MD simulations. Solvent effects in the static QM and MM calculations were mimicked using implicit solvent models (COSMO and Poisson-Boltzmann, respectively). The comparison of QM and MM geometries and energies revealed that the two methodologies provide qualitatively consistent results in most of the cases. Even though we found some differences, these were insufficient to indicate any systematic corrections of the RNA FF terms that could improve the performance of classical MD in simulating tetranucleotides. On the basis of these findings, we inferred that the overpopulation of intercalated conformers in the MD simulations of RNA tetramers, which were not observed experimentally, might be predominantly caused by imbalanced water-solvent and water-water interactions. Apart from the large-scale QM calculations performed to assess the performance of the AMBER FF, a representative spectrum of faster QM methods was tested.

Quantum chemical benchmark study on 46 RNA backbone families using a dinucleotide unit

We have created a benchmark set of quantum chemical structure-energy data denoted as UpU46, which consists of 46 uracil dinucleotides (UpU), representing all known 46 RNA backbone conformational families. Penalty-function-based restrained optimizations with COSMO TPSS-D3/def2-TZVP ensure a balance between keeping the target conformation and geometry relaxation. The backbone geometries are close to the clustering-means of their respective RNA bioinformatics family classification. High-level wave function methods (DLPNO-CCSD(T) as reference) and a wide-range of dispersion-corrected or inclusive DFT methods (DFT-D3, VV10, LC-BOP-LRD, M06-2X, M11, and more) are used to evaluate the conformational energies. The results are compared to the Amber RNA bsc0χOL3 force field. Most dispersion-corrected DFT methods surpass the Amber force field significantly in accuracy and yield mean absolute deviations (MADs) for relative conformational energies of ∼0.4-0.6 kcal/mol. Double-hybrid density functionals represent the most accurate class of density functionals. Low-cost quantum chemical methods such as PM6-D3H+, HF-3c, DFTB3-D3, as well as small basis set calculations corrected for basis set superposition errors (BSSEs) by the gCP procedure are also tested. Unfortunately, the presently available low-cost methods are struggling to describe the UpU conformational energies with satisfactory accuracy. The UpU46 benchmark is an ideal test for benchmarking and development of fast methods to describe nucleic acids, including force fields.

An introduction to recurrent nucleotide interactions in RNA

RNA secondary structure diagrams familiar to molecular biologists summarize at a glance the folding of RNA chains to form Watson–Crick paired double helices. However, they can be misleading: First of all, they imply that the nucleotides in loops and linker segments, which can amount to 35% to 50% of a structured RNA, do not significantly interact with other nucleotides. Secondly, they give the impression that RNA molecules are loosely organized in three-dimensional (3D) space. In fact, structured RNAs are compactly folded as a result of numerous long-range, sequence-specific interactions, many of which involve loop or linker nucleotides. Here, we provide an introduction for students and researchers of RNA on the types, prevalence, and sequence variations of inter-nucleotide interactions that structure and stabilize RNA 3D motifs and architectures, using Escherichia coli (E. coli) 16S ribosomal RNA as a concrete example. The picture that emerges is that almost all nucleotides in structured RNA molecules, including those in nominally single-stranded loop or linker regions, form specific interactions that stabilize functional structures or mediate interactions with other molecules. The small number of noninteracting, ‘looped-out’ nucleotides make it possible for the RNA chain to form sharp turns. Base-pairing is the most specific interaction in RNA as it involves edge-to-edge hydrogen bonding (H-bonding) of the bases. Non-Watson–Crick base pairs are a significant fraction (30% or more) of base pairs in structured RNAs.

Identifying novel sequence variants of RNA 3D motifs

Predicting RNA 3D structure from sequence is a major challenge in biophysics. An important sub-goal is accurately identifying recurrent 3D motifs from RNA internal and hairpin loop sequences extracted from secondary structure (2D) diagrams. We have developed and validated new probabilistic models for 3D motif sequences based on hybrid Stochastic Context-Free Grammars and Markov Random Fields (SCFG/MRF). The SCFG/MRF models are constructed using atomic-resolution RNA 3D structures. To parameterize each model, we use all instances of each motif found in the RNA 3D Motif Atlas and annotations of pairwise nucleotide interactions generated by the FR3D software. Isostericity relations between non-Watson–Crick basepairs are used in scoring sequence variants. SCFG techniques model nested pairs and insertions, while MRF ideas handle crossing interactions and base triples. We use test sets of randomly-generated sequences to set acceptance and rejection thresholds for each motif group and thus control the false positive rate. Validation was carried out by comparing results for four motif groups to RMDetect. The software developed for sequence scoring (JAR3D) is structured to automatically incorporate new motifs as they accumulate in the RNA 3D Motif Atlas when new structures are solved and is available free for download.

RNABP COGEST: a resource for investigating functional RNAs

Structural bioinformatics of RNA has evolved mainly in response to the rapidly accumulating evidence that non-(protein)-coding RNAs (ncRNAs) play critical roles in gene regulation and development. The structures and functions of most ncRNAs are however still unknown. Most of the available RNA structural databases rely heavily on known 3D structures, and contextually correlate base pairing geometry with actual 3D RNA structures. None of the databases provide any direct information about stabilization energies. However, the intrinsic interaction energies of constituent base pairs can provide significant insights into their roles in the overall dynamics of RNA motifs and structures. Quantum mechanical (QM) computations provide the only approach toward their accurate quantification and characterization. ‘RNA Base Pair Count, Geometry and Stability’ (http://bioinf.iiit.ac.in/RNABPCOGEST) brings together information, extracted from literature data, regarding occurrence frequency, experimental and quantum chemically optimized geometries, and computed interaction energies, for non-canonical base pairs observed in a non-redundant dataset of functional RNA structures. The database is designed to enable the QM community, on the one hand, to identify appropriate biologically relevant model systems and also enable the biology community to easily sift through diverse computational results to gain theoretical insights which could promote hypothesis driven biological research.

Database URL:http://bioinf.iiit.ac.in/RNABPCOGEST

The role of N7 protonation of guanine in determining the structure, stability and function of RNA base pairs

Ab initio computations and bioinformatics studies reveal that stabilization of some important RNA structural motifs might involve N7 protonation of guanine.

Towards biochemically relevant QM computations on nucleic acids: controlled electronic structure geometry optimization of nucleic acid structural motifs using penalty restraint functions

Recent developments in dispersion-corrected density functional theory methods allow for the first time the description of large fragments of nucleic acids (hundreds of atoms) with an accuracy clearly surpassing the accuracy of common biomolecular force fields. Such calculations can significantly improve the description of the potential energy surface of nucleic acid molecules, which may be useful for studies of molecular interactions and conformational preferences of nucleic acids, as well as verification and parameterization of other methods. The first of such studies, however, demonstrated that successful applications of accurate QM calculations to larger nucleic acid building blocks are hampered by difficulties in obtaining geometries that are biochemically relevant and are not biased by non-native structural features. We present an approach that can greatly facilitate large-scale QM studies on nucleic acids, namely electronic structure geometry optimization of nucleic acid fragments utilizing a penalty function to restrain key internal coordinates with a specific focus on the torsional backbone angles. This work explores the viability of these restraint optimizations for DFT-D3, PM6-D3H and HF-3c optimizations on a set of examples (a UpA dinucleotide, a DNA G-quadruplex and a B-DNA fragment). Evaluation of different penalty function strengths reveals only a minor system-dependency and reasonable restraint values range from 0.01 to 0.05 Eh rad(-2) for the backbone torsions. Restraints are crucial to perform the QM calculations on biochemically relevant conformations in implicit solvation and gas phase geometry optimizations. The reasons for using restrained instead of constrained or unconstrained optimizations are explained and an open-source external optimizer is provided.

Energies and 2'-Hydroxyl Group Orientations of RNA Backbone Conformations. Benchmark CCSD(T)/CBS Database, Electronic Analysis, and Assessment of DFT Methods and MD Simulations

Sugar-phosphate backbone is an electronically complex molecular segment imparting RNA molecules high flexibility and architectonic heterogeneity necessary for their biological functions. The structural variability of RNA molecules is amplified by the presence of the 2'-hydroxyl group, capable of forming multitude of intra- and intermolecular interactions. Bioinformatics studies based on X-ray structure database revealed that RNA backbone samples at least 46 substates known as rotameric families. The present study provides a comprehensive analysis of RNA backbone conformational preferences and 2'-hydroxyl group orientations. First, we create a benchmark database of estimated CCSD(T)/CBS relative energies of all rotameric families and test performance of dispersion-corrected DFT-D3 methods and molecular mechanics in vacuum and in continuum solvent. The performance of the DFT-D3 methods is in general quite satisfactory. The B-LYP-D3 method provides the best trade-off between accuracy and computational demands. B3-LYP-D3 slightly outperforms the new PW6B95-D3 and MPW1B95-D3 and is the second most accurate density functional of the study. The best agreement with CCSD(T)/CBS is provided by DSD-B-LYP-D3 double-hybrid functional, although its large-scale applications may be limited by high computational costs. Molecular mechanics does not reproduce the fine energy differences between the RNA backbone substates. We also demonstrate that the differences in the magnitude of the hyperconjugation effect do not correlate with the energy ranking of the backbone conformations. Further, we investigated the 2'-hydroxyl group orientation preferences. For all families, we conducted a QM and MM hydroxyl group rigid scan in gas phase and solvent. We then carried out set of explicit solvent MD simulations of folded RNAs and analyze 2'-hydroxyl group orientations of different backbone families in MD. The solvent energy profiles determined primarily by the sugar pucker match well with the distribution data derived from the simulations. The QM and MM energy profiles predict the same 2'-hydroxyl group orientation preferences. Finally, we demonstrate that the high energy of unfavorable and rarely sampled 2'-hydroxyl group orientations can be attributed to clashes between occupied orbitals.

Isosteric and nonisosteric base pairs in RNA motifs: molecular dynamics and bioinformatics study of the sarcin-ricin internal loop

The sarcin-ricin RNA motif (SR motif) is one of the most prominent recurrent RNA building blocks that occurs in many different RNA contexts and folds autonomously, that is, in a context-independent manner. In this study, we combined bioinformatics analysis with explicit-solvent molecular dynamics (MD) simulations to better understand the relation between the RNA sequence and the evolutionary patterns of the SR motif. A SHAPE probing experiment was also performed to confirm the fidelity of the MD simulations. We identified 57 instances of the SR motif in a nonredundant subset of the RNA X-ray structure database and analyzed their base pairing, base-phosphate, and backbone-backbone interactions. We extracted sequences aligned to these instances from large rRNA alignments to determine the frequency of occurrence for different sequence variants. We then used a simple scoring scheme based on isostericity to suggest 10 sequence variants with a highly variable expected degree of compatibility with the SR motif 3D structure. We carried out MD simulations of SR motifs with these base substitutions. Nonisosteric base substitutions led to unstable structures, but so did isosteric substitutions which were unable to make key base-phosphate interactions. The MD technique explains why some potentially isosteric SR motifs are not realized during evolution. We also found that the inability to form stable cWW geometry is an important factor in the case of the first base pair of the flexible region of the SR motif. A comparison of structural, bioinformatics, SHAPE probing, and MD simulation data reveals that explicit solvent MD simulations neatly reflect the viability of different sequence variants of the SR motif. Thus, MD simulations can efficiently complement bioinformatics tools in studies of conservation patterns of RNA motifs and provide atomistic insight into the role of their different signature interactions.

Computational approaches to predicting the impact of novel bases on RNA structure and stability

The use of computational modeling techniques to gain insight into nucleobase interactions has been a challenging endeavor to date. Accurate treatment requires the tackling of many challenges but also holds the promise of great rewards. The development of effective computational approaches to predict the binding affinities of nucleobases and analogues can, for example, streamline the process of developing novel nucleobase modifications, which should facilitate the development of new RNAi-based therapeutics. This brief review focuses on available computational approaches to predicting base pairing affinity in RNA-based contexts such as nucleobase-nucleobase interactions in duplexes and nucleobase-protein interactions. The challenges associated with such modeling along with potential future directions for the field are highlighted.

Higher order structural effects stabilizing the reverse Watson-Crick Guanine-Cytosine base pair in functional RNAs

The G:C reverse Watson-Crick (W:W trans) base pair, also known as Levitt base pair in the context of tRNAs, is a structurally and functionally important base pair that contributes to tertiary interactions joining distant domains in functional RNA molecules and also participates in metabolite binding in riboswitches. We previously indicated that the isolated G:C W:W trans base pair is a rather unstable geometry, and that dicationic metal binding to the Guanine base or posttranscriptional modification of the Guanine can increase its stability. Herein, we extend our survey and report on other H-bonding interactions that can increase the stability of this base pair. To this aim, we performed a bioinformatics search of the PDB to locate all the occurencies of G:C trans base pairs. Interestingly, 66% of the G:C trans base pairs in the PDB are engaged in additional H-bonding interactions with other bases, the RNA backbone or structured water molecules. High level quantum mechanical calculations on a data set of representative crystal structures were performed to shed light on the structural stability and energetics of the various crystallographic motifs. This analysis was extended to the binding of the preQ1 metabolite to a preQ1-II riboswitch.

Molecular modeling of nucleic acid structure: energy and sampling

An overview of computer simulation techniques as applied to nucleic acid systems is presented. This unit discusses methods used to treat energy and to sample representative configurations. Emphasis is placed on molecular mechanics and empirical force fields.

Configuration of wobble base pairs having pyrimidines as anticodon wobble bases: significance for codon degeneracy

Degeneracy of the genetic code was attributed by Crick to imprecise hydrogen-bonded base-pairing at the wobble position during codon-anticodon pairing. The Crick wobble rules define but do not explain the RNA base pair combinations allowed at this position. We select six pyrimidine bases functioning as anticodon wobble bases (AWBs) to study their H-bonded pairing properties with the four major RNA bases using density functional theory at the B3LYP/6-31G(d,p) level. This is done to assess the extent to which the configuration of a solitary RNA wobble base pair may in itself determine specificity and degeneracy of the genetic code by allowing or disallowing the given base pair during codon-anticodon pairing. Calculated values of select configuration markers for the base pairs screen well between allowed and disallowed base pairs for most cases examined here, where the base pair width emerges as an important factor. A few allowed wobble pairs invoke the involvement of RNA nucleoside conformation, as well as involvement of the exocyclic substituent in H-bonding. This study, however, cannot explain the disallowed status of the Ura⋯Gua wobble pair on the basis of configuration alone. Explanation of the allowed status of the V⋯Ura pair requires further study on the mediatory role of water molecules. Apart from these two cases, these computational results are sufficient, on the basis of base pair configuration alone, to account for the specificity and degeneracy of the genetic code for all known cases of codon-anticodon pairing which involve the pyrimidine AWBs studied here.

RNA structure and dynamics: a base pairing perspective

RNA is now known to possess various structural, regulatory and enzymatic functions for survival of cellular organisms. Functional RNA structures are generally created by three-dimensional organization of small structural motifs, formed by base pairing between self-complementary sequences from different parts of the RNA chain. In addition to the canonical Watson-Crick or wobble base pairs, several non-canonical base pairs are found to be crucial to the structural organization of RNA molecules. They appear within different structural motifs and are found to stabilize the molecule through long-range intra-molecular interactions between basic structural motifs like double helices and loops. These base pairs also impart functional variation to the minor groove of A-form RNA helices, thus forming anchoring site for metabolites and ligands. Non-canonical base pairs are formed by edge-to-edge hydrogen bonding interactions between the bases. A large number of theoretical studies have been done to detect and analyze these non-canonical base pairs within crystal or NMR derived structures of different functional RNA. Theoretical studies of these isolated base pairs using ab initio quantum chemical methods as well as molecular dynamics simulations of larger fragments have also established that many of these non-canonical base pairs are as stable as the canonical Watson-Crick base pairs. This review focuses on the various structural aspects of non-canonical base pairs in the organization of RNA molecules and the possible applications of these base pairs in predicting RNA structures with more accuracy.

Relative stability of different DNA guanine quadruplex stem topologies derived using large-scale quantum-chemical computations

We provide theoretical predictions of the intrinsic stability of different arrangements of guanine quadruplex (G-DNA) stems. Most computational studies of nucleic acids have applied Molecular Mechanics (MM) approaches using simple pairwise-additive force fields. The principle limitation of such calculations is the highly approximate nature of the force fields. In this study, we for the first time apply accurate QM computations (DFT-D3 with large atomic orbital basis sets) to essentially complete DNA building blocks, seven different folds of the cation-stabilized two-quartet G-DNA stem, each having more than 250 atoms. The solvent effects are approximated by COSMO continuum solvent. We reveal sizable differences between MM and QM descriptions of relative energies of different G-DNA stems, which apparently reflect approximations of the DNA force field. Using the QM energy data, we propose correction to earlier free energy estimates of relative stabilities of different parallel, hybrid, and antiparallel G-stem folds based on classical simulations. The new energy ranking visibly improves the agreement between theory and experiment. We predict the 5'-anti-anti-3' GpG dinucleotide step to be the most stable one, closely followed by the 5'-syn-anti-3' step. The results are in good agreement with known experimental structures of 2-, 3-, and 4-quartet G-DNA stems. Besides providing specific results for G-DNA, our study highlights basic limitations of force field modeling of nucleic acids. Although QM computations have their own limitations, mainly the lack of conformational sampling and the approximate description of the solvent, they can substantially improve the quality of calculations currently relying exclusively on force fields.

How to understand quantum chemical computations on DNA and RNA systems? A practical guide for non-specialists

In this review primarily written for non-experts we explain basic methodological aspects and interpretation of modern quantum chemical (QM) computations applied to nucleic acids. We introduce current reference QM computations on small model systems consisting of dozens of atoms. Then we comment on recent advance of fast and accurate dispersion-corrected density functional theory methods, which will allow computations of small but complete nucleic acids building blocks in the near future. The qualitative difference between QM and molecular mechanics (MM, force field) computations is discussed. We also explain relation of QM and molecular simulation computations to experiments.

The DNA and RNA sugar-phosphate backbone emerges as the key player. An overview of quantum-chemical, structural biology and simulation studies

Knowledge of geometrical and physico-chemical properties of the sugar-phosphate backbone substantially contributes to the comprehension of the structural dynamics, function and evolution of nucleic acids. We provide a side by side overview of structural biology/bioinformatics, quantum chemical and molecular mechanical/simulation studies of the nucleic acids backbone. We highlight main features, advantages and limitations of these techniques, with a special emphasis given to their synergy. The present status of the research is then illustrated by selected examples which include classification of DNA and RNA backbone families, benchmark structure-energy quantum chemical calculations, parameterization of the dihedral space of simulation force fields, incorporation of arsenate into DNA, sugar-phosphate backbone self-cleavage in small RNA enzymes, and intricate geometries of the backbone in recurrent RNA building blocks. Although not apparent from the current literature showing limited overlaps between the QM, simulation and bioinformatics studies of the nucleic acids backbone, there in fact should be a major cooperative interaction between these three approaches in studies of the sugar-phosphate backbone.

Structure and energy of non-canonical basepairs: comparison of various computational chemistry methods with crystallographic ensembles

Different types of non-canonical basepairs, in addition to the Watson-Crick ones, are observed quite frequently in RNA. Their importance in the three dimensional structure is not fully understood, but their various roles have been proposed by different groups. We have analyzed the energetics and geometry of 32 most frequently observed basepairs in the functional RNA crystal structures using different popular empirical, semi-empirical and ab initio quantum chemical methods and compared their optimized geometry with the crystal data. These basepairs are classified into three categories: polar, non-polar and sugar-mediated, depending on the types of atoms involved in hydrogen bonding. In case of polar basepairs, most of the methods give rise to optimized structures close to their initial geometry. The interaction energies also follow similar trends, with the polar ones having more attractive interaction energies. Some of the C-H...O/N hydrogen bond mediated non-polar basepairs are also found to be significantly stable in terms of their interaction energy values. Few polar basepairs, having amino or carboxyl groups not hydrogen bonded to anything, such as G:G H:W C, show large flexibility. Most of the non-polar basepairs, except A:G s:s T and A:G w:s C, are found to be stable; indicating C-H...O/N interaction also plays a prominent role in stabilizing the basepairs. The sugar mediated basepairs show variability in their structures, due to the involvement of flexible ribose sugar. These presumably indicate that the most of the polar basepairs along with few non-polar ones act as seed for RNA folding while few may act as some conformational switch in the RNA.

Squaring theory with practice in RNA design

Ribonucleic acid (RNA) design offers unique opportunities for engineering genetic networks and nanostructures that self-assemble within living cells. Recent years have seen the creation of increasingly complex RNA devices, including proof-of-concept applications for in vivo three-dimensional scaffolding, imaging, computing, and control of biological behaviors. Expert intuition and simple design rules--the stability of double helices, the modularity of noncanonical RNA motifs, and geometric closure--have enabled these successful applications. Going beyond heuristics, emerging algorithms may enable automated design of RNAs with nucleotide-level accuracy but, as illustrated on a recent RNA square design, are not yet fully predictive. Looking ahead, technological advances in RNA synthesis and interrogation are poised to radically accelerate the discovery and stringent testing of design methods.

Understanding the role of base stacking in nucleic acids. MD and QM analysis of tandem GA base pairs in RNA duplexes

Preceding NMR experiments show that the conformation of tandem GA base pairs, an important recurrent non-canonical building block in RNA duplexes, is context dependent. The GA base pairs adopt "sheared" N3(G)-N6(A), N2(G)-N7(A) geometry in the r(CGAG)(2) and r(iGGAiC)(2) contexts while switching to "imino" N1(G)-N1(A), O6(G)-N6(A) geometry in the r(GGAC)(2) and r(iCGAiG)(2) contexts (iC and iG stand for isocytosine and isoguanine, respectively). As base stacking is likely to be one of the key sources of the context dependence of the conformation of GA base pairs, we calculated base stacking energies in duplexes containing such base pairs, to see if this dependence can be predicted by stacking energy calculations. When investigating the context dependence of the GA geometry two different conformations of the same duplex were compared (imino vs. sheared). The geometries were generated via explicit solvent MD simulations of the respective RNA duplexes, while the subsequent QM energy calculations focused on base stacking interactions of the four internal base pairs. Geometrical relaxation of nucleobase atoms prior to the stacking energy computations has a non-negligible effect on the results. The stacking energies were derived at the DFT-D/6-311++G(3df,3pd) level. We show a rather good correspondence between the intrinsic gas-phase stacking energies and the NMR-determined GA geometries. The conformation with more favorable gas-phase stacking is in most cases the one observed in experiments. This correlation is not improved when including solvent effects via the COSMO method. On the other side, the stacking calculations do not predict the relative thermodynamic stability of duplex formation for different sequences.

Testing the nearest neighbor model for canonical RNA base pairs: revision of GU parameters

Thermodynamic parameters for GU pairs are important for predicting the secondary structures of RNA and for finding genomic sequences that code for structured RNA. Optical melting curves were measured for 29 RNA duplexes with GU pairs to improve nearest neighbor parameters for predicting stabilities of helixes. The updated model eliminates a prior penalty assumed for terminal GU pairs. Six additional duplexes with the 5'GG/3'UU motif were added to the single representation in the previous database. This revises the ΔG°(37) for the 5'GG/3'UU motif from an unfavorable 0.5 kcal/mol to a favorable -0.2 kcal/mol. Similarly, the ΔG°(37) for the 5'UG/3'GU motif changes from 0.3 to -0.6 kcal/mol. The correlation coefficients between predicted and experimental ΔG°(37), ΔH°, and ΔS° for the expanded database are 0.95, 0.89, and 0.87, respectively. The results should improve predictions of RNA secondary structure.

Role of wobble base pair geometry for codon degeneracy: purine-type bases at the anticodon wobble position

Codon degeneracy is a key feature of the genetic code, explained by Crick (J Mol Biol 19:548-555, 1966) in terms of imprecision of base pairing at the codon third position (the wobble position) of the codon-anticodon duplex. The Crick wobble rules define, but do not explain, which base pairs are allowed/disallowed at the wobble position of this duplex. This work examines whether the H-bonded configurations of solitary RNA base pairs can in themselves help decide which base pairs are allowed at the wobble position during codon-anticodon pairing. Taking the purine-type bases guanine, hypoxanthine, queuine and adenine as anticodon wobble bases, H-bonded pairing energies and optimized configurations of numerous RNA base pairs are calculated in gas and modeled aqueous phase at the B3LYP/6-31 G(d,p) level. Calculated descriptors of alignment of these solitary base pairs are able to screen between allowed and disallowed base pairs for all cases studied here, except two cases which invoke base-sugar interactions in the codon wobble nucleoside. The exclusion of adenine from the anticodon wobble position cannot be explained on the basis of pairing facility or base pair geometry. These DFT results thus account for the specificity and degeneracy of the genetic code for all cases involving guanine, hypoxanthine and queuine as anticodon wobble bases.

Understanding the Sequence Preference of Recurrent RNA Building Blocks using Quantum Chemistry: The Intrastrand RNA Dinucleotide Platform

Folded RNA molecules are shaped by an astonishing variety of highly conserved noncanonical molecular interactions and backbone topologies. The dinucleotide platform is a widespread recurrent RNA modular building submotif formed by the side-by-side pairing of bases from two consecutive nucleotides within a single strand, with highly specific sequence preferences. This unique arrangement of bases is cemented by an intricate network of noncanonical hydrogen bonds and facilitated by a distinctive backbone topology. The present study investigates the gas-phase intrinsic stabilities of the three most common RNA dinucleotide platforms - 5'-GpU-3', ApA, and UpC - via state-of-the-art quantum-chemical (QM) techniques. The mean stability of base-base interactions decreases with sequence in the order GpU > ApA > UpC. Bader's atoms-in-molecules analysis reveals that the N2(G)…O4(U) hydrogen bond of the GpU platform is stronger than the corresponding hydrogen bonds in the other two platforms. The mixed-pucker sugar-phosphate backbone conformation found in most GpU platforms, in which the 5'-ribose sugar (G) is in the C2'-endo form and the 3'-sugar (U) in the C3'-endo form, is intrinsically more stable than the standard A-RNA backbone arrangement, partially as a result of a favorable O2'…O2P intra-platform interaction. Our results thus validate the hypothesis of Lu et al. (Lu Xiang-Jun, et al. Nucleic Acids Res. 2010, 38, 4868-4876), that the superior stability of GpU platforms is partially mediated by the strong O2'…O2P hydrogen bond. In contrast, ApA and especially UpC platform-compatible backbone conformations are rather diverse and do not display any characteristic structural features. The average stabilities of ApA and UpC derived backbone conformers are also lower than those of GpU platforms. Thus, the observed structural and evolutionary patterns of the dinucleotide platforms can be accounted for, to a large extent, by their intrinsic properties as described by modern QM calculations. In contrast, we show that the dinucleotide platform is not properly described in the course of atomistic explicit-solvent simulations. Our work also gives methodological insights into QM calculations of experimental RNA backbone geometries. Such calculations are inherently complicated by rather large data and refinement uncertainties in the available RNA experimental structures, which often preclude reliable energy computations.

Theoretical calculation of the OH vibrational overtone spectra of 1,5-pentanediol and 1,6-hexanediol

It is well-known that intramolecular hydrogen bonding affects the relative energetics of conformers, as well as the OH stretching peak positions, intensities, and width. In this study we simulated the Δv(OH) = 3, 4 overtone spectra of 1,5-pentanediol (PeD) and 1,6-hexanediol (HD) using the peak positions, intensities, and width calculated from the B3LYP/6-31+G(d,p) method. Furthermore, room temperature free energy calculations were performed using B3LYP/6-31+G(d,p) MP2/6-31+G(d,p), and MP2/6-311++G(3df,3pd) to obtain the relative population of the conformers. From the calculation of 109 and 381 distinct conformers for PeD and HD, respectively, we find that for these long chain diols the intramolecular hydrogen bonded conformers are not the most dominant conformation at room temperature. This is in stark contrast with shorter chain diols such as ethylene glycol for which the hydrogen bonded conformer dominates the population at room temperature. On the other hand, we found that the correlation between the hydrogen bonded OH red shift versus the homogeneous width, Γ = 0.0155(Δω)(1.36), which was derived for shorter chain diols, is valid even for these longer chain diols. We also showed that the intramolecular hydrogen bonded OH initially decays through the CCOH torsion and COH bending mode no matter how long the alkanediol chain length is for 1,n-alkanediols for n up to 6.

A-minor tertiary interactions in RNA kink-turns. Molecular dynamics and quantum chemical analysis

The RNA kink-turn is an important recurrent RNA motif, an internal loop with characteristic consensus sequence forming highly conserved three-dimensional structure. Functional arrangement of RNA kink-turns shows a sharp bend in the phosphodiester backbone. Among other signature interactions, kink-turns form A-minor interaction between their two stems. Most kink-turns possess extended A-minor I (A-I) interaction where adenine of the second A•G base pair of the NC-stem interacts with the first canonical pair of the C-stem (i.e., the receptor pair) via trans-sugar-edge/sugar-edge (tSS) and cis-sugar-edge/sugar-edge (cSS) interactions. The remaining kink-turns have less compact A-minor 0 (A-0) interaction with just one tSS contact. We show that kink-turns with A-I in ribosomal X-ray structures keep G═C receptor base pair during evolution while the inverted pair (C═G) is not realized. In contrast, kink-turns with A-0 in the observed structures alternate G═C and C═G base pairs in sequences. We carried out an extended set (~5 μs) of explicit-solvent molecular dynamics simulations of kink-turns to rationalize this structural/evolutionary pattern. The simulations were done using a net-neutral Na(+) cation atmosphere (with ~0.25 M cation concentration) supplemented by simulations with either excess salt KCl atmosphere or inclusion of Mg(2+). The results do not seem to depend on the treatment of ions. The simulations started with X-ray structures of several kink-turns while we tested the response of the simulated system to base substitutions, modest structural perturbations and constraints. The trends seen in the simulations reveal that the A-I/G═C arrangement is preferred over all three other structures. The A-I/C═G triple appears structurally entirely unstable, consistent with the covariation patterns seen during the evolution. The A-0 arrangements tend to shift toward the A-I pattern in simulations, which suggests that formation of the A-0 interaction is likely supported by the surrounding protein and RNA molecules. A-0 may also be stabilized by additional kink-turn nucleotides not belonging to the kink-turn consensus, as shown for the kink-turn from ribosomal Helix 15. Quantum-chemical calculations on all four A-minor triples suggest that there is a different balance of electrostatic and dispersion stabilization in the A-I/G═C and A-I/C═G triples, which may explain different behavior of these otherwise isosteric triples in the context of kink-turns.

Noncanonical hydrogen bonding in nucleic acids. Benchmark evaluation of key base-phosphate interactions in folded RNA molecules using quantum-chemical calculations and molecular dynamics simulations

RNA molecules are stabilized by a wide range of noncanonical interactions that are not present in DNA. Among them, the recently classified base-phosphate (BPh) interactions belong to the most important ones. Twelve percent of nucleotides in the ribosomal crystal structures are involved in BPh interactions. BPh interactions are highly conserved and provide major constraints on RNA sequence evolution. Here we provide assessment of the energetics of BPh interactions using MP2 computations extrapolated to the complete basis set of atomic orbitals and corrected for higher-order electron correlation effects. The reference computations are compared with DFT-D and DFT-D3 approaches, the SAPT method, and the molecular mechanics force field. The computations, besides providing the basic benchmark for the BPh interactions, allow some refinements of the original classification, including identification of some potential doubly bonded BPh patterns. The reference computations are followed by analysis of some larger RNA fragments that consider the context of the BPh interactions. The computations demonstrate the complexity of interaction patterns utilizing the BPh interactions in real RNA structures. The BPh interactions are often involved in intricate interaction networks. We studied BPh interactions of protonated adenine that can contribute to catalysis of hairpin ribozyme, the key BPh interaction in the S-turn motif of the sarcin-ricin loop, which may predetermine the S-turn topology and complex BPh patterns from the glmS riboswitch. Finally, the structural stability of BPh interactions in explicit solvent molecular dynamics simulations is assessed. The simulations well preserve key BPh interactions and allow dissection of structurally/functionally important water-meditated BPh bridges, which could not be considered in earlier bioinformatics classification of BPh interactions.

Benchmarking AMBER force fields for RNA: comparisons to NMR spectra for single-stranded r(GACC) are improved by revised χ torsions

Accurately modeling unpaired regions of RNA is important for predicting structure, dynamics, and thermodynamics of folded RNA. Comparisons between NMR data and molecular dynamics simulations provide a test of force fields used for modeling. Here, NMR spectroscopy, including NOESY, (1)H-(31)P HETCOR, DQF-COSY, and TOCSY, was used to determine conformational preferences for single-stranded GACC RNA. The spectra are consistent with a conformational ensemble containing major and minor A-form-like structures. In a series of 50 ns molecular dynamics (MD) simulations with the AMBER99 force field in explicit solvent, initial A-form-like structures rapidly evolve to disordered conformations. A set of 50 ns simulations with revised χ torsions (AMBER99χ force field) gives two primary conformations, consistent with the NMR spectra. A single 1.9 μs MD simulation with the AMBER99χ force field showed that the major and minor conformations are retained for almost 68% of the time in the first 700 ns, with multiple transformations from A-form to non-A-form conformations. For the rest of the simulation, random-coil structures and a stable non-A-form conformation inconsistent with NMR spectra were seen. Evidently, the AMBER99χ force field improves structural predictions for single-stranded GACC RNA compared to the AMBER99 force field, but further force field improvements are needed.

Molecular modeling of nucleic acid structure: energy and sampling

An overview of computer simulation techniques as applied to nucleic acid systems is presented. This unit expands an accompanying overview unit (UNIT 7.5) by discussing methods used to treat the energy and sample representative configurations. Emphasis is placed on molecular mechanics and empirical force fields.