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

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


■ INTRODUCTION
The gene expression process inevitably involves formation of mixed RNA and DNA duplexes (hybrids) during the transcription, with the newly synthesized RNA strand temporarily base paired to its DNA template. 1The opposite process occurs during the reverse transcription, when the DNA is synthesized based on the RNA template. 2The combination of ubiquity and relatively short lifetimes makes hybrids an important but less studied component of cellular biology.−7 Structurally, the hybrids form classical right-handed double helices, but the major structural biology methods somewhat disagree on the details of their geometries.−20 Overall, there is a notable degree of ambiguity in the experimental structural characterization of DNA/RNA hybrids, a situation which could be greatly improved by well-executed computational studies.
When performing molecular dynamics (MD) simulations of the hybrids, the choice of the nucleic acid force fields (f fs) is not trivial.−23 This leads to a question of which combination of nucleic acid f fs should be utilized for hybrids to obtain optimal performance consistent with the experiments.In some cases, both RNA and DNA f f parameters are developed by a single research group, 22,24−26 making their combined use logical.−35 In those cases, the choice of the second nucleic acids f f is arbitrary and difficult to justify, especially when no specific recommendation is given by the f f developers and the hybrids are seldom among the systems used for benchmarking.−48 Here, we present a benchmark study that evaluates the performance of several modern AMBER nucleic acid f fs for MD simulations of DNA/RNA hybrids.We also test the CHARMM36 49,50 and the latest Drude 51−53 and Amoeba 54,55 polarizable ffs.To observe the helical form changes within the same sequence context, we utilize the structure of the Dickerson−Drew B-DNA dodecamer 56 and its modeled A-RNA and hybrid sequence equivalents.As an experimental reference for the hybrid structure, we use the polypurine tract (PPT) hybrid duplex, 13 allowing comparison of the base pair and helical parameters predicted by the f fs with experimental X-ray crystallography values.We find that none of the tested f fs are able to fully reproduce the hybrid structure, although the performance of some f fs may be sufficient for many specific applications.A significant challenge arises from the inability of many f fs to sample the experimental C3′-endo pucker of the DNA strand, resulting in shifts of all major helical parameters toward more "B-form-like" geometries, even in sequences where experimental data clearly indicate dominance of the Aform.On the other hand, the simulations correctly reproduce the anisotropic buckle of the hybrids.The data presented herein offer insights into the MD simulation description of the hybrids, complement existing experimental data and outline potential opportunities for further f f adjustments.Furthermore, the results can guide the selection of utilizable f f combinations for future MD studies of DNA/RNA hybrids and their complexes with proteins.

