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.


■ 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 f fs.
−18 Despite the suboptimal performance for certain motifs, mainly the short RNA single strands, 14,16,19−21 the AMBER OL3 f f from 2010 7 still remains the state-of-the-art for RNA simulations and safest option to start with. 22It has been shown in the past decade that the OL3 performance can be partially tweaked by optimizing simulation settings and by adding new f f terms that are orthogonal to the current f f 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. 29n 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 G S+1 (C8−H8)•••G S+1 (O5′) interaction resulted in the spurious flipping of the G S+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 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). 31The attempt to adjust the 0BPh interaction in the structural context of UUCG TL via nonbonded fix (NBfix) stabilized the native G S+1 phosphate conformation 30 and indicated that it could be beneficial for the stability of native states of both GAGA and UUCG TLs. 32n 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 (NBfix 0BPh modification) as a general improvement of the AMBER OL3 RNA f f (and likely other f f variants using the original AMBER nonbonded parameters).We first used the reweighting approach to optimize settings for NBfix 0BPh modification.Reweighting methods are one of the strongest tools of f f 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 NBfix 0BPh modification (with a similar setting as those for RNA nucleotides) on DNA duplex and guanine quadruplexes (GQs).We show that the NBfix 0BPh 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 NBfix 0BPh 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.

Details about AMBER RNA and DNA f fs and Their
Adjustments.We used standard AMBER OL3 (known also as χOL3) 7,34−36 and AMBER OL15 37 f fs for RNA and DNA simulations, respectively.OL3 RNA f f 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,24This RNA ff version is abbreviated as OL3 CP 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 OL3 CP f f 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. 14We 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. 29For some specific tests, we also used standard OL3 (i.e., without phosphate modification) and older f f 99bsc0 36 RNA f fs.
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 OL3 CP and DNA OL15 f fs.Namely, we reduced the minimum-energy distance of the Lennard−Jones potential (i.e., R i,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 NBfix 0BPh ; 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 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.
in the RNA OL3 CP , whereas the OS atom type involves O5′, O3′, and also O4' atoms from deoxyribose sugar in the DNA OL15 f f.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 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,60REST2 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). 61Further details about REST2 settings can be found elsewhere. 14One 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. 16The 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,63T-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 R i,j values of the Lennard−Jones potential for the NBfix 0BPh 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.−70 We analyzed separately four NMR observables, i.e., (i) backbone 3 J scalar couplings, (ii) sugar 3 J 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. 3 J scalar couplings were calculated via Karplus relationships; NOEs and uNOEs were obtained as averages over the N samples, i.e., u ( )NOE ( ) ; 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 The combination of all those analyzed NMR observables (calculated as weighted arithmetic mean) provided the total χ 2 value for each REST2 simulation: , where χ i 2 and n i are individual χ 2 values and the number of observables, respectively, corresponding to the particular NMR observables.

Journal of Chemical Theory and Computation
The lower the total χ 2 value 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°f old ). 28,32The reference native structure of GAGA TL was taken from our previous work. 14The ε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°f old estimations and convergence).

