Molecular dynamics simulations of RNA: an in silico single molecule approach

Influence of ion and hydration atmospheres on RNA structure and dynamics: insights from advanced theoretical and computational methods

RNA, a highly charged biopolymer composed of negatively charged phosphate groups, defies electrostatic repulsion to adopt well-defined, compact structures. Hence, the presence of positively charged metal ions is crucial not only for RNA's charge neutralization, but they also coherently decorate the ion atmosphere of RNA to stabilize its compact fold. This feature article elucidates various modes of close RNA-ion interactions, with a special emphasis on Mg2+ as an outer-sphere and inner-sphere ion. Through examples, we highlight how inner-sphere chelated Mg2+ stabilizes RNA pseudoknots, while outer-sphere ions can also exert long-range electrostatic interactions, inducing groove narrowing, coaxial helical stacking, and RNA ring formation. In addition to investigating the RNA's ion environment, we note that the RNA's hydration environment is relatively underexplored. Our study delves into its profound interplay with the structural dynamics of RNA, employing state-of-the-art atomistic simulation techniques. Through examples, we illustrate how specific ions and water molecules are associated with RNA functions, leveraging atomistic simulations to identify preferential ion binding and hydration sites. However, understanding their impact(s) on the RNA structure remains challenging due to the involvement of large length and long time scales associated with RNA's dynamic nature. Nevertheless, our contributions and recent advances in coarse-grained simulation techniques offer insights into large-scale structural changes dynamically linked to the RNA ion atmosphere. In this connection, we also review how different cutting-edge computational simulation methods provide a microscopic lens into the influence of ions and hydration on RNA structure and dynamics, elucidating distinct ion atmospheric components and specific hydration layers and their individual and collective impacts.

Computer-Assisted Design and Characterization of RNA Nanostructures

Molecular dynamics (MD) simulations can aid in the design and characterization of RNA nanomaterials, providing details about structural and dynamical properties as a function of sequence and environment. Here, we describe how to perform explicit and implicit solvent all-atom MD simulations for RNA nanoring systems.

Dynamic effects of the spine of hydrated magnesium on viral RNA pseudoknot structure

In the cellular environment, a viral RNA Pseudoknot (PK) structure is responsive to its prevailing ion atmosphere created by a mixture of monovalent and divalent cations. We investigate the influence of such a mixed-salt environment on RNA-PK structure at an atomic resolution through three sets of 1.5 μs explicit solvent molecular dynamics (MD) simulations and also by building a dynamic counterion-condensation (DCC) model at varying divalent Mg²⁺ concentrations. The DCC model includes explicit interaction of the ligand and adjacent chelated Mg²⁺ by extending the recently developed generalized Manning condensation model. Its implementation within an all-atom structure-based molecular dynamics framework bolsters its opportunity to explore large-length scale and long-timescale phenomena associated with complex RNA systems immersed in its dynamic ion environment. In the present case of RNA-PK, both explicit MD and DCC simulations reveal a spine of hydrated-Mg²⁺ to induce stem-I and stem-II closure where the minor groove between these stems is akin to breathing. Mg²⁺ mediated minor-groove narrowing is coupled with local base-flipping dynamics of a base triple and quadruple, changing the stem structure of such RNA. Contrary to the conversational view of the indispensable need for Mg²⁺ for the tertiary structure of RNA, the study warns about the plausible detrimental effect of specific Mg²⁺-phosphate interactions on the RNA-PK structure beyond a certain concentration of Mg²⁺.

Petkov P, Milenov TI, Vayssilov GN. Sodium and Magnesium Ion Location at the Backbone and at the Nucleobase of RNA: Ab Initio Molecular Dynamics in Water Solution

The interactions between Na⁺ or Mg²⁺ ions with different parts of single-stranded RNA molecules, namely, the oxygen atoms from the phosphate groups or the guanine base, in water solution have been studied using first-principles molecular dynamics. Sodium ions were found to be much more mobile than Mg²⁺ ions and readily underwent transitions between a state directly bonded to RNA oxygen atoms and a completely solvated state. The inner solvation shell of Na⁺ ions fluctuated stochastically at a femtosecond timescale coordinating on average 5 oxygen atoms for bonded Na⁺ ions and 5.5 oxygen atoms for solvated Na⁺ ions. In contrast, the inner solvation shell of Mg²⁺ ions was stable in both RNA-bonded and completely solvated states. In both cases, Mg²⁺ ions coordinated 6 oxygen atoms from the inner solvation shell. Consistent with their stable solvation shells, Mg²⁺ ions were more effective than Na⁺ ions in stabilizing the RNA backbone conformation. The exclusion zones between the first and second solvation shells, solvation shell widths, and angles for binding to carbonyl oxygen of guanine for solvated Na⁺ or Mg²⁺ ions exhibited a number of quantitative differences when compared with RNA crystallographic data. The presented results support the distinct capacity of Mg²⁺ ions to support the RNA structure not only in the crystal phase but also in the dynamic water environment both on the side of the phosphate moiety and on the side of the nucleobase.

Algorithm to minimize MPI communications in the parallelized fast multipole method combined with molecular dynamics calculations

In the era of exascale supercomputers, large-scale, and long-time molecular dynamics (MD) calculations are expected to make breakthroughs in various fields of science and technology. Here, we propose a new algorithm to improve the parallelization performance of message passing interface (MPI)-communication in the MPI-parallelized fast multipole method (FMM) combined with MD calculations under three-dimensional periodic boundary conditions. Our approach enables a drastic reduction in the amount of communication data, including the atomic coordinates and multipole coefficients, both of which are required to calculate the electrostatic interaction by using the FMM. In communications of multipole coefficients, the reduction rate of communication data in the new algorithm relative to the amount of data in the conventional one increases as both the number of FMM levels and the number of MPI processes increase. The aforementioned rate increase could exceed 50% as the number of MPI processes becomes larger for very large systems. The proposed algorithm, named the minimum-transferred data (MTD) method, should enable large-scale and long-time MD calculations to be calculated efficiently, under the condition of massive MPI-parallelization on exascale supercomputers.

