The mechanism of RNA base fraying: Molecular dynamics simulations analyzed with core-set Markov state models

Simulation-Guided Conformational Space Exploration to Assess Reactive Conformations of a Ribozyme

Self-splicing ribozymes are small ribonucleic acid (RNA) enzymes that catalyze their own cleavage through a transphosphoesterification reaction. While this process is involved in some specific steps of viral RNA replication and splicing, it is also of importance in the context of the (putative) first autocatalytic RNA-based systems that could have preceded the emergence of modern life. The uncatalyzed phosphoester bond formation is thermodynamically very unfavorable, and many experimental studies have focused on understanding the molecular features of catalysis in these ribozymes. However, chemical reaction paths are short-lived and not easily characterized by experimental approaches, so molecular simulation approaches appear as an ideal tool to unveil the molecular details of the reaction. Here, we focus on the model hairpin ribozyme. We show that identifying a relevant initial conformation for reactivity studies, which is frequently overlooked in mixed quantum-classical studies that predominantly concentrate on the chemical reaction itself, can be highly challenging. These challenges stem from limitations in both available experimental structures (which are chemically altered to prevent self-cleavage) and the accuracy of force fields, together with the necessity for comprehensive sampling. We show that molecular dynamics simulations, combined with extensive conformational phase space exploration with Hamiltonian replica-exchange simulations, enable us to characterize the relevant conformational basins of the minimal hairpin ribozyme in the ligated state prior to self-cleavage. We find that what is usually considered a canonical reactive conformation with active site geometries and hydrogen-bond patterns that are optimal for the addition-elimination reaction with general acid/general base catalysis is metastable and only marginally populated. The thermodynamically stable conformation appears to be consistent with the expectations of a mechanism that does not require the direct participation of ribozyme residues in the reaction. While these observations may suffer from forcefield inaccuracies, all investigated forcefields lead to the same conclusions upon proper sampling, contrasting with previous investigations on shorter timescales suggesting that at least one reparametrization of the Amber99 forcefield allowed to stabilize aligned active site conformations. Our study demonstrates that identifying the most pertinent reactant state conformation holds equal importance alongside the accurate determination of the thermodynamics and kinetics of the chemical steps of the reaction.

Atomistic simulations of RNA duplex thermal denaturation: Sequence- and forcefield-dependence

Double-stranded RNA is the end-product of template-based replication, and is also the functional state of some biological RNAs. Similarly to proteins and DNA, they can be denatured by temperature, with important physiological and technological implications. Here, we use an in silico strategy to probe the thermal denaturation of RNA duplexes. Following previous results that were obtained on a few different duplexes, and which nuanced the canonical 2-state picture of nucleic acid denaturation, we here specifically address three different aspects that greatly improve our description of the temperature-induced dsRNA separation. First, we investigate the effect of the spatial distribution of weak and strong base-pairs among the duplex sequence. We show that the deviations from the two-state dehybridization mechanism are more pronounced when a strong core is flanked with weak extremities, while duplexes with a weak core but strong extremities exhibit a two-state behavior, which can be explained by the key role played by base fraying. This was later verified by generating artificial hairpin or circular states containing one or two locked duplex extremities, which results in an important reinforcement of the entire HB structure of the duplex and higher melting temperatures. Finally, we demonstrate that our results are little sensitive to the employed combination of RNA and water forcefields. The trends in thermal stability among the different sequences as well as the observed unfolding mechanisms (and the deviations from a two-state scenario) remain the same regardless of the employed atomistic models. However, our study points to possible limitations of recent reparametrizations of the Amber RNA forcefield, which sometimes results in duplexes that readily denature under ambient conditions, in contradiction with available experimental results.

Mechanism of Ligand Binding to Theophylline RNA Aptamer

