Search

In computational physics, chemistry, and biology, the implementation of new techniques in a shared and open source software lowers barriers to entry and promotes rapid scientific progress. However, effectively training new software users presents several challenges. Common methods like direct knowledge transfer and in-person workshops are limited in reach and comprehensiveness. Furthermore, while the COVID-19 pandemic highlighted the benefits of online training, traditional online tutorials can quickly become outdated and may not cover all the software's functionalities. To address these issues, here we introduce ``PLUMED Tutorials'', a collaborative model for developing, sharing, and updating online tutorials. This initiative utilizes repository management and continuous integration to ensure compatibility with software updates. Moreover, the tutorials are interconnected to form a structured learning path and are enriched with automatic annotations to provide broader context. This paper illustrates the development, features, and advantages of PLUMED Tutorials, aiming to foster an open community for creating and sharing educational resources., Comment: 26 pages, 5 figures

Molecular dynamics (MD) simulations play a crucial role in resolving the underlying conformational dynamics of molecular systems. However, their capability to correctly reproduce and predict dynamics in agreement with experiments is limited by the accuracy of the force-field model. This capability can be improved by refining the structural ensembles or the force-field parameters. Furthermore, discrepancies with experimental data can be due to imprecise forward models, namely, functions mapping simulated structures to experimental observables. Here, we introduce MDRefine, a Python package aimed at implementing the refinement of the ensemble, the force-field and/or the forward model by comparing MD-generated trajectories with experimental data. The software consists of several tools that can be employed separately from each other or combined together in different ways, providing a seamless interpolation between these three different types of refinement. We use some benchmark cases to show that the combined approach is superior to separately applied refinements. Source code, documentation and examples are freely available at https://pypi.org/project/MDRefine and https://github.com/bussilab/MDRefine., Comment: 39 pages, 20 figures

RNA molecules exhibit various biological functions intrinsically dependent on their diverse ecosystem of highly flexible structures. This flexibility arises from complex hydrogen-bonding networks defined by canonical and non-canonical base pairs that require protonation events to stabilize or perturb these interactions. Constant pH molecular dynamics (CpHMD) methods provide a reliable framework to explore the conformational and protonation space of dynamic structures and for robust calculations of pH-dependent properties, such as the pK$_\mathrm{a}$ of titrable sites. Despite growing biological evidence concerning pH regulation of certain motifs and in biotechnological applications, pH-sensitive in silico methods have rarely been applied to nucleic acids. In this work, we extended the stochastic titration CpHMD method to include RNA parameters from the standard $\chi$OL3 AMBER force field and highlighted its capability to depict titration events of nucleotides in single-stranded RNAs. We validated the method using trimers and pentamers with a single central titrable site while integrating a well-tempered metadynamics approach into the st-CpHMD methodology (CpH-MetaD) using PLUMED. This approach enhanced the convergence of the conformational landscape and enabled more efficient sampling of protonation-conformation coupling. Our pK$_\mathrm{a}$ estimates agree with experimental data, validating the method's ability to reproduce electrostatic changes around a titrable nucleobase in single-stranded RNA. These findings provided molecular insight into intramolecular phenomena, such as nucleobase stacking and phosphate interactions, that dictate the experimentally observed pK$_\mathrm{a}$ shifts between different strands. Overall, this work validates both the st-CpHMD and the metadynamics integration as reliable tools for studying biologically relevant RNA systems.