■ METHODS
Selection of Initial Structures.We have used the X-ray structure of the B-DNA dodecamer (Dickerson−Drew dodecamer (DD); PDB: 1BNA) 56 (Figure 1) as a starting structure for simulations of the B-DNA duplex (DD_DNA).A-RNA (DD_RNA) and hybrid (DD_hybrid) duplexes of identical sequence were prepared with the Nucleic Acid Builder. 57Thymines were replaced with uracils for the RNA strands and the initial helical geometry of the modeled duplexes corresponded to the A-form.We also prepared hybrid duplex structure in B-form to verify the convergence of helical forms in simulations.In addition to the DD structures, simulations were also performed using the X-ray structure of the polypurine tract hybrid duplex (PPT) (Figure 1). 13Lastly, to explore performance of a sequence radically different from PPT, some simulations were also performed with modeled hybrid duplex possessing 0% deoxypyrimidine content (low_-dPyr), with randomized RNA strand sequence of 5′-CCUCUCUCUCCC-3′.
System Building and Simulation Protocol − the AMBER and CHARMM36 f fs.From the AMBER family of f fs, we tested the OL15, 58 OL21, 24 bsc1, 28 DES-Amber, 26 and Tumuc1 29 f fs for DNA and OL3, 22 ROC-RNA, 27 and DES-Amber 26 f fs for RNA.The tLeap program of AMBER 22 59 was used to generate initial files for all ffs except DES-Amber (see below) and simulations were subsequently performed in AMBER22 using the pmemd.MPI and pmemd.cuda 60programs for equilibration 61 and production simulations, respectively.The duplexes were placed in an octahedral box of SPC/E 62 water molecules with minimal distance of 12 Å between the solute and the box border.KCl ions parametrized by Joung and Cheatham 63 were added to neutralize the systems and to establish excess-salt ion concentration of 0.15 M. The production simulations were run for 3 μs and three independent trajectories were obtained for each system and ff combination.We used SHAKE 64 and hydrogen mass repartitioning 65 to allow a 4 fs integration step.Long-range electrostatics were treated with particle mesh Ewald 66 and periodic boundary conditions were applied.The distance cutoff for Lennard-Jones interactions was set to 9 Å.The production simulations were performed in a constant pressure ensemble with pressure and temperature being regulated with Monte Carlo barostat and Langevin thermostat, 67 respectively.In all simulations, we have applied 1 kcal/mol structure-specific HBfix (sHBfix) potential 68 between donors and acceptors to stabilize the Watson−Crick H-bonds in the terminal base pairs, thus essentially eliminating the terminal base pair fraying.Although a minor population of frayed terminal base pairs is to be realistically expected, its estimated thermodynamic convergence in MD simulations is far beyond 10 μs. 69This would lead to major sampling uncertainty and complicate the f f comparison.The problem has been completely circumvented by increasing stability of the terminal base pairs.
As the latest DES-Amber parameters for nucleic acids were unavailable for AMBER at the time, we have used Gromacs2018 70 to perform all DES-Amber 26 simulations, utilizing the parameters published in ref 26 and provided by the authors on our request.The recommended TIP4P-D water model 71 along with CHARMM22 ions 71 was utilized in all DES-Amber simulations.Note that the Gromacs DES-Amber parameters 26 for DNA and RNA as provided by the authors contain overlapping atomic types between the two biopolymers and therefore do not natively allow hybrid simulations.To circumvent this, we modified all the DNA atomic types and the associated parameters by adding the "D" prefix to globally distinguish them from the RNA-specific parameters.
We have likewise used Gromacs2018 to perform all the CHARMM36 49,50 simulations, utilizing the CHARMMmodified TIP3P water molecules and CHARMM36 KCl parameters for ions.Due to the different molecular dynamics engine, the simulation protocol of the DES-Amber and CHARMM36 simulations somewhat differed from the ones used with AMBER22.Namely, the simulations utilized PLUMED 2.5 72 to implement the sHBfix for the terminal base pairs (see above). 73The simulations were also performed in a rhombic dodecahedral box.Bonds involving hydrogens were constrained using the LINCS algorithm. 74The cutoff distance for the direct space summation of the electrostatic interactions was 12 Å and the simulations were performed at 298 K using the stochastic velocity rescale thermostat. 75The other settings, equilibration protocols, and system building choices were the same as in the AMBER22 simulations.
System Building and Simulation Protocol − DRUDE and AMOEBA Polarizable f fs.For the simulations with CHARMM Drude19 [51][52][53]76 and AMOEBA18 54,55 polarizable f fs, the systems were prebuilt in AMBER22 followed by minimization and equilibration using the OL21/OL3 f fs.In the case of Drude, the resulting structures were transformed into the Drude polarizable model using the CHARMM software (version 44b1). 77 During te conversion process, Drude particles were introduced for all heavy atoms and the lone pairs associated with each hydrogen acceptor.The SPC/E water molecules were converted into the polarizable SWM4-NDP model.78 Each system was then subjected to minimization and equilibration using the NAMD 2.13 package.79 First, we performed 1000-step minimization of the waters and ions while keeping the nucleic acid atoms restrained by a harmonic 500 kcal/mol/Å 2 positional restraint.Keeping the duplexes restrained, this was followed by 500 ps-long MD simulation under constant pressure conditions, to relax the overall density of the systems.Afterward, the duplexes were relaxed through several minimization cycles, gradually reducing the strength of the positional restraints placed on the nucleic acid atoms in a manner outlined in ref 80. Finally, the systems were heated under constant volume conditions for 100 ps and then equilibrated again under constant pressure conditions for additional 100 ps.PLUMED 2.7 72 was used to implement the sHBfix for the terminal base pairs.The obtained structures were then utilized for three separate 1 μs-long production simulations at 298 K in OpenMM 8.0.81 Drude Langevin integrator 82,83 was utilized with a time step of 1 fs.The pressure was maintained at 1 bar utilizing the Monte Carlo barostat.84 The covalent bonds involving hydrogens were kept rigid using the SHAKE 64 and SETTLE 85 algorithms for solute and waters, respectively.A constraint of 0.2 Å was applied to limit the length of Drude-nuclei bonds.Electrostatic interactions were treated using the particle mesh Ewald method (PME) 86 with a 12 Å cutoff for the real space term.Nonbonded interactions were truncated at 12 Å, using a switching function from 10 to 12 Å.
The AMOEBA simulations were performed in GPUaccelerated Tinker version. 87The prebuilt and equilibrated initial structures (see above) were converted into xyz coordinates and minimized in 10 000 steps using the steepest descent method.This was followed by heating the systems to 300 K and equilibrating the pressure to 1 bar before the production simulations were run for 1 μs in constant volume ensemble.The RESPA integrator was used with the integration step of 1 and 2 fs for minimization/equilibration and production simulations, respectively.Stochastic velocity rescale thermostat 75 and Monte Carlo barostat were used to maintain temperature and pressure, respectively.The applied real-space cutoff for electrostatics and van der Waals was 7 and 12 Å, respectively.All other settings were set to default.As Tinker did not at the time allow use of sHBfix or patching with Plumed code, we instead stabilized the H-bonds of the terminal base pairs by applying a linearly growing penalty function, with a force constant of 10 kcal/mol/Å for hydrogenacceptor distances greater than 2.5 Å.
Analyses.Cpptraj 88 and VMD 89 were used to analyze and visualize the trajectories.Molecular figures and graphs were prepared with Povray and gnuplot, respectively.−92 We used the Altona et al. definition to calculate the sugar pucker. 93The terminal base pairs as well as the base pair steps between the terminal and subterminal base pairs were not included in the analyses.For the DD structures, we evaluated the base pair parameters as separate averages of the AT(U) and GC base pairs.For the other sequences, averages of all nonterminal base pairs were calculated.For the CHARMM36 simulations, simulation frames with disrupted nonterminal base pairs were excluded from the analyses (see below).The histogram analyses were performed in cpptraj with bin sizes ranging from 0.1 to 5 depending on the range of absolute values observed for the individual parameters (see Supporting Information).All histograms were normalized by cpptraj so that the sum of the bins equals one and the presented graphs are shown with cubic spline smoothing applied by Gnuplot.Unless specified otherwise, the presented The DNA nucleotide puckers were restrained to C3′-endo region (pseudorotation values of −10°to 40°) by five flat-well dihedral restraint potential functions applied on each 2-deoxyribose (see the Supporting Information).b The Li and Merz parameters 94 for KCl ions were utilized.
The initial structure of the hybrid corresponded to the B-form helix (see Methods).d The χ dihedral potentials of the OL15 DNA f f were replaced with the RNA χ dihedral potentials from OL3 f f.See ref 95 for more details.e 2 kcal/mol stabilizing sHBfix potential was applied on every base pairing H-bond.
values correspond to combined simulation ensembles of the three simulations performed for each system and f f combination.The first 500 ns of each production simulation were considered an extended equilibration period and were not included in the analyses.All the trajectories were visually inspected.
Convergence of the Helical Forms.We have executed at least three independent simulations for each system and f f combination (Table 1).Comparison of the independent simulations revealed average RMSD between the ensembles of less than 0.1 Å, suggesting excellent convergence.In fact, we universally observed the systems approaching convergence of their helical parameters in ∼200−300 ns, with B-DNA and hybrids taking longer than the A-RNA due to sampling of the BI/BII populations in the DNA strand(s).Convergence of the hybrid's helical form was of special concern as in principle it can be initially modeled as either the A-form (utilized in most of our simulations; see above) or B-form, with the initial The DNA nucleotide puckers were restrained to C3′-endo region (pseudorotation values of −10°to 40°) by five flat-well dihedral restraint potential functions applied on each 2-deoxyribose (see the Supporting Information).b The Li and Merz parameters 94 for KCl ions were utilized.