Studying RNA-ligand interactions and quantifying their binding thermodynamics and kinetics are of particular relevance in the field of drug discovery. Here, we combined biochemical binding assays and accelerated molecular simulations to investigate ligand binding and dissociation in RNA using the theophylline-binding RNA as a model system. All-atom simulations using a Ligand Gaussian accelerated Molecular Dynamics method (LiGaMD) have captured repetitive binding and dissociation of theophylline and caffeine to RNA. Theophylline's binding free energy and kinetic rate constants align with our experimental data, while caffeine's binding affinity is over 10,000 times weaker, and its kinetics could not be determined. LiGaMD simulations allowed us to identify distinct low-energy conformations and multiple ligand binding pathways to RNA. Simulations revealed a "conformational selection" mechanism for ligand binding to the flexible RNA aptamer, which provides important mechanistic insights into ligand binding to the theophylline-binding model. Our findings suggest that compound docking using a structural ensemble of representative RNA conformations would be necessary for structure-based drug design of flexible RNA.

Toward a Molecular Mechanism of Complementary RNA Duplexes Denaturation

RNA duplexes are relatively rare but play very important biological roles. As an end-product of template-based RNA replication, they also have key implications for hypothetical primitive forms of life. Unless they are specifically separated by enzymes, these duplexes denature upon a temperature increase. However, mechanistic and kinetic aspects of RNA (and DNA) duplex thermal denaturation remain unclear at the microscopic level. We propose an in silico strategy that probes the thermal denaturation of RNA duplexes and allows for an extensive conformational space exploration along a wide temperature range with atomistic precision. We show that this approach first accounts for the strong sequence and length dependence of the duplexes melting temperature, reproducing the trends seen in the experiments and predicted by nearest-neighbor models. The simulations are then instrumental at providing a molecular picture of the temperature-induced strand separation. The textbook canonical "all-or-nothing" two-state model, very much inspired by the protein folding mechanism, can be nuanced. We demonstrate that a temperature increase leads to significantly distorted but stable structures with extensive base-fraying at the extremities, and that the fully formed duplexes typically do not form around melting. The duplex separation therefore appears as much more gradual than commonly thought.

RNA:RNA interaction in ternary complexes resolved by chemical probing

RNA regulation can be performed by a second targeting RNA molecule, such as in the microRNA regulation mechanism. Selective 2'-hydroxyl acylation analyzed by primer extension (SHAPE) probes the structure of RNA molecules and can resolve RNA:protein interactions, but RNA:RNA interactions have not yet been addressed with this technique. Here, we apply SHAPE to investigate RNA-mediated binding processes in RNA:RNA and RNA:RNA-RBP complexes. We use RNA:RNA binding by SHAPE (RABS) to investigate microRNA-34a (miR-34a) binding its mRNA target, the silent information regulator 1 (mSIRT1), both with and without the Argonaute protein, constituting the RNA-induced silencing complex (RISC). We show that the seed of the mRNA target must be bound to the microRNA loaded into RISC to enable further binding of the compensatory region by RISC, while the naked miR-34a is able to bind the compensatory region without seed interaction. The method presented here provides complementary structural evidence for the commonly performed luciferase-assay-based evaluation of microRNA binding-site efficiency and specificity on the mRNA target site and could therefore be used in conjunction with it. The method can be applied to any nucleic acid-mediated RNA- or RBP-binding process, such as splicing, antisense RNA binding, or regulation by RISC, providing important insight into the targeted RNA structure.

Thermodynamics and kinetics of an A-U RNA base pair under force studied by molecular dynamics simulations

Mechanical force has been widely used to study RNA folding and unfolding. Understanding how the force affects the opening and closing of a single base pair, which is a basic step for RNA folding and unfolding and a fundamental behavior in some important biological activities, is crucial to understanding the mechanism of RNA folding and unfolding under mechanical force. In this work, we investigated the opening and closing process of an RNA base pair under mechanical force with constant-force stretching molecular dynamics simulations. It was found that high mechanical force results in overstretching, and the open state is a high-energy state. The enthalpy and entropy change of the base-pair opening-closing transition were obtained and the results at low forces were in good agreement with the nearest-neighbor model. The temperature and force dependence of the opening and closing rates were also obtained. The position of the transition state for the base-pair opening-closing transition under mechanical force was determined. The free energy barrier of opening a base pair without force is the enthalpy increase, and the work done by the force from the closed state to the transition state decreases the barrier and increases the opening rate. The free energy barrier of closing the base pair without force results from the entropy loss, and the work done by the force from the open state to the transition state increases the barrier and decreases the closing rate. The transition rates are strongly dependent on the temperature and force, while the transition path times are weakly dependent on force and temperature.

