Search

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.

Phosphorothioates (PTs) are important chemical modifications of the RNA backbone where a single nonbridging oxygen of the phosphate is replaced with a sulfur atom. PT can stabilize RNAs by protecting them from hydrolysis and is commonly used as a tool to explore their function. It is, however, unclear what basic physical effects PT has on RNA stability and electronic structure. Here, we present molecular dynamics (MD) simulations, quantum mechanical (QM) calculations, and NMR spectroscopy measurements, exploring the effects of PT modifications in the structural context of the neomycin-sensing riboswitch (NSR). The NSR is the smallest biologically functional riboswitch with a well-defined structure stabilized by a U-turn motif. Three of the signature interactions of the U-turn: an H-bond, an anion-π interaction, and a potassium binding site; are formed by RNA phosphates, making the NSR an ideal model for studying how PT affects RNA structure and dynamics. By comparing with high-level QM calculations, we reveal the distinct physical properties of the individual interactions facilitated by the PT. The sulfur substitution, besides weakening the direct H-bond interaction, reduces the directionality of H-bonding while increasing its dispersion and induction components. It also reduces the induction and increases the dispersion component of the anion-π stacking. The sulfur force-field parameters commonly employed in the literature do not reflect these distinctions, leading to the unsatisfactory description of PT in simulations of the NSR. We show that it is not possible to accurately describe the PT interactions using one universal set of van der Waals sulfur parameters and provide suggestions for improving the force-field performance.

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.

We present the first, to the best of our knowledge, ab initio molecular dynamics (AIMD) investigation on three aqueous solutions where an abasic cyclic nucleotide model is solvated in the presence of distinct cations (i.e., Na + , K + and Mg 2+ ). We elucidate the typical modalities of interaction between those ionic species and the nucleotide moiety by first-principles numerical simulations, starting from an inner-shell binding configuration on a time scale of 100 ps (total simulation time of ∼600 ps). Whereas the strong "structure-maker" Mg 2+ is permanently bound to one of the two oxygen atoms of the phosphate group of the nucleotide model, Na + and K + show binding times τ b of 65 ps and 10-15 ps, respectively, thus reflecting their chemical nature in aqueous solutions. Furthermore, we qualitatively relate these findings to approximate free-energy barriers of the cations' unbinding obtained by means of exploratory well-tempered metadynamics. With the aim of shedding light on the features of commonly employed force-fields (FFs), classical MD simulations (almost 200 trajectories with a total simulation time of ∼18 μs) using the biomolecular AMBER FF are also reported. By choosing several combinations of the parametrization for the water environment (i.e., TIP3P, SPC/E and OPC) and cations (i.e., Joung-Cheatham, Li-Merz 12-6 and Li-Merz 12-6-4), we found significant differences in the radial distribution functions and residence times compared to the ab initio results. The Na + and K + ions wrongly show quasi-identical radial distribution functions and the Li & Merz 12-6-4 Lennard-Jones parameters for Mg 2+ were found to be essential in quickly reaching the binding state consistent with AIMD.

Although atomistic explicit-solvent Molecular Dynamics (MD) is a popular tool to study protein-RNA recognition, satisfactory MD description of protein-RNA complexes is not always achieved. Unfortunately, it is often difficult to separate MD simulation instabilities primarily caused by the simple point-charge molecular mechanics (MM) force fields from problems related to the notorious uncertainties in the starting structures. Herein, we report a series of large-scale QM/MM calculations on the U1A protein-RNA complex. This experimentally well-characterized system has an intricate protein-RNA interface, which is very unstable in MD simulations. The QM/MM calculations identify several H-bonds poorly described by the MM method and thus indicate the sources of instabilities of the U1A interface in MD simulations. The results suggest that advanced QM/MM computations could be used to indirectly rationalize problems seen in MM-based MD simulations of protein-RNA complexes. As the most accurate QM method, we employ the computationally demanding meta-GGA density functional TPSS-D3(BJ)/def2-TZVP level of theory. Because considerably faster methods would be needed to extend sampling and to study even larger protein-RNA interfaces, a set of low-cost QM/MM methods is compared to the TPSS-D3(BJ)/def2-TZVP data. The PBEh-3c and B97-3c density functional composite methods appear to be suitable for protein-RNA interfaces. In contrast, HF-3c and the tight-binding Hamiltonians DFTB3-D3 and GFN-xTB perform unsatisfactorily and do not provide any advantage over the MM description. These conclusions are supported also by similar analysis of a simple HutP protein-RNA interface, which is well-described by MD with the exception of just one H-bond. Some other methodological aspects of QM/MM calculations on protein-RNA interfaces are discussed.