Molecular dynamics simulation of the binding process of ligands to the add adenine riboswitch aptamer

Riboswitches are RNA-structured elements that modulate gene expression through changing their conformation in response to specific metabolite binding. However, the regulation mechanisms of riboswitches by ligand binding are still not well understood. At present two possible ligand-binding mechanisms have been proposed: conformational selection and induced fit. Based on explicit-solvent molecular dynamics (MD) simulations, we have studied the process of the binding of ligands (adenines) to add adenine riboswitch aptamer (AARA) in detail. Our results show that the relative high flexibility of the junction J23 of AARA allows the ligands to be captured by the binding pocket of AARA in a near-native state, which may be driven by hydrophobic and base-stacking interactions. In addition, the binding of a ligand makes the stem P1 and the junction J23 of AARA more stable than in the absence of the ligand. Generally, our analyses show that the ligand-binding process of the add adenine riboswitch may be partially embodied by a conformational selection mechanism.

Recent Advances in Machine Learning Based Prediction of RNA-protein Interactions

The interactions between RNAs and proteins play critical roles in many biological processes. Therefore, characterizing these interactions becomes critical for mechanistic, biomedical, and clinical studies. Many experimental methods can be used to determine RNA-protein interactions in multiple aspects. However, due to the facts that RNA-protein interactions are tissuespecific and condition-specific, as well as these interactions are weak and frequently compete with each other, those experimental techniques can not be made full use of to discover the complete spectrum of RNA-protein interactions. To moderate these issues, continuous efforts have been devoted to developing high quality computational techniques to study the interactions between RNAs and proteins. Many important progresses have been achieved with the application of novel techniques and strategies, such as machine learning techniques. Especially, with the development and application of CLIP techniques, more and more experimental data on RNA-protein interaction under specific biological conditions are available. These CLIP data altogether provide a rich source for developing advanced machine learning predictors. In this review, recent progresses on computational predictors for RNA-protein interaction were summarized in the following aspects: dataset, prediction strategies, and input features. Possible future developments were also discussed at the end of the review.

Protein Interaction with Charged Macromolecules: From Model Polymers to Unfolded Proteins and Post-Translational Modifications

Interaction of proteins with charged macromolecules is involved in many processes in cells. Firstly, there are many naturally occurred charged polymers such as DNA and RNA, polyphosphates, sulfated glycosaminoglycans, etc., as well as pronouncedly charged proteins such as histones or actin. Electrostatic interactions are also important for "generic" proteins, which are not generally considered as polyanions or polycations. Finally, protein behavior can be altered due to post-translational modifications such as phosphorylation, sulfation, and glycation, which change a local charge of the protein region. Herein we review molecular modeling for the investigation of such interactions, from model polyanions and polycations to unfolded proteins. We will show that electrostatic interactions are ubiquitous, and molecular dynamics simulations provide an outstanding opportunity to look inside binding and reveal the contribution of electrostatic interactions. Since a molecular dynamics simulation is only a model, we will comprehensively consider its relationship with the experimental data.

Dynamics of metal ions around an RNA molecule

Noncoding RNA molecules take part in many biological processes, while metal ions play crucial roles in helping RNAs to perform their functions. However, the statics and dynamics of these metal ions around RNA molecules are still not well understood. In this work, we report a detailed molecular dynamics study of the type-I preQ_{1}-bound riboswitch aptamer domain (PRAD) at different ionic conditions (K^{+}, Na^{+}, and Mg^{2+}). The results show that the structural properties and flexibility of the PRAD molecule greatly influence the distributions and dynamics of metal ions around it. Simultaneously, Na^{+} ions show a stronger competitiveness with Mg^{2+} ions than K^{+} ions, and the three types of metal ions have different modes of interaction with the RNA molecule. Furthermore, we have also investigated specific binding sites of metal ions on the PRAD molecule and found that the dynamics and hydration structures of metal ions located at the ion-binding sites were obviously affected by the RNA structure near these ion-binding sites. These results may be useful to understand the role of the metal ions in noncoding RNA functions.

RNA Structural Differentiation: Opportunities with Pattern Recognition

Our awareness and appreciation of the many regulatory roles of RNA have dramatically increased in the past decade. This understanding, in addition to the impact of RNA in many disease states, has renewed interest in developing selective RNA-targeted small molecule probes. However, the fundamental guiding principles in RNA molecular recognition that could accelerate these efforts remain elusive. While high-resolution structural characterization can provide invaluable insight, examples of well-characterized RNA structures, not to mention small molecule:RNA complexes, remain limited. This Perspective provides an overview of the current techniques used to understand RNA molecular recognition when high-resolution structural information is unavailable. We will place particular emphasis on a new method, pattern recognition of RNA with small molecules (PRRSM), that provides rapid insight into critical components of RNA recognition and differentiation by small molecules as well as into RNA structural features.

Atomistic Analysis of ToxN and ToxI Complex Unbinding Mechanism

ToxIN is a triangular structure formed by three protein toxins (ToxNs) and three specific noncoding RNA antitoxins (ToxIs). To respond to stimuli, ToxI is preferentially degraded, releasing the ToxN. Thus, the dynamic character is essential in the normal function interactions between ToxN and ToxI. Here, equilibrated molecular dynamics (MD) simulations were performed to study the stability of ToxN and ToxI. The results indicate that ToxI adjusts the conformation of 3' and 5' termini to bind to ToxN. Steered molecular dynamics (SMD) simulations combined with the recently developed thermodynamic integration in 3nD (TI3nD) method were carried out to investigate ToxN unbinding from the ToxIN complex. The potentials of mean force (PMFs) and atomistic pictures suggest the unbinding mechanism as follows: (1) dissociation of the 5' terminus from ToxN, (2) missing the interactions involved in the 3' terminus of ToxI without three nucleotides (G31, A32, and A33), (3) starting to unfold for ToxI, (4) leaving the binding package of ToxN for three nucleotides of ToxI, (5) unfolding of ToxI. This work provides information on the structure-function relationship at the atomistic level, which is helpful for designing new potent antibacterial drugs in the future.