An overview of structural approaches to study therapeutic RNAs

RNAs provide considerable opportunities as therapeutic agent to expand the plethora of classical therapeutic targets, from extracellular and surface proteins to intracellular nucleic acids and its regulators, in a wide range of diseases. RNA versatility can be exploited to recognize cell types, perform cell therapy, and develop new vaccine classes. Therapeutic RNAs (aptamers, antisense nucleotides, siRNA, miRNA, mRNA and CRISPR-Cas9) can modulate or induce protein expression, inhibit molecular interactions, achieve genome editing as well as exon-skipping. A common RNA thread, which makes it very promising for therapeutic applications, is its structure, flexibility, and binding specificity. Moreover, RNA displays peculiar structural plasticity compared to proteins as well as to DNA. Here we summarize the recent advances and applications of therapeutic RNAs, and the experimental and computational methods to analyze their structure, by biophysical techniques (liquid-state NMR, scattering, reactivity, and computational simulations), with a focus on dynamic and flexibility aspects and to binding analysis. This will provide insights on the currently available RNA therapeutic applications and on the best techniques to evaluate its dynamics and reactivity.

DADApy: Distance-based analysis of data-manifolds in Python

DADApy is a Python software package for analyzing and characterizing high-dimensional data manifolds. It provides methods for estimating the intrinsic dimension and the probability density, for performing density-based clustering, and for comparing different distance metrics. We review the main functionalities of the package and exemplify its usage in a synthetic dataset and in a real-world application. DADApy is freely available under the open-source Apache 2.0 license.

A Computational Study of RNA Tetraloop Thermodynamics, Including Misfolded States

An important characteristic of RNA folding is the adoption of alternative configurations of similar stability, often referred to as misfolded configurations. These configurations are considered to compete with correctly folded configurations, although their rigorous thermodynamic and structural characterization remains elusive. Tetraloop motifs found in large ribozymes are ideal systems for an atomistically detailed computational quantification of folding free energy landscapes and the structural characterization of their constituent free energy basins, including nonnative states. In this work, we studied a group of closely related 10-mer tetraloops using a combined parallel tempering and metadynamics technique that allows a reliable sampling of the free energy landscapes, requiring only knowledge that the stem folds into a canonical A-RNA configuration. We isolated and analyzed unfolded, folded, and misfolded populations that correspond to different free energy basins. We identified a distinct misfolded state that has a stability very close to that of the correctly folded state. This misfolded state contains a predominant population that shares the same structural features across all tetraloops studied here and lacks the noncanonical A-G base pair in its loop portion. Further analysis performed with biased trajectories showed that although this competitive misfolded state is not an essential intermediate, it is visited in most of the transitions from unfolded to correctly folded states. Moreover, the tetraloops can transition from this misfolded state to the correctly folded state without requiring extensive unfolding.

Determining Sequence-Dependent DNA Oligonucleotide Hybridization and Dehybridization Mechanisms Using Coarse-Grained Molecular Simulation, Markov State Models, and Infrared Spectroscopy