■ RESULTS AND DISCUSSION
Reweighting Reveals Optimal Settings for the NBfix 0BPh Modification.The reweighting algorithm 66 allows one to efficiently re-evaluate results of MD simulations under the assumption of a modified f f parametrization without the necessity to perform new simulations. 27,32,33Here, we reweighted the NBfix 0BPh settings for a set of enhanced sampling simulations of RNA TNs and TLs.More specifically, we used results from OL3 CP REST2 simulations with gHBfix19 and tHBfix20 potentials of r(AAAA), r(CCCC), and r(CAAU) TNs 29 and from folding OL3 CP ST-MetaD simulations with gHBfix19 of GAGA and UUCG TLs. 32We monitored changes in total χ 2 values (and its components) and ΔG°f old energies for TNs and TLs, respectively.Each NBfix 0BPh setting was defined by scanning R i,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.48intervals for purine and pyrimidine nucleotides, respectively, using 0.025 Å steps.
Results from TN simulations indicate that a 0.25 Å decrease of R i,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 R i,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 R i,j values for the −H8•••O5′− atom pairs of purine nucleotides is the dominant contributor for shifting ΔG°f old 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 R i,j values for both purines and pyrimidines simultaneously (Figure 2A).
Modifications of R i,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 R i,j value for only pyrimidine residues indicates stabilization; ΔG°f old is decreasing from ∼5.4 kcal/mol (default R i,j value) up to ∼3.5 kcal/mol (final and lowest tested R i,j ) (Figure 2B).In contrast, decreased R i,j values for sole purines indicate destabilization (ΔG°f old is increasing from ∼5.4 kcal/mol up to ∼6.0 kcal/mol, Figure 2B).Such a counterintuitive effect of the decreased R i,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 (G L4 ) is sampling noncanonical anti orientation of the χ dihedral (syn states are required for the native conformation, and those are ).The reweighting data are still suggesting that the UUCG native state is significantly destabilized by the ff even after the application of the NBfix 0BPh modification.This is not surprising because multiple f f 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. 28ot affected by the NBfix 0BPh modification; see Figure 1 and ref 32 for further discussion).Most importantly, simultaneous modification of R i,j values for both purines and pyrimidines, i.e., the NBfix 0BPh 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 R i,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°f old energy with decreasing R i,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 R i,j parameters that would be enough to relax the steric clash of −H8•••O5′− and −H6•••O5′− atom pairs observed in the original OL3 CP f f.We found that a 0.25 Å decrease of the particular R i,j parameter for both purines and pyrimidines appears to be the optimal choice.NBfix 0BPh Modification of the AMBER OL3 CP RNA f f Improves the Structural Description of Single-Stranded RNA Motifs.We used the optimized setting for the NBfix 0BPh 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,29We used OL3 CP f f 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 2 for TNs and HNs, respectively.
The NBfix 0BPh modification further fine-tunes especially the r(CAAU) and r(AAAA) ensembles, which now agree excellently with NMR data (total χ 2 value below 1; Table 1).The NBfix 0BPh modification indirectly stabilizes the anti state of A residues and thus resolves the previously reported excessive sampling of syn states. 29,74We also observed that NBfix 0BPh increased the population of canonical A-form like states in Simulations used a combination of OL3 CP RNA f f 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 χ 2 values were obtained by comparing calculated and experimental backbone 3 J scalar couplings, sugar 3 J 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−70 Application of NBfix 0BPh 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 χ 2 values; 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,75Hence, further adjustment of AMBER RNA f fs 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. 17It 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. 17STAfix, 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 NBfix 0BPh 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,76On 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. 16r(UCAAUC) REST2 simulations performed here show that the NBfix 0BPh 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 χ 2 values 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 NBfix 0BPh 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 NBfix 0BPh 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 OL3 CP 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. 16We observed that calculated χ 2 values are almost an order of magnitude higher than those from TN simulations and the other r(UCAAUC) HN (Tables 1 and 2).Inspection of individual χ 2 components (see Methods) revealed that the largest deviations between experimentally measured and predicted results from simulations come from sugar 3 J scalar couplings.Experimental data suggest that each residue is sampling a mixture of C3′-endo and C2′-endo sugar puckers. 16Indeed, 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. 16However, the formation of conformers with the G5C2 base pair is not supported. 16n 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 OL3 CP f f 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 NBfix 0BPh Modification.We performed ST-MetaD simulation of the r(gcGAGAgc) TL to explicitly test the effect of the NBfix 0BPh modification on the folding of this small 8-mer RNA TL motif.We obtained a ΔG°f old 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°f old = 0.1 ± 0.2 kcal/mol), 32 indicating that the NBfix 0BPh modification could increase the population of the GAGA native state by ∼24% at 298 K.
MD Simulations of Common NA Motifs with the NBfix 0BPh Modification Did Not Reveal Any Side Effects.Despite the clear benefits of the NBfix 0BPh 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 f f setup just without NBfix 0BPh ; see Methods).We observed that the NBfix 0BPh modification did not introduce any side effects in the structural description of these motifs (full details are in the Supporting Information).
The NBfix 0BPh Modification Disfavors Spurious Ladderlike Structures.We showed in a previous section that NBfix 0BPh 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 NBfix 0BPh modestly reduces the inclination of the A-form duplex (see the Supporting Information).Therefore, we also tested the possible effect of NBfix 0BPh 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. 7Thus, we tested NBfix 0BPh on a longer 1QC0 duplex 42 (i.e., the canonical decamer with r(GCACCGUUGG) sequence, see Methods) and short r(CGCG) 2 duplex known for a quick (∼10 ns) transition from the canonical A-form to the spurious ladderlike structure in obsolete AMBER f f versions lacking the OL3 correction. 79We compared the behavior of two f fs, i.e., the outdated f f 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 NBfix 0BPh 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) 2 simulations 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 f f 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 f f, 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 NBfix 0BPh 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) 2 duplex 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.NBfix 0BPh alone did not eliminate the ladder, at least on the 1 μs time scale (Table S4 in Supporting Information).
In summary, the NBfix 0BPh 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 NBfix 0BPh modification.The OL3 f f can be safely combined with NBfix 0BPh , 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 f f for atoms involved in the 0BPh intranucleotide molecular interactions, i.e., the NBfix 0BPh modification.This was motivated by the recently identified spurious steric clash in the sugar−phosphate backbone associated with the 0BPh interaction. 30We

