Comprehensive Assessment of Force-Field Performance in Molecular Dynamics Simulations of DNA/RNA Hybrid Duplexes.

Mixed double helices formed by RNA and DNA strands, commonly referred to as hybrid duplexes or hybrids, are essential in biological processes like transcription and reverse transcription. They are also important for their applications in CRISPR gene editing and nanotechnology. Yet, despite their significance, the hybrid duplexes have been seldom modeled by atomistic molecular dynamics methodology, and there is no benchmark study systematically assessing the force-field performance. Here, we present an extensive benchmark study of polypurine tract (PPT) and Dickerson-Drew dodecamer hybrid duplexes using contemporary and commonly utilized pairwise additive and polarizable nucleic acid force fields. Our findings indicate that none of the available force-field choices accurately reproduces all the characteristic structural details of the hybrid duplexes. The AMBER force fields are unable to populate the C3'-endo (north) pucker of the DNA strand and underestimate inclination. The CHARMM force field accurately describes the C3'-endo pucker and inclination but shows base pair instability. The polarizable force fields struggle with accurately reproducing the helical parameters. Some force-field combinations even demonstrate a discernible conflict between the RNA and DNA parameters. In this work, we offer a candid assessment of the force-field performance for mixed DNA/RNA duplexes. We provide guidance on selecting utilizable force-field combinations and also highlight potential pitfalls and best practices for obtaining optimal performance.

Mechanical Stability and Unfolding Pathways of Parallel Tetrameric G-Quadruplexes Probed by Pulling Simulations.

Guanine quadruplex (GQ) is a noncanonical nucleic acid structure formed by guanine-rich DNA and RNA sequences. Folding of GQs is a complex process, where several aspects remain elusive, despite being important for understanding structure formation and biological functions of GQs. Pulling experiments are a common tool for acquiring insights into the folding landscape of GQs. Herein, we applied a computational pulling strategy─steered molecular dynamics (SMD) simulations─in combination with standard molecular dynamics (MD) simulations to explore the unfolding landscapes of tetrameric parallel GQs. We identified anisotropic properties of elastic conformational changes, unfolding transitions, and GQ mechanical stabilities. Using a special set of structural parameters, we found that the vertical component of pulling force (perpendicular to the average G-quartet plane) plays a significant role in disrupting GQ structures and weakening their mechanical stabilities. We demonstrated that the magnitude of the vertical force component depends on the pulling anchor positions and the number of G-quartets. Typical unfolding transitions for tetrameric parallel GQs involve base unzipping, opening of the G-stem, strand slippage, and rotation to cross-like structures. The unzipping was detected as the first and dominant unfolding event, and it usually started at the 3'-end. Furthermore, results from both SMD and standard MD simulations indicate that partial spiral conformations serve as a transient ensemble during the (un)folding of GQs.

Molecular dynamics simulations reveal the parallel stranded d(GGGA)3GGG DNA quadruplex folds via multiple paths from a coil-like ensemble.

G-quadruplexes (G4s) are non-canonical nucleic acid structures that fold through complex processes. Characterization of the G4 folding landscape may help to elucidate biological roles of G4s but is challenging both experimentally and computationally. Here, we achieved complete folding of a three-quartet parallel DNA G4 with (GGGA)3GGG sequence using all-atom explicit-solvent enhanced-sampling molecular dynamics (MD) simulations. The simulations suggested early formation of guanine stacks in the G-tracts, which behave as semi-rigid blocks in the folding process. The folding continues via the formation of a collapsed compact coil-like ensemble. Structuring of the G4 from the coil then proceeds via various cross-like, hairpin, slip-stranded and two-quartet ensembles and can bypass the G-triplex structure. Folding of the parallel G4 does not appear to involve any salient intermediates and is a multi-pathway process. We also carried out an extended set of simulations of parallel G-hairpins. While parallel G-hairpins are extremely unstable when isolated, they are more stable inside the coil structure. On the methodology side, we show that the AMBER DNA force field predicts the folded G4 to be less stable than the unfolded ensemble, uncovering substantial force-field issues. Overall, we provide unique atomistic insights into the folding landscape of parallel-stranded G4 but also reveal limitations of current state-of-the-art MD techniques.

