Recent modifications and improvements to standard nucleic acid force fields have attempted to fix problems and issues that have been observed as longer timescale simulations have become routine. Although previous work has shown the ability to fold the UUCG stem-loop structure, until now no group has attempted to quantify the performance of current force fields using highly converged structural populations of the tetraloop conformational ensemble. In this study, we report the use of multiple independent sets of multidimensional replica exchange molecular dynamics (M-REMD) simulations with different initial conditions to generate well-converged conformational ensembles for the tetranucleotides r(GACC) and r(CCCC), as well as the larger UUCG tetraloop motif. By generating what is to our knowledge the most complete RNA structure ensembles reported to date for these systems, we remove the coupling between force field errors and errors due to incomplete sampling, providing a comprehensive comparison between current top-performing MD force fields for RNA. Of the RNA force fields tested in this study, none demonstrate the ability to correctly identify the most thermodynamically stable structure for all three systems. We discuss the deficiencies present in each potential function and suggest areas where improvements can be made. The results imply that although "short" (nsec-μsec timescale) simulations may stay close to their respective experimental structures and may well reproduce experimental observables, inevitably the current force fields will populate alternative incorrect structures that are more stable than those observed via experiment.

We compare the performance of two different RNA force fields in four water models in simulating the conformational ensembles r(GACC) and r(CCCC). With the increased sampling facilitated by multidimensional replica exchange molecular dynamics (M-REMD), populations are compared to NMR data to evaluate force field reliability. The combination of AMBER ff12 with vdW(bb) modifications and the OPC water model produces results in quantitative agreement with the NMR ensemble that have eluded us to date.

A necessary step to properly assess and validate the performance of force fields for biomolecules is to exhaustively sample the accessible conformational space, which is challenging for large RNA structures. Given questions regarding the reliability of modeling RNA structure and dynamics with current methods, we have begun to use RNA tetranucleotides to evaluate force fields. These systems, though small, display considerable conformational variability and complete sampling with standard simulation methods remains challenging. Here we compare and discuss the performance of known variations of replica exchange molecular dynamics (REMD) methods, specifically temperature REMD (T-REMD), Hamiltonian REMD (H-REMD), and multidimensional REMD (M-REMD) methods, which have been implemented in Amber’s accelerated GPU code. Using two independent simulations, we show that M-REMD not only makes very efficient use of emerging large-scale GPU clusters, like Blue Waters at the University of Illinois, but also is critically important in generating the converged ensemble more efficiently than either T-REMD or H-REMD. With 57.6 μs aggregate sampling of a conformational ensemble with M-REMD methods, the populations can be compared to NMR data to evaluate force field reliability and further understand how putative changes to the force field may alter populations to be in more consistent agreement with experiment.

Many problems studied via molecular dynamics require accurate estimates of various thermodynamic properties, such as the free energies of different states of a system, which in turn requires well-converged sampling of the ensemble of possible structures. Enhanced sampling techniques are often applied to provide faster convergence than is possible with traditional molecular dynamics simulations. Hamiltonian replica exchange molecular dynamics (H-REMD) is a particularly attractive method, as it allows the incorporation of a variety of enhanced sampling techniques through modifications to the various Hamiltonians. In this work, we study the enhanced sampling of the RNA tetranucleotide r(GACC) provided by H-REMD combined with accelerated molecular dynamics (aMD), where a boosting potential is applied to torsions, and compare this to the enhanced sampling provided by H-REMD in which torsion potential barrier heights are scaled down to lower force constants. We show that H-REMD and multidimensional REMD (M-REMD) combined with aMD does indeed enhance sampling for r(GACC), and that the addition of the temperature dimension in the M-REMD simulations is necessary to efficiently sample rare conformations. Interestingly, we find that the rate of convergence can be improved in a single H-REMD dimension by simply increasing the number of replicas from 8 to 24 without increasing the maximum level of bias. The results also indicate that factors beyond replica spacing, such as round trip times and time spent at each replica, must be considered in order to achieve optimal sampling efficiency.

Base flipping is a common strategy utilized by many enzymes to gain access to the functional groups of nucleic acid bases in duplex DNA which are otherwise protected by the DNA backbone and hydrogen bonding with their partner bases. Several X-ray crystallography studies have revealed flipped conformations of nucleotides bound to enzymes. However, little is known about the base flipping process itself, even less about the role of the enzymes. Computational studies have used umbrella sampling to elicit the free energy profile of the base flipping process using a pseudo-dihedral angle to represent the reaction coordinate. In this study we have used an unrestrained trajectory in which a flipped base spontaneously reinserted into the helix in order to evaluate and improve the previously defined pseudo-dihedral angle. Our modified pseudo-dihedral angles use a new atom selection to improve the numerical stability of the restraints and also provide better correlation to the extent of flipping observed in simulations. Furthermore, based on the comparison of PMFs generated using different reaction coordinates, we observed that the shape of a flipping PMF profile is strongly dependent on the definition of the reaction coordinate, even for the same data set.

The nudged elastic band method (NEB) can be used to find a minimum energy path between two given starting structures. This method has been available in the standard release of the Amber9 and Amber10 suite of programs. In this paper a novel implementation of this method will be discussed, in which the nudged elastic band method is applied to only a specific, user-defined subset of atoms in a particular system, returning comparable results and minimum energy pathways as the standard implementation for an alanine dipeptide test system. This allows incorporation of explicit solvent with simulated systems, which may be preferred in many cases to an implicit solvent model. From a computational standpoint, this implementation of NEB also reduces the communication overhead inside the code, resulting in better performance for larger systems.

Base eversion is a fundamental process in the biochemistry of nucleic acids, allowing proteins engaged in DNA repair and epigenetic modifications to access target bases in DNA. Crystal structures reveal endpoints of these processes, but not the pathways involved in the dynamic process of base recognition. To elucidate the pathway taken by 8-oxoguanine during base excision repair by Fpg, we calculated free energy surfaces during eversion of the damaged base through the major and minor grooves. The minor groove pathway and free energy barrier (6–7 kcal/mol) are consistent with previously reported results.1 However, eversion of 8-oxoG through the major groove encounters a significantly lower barrier (3–4 kcal/mol) more consistent with experimentally determined2 rates of enzymatic sliding during lesion search. Major groove eversion has been suggested for other glycosylases, suggesting that in addition to function, dynamics of base eversion may also be conserved.

Formamidopyrimidine-DNA glycosylase (Fpg) excises 8-oxoguanine (oxoG) from DNA but ignores normal guanine. We combined molecular dynamics simulation and stopped-flow kinetics with fluorescence detection to track the events in the recognition of oxoG by Fpg and its mutants with a key phenylalanine residue, which intercalates next to the damaged base, changed to either alanine (F110A) or fluorescent reporter tryptophan (F110W). Guanine was sampled by Fpg, as evident from the F110W stopped-flow traces, but less extensively than oxoG. The wedgeless F110A enzyme could bend DNA but failed to proceed further in oxoG recognition. Modeling of the base eversion with energy decomposition suggested that the wedge destabilizes the intrahelical base primarily through buckling both surrounding base pairs. Replacement of oxoG with abasic (AP) site rescued the activity, and calculations suggested that wedge insertion is not required for AP site destabilization and eversion. Our results suggest that Fpg, and possibly other DNA glycosylases, convert part of the binding energy into active destabilization of their substrates, using the energy differences between normal and damaged bases for fast substrate discrimination.