The computational efficiency of approximate quantum mechanical methods allows their use for the construction of multidimensional reaction free energy profiles. It has recently been demonstrated that quantum models based on the neglect of diatomic differential overlap (NNDO) approximation have difficulty modeling deoxyribose and ribose sugar ring puckers and thus limit their predictive value in the study of RNA and DNA systems. A method has been introduced in our previous work to improve the description of the sugar puckering conformational landscape that uses a multidimensional B-spline correction map (BMAP correction) for systems involving intrinsically coupled torsion angles. This method greatly improved the adiabatic potential energy surface profiles of DNA and RNA sugar rings relative to high-level ab initio methods even for highly problematic NDDO-based models. In the present work, a BMAP correction is developed, implemented, and tested in molecular dynamics simulations using the AM1/d-PhoT semiempirical Hamiltonian for biological phosphoryl transfer reactions. Results are presented for gas-phase adiabatic potential energy surfaces of RNA transesterification model reactions and condensed-phase QM/MM free energy surfaces for nonenzymatic and RNase A-catalyzed transesterification reactions. The results show that the BMAP correction is stable, efficient, and leads to improvement in both the potential energy and free energy profiles for the reactions studied, as compared with ab initio and experimental reference data. Exploration of the effect of the size of the quantum mechanical region indicates the best agreement with experimental reaction barriers occurs when the full CpA dinucleotide substrate is treated quantum mechanically with the sugar pucker correction.

AbstractRecently, a crystal structure has been reported of a new catalytic RNA, the TS ribozyme, that has been identified through comparative genomics and is believed to be a metalloribozyme having novel mechanistic features. Although this data provides invaluable structural information, analysis suggests a conformational change is required to arrive at a catalytically relevant state. We report results of molecular simulations that predict a spontaneous local rearrangement of the active site, leading to solution structures consistent with available functional data and providing competing mechanistic hypotheses that can be experimentally tested. The two competing hypotheses differ in the proposed identity of the catalytic general acid: either a water molecule coordinating a Mg2+ ion bound at the Watson–Crick edge of residue C7, or the N3 position of residue C7 itself.

We present results from molecular dynamics simulations and free energy calculations of the twister ribozyme at different stages along the reaction path to gain insight into its mechanism. The results, together with recent biochemical experiments, provide support for a mechanism involving general-acid catalysis by a conserved adenine residue in the active site. Although adenine has been previously implicated as a general acid acting through the N1 position in other ribozymes such as the hairpin and VS ribozymes, in the twister ribozyme there may be a twist. Biochemical experiments suggest that general acid catalysis may occur through the N3 position, which has never before been implicated in this role; however, currently, there is a lack of a detailed structural model for the active state of the twister ribozyme in solution that is consistent with these and other experiments. Simulations in a crystalline environment reported here are consistent with X-ray crystallographic data, and suggest that crystal packing contacts trap the RNA in an inactive conformation with U-1 in an extruded state that is incompatible with an in-line attack to the scissile phosphate. Simulations in solution, on the other hand, reveal this region to be dynamic and able to adopt a conformation where U-1 is stacked with G33. In this state, the nucleophile is in line with the scissile phosphate, and the N1 position of G33 and N3 position of A1 are poised to act as a general base and acid, respectively, as supported by mutational experiments. Free energy calculations further predict the electrostatic environment causes a shift of the microscopic pKa at the N3 position of A1 toward neutrality by approximately 5 pKa units. These results offer a unified interpretation of a broad range of currently available experimental data that points to a novel mode of general acid catalysis through the N3 position of an adenine nucleobase, thus expanding the repertoire of known mechanistic strategies employed by small nucleolytic ribozymes.