A robust understanding of the sequence-dependent thermodynamics of DNA hybridization has enabled rapid advances in DNA nanotechnology. A fundamental understanding of the sequence-dependent kinetics and mechanisms of hybridization and dehybridization remains comparatively underdeveloped. In this work, we establish new understanding of the sequence-dependent hybridization/dehybridization kinetics and mechanism within a family of self-complementary pairs of 10-mer DNA oligomers by integrating coarse-grained molecular simulation, machine learning of the slow dynamical modes, data-driven inference of long-time kinetic models, and experimental temperature-jump infrared spectroscopy. For a repetitive ATATATATAT sequence, we resolve a rugged dynamical landscape comprising multiple metastable states, numerous competing hybridization/dehybridization pathways, and a spectrum of dynamical relaxations. Introduction of a G:C pair at the terminus (GATATATATC) or center (ATATGCATAT) of the sequence reduces the ruggedness of the dynamics landscape by eliminating a number of metastable states and reducing the number of competing dynamical pathways. Only by introducing a G:C pair midway between the terminus and the center to maximally disrupt the repetitive nature of the sequence (ATGATATCAT) do we recover a canonical "all-or-nothing" two-state model of hybridization/dehybridization with no intermediate metastable states. Our results establish new understanding of the dynamical richness of sequence-dependent kinetics and mechanisms of DNA hybridization/dehybridization by furnishing quantitative and predictive kinetic models of the dynamical transition network between metastable states, present a molecular basis with which to understand experimental temperature jump data, and furnish foundational design rules by which to rationally engineer the kinetics and pathways of DNA association and dissociation for DNA nanotechnology applications.

Unsupervised Learning Methods for Molecular Simulation Data

Unsupervised learning is becoming an essential tool to analyze the increasingly large amounts of data produced by atomistic and molecular simulations, in material science, solid state physics, biophysics, and biochemistry. In this Review, we provide a comprehensive overview of the methods of unsupervised learning that have been most commonly used to investigate simulation data and indicate likely directions for further developments in the field. In particular, we discuss feature representation of molecular systems and present state-of-the-art algorithms of dimensionality reduction, density estimation, and clustering, and kinetic models. We divide our discussion into self-contained sections, each discussing a specific method. In each section, we briefly touch upon the mathematical and algorithmic foundations of the method, highlight its strengths and limitations, and describe the specific ways in which it has been used-or can be used-to analyze molecular simulation data.

Chelated Magnesium Logic Gate Regulates Riboswitch Pseudoknot Formation

Magnesium plays a critical role in the structure, dynamics, and function of RNA. The precise microscopic effect of chelated magnesium on RNA structure is yet to be explored. Magnesium is known to act through its diffuse cloud around RNA, through the outer sphere (water-mediated), inner sphere, and often chelated ion-mediated interactions. A mechanism is proposed for the role of experimentally discovered site-specific chelated magnesium ions on the conformational dynamics of SAM-I riboswitch aptamers in bacteria. This mechanism is observed with atomistic simulations performed in a physiological mixed salt environment at a high temperature. The simulations were validated with phosphorothioate interference mapping experiments that help to identify crucial inner-sphere Mg²⁺ sites prescribing an appropriate initial distribution of inner- and outer-sphere magnesium ions to maintain a physiological ion concentration of monovalent and divalent salts. A concerted role of two chelated magnesium ions is newly discovered since the presence of both supports the formation of the pseudoknot. This constitutes a logical AND gate. The absence of any of these magnesium ions instigates the dissociation of long-range pseudoknot interaction exposing the inner core of the RNA. A base triple is the epicenter of the magnesium chelation effect. It allosterically controls RNA pseudoknot by bolstering the direct effect of magnesium chelation in protecting the functional fold of RNA to control ON and OFF transcription switching.

Targeting structural features of viral genomes with a nano-sized supramolecular drug

RNA targeting is an exciting frontier for drug design. Intriguing targets include functional RNA structures in structurally-conserved untranslated regions (UTRs) of many lethal viruses. However, computational docking screens, valuable in protein structure targeting, fail for inherently flexible RNA. Herein we harness MD simulations with Markov state modeling to enable nanosize metallo-supramolecular cylinders to explore the dynamic RNA conformational landscape of HIV-1 TAR untranslated region RNA (representative for many viruses) replicating experimental observations. These cylinders are exciting as they have unprecedented nucleic acid binding and are the first supramolecular helicates shown to have anti-viral activity in cellulo: the approach developed in this study provides additional new insight about how such viral UTR structures might be targeted with the cylinder binding into the heart of an RNA-bulge cavity, how that reduces the conformational flexibility of the RNA and molecular details of the insertion mechanism. The approach and understanding developed represents a new roadmap for design of supramolecular drugs to target RNA structural motifs across biology and nucleic acid nanoscience.