c
The initial structure of the hybrid corresponded to the B-form helix (see Methods).d The χ dihedral potentials of the OL15 DNA f f were replaced with the RNA χ dihedral potentials from OL3 f f.See ref 95 for more details.e 2 kcal/mol stabilizing sHBfix potential was applied on every base pairing H-bond.f Average values observed in the PPT experimental structure. 13See Supporting Information Table S1 for values observed for other experimental structures of the hybrid duplexes.structure potentially affecting the results.However, by comparing the average structures obtained from simulations started from the two helical forms, we observed the RMSD of the second halves of the simulation ensembles of ∼0.1 Å.The resulting helical parameters were also very similar (Table 2), showing differences on par with the negligible differences observed among parallel simulations started from the same helical form.This confirms that the helical form of the starting structure is not affecting the results.We have likewise monitored the time development of all the individual helical parameters, dihedrals and sugar puckers described in the text, observing extensive and frequent fluctuations throughout the simulations with all f fs.With the exception of the terminal base pairs (which were not included in the analyses), we observed largely uniform distribution of the parameters along the DNA and RNA strands.In other words, there was no significant tendency for individual nucleotides to populate some specific states, showing only minor deviations from the average values calculated for the entire strand.In conclusion, we suggest that our simulations were long enough to capture all the qualitative trends associated with the individual systems and f fs as well as to achieve a very good quantitative convergence.