Journal of Chemical Theory and Computation
optimized NBfix 0BPh settings by the reweighting method, which is an efficient approach to assess the effects of f f changes without the necessity to perform actual extensive simulations.Reweighted data indicate that the 0.25 Å decrease of pairwise R i,j parameters for both −H8•••O5′− and −H6•••O5′− atom pairs is the optimal adjustment.The NBfix 0BPh 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 NBfix 0BPh modification provided better agreement with the available experimental data.We suggest that the NBfix 0BPh modification should be used in all future simulations of these small RNA oligonucleotides because it directly corrects a real imbalance of the AMBER f f nonbonded terms.Elimination of its consequences by other f f terms, such as dihedral potential reparameterizations or gHBfix, is thus not desirable.We also did not find any undesired side effects of the NBfix 0BPh modification in standard simulations of larger and more complicated RNA motifs with noncanonical interactions.The data for canonical duplexes indicate that NBfix 0BPh may help to prevent A-RNA distortions into the ladderlike structure, although the OL3 still remains essential in preventing this artifact.
Hence, the NBfix 0BPh modification could be considered as a general addition to the OL3 RNA f f (and likely also to other f f variants using the original AMBER nonbonded parameters).NBfix 0BPh 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 NBfix 0BPh modification on top of the AMBER OL3 f f, 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, NBfix 0BPh modification corrects a moderate but evident imbalance in the repulsion term of the AMBER f f nonbonded terms for the 0BPh interaction, which affects the accuracy mainly in case of the simulations of small RNA oligonucleotides.As the NBfix 0BPh 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 f f from the AMBER family.

■ ASSOCIATED CONTENT
* sı Supporting Information

Figure 2 .
Figure 2. Effects of different setting of NBfix 0BPh modifications on the stability of the native state of RNA TLs.Folding free energies (ΔG°f old ) were estimated by reweighting of GAGA (A) and UUCG (B) trajectories from OL3 CP ST-MetaD simulations with the gHBfix19 potential.Plots show the dependence of ΔG°f old energies on modified vdW parameters for purines and pyrimidines, i.e., changes of R i,j values for −H8•••O5′− (cyan line) and −H6•••O5′− (magenta line) atom pairs from the original (standard AMBER f f) values.The black line presents the outcome, where R i,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 NBfix 0BPh modification (0.25 Å decrease of pairwise terms; green vertical line) stabilizes native states of both TLs by providing lower ΔG°f old energies, which are in better agreement with experiments (reported ΔG°f old 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 NBfix 0BPh modification.This is not surprising because multiple f f 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

Table 1 .
REST2 Simulations of RNA TNs and Their Comparison with Experiments a

Table 2 .
REST2 Simulations of RNA HNs and Their Comparison with Experiments aSimulations used a combination of OL3 CP RNA f f 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.values were obtained by comparing calculated and experimental backbone 3 J scalar couplings, sugar 3 J 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).
a b χ 2