Pressure control using stochastic cell rescaling

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.

Analyzing and biasing simulations with PLUMED

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.

Assessing the accuracy of direct-coupling analysis for RNA contact prediction

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.

Using the Maximum Entropy Principle to Combine Simulations and Solution Experiments

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.

A nucleobase-centered coarse-grained representation for structure prediction of RNA motifs

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.

Effects and limitations of a nucleobase-driven backmapping procedure for nucleic acids using steered Molecular Dynamics

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.

Exploring RNA structure and dynamics through enhanced sampling simulations

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.

Understanding In-line Probing Experiments by Modeling Cleavage of Non-reactive RNA Nucleotides

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.

Unravelling Mg$^{2+}$-RNA binding with atomistic molecular dynamics

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.

Folding of guanine quadruplex molecules -- funnel-like mechanism or kinetic partitioning? An overview from MD simulation studies

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