The HutP protein from B. subtilis regulates histidine metabolism by interacting with an antiterminator mRNA hairpin in response to the binding of l-histidine and Mg 2+ . We studied the functional ligand-bound HutP hexamer complexed with two mRNAs using all-atom microsecond-scale explicit-solvent MD simulations performed with the Amber force fields. The experimentally observed protein-RNA interface exhibited good structural stability in the simulations with the exception of some fluctuations in an unusual adenine-threonine interaction involving two closely spaced H-bonds. We further investigated this interaction by comparing QM/MM and MM optimizations, using the QM region comprising almost 350 atoms described at the DFT-D3 level. The QM/MM method clearly improved the adenine-threonine interaction compared to MM, especially when the X-H bond lengths were frozen during the MM optimization to mimic the use of SHAKE in the MD simulations. Thus, both the MM approximation and the use of SHAKE can compromise the description of H-bonds at protein-RNA interfaces. The simulations also revealed a notable Mg 2+ -parameter dependence in the behavior of the ligand-binding pocket (LBP). With the SPC/E water model, the 12-6 Åqvist and Li&Merz parameters provided an entirely stable LBP structure, but the 12-6 Allnér and 12-6-4 Li&Merz parametrizations resulted in a progressive loss of direct nitrogen-Mg 2+ LBP coordination. The Allnér and Li&Merz 12-6 parametrizations were also tested with the TIP3P water model; the LBP was destabilized in both cases. This illustrates the difficulty of consistently describing different Mg 2+ interactions using nonpolarizable force fields. Overall, the simulations support the hypothesis that HutP protein becomes fully structured upon ligand binding. Subsequent RNA binding does not affect the protein structure, in keeping with the mechanism inferred from experimental structures.

Recent developments in dispersion-corrected density functional theory methods allow for the first time the description of large fragments of nucleic acids (hundreds of atoms) with an accuracy clearly surpassing the accuracy of common biomolecular force fields. Such calculations can significantly improve the description of the potential energy surface of nucleic acid molecules, which may be useful for studies of molecular interactions and conformational preferences of nucleic acids, as well as verification and parameterization of other methods. The first of such studies, however, demonstrated that successful applications of accurate QM calculations to larger nucleic acid building blocks are hampered by difficulties in obtaining geometries that are biochemically relevant and are not biased by non-native structural features. We present an approach that can greatly facilitate large-scale QM studies on nucleic acids, namely electronic structure geometry optimization of nucleic acid fragments utilizing a penalty function to restrain key internal coordinates with a specific focus on the torsional backbone angles. This work explores the viability of these restraint optimizations for DFT-D3, PM6-D3H and HF-3c optimizations on a set of examples (a UpA dinucleotide, a DNA G-quadruplex and a B-DNA fragment). Evaluation of different penalty function strengths reveals only a minor system-dependency and reasonable restraint values range from 0.01 to 0.05 Eh rad(-2) for the backbone torsions. Restraints are crucial to perform the QM calculations on biochemically relevant conformations in implicit solvation and gas phase geometry optimizations. The reasons for using restrained instead of constrained or unconstrained optimizations are explained and an open-source external optimizer is provided.