RNA enzymes serve as a potentially powerful platform from which to design catalysts and engineer new biotechnology. A fundamental understanding of these systems provides insight to guide design. The hepatitis delta virus ribozyme (HDVr) is a small, self-cleaving RNA motif widely distributed in nature, which has served as a paradigm for understanding the basic principles of RNA catalysis. Nevertheless, questions remain regarding the precise roles of divalent metal ions and key nucleotides in catalysis. In an effort to establish a reaction mechanism model consistent with available experimental data, we utilize molecular dynamics simulations to explore different conformations and metal ion binding modes along the HDVr reaction path. Building upon recent crystallographic data, our results provide a dynamic model of the HDVr reaction mechanism involving a conformational switch between multiple noncanonical G25:U20 base pair conformations in the active site. These local nucleobase dynamics play an important role in catalysis by modulating the metal binding environments of two Mg2+ ions that support catalysis at different steps of the reaction pathway. The first ion plays a structural role by inducing a base pair flip necessary to obtain the catalytic fold in which C75 moves towards to the scissile phosphate in the active site. Ejection of this ion then permits a second ion to bind elsewhere in the active site and facilitate nucleophile activation. The simulations collectively describe a mechanistic scenario that is consistent with currently available experimental data from crystallography, phosphorothioate substitutions, and chemical probing studies. Avenues for further experimental verification are suggested.Keywords: conformation switching; HDV; HDVr; metal ions; QM/MM; reaction mechanism; ribozyme

A new approach for performing Particle Mesh Ewald in ab initio quantum mechanical/molecular mechanical (QM/MM) simulations with extended atomic orbital basis sets is presented. The new approach, the Ambient-Potential Composite Ewald (CEw) method, does not perform the QM/MM interaction with Mulliken charges nor electrostatically fit charges. Instead the nuclei and electron density interact directly with the MM environment, but in a manner that avoids the use of dense Fourier transform grids. By performing the electrostatics with the underlying QM density, the CEw method avoids self-consistent field instabilities that have been encountered with simple charge mapping procedures. Potential of mean force (PMF) profiles of the p-nitrophenyl phosphate dissociation reaction in explicit solvent are computed from PBE0/6-31G* QM/MM molecular dynamics simulations with various electrostatic protocols. The CEw profiles are shown to be stable with respect to real-space Ewald cutoff, whereas the PMFs computed from truncated and switched electrostatics produce artifacts. PBE0/6-311G**, AM1/d-PhoT, and DFTB2 QM/MM simulations are performed to generate two-dimensional PMF profiles of the phosphoryl transesterification reactions with ethoxide and phenoxide leaving groups. The semiempirical models incorrectly produce a concerted ethoxide mechanism, whereas PBE0 correctly produces a stepwise mechanism. The ab initio reaction barriers agree more closely to experiment than the semiempirical models. The failure of Mulliken-charge QM/MM-Ewald is analyzed.

In drug discovery, protonation states and tautomerization are easily overlooked. Through a Merck–Rutgers collaboration, this paper re-examined the initial settings and preparations for the Thermodynamic Integration (TI) calculation in AMBER Free-Energy Workflows, demonstrating the value of careful consideration of ligand protonation and tautomer state. Finally, promising results comparing AMBER TI and Schrödinger FEP+ are shown that should encourage others to explore the value of TI in routine Structure-based Drug Design.

Divalent metal ions are of fundamental importance to the function and folding of nucleic acids. Divalent metal ion–nucleic acid interactions are complex in nature and include both territorial and site specific binding. Commonly employed nonbonded divalent ion models, however, are often parametrized against bulk ion properties and are subsequently utilized in biomolecular simulations without considering any data related to interactions at specific nucleic acid sites. Previously, we assessed the ability of 17 different nonbonded Mg2+ ion models to reproduce different properties of Mg2+ in aqueous solution including radial distribution functions, solvation free energies, water exchange rates, and translational diffusion coefficients. In the present work, we depart from the recently developed 12–6–4 potential models for divalent metal ions developed by Li and Merz and tune the pairwise parameters for Mg2+, Mn2+, Zn2+, and Cd2+ binding dimethyl phosphate, adenosine, and guanosine in order to reproduce experimental site specific binding free energies derived from potentiometric pH titration data. We further apply these parameters to investigate a metal ion migration previously proposed to occur during the catalytic reaction of the hammerhead ribozyme. The new parameters are shown to be accurate and balanced for nucleic acid binding in comparison with available experimental data and provide an important tool for molecular dynamics and free energy simulations of nucleic acids where these ions may exhibit different binding modes.