Simple Adjustment of Intranucleotide Base-Phosphate Interaction in the OL3 AMBER Force Field Improves RNA Simulations.

Molecular dynamics (MD) simulations represent an established tool to study RNA molecules. The outcome of MD studies depends, however, on the quality of the force field (ff). Here we suggest a correction for the widely used AMBER OL3 ff by adding a simple adjustment of the nonbonded parameters. The reparameterization of the Lennard-Jones potential for the -H8···O5'- and -H6···O5'- atom pairs addresses an intranucleotide steric clash occurring in the type 0 base-phosphate interaction (0BPh). The nonbonded fix (NBfix) modification of 0BPh interactions (NBfix0BPh modification) was tuned via a reweighting approach and subsequently tested using an extensive set of standard and enhanced sampling simulations of both unstructured and folded RNA motifs. The modification corrects minor but visible intranucleotide clash for the anti nucleobase conformation. We observed that structural ensembles of small RNA benchmark motifs simulated with the NBfix0BPh modification provide better agreement with experiments. No side effects of the modification were observed in standard simulations of larger structured RNA motifs. We suggest that the combination of OL3 RNA ff and NBfix0BPh modification is a viable option to improve RNA MD simulations.

Complexity of Guanine Quadruplex Unfolding Pathways Revealed by Atomistic Pulling Simulations.

Guanine quadruplexes (GQs) are non-canonical nucleic acid structures involved in many biological processes. GQs formed in single-stranded regions often need to be unwound by cellular machinery, so their mechanochemical properties are important. Here, we performed steered molecular dynamics simulations of human telomeric GQs to study their unfolding. We examined four pulling regimes, including a very slow setup with pulling velocity and force load accessible to high-speed atomic force microscopy. We identified multiple factors affecting the unfolding mechanism, i.e.,: (i) the more the direction of force was perpendicular to the GQ channel axis (determined by GQ topology), the more the base unzipping mechanism happened, (ii) the more parallel the direction of force was, GQ opening and cross-like GQs were more likely to occur, (iii) strand slippage mechanism was possible for GQs with an all-anti pattern in a strand, and (iv) slower pulling velocity led to richer structural dynamics with sampling of more intermediates and partial refolding events. We also identified that a GQ may eventually unfold after a force drop under forces smaller than those that the GQ withstood before the drop. Finally, we found out that different unfolding intermediates could have very similar chain end-to-end distances, which reveals some limitations of structural interpretations of single-molecule spectroscopic data.

Atomistic Picture of Opening-Closing Dynamics of DNA Holliday Junction Obtained by Molecular Simulations.

Holliday junction (HJ) is a noncanonical four-way DNA structure with a prominent role in DNA repair, recombination, and DNA nanotechnology. By rearranging its four arms, HJ can adopt either closed or open state. With enzymes typically recognizing only a single state, acquiring detailed knowledge of the rearrangement process is an important step toward fully understanding the biological function of HJs. Here, we carried out standard all-atom molecular dynamics (MD) simulations of the spontaneous opening-closing transitions, which revealed complex conformational transitions of HJs with an involvement of previously unconsidered "half-closed" intermediates. Detailed free-energy landscapes of the transitions were obtained by sophisticated enhanced sampling simulations. Because the force field overstabilizes the closed conformation of HJs, we developed a system-specific modification which for the first time allows the observation of spontaneous opening-closing HJ transitions in unbiased MD simulations and opens the possibilities for more accurate HJ computational studies of biological processes and nanomaterials.

Sensitivity of the RNA Structure to Ion Conditions as Probed by Molecular Dynamics Simulations of Common Canonical RNA Duplexes.

RNA molecules play a key role in countless biochemical processes. RNA interactions, which are of highly diverse nature, are determined by the fact that RNA is a highly negatively charged polyelectrolyte, which leads to intimate interactions with an ion atmosphere. Although RNA molecules are formally single-stranded, canonical (Watson-Crick) duplexes are key components of folded RNAs. A double-stranded (ds) RNA is also important for the design of RNA-based nanostructures and assemblies. Despite the fact that the description of canonical dsRNA is considered the least problematic part of RNA modeling, the imperfect shape and flexibility of dsRNA can lead to imbalances in the simulations of larger RNAs and RNA-containing assemblies. We present a comprehensive set of molecular dynamics (MD) simulations of four canonical A-RNA duplexes. Our focus was directed toward the characterization of the influence of varying ion concentrations and of the size of the solvation box. We compared several water models and four RNA force fields. The simulations showed that the A-RNA shape was most sensitive to the RNA force field, with some force fields leading to a reduced inclination of the A-RNA duplexes. The ions and water models played a minor role. The effect of the box size was negligible, and even boxes with a small fraction of the bulk solvent outside the RNA hydration sphere were sufficient for the simulation of the dsRNA.

