Simple adjustment of intranucleotide base-phosphate interaction in the OL3 AMBER force field improves RNA simulations

Simple Adjustment of Intranucleotide Base-Phosphate Interaction in the OL3 AMBER Force Field Improves RNA Simulations Vojtech Mlynsky,*Petra Kuhrová, Petr Stadlbauer, Miroslav Krepl, Michal Otyepka, Pavel Banás,* and JiríSponer* Cite This: J. Chem. Theory Comput. 2023, 19, 8423−8433 Read Online ACCESS Metrics & More Article Recommendations * sı Supporting Information ABSTRACT: 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. ■INTRODUCTION Insights into RNA structural dynamics at an atomistic description are essential for the understanding of biomolecular motions and processes. Molecular dynamics (MD) simulation is an established theoretical tool that benefits from the ability to overcome experimental limits when finding links between the RNA structure, dynamics, and function. 1−5 Compact folded RNA molecules are typically well-described by modern empirical potentials (force fields, ffs) on a sub-microsecond time scale when starting simulations from established experimental structures. The general usage, applicability, and reproducibility of MD simulations are continuously increasing. However, the MD simulation studies remain limited by the approximative nature of the ffs. A number of RNA ffs are presently available to carry out MD simulations of RNA systems, 6−11 but none of them is flawless. 5,12−18 Despite the suboptimal performance for certain motifs, mainly the short RNA single strands, 14,16,19−21 the AMBER OL3 ff from 2010 7 still remains the state-of-the-art for RNA simulations and safest option to start with. 22 It has been shown in the past decade that the OL3 performance can be partially tweaked by optimizing simulation settings and by adding new ff terms that are orthogonal to the current ff terms; namely, (i) an update of parameters for phosphate oxygens by Steinbrecher and co-workers 23 was shown to improve the general simulation outcome, 24,25 (ii) the combination of OL3 with the four-point OPC water model 26 revealed benefits for structural description of short RNA single-strands (tetranucleotides, TNs) and tetraloops (TLs), 20,25,27 (iii) controllable fine-tuning of specific pairwise H-bond interactions via an external general H-bond fix potential (gHBfix) corrected structural ensembles of several RNA motifs, 13,14,28 and (iv) specific adjustment of interactions formed by terminal nucleotides via tHBfix further improved the agreement with experiments for RNA TNs. 29 In our recent work, 30 we spotted another minor but visible imbalance of the AMBER ff connected with the imprecise description of RNA intranucleotide interactions when van der Waals (vdW) clash in GS+1(C8−H8)···GS+1(O5′) interaction resulted in the spurious flipping of the GS+1 sugar−phosphate backbone and population of non-native structures of UUCG TL. The identified vdW clash in UUCG TL 30 is related to the weak −CH···O−H-bond between −C8H8 and −C6H6 Received: September 9, 2023 Revised: October 18, 2023 Accepted: October 23, 2023 Published: November 9, 2023 Articlepubs.acs.org/JCTC © 2023 The Authors. Published by American Chemical Society 8423 https://doi.org/10.1021/acs.jctc.3c00990 J. Chem. Theory Comput. 2023, 19, 8423−8433 This article is licensed under CC-BY 4.0 Downloaded via TECHL UNIV OF OSTRAVA on June 10, 2024 at 07:21:18 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. groups of purines and pyrimidines, respectively, and bridging O5′phosphate oxygen from the same residue, which was categorized as base-phosphate interaction type 0 (0BPh, Figure 1). 31 The attempt to adjust the 0BPh interaction in the structural context of UUCG TL via nonbonded fix (NBfix) stabilized the native GS+1 phosphate conformation 30 and indicated that it could be beneficial for the stability of native states of both GAGA and UUCG TLs. 32 In this work, we performed an extensive set of standard and enhanced sampling MD simulations on diverse RNA systems and considered the adjustment of 0BPh interaction via NBfix (NBfix0BPh modification) as a general improvement of the AMBER OL3 RNA ff (and likely other ff variants using the original AMBER nonbonded parameters). We first used the reweighting approach to optimize settings for NBfix0BPh modification. Reweighting methods are one of the strongest tools of ff development 27,33 and can re-evaluate the results of MD simulations under the assumption of a modified ff parametrization without repeating the entire simulation. Reweighting allowed us to identify the most suitable modified parameters for both purines and pyrimidines, which were subsequently validated on benchmark systems through largescale MD simulations. The testing set contained singlestranded RNAs, duplexes, TLs, and some examples of more structured RNA motifs with noncanonical interactions. We also investigated possible effects of the NBfix0BPh modification (with a similar setting as those for RNA nucleotides) on DNA duplex and guanine quadruplexes (GQs). We show that the NBfix0BPh modification improves the structural behavior of small RNA motifs by increasing the agreement with experimental data sets. The stability of canonical A-RNA duplexes was also improved, whereas other bigger and more structured RNA motifs were not affected by the NBfix0BPh modification on the affordable simulation time scale. In addition, this paper presents an extensive set of new simulations with two versions of the external gHBfix potential and supports previous claims 14,28 that gHBfix is beneficial for proper structural description of RNA motifs. ■METHODS Details about AMBER RNA and DNA ffs and Their Adjustments. We used standard AMBER OL3 (known also as χOL3) 7,34−36 and AMBER OL15 37 ffs for RNA and DNA simulations, respectively. OL3 RNA ff was further adjusted by the van der Waals (vdW) modification of phosphate oxygens developed by Steinbrecher et al., 23 where the affected dihedrals were adjusted as described elsewhere. 13,24 This RNA ff version is abbreviated as OL3CP henceforth, and the AMBER library file can be found in the Supporting Information of ref 13. In addition, we used gHBfix19 potential, which is the correction for the OL3CP ff from 2019, where all −NH···N−base−base interactions are strengthened by 1.0 kcal/mol and all −OH··· bO−and −OH···nbO−sugar−phosphate interactions are weakened by 0.5 kcal/mol. 14 We also tested the latest optimized gHBfix version from 2021 (gHBfix21), where all RNA H-bond donor···H-bond acceptor interactions are modified; i.e., base donor−base acceptor, base donor−sugar acceptor, base donor−phosphate acceptor, sugar donor−base acceptor, sugar donor−sugar acceptor, and sugar donor− phosphate acceptor are adjusted specifically (see ref 28 for a full description). Simulations of RNA single strands were also performed with the tHBfix20 potential, where additional correction is added to interactions formed by terminal residues. 29 For some specific tests, we also used standard OL3 (i.e., without phosphate modification) and older ff 99bsc0 36 RNA ffs. On top of that, we modified the pairwise vdW parameters via breakage of the combination (mixing) rules via the nonbonded fix (NBfix) approach 38 for atoms involved in 0BPh intranucleotide interactions for both RNA OL3CP and DNA OL15 ffs. Namely, we reduced the minimum-energy distance of the Lennard−Jones potential (i.e., Ri,j parameter) for the −H8··· O5′− and −H6···O5′− pairs, i.e., between H5−OR/OS and H4−OR/OS atom types, by 0.25 Å to 2.8808 and 2.9308 Å for purine and pyrimidine nucleotides, respectively (modification labeled as NBfix0BPh; see Table S1 in Supporting Information for the list of original and modified vdW parameters). The OR atom type belongs to O5′and O3′bridging phosphate oxygens Figure 1. Definition and description of the weak −CH···O−H-bond that is classified as intranucleotide base-phosphate interaction type 0 (0BPh). 31 Panels illustrate examples of 0BPh interactions with established −C6H6···O5′− and −C8H8···O5′− H-bonds for pyrimidines (C3, panel A) and purines (G10, panel B), respectively. The 0BPh is not formed in syn orientation of the nucleobase due to significantly higher −H6/H8···O5′− distances (G9, panel C). Snapshots of all three nucleotides were taken from the NMR structure of UUCG TL (PDB ID 2KOC), 39 and measured distances between H6/H8 and O5′atoms are shown as averages (with errors as standard deviations) over 20 deposited structures. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c00990 J. Chem. Theory Comput. 2023, 19, 8423−8433 8424 in the RNA OL3CP, whereas the OS atom type involves O5′, O3′, and also O4’ atoms from deoxyribose sugar in the DNA OL15 ff. Depths of the potential well (εi,j parameters) were kept at a default value of 0.0505 kcal/mol. In other words, we attempted to decrease the repulsion between H8 and H6 atoms of all purine and pyrimidine bases and O5′oxygens of phosphates. We note that O3′and O4’ atoms have higher distances from H8 and H6 in comparison with O5′atoms, and thus, modified Lennard−Jones potentials for H8/H6···O3′/ O4’ atom pairs are supposed to have marginal effects. Starting Structures and Simulation Setup. Initial coordinates of r(AAAA), r(CAAU), r(CCCC), r(GACC), and r(UUUU) tetranucleotides (TNs); r(UCAAUC) and r(UCUCGU) hexanucleotides (HNs); and r(gcGAGAgc) 8mer tetraloop (GAGA TL) were prepared using the Nucleic Acid Builder of AmberTools14 40 as one strand of an A-form duplex. RNA single strands were solvated using a cubic box of the OPC water 26 with a minimum distance between box walls and solute of 12 Å, yielding ∼2200 water molecules added (∼40 ×40 ×40 Å3box size), ∼4100 water molecules added (∼50 ×50 ×50 Å3box size), and ∼7600 water molecules added (∼60 ×60 ×60 Å3box size) for TNs, HNs, and the GAGA TL, respectively. Enhanced sampling MD simulations were performed in ∼0.15 M KCl salt excess (TNs and HNs) and ∼1.0 M KCl salt excess (GAGA TL) using Joung− Cheatham (JC) ionic parameters 41 for the TIP4P-EW water. Starting topologies and coordinates of two RNA duplexes (PDB ID 1QC0; 42 only 10 canonical base pairs were considered, and PDB ID 1RNA 43 ), RNA sarcin−ricin loop (SRL, PDB ID 3DW4; 44 residues 2649−2671), RNA kink-turn 7 (Kt-7, PDB ID 1S72; 45 residues 76−83 and 91−101), DNA duplex (PDB ID 1BNA 46 ), and two guanine quadruplexes (GQs) formed by the human telomeric sequence (parallelstranded GQ, PDB ID 1KF1 47 and (3 + 1) hybrid GQ, PDB ID 2GKU 48 ) were prepared from particular experimental structures by using the tLEaP module of the AMBER 16 program package 49 (see the Supporting Information of ref 14 for details about structure preparation). A short r(CGCG)2 duplex was prepared using the Nucleic Acid Builder of AmberTools14. Standard MD simulations were carried out in a cubic box of OPC 26 and SPC/E 50 water models (for RNA and DNA simulations, respectively) with a minimum distance between box walls and solute of 10 Å and with ∼0.15 M KCl salt using the JC ionic parameters. The TIP3P water model 51 was also used for some specific tests involving RNA duplexes (Table S2 in the Supporting Information). All MD simulations were run at T= 298 K with the hydrogen mass repartitioning 52 allowing a 4 fs integration time step (see the Supporting Information of ref 14 for other details about minimization and equilibration protocols). Standard MD simulations were run in AMBER18, 53 whereas both AMBER18 and GROMACS2018 54 were used for enhanced sampling simulations. PARMED 55 was used to convert AMBER topologies and coordinates into GROMACS inputs. Enhanced Sampling Simulations. We used two different enhanced sampling schemes, i.e., a standard replica exchange solute tempering (REST2) protocol 56 and well-tempered metadynamics 57−59 (MetaD) in combination with the REST2 method (ST-MetaD). 32,60 REST2 simulations were performed at 298 K (the reference replica) with 8 and 12 replicas for TNs and HNs, respectively. The scaling factor (λ) values ranged from 1 to 0.601700871 and to 0.59984 for 8 and 12 replicas, respectively. Those values were chosen to maintain an exchange rate above 20%. The effective solute temperature ranged from 298 to ∼500 K. REST2 simulations were performed with the AMBER GPU MD simulation engine (pmemd.cuda). 61 Further details about REST2 settings can be found elsewhere. 14 One test simulation of r(UCUCGU) HN was performed at 275 K temperature corresponding to the experimental conditions where the biggest number of NMR signals was obtained. 16 The same λvalues were applied for scaling, which resulted in an effective solute temperature range from 275 to ∼460 K. We note that the 275 and 298 K REST2 simulations revealed comparable results considering the limits of sampling. ST-MetaD simulations of GAGA TL were performed with 12 replicas starting from unfolded single strands and were simulated in the effective temperature range of 298−497 K for 5μs per replica. The average acceptance rate was ∼30%. The εRMSD metric 62 was used as a biased collective variable. 32,63 ST-MetaD simulations were carried out using a GPU-capable version of GROMACS2018 54 in combination with PLUMED 2.5 64,65 (see ref 32 for further details about ST-MetaD settings). Besides newly performed MD simulations, we also used some trajectories from our previous works (see Table S2 for a full list of standard as well as enhanced sampling simulations). The reweighting algorithm, 66 which enables fast reevaluation of the results from MD simulations, was used on the range of modified Ri,j values of the Lennard−Jones potential for the NBfix0BPh modification in the attempt to find the optimal setup. Reweighting was performed on trajectories from both REST2 and ST-MetaD simulations using snapshots only from the reference replica (the lowest REST2 replica corresponding to 298 K; see ref 32 for further details about the application of the reweighting approach). Conformational Analysis. The most populated conformations from TNs and HNs structural ensembles were identified by clustering, which is based on an algorithm introduced by Rodriguez and Laio 67 in combination with the εRMSD metric 62 (see ref 29 for more details). Comparison between MD and NMR Data. The conformational ensembles of TNs and HNs obtained from REST2 simulations were compared with previously published NMR experiments. 16,21,68−70 We analyzed separately four NMR observables, i.e., (i) backbone 3J scalar couplings, (ii) sugar 3J scalar couplings, (iii) nuclear Overhauser effect intensities (NOEs), and (iv) the absence of specific peaks in NOE spectroscopy (uNOEs). In addition, we also considered ambiguous NOEs (ambNOEs) resulting from the sum of overlapping peaks 16,21 as the fifth individual NMR component for HNs. 3J scalar couplings were calculated via Karplus relationships; NOEs and uNOEs were obtained as averages over the Nsamples, i.e., u( )NOE ( ) r N CALC 1/6 i N i 6 = ; and ambNOEs were calculated by summing the contribution from either two or three nuclei pairs and again averaged over the N samples, i.e., ambNOE r r r N CALC (1 2 (3)) 1/6 i N i i i 6 6 6 = + + i k j j jy { z z z . The combination of all those analyzed NMR observables (calculated as weighted arithmetic mean) provided the total χ2value for each REST2 simulation: n n total 2i n ii i n i 2 = , where χi 2 and niare individual χ2values and the number of observables, respectively, corresponding to the particular NMR observables. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c00990 J. Chem. Theory Comput. 2023, 19, 8423−8433 8425 The lower the total χ2value is, the better is the agreement between the experiment and REST2 simulation (for detailed explanations, see, e.g., refs 27,29) ST-MetaD simulation of the GAGA TL provided populations of the native structure and other conformations, which can be used for estimation of the folding free energy balance (ΔG°fold). 28,32 The reference native structure of GAGA TL was taken from our previous work. 14 The εRMSD threshold separating the folded and un(mis)folded states was set at a value of 0.7 (see ref 32 for details about ΔG°fold estimations and convergence). ■RESULTS AND DISCUSSION Reweighting Reveals Optimal Settings for the NBfix0BPh Modification. The reweighting algorithm 66 allows one to efficiently re-evaluate results of MD simulations under the assumption of a modified ff parametrization without the necessity to perform new simulations. 27,32,33 Here, we reweighted the NBfix0BPh settings for a set of enhanced sampling simulations of RNA TNs and TLs. More specifically, we used results from OL3CP REST2 simulations with gHBfix19 and tHBfix20 potentials of r(AAAA), r(CCCC), and r(CAAU) TNs 29 and from folding OL3CP ST-MetaD simulations with gHBfix19 of GAGA and UUCG TLs. 32 We monitored changes in total χ2values (and its components) and ΔG°fold energies for TNs and TLs, respectively. Each NBfix0BPh setting was defined by scanning Ri,j values of the Lennard−Jones potential for the −H8···O5′− and −H6···O5′− atom pairs from 2.33− 3.33 and 2.48−3.48 intervals for purine and pyrimidine nucleotides, respectively, using 0.025 Å steps. Results from TN simulations indicate that a 0.25 Å decrease of Ri,j parameters for both purine and pyrimidine nucleotides is optimal (Figures S1−S5 in the Supporting Information). Reweighted data from r(CCCC) and r(CAAU) simulations with the sole gHBfix19 suggest an even larger decrease (than the finally selected 0.25 Å) of the Ri,j value for pyrimidines, i.e., −H6···O5′−atom pairs, but that is contradicted by data from simulations with combined gHBfix19 + tHBfix20 potentials (those giving much better agreement with experiments; 29 see Figures S2, S4 and S5 in Supporting Information). Data from the reweighting of ST-MetaD TL simulations support results from TNs. Reweighting of the GAGA TL simulation shows that the adjustment of Ri,j values for the −H8···O5′− atom pairs of purine nucleotides is the dominant contributor for shifting ΔG°fold energy (Figure 2A). This is not surprising because six out of eight nucleotides from the r(gcGAGAgc) TL are purines. Data even suggested that a larger decrease would be vital, but that was not confirmed when reweighting Ri,j values for both purines and pyrimidines simultaneously (Figure 2A). Modifications of Ri,j parameters for both purines and pyrimidines significantly affect the stability of the UUCG TL native state, but in the opposite direction (Figure 2B). More specifically, a decrease of the Ri,j value for only pyrimidine residues indicates stabilization; ΔG°fold is decreasing from ∼5.4 kcal/mol (default Ri,j value) up to ∼3.5 kcal/mol (final and lowest tested Ri,j) (Figure 2B). In contrast, decreased Ri,j values for sole purines indicate destabilization (ΔG°fold is increasing from ∼5.4 kcal/mol up to ∼6.0 kcal/mol, Figure 2B). Such a counterintuitive effect of the decreased Ri,j values for the −H8···O5′− atom pairs of purines, i.e., small destabilization of the UUCG native state, is caused by the stabilization of misfolded states, where the key G residue in the loop (GL4) is sampling noncanonical anti orientation of the χdihedral (syn states are required for the native conformation, and those are Figure 2. Effects of different setting of NBfix0BPh modifications on the stability of the native state of RNA TLs. Folding free energies (ΔG°fold) were estimated by reweighting of GAGA (A) and UUCG (B) trajectories from OL3CP ST-MetaD simulations with the gHBfix19 potential. Plots show the dependence of ΔG°fold energies on modified vdW parameters for purines and pyrimidines, i.e., changes of Ri,j values for −H8···O5′− (cyan line) and −H6···O5′− (magenta line) atom pairs from the original (standard AMBER ff) values. The black line presents the outcome, where Ri,j values for both purines and pyrimidines were modified simultaneously by the same value (see Methods for details). The gray vertical line indicates the position of the original (unmodified) values. The NBfix0BPh modification (0.25 Å decrease of pairwise terms; green vertical line) stabilizes native states of both TLs by providing lower ΔG°fold energies, which are in better agreement with experiments (reported ΔG°fold energies are between −0.7 and −1.3 kcal/mol for both TLs 71−73 ). The reweighting data are still suggesting that the UUCG native state is significantly destabilized by the ff even after the application of the NBfix0BPh modification. This is not surprising because multiple ff disbalances affect the stability of UUCG TL (see ref 30 for details). We note that the latest optimized gHBfix version (gHBfix21) achieved a more significant improvement of free energy of the UUCG TL compared to gHBfix19. 28 Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c00990 J. Chem. Theory Comput. 2023, 19, 8423−8433 8426 not affected by the NBfix0BPh modification; see Figure 1 and ref 32 for further discussion). Most importantly, simultaneous modification of Ri,j values for both purines and pyrimidines, i.e., the NBfix0BPh modification, revealed substantial stabilization of the UUCG native state; i.e., modified pyrimidines provide a dominant effect for the UUCG TL (Figure 2B). In summary, reweighted data from both TN and TL motifs showed that decreased Ri,j parameters for both −H8···O5′− and −H6···O5′− atom pairs provided better agreement with experiments. Effective usage of reweighting requires a broad exploration of the conformational space, which we believe was achieved in trajectories from REST2 and ST-MetaD simulations of TNs and TLs, respectively. The accuracy of reweighted results decreases with increasing deviation of the parameter from original values (see, e.g., the continuous decrease of ΔG°fold energy with decreasing Ri,j parameter for pyrimidines on Figure 2B). Those large parameter changes often result in problems (e.g., clashes) in subsequent simulations with modified parameters as they drive the system toward states that were not sampled in original trajectories used for reweighting. Hence, we aimed for a reasonably small change of Ri,j parameters that would be enough to relax the steric clash of −H8···O5′− and −H6···O5′− atom pairs observed in the original OL3CP ff. We found that a 0.25 Å decrease of the particular Ri,j parameter for both purines and pyrimidines appears to be the optimal choice. NBfix0BPh Modification of the AMBER OL3CP RNA ff Improves the Structural Description of Single-Stranded RNA Motifs. We used the optimized setting for the NBfix0BPh modification (0.25 Å decrease) and explicitly tested its effect on the structural description of five RNA TNs and two HNs, i.e., benchmark set of motifs with available NMR data, in a large set of enhanced sampling simulations. We again opted for the REST2 protocol, which was shown to provide ensembles with sufficient convergence for those single-stranded motifs. 14,29 We used OL3CP ff combined with gHBfix19, gHBfix19 + tHBfix20, and gHBfix21 + tHBfix20 potentials, which were shown to improve the structural behavior of the short RNA single strand motifs (see Methods and Table S2 in Supporting Information). The data are summarized in Tables 1 and 2for TNs and HNs, respectively. The NBfix0BPh modification further fine-tunes especially the r(CAAU) and r(AAAA) ensembles, which now agree excellently with NMR data (total χ2value below 1; Table 1). The NBfix0BPh modification indirectly stabilizes the anti state of A residues and thus resolves the previously reported excessive sampling of syn states. 29,74 We also observed that NBfix0BPh increased the population of canonical A-form like states in Table 1. REST2 Simulations of RNA TNs and Their Comparison with Experiments a motif gHBfix version NBfix0BPh χ2(3J-backbone, 3J-sugar, NOE, uNOE (#of violations)/total) b clustering (A-form/loop/bulged-out/ intercalated; in %) #of clusters r(GACC) c gHBfix19 no 0.39, 0.50, 1.13, 0.34 (3)/0.40 ∼74/∼7/∼5/∼2 5 r(GACC) c gHBfix19_tHBfix20 no 0.34, 0.37, 1.11, 0.15 (3)/0.23 ∼81/∼5/∼4/∼0 4 r(GACC) gHBfix21_tHBfix20 no 0.28, 0.39, 2.62, 0.02 (1)/0.20 ∼86/∼2/∼4/∼0 3 r(GACC) gHBfix19 yes 0.29, 0.33, 1.23, 0.00 (0)/0.10 ∼91/∼1/∼3/∼0 4 r(GACC) gHBfix19_tHBfix20 yes 0.29, 0.34, 1.34, 0.00 (1)/0.11 ∼90/∼2/∼4/∼0 4 r(GACC) gHBfix21_tHBfix20 yes 0.28, 0.39, 2.70, 0.00 (0)/0.19 ∼93/∼2/∼2/∼0 3 r(CAAU) c gHBfix19 no 0.73, 1.14, 1.81, 6.07 (28)/5.37 ∼33/∼0/∼3/∼51 5 r(CAAU) c gHBfix19 + tHBfix20 no 0.49, 0.80, 1.40, 2.22 (29)/2.05 ∼76/∼0/∼6/∼3 8 r(CAAU) gHBfix21 + tHBfix20 no 0.47, 1.55, 1.63, 3.78 (36)/3.41 ∼73/∼0/∼9/∼1 5 r(CAAU) gHBfix19 yes 0.55, 0.92, 1.28, 4.88 (21)/4.30 ∼49/∼0/∼0/∼44 3 r(CAAU) gHBfix19 + tHBfix20 yes 0.33, 0.56, 0.89, 0.32 (12)/0.37 ∼89/∼0/∼0/∼5 3 r(CAAU) gHBfix21 + tHBfix20 yes 0.30, 0.51, 1.02, 0.65 (12)/0.66 ∼85/∼0/∼6/∼0 5 r(AAAA) c gHBfix19 no 0.66, 1.23, 1.08, 1.14 (12)/1.11 ∼30/∼0/∼6/∼10 17 r(AAAA) c gHBfix19 + tHBfix20 no 0.68, 1.48, 1.28, 1.06 (11)/1.08 ∼32/∼0/∼5/∼1 18 r(AAAA) gHBfix21 + tHBfix20 no 0.61, 1.84, 1.07, 0.82 (6)/0.87 ∼38/∼0/∼8/∼0 15 r(AAAA) gHBfix19 yes 0.44, 0.83, 1.10, 0.38 (5)/0.48 ∼66/∼0/∼3/∼12 7 r(AAAA) gHBfix19 + tHBfix20 yes 0.42, 0.78, 1.14, 0.02 (2)/0.20 ∼72/∼0/∼3/∼2 8 r(AAAA) gHBfix21 + tHBfix20 yes 0.35, 0.66, 1.21, 0.12 (2)/0.28 ∼70/∼0/∼8/∼0 5 r(CCCC) c gHBfix19 no 0.44, 0.68, 1.54, 3.22 (11)/2.83 ∼55/∼0/∼0/∼40 2 r(CCCC) c gHBfix19 + tHBfix20 no 0.25, 0.38, 1.22, 0.50 (6)/0.55 ∼85/∼0/∼2/∼8 5 r(CCCC) gHBfix21 + tHBfix20 no 0.22, 0.35, 3.38, 0.10 (3)/0.41 ∼86/∼0/∼5/∼3 4 r(CCCC) gHBfix19 yes 0.28, 0.38, 1.28, 1.87 (5)/1.68 ∼76/∼0/∼0/∼23 2 r(CCCC) gHBfix19 + tHBfix20 yes 0.20, 0.30, 1.41, 0.12 (3)/0.25 ∼93/∼0/∼0/∼4 2 r(CCCC) gHBfix21 + tHBfix20 yes 0.18, 0.27, 3.92, 0.02 (1)/0.39 ∼96/∼0/∼0/∼0 1 r(UUUU) c gHBfix19 no 0.53, 0.87, 3.53, 1.65 (19)/1.63 ∼5/∼0/∼0/∼3 10 r(UUUU) c gHBfix19 + tHBfix20 no 0.52, 0.97, 3.46, 1.86 (17)/1.81 ∼14/∼0/∼0/∼0 10 r(UUUU) gHBfix21 + tHBfix20 no 0.47, 1.69, 3.58, 1.47 (11)/1.49 ∼8/∼25/∼0/∼0 9 r(UUUU) gHBfix19 yes 0.48, 0.54, 4.50, 1.72 (19)/1.70 ∼36/∼0/∼2/∼2 9 r(UUUU) gHBfix19 + tHBfix20 yes 0.49, 0.60, 4.55, 2.04 (18)/1.99 ∼29/∼0/∼5/∼0 9 r(UUUU) gHBfix21 + tHBfix20 yes 0.46, 1.06, 4.44, 1.71 (12)/1.71 ∼24/∼27/∼1/∼0 11 a Simulations used a combination of OL3CP RNA ff with either gHBfix19 or gHBfix21 potentials and sometimes also with tHBfix20 (see Methods). REST2 simulations used eight replicas and were run for 10 μs (the last 7 μs was used for data analysis). Experimental data are taken from refs 68−70. b χ2values were obtained by comparing calculated and experimental backbone 3J scalar couplings, sugar 3J scalar couplings, nuclear Overhauser effect intensities (NOEs), and the absence of specific peaks in NOE spectroscopy (uNOEs, the number of violations is reported). See Methods for details. c Data taken from previous papers. 14,29 Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c00990 J. Chem. Theory Comput. 2023, 19, 8423−8433 8427 structural ensembles of all simulated TNs (Table 1), which is expected to be the main conformer based on available experimental data. 21,68−70 Application of NBfix0BPh did not provide better agreement with the experiment for r(UUUU) TN. We did observe an increase (∼20%) in the population of the canonical A-form like structures (as for other TNs), but it was associated with a slight increase of total χ2values; i.e., the agreement between simulations and experiment was modestly worsened for all tested MD setups (Table 1). Experimental data show that r(UUUU) is more dynamical than other TNs and should sample unstacked/disordered states. 68,75 Hence, further adjustment of AMBER RNA ffs is probably required potentially including revision of stacking interactions 74 to improve the behavior of poly-U sequences in MD simulations. Recently, modification of the AMBER potential (named STAfix) by major weakening of intramolecular stacking was employed to study the spontaneous binding of RNA single strands to RNA recognition motif proteins. 17 It efficiently eliminated majority of previously identified spurious r(UUUU) 1-3_2-4 stacked conformers 29 and supported the sampling of unstructured states in structural ensembles, which resulted in excellent agreement with the NMR data. 17 STAfix, however, is not suitable for simulations of folded RNAs; it is a specialpurpose modification that has been prepared for binding of unstructured single-stranded RNAs to proteins. We also note that we have significantly fewer NMR signals for r(UUUU) than for other TNs, which complicate predictions from the experiment. The NBfix0BPh modification improves the structural behavior of two r(UCAAUC) and r(UCUCGU) HN motifs. It has been shown previously by both experiments and simulations that r(UCAAUC) should sample predominantly structures around the canonical A-form state. 21,76 On the other hand, the presence of the G residue in r(UCUCGU) HN significantly affects the conformational space, and this HN motif is much more dynamical than shorter RNA TNs. 16 r(UCAAUC) REST2 simulations performed here show that the NBfix0BPh modification slightly enhanced the population of canonical Aform like states for r(UCAAUC) HN, and results are in good agreement with the experiment 21 (total χ2values close to 1; Table 2). r(UCUCGU) REST2 simulations revealed disagreement with the experiment, with significantly higher total χ2 values than obtained for other tested single-stranded motifs. As the NBfix0BPh modification did not appear to significantly affect r(UCUCGU) structural ensembles, we performed simulations with only four different setups for this motif (Table 2). The gHBfix21 potential provided slightly worse agreement with the experiment 16 than gHBfix19 as it increased the population of looplike structures. We also performed one REST2 simulation at a lower temperature (where the biggest number of NMR signals was obtained; see Methods), but it did not improve the agreement (Table 2). We made some initial attempts to identify the source of this huge discrepancy between experiment and theory (see the next section). In summary, the NBfix0BPh modification positively affects structural ensembles of single-stranded RNA motifs. Most importantly, it resolves the previously reported excessive sampling of syn states for A residues 29,74 by stabilizing anti states. r(UUUU) TN and especially r(UCUCGU) HN motifs remain challenging for the AMBER OL3CP ff even with the additional refinements. Despite some uncertainty in experimental data sets for those motifs, it is evident that a larger reparameterization than a simple adjustment of pairwise terms for one interaction would be required to improve the structural behavior of those motifs. Taking all of the available data together clearly demonstrates that we remain far from having flawless ff for RNA molecules. Measured NMR signals for r(UCUCGU) HN Are from Multiple Structures. r(UCUCGU) REST2 simulations revealed a striking disagreement between MD predictions and the experimental data. 16 We observed that calculated χ2 values are almost an order of magnitude higher than those from TN simulations and the other r(UCAAUC) HN (Tables 1and 2). Inspection of individual χ2components (see Methods) revealed that the largest deviations between Table 2. REST2 Simulations of RNA HNs and Their Comparison with Experiments a motif gHBfix version NBfix0BPh χ2(3J-backbone, 3J-sugar, NOE, ambNOE, uNOE (#of violations)/total) b clustering (A-form/loop/bulged-out/ intercalated; in %) #of clusters r(UCAAUC) gHBfix19 no 0.52, 3.48, 1.44, 1.92, 1.73 (27)/1.70 ∼26/∼0/∼24/∼2 11 r(UCAAUC) gHBfix19 + tHBfix20 no 0.54, 3.34, 1.40, 1.95, 1.46 (24)/1.46 ∼21/∼0/∼32/∼2 13 r(UCAAUC) gHBfix21 + tHBfix20 no 0.45, 5.58, 1.96, 1.61, 2.27 (33)/2.22 ∼7/∼12/∼49/∼0 13 r(UCAAUC) gHBfix19 yes 0.39, 1.95, 0.91, 0.98, 1.21 (17)/1.16 ∼39/∼0/∼51/∼2 6 r(UCAAUC) gHBfix19 + tHBfix20 yes 0.38, 1.80, 0.91, 0.98, 1.19 (16)/1.15 ∼41/∼0/∼36/∼2 9 r(UCAAUC) gHBfix21 + tHBfix20 yes 0.33, 2.85, 1.15, 1.14, 1.69 (21)/1.61 ∼28/∼0/∼57/∼0 7 r(UCUCGU) gHBfix21 + tHBfix20 no 4.02, 32.87, 7.01, 4.18, 6.25 (1)/11.60 ∼3/∼79/∼3/∼0 6 r(UCUCGU) gHBfix19 + tHBfix20 yes 3.38, 17.54, 2.83, 3.38, 2.25 (1)/6.05 ∼34/∼48/∼12/∼0 8 r(UCUCGU) c gHBfix19 + tHBfix20 yes 3.62, 20.94, 3.79, 4.00, 4.00 (1)/7.39 ∼33/∼57/∼4/∼0 8 r(UCUCGU) gHBfix21 + tHBfix20 yes 3.88, 29.12, 5.90, 6.66, 6.25 (1)/10.43 ∼8/∼79/∼3/∼0 8 a Simulations used a combination of OL3CP RNA ff with either gHBfix19 or gHBfix21 potentials and sometimes also with tHBfix20 (see Methods). REST2 simulations used 12 replicas and were run for 10 μs (the last 7 μs was used for data analysis). Experimental data are taken from refs 16,21. b χ2values were obtained by comparing calculated and experimental backbone 3J scalar couplings, sugar 3J scalar couplings, NOEs, ambiguous NOEs (ambNOEs), and uNOEs (the number of violations is also reported). See Methods for details. c REST2 simulation run with unbiased replica shifted from 298 to 275 K (to mimic the experimental temperature). Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c00990 J. Chem. Theory Comput. 2023, 19, 8423−8433 8428 experimentally measured and predicted results from simulations come from sugar 3J scalar couplings. Experimental data suggest that each residue is sampling a mixture of C3′-endo and C2′-endo sugar puckers. 16 Indeed, MD is able to sample transitions between these most common RNA pucker conformations, but overall populations are shifted from those measured experimentally for all residues (with the G5 residue showing the largest deviation; see Table S3 in Supporting Information). Analysis of most populated conformers during REST2 simulation showed that besides the canonical A-formlike states, MD favors formation of a looplike conformer with formation of G5C2 base pair (Figure S6 in Supporting Information). The experiment shows much weaker NMR Aform signals than for TN and r(UCAAUC) HN motifs, which indicates the presence of another state. 16 However, the formation of conformers with the G5C2 base pair is not supported. 16 In summary, r(UCUCGU) NMR data are composed of different conformers, but those sampled during REST2 simulations (besides the canonical A-form-like states) are apparently different from those detected in the experiment. We speculate that a different preference for the C2′-endo sugar pucker in the AMBER OL3CP ff simulation could also (to some extent) affect the structural ensemble of the r- (UCUCGU) HN. This is out of the scope of this work and will be addressed in the future. The Native State of GAGA TL Is Stabilized by the NBfix0BPh Modification. We performed ST-MetaD simulation of the r(gcGAGAgc) TL to explicitly test the effect of the NBfix0BPh modification on the folding of this small 8-mer RNA TL motif. We obtained a ΔG°fold energy of −0.5 ±0.6 kcal/mol (corresponding to 68.3 ±20.0% of the population of the native state). The result is in an excellent agreement with data from the reweighting approach (see Figure 2 and ref 32). Despite the fact that the “true” uncertainty in the convergence of ST-MetaD simulations is probably larger that the statistical error (the results should be ideally derived by a series of independent simulations 32,77,78 ), the agreement between predicted data from reweighting and the actual simulation result is encouraging. The GAGA native state is stabilized by ∼0.6 kcal/mol in comparison with the control simulation (ΔG°fold = 0.1 ±0.2 kcal/mol), 32 indicating that the NBfix0BPh modification could increase the population of the GAGA native state by ∼24% at 298 K. MD Simulations of Common NA Motifs with the NBfix0BPh Modification Did Not Reveal Any Side Effects. Despite the clear benefits of the NBfix0BPh modification on the structural description of small RNA motifs, one should test its effects on bigger and structured RNA and DNA motifs before any claims are made about its general applicability. Hence, we performed standard MD simulations of systems commonly used for NA ff validation, i.e., RNA and DNA duplexes, RNA SRL and Kt-7 motifs, and two DNA GQs. We compared the results with the control simulation set (i.e., the same ff setup just without NBfix0BPh; see Methods). We observed that the NBfix0BPh modification did not introduce any side effects in the structural description of these motifs (full details are in the Supporting Information). The NBfix0BPh Modification Disfavors Spurious Ladderlike Structures. We showed in a previous section that NBfix0BPh supports canonical A-form backbone conformation in short RNA single strands. In addition, simulations of RNA duplexes revealed that, similarly to the OL3 correction, 79 NBfix0BPh modestly reduces the inclination of the A-form duplex (see the Supporting Information). Therefore, we also tested the possible effect of NBfix0BPh on the potential formation of a spurious deformed duplex structure known as the ladderlike RNA structure, 80 which is associated with the population of high-anti χangle in RNA duplexes. Elimination of the ladderlike structure has been the key improvement introduced by the OL3 parametrization. 7 Thus, we tested NBfix0BPh on a longer 1QC0 duplex 42 (i.e., the canonical decamer with r(GCACCGUUGG) sequence, see Methods) and short r(CGCG)2duplex known for a quick (∼10 ns) transition from the canonical A-form to the spurious ladderlike structure in obsolete AMBER ff versions lacking the OL3 correction. 79 We compared the behavior of two ffs, i.e., the outdated ff 99bsc0 36 and the currently recommended OL3 ff, in combination with two different water models (TIP3P 51 and OPC 26 ) and explicitly probed possible effects of adding the NBfix0BPh modification. As a result, we compared the behavior of eight different simulation setups (Table S4 in Supporting Information) on both duplexes. The initial set of r(CGCG)2 simulations revealed a high tendency to end fraying, i.e., loss of terminal base pairs, which significantly affected outcomes of the simulations (Table S4 in Supporting Information). Therefore, we used the gHBfix21 potential 28 in the next set of r(CGCG)2simulations to stabilize base−base interactions and suppress the extensive base pair fraying. We observed that a spontaneous transition into the ladder structure occurred only in the ff 99bsc0 + TIP3P setup for both short and long duplexes (Table S4 in the Supporting Information). Thus, besides the known reparameterization of χdihedrals of nucleosides in the OL3 ff, 7 there might be other factors affecting the propensity of transitions of canonical duplexes to ladders, namely, (i) the four-point OPC water model and (ii) the NBfix0BPh modification (effective removal of the clash in the A-RNA sugar−phosphate backbone 30 ). Subsequently, we started simulations from the ladder conformation of the r(CGCG)2duplex and observed that only setups with OL3 ff were able to undergo successful transitions back to the canonical A-form structures on the 1 μs time scale; the ladders were corrected in all four OL3 simulations on time scales ∼100−250 ns. NBfix0BPh alone did not eliminate the ladder, at least on the 1 μs time scale (Table S4 in Supporting Information). In summary, the NBfix0BPh modification is a vital addition to the standard AMBER OL3 RNA ff. It can further stabilize canonical A-form duplexes and prevent transitions to artificial ladders. Results from duplex simulations revealed that anti/ high-anti χimbalance of nucleosides (corrected by the crucial OL3 parametrization 7 ) was not the only driving factor for spurious transitions into ladderlike structures. The clash in the backbone remained 30 even after applying OL3 correction, and this clash can be efficiently removed by the NBfix0BPh modification. The OL3 ff can be safely combined with NBfix0BPh, and this combination further reduces the likelihood of sampling various spurious structures. ■CONCLUDING REMARKS We presented the adjustment of pairwise Lennard−Jones parameters of the OL3 AMBER RNA ff for atoms involved in the 0BPh intranucleotide molecular interactions, i.e., the NBfix0BPh modification. This was motivated by the recently identified spurious steric clash in the sugar−phosphate backbone associated with the 0BPh interaction. 30 We Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c00990 J. Chem. Theory Comput. 2023, 19, 8423−8433 8429 optimized NBfix0BPh settings by the reweighting method, which is an efficient approach to assess the effects of ff changes without the necessity to perform actual extensive simulations. Reweighted data indicate that the 0.25 Å decrease of pairwise Ri,j parameters for both −H8···O5′− and −H6···O5′− atom pairs is the optimal adjustment. The NBfix0BPh modification was further verified by an extensive set of new enhanced sampling simulations for a contemporary recommended set of benchmark RNA systems including tetranucleotides, hexanucleotides, and tetraloops. The simulations with the NBfix0BPh modification provided better agreement with the available experimental data. We suggest that the NBfix0BPh modification should be used in all future simulations of these small RNA oligonucleotides because it directly corrects a real imbalance of the AMBER ff nonbonded terms. Elimination of its consequences by other ff terms, such as dihedral potential reparameterizations or gHBfix, is thus not desirable. We also did not find any undesired side effects of the NBfix0BPh modification in standard simulations of larger and more complicated RNA motifs with noncanonical interactions. The data for canonical duplexes indicate that NBfix0BPh may help to prevent A-RNA distortions into the ladderlike structure, although the OL3 still remains essential in preventing this artifact. Hence, the NBfix0BPh modification could be considered as a general addition to the OL3 RNA ff (and likely also to other ff variants using the original AMBER nonbonded parameters). NBfix0BPh eliminates artificial steric clash within nucleotides sampling the anti nucleobase conformation, which improves the description of syn/anti balance. Considering the results shown here, we suggest the utilization of the NBfix0BPh modification on top of the AMBER OL3 ff, 7 modified parameters for phosphates by Steinbrecher et al., 23 OPC water model, 26 and the gHBfix21 potential 28 for MD simulations of RNA especially for small motifs such as short single strands and tetraloops. In summary, NBfix0BPh modification corrects a moderate but evident imbalance in the repulsion term of the AMBER ff nonbonded terms for the 0BPh interaction, which affects the accuracy mainly in case of the simulations of small RNA oligonucleotides. As the NBfix0BPh modification selectively targets only the solute−solute van der Waals term related to the 0BPh interaction, it is robust and does not produce any side effects. Thus, it could be used in combination with any ff from the AMBER family. ■ASSOCIATED CONTENT * sı Supporting Information The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.3c00990. Details about standard MD simulations of various RNA and DNA motifs; supporting tables and figures (PDF) Coordinates of the starting systems (ZIP) ■AUTHOR INFORMATION Corresponding Authors Vojte ch Mlynsky−Institute of Biophysics of the Czech Academy of Sciences, Brno 612 00, Czech Republic; orcid.org/0000-0003-2769-1553; Email: mlynsky@ ibp.cz Pavel Banás−Czech Advanced Technology and Research Institute, CATRIN, Olomouc 779 00, Czech Republic; Institute of Biophysics of the Czech Academy of Sciences, Brno 612 00, Czech Republic; IT4Innovations, VSB−Technical University of Ostrava, Ostrava-Poruba 708 00, Czech Republic; orcid.org/0000-0002-7137-8225; Email: [email protected] Jir í S  poner −Institute of Biophysics of the Czech Academy of Sciences, Brno 612 00, Czech Republic; orcid.org/00000001-6558-6186; Email: [email protected] Authors Petra Kuhrová −Institute of Biophysics of the Czech Academy of Sciences, Brno 612 00, Czech Republic; Czech Advanced Technology and Research Institute, CATRIN, Olomouc 779 00, Czech Republic; orcid.org/0000-0003-1593-5282 Petr Stadlbauer −Institute of Biophysics of the Czech Academy of Sciences, Brno 612 00, Czech Republic; Czech Advanced Technology and Research Institute, CATRIN, Olomouc 779 00, Czech Republic; orcid.org/0000-00025470-8376 Miroslav Krepl −Institute of Biophysics of the Czech Academy of Sciences, Brno 612 00, Czech Republic; Czech Advanced Technology and Research Institute, CATRIN, Olomouc 779 00, Czech Republic; orcid.org/0000-0002-9833-4281 Michal Otyepka −Czech Advanced Technology and Research Institute, CATRIN, Olomouc 779 00, Czech Republic; IT4Innovations, VSB−Technical University of Ostrava, Ostrava-Poruba 708 00, Czech Republic; orcid.org/00000002-1066-5677 Complete contact information is available at: https://pubs.acs.org/10.1021/acs.jctc.3c00990 Notes The authors declare no competing financial interest. ■ACKNOWLEDGMENTS This work was supported by the Czech Science Foundation to V.M., P.S., M.K., and J.S. (grant 23-05639S). This research also received the support of EXA4MIND, a European Union′s Horizon Europe Research and Innovation Programme under grant agreement no. 101092944 (M.O. and P.B.). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the granting authority can be held responsible for them. M.O. and P.B. were also supported by ERDF/ESF project TECHSCALE (CZ.02.01.01/00/22_008/0004587). ■REFERENCES (1) Cheatham, T. E., 3rd; Case, D. A Twenty-five years of nucleic acid simulations. Biopolymers 2013,99 (12), 969−977. (2) Sponer, J.; Banas, P.; Jurecka, P.; Zgarbova, M.; Kuhrova, P.; Havrila, M.; Krepl, M.; Stadlbauer, P.; Otyepka, M. Molecular Dynamics Simulations of Nucleic Acids. From Tetranucleotides to the Ribosome. J. Phys. Chem. Lett. 2014,5(10), 1771−1782. (3) Vangaveti, S.; Ranganathan, S. V.; Chen, A. A. Advances in RNA molecular dynamics: a simulator’s guide to RNA force fields. Wiley Interdiscip. Rev.: RNA 2017,8(2), 1396 DOI: 10.1002/wrna.1396. (4) Smith, L. G.; Zhao, J.; Mathews, D. H.; Turner, D. H. Physicsbased all-atom modeling of RNA energetics and structure. Wiley Interdiscip. Rev.: RNA 2017,8(5), No. e1422, DOI: 10.1002/ wrna.1422. (5) Sponer, J.; Bussi, G.; Krepl, M.; Banas, P.; Bottaro, S.; Cunha, R. A.; Gil-Ley, A.; Pinamonti, G.; Poblete, S.; Jurecka, P.; et al. RNA Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c00990 J. Chem. Theory Comput. 2023, 19, 8423−8433 8430 Structural Dynamics As Captured by Molecular Simulations: A Comprehensive Overview. Chem. Rev. 2018,118 (8), 4177−4338. (6) Denning, E. J.; Priyakumar, U. D.; Nilsson, L.; Mackerell, A. D. Impact of 2’-Hydroxyl Sampling on the Conformational Properties of RNA: Update of the CHARMM All-Atom Additive Force Field for RNA. J. Comput. Chem. 2011,32 (9), 1929−1943. (7) Zgarbova, M.; Otyepka, M.; Sponer, J.; Mladek, A.; Banas, P.; Cheatham, T. E.; Jurecka, P. Refinement of the Cornell et al. Nucleic Acids Force Field Based on Reference Quantum Chemical Calculations of Glycosidic Torsion Profiles. J. Chem. Theory Comput 2011,7(9), 2886−2902. (8) Chen, A. A.; Garcia, A. E. High-resolution reversible folding of hyperstable RNA tetraloops using molecular dynamics simulations. Proc. Natl. Acad. Sci. U. S. A. 2013,110 (42), 16820−16825. (9) Yang, C.; Lim, M.; Kim, E.; Pak, Y. Predicting RNA Structures via a Simple van der Waals Correction to an All-Atom Force Field. J. Chem. Theory Comput 2017,13 (2), 395−399. (10) Aytenfisu, A. H.; Spasic, A.; Grossfield, A.; Stern, H. A.; Mathews, D. H. Revised RNA Dihedral Parameters for the Amber Force Field Improve RNA Molecular Dynamics. J. Chem. Theory Comput 2017,13 (2), 900−915. (11) Tan, D.; Piana, S.; Dirks, R. M.; Shaw, D. E. RNA force field with accuracy comparable to state-of-the-art protein force fields. Proc. Natl. Acad. Sci. U. S. A. 2018,115 (7), E1346−E1355. (12) Bergonzo, C.; Henriksen, N. M.; Roe, D. R.; Cheatham, T. E., 3rd Highly sampled tetranucleotide and tetraloop motifs enable evaluation of common RNA force fields. RNA 2015,21 (9), 1578− 1590. (13) Kuhrova, P.; Best, R. B.; Bottaro, S.; Bussi, G.; Sponer, J.; Otyepka, M.; Banas, P. Computer Folding of RNA Tetraloops: Identification of Key Force Field Deficiencies. J. Chem. Theory Comput 2016,12 (9), 4534−4548. (14) Kuhrova, P.; Mlynsky, V.; Zgarbova, M.; Krepl, M.; Bussi, G.; Best, R. B.; Otyepka, M.; Sponer, J.; Banas, P. Improving the Performance of the Amber RNA Force Field by Tuning the Hydrogen-Bonding Interactions. J. Chem. Theory Comput 2019,15 (5), 3288−3305. (15) Kuhrova, P.; Mlynsky, V.; Zgarbova, M.; Krepl, M.; Bussi, G.; Best, R. B.; Otyepka, M.; Sponer, J.; Banas, P. Correction to ″Improving the Performance of the Amber RNA Force Field by Tuning the Hydrogen-Bonding Interactions″.J. Chem. Theory Comput 2020,16 (1), 818−819. (16) Zhao, J.; Kennedy, S. D.; Turner, D. H. Nuclear Magnetic Resonance Spectra and AMBER OL3 and ROC-RNA Simulations of UCUCGU Reveal Force Field Strengths and Weaknesses for SingleStranded RNA. J. Chem. Theory Comput 2022,18 (2), 1241−1254. (17) Krepl, M.; Pokorna, P.; Mlynsky, V.; Stadlbauer, P.; Sponer, J. Spontaneous binding of single-stranded RNAs to RRM proteins visualized by unbiased atomistic simulations with a rescaled RNA force field. Nucleic Acids Res. 2022,50 (21), 12480−12496. (18) Kuhrova, P.; Mlynsky, V.; Otyepka, M.; Sponer, J.; Banas, P. Sensitivity of the RNA Structure to Ion Conditions as Probed by Molecular Dynamics Simulations of Common Canonical RNA Duplexes. J. Chem. Inf Model 2023,63 (7), 2133−2146. (19) Cesari, A.; Gil-Ley, A.; Bussi, G. Combining Simulations and Solution Experiments as a Paradigm for RNA Force Field Refinement. J. Chem. Theory Comput 2016,12 (12), 6192−6200. (20) Cesari, A.; Bottaro, S.; Lindorff-Larsen, K.; Banas, P.; Sponer, J.; Bussi, G. Fitting Corrections to an RNA Force Field Using Experimental Data. J. Chem. Theory Comput 2019,15 (6), 3425− 3431. (21) Zhao, J.; Kennedy, S. D.; Berger, K. D.; Turner, D. H. Nuclear Magnetic Resonance of Single-Stranded RNAs and DNAs of CAAU and UCAAUC as Benchmarks for Molecular Dynamics Simulations. J. Chem. Theory Comput 2020,16 (3), 1968−1984. (22) Case, D. A.; Aktulga, H. M.; Belfon, K.; Ben-Shalom, I. Y.; Berryman, J. T.; Brozell, S. R.; Cerutti, D. S. T.E.; Cheatham, I.; Cisneros, G. A.; Cruzeiro, V. W. D.; et al. AMBER 2023; University of California, San Francisco;2023. (23) Steinbrecher, T.; Latzer, J.; Case, D. A. Revised AMBER Parameters for Bioorganic Phosphates. J. Chem. Theory Comput 2012, 8(11), 4405−4412. (24) Mlynsky, V.; Kuhrova, P.; Zgarbova, M.; Jurecka, P.; Walter, N. G.; Otyepka, M.; Sponer, J.; Banas, P. Reactive conformation of the active site in the hairpin ribozyme achieved by molecular dynamics simulations with epsilon/zeta force field reparametrizations. J. Phys. Chem. B 2015,119 (11), 4220−4229. (25) Bergonzo, C.; Cheatham, T. E., 3rd. Improved Force Field Parameters Lead to a Better Description of RNA Structure. J. Chem. Theory Comput 2015,11 (9), 3969−3972. (26) Izadi, S.; Anandakrishnan, R.; Onufriev, A. V. Building Water Models: A Different Approach. J. Phys. Chem. Lett. 2014,5(21), 3863−3871. (27) Bottaro, S.; Bussi, G.; Kennedy, S. D.; Turner, D. H.; LindorffLarsen, K. Conformational ensembles of RNA oligonucleotides from integrating NMR and molecular simulations. Sci. Adv. 2018,4(5), No. eaar8521. (28) Frohlking, T.; Mlynsky, V.; Janecek, M.; Kuhrova, P.; Krepl, M.; Banas, P.; Sponer, J.; Bussi, G. Automatic Learning of HydrogenBond Fixes in the AMBER RNA Force Field. J. Chem. Theory Comput 2022,18 (7), 4490−4502. (29) Mlynsky, V.; Kuhrova, P.; Kuhr, T.; Otyepka, M.; Bussi, G.; Banas, P.; Sponer, J. Fine-Tuning of the AMBER RNA Force Field with a New Term Adjusting Interactions of Terminal Nucleotides. J. Chem. Theory Comput 2020,16 (6), 3936−3946. (30) Mrazikova, K.; Mlynsky, V.; Kuhrova, P.; Pokorna, P.; Kruse, H.; Krepl, M.; Otyepka, M.; Banas, P.; Sponer, J. UUCG RNA Tetraloop as a Formidable Force-Field Challenge for MD Simulations. J. Chem. Theory Comput 2020,16 (12), 7601−7617. (31) Zirbel, C. L.; Sponer, J. E.; Sponer, J.; Stombaugh, J.; Leontis, N. B. Classification and energetics of the base-phosphate interactions in RNA. Nucleic Acids Res. 2009,37 (15), 4898−4918. (32) Mlynsky, V.; Janecek, M.; Kuhrova, P.; Frohlking, T.; Otyepka, M.; Bussi, G.; Banas, P.; Sponer, J. Toward Convergence in Folding Simulations of RNA Tetraloops: Comparison of Enhanced Sampling Techniques and Effects of Force Field Modifications. J. Chem. Theory Comput 2022,18 (4), 2642−2656. (33) Kofinger, J.; Rozycki, B.; Hummer, G. Inferring Structural Ensembles of Flexible and Dynamic Macromolecules Using Bayesian, Maximum Entropy, and Minimal-Ensemble Refinement Methods. Methods Mol. Biol. 2019,2022, 341−352. (34) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules (vol 117, pg 5179, 1995). J. Am. Chem. Soc. 1996,118 (9), 2309−2309. (35) Wang, J. M.; Cieplak, P.; Kollman, P. A. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comput. Chem. 2000,21 (12), 1049−1074. (36) Perez, A.; Marchan, I.; Svozil, D.; Sponer, J.; Cheatham, T. E.; Laughton, C. A.; Orozco, M. Refinement of the AMBER Force Field for Nucleic Acids: Improving the Description of Alpha/Gamma Conformers. Biophys. J. 2007,92 (11), 3817−3829. (37) Zgarbova, M.; Sponer, J.; Otyepka, M.; Cheatham, T. E.; Galindo-Murillo, R.; Jurecka, P. Refinement of the Sugar−Phosphate Backbone Torsion Beta for AMBER Force Fields Improves the Description of Zand B-DNA. J. Chem. Theory Comput 2015,11 (12), 5723−5736. (38) Yoo, J.; Aksimentiev, A. New tricks for old dogs: improving the accuracy of biomolecular force fields by pair-specific corrections to non-bonded interactions. Phys. Chem. Chem. Phys. 2018,20 (13), 8432−8449. (39) Nozinovic, S.; Furtig, B.; Jonker, H. R.; Richter, C.; Schwalbe, H. High-resolution NMR structure of an RNA model system: the 14mer cUUCGg tetraloop hairpin RNA. Nucleic Acids Res. 2010,38 (2), 683−694. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c00990 J. Chem. Theory Comput. 2023, 19, 8423−8433 8431