Volume-scaled common nearest neighbor clustering algorithm with free-energy hierarchy

The combination of Markov state modeling (MSM) and molecular dynamics (MD) simulations has been shown in recent years to be a valuable approach to unravel the slow processes of molecular systems with increasing complexity. While the algorithms for intermediate steps in the MSM workflow such as featurization and dimensionality reduction have been specifically adapted to MD datasets, conventional clustering methods are generally applied to the discretization step. This work adds to recent efforts to develop specialized density-based clustering algorithms for the Boltzmann-weighted data from MD simulations. We introduce the volume-scaled common nearest neighbor (vs-CNN) clustering that is an adapted version of the common nearest neighbor (CNN) algorithm. A major advantage of the proposed algorithm is that the introduced density-based criterion directly links to a free-energy notion via Boltzmann inversion. Such a free-energy perspective allows a straightforward hierarchical scheme to identify conformational clusters at different levels of a generally rugged free-energy landscape of complex molecular systems.

Salt dependent mesoscopic model for RNA at multiple strand concentrations

Mesoscopic models can be used for the description of the thermodynamic properties of RNA duplexes. With the use of experimental melting temperatures, its parametrization can provide important insights into its hydrogen bonds and stacking interactions as has been done for high sodium concentrations. However, the RNA parametrization for lower salt concentrations is still missing due to the limited amount of published melting temperature data. While the Peyrard-Bishop (PB) parametrization was found to be largely independent of strand concentrations, it requires that all temperatures are provided at the same strand concentrations. Here we adapted the PB model to handle multiple strand concentrations and in this way we were able to make use of an experimental set of temperatures to model the hydrogen bond and stacking interactions at low and intermediate sodium concentrations. For the parametrizations we make a distinction between terminal and internal base pairs, and the resulting potentials were qualitatively similar as we obtained previously for DNA. The main difference from DNA parameters, was the Morse potentials at low sodium concentrations for terminal r(AU) which is stronger than d(AT), suggesting higher hydrogen bond strength.

Free Energy Landscape and Conformational Kinetics of Hoogsteen Base Pairing in DNA vs. RNA

Genetic information is encoded in the DNA double helix, which, in its physiological milieu, is characterized by the iconical Watson-Crick nucleo-base pairing. Recent NMR relaxation experiments revealed the transient presence of an alternative, Hoogsteen (HG) base pairing pattern in naked DNA duplexes, and estimated its relative stability and lifetime. In contrast with DNA, such structures were not observed in RNA duplexes. Understanding HG base pairing is important because the underlying "breathing" motion between the two conformations can significantly modulate protein binding. However, a detailed mechanistic insight into the transition pathways and kinetics is still missing. We performed enhanced sampling simulation (with combined metadynamics and adaptive force-bias method) and Markov state modeling to obtain accurate free energy, kinetics, and the intermediates in the transition pathway between Watson-Crick and HG base pairs for both naked B-DNA and A-RNA duplexes. The Markov state model constructed from our unbiased MD simulation data revealed previously unknown complex extrahelical intermediates in the seemingly simple process of base flipping in B-DNA. Extending our calculation to A-RNA, for which HG base pairing is not observed experimentally, resulted in relatively unstable, single-hydrogen-bonded, distorted Hoogsteen-like bases. Unlike B-DNA, the transition pathway primarily involved base paired and intrahelical intermediates with transition timescales much longer than that of B-DNA. The seemingly obvious flip-over reaction coordinate (i.e., the glycosidic torsion angle) is unable to resolve the intermediates. Instead, a multidimensional picture involving backbone dihedral angles and distance between hydrogen bond donor and acceptor atoms is required to gain insight into the molecular mechanism.