■ RESULTS AND DISCUSSION
Below, we present an extensive set of MD simulations involving A-RNA, B-DNA, and hybrid DNA/RNA duplexes (126 independent simulations in total with a cumulative time of 330 μs, see Table 1).Average simulation values of the helical parameters relevant for the study are given in Table 2 and Supporting Information Table S2.To summarize, by simulating A-RNA, B-DNA, and hybrid helices of the same sequence (Figure 1 and Table 1), we were able to capture the changes in helical form balance among the three duplexes (Table 2).The details differed significantly depending on the utilized f fs, but in general, the individual helical parameters of the hybrid duplexes were either closer but not identical to the A-form values (inclination, helical twist, propeller twist, minor groove width), near midway between values characteristic for A-and B-forms (x-displacement, roll, slide, helical rise), or they possessed unique values divergent from both forms (buckle, tilt, tip) (Table 2 and Supporting Information Table S2).By comparing our simulations with the helical values observed experimentally for the PPT, we then evaluated how accurately the individual f f combinations reproduce this hybrid structure (Table 2 and Supporting Information Table S2).The following text details the key observations derived from our analyses.
AMBER Family of f fs Did Not Reproduce the C3′-Endo Pucker of the DNA Strand of the Hybrid Duplexes.In all hybrid simulations using the nonpolarizable AMBER f fs, we observed the DNA nucleotides possessing exclusively C2′-endo puckers (∼135−180°) (Figure 2) along with N-glycosidic χ dihedrals in the high-anti region (∼250°), both characteristic features of the B-form duplexes.This simulation result stands in stark contrast to experiments that instead indicate C3′-endo pucker (∼−10−40°) and χ anti (∼200°), exclusively in the case of X-ray structures and as a population mixture of C2′/C3′-endo puckers in NMR data.The simulations showed slightly lower pucker values for the DNA in hybrids than in pure DNA, but still well within the C2′-endo region.We suggest these observations are related to the known difficulties with addressing the A-DNA/B-DNA balance and transitions in MD simulations using the AMBER f fs, where the absence of the C3′-endo pucker for DNA is also problematic. 95,96ccording to the solution NMR measurements, the sugar pucker of the DNA nucleotides in hybrids exists in a state of C2′/C3′-endo dynamic equilibrium. 14,15,97Therefore, we conclude the complete inaccessibility of the C3′-endo pucker for the DNA strand in our hybrid simulations using all the AMBER ffs points to a significant imbalance.The absolute preference for C2′-endo leads to majority of the helical parameters shifting closer to the B-form (Table 2 and Supporting Information Table S2).Most significant is the reduction of inclination and roll, which in the case of the PPT structure clearly contradicts the experimental data (Table 2).
The inclination and roll are interrelated parameters in global helical and local base pair steps coordinate systems, respectively. 98,99We also noticed that the backbone dihedral angles of the hybrids correspond closely to the pure duplexes.In other words, the dihedrals observed for the hybrid RNA strand were virtually identical to the A-RNA duplex while the hybrid DNA strand was nearly identical to the B-DNA duplex, showing no discernible adaptation to the unique structural environment of the hybrids (Figure 3).The C3′-endo pucker of the hybrid DNA strand can be safely enforced by restraints (Table 1 and Supporting Information Figure S1) without undesired side-effects on the duplex stability.It was tested for the OL21/OL3 f f combination where it leads to significant shifts of the helical parameters toward the A-form (Table 2).We suggest that restraining the DNA sugar pucker in hybrids to C3′-endo might be helpful in simulations of hybrids using the AMBER f fs when such pucker is required, e.g., for correct interaction geometry with a protein or when detailing some catalytic mechanisms.In contrast, replacing the DNA χ dihedral potentials in OL15 with the RNA potentials from OL3 f f, as recently tested in ref 95 to improve the A/B-form balance of the DNA simulations, provided negligible improvements for the hybrid (Supporting Information Figure S1 and Table 2).We also note the absence of C3′-endo DNA pucker in AMBER simulations might not always be detrimental for the helical parameters as we observed quite high inclination in the low_dPyr simulations (see Methods and Table 2), despite the complete lack of the C3′-endo population on the DNA strand (Supporting Information Figure S1).However, the observed inclination could still be underestimated as we lack experimental reference for the low_dPyr structure (see Methods).It nevertheless hints at potential sequence-specific effects on f f performance for the hybrids, at least for selected helical parameters.However, due to the sheer scale of the imbalance, we do not expect the C3′-endo DNA pucker to be significantly populated for any hybrid structures in simulations using AMBER f fs.