Insights into the effect of molecular crowding on the structure, interactions and functions of siRNA-PAZ complex through molecular dynamics studies

siRNA molecules are well known to be involved in the post-transcriptional regulation of gene expression and play a key role in understanding the intricacies of eukaryotic gene regulation. While it is widely known that 3' end of siRNA binds to the PAZ domain of Argonaute proteins, it remains unclear whether the molecular crowding facilitates or hinders the overall siRNA-protein interactions during RNA interference. The biological interaction networks controlling the cellular functions of any biological cell may behave very differently in crowded environment as compared to the dilute conditions. Therefore, it is of interest to study the siRNA-protein interactions under more physiologically relevant conditions. In our previous work, we studied the role of hydrogen bond interactions and water network interactions towards the structural integrity of siRNA-PAZ complex. We also described the motions relevant for the functioning of the complex and analyzed the biphasic interaction of the 3' end of siRNA within the PAZ domain under aqueous condition. In the present work, we studied the dynamics of siRNA-PAZ complex in the presence of water-soluble crowding agent. We observed significant changes in interactions and dynamics of siRNA-PAZ complex in the presence of crowder. The conserved H-bond interactions were destabilized by ≈12%, while interfacial water networks were stabilized in the presence of crowder. The analysis of siRNA-PAZ dynamics revealed the stabilizing role of ARG52, ARG172, LYS53 and LYS173 in siRNA-PAZ complex. Interestingly, despite increase in flexibility as measured by RMSD, crowding stabilized top modes. Communicated by Ramaswamy H. Sarma.

Above the Epitranscriptome: RNA Modifications and Stem Cell Identity

Sequence databases and transcriptome-wide mapping have revealed different reversible and dynamic chemical modifications of the nitrogen bases of RNA molecules. Modifications occur in coding RNAs and noncoding-RNAs post-transcriptionally and they can influence the RNA structure, metabolism, and function. The result is the expansion of the variety of the transcriptome. In fact, depending on the type of modification, RNA molecules enter into a specific program exerting the role of the player or/and the target in biological and pathological processes. Many research groups are exploring the role of RNA modifications (alias epitranscriptome) in cell proliferation, survival, and in more specialized activities. More recently, the role of RNA modifications has been also explored in stem cell biology. Our understanding in this context is still in its infancy. Available evidence addresses the role of RNA modifications in self-renewal, commitment, and differentiation processes of stem cells. In this review, we will focus on five epitranscriptomic marks: N6-methyladenosine, N1-methyladenosine, 5-methylcytosine, Pseudouridine (Ψ) and Adenosine-to-Inosine editing. We will provide insights into the function and the distribution of these chemical modifications in coding RNAs and noncoding-RNAs. Mainly, we will emphasize the role of epitranscriptomic mechanisms in the biology of naïve, primed, embryonic, adult, and cancer stem cells.

Intrinsic Base-Pair Rearrangement in the Hairpin Ribozyme Directs RNA Conformational Sampling and Tertiary Interface Formation

Dynamic fluctuations in RNA structure enable conformational changes that are required for catalysis and recognition. In the hairpin ribozyme, the catalytically active structure is formed as an intricate tertiary interface between two RNA internal loops. Substantial alterations in the structure of each loop are observed upon interface formation, or docking. The very slow on-rate for this relatively tight interaction has led us to hypothesize a double conformational capture mechanism for RNA-RNA recognition. We used extensive molecular dynamics simulations to assess conformational sampling in the undocked form of the loop domain containing the scissile phosphate (loop A). We observed several major accessible conformations with distinctive patterns of hydrogen bonding and base stacking interactions in the active-site internal loop. Several important conformational features characteristic of the docked state were observed in well-populated substates, consistent with the kinetic sampling of docking-competent states by isolated loop A. Our observations suggest a hybrid or multistage binding mechanism, in which initial conformational selection of a docking-competent state is followed by induced-fit adjustment to an in-line, chemically reactive state only after formation of the initial complex with loop B.

Coarse-grained modeling of RNA 3D structure

Functional RNA molecules depend on three-dimensional (3D) structures to carry out their tasks within the cell. Understanding how these molecules interact to carry out their biological roles requires a detailed knowledge of RNA 3D structure and dynamics as well as thermodynamics, which strongly governs the folding of RNA and RNA-RNA interactions as well as a host of other interactions within the cellular environment. Experimental determination of these properties is difficult, and various computational methods have been developed to model the folding of RNA 3D structures and their interactions with other molecules. However, computational methods also have their limitations, especially when the biological effects demand computation of the dynamics beyond a few hundred nanoseconds. For the researcher confronted with such challenges, a more amenable approach is to resort to coarse-grained modeling to reduce the number of data points and computational demand to a more tractable size, while sacrificing as little critical information as possible. This review presents an introduction to the topic of coarse-grained modeling of RNA 3D structures and dynamics, covering both high- and low-resolution strategies. We discuss how physics-based approaches compare with knowledge based methods that rely on databases of information. In the course of this review, we discuss important aspects in the reasoning process behind building different models and the goals and pitfalls that can result.

Identification of possible siRNA molecules for TDP43 mutants causing amyotrophic lateral sclerosis: In silico design and molecular dynamics study