Spontaneous binding of single-stranded RNAs to RRM proteins visualized by unbiased atomistic simulations with a rescaled RNA force field.

Recognition of single-stranded RNA (ssRNA) by RNA recognition motif (RRM) domains is an important class of protein-RNA interactions. Many such complexes were characterized using nuclear magnetic resonance (NMR) and/or X-ray crystallography techniques, revealing ensemble-averaged pictures of the bound states. However, it is becoming widely accepted that better understanding of protein-RNA interactions would be obtained from ensemble descriptions. Indeed, earlier molecular dynamics simulations of bound states indicated visible dynamics at the RNA-RRM interfaces. Here, we report the first atomistic simulation study of spontaneous binding of short RNA sequences to RRM domains of HuR and SRSF1 proteins. Using a millisecond-scale aggregate ensemble of unbiased simulations, we were able to observe a few dozen binding events. HuR RRM3 utilizes a pre-binding state to navigate the RNA sequence to its partially disordered bound state and then to dynamically scan its different binding registers. SRSF1 RRM2 binding is more straightforward but still multiple-pathway. The present study necessitated development of a goal-specific force field modification, scaling down the intramolecular van der Waals interactions of the RNA which also improves description of the RNA-RRM bound state. Our study opens up a new avenue for large-scale atomistic investigations of binding landscapes of protein-RNA complexes, and future perspectives of such research are discussed.

Multiscale Modeling of Phosphate···π Contacts in RNA U-Turns Exposes Differences between Quantum-Chemical and AMBER Force Field Descriptions.

Phosphate···π, also called anion···π, contacts occur between nucleobases and anionic phosphate oxygens (OP2) in r(GNRA) and r(UNNN) U-turn motifs (N = A,G,C,U; R = A,G). These contacts were investigated using state-of-the-art quantum-chemical methods (QM) to characterize their physicochemical properties and to serve as a reference to evaluate AMBER force field (AFF) performance. We found that phosphate···π interaction energies calculated with the AFF for dimethyl phosphate···nucleobase model systems are less stabilizing in comparison with double-hybrid DFT and that minimum contact distances are larger for all nucleobases. These distance stretches are also observed in large-scale AFF vs QM/MM computations and classical molecular dynamics (MD) simulations on several r(gcGNRAgc) tetraloop hairpins when compared to experimental data extracted from X-ray/cryo-EM structures (res. ≤ 2.5 Å) using the WebFR3D bioinformatic tool. MD simulations further revealed shifted OP2/nucleobase positions. We propose that discrepancies between the QM and AFF result from a combination of missing polarization in the AFF combined with too large AFF Lennard-Jones (LJ) radii of nucleobase carbon atoms in addition to an exaggerated short-range repulsion of the r-12 LJ repulsive term. We compared these results with earlier data gathered on lone pair···π contacts in CpG Z-steps occurring in r(UNCG) tetraloops. In both instances, charge transfer calculations do not support any significant n → π* donation effects. We also investigated thiophosphate···π contacts that showed reduced stabilizing interaction energies when compared to phosphate···π contacts. Thus, we challenge suggestions that the experimentally observed enhanced thermodynamic stability of phosphorothioated r(GNRA) tetraloops can be explained by larger London dispersion.

Automatic Learning of Hydrogen-Bond Fixes in the AMBER RNA Force Field.

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. Here, we 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 gHBfixopt. 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 the 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.

Toward Convergence in Folding Simulations of RNA Tetraloops: Comparison of Enhanced Sampling Techniques and Effects of Force Field Modifications.