CHARMM36 f f Reproduces the C3′-Endo Pucker of the Hybrids but Has Trouble to Keep Stable Base
Pairing.The simulations using the CHARMM36 f f accurately reproduced the experimental PPT hybrid structure geometry, including the C3′-endo pucker of the DNA strand (Figure 2) and the inclination (Table 2).Contrasting the AMBER f fs, we also observed the DNA strand of the hybrid mostly possessing dihedrals typical for pure A-RNA (Figure 3), rationalizing the shift of the helical parameters toward the A-form (Table 2).The CHARMM36 performance we observed for the hybrids is very good and in agreement with the previous studies utilizing the CHARMM ff. 37,38,100However, we also noticed large instabilities of the base pairing in all our CHARMM36 simulations (Supporting Information Figure S2 and Table S3).Namely, several consecutive base pairs near the terminus were temporarily or permanently disrupted in CHARMM36 simulations (both hybrids and pure duplexes).These large instabilities occurred despite the terminal base pairs being supported with the sHBfix (see Methods), which was sufficient to almost entirely prevent terminal base pair fraying in simulations using the AMBER f fs.While some degree of fraying is to be expected, 69 we suggest that such large scale instability on a microsecond time scale (up to 30% of simulation frames were affected by this), including the thermodynamically more stable GC/CG base pairs (Supporting Information Figure S2 and Table S3), likely reflects imbalance within the CHARMM36 f f rather than genuine biomolecular dynamics.
A possible way to overcome this problem and still benefit from the satisfactory performance of CHARMM36 f f for hybrids is to apply the sHBfix potentials 68 to increase the stability of every base pair and not just the terminal ones (Table 1).Our CHARMM36 simulations where we have applied sHBfix in such a fashion revealed negligible instability of the base pairs while fully maintaining the quite accurate description of the hybrid structure (Table 2).The advantage of stabilizing the base pairing with sHBfix instead of standard distance restraints is that, depending on the strength of the sHBfix, temporary disruptions of the base pairing and breathing motions are still allowed to happen albeit with lower frequency, leading to a more natural dynamics. 23olarizable f fs Can Populate C3′-Endo Pucker for DNA Nucleotides but Significantly Underestimate the Inclination in Both DD and PPT Hybrids.The polarizable DRUDE and AMOEBA f fs sample a minor population of the C3′-endo pucker in the DNA strand of the hybrid systems (Figure 2).Notably, in the DD_hybrid simulations using the AMOEBA f f, the C3′-endo population was even dominant (Supporting Information Figure S1).However, surprisingly, the population of the DNA C3′-endo pucker did not translate into increased base pair inclination of the hybrids (Table 2).In fact, the inclination was generally lower than with the AMBER family of f fs where the C3′-endo pucker for DNA was completely absent and where higher inclination could be achieved by restraining the C3′-endo pucker (see above).The other helical parameters were also farther away from the experimental values in comparison with the simulations using the nonpolarizable f fs.We noticed the polarizable and nonpolarizable f fs populating different values of the β backbone dihedral (Supporting Information Figure S3 and Table S4), which might be one of the reasons behind the suboptimal reproduction of the helical parameters.The ability to reproduce the C3′-endo pucker of the DNA strand is promising, but it is also evident the new polarizable f fs require an extensive testing over a broad set of systems, while some adjustments of the existing versions might be necessary. 101dditional Comments on Specific f fs.We have also observed some problems with specific f fs and their combinations (Table 1).Namely, the ROC-RNA f f reversibly but significantly populated the α(trans)/γ(trans) dihedral region of the RNA backbone (Figure 4a), both in pure RNA and in hybrids.Such a large population of this dihedral angle combination was not seen for any other f fs.It might have contributed to the low inclination and drift toward the B-form observed for the hybrids with this f f (Table 2), as by excluding the frames with α(trans)/γ(trans) states from the ensemble, the inclination increased.Second, while the Tumuc1 f f performed well for the B-DNA simulations, when combined with either OL3 or ROC-RNA, we observed a visible drift toward the B-form for the hybrids, once again reducing both the inclination and x-displacement below their experimental values (Table 2 and Figure 2).Lastly, the DES-Amber f f revealed a complete lack of backbone BII-states in pure B-DNA, as admitted also in the original paper. 26Consequently, the BII-states were absent also for the DNA strand in the hybrids (Figure 4b).The DES-Amber f f also revealed an unusual population of non-native pucker values in the C4′exo/O4′-endo region for DNA nucleotides in both B-DNA and hybrids (Figure 4c).This may have shifted some of the hybrid helical parameters closer to the experimental values as the pucker value was lowered closer to the experimentally indicated C3′-endo state, but did not reach it (Figure 4c and Table 2).In any case, C4′-exo/O4′-endo puckers are nonnative for B-DNA and their presence in simulations is likely undesirable for both B-DNA and the hybrids.
In general, we also recommend careful verification that the hybrid systems are correctly built as some f fs (e.g., DES-Amber f f, see Methods) can contain overlapping atomic types for the DNA and RNA nucleotides.A good practice is to always check that the same bonded potential energies of the individual strands are obtained for systems built with lone DNA and RNA strands and a single set of parameters loaded, and once again when the same molecules are part of the hybrid.Disagreement between the calculated energies indicates an overlap between the parameters.
Water Model and Ion Parameters Have Minimal Influence.We have utilized water models recommended by the f f developers, whose recommendations are given in the original papers.In other cases, we have used the SPC/E water model, 62 which is well tested for simulations of nucleic acids in general.With certain f f combinations, we have also tested the OPC water model 102 for the hybrid simulations (see Methods and Table 1).The results of the present simulations indicated only a minor sensitivity to the water model used, with the solute f f proving to be the decisive factor.Specifically, the hybrid simulations using the OPC water revealed a minor but systematic drop in inclination compared to the SPC/E simulations (Table 2 and Supporting Information Figure S4).While this further exacerbates the undesirable reduction of inclination we observe for the hybrids, it is a very minor effect.Although we did not test the water models systematically, we tentatively suggest that for simulations of hybrid duplexes, all the water models commonly employed for simulations of nucleic acids (e.g., TIP3P, SPC/E and OPC) are reasonable choices.The results are only moderately sensitive to the water model choice, echoing similar conclusions made for simulations of A-RNA duplexes, 80,99,103 but contrasting the rather significant solvent effects observed for the RNA tetranucleotides or G-quadruplexes. 104,105Therefore, it might be best to follow the recommendation made by the authors of the particular nucleic acid f f, provided any are given.Obviously, this can lead to incompatible water model requirements for some DNA and RNA f f combinations.For instance, the DES-Amber ff utilizes the TIP4P-D water model along with the CHARMM22 ions, which is however an untested setting with all the other AMBER f fs.Lastly, with certain f f combinations, we have also tested the Li and Merz ion parameters 94 for the hybrid simulations, observing that the choice of ion parameters has even smaller influence than the water model (Table 2 and Supporting Information Figure S4).
Hybrid Duplexes Display Buckle Anisotropy.An intriguing difference between the pure duplexes and hybrids was observed for the buckle, which describes the V-shaped deformation of the base pairs (Figure 5a).Namely, the distribution of the buckle peaked around 0°for both A-RNA and B-DNA, meaning that on average, the base pairs did not have a preferred direction for the V-shaped deformation as they fluctuated in both directions equally.In contrast, the buckle peaked at ∼−5°for the hybrids, making the base pairs

