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.

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. 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-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 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 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.

The YTH domain of YTHDC1 belongs to a class of protein "readers", recognizing the N6methyladenosine (m 6 A) chemical modification in mRNA. Static ensemble-averaged structures revealed details of N6-methyl recognition via a conserved aromatic cage. Here, we performed molecular dynamics (MD) simulations along with nuclear magnetic resonance (NMR) and isothermal titration calorimetry (ITC) to examine how dynamics and solvent interactions contribute to the m 6 A recognition and negative selectivity towards unmethylated substrate. The structured water molecules surrounding the bound RNA and the methylated substrate's ability to exclude bulk water molecules contribute to the YTH domain's preference for m 6 A. Intrusions of bulk water deep into the binding pocket disrupt binding of unmethylated adenosine. The YTHDC1's preference for the 5′-Gm 6 A-3′ motif is partially facilitated by a network of water-mediated interactions between the 2-amino group of the guanosine and residues in the m 6 A binding pocket. The 5′-Im 6 A-3′ (where I is inosine) motif can be recognized too but disruption of the water network lowers affinity. The D479A mutant also disrupts the water network and destabilizes m 6 A binding. Our interdisciplinary study of YTHDC1 protein/RNA complex reveals an unusual physical mechanism by which solvent interactions contributes towards m 6 A recognition.