This letter illustrates the opinion of the molecular dynamics (MD) community on the need to adopt a new FAIR paradigm for the use of molecular simulations. It highlights the necessity of a collaborative effort to create, establish, and sustain a database that allows findability, accessibility, interoperability, and reusability of molecular dynamics simulation data. Such a development would democratize the field and significantly improve the impact of MD simulations on life science research. This will transform our working paradigm, pushing the field to a new frontier. We invite you to support our initiative at the MDDB community (https://mddbr.eu/community/)

N6-methyladenosine (m6A) is a prevalent RNA post-transcriptional modification that plays crucial roles in RNA stability, structural dynamics, and interactions with proteins. The YT521-B (YTH) family of proteins, which are notable m6A readers, function through their highly conserved YTH domain. Recent structural investigations and molecular dynamics (MD) simulations have shed light on the recognition mechanism of m6A by the YTHDC1 protein. Despite advancements, using MD to predict the stabilization induced by m6A on the free energy of binding between RNA and YTH proteins remains challenging, due to inaccuracy of the employed force field and limited sampling. For instance, simulations often fail to sufficiently capture the hydration dynamics of the binding pocket. This study addresses these challenges through an innovative methodology that integrates metadynamics, alchemical simulations, and force-field refinement. Importantly, our research identifies hydration of the binding pocket as giving only a minor contribution to the binding free energy and emphasizes the critical importance of precisely tuning force-field parameters to experimental data. By employing a fitting strategy built on alchemical calculations, we refine the m6A partial charges parameters, thereby enabling the simultaneous reproduction of N6 methylation on both the protein binding free energy and the thermodynamic stability of nine RNA duplexes. Our findings underscore the sensitivity of binding free energies to partial charges, highlighting the necessity for thorough parameterization and validation against experimental observations across a range of structural contexts.

Ribonucleic acids (RNA) are unique in that they can store genetic information, replicate and perform catalysis. Importantly, RNA molecules are highly dynamic, and thus determining the ensemble of conformations that they populate is crucial not only to elucidate their biological functions, but also for their potential use as therapeutic targets. Computational and experimental techniques provide complementary views on RNA dynamics, and their integration is fundamental to improve the accuracy of computations and the resolution of experiments. Recent exciting developments in this field, were discussed at the CECAM workshop ``RNA dynamics from experimental and computational approaches'', in Paris, June 26-28, 2023. This report outlines key `take-home' messages that emerged during this workshop from the presentations and discussions.

A novel method combining the ensemble refinement by maximum entropy principle and the force field fitting approach is presented. Its formulation allows to continuously interpolate in between these two methods, which can thus be interpreted as two limiting cases. A cross-validation procedure enables to correctly assess the relative weight of both of them, distinguishing scenarios where the combined approach is meaningful from those in which either ensemble refinement or force field fitting separately prevails. The efficacy of their combination is examined for a realistic case study of RNA oligomers. Within the new scheme, molecular dynamics simulations are integrated with experimental data provided by nuclear-magnetic-resonance measures. We show that force field corrections are in general superior when applied to the appropriate force field terms, but are automatically discarded by the method when applied to inappropriate force field terms.

A novel method combining maximum entropy principle, the Bayesian-inference of ensembles approach, and the optimization of empirical forward models is presented. Here we focus on the Karplus parameters for RNA systems, which relate the dihedral angles of $\gamma$, $\beta$, and the dihedrals in the sugar ring to the corresponding $^3J$-coupling signal between coupling protons. Extensive molecular simulations are performed on a set of RNA tetramers and hexamers and combined with available nucleic-magnetic-resonance data. Within the new framework, the sampled structural dynamics can be reweighted to match experimental data while the error arising from inaccuracies in the forward models can be corrected simultaneously and consequently does not leak into the reweighted ensemble. Carefully crafted cross-validation procedure and regularization terms enable obtaining transferable Karplus parameters. Our approach identifies the optimal regularization strength and new sets of Karplus parameters balancing good agreement between simulations and experiments with minimal changes to the original ensemble., Comment: Submitted to journal; added zenodo link; replaced fig. 3 with correct one

Post-transcriptional modifications are crucial for RNA function, with roles ranging from the stabilization of functional RNA structures to modulation of RNA--protein interactions. Additionally, artificially modified RNAs have been suggested as optimal oligonucleotides for therapeutic purposes. The impact of chemical modifications on secondary structure has been rationalized for some of the most common modifications. However, the characterization of how the modifications affect the three-dimensional RNA structure and dynamics and its capability to bind proteins is still highly challenging. Molecular dynamics simulations, coupled with enhanced sampling methods and integration of experimental data, provide a direct access to RNA structural dynamics. In the context of RNA chemical modifications, alchemical simulations where a wild type nucleotide is converted to a modified one are particularly common. In this Chapter, we review recent molecular dynamics studies of modified ribonucleotides. We discuss the technical aspects of the reviewed works, including the employed force fields, enhanced sampling methods, and alchemical methods, in a way that is accessible to experimentalists. Finally, we provide our perspective on this quickly growing field of research. The goal of this Chapter is to provide a guide for experimentalists to understand molecular dynamics works and, at the same time, give to molecular dynamics experts a solid review of published articles that will be a useful starting point for new research., Comment: Submitted as a chapter for the book "RNA Structure and Function", series "RNA Technologies", published by Springer

Chemical probing experiments such as SHAPE are routinely used to probe RNA molecules. In this work, we use atomistic molecular dynamics simulations to test the hypothesis that binding of RNA with SHAPE reagents is affected by cooperative effects leading to an observed reactivity that is dependent on the reagent concentration. We develop a general technique that enables the calculation of the affinity for arbitrary molecules as a function of their concentration in the grand-canonical ensemble. Our simulations of an RNA structural motif suggest that, at the concentration typically used in SHAPE experiments, cooperative binding would lead to a measurable concentration-dependent reactivity. We also provide a qualitative validation of this statement by analyzing a new set of experiments collected at different reagent concentrations.

Conformational dynamics is crucial for ribonucleic acid (RNA) function. Techniques such as nuclear magnetic resonance, cryo-electron microscopy, small- and wide-angle X-ray scattering, chemical probing, single-molecule F\"orster resonance energy transfer or even thermal or mechanical denaturation experiments probe RNA dynamics at different time and space resolutions. Their combination with accurate atomistic molecular dynamics (MD) simulations paves the way for quantitative and detailed studies of RNA dynamics. First, experiments provide a quantitative validation tool for MD simulations. Second, available data can be used to refine simulated structural ensembles to match experiments. Finally, comparison with experiments allows for improving MD force fields that are transferable to new systems for which data is not available. Here we review the recent literature and provide our perspective on this field., Comment: Submitted to journal

Performing alchemical transformations, in which one molecular system is nonphysically changed to another system, is a popular approach adopted in performing free energy calculations associated with various biophysical processes, such as protein-ligand binding or the transfer of a molecule between environments. While the sampling of alchemical intermediate states in either parallel (e.g. Hamiltonian replica exchange) or serial manner (e.g. expanded ensemble) can bridge the high-probability regions in the configurational space between two end states of interest, alchemical methods can fail in scenarios where the most important slow degrees of freedom in the configurational space are in large part orthogonal to the alchemical variable, or if the system gets trapped in a deep basin extending in both the configurational and alchemical space. To alleviate these issues, we propose to use alchemical variables as an additional dimension in metadynamics, augmenting the ability both to sample collective variables and to enhance sampling in free energy calculations. In this study, we validate our implementation of alchemical metadynamics in PLUMED with test systems and alchemical processes with varying complexities and dimensions of collective variable space, including the interconversion between the torsional metastable states of a toy system and the methylation of a nucleoside both in the isolated form and in a duplex. We show that multi-dimensional alchemical metadynamics can address the challenges mentioned above and further accelerate sampling by introducing configurational collective variables. The method can trivially be combined with other metadynamics-based algorithms implemented in PLUMED. The necessary PLUMED code changes have already been released for general use in PLUMED 2.8.

Post-transcriptional modifications are crucial for RNA function and can affect its structure and dynamics. Force-field based classical molecular dynamics simulations are a fundamental tool to characterize biomolecular dynamics and their application to RNA is flourishing. Here we show that the set of force-field parameters for N$^6$-methyladenosine (m$^6$A) developed for the commonly used AMBER force field does not reproduce duplex denaturation experiments and, specifically, cannot be used to describe both paired and unpaired states. Then we use reweighting techniques to derive new parameters matching available experimental data. The resulting force field can be used to properly describe paired and unpaired m$^6$A in both syn and anti conformation, and thus opens the way to the use of molecular simulations to investigate the effects of N6 methylations on RNA structural dynamics., Comment: Supporting information available as ancillary material. Version to be submitted to journal

The capability of current force fields to reproduce RNA structural dynamics is limited. Several methods have been developed to take advantage of experimental data in order to enforce agreement with experiments. We herein extend an existing framework, which allows arbitrarily chosen force-field correction terms to be fitted by quantification of the discrepancy between observables back-calculated from simulation and corresponding experiments. We apply a robust regularization protocol to avoid overfitting, and additionally introduce and compare a number of different regularization strategies, namely L1-, L2-, Kish Size-, Relative Kish Size- and Relative Entropy-penalties. The training set includes a GACC tetramer as well as more challenging systems, namely gcGAGAgc and gcUUCGgc RNA tetraloops. Specific intramolecular hydrogen bonds in the AMBER RNA force field are corrected with automatically determined parameters that we call gHBfix$_{opt}$. A validation involving a separate simulation of a system present in the training set (gcUUCGgc) and new systems not seen during training (CAAU and UUUU tetramers) displays improvements regarding native population of the tetraloop as well as good agreement with NMR-experiments for tetramers when using the new parameters. Then we simulate folded RNAs (a kink-turn and L1 stalk rRNA) including hydrogen bond types not sufficiently present in the training set. This allows a final modification of the parameter set which is named gHBfix21 and is suggested to be applicable to a wider range of RNA systems., Comment: Supporting information included in ancillary files

Small-angle X-ray scattering (SAXS) experiments are increasingly used to probe RNA structure. A number of \emph{forward models} that relate measured SAXS intensities and structural features, and that are suitable to model either explicit-solvent effects or solute dynamics, have been proposed in the past years. Here we introduce an approach that integrates atomistic molecular dynamics simulations and SAXS experiments to reconstruct RNA structural ensembles while simultaneously accounting for both RNA conformational dynamics and explicit-solvent effects. Our protocol exploits SAXS pure-solute forward models and enhanced sampling methods to sample an heterogenous ensemble of structures, with no information towards the experiments provided on-the-fly. The generated structural ensemble is then reweighted through the maximum entropy principle so as to match reference SAXS experimental data at multiple ionic conditions. Importantly, accurate explicit-solvent forward models are used at this reweighting stage. We apply this framework to the GTPase-associated center, a relevant RNA molecule involved in protein translation, in order to elucidate its ion-dependent conformational ensembles. We show that (a) both solvent and dynamics are crucial to reproduce experimental SAXS data and (b) the resulting dynamical ensembles contain an ion-dependent fraction of extended structures., Comment: Supporting information included in ancillary files. This version includes corrections implemented after receiving feedbacks from the community

Molecular dynamics simulations require barostats to be performed at constant pressure. The usual recipe is to employ the Berendsen barostat first, which displays a first-order volume relaxation efficient in equilibration but results in incorrect volume fluctuations, followed by a second order or Monte Carlo barostat for production runs. In this paper, we introduce stochastic cell rescaling, a first-order barostat that samples the correct volume fluctuations by including a suitable noise term. The algorithm is shown to report volume fluctuations compatible with the isobaric ensemble and its anisotropic variant is tested on a membrane simulation. Stochastic cell rescaling can be straightforwardly implemented in existing codes and can be used effectively both in equilibration and in production phases., Comment: Supporting information included in ancillary files. The article has been accepted by J. Chem. Phys

Biomolecular force fields have been traditionally derived based on a mixture of reference quantum chemistry data and experimental information obtained on small fragments. However, the possibility to run extensive molecular dynamics simulations on larger systems achieving ergodic sampling is paving the way to directly using such simulations along with solution experiments obtained on macromolecular systems. Recently, a number of methods have been introduced to automatize this approach. Here we review these methods, highlight their relationship with machine learning methods, and discuss the open challenges in the field., Comment: Supporting information included in ancillary files. The article has been accepted by J. Chem. Phys

RNA function crucially depends on its structure. Thermodynamic models currently used for secondary structure prediction rely on computing the partition function of folding ensembles, and can thus estimate minimum free-energy structures and ensemble populations. These models sometimes fail in identifying native structures unless complemented by auxiliary experimental data. Here, we build a set of models that combine thermodynamic parameters, chemical probing data (DMS, SHAPE), and co-evolutionary data (Direct Coupling Analysis, DCA) through a network that outputs perturbations to the ensemble free energy. Perturbations are trained to increase the ensemble populations of a representative set of known native RNA structures. In the chemical probing nodes of the network, a convolutional window combines neighboring reactivities, enlightening their structural information content and the contribution of local conformational ensembles. Regularization is used to limit overfitting and improve transferability. The most transferable model is selected through a cross-validation strategy that estimates the performance of models on systems on which they are not trained. With the selected model we obtain increased ensemble populations for native structures and more accurate predictions in an independent validation set. The flexibility of the approach allows the model to be easily retrained and adapted to incorporate arbitrary experimental information., Comment: Supporting Information are reported in a Jupyter notebook, in both format .pdf (read-only) and .ipynb (interactive and editable). Full data required to run the .ipynb notebook can be found at the GitHub repository reported in the notebook

This chapter discusses how the PLUMED plugin for molecular dynamics can be used to analyze and bias molecular dynamics trajectories. The chapter begins by introducing the notion of a collective variable and by then explaining how the free energy can be computed as a function of one or more collective variables. A number of practical issues mostly around periodic boundary conditions that arise when these types of calculations are performed using PLUMED are then discussed. Later parts of the chapter discuss how PLUMED can be used to perform enhanced sampling simulations that introduce simulation biases or multiple replicas of the system and Monte Carlo exchanges between these replicas. This section is then followed by a discussion on how free-energy surfaces and associated error bars can be extracted from such simulations by using weighted histogram and block averaging techniques.

Many non-coding RNAs are known to play a role in the cell directly linked to their structure. Structure prediction based on the sole sequence is however a challenging task. On the other hand, thanks to the low cost of sequencing technologies, a very large number of homologous sequences are becoming available for many RNA families. In the protein community, it has emerged in the last decade the idea of exploiting the covariance of mutations within a family to predict the protein structure using the direct-coupling-analysis (DCA) method. The application of DCA to RNA systems has been limited so far. We here perform an assessment of the DCA method on 17 riboswitch families, comparing it with the commonly used mutual information analysis and with state-of-the-art R-scape covariance method. We also compare different flavors of DCA, including mean-field, pseudo-likelihood, and a proposed stochastic procedure (Boltzmann learning) for solving exactly the DCA inverse problem. Boltzmann learning outperforms the other methods in predicting contacts observed in high resolution crystal structures., Comment: Supporting information included in ancillary files

The process of RNA base fraying (i.e. the transient opening of the termini of a helix) is involved in many aspects of RNA dynamics. We here use molecular dynamics simulations and Markov state models to characterize the kinetics of RNA fraying and its sequence and direction dependence. In particular, we first introduce a method for determining biomolecular dynamics employing core-set Markov state models constructed using an advanced clustering technique. The method is validated on previously reported simulations. We then use the method to analyze extensive trajectories for four different RNA model duplexes. Results obtained using D. E. Shaw research and AMBER force fields are compared and discussed in detail, and show a non-trivial interplay between the stability of intermediate states and the overall fraying kinetics., Comment: Supporting information included in ancillary files

Molecular dynamics (MD) simulations allow investigating the structural dynamics of biomolecular systems with unrivaled time and space resolution. However, in order to compensate for the inaccuracies of the utilized empirical force fields, it is becoming common to integrate MD simulations with experimental data obtained from ensemble measurements. We here review the approaches that can be used to combine MD and experiment under the guidance of the maximum entropy principle. We mostly focus on methods based on Lagrangian multipliers, either implemented as reweighting of existing simulations or through an on-the-fly optimization. We discuss how errors in the experimental data can be modeled and accounted for. Finally, we use simple model systems to illustrate the typical difficulties arising when applying these methods., Comment: Submitted to Computation

We introduce the SPlit-and-conQueR (SPQR) model, a coarse-grained representation of RNA designed for structure prediction and refinement. In our approach, the representation of a nucleotide consists of a point particle for the phosphate group and an anisotropic particle for the nucleoside. The interactions are, in principle, knowledge-based potentials inspired by the ESCORE function, a base-centered scoring function. However, a special treatment is given to base-pairing interactions and certain geometrical conformations which are lost in a raw knowledge-base model. This results in a representation able to describe planar canonical and non-canonical base pairs and base-phosphate interactions and to distinguish sugar puckers and glycosidic torsion conformations. The model is applied to the folding of several structures, including duplexes with internal loops of non-canonical base pairs, tetraloops, junctions and a pseudoknot. For the majority of these systems, experimental structures are correctly predicted at the level of individual contacts. We also propose a method for efficiently reintroducing atomistic detail from the coarse-grained representation., Comment: Accepted for publication on Nucleic Acids Research

Coarse-grained models can be of great help to address the problem of structure prediction in nucleic acids. On one hand they can make the prediction more efficient, while on the other hand, they can also help to identify the essential degrees of freedom and interactions for the description of a number of structures. With the aim to provide an all-atom representation in an explicit solvent to the predictions of our SPlit and conQueR (SPQR) coarse-grained model of RNA, we recently introduced a backmapping procedure which enforces the predicted structure into an atomistic one by means of steered Molecular Dynamics. These simulations minimize the ERMSD, a particular metric which deals exclusively with the relative arrangement of nucleobases, between the atomistic representation and the target structure. In this paper, we explore the effects of this approach on the resulting interaction networks and backbone conformations by applying it on a set of fragments using as a target their native structure. We find that the geometry of the target structures can be reliably recovered, with limitations in the regions with unpaired bases such as bulges. In addition, we observe that the folding pathway can also change depending on the parameters used in the definition of the ERMSD and the use of other metrics such as the RMSD., Comment: Submitted

RNA function is intimately related to its structural dynamics. Molecular dynamics simulations are useful for exploring biomolecular flexibility but are severely limited by the accessible timescale. Enhanced sampling methods allow this timescale to be effectively extended in order to probe biologically-relevant conformational changes and chemical reactions. Here, we review the role of enhanced sampling techniques in the study of RNA systems. We discuss the challenges and promises associated with the application of these methods to force-field validation, exploration of conformational landscapes and ion/ligand-RNA interactions, as well as catalytic pathways. Important technical aspects of these methods, such as the choice of the biased collective variables and the analysis of multi-replica simulations, are examined in detail. Finally, a perspective on the role of these methods in the characterization of RNA dynamics is provided., Comment: Accepted for publication on Current Opinion in Structural Biology

Ribonucleic acid (RNA) is involved in many regulatory and catalytic processes in the cell. The function of any RNA molecule is intimately related with its structure. In-line probing experiments provide valuable structural datasets for a variety of RNAs and are used to characterize conformational changes in riboswitches. However, the structural determinants that lead to differential reactivities in unpaired nucleotides have not been investigated yet. In this work we used a combination of theoretical approaches, i.e., classical molecular dynamics simulations, multiscale quantum mechanical/molecular mechanical calculations, and enhanced sampling techniques in order to compute and interpret the differential reactivity of individual residues in several RNA motifs including members of the most important GNRA and UNCG tetraloop families. Simulations on the multi ns timescale are required to converge the related free-energy landscapes. The results for uGAAAg and cUUCGg tetraloops and double helices are compared with available data from in-line probing experiments and show that the introduced technique is able to distinguish between nucleotides of the uGAAAg tetraloop based on their structural predispositions towards phosphodiester backbone cleavage. For the cUUCGg tetraloop, more advanced ab initio calculations would be required. This study is the first attempt to computationally classify chemical probing experiments and paves the way for an identification of tertiary structures based on the measured reactivity of non-reactive nucleotides.

Interaction with divalent cations is of paramount importance for RNA structural stability and function. We here report a detailed molecular dynamics study of all the possible binding sites for Mg$^{2+}$ on a RNA duplex, including both direct (inner sphere) and indirect (outer sphere) binding. In order to tackle sampling issues, we develop a modified version of bias-exchange metadynamics which allows us to simultaneously compute affinities with previously unreported statistical accuracy. Results correctly reproduce trends observed in crystallographic databases. Based on this, we simulate a carefully chosen set of models that allows us to quantify the effects of competition with monovalent cations, RNA flexibility, and RNA hybridization. Our simulations reproduce the decrease and increase of Mg$^{2+}$ affinity due to ion competition and hybridization respectively, and predict that RNA flexibility has a site dependent effect. This suggests a non trivial interplay between RNA conformational entropy and divalent cation binding.

Nowadays different experimental techniques, such as single molecule or relaxation experiments, can provide dynamic properties of biomolecular systems, but the amount of detail obtainable with these methods is often limited in terms of time or spatial resolution. Here we use state-of-the-art computational techniques, namely atomistic molecular dynamics and Markov state models, to provide insight into the rapid dynamics of short RNA oligonucleotides, in order to elucidate the kinetics of stacking interactions. Analysis of multiple microsecond-long simulations indicates that the main relaxation modes of such molecules can consist of transitions between alternative folded states, rather than between random coils and native structures. After properly removing structures that are artificially stabilized by known inaccuracies of the current RNA AMBER force field, the kinetic properties predicted are consistent with the timescales of previously reported relaxation experiments., Comment: Accepted for publication on J. Chem. Theory Comput

Background: Guanine quadruplexes (GQs) play vital roles in many cellular processes and are of much interest as drug targets. In contrast to the availability of many structural studies, there is still limited knowledge on GQ folding. Scope of review: We review recent molecular dynamics (MD) simulation studies of the folding of GQs, with an emphasis paid to the human telomeric DNA GQ. We explain the basic principles and limitations of all types of MD methods used to study unfolding and folding in a way accessible to non-specialists. We discuss the potential role of G-hairpin, G-triplex and alternative GQ intermediates in the folding process. We argue that, in general, folding of GQs is fundamentally different from funneled folding of small fast-folding proteins, and can be best described by a kinetic partitioning (KP) mechanism. KP is a competition between at least two (but often many) well-separated and structurally different conformational ensembles. Major conclusions: The KP mechanism is the only plausible way to explain experiments reporting long time-scales of GQ folding and the existence of long-lived sub-states. A significant part of the natural partitioning of the free energy landscape of GQs comes from the ability of the GQ-forming sequences to populate a large number of syn-anti patterns in their G-tracts. The extreme complexity of the KP of GQs typically prevents an appropriate description of the folding landscape using just a few order parameters or collective variables. General significance: We reconcile available computational and experimental studies of GQ folding and formulate basic principles characterizing GQ folding landscapes, Comment: \c{opyright} 2016. This manuscript version is made available under the CC-BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-nd/4.0/

Recent computational efforts have shown that the current potential energy models used in molecular dynamics are not accurate enough to describe the conformational ensemble of RNA oligomers and suggest that molecular dynamics should be complemented with experimental data. We here propose a scheme based on the maximum entropy principle to combine simulations with bulk experiments. In the proposed scheme the noise arising from both the measurements and the forward models used to back calculate the experimental observables is explicitly taken into account. The method is tested on RNA nucleosides and is then used to construct chemically consistent corrections to the Amber RNA force field that allow a large set of experimental data on nucleosides and dinucleosides to be correctly reproduced. The transferability of these corrections is assessed against independent data on tetranucleotides and displays a previously unreported agreement with experiments. This procedure can be applied to enforce multiple experimental data on multiple systems in a self-consistent framework thus suggesting a new paradigm for force field refinement., Comment: Reprinted (adapted) with permission from J Chem Theory Comput doi: 10.1021/acs.jctc.6b00944. Copyright (2016) American Chemical Society

The computer-aided folding of biomolecules, particularly RNAs, is one of the most difficult challenges in computational structural biology. RNA tetraloops are fundamental RNA motifs playing key roles in RNA folding and RNA-RNA and RNA-protein interactions. Although state-of-the-art Molecular Dynamics (MD) force fields correctly describe the native state of these tetraloops as a stable free-energy basin on the microsecond time scale, enhanced sampling techniques reveal that the native state is not the global free energy minimum, suggesting yet unidentified significant imbalances in the force fields. Here, we tested our ability to fold the RNA tetraloops in various force fields and simulation settings. We employed three different enhanced sampling techniques, namely, temperature replica exchange MD (T-REMD), replica exchange with solute tempering (REST2), and well-tempered metadynamics (WT-MetaD). We aimed to separate problems caused by limited sampling from those due to force-field inaccuracies. We found that none of the contemporary force fields is able to correctly describe folding of the 5'-GAGA-3' tetraloop over a range of simulation conditions. We thus aimed to identify which terms of the force field are responsible for this poor description of TL folding. We showed that at least two different imbalances contribute to this behavior, namely, overstabilization of base-phosphate and/or sugar-phosphate interactions and underestimated stability of the hydrogen bonding interaction in base pairing. The first artifact stabilizes the unfolded ensemble, while the second one destabilizes the folded state. The former problem might be partially alleviated by reparametrization of the van der Waals parameters of the phosphate oxygens suggested by Case et al., while in order to overcome the latter effect we suggest local potentials to better capture hydrogen bonding interactions., Comment: Reprinted with permission from J. Chem. Theory Comput. 2016, 12 (9), pp 4534-4548. Copyright 2016 American Chemical Society

We report the folding thermodynamics of ccUUCGgg and ccGAGAgg RNA tetraloops using atomistic molecular dynamics simulations. We obtain a previously unreported estimation of the folding free energy using parallel tempering in combination with well-tempered metadynamics. A key ingredient is the use of a recently developed metric distance, eRMSD, as a biased collective variable. We find that the native fold of both tetraloops is not the global free energy minimum using the Amber\c{hi}OL3 force field. The estimated folding free energies are 30.2kJ/mol for UUCG and 7.5 kJ/mol for GAGA, in striking disagreement with experimental data. We evaluate the viability of all possible one-dimensional backbone force field corrections. We find that disfavoring the gauche+ region of {\alpha} and {\zeta} angles consistently improves the existing force field. The level of accuracy achieved with these corrections, however, cannot be considered sufficient by judging on the basis of available thermodynamic data and solution experiments., Comment: Journal of Physical Chemistry Letters (2016)

The computational study of conformational transitions in nucleic acids still faces many challenges. For example, in the case of single stranded RNA tetranucleotides, agreement between simulations and experiments is not satisfactory due to inaccuracies in the force fields commonly used in molecular dynamics simulations. We here use experimental data collected from high-resolution X-ray structures to attempt an improvement of the latest version of the AMBER force field. A modified metadynamics algorithm is used to calculate correcting potentials designed to enforce experimental distributions of backbone torsion angles. Replica-exchange simulations of tetranucleotides including these correcting potentials show significantly better agreement with independent solution experiments for the oligonucleotides containing pyrimidine bases. Although the proposed corrections do not seem to be portable to generic RNA systems, the simulations revealed the importance of the alpha and beta backbone angles on the modulation of the RNA conformational ensemble. The correction protocol presented here suggests a systematic procedure for force-field refinement.

We introduce a method for predicting RNA folding pathways, with an application to the most important RNA tetraloops. The method is based on the idea that ensembles of three-dimensional fragments extracted from high-resolution crystal structures are heterogeneous enough to describe metastable as well as intermediate states. These ensembles are first validated by performing a quantitative comparison against available solution NMR data of a set of RNA tetranucleotides. Notably, the agreement is better with respect to the one obtained by comparing NMR with extensive all-atom molecular dynamics simulations. We then propose a procedure based on diffusion maps and Markov models that makes it possible to obtain reaction pathways and their relative probabilities from fragment ensembles. This approach is applied to study the helix-to-loop folding pathway of all the tetraloops from the GNRA and UNCG families. The results give detailed insights into the folding mechanism that are compatible with available experimental data and clarify the role of intermediate states observed in previous simulation studies. The method is computationally inexpensive and can be used to study arbitrary conformational transitions., Comment: Accepted for publication on Nucleic Acids Research

DNA G-hairpins are potential key structures participating in folding of human telomeric guanine quadruplexes (GQ). We examined their properties by standard MD simulations starting from the folded state and long T-REMD starting from the unfolded state, accumulating ~130 {\mu}s of atomistic simulations. Antiparallel G-hairpins should spontaneously form in all stages of the folding to support lateral and diagonal loops, with sub-{\mu}s scale rearrangements between them. We found no clear predisposition for direct folding into specific GQ topologies with specific syn/anti patterns. Our key prediction stemming from the T-REMD is that an ideal unfolded ensemble of the full GQ sequence populates all 4096 syn/anti combinations of its four G-stretches. The simulations can propose idealized folding pathways but we explain that such few-state pathways may be misleading. In the context of the available experimental data, the simulations strongly suggest that the GQ folding could be best understood by the kinetic partitioning mechanism with a set of deep competing minima on the folding landscape, with only a small fraction of molecules directly folding to the native fold. The landscape should further include nonspecific collapse processes where the molecules move via diffusion and consecutive random rare transitions, which could, e.g., structure the propeller loops., Comment: This article has been accepted for publication in Nucleic Acids Research Published by Oxford University Press

Non structural protein 3 (NS3) helicase from hepatitis C virus is an enzyme that unwinds and translocates along nucleic acids with an ATP-dependent mechanism and has a key role in the replication of the viral RNA. An inchworm-like mechanism for translocation has been proposed based on crystal structures and single molecule experiments. We here perform atomistic molecular dynamics in explicit solvent on the microsecond time scale of the available experimental structures. We also construct and simulate putative intermediates for the translocation process, and we perform non-equilibrium targeted simulations to estimate their relative stability. For each of the simulated structures we carefully characterize the available conformational space, the ligand binding pocket, and the RNA binding cleft. The analysis of the hydrogen bond network and of the non-equilibrium trajectories indicates an ATP-dependent stabilization of one of the protein conformers. Additionally, enthalpy calculations suggest that entropic effects might be crucial for the stabilization of the experimentally observed structures., Comment: This article has been accepted for publication in Nucleic Acids Research Published by Oxford University Press

Elastic network models (ENMs) are valuable and efficient tools for characterizing the collective internal dynamics of proteins based on the knowledge of their native structures. The increasing evidence that the biological functionality of RNAs is often linked to their innate internal motions, poses the question of whether ENM approaches can be successfully extended to this class of biomolecules. This issue is tackled here by considering various families of elastic networks of increasing complexity applied to a representative set of RNAs. The fluctuations predicted by the alternative ENMs are stringently validated by comparison against extensive molecular dynamics simulations and SHAPE experiments. We find that simulations and experimental data are systematically best reproduced by either an all-atom or a three-beads-per-nucleotide representation (sugar-base-phosphate), with the latter arguably providing the best balance of accuracy and computational complexity., Comment: This article has been accepted for publication in Nucleic Acids Research Published by Oxford University Press

Riboswitches are cis-acting regulatory RNA elements prevalently located in the leader sequences of bacterial mRNA. An adenine sensing riboswitch cis-regulates adeninosine deaminase gene (add) in Vibrio vulnificus. The structural mechanism regulating its conformational changes upon ligand binding mostly remains to be elucidated. In this open framework it has been suggested that the ligand stabilizes the interaction of the distal "kissing loop" complex. Using accurate full-atom molecular dynamics with explicit solvent in combination with enhanced sampling techniques and advanced analysis methods it could be possible to provide a more detailed perspective on the formation of these tertiary contacts. In this work, we used umbrella sampling simulations to study the thermodynamics of the kissing loop complex in the presence and in the absence of the cognate ligand. We enforced the breaking/formation of the loop-loop interaction restraining the distance between the two loops. We also assessed the convergence of the results by using two alternative initialization protocols. A structural analysis was performed using a novel approach to analyze base contacts. Our simulations qualitatively indicated that the ligand could stabilize the kissing loop complex. We also compared with previously published simulation studies. Kissing complex stabilization given by the ligand was compatible with available experimental data. However, the dependence of its value on the initialization protocol of the umbrella sampling simulations posed some questions on the quantitative interpretation of the results and called for better converged enhanced sampling simulations., Comment: Accepted for publication on BMC Bioinformatics

RNA is a fundamental class of biomolecules that mediate a large variety of molecular processes within the cell. Computational algorithms can be of great help in the understanding of RNA structure-function relationship. One of the main challenges in this field is the development of structure-prediction algorithms, which aim at the prediction of the three-dimensional (3D) native fold from the sole knowledge of the sequence. In a recent paper, we have introduced a scoring function for RNA structure prediction. Here, we analyze in detail the performance of the method, we underline strengths and shortcomings, and we discuss the results with respect to state-of-the-art techniques. These observations provide a starting point for improving current methodologies, thus paving the way to the advances of more accurate approaches for RNA 3D structure prediction., Comment: Accepted for publication on RNA & Disease

The computational study of conformational transitions in RNA and proteins with atomistic molecular dynamics often requires suitable enhanced sampling techniques. We here introduce a novel method where concurrent metadynamics are integrated in a Hamiltonian replica-exchange scheme. The ladder of replicas is built with different strength of the bias potential exploiting the tunability of well-tempered metadynamics. Using this method, free-energy barriers of individual collective variables are significantly reduced compared with simple force-field scaling. The introduced methodology is flexible and allows adaptive bias potentials to be self-consistently constructed for a large number of simple collective variables, such as distances and dihedral angles. The method is tested on alanine dipeptide and applied to the difficult problem of conformational sampling in a tetranucleotide., Comment: Accepted for publication on J. Chem. Theory Comput

Many recently introduced enhanced sampling techniques are based on biasing coarse descriptors (collective variables) of a molecular system on the fly. Sometimes the calculation of such collective variables is expensive and becomes a bottleneck in molecular dynamics simulations. An algorithm to treat smooth biasing forces within a multiple time step framework is here discussed. The implementation is simple and allows a speed up when expensive collective variables are employed. The gain can be substantial when using massively parallel or GPU-based molecular dynamics software. Moreover, a theoretical framework to assess the sampling accuracy is introduced, which can be used to assess the choice of the integration time step in both single and multiple time step biased simulations., Comment: Accepted for publication on J. Chem. Theory Comput

The intricate network of interactions observed in RNA three-dimensional structures is often described in terms of a multitude of geometrical properties, including helical parameters, base pairing/stacking, hydrogen bonding and backbone conformation. We show that a simple molecular representation consisting in one oriented bead per nucleotide can account for the fundamental structural properties of RNA. In this framework, canonical Watson-Crick, non-Watson-Crick base-pairing and base-stacking interactions can be unambiguously identified within a well-defined interaction shell. We validate this representation by performing two independent, complementary tests. First, we use it to construct a sequence-independent, knowledge-based scoring function for RNA structural prediction, which compares favorably to fully atomistic, state-of-the-art techniques. Second, we define a metric to measure deviation between RNA structures that directly reports on the differences in the base-base interaction network. The effectiveness of this metric is tested with respect to the ability to discriminate between structurally and kinetically distant RNA conformations, performing better compared to standard techniques. Taken together, our results suggest that this minimalist, nucleobase-centric representation captures the main interactions that are relevant for describing RNA structure and dynamics., Comment: Accepted for publication on Nucleic Acids Research