Journal of Chemical Theory and Computation
buckled in the 3′ direction of the RNA strand on average (Table 2 and Figure 5).Examination of the helical parameters in nine different hybrid duplexes deposited in the PDB database revealed the buckle anisotropy to be a universal structural feature for the hybrids (Supporting Information Table S1).All tested f fs managed to qualitatively reproduce the buckle anisotropy, but the buckle value significantly varied among the individual f f combinations (Table 2).While seemingly a genuine structural effect, we suspect that some f fs could be exaggerating it in order to relieve potential conflicts between the RNA and DNA f f parameters.Lastly, note that the tilt and tip parameters were similarly shifted for the hybrids as the buckle was (Supporting Information Table S2).

■ CONCLUSIONS
In this study, we present a comprehensive MD simulation benchmark focused on hybrid DNA/RNA duplexes, utilizing modern, state-of-the-art nucleic acid f fs.Our findings suggest that none of the tested f fs provide a fully satisfactory performance.The OL-family of AMBER f fs (OL21, OL15 for DNA, and OL3 for RNA), the bsc1 (paired with OL3 for RNA) or CHARMM36, currently represents the most promising albeit far from perfect choices for simulating the hybrids.The universal limitation of all tested AMBER f fs was the inability to reproduce the C3′-endo pucker of the DNA strand in the hybrids, which is accompanied by reduced helical inclination.This appears to be primarily the issue of the DNA f f parameters.We show that when utilizing the AMBER f fs for hybrid simulations, it is straightforwardly possible to restrain the C3′-endo pucker for the DNA nucleotides, which greatly improves the inclination.The CHARMM36 f f provided satisfactory sugar pucker distributions and inclinations, but the f f severely underestimates stability of the base pairs in both hybrids and pure duplexes.It leads to large populations of structures with several disrupted base pairs, not only the terminal ones.We demonstrate that this issue can be resolved by reinforcing stability of the base pairs using sHBfix, which significantly improved simulation performance with CHARMM36.With the help of sHBfix, CHARMM36 provides the best results for the hybrid duplexes.The latest versions of polarizable DRUDE and AMOEBA f fs demonstrated a promising ability to populate the C3′-endo pucker as a minor state and to describe the pucker transitions.However, their inclination was even lower than with the nonpolarizable f fs, possibly requiring further f f refinement, perhaps of the backbone dihedral potentials.In summary, the results show that accurate simulations of hybrid DNA/RNA duplexes are challenging for contemporary f fs.However, we underline that effectiveness of the nucleic acids f fs is highly dependent on the specific molecular systems being studied.Therefore, different simulation performance could potentially be observed for DNA/RNA structures other than hybrid duplexes (not explored in this work).However, in the specific case of the hybrid duplexes, it is apparent that additional parametrization efforts would have to be made to develop a perfect f f.At the very least, DNA/RNA hybrids represent a straightforward test system that should be considered in future studies aimed at ff refinement, especially of the DNA parameters.Still, satisfactory performance for the hybrid simulations can be achieved with current generation of f fs, provided the authors are cognizant of their specific caveats and apply the appropriate workarounds, as outlined in this work.