The DNA binding protein, TDP43 is a major protein involved in amyotrophic lateral sclerosis and other neurological disorders such as frontotemporal dementia, Alzheimer disease, etc. In the present study, we have designed possible siRNAs for the glycine rich region of tardbp mutants causing ALS disorder based on a systematic theoretical approach including (i) identification of respective codons for all mutants (reported at the protein level) based on both minimum free energy and probabilistic approaches, (ii) rational design of siRNA, (iii) secondary structure analysis for the target accessibility of siRNA, (iii) determination of the ability of siRNA to interact with mRNA and the formation/stability of duplex via molecular dynamics study for a period of 15ns and (iv) characterization of mRNA-siRNA duplex stability based on thermo-physical analysis. The stable GC-rich siRNA expressed strong binding affinity towards mRNA and forms stable duplex in A-form. The linear dependence between the thermo-physical parameters such as Tm, GC content and binding free energy revealed the ability of the identified siRNAs to interact with mRNA in comparable to that of the experimentally reported siRNAs. Hence, this present study proposes few siRNAs as the possible gene silencing agents in RNAi therapy based on the in silico approach.

Sodium and Potassium Interactions with Nucleic Acids

Metal ions are essential cofactors for the structure and functions of nucleic acids. Yet, the early discovery in the 70s of the crucial role of Mg(2+) in stabilizing tRNA structures has occulted for a long time the importance of monovalent cations. Renewed interest in these ions was brought in the late 90s by the discovery of specific potassium metal ions in the core of a group I intron. Their importance in nucleic acid folding and catalytic activity is now well established. However, detection of K(+) and Na(+) ions is notoriously problematic and the question about their specificity is recurrent. Here we review the different methods that can be used to detect K(+) and Na(+) ions in nucleic acid structures such as X-ray crystallography, nuclear magnetic resonance or molecular dynamics simulations. We also discuss specific versus non-specific binding to different structures through various examples.

Prediction of solution properties and dynamics of RNAs by means of Brownian dynamics simulation of coarse-grained models: Ribosomal 5S RNA and phenylalanine transfer RNA

Background

The possibility of validating biological macromolecules with locally disordered domains like RNA against solution properties is helpful to understand their function. In this work, we present a computational scheme for predicting global properties and mimicking the internal dynamics of RNA molecules in solution. A simple coarse-grained model with one bead per nucleotide and two types of intra-molecular interactions (elastic interactions and excluded volume interactions) is used to represent the RNA chain. The elastic interactions are modeled by a set of Hooke springs that form a minimalist elastic network. The Brownian dynamics technique is employed to simulate the time evolution of the RNA conformations.

Results

That scheme is applied to the 5S ribosomal RNA of E. Coli and the yeast phenylalanine transfer RNA. From the Brownian trajectory, several solution properties (radius of gyration, translational diffusion coefficient, and a rotational relaxation time) are calculated. For the case of yeast phenylalanine transfer RNA, the time evolution and the probability distribution of the inter-arm angle is also computed.

Conclusions

The general good agreement between our results and some experimental data indicates that the model is able to capture the tertiary structure of RNA in solution. Our simulation results also compare quite well with other numerical data. An advantage of the scheme described here is the possibility of visualizing the real time macromolecular dynamics.

Electronic supplementary material

The online version of this article (doi:10.1186/s13628-015-0025-7) contains supplementary material, which is available to authorized users.

Wobble pairs of the HDV ribozyme play specific roles in stabilization of active site dynamics