Atomistic molecular dynamics simulations represent an established technique for investigation of RNA structural dynamics. Despite continuous development, contemporary RNA simulations still suffer from suboptimal accuracy of empirical potentials (force fields, ffs) and sampling limitations. Development of efficient enhanced sampling techniques is important for two reasons. First, they allow us to overcome the sampling limitations, and second, they can be used to quantify ff imbalances provided they reach a sufficient convergence. Here, we study two RNA tetraloops (TLs), namely the GAGA and UUCG motifs. We perform extensive folding simulations and calculate folding free energies (ΔGfold°) with the aim to compare different enhanced sampling techniques and to test several modifications of the nonbonded terms extending the AMBER OL3 RNA ff. We demonstrate that replica-exchange solute tempering (REST2) simulations with 12-16 replicas do not show any sign of convergence even when extended to a timescale of 120 µs per replica. However, the combination of REST2 with well-tempered metadynamics (ST-MetaD) achieves good convergence on a timescale of 5-10 µs per replica, improving the sampling efficiency by at least 2 orders of magnitude. Effects of ff modifications on ΔGfold° energies were initially explored by the reweighting approach and then validated by new simulations. We tested several manually prepared variants of the gHBfix potential which improve stability of the native state of both TLs by ∼2 kcal/mol. This is sufficient to conveniently stabilize the folded GAGA TL while the UUCG TL still remains under-stabilized. Appropriate adjustment of van der Waals parameters for C-H···O5' base-phosphate interaction may further stabilize the native states of both TLs by ∼0.6 kcal/mol.

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

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

W-RESP: Well-Restrained Electrostatic Potential-Derived Charges. Revisiting the Charge Derivation Model.

Representation of electrostatic interactions by a Coulombic pairwise potential between atom-centered partial charges is a fundamental and crucial part of empirical force fields used in classical molecular dynamics simulations. The broad success of the AMBER force-field family originates mainly from the restrained electrostatic potential (RESP) charge model, which derives partial charges to reproduce the electrostatic field around the molecules. However, the description of the electrostatic potential around molecules by standard RESP may be biased for some types of molecules. In this study, we modified the RESP charge derivation model to improve its description of the electrostatic potential around molecules and thus electrostatic interactions in the force field. In particular, we reoptimized the atomic radii for definition of the grid points around the molecule, redesigned the restraining scheme, and included extra point (EP) charges. The RESP fitting was significantly improved for aromatic heterocyclic molecules. Thus, the suggested W-RESP(-EP) charge derivation model shows some potential for improving the performance of the nucleic acid force fields, for which the poor description of nonbonded interactions, such as the underestimated stability of base pairing, is well-established. We also report some preliminary simulation tests (around 1 ms of simulation data) on A-RNA duplexes, tetranucleotides, and tetraloops. The simulations reveal no adverse effects, while the description of base-pairing interactions might be improved. The new charges can thus be used in future attempts to improve the nucleic acid simulation force fields, in combination with reparametrization of the other terms.

UUCG RNA Tetraloop as a Formidable Force-Field Challenge for MD Simulations.