Dynamical matrix propagator scheme for large-scale proton dynamics simulations

We derive a matrix formalism for the simulation of long range proton dynamics for extended systems and timescales. On the basis of an ab initio molecular dynamics simulation, we construct a Markov chain, which allows us to store the entire proton dynamics in an M × M transition matrix (where M is the number of oxygen atoms). In this article, we start from common topology features of the hydrogen bond network of good proton conductors and utilize them as constituent constraints of our dynamic model. We present a thorough mathematical derivation of our approach and verify its uniqueness and correct asymptotic behavior. We propagate the proton distribution by means of transition matrices, which contain kinetic data from both ultra-short (sub-ps) and intermediate (ps) timescales. This concept allows us to keep the most relevant features from the microscopic level while effectively reaching larger time and length scales. We demonstrate the applicability of the transition matrices for the description of proton conduction trends in proton exchange membrane materials.

On the Stability of the Water-Soluble Chlorophyll-Binding Protein (WSCP) Studied by Molecular Dynamics Simulations

The water-soluble chlorophyll-binding protein (WSCP) is assumed to be not a part of the photosynthetic process. Applying molecular dynamics (MD) simulations, we aimed to obtain insight into the exceptional stability of WSCP. We analyzed dynamical features such as the hydrogen bond network, flexibility, and force distributions. The WSCP structure contains two cysteines at the interfaces of every protein chain, which are in close contact with the cysteines of the other dimer. We tested if a connection of these cysteines between different protein chains influences the dynamical behavior to investigate any influences on the thermal stability. We find that the hydrogen bond network is very stable regardless of the presence or absence of the hypothetical disulfide bridges and/or the chlorophyll units. Furthermore, it is found that the phytyl chains of the chlorophyll units are extremely flexible, much more than what is seen in crystal structures. Nonetheless, they seem to protect a photochemically active site of the chlorophylls over the complete simulation time. Finally, we also find that a cavity in the chlorophyll-surrounding sheath exists, which may allow access for individual small molecules to the core of WSCP.

Asymmetric base-pair opening drives helicase unwinding dynamics

The opening of a Watson-Crick double helix is required for crucial cellular processes, including replication, repair, and transcription. It has long been assumed that RNA or DNA base pairs are broken by the concerted symmetric movement of complementary nucleobases. By analyzing thousands of base-pair opening and closing events from molecular simulations, here, we uncover a systematic stepwise process driven by the asymmetric flipping-out probability of paired nucleobases. We demonstrate experimentally that such asymmetry strongly biases the unwinding efficiency of DNA helicases toward substrates that bear highly dynamic nucleobases, such as pyrimidines, on the displaced strand. Duplex substrates with identical thermodynamic stability are thus shown to be more easily unwound from one side than the other, in a quantifiable and predictable manner. Our results indicate a possible layer of gene regulation coded in the direction-dependent unwindability of the double helix.

Kinetic Mechanism of RNA Helix-Terminal Basepairing-A Kinetic Minima Network Analysis

RNA functions are often kinetically controlled. The folding kinetics of RNAs involves global structural changes and local nucleotide movement, such as base flipping. The most elementary step in RNA folding is the closing and opening of a basepair. By integrating molecular dynamics simulation, master equation, and kinetic Monte Carlo simulation, we investigate the kinetics mechanism of RNA helix-terminal basepairing. The study reveals a six-state folding scheme with three dominant folding pathways of tens, hundreds, and thousands of nanoseconds of folding timescales, respectively. The overall kinetics is rate limited by the detrapping of a misfolded state with the overall folding time of 10^-5 s. Moreover, the analysis examines the different roles of the various driving forces, such as the basepairing and stacking interactions and the ion binding/dissociation effects on structural changes. The results may provide useful insights for developing a basepair opening/closing rate model and further kinetics models of large RNAs.