The hepatitis delta virus (HDV) is the only known human pathogen whose genome contains a catalytic RNA motif (ribozyme). The overall architecture of the HDV ribozyme is that of a double-nested pseudoknot, with two GU pairs flanking the active site. Although extensive studies have shown that mutation of either wobble results in decreased catalytic activity, little work has focused on linking these mutations to specific structural effects on catalytic fitness. Here we use molecular dynamics simulations based on an activated structure to probe the active site dynamics as a result of wobble pair mutations. In both wild-type and mutant ribozymes, the in-line fitness of the active site (as a measure of catalytic proficiency) strongly depends on the presence of a C75(N3H3+)N1(O5') hydrogen bond, which positions C75 as the general acid for the reaction. Our mutational analyses show that each GU wobble supports catalytically fit conformations in distinct ways; the reverse G25U20 wobble promotes high in-line fitness, high occupancy of the C75(N3H3+)G1(O5') general-acid hydrogen bond and stabilization of the G1U37 wobble, while the G1U37 wobble acts more locally by stabilizing high in-line fitness and the C75(N3H3+)G1(O5') hydrogen bond. We also find that stable type I A-minor and P1.1 hydrogen bonding above and below the active site, respectively, prevent local structural disorder from spreading and disrupting global conformation. Taken together, our results define specific, often redundant architectural roles for several structural motifs of the HDV ribozyme active site, expanding the known roles of these motifs within all HDV-like ribozymes and other structured RNAs.

Molecular view of ligands specificity for CAG repeats in anti-Huntington therapy

Huntington's disease is a fatal and devastating neurodegenerative genetic disorder for which there is currently no cure. It is characterized by Huntingtin protein's mRNA transcripts with 36 or more CAG repeats. Inhibiting the formation of pathological complexes between these expanded transcripts and target proteins may be a valuable strategy against the disease. Yet, the rational design of molecules specifically targeting the expanded CAG repeats is limited by the lack of structural information. Here, we use well-tempered metadynamics-based free energy calculations to investigate pose and affinity of two ligands targeting CAG repeats for which affinities have been previously measured. The first consists of two 4-guanidinophenyl rings linked by an ester group. It is the most potent ligand identified so far, with Kd = 60(30) nM. The second consists of a 4-phenyl dihydroimidazole and 4-1H-indole dihydroimidazole connected by a C-C bond (Kd = 700(80) nM). Our calculations reproduce the experimental affinities and uncover the recognition pattern between ligands' and their RNA target. They also provide a molecular basis for the markedly different affinity of the two ligands for CAG repeats as observed experimentally. These findings may pave the way for a structure-based hit-to-lead optimization to further improve ligand selectivity toward CAG repeat-containing mRNAs.

Conformational dynamics of bacterial and human cytoplasmic models of the ribosomal A-site

The aminoacyl-tRNA binding site (A-site) is located in helix 44 of small ribosomal subunit. The mobile adenines 1492 and 1493 (Escherichia coli numbering), forming the A-site bulge, act as a functional switch that ensures mRNA decoding accuracy. Structural data on the oligonucleotide models mimicking the ribosomal A-site with sequences corresponding to bacterial and human cytoplasmic sites confirm that this RNA motif forms also without the ribosome context. We performed all-atom molecular dynamics simulations of these crystallographic A-site models to compare their conformational properties. We found that the human A-site bulge is more internally flexible than the bacterial one and has different base pairing preferences, which result in the overall different shapes of these bulges and cation density distributions. Also, in the human A-site model we observed repetitive destacking of A1492, while A1493 was more stably paired than in the bacterial variant. Based on the dynamics of the A-sites we suggest why aminoglycoside antibiotics, which target the bacterial A-site, have lower binding affinities and anti-translational activities toward the human variant.

Computational methods toward accurate RNA structure prediction using coarse-grained and all-atom models

Computational methods can provide significant insights into RNA structure and dynamics, bridging the gap in our understanding of the relationship between structure and biological function. Simulations enrich and enhance our understanding of data derived on the bench, as well as provide feasible alternatives to costly or technically challenging experiments. Coarse-grained computational models of RNA are especially important in this regard, as they allow analysis of events occurring in timescales relevant to RNA biological function, which are inaccessible through experimental methods alone. We have developed a three-bead coarse-grained model of RNA for discrete molecular dynamics simulations. This model is efficient in de novo prediction of short RNA tertiary structure, starting from RNA primary sequences of less than 50 nucleotides. To complement this model, we have incorporated additional base-pairing constraints and have developed a bias potential reliant on data obtained from hydroxyl probing experiments that guide RNA folding to its correct state. By introducing experimentally derived constraints to our computer simulations, we are able to make reliable predictions of RNA tertiary structures up to a few hundred nucleotides. Our refined model exemplifies a valuable benefit achieved through integration of computation and experimental methods.

Multiscale methods for computational RNA enzymology

RNA catalysis is of fundamental importance to biology and yet remains ill-understood due to its complex nature. The multidimensional "problem space" of RNA catalysis includes both local and global conformational rearrangements, changes in the ion atmosphere around nucleic acids and metal ion binding, dependence on potentially correlated protonation states of key residues, and bond breaking/forming in the chemical steps of the reaction. The goal of this chapter is to summarize and apply multiscale modeling methods in an effort to target the different parts of the RNA catalysis problem space while also addressing the limitations and pitfalls of these methods. Classical molecular dynamics simulations, reference interaction site model calculations, constant pH molecular dynamics (CpHMD) simulations, Hamiltonian replica exchange molecular dynamics, and quantum mechanical/molecular mechanical simulations will be discussed in the context of the study of RNA backbone cleavage transesterification. This reaction is catalyzed by both RNA and protein enzymes, and here we examine the different mechanistic strategies taken by the hepatitis delta virus ribozyme and RNase A.

The purine riboswitch as a model system for exploring RNA biology and chemistry

Over the past decade the purine riboswitch, and in particular its nucleobase-binding aptamer domain, has emerged as an important model system for exploring various aspects of RNA structure and function. Its relatively small size, structural simplicity and readily observable activity enable application of a wide variety of experimental approaches towards the study of this RNA. These analyses have yielded important insights into small molecule recognition, co-transcriptional folding and secondary structural switching, and conformational dynamics that serve as a paradigm for other RNAs. In this article, the current state of understanding of the purine riboswitch family and how this growing knowledge base is starting to be exploited in the creation of novel RNA devices are examined. This article is part of a Special Issue entitled: Riboswitches.

Interconversion between parallel and antiparallel conformations of a 4H RNA junction in domain 3 of foot-and-mouth disease virus IRES captured by dynamics simulations

RNA junctions are common secondary structural elements present in a wide range of RNA species. They play crucial roles in directing the overall folding of RNA molecules as well as in a variety of biological functions. In particular, there has been great interest in the dynamics of RNA junctions, including conformational pathways of fully base-paired 4-way (4H) RNA junctions. In such constructs, all nucleotides participate in one of the four double-stranded stem regions, with no connecting loops. Dynamical aspects of these 4H RNAs are interesting because frequent interchanges between parallel and antiparallel conformations are thought to occur without binding of other factors. Gel electrophoresis and single-molecule fluorescence resonance energy transfer experiments have suggested two possible pathways: one involves a helical rearrangement via disruption of coaxial stacking, and the other occurs by a rotation between the helical axes of coaxially stacked conformers. Employing molecular dynamics simulations, we explore this conformational variability in a 4H junction derived from domain 3 of the foot-and-mouth disease virus internal ribosome entry site (IRES); this junction contains highly conserved motifs for RNA-RNA and RNA-protein interactions, important for IRES activity. Our simulations capture transitions of the 4H junction between parallel and antiparallel conformations. The interconversion is virtually barrier-free and occurs via a rotation between the axes of coaxially stacked helices with a transient perpendicular intermediate. We characterize this transition, with various interhelical orientations, by pseudodihedral angle and interhelical distance measures. The high flexibility of the junction, as also demonstrated experimentally, is suitable for IRES activity. Because foot-and-mouth disease virus IRES structure depends on long-range interactions involving domain 3, the perpendicular intermediate, which maintains coaxial stacking of helices and thereby consensus primary and secondary structure information, may be beneficial for guiding the overall organization of the RNA system in domain 3.

Disparate HDV ribozyme crystal structures represent intermediates on a rugged free-energy landscape

The hepatitis delta virus (HDV) ribozyme is a member of the class of small, self-cleaving catalytic RNAs found in a wide range of genomes from HDV to human. Both pre- and post-catalysis (precursor and product) crystal structures of the cis-acting genomic HDV ribozyme have been determined. These structures, together with extensive solution probing, have suggested that a significant conformational change accompanies catalysis. A recent crystal structure of a trans-acting precursor, obtained at low pH and by molecular replacement from the previous product conformation, conforms to the product, raising the possibility that it represents an activated conformer past the conformational change. Here, using fluorescence resonance energy transfer (FRET), we discovered that cleavage of this ribozyme at physiological pH is accompanied by a structural lengthening in magnitude comparable to previous trans-acting HDV ribozymes. Conformational heterogeneity observed by FRET in solution appears to have been removed upon crystallization. Analysis of a total of 1.8 µsec of molecular dynamics (MD) simulations showed that the crystallographically unresolved cleavage site conformation is likely correctly modeled after the hammerhead ribozyme, but that crystal contacts and the removal of several 2'-oxygens near the scissile phosphate compromise catalytic in-line fitness. A cis-acting version of the ribozyme exhibits a more dynamic active site, while a G-1 residue upstream of the scissile phosphate favors poor fitness, allowing us to rationalize corresponding changes in catalytic activity. Based on these data, we propose that the available crystal structures of the HDV ribozyme represent intermediates on an overall rugged RNA folding free-energy landscape.

Dynamics of riboswitches: Molecular simulations

Riboswitch RNAs play key roles in bacterial metabolism and represent a promising new class of antibiotic targets for treatment of infectious disease. While many studies of riboswitches have been performed, the exact mechanism of riboswitch operation is still not fully understood at the atomistic level of detail. Molecular dynamics simulations are useful for interpreting existing experimental data and producing predictions for new experiments. Here, a wide range of computational studies on riboswitches is reviewed. By elucidating the key principles of riboswitch operation, computation may aid in the effort to design more specific antibiotics with affinities greater than those of the native ligand. Such a detailed understanding may be required to improve efficacy and reduce side effects. These studies are laying the groundwork for understanding the action mechanism of new compounds that inhibit riboswitch activity. Future directions such as magnesium effects, large-scale conformational changes, expression platforms and co-transcriptional folding are also discussed. This article is part of a Special Issue entitled: Riboswitches.

Cooperative gene regulation by microRNA pairs and their identification using a computational workflow

MicroRNAs (miRNAs) are an integral part of gene regulation at the post-transcriptional level. Recently, it has been shown that pairs of miRNAs can repress the translation of a target mRNA in a cooperative manner, which leads to an enhanced effectiveness and specificity in target repression. However, it remains unclear which miRNA pairs can synergize and which genes are target of cooperative miRNA regulation. In this paper, we present a computational workflow for the prediction and analysis of cooperating miRNAs and their mutual target genes, which we refer to as RNA triplexes. The workflow integrates methods of miRNA target prediction; triplex structure analysis; molecular dynamics simulations and mathematical modeling for a reliable prediction of functional RNA triplexes and target repression efficiency. In a case study we analyzed the human genome and identified several thousand targets of cooperative gene regulation. Our results suggest that miRNA cooperativity is a frequent mechanism for an enhanced target repression by pairs of miRNAs facilitating distinctive and fine-tuned target gene expression patterns. Human RNA triplexes predicted and characterized in this study are organized in a web resource at www.sbi.uni-rostock.de/triplexrna/.

Molecular Dynamics Simulations of Nucleic Acids. From Tetranucleotides to the Ribosome

We present a brief overview of explicit solvent molecular dynamics (MD) simulations of nucleic acids. We explain physical chemistry limitations of the simulations, namely, the molecular mechanics (MM) force field (FF) approximation and limited time scale. Further, we discuss relations and differences between simulations and experiments, compare standard and enhanced sampling simulations, discuss the role of starting structures, comment on different versions of nucleic acid FFs, and relate MM computations with contemporary quantum chemistry. Despite its limitations, we show that MD is a powerful technique for studying the structural dynamics of nucleic acids with a fast growing potential that substantially complements experimental results and aids their interpretation.

A novel implicit solvent model for simulating the molecular dynamics of RNA

Although molecular dynamics simulations can be accelerated by more than an order of magnitude by implicitly describing the influence of the solvent with a continuum model, most currently available implicit solvent simulations cannot robustly simulate the structure and dynamics of nucleic acids. The difficulties become exacerbated especially for RNAs, suggesting the presence of serious physical flaws in the prior continuum models for the influence of the solvent and counter ions on the nucleic acids. We present a novel, to our knowledge, implicit solvent model for simulating nucleic acids by combining the Langevin-Debye model and the Poisson-Boltzmann equation to provide a better estimate of the electrostatic screening of both the water and counter ions. Tests of the model involve comparisons of implicit and explicit solvent simulations for three RNA targets with 20, 29, and 75 nucleotides. The model provides reasonable agreement with explicit solvent simulations, and directions for future improvement are noted.

Multidimensional Replica Exchange Molecular Dynamics Yields a Converged Ensemble of an RNA Tetranucleotide

A necessary step to properly assess and validate the performance of force fields for biomolecules is to exhaustively sample the accessible conformational space, which is challenging for large RNA structures. Given questions regarding the reliability of modeling RNA structure and dynamics with current methods, we have begun to use RNA tetranucleotides to evaluate force fields. These systems, though small, display considerable conformational variability and complete sampling with standard simulation methods remains challenging. Here we compare and discuss the performance of known variations of replica exchange molecular dynamics (REMD) methods, specifically temperature REMD (T-REMD), Hamiltonian REMD (H-REMD), and multidimensional REMD (M-REMD) methods, which have been implemented in Amber’s accelerated GPU code. Using two independent simulations, we show that M-REMD not only makes very efficient use of emerging large-scale GPU clusters, like Blue Waters at the University of Illinois, but also is critically important in generating the converged ensemble more efficiently than either T-REMD or H-REMD. With 57.6 μs aggregate sampling of a conformational ensemble with M-REMD methods, the populations can be compared to NMR data to evaluate force field reliability and further understand how putative changes to the force field may alter populations to be in more consistent agreement with experiment.

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.

Novel features for identifying A-minors in three-dimensional RNA molecules

RNA tertiary interactions or tertiary motifs are conserved structural patterns formed by pairwise interactions between nucleotides. They include base-pairing, base-stacking, and base-phosphate interactions. A-minor motifs are the most common tertiary interactions in the large ribosomal subunit. The A-minor motif is a nucleotide triple in which minor groove edges of an adenine base are inserted into the minor groove of neighboring helices, leading to interaction with a stabilizing base pair. We propose here novel features for identifying and predicting A-minor motifs in a given three-dimensional RNA molecule. By utilizing the features together with machine learning algorithms including random forests and support vector machines, we show experimentally that our approach is capable of predicting A-minor motifs in the given RNA molecule effectively, demonstrating the usefulness of the proposed approach. The techniques developed from this work will be useful for molecular biologists and biochemists to analyze RNA tertiary motifs, specifically A-minor interactions.

Effects of hypoxanthine substitution in peptide nucleic acids targeting KRAS2 oncogenic mRNA molecules: theory and experiment

Genetic disorders can arise from single base substitutions in a single gene. A single base substitution for wild type guanine in the twelfth codon of KRAS2 mRNA occurs frequently to initiate lung, pancreatic, and colon cancer. We have observed single base mismatch specificity in radioimaging of mutant KRAS2 mRNA in tumors in mice by in vivo hybridization with radiolabeled peptide nucleic acid (PNA) dodecamers. We hypothesized that multimutant specificity could be achieved with a PNA dodecamer incorporating hypoxanthine, which can form Watson-Crick base pairs with adenine, cytosine, thymine, and uracil. Using molecular dynamics simulations and free energy calculations, we show that hypoxanthine substitutions in PNAs are tolerated in KRAS2 RNA:PNA duplexes where wild type guanine is replaced by mutant uracil or adenine in RNA. To validate our predictions, we synthesized PNA dodecamers with hypoxanthine, and then measured the thermal stability of RNA:PNA duplexes. Circular dichroism thermal melting results showed that hypoxanthine-containing PNAs are more stable in duplexes where hypoxanthine-adenine and hypoxanthine-uracil base pairs are formed than single mismatch duplexes or duplexes containing hypoxanthine-guanine opposition.

Biomolecular simulation: a computational microscope for molecular biology

Molecular dynamics simulations capture the behavior of biological macromolecules in full atomic detail, but their computational demands, combined with the challenge of appropriately modeling the relevant physics, have historically restricted their length and accuracy. Dramatic recent improvements in achievable simulation speed and the underlying physical models have enabled atomic-level simulations on timescales as long as milliseconds that capture key biochemical processes such as protein folding, drug binding, membrane transport, and the conformational changes critical to protein function. Such simulation may serve as a computational microscope, revealing biomolecular mechanisms at spatial and temporal scales that are difficult to observe experimentally. We describe the rapidly evolving state of the art for atomic-level biomolecular simulation, illustrate the types of biological discoveries that can now be made through simulation, and discuss challenges motivating continued innovation in this field.

RNA/Peptide Binding Driven by Electrostatics-Insight from Bidirectional Pulling Simulations

RNA/protein interactions play crucial roles in controlling gene expression. They are becoming important targets for pharmaceutical applications. Due to RNA flexibility and to the strength of electrostatic interactions, standard docking methods are insufficient. We here present a computational method which allows studying the binding of RNA molecules and charged peptides with atomistic, explicit-solvent molecular dynamics. In our method, a suitable estimate of the electrostatic interaction is used as an order parameter (collective variable) which is then accelerated using bidirectional pulling simulations. Since the electrostatic interaction is only used to enhance the sampling, the approximations used to compute it do not affect the final accuracy. The method is employed to characterize the binding of TAR RNA from HIV-1 and a small cyclic peptide. Our simulation protocol allows blindly predicting the binding pocket and pose as well as the binding affinity. The method is general and could be applied to study other electrostatics-driven binding events.

Bridging the gap between theory and experiment to derive a detailed understanding of hammerhead ribozyme catalysis

Herein we summarize our progress toward the understanding of hammerhead ribozyme (HHR) catalysis through a multiscale simulation strategy. Simulation results collectively paint a picture of HHR catalysis: HHR first folds to form an electronegative active site pocket to recruit a threshold occupation of cationic charges, either a Mg(2+) ion or multiple monovalent cations. Catalytically active conformations that have good in-line fitness are supported by specific metal ion coordination patterns that involve either a bridging Mg(2+) ion or multiple Na(+) ions, one of which is also in a bridging coordination pattern. In the case of a single Mg(2+) ion bound in the active site, the Mg(2+) ion undergoes a migration that is coupled with deprotonation of the nucleophile (C17:O2'). As the reaction proceeds, the Mg(2+) ion stabilizes the accumulating charge of the leaving group and significantly increases the general acid ability of G8:O2'. Further computational mutagenesis simulations suggest that the disruptions due to mutations may severely impact HHR catalysis at different stages of the reaction. Catalytic mechanisms supported by the simulation results are consistent with available structural and biochemical experiments, and together they advance our understanding of HHR catalysis.

Candidate RNA structures for domain 3 of the foot-and-mouth-disease virus internal ribosome entry site

The foot-and-mouth-disease virus (FMDV) utilizes non-canonical translation initiation for viral protein synthesis, by forming a specific RNA structure called internal ribosome entry site (IRES). Domain 3 in FMDV IRES is phylogenetically conserved and highly structured; it contains four-way junctions where intramolecular RNA–RNA interactions serve as a scaffold for the RNA to fold for efficient IRES activity. Although the 3D structure of domain 3 is crucial to exploring and deciphering the initiation mechanism of translation, little is known. Here, we employ a combination of various modeling approaches to propose candidate tertiary structures for the apical region of domain 3, thought to be crucial for IRES function. We begin by modeling junction topology candidates and build atomic 3D models consistent with available experimental data. We then investigate each of the four candidate 3D structures by molecular dynamics simulations to determine the most energetically favorable configurations and to analyze specific tertiary interactions. Only one model emerges as viable containing not only the specific binding site for the GNRA tetraloop but also helical arrangements which enhance the stability of domain 3. These collective findings, together with available experimental data, suggest a plausible theoretical tertiary structure of the apical region in FMDV IRES domain 3.

Influence of thiostrepton binding on the ribosomal GTPase associated region characterized by molecular dynamics simulation

The thiostrepton antibiotic inhibits bacterial protein synthesis by binding to a cleft formed by the ribosomal protein L11 and 23S's rRNA helices 43-44 on the 70S ribosome. It was proposed from crystal structures that the ligand restricts L11's N-terminal movement and thus prevents proper translation factor binding. An exact understanding of thiostrepton's impact on the binding site's dynamics at atomistic resolution is still missing. Here we report an all-atom molecular dynamics simulations of the binary L11·rRNA and the ternary L11·rRNA·thiostrepton complex (rRNA = helices 43-44). We demonstrate that thiostrepton directly impacts the binding site's atomic and biomacromolecular dynamics.

Molecular dynamics re-refinement of two different small RNA loop structures using the original NMR data suggest a common structure

Restrained molecular dynamics simulations are a robust, though perhaps underused, tool for the end-stage refinement of biomolecular structures. We demonstrate their utility-using modern simulation protocols, optimized force fields, and inclusion of explicit solvent and mobile counterions-by re-investigating the solution structures of two RNA hairpins that had previously been refined using conventional techniques. The structures, both domain 5 group II intron ribozymes from yeast ai5γ and Pylaiella littoralis, share a nearly identical primary sequence yet the published 3D structures appear quite different. Relatively long restrained MD simulations using the original NMR restraint data identified the presence of a small set of violated distance restraints in one structure and a possibly incorrect trapped bulge nucleotide conformation in the other structure. The removal of problematic distance restraints and the addition of a heating step yielded representative ensembles with very similar 3D structures and much lower pairwise RMSD values. Analysis of ion density during the restrained simulations helped to explain chemical shift perturbation data published previously. These results suggest that restrained MD simulations, with proper caution, can be used to "update" older structures or aid in the refinement of new structures that lack sufficient experimental data to produce a high quality result. Notable cautions include the need for sufficient sampling, awareness of potential force field bias (such as small angle deviations with the current AMBER force fields), and a proper balance between the various restraint weights.

Counterion Redistribution upon Binding of a Tat-Protein Mimic to HIV-1 TAR RNA

Binding of proteins and small molecules to RNA involves many electrostatic interactions, which may alter the distribution of ions around the RNA molecule. Here, we use molecular dynamics simulations to investigate how binding of a cyclic peptide mimic of the HIV-1 Tat protein affects the ionic distribution around the HIV-1 TAR RNA element. The calculations reproduce the structural properties observed in NMR studies of TAR and its complex. They also provide insight into the rearrangement of counterions during the molecular recognition events leading to the formation of the protein/RNA complex.

Revision of AMBER Torsional Parameters for RNA Improves Free Energy Predictions for Tetramer Duplexes with GC and iGiC Base Pairs

All-atom force fields are important for predicting thermodynamic, structural, and dynamic properties of RNA. In this paper, results are reported for thermodynamic integration calculations of free energy differences of duplex formation when CG pairs in the RNA duplexes r(CCGG)(2), r(GGCC)(2), r(GCGC)(2), and r(CGCG)(2) are replaced by isocytidine-isoguanosine (iCiG) pairs. Agreement with experiment was improved when ε/ζ, α/γ, β, and χ torsional parameters in the AMBER99 force field were revised on the basis of quantum mechanical calculations. The revised force field, AMBER99TOR, brings free energy difference predictions to within 1.3, 1.4, 2.3, and 2.6 kcal/mol at 300 K, respectively, compared to experimental results for the thermodynamic cycles of CCGG → iCiCiGiG, GGCC → iGiGiCiC, GCGC → iGiCiGiC, and CGCG → iCiGiCiG. In contrast, unmodified AMBER99 predictions for GGCC → iGiGiCiC and GCGC → iGiCiGiC differ from experiment by 11.7 and 12.6 kcal/mol, respectively. In order to test the dynamic stability of the above duplexes with AMBER99TOR, four individual 50 ns molecular dynamics (MD) simulations in explicit solvent were run. All except r(CCGG)(2) retained A-form conformation for ≥82% of the time. This is consistent with NMR spectra of r(iGiGiCiC)(2), which reveal an A-form conformation. In MD simulations, r(CCGG)(2) retained A-form conformation 52% of the time, suggesting that its terminal base pairs may fray. The results indicate that revised backbone parameters improve predictions of RNA properties and that comparisons to measured sequence dependent thermodynamics provide useful benchmarks for testing force fields and computational methods.