Explicit solvent atomistic molecular dynamics (MD) simulations represent an established technique to study structural dynamics of RNA molecules and an important complement for diverse experimental methods. However, performance of molecular mechanical (MM) force fields (ff's) remains far from satisfactory even after decades of development, as apparent from a problematic structural description of some important RNA motifs. Actually, some of the smallest RNA molecules belong to the most challenging systems for MD simulations and, among them, the UUCG tetraloop is saliently difficult. We report a detailed analysis of UUCG MD simulations, depicting the sequence of events leading to the loss of the UUCG native state during MD simulations. The total amount of MD simulation data analyzed in this work is close to 1.3 ms. We identify molecular interactions, backbone conformations, and substates that are involved in the process. Then, we unravel specific ff deficiencies using diverse quantum mechanical/molecular mechanical (QM/MM) and QM calculations. Comparison between the MM and QM methods shows discrepancies in the description of the 5'-flanking phosphate moiety and both signature sugar-base interactions. Our work indicates that poor behavior of the UUCG tetraloop in simulations is a complex issue that cannot be attributed to one dominant and straightforwardly correctable factor. Instead, there is a concerted effect of multiple ff inaccuracies that are coupled and amplifying each other. We attempted to improve the simulation behavior by some carefully tailored interventions, but the results were still far from satisfactory, underlying the difficulties in development of accurate nucleic acid ff's.

Fine-Tuning of the AMBER RNA Force Field with a New Term Adjusting Interactions of Terminal Nucleotides.

Determination of RNA structural-dynamic properties is challenging for experimental methods. Thus, atomistic molecular dynamics (MD) simulations represent a helpful technique complementary to experiments. However, contemporary MD methods still suffer from limitations of force fields (ffs), including imbalances in the nonbonded ff terms. We have recently demonstrated that some improvement of state-of-the-art AMBER RNA ff can be achieved by adding a new term for H-bonding called gHBfix, which increases tuning flexibility and reduces risk of side-effects. Still, the first gHBfix version did not fully correct simulations of short RNA tetranucleotides (TNs). TNs are key benchmark systems due to availability of unique NMR data, although giving too much weight on improving TN simulations can easily lead to overfitting to A-form RNA. Here we combine the gHBfix version with another term called tHBfix, which separately treats H-bond interactions formed by terminal nucleotides. This allows to refine simulations of RNA TNs without affecting simulations of other RNAs. The approach is in line with adopted strategy of current RNA ffs, where the terminal nucleotides possess different parameters for terminal atoms than the internal nucleotides. Combination of gHBfix with tHBfix significantly improves the behavior of RNA TNs during well-converged enhanced-sampling simulations using replica exchange with solute tempering. TNs mostly populate canonical A-form like states while spurious intercalated structures are largely suppressed. Still, simulations of r(AAAA) and r(UUUU) TNs show some residual discrepancies with primary NMR data which suggests that future tuning of some other ff terms might be useful. Nevertheless, the tHBfix has a clear potential to improve modeling of key biochemical processes, where interactions of RNA single stranded ends are involved.

Improving the Performance of the Amber RNA Force Field by Tuning the Hydrogen-Bonding Interactions.

Molecular dynamics (MD) simulations became a leading tool for investigation of structural dynamics of nucleic acids. Despite recent efforts to improve the empirical potentials (force fields, ffs), RNA ffs have persisting deficiencies, which hamper their utilization in quantitatively accurate simulations. Previous studies have shown that at least two salient problems contribute to difficulties in the description of free-energy landscapes of small RNA motifs: (i) excessive stabilization of the unfolded single-stranded RNA ensemble by intramolecular base-phosphate and sugar-phosphate interactions and (ii) destabilization of the native folded state by underestimation of stability of base pairing. Here, we introduce a general ff term (gHBfix) that can selectively fine-tune nonbonding interaction terms in RNA ffs, in particular, the H bonds. The gHBfix potential affects the pairwise interactions between all possible pairs of the specific atom types, while all other interactions remain intact; i.e., it is not a structure-based model. In order to probe the ability of the gHBfix potential to refine the ff nonbonded terms, we performed an extensive set of folding simulations of RNA tetranucleotides and tetraloops. On the basis of these data, we propose particular gHBfix parameters to modify the AMBER RNA ff. The suggested parametrization significantly improves the agreement between experimental data and the simulation conformational ensembles, although our current ff version still remains far from being flawless. While attempts to tune the RNA ffs by conventional reparametrizations of dihedral potentials or nonbonded terms can lead to major undesired side effects, as we demonstrate for some recently published ffs, gHBfix has a clear promising potential to improve the ff performance while avoiding introduction of major new imbalances.

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.

Molecular Dynamics Simulations Reveal an Interplay between SHAPE Reagent Binding and RNA Flexibility.

The function of RNA molecules usually depends on their overall fold and on the presence of specific structural motifs. Chemical probing methods are routinely used in combination with nearest-neighbor models to determine RNA secondary structure. Among the available methods, SHAPE is relevant due to its capability to probe all RNA nucleotides and the possibility to be used in vivo. However, the structural determinants for SHAPE reactivity and its mechanism of reaction are still unclear. Here molecular dynamics simulations and enhanced sampling techniques are used to predict the accessibility of nucleotide analogs and larger RNA structural motifs to SHAPE reagents. We show that local RNA reconformations are crucial in allowing reagents to reach the 2'-OH group of a particular nucleotide and that sugar pucker is a major structural factor influencing SHAPE reactivity.