■ ASSOCIATED CONTENT
* sı Supporting Information

Figure 1 .
Figure 1.Dickerson−Drew dodecamer (DD) and polypurine tract hybrid duplex (PPT) structures and sequences.The DNA and RNA sequence letters are colored blue and red, respectively.

Figure 2 .
Figure 2. Histograms of the DNA sugar puckers, inclination and x-displacement in MD simulations of the PPT structure.Example of nucleotides possessing C3′-endo and C2′-endo puckers is shown at top left corner.The individual f f combinations are color-coded in the graphs according to the legend in the bottom right corner.Black dots on the x-axes indicate average experimental values.

Figure 3 .
Figure 3. Histograms of the dihedral angles in MD simulations of the DD structures.Distribution of the ζ (zeta) dihedral angle in pure DNA and RNA and the hybrid DD structures, using the OL21/OL3 (top) and CHARMM36 f fs (bottom).The visualized data sets are representative of the development observed also for the other dihedrals.The top graph is representative of a performance seen also with the other AMBER ffs.

Figure 4 .
Figure 4. Histograms of the sugar puckers and dihedral angles in MD simulations of the DD structures.(a) α and γ dihedrals observed for the RNA strand in simulations of the DD_RNA and DD_hybrid structures using the ROC-RNA f f (combined with OL21 f f for the hybrid).(b) ζ dihedral angle in DD_DNA and DD_hybrid structures, using the OL21/OL3 and DES-Amber f fs.Correlated behavior was observed for the ε dihedral as it pertains to the BI/BII B-form substates.(c) Sugar puckers observed for the DD structures using the DES-Amber f f.

Figure 5 .
Figure 5. Anisotropic base pair buckle observed in the hybrid DNA/ RNA experimental structures and simulations.(a) 2D scheme of the anisotropic buckle of the hybrids.Note that the extent of the Vshaped deformation of the base pairs is greatly exaggerated for illustrative purposes.(b) Distribution of the buckle in MD simulations of the PPT structure.Black dot on the x-axis indicates the average experimental value.

Table 1 .
List of Simulations

Table 2 .
Selected Helical Parameters Observed in the MD Simulations