Machine learning potentials (MLPs) have revolutionized molecular simulation by providing efficient and accurate models for predicting atomic interactions. MLPs continue to advance and have had profound impact in applications that include drug discovery, enzyme catalysis, and materials design. The current landscape of MLP software presents challenges due to the limited interoperability between packages, which can lead to inconsistent benchmarking practices and necessitates separate interfaces with molecular dynamics (MD) software. To address these issues, we present DeePMD-GNN, a plugin for the DeePMD-kit framework that extends its capabilities to support external graph neural network (GNN) potentials.DeePMD-GNN enables the seamless integration of popular GNN-based models, such as NequIP and MACE, within the DeePMD-kit ecosystem. Furthermore, the new software infrastructure allows GNN models to be used within combined quantum mechanical/molecular mechanical (QM/MM) applications using the range corrected ΔMLP formalism.We demonstrate the application of DeePMD-GNN by performing benchmark calculations of NequIP, MACE, and DPA-2 models developed under consistent training conditions to ensure fair comparison.

MTR1 is an in vitro-selected alkyl transferase ribozyme that transfers an alkyl group from O6-alkylguanine to N1 of the target adenine in the RNA substrate (A63). The structure of the ribozyme suggested a mechanism in which a cytosine (C10) acts as a general acid to protonate O6-alkylguanine N1. Here, we have analyzed the role of the C10 general acid and the A63 nucleophile by atomic mutagenesis and computation. C10 was substituted by n1c and n1c, c5n variants. The n1c variant has an elevated pKa (11.4 as the free nucleotide) and leads to a 104-fold lower activity that is pH-independent. Addition of the second c5n substitution with a lower pKa restored both the rate and pH dependence of alkyl transfer. Quantum mechanical calculations indicate that protonation of O6-alkylguanine lowers the barrier to alkyl transfer and that there is a significantly elevated barrier to proton transfer for the n1c single substitution. The calculated pKa values are in good agreement with the apparent values from measured rates. Increasing the pKa of the nucleophile by A63 n7c substitution led to a 6-fold higher rate. The increased reactivity of the nucleophile corresponds to a βnuc of ∼0.5, indicating significant C–N bond formation in the transition state. Taken together, these results are consistent with a two-step mechanism comprising protonation of the O6-alkylguanine followed by alkyl transfer.

Accurate in silico predictions of how strongly small molecules bind to proteins, such as those afforded by relative binding free energy (RBFE) calculations, can greatly increase the efficiency of the hit-to-lead and lead optimization stages of the drug discovery process. The success of such calculations, however, relies heavily on their precision. Here, we show that a recently developed alchemical enhanced sampling (ACES) approach can consistently improve the precision of RBFE calculations on a large and diverse set of proteins and small molecule ligands. The addition of ACES to conventional RBFE calculations lowered the average hysteresis by over 35% (0.3–0.4 kcal/mol) and the average replicate spread by over 25% (0.2–0.3 kcal/mol) across a set of 10 protein targets and 213 small molecules while maintaining similar or improved accuracy. We show in atomic detail how ACES improved convergence of several representative RBFE calculations through enhancing the sampling of important slowly transitioning ligand degrees of freedom.

Antimicrobial resistance in bacteria often arises from their ability to actively identify and expel toxic compounds. The bacterium strain Pseudomonas putida DOT-T1E utilizes its TtgABC efflux pump to confer robust resistance against antibiotics, flavonoids, and organic solvents. This resistance mechanism is intricately regulated at the transcriptional level by the TtgR protein. Through molecular dynamics and alchemical free energy simulations, we systematically examine the binding of seven flavonoids and their derivatives with the TtgR transcriptional regulator. Our simulations reveal distinct binding geometries and free energies for the flavonoids in the active site of the protein, which are driven by a range of noncovalent forces encompassing van der Waals, electrostatic, and hydrogen bonding interactions. The interplay of molecular structures, substituent patterns, and intermolecular interactions effectively stabilizes the bound flavonoids, confining their movements within the TtgR binding pocket. These findings yield valuable insights into the molecular determinants that govern ligand recognition in TtgR and shed light on the mechanism of antimicrobial resistance in P. putida DOT-T1E.

We present software infrastructure for the design and testing of new quantum mechanical/molecular mechanical and machine-learning potential (QM/MM−ΔMLP) force fields for a wide range of applications. The software integrates Amber’s molecular dynamics simulation capabilities with fast, approximate quantum models in the xtb package and machine-learning potential corrections in DeePMD-kit. The xtb package implements the recently developed density-functional tight-binding QM models with multipolar electrostatics and density-dependent dispersion (GFN2-xTB), and the interface with Amber enables their use in periodic boundary QM/MM simulations with linear-scaling QM/MM particle-mesh Ewald electrostatics. The accuracy of the semiempirical models is enhanced by including machine-learning correction potentials (ΔMLPs) enabled through an interface with the DeePMD-kit software. The goal of this paper is to present and validate the implementation of this software infrastructure in molecular dynamics and free energy simulations. The utility of the new infrastructure is demonstrated in proof-of-concept example applications. The software elements presented here are open source and freely available. Their interface provides a powerful enabling technology for the design of new QM/MM−ΔMLP models for studying a wide range of problems, including biomolecular reactivity and protein–ligand binding.

We report the development and testing of new integrated cyberinfrastructure for performing free energy simulations with generalized hybrid quantum mechanical/molecular mechanical (QM/MM) and machine learning potentials (MLPs) in Amber. The Sander molecular dynamics program has been extended to leverage fast, density-functional tight-binding models implemented in the DFTB+ and xTB packages, and an interface to the DeePMD-kit software enables the use of MLPs. The software is integrated through application program interfaces that circumvent the need to perform “system calls” and enable the incorporation of long-range Ewald electrostatics into the external software’s self-consistent field procedure. The infrastructure provides access to QM/MM models that may serve as the foundation for QM/MM–ΔMLP potentials, which supplement the semiempirical QM/MM model with a MLP correction trained to reproduce ab initio QM/MM energies and forces. Efficient optimization of minimum free energy pathways is enabled through a new surface-accelerated finite-temperature string method implemented in the FE-ToolKit package. Furthermore, we interfaced Sander with the i-PI software by implementing the socket communication protocol used in the i-PI client–server model. The new interface with i-PI allows for the treatment of nuclear quantum effects with semiempirical QM/MM–ΔMLP models. The modular interoperable software is demonstrated on proton transfer reactions in guanine-thymine mispairs in a B-form deoxyribonucleic acid helix. The current work represents a considerable advance in the development of modular software for performing free energy simulations of chemical reactions that are important in a wide range of applications.

Rare tautomeric forms of nucleobases can lead to Watson–Crick-like (WC-like) mispairs in DNA, but the process of proton transfer is fast and difficult to detect experimentally. NMR studies show evidence for the existence of short-time WC-like guanine–thymine (G-T) mispairs; however, the mechanism of proton transfer and the degree to which nuclear quantum effects play a role are unclear. We use a B-DNA helix exhibiting a wGT mispair as a model system to study tautomerization reactions. We perform ab initio (PBE0/6-31G*) quantum mechanical/molecular mechanical (QM/MM) simulations to examine the free energy surface for tautomerization. We demonstrate that while the ab initio QM/MM simulations are accurate, considerable sampling is required to achieve high precision in the free energy barriers. To address this problem, we develop a QM/MM machine learning potential correction (QM/MM-ΔMLP) that is able to improve the computational efficiency, greatly extend the accessible time scales of the simulations, and enable practical application of path integral molecular dynamics to examine nuclear quantum effects. We find that the inclusion of nuclear quantum effects has only a modest effect on the mechanistic pathway but leads to a considerable lowering of the free energy barrier for the GT*⇌G*T equilibrium. Our results enable a rationalization of observed experimental data and the prediction of populations of rare tautomeric forms of nucleobases and rates of their interconversion in B-DNA.

An alchemical enhanced sampling (ACES) method has recently been introduced to facilitate importance sampling in free energy simulations. The method achieves enhanced sampling from Hamiltonian replica exchange within a dual topology framework while utilizing new smoothstep softcore potentials. A common sampling problem encountered in lead optimization is the functionalization of aromatic rings that exhibit distinct conformational preferences when interacting with the protein. It is difficult to converge the distribution of ring conformations due to the long time scale of ring flipping events; however, the ACES method addresses this issue by modeling the syn and anti ring conformations within a dual topology. ACES thereby samples the conformer distributions by alchemically tunneling between states, as opposed to traversing a physical pathway with a high rotational barrier. We demonstrate the use of ACES to overcome conformational sampling issues involving ring flipping in ML300-derived noncovalent inhibitors of SARS-CoV-2 Main Protease (Mpro). The demonstrations explore how the use of replica exchange and the choice of softcore selection affects the convergence of the ring conformation distributions. Furthermore, we examine how the accuracy of the calculated free energies is affected by the degree of phase space overlap (PSO) between adjacent states (i.e., between neighboring λ-windows) and the Hamiltonian replica exchange acceptance ratios. Both of these factors are sensitive to the spacing between the intermediate states. We introduce a new method for choosing a schedule of λ values. The method analyzes short “burn-in” simulations to construct a 2D map of the nonlocal PSO. The schedule is obtained by optimizing an alchemical pathway on the 2D map that equalizes the PSO between the λ intervals. The optimized phase space overlap λ-spacing method (Opt-PSO) leads to more numerous end-to-end single passes and round trips due to the correlation between PSO and Hamiltonian replica exchange acceptance ratios. The improved exchange statistics enhance the efficiency of ACES method. The method has been implemented into the FE-ToolKit software package, which is freely available.

We present a surface-accelerated string method (SASM) to efficiently optimize low-dimensional reaction pathways from the sampling performed with expensive quantum mechanical/molecular mechanical (QM/MM) Hamiltonians. The SASM accelerates the convergence of the path using the aggregate sampling obtained from the current and previous string iterations, whereas approaches like the string method in collective variables (SMCV) or the modified string method in collective variables (MSMCV) update the path only from the sampling obtained from the current iteration. Furthermore, the SASM decouples the number of images used to perform sampling from the number of synthetic images used to represent the path. The path is optimized on the current best estimate of the free energy surface obtained from all available sampling, and the proposed set of new simulations is not restricted to being located along the optimized path. Instead, the umbrella potential placement is chosen to extend the range of the free energy surface and improve the quality of the free energy estimates near the path. In this manner, the SASM is shown to improve the exploration for a minimum free energy pathway in regions where the free energy surface is relatively flat. Furthermore, it improves the quality of the free energy profile when the string is discretized with too few images. We compare the SASM, SMCV, and MSMCV using 3 QM/MM applications: a ribozyme methyltransferase reaction using 2 reaction coordinates, the 2′-O-transphosphorylation reaction of Hammerhead ribozyme using 3 reaction coordinates, and a tautomeric reaction in B-DNA using 5 reaction coordinates. We show that SASM converges the paths using roughly 3 times less sampling than the SMCV and MSMCV methods. All three algorithms have been implemented in the FE-ToolKit package made freely available.

AmberTools is a free and open-source collection of programs used to set up, run, and analyze molecular simulations. The newer features contained within AmberTools23 are briefly described in this Application note.

This Perspective provides a contextual explanation of the current state-of-the-art alchemical free energy methods and their role in drug discovery as well as highlights select emerging technologies. The narrative attempts to answer basic questions about what goes on “under the hood” in free energy simulations and provide general guidelines for how to run simulations and analyze the results. It is the hope that this work will provide a valuable introduction to students and scientists in the field.

DeePMD-kit is a powerful open-source software package that facilitates molecular dynamics simulations using machine learning potentials known as Deep Potential (DP) models. This package, which was released in 2017, has been widely used in the fields of physics, chemistry, biology, and material science for studying atomistic systems. The current version of DeePMD-kit offers numerous advanced features, such as DeepPot-SE, attention-based and hybrid descriptors, the ability to fit tensile properties, type embedding, model deviation, DP-range correction, DP long range, graphics processing unit support for customized operators, model compression, non-von Neumann molecular dynamics, and improved usability, including documentation, compiled binary packages, graphical user interfaces, and application programming interfaces. This article presents an overview of the current major version of the DeePMD-kit package, highlighting its features and technical details. Additionally, this article presents a comprehensive procedure for conducting molecular dynamics as a representative application, benchmarks the accuracy and efficiency of different models, and discusses ongoing developments.

Pistol ribozyme (Psr) is a distinct class of small endonucleolytic ribozymes, which are important experimental systems for defining fundamental principles of RNA catalysis and designing valuable tools in biotechnology. High-resolution structures of Psr, extensive structure–function studies, and computation support a mechanism involving one or more catalytic guanosine nucleobases acting as a general base and divalent metal ion-bound water acting as an acid to catalyze RNA 2′-O-transphosphorylation. Yet, for a wide range of pH and metal ion concentrations, the rate of Psr catalysis is too fast to measure manually and the reaction steps that limit catalysis are not well understood. Here, we use stopped-flow fluorescence spectroscopy to evaluate Psr temperature dependence, solvent H/D isotope effects, and divalent metal ion affinity and specificity unconstrained by limitations due to fast kinetics. The results show that Psr catalysis is characterized by small apparent activation enthalpy and entropy changes and minimal transition state H/D fractionation, suggesting that one or more pre-equilibrium steps rather than chemistry is rate limiting. Quantitative analyses of divalent ion dependence confirm that metal aquo ion pKa correlates with higher rates of catalysis independent of differences in ion binding affinity. However, ambiguity regarding the rate-limiting step and similar correlation with related attributes such as ionic radius and hydration free energy complicate a definitive mechanistic interpretation. These new data provide a framework for further interrogation of Psr transition state stabilization and show how thermal instability, metal ion insolubility at optimal pH, and pre-equilibrium steps such as ion binding and folding limit the catalytic power of Psr suggesting potential strategies for further optimization.

We use the modified Bigeleisen–Mayer equation to compute kinetic isotope effect values for non-enzymatic phosphoryl transfer reactions from classical and path integral molecular dynamics umbrella sampling. The modified form of the Bigeleisen–Mayer equation consists of a ratio of imaginary mode vibrational frequencies and a contribution arising from the isotopic substitution’s effect on the activation free energy, which can be computed from path integral simulation. In the present study, we describe a practical method for estimating the frequency ratio correction directly from umbrella sampling in a manner that does not require normal mode analysis of many geometry optimized structures. Instead, the method relates the frequency ratio to the change in the mass weighted coordinate representation of the minimum free energy path at the transition state induced by isotopic substitution. The method is applied to the calculation of ^16/18O and ^32/34S primary kinetic isotope effect values for six non-enzymatic phosphoryl transfer reactions. We demonstrate that the results are consistent with the analysis of geometry optimized transition state ensembles using the traditional Bigeleisen–Mayer equation. The method thus presents a new practical tool to enable facile calculation of kinetic isotope effect values for complex chemical reactions in the condensed phase.

A methyltransferase ribozyme (MTR1) was selected in vitro to catalyze alkyl transfer from exogenous O6-methylguanine (O6mG) to a target adenine N1, and recently, high-resolution crystal structures have become available. We use a combination of classical molecular dynamics, ab initio quantum mechanical/molecular mechanical (QM/MM) and alchemical free energy (AFE) simulations to elucidate the atomic-level solution mechanism of MTR1. Simulations identify an active reactant state involving protonation of C10 that hydrogen bonds with O6mG:N1. The deduced mechanism involves a stepwise mechanism with two transition states corresponding to proton transfer from C10:N3 to O6mG:N1 and rate-controlling methyl transfer (19.4  kcal·mol−1 barrier). AFE simulations predict the pKa for C10 to be 6.3, close to the experimental apparent pKa of 6.2, further implicating it as a critical general acid. The intrinsic rate derived from QM/MM simulations, together with pKa calculations, enables us to predict an activity–pH profile that agrees well with experiment. The insights gained provide further support for a putative RNA world and establish new design principles for RNA-based biochemical tools.

Modern semiempirical electronic structure methods have considerable promise in drug discovery as universal "force fields" that can reliably model biological and drug-like molecules. Herein, we compare the performance of several NDDO-based semiempirical (MNDO/d, AM1, PM6 and ODM2), density-functional tight-binding based (DFTB3, GFN1-xTB and GFN2-xTB) models with pure machine learning potentials (ANI-1x and ANI-2x) and hybrid quantum mechanical/machine learning potentials (AIQM1 and QDπ) for a wide range of data computed at a consistent ωB97X/6-31G* level of theory (as in the ANI-1x database). This data includes conformational energies, intermolecular interactions, tautomers, and protonation states. Additional comparisons are made to a set of natural and synthetic nucleic acids from the artificially expanded genetic information system (AEGIS). This dataset has important implications in the design of new biotechnology and therapeutics. Finally, weexamine acid/base chemistry relevant for RNA cleavage reactions catalyzed by small nucleolytic ribozymes and ribonucleases. Overall, the recently developed QDπ model performs exceptionally well across all datasets, having especially high accuracy for tautomers and protonation states relevant to drug discovery.

RNA strand cleavage by 2′-O-transphosphorylation is catalyzed not only by numerous nucleolytic RNA enzymes (ribozymes) but also by hydroxide or hydronium ions. In experiments, both cleavage of the 5′-linked nucleoside and isomerization between 3′,5′- and 2′,5′-phosphodiesters occur under acidic conditions, while only the cleavage reaction is observed under basic conditions. An ab initio path-integral approach for simulating kinetic isotope effects is used to reveal the reaction mechanisms for RNA cleavage and isomerization reactions under acidic conditions. Moreover, the proposed mechanisms can also be combined through the experimental pH-rate profiles.

Ribonucleases and small nucleolytic ribozymes are both able to catalyze RNA strand cleavage through 2′-O-transphosphorylation, provoking the question of whether protein and RNA enzymes facilitate mechanisms that pass through the same or distinct transition states. Here, we report the primary and secondary 18O kinetic isotope effects for hepatitis delta virus ribozyme catalysis that reveal a dissociative, metaphosphate-like transition state in stark contrast to the late, associative transition states observed for reactions catalyzed by specific base, Zn2+ ions, or ribonuclease A. This new information provides evidence for a discrete ribozyme active site design that modulates the RNA cleavage pathway to pass through an altered transition state.

A new direction has emerged in molecular simulations in recent years, where potential energy surfaces (PES) are constructed using machine learning (ML) methods. These ML models, combining the accuracy of quantum mechanical models and the efficiency of empirical atomic potential models, have been demonstrated by many studies to have extensive application prospects. This chapter introduces a recently developed ML model, Deep Potential (DP), and the corresponding package, DeePMD-kit. First, we present the basic theory of the DP method. Then, we show how to train and test a DP model for a gas-phase methane molecule using the DeePMD-kit package. Next, we introduce some recent progress on simulations of biomolecular processes by integrating the DeePMD-kit with the AMBER molecular simulation software suite. Finally, we provide a supplement on points that require further explanation.

We report QDπ-v1.0 for modeling the internal energy of drug molecules containing H, C, N, and O atoms. The QDπ model is in the form of a quantum mechanical/machine learning potential correction (QM/Δ-MLP) that uses a fast third-order self-consistent density-functional tight-binding (DFTB3/3OB) model that is corrected to a quantitatively high-level of accuracy through a deep-learning potential (DeepPot-SE). The model has the advantage that it is able to properly treat electrostatic interactions and handle changes in charge/protonation states. The model is trained against reference data computed at the ωB97X/6-31G* level (as in the ANI-1x data set) and compared to several other approximate semiempirical and machine learning potentials (ANI-1x, ANI-2x, DFTB3, MNDO/d, AM1, PM6, GFN1-xTB, and GFN2-xTB). The QDπ model is demonstrated to be accurate for a wide range of intra- and intermolecular interactions (despite its intended use as an internal energy model) and has shown to perform exceptionally well for relative protonation/deprotonation energies and tautomers. An example application to model reactions involved in RNA strand cleavage catalyzed by protein and nucleic acid enzymes illustrates QDπ has average errors less than 0.5 kcal/mol, whereas the other models compared have errors over an order of magnitude greater. Taken together, this makes QDπ highly attractive as a potential force field model for drug discovery.

We present an alchemical enhanced sampling (ACES) method implemented in the GPU-accelerated AMBER free energy MD engine. The methods hinges on the creation of an “enhanced sampling state” by reducing or eliminating selected potential energy terms and interactions that lead to kinetic traps and conformational barriers while maintaining those terms that curtail the need to otherwise sample large volumes of phase space. For example, the enhanced sampling state might involve transforming regions of a ligand and/or protein side chain into a noninteracting “dummy state” with internal electrostatics and torsion angle terms turned off. The enhanced sampling state is connected to a real-state end point through a Hamiltonian replica exchange (HREMD) framework that is facilitated by newly developed alchemical transformation pathways and smoothstep softcore potentials. This creates a counterdiffusion of real and enhanced-sampling states along the HREMD network. The effect of a differential response of the environment to the real and enhanced-sampling states is minimized by leveraging the dual topology framework in AMBER to construct a counterbalancing HREMD network in the opposite alchemical direction with the same (or similar) real and enhanced sampling states at inverted end points. The method has been demonstrated in a series of test cases of increasing complexity where traditional MD, and in several cases alternative REST2-like enhanced sampling methods, are shown to fail. The hydration free energy for acetic acid was shown to be independent of the starting conformation, and the values for four additional edge case molecules from the FreeSolv database were shown to have a significantly closer agreement with experiment using ACES. The method was further able to handle different rotamer states in a Cdk2 ligand identified as fractionally occupied in crystal structures. Finally, ACES was applied to T4-lysozyme and demonstrated that the side chain distribution of V111χ1 could be reliably reproduced for the apo state, bound to p-xylene, and in p-xylene→ benzene transformations. In these cases, the ACES method is shown to be highly robust and superior to a REST2-like enhanced sampling implementation alone.

We develop a framework for the design of optimized alchemical transformation pathways in free energy simulations using nonlinear mixing and a new functional form for so-called “softcore” potentials. We describe the implementation and testing of this framework in the GPU-accelerated AMBER software suite. The new optimized alchemical transformation pathways integrate a number of important features, including (1) the use of smoothstep functions to stabilize behavior near the transformation end points, (2) consistent power scaling of Coulomb and Lennard-Jones (LJ) interactions with unitless control parameters to maintain balance of electrostatic attractions and exchange repulsions, (3) pairwise form based on the LJ contact radius for the effective interaction distance with separation-shifted scaling, and (4) rigorous smoothing of the potential at the nonbonded cutoff boundary. The new softcore potential form is combined with smoothly transforming nonlinear λ weights for mixing specific potential energy terms, along with flexible λ-scheduling features, to enable robust and stable alchemical transformation pathways. The resulting pathways are demonstrated and tested, and shown to be superior to the traditional methods in terms of numerical stability and minimal variance of the free energy estimates for all cases considered. The framework presented here can be used to design new alchemical enhanced sampling methods, and leveraged in robust free energy workflows for large ligand data sets.

We report an automated workflow for production free-energy simulation setup and analysis (ProFESSA) using the GPU-accelerated AMBER free-energy engine with enhanced sampling features and analysis tools, part of the AMBER Drug Discovery Boost package that has been integrated into the AMBER22 release. The workflow establishes a flexible, end-to-end pipeline for performing alchemical free-energy simulations that brings to bear technologies, including new enhanced sampling features and analysis tools, to practical drug discovery problems. ProFESSA provides the user with top-level control of large sets of free-energy calculations and offers access to the following key functionalities: (1) automated setup of file infrastructure; (2) enhanced conformational and alchemical sampling with the ACES method; and (3) network-wide free-energy analysis with the optional imposition of cycle closure and experimental constraints. The workflow is applied to perform absolute and relative solvation free-energy and relative ligand–protein binding free-energy calculations using different atom-mapping procedures. Results demonstrate that the workflow is internally consistent and highly robust. Further, the application of a new network-wide Lagrange multiplier constraint analysis that imposes key experimental constraints substantially improves binding free-energy predictions.

We describe the generalized weighted thermodynamic perturbation (gwTP) method for estimating the free energy surface of an expensive “high-level” potential energy function from the umbrella sampling performed with multiple inexpensive “low-level” reference potentials. The gwTP method is a generalization of the weighted thermodynamic perturbation (wTP) method developed by Li and co-workers [J. Chem. Theory Comput. 2018, 14, 5583–5596] that uses a single “low-level” reference potential. The gwTP method offers new possibilities in model design whereby the sampling generated from several low-level potentials may be combined (e.g., specific reaction parameter models that might have variable accuracy at different stages of a multistep reaction). The gwTP method is especially well suited for use with machine learning potentials (MLPs) that are trained against computationally expensive ab initio quantum mechanical/molecular mechanical (QM/MM) energies and forces using active learning procedures that naturally produce multiple distinct neural network potentials. Simulations can be performed with greater sampling using the fast MLPs and then corrected to the ab initio level using gwTP. The capabilities of the gwTP method are demonstrated by creating reference potentials based on the MNDO/d and DFTB2/MIO semiempirical models supplemented with the “range-corrected deep potential” (DPRc). The DPRc parameters are trained to ab initio QM/MM data, and the potentials are used to calculate the free energy surface of stepwise mechanisms for nonenzymatic RNA 2′-O-transesterification model reactions. The extended sampling made possible by the reference potentials allows one to identify unequilibrated portions of the simulations that are not always evident from the short time scale commonly used with ab initio QM/MM potentials. We show that the reference potential approach can yield more accurate ab initio free energy predictions than the wTP method or what can be reasonably afforded from explicit ab initio QM/MM sampling.

With the rise in implementation of peer-led learning in higher education, the interactions between peer instructors and their students warrant further investigation as an increasingly critical component of student learning. In this work, Teaching Interns (TIs) are undergraduate peer instructors that lead supplemental learning sessions in General Chemistry. Each week, TIs participate in pedagogy training and complete written reflections on their learning sessions. For this multicase study, six TIs were observed in their office hours over varying time periods. A qualitative approach was taken to analyze their verbal behaviors and the extent to which those behaviors matched their beliefs about teaching. Specifically, discourse analysis allowed for the characterization of the interactions between TIs and students, while analysis of TIs’ weekly written reflections provided insight into their teaching beliefs and perception of their own teaching sessions. The results presented here suggest that even at the start of the program, TIs hold some productive beliefs about teaching, though these beliefs were not always evident in their interactions with students. Over time, TIs generally shifted toward more student-centered discourse and honed their abilities to convey or elicit deeper knowledge among their students. Further, evidence from the TIs’ reflections suggest that they became better at self-monitoring their own teaching behaviors, shrinking the gap between their practices and espoused beliefs about teaching, and that they turned their focus toward student learning versus simply managing their sessions. Taken together, this work provides additional support for the further development and study of peer instruction programs.

Electrostatic interactions are fundamental to RNA structure and function, and intimately influenced by solvation and the ion atmosphere. RNA enzymes, or ribozymes, are catalytic RNAs that are able to enhance reaction rates over a million-fold, despite having only a limited repertoire of building blocks and available set of chemical functional groups. Ribozyme active sites usually occur at junctions where negatively charged helices come together, and in many cases leverage this strained electrostatic environment to recruit metal ions in solution that can assist in catalysis. Similar strategies have been implicated in related artificially engineered DNA enzymes. Herein, we apply Poisson–Boltzmann, 3D-RISM, and molecular simulations to study a set of metal-dependent small self-cleaving ribozymes (hammerhead, pistol, and Varkud satellite) as well as an artificially engineered DNAzyme (8–17) to examine electrostatic features and their relation to the recruitment of monovalent and divalent metal ions important for activity. We examine several fundamental roles for these ions that include: (1) structural integrity of the catalytically active state, (2) pKa tuning of residues involved in acid–base catalysis, and (3) direct electrostatic stabilization of the transition state via Lewis acid catalysis. Taken together, these examples demonstrate how RNA electrostatics orchestrates the site-specific and territorial binding of metal ions to play important roles in catalysis.

DNA polymerases have evolved to feature a highly conserved activity across the tree of life: formation of, without exception, internucleotidyl O–P linkages. Can this linkage selectivity be overcome by design to produce xenonucleic acids? Here, we report that the structure-guided redesign of an archaeal DNA polymerase, 9°N, exhibits a new activity undetectable in the wild-type enzyme: catalyzing the formation of internucleotidyl N–P linkages using 3′-NH2-ddNTPs. Replacing a metal-binding aspartate in the 9°N active site with asparagine was key to the emergence of this unnatural enzyme activity. MD simulations provided insights into how a single substitution enhances the productive positioning of a 3′-amino nucleophile in the active site. Further remodeling of the protein–nucleic acid interface in the finger subdomain yielded a quadruple-mutant variant (9°N-NRQS) displaying DNA-dependent NP-DNA polymerase activity. In addition, the engineered promiscuity of 9°N-NRQS was leveraged for one-pot synthesis of DNA─NP-DNA copolymers. This work sheds light on the molecular basis of substrate fidelity and latent promiscuity in enzymes.

We present a fast, accurate, and robust approach for determination of free energy profiles and kinetic isotope effects for RNA 2′-O-transphosphorylation reactions with inclusion of nuclear quantum effects. We apply a deep potential range correction (DPRc) for combined quantum mechanical/molecular mechanical (QM/MM) simulations of reactions in the condensed phase. The method uses the second-order density-functional tight-binding method (DFTB2) as a fast, approximate base QM model. The DPRc model modifies the DFTB2 QM interactions and applies short-range corrections to the QM/MM interactions to reproduce ab initio DFT (PBE0/6-31G*) QM/MM energies and forces. The DPRc thus enables both QM and QM/MM interactions to be tuned to high accuracy, and the QM/MM corrections are designed to smoothly vanish at a specified cutoff boundary (6 Å in the present work). The computational speed-up afforded by the QM/MM+DPRc model enables free energy profiles to be calculated that include rigorous long-range QM/MM interactions under periodic boundary conditions and nuclear quantum effects through a path integral approach using a new interface between the AMBER and i-PI software. The approach is demonstrated through the calculation of free energy profiles of a native RNA cleavage model reaction and reactions involving thio-substitutions, which are important experimental probes of the mechanism. The DFTB2+DPRc QM/MM free energy surfaces agree very closely with the PBE0/6-31G* QM/MM results, and it is vastly superior to the DFTB2 QM/MM surfaces with and without weighted thermodynamic perturbation corrections. 18O and 34S primary kinetic isotope effects are compared, and the influence of nuclear quantum effects on the free energy profiles is examined.

The pistol ribozyme (Psr) is one among the most recently discovered classes of small nucleolytic ribozymes that catalyze site-specific RNA self-cleavage through 2′-O-transphosphorylation. The Psr contains a conserved guanine (G40) that in crystal structures is in a position suggesting it plays the role of the general base to abstract a proton from the nucleophile to activate the reaction. Although some functional data is consistent with this mechanistic role, a notable exception is 2-aminopurine (2AP) substitution which has no effect on the rate, unlike similar substitutions across other so-called “G + M” and “G + A” ribozyme classes. Herein we postulate that an alternate conserved guanine, G42, is the primary general base, and provide evidence from molecular simulations that the active site of Psr can undergo local refolding into a structure that is consistent with the common “L-platform/L-scaffold” architecture identified in G + M and G + A ribozyme classes with Psr currently the notable exception. We summarize the key currently available experimental data and present new classical and combined quantum mechanical/molecular mechanical simulation results that collectively suggest a new hypothesis. We hypothesize that there are two available catalytic pathways supported by different conformational states connected by a local refolding of the active site: (1) a primary pathway with an active site architecture aligned with the L-platform/L-scaffold framework where G42 acts as a general base, and (2) a secondary pathway with the crystallographic active site architecture where G40 acts as a general base. We go on to make several experimentally testable predictions, and suggest specific experiments that would ultimately bring closure to the mystery as to “who stole the proton in the pistol ribozyme?”.

We report a new online suite of tools that enables inquiry-based active-learning activities to develop students’ representational competence about atomic orbitals. Orbital Explorer is a Web site for the visualization and interactive investigation of atomic orbital properties. Orbital Explorer contains two integrated tools, namely, Atomic Orbital Explorer, which enables one to visualize and interrogate individual atomic orbitals, and Orbital RDF Comparison, which enables one to make a more detailed quantitative comparison of orbital energies and properties of orbital radial distribution functions (RDFs). In addition, we present an original chemistry educational gamification design, BingOrbital, constructed in a format resembling Bingo (American version). The game aims to reinforce the recognition of atomic orbitals based on the RDF and three-dimensional isosurface and has been applied as an engaging retrieval practice tool. A companion set of example activities that use the Orbital Explorer and BingOrbital game have been presented in another article.

Atomic orbitals represent an essential construct used to develop chemical bonding models, upon which other more advanced chemistry topics are built. In this article, we share a series of active-learning activities and a gamified approach to develop students’ representational competence about atomic orbitals and to engage students in learning the properties of atomic orbitals. These properties are essential for understanding an array of fundamental concepts such as penetration and shielding, relationships such as periodic trends, and models used to describe chemical bonding. The activities employ an inquiry-based approach to engage students in exploring the relationship between atomic orbitals’ spatial properties and quantum numbers. The activities guide students to collect data to verify periodic trends and construct electronic configurations. The activities utilize Orbital Explorer Web site for visualization, comparison, and analysis of atomic orbitals. The Orbital Explorer Web site is described in a related Technology Report. The activities and the game are suitable to be conducted in both in-person and remote-teaching settings.

Complete understanding of most, if not all chemical processes requires at its very core the knowledge of the underlying free-energy change. In computer-aided drug design, for instance, such processes as binding of a drug to a protein or its spontaneous partitioning across the cell membrane cannot be predicted reliably without considering how the associated free energy varies. Owing to relentless theoretical developments, which have benefited from ever-growing computational resources, free-energy calculations leaning on statistical-mechanics simulations are now part of the arsenal of robust and well-characterized modeling tools. However, as will be explained below and touched upon throughout the chapters of this book, it is still challenging to obtain accurate and reliable free-energy predictions for biomolecules due to the many nuances in the system setup and the unknown unknowns such as whether a given simulation is globally converged or only locally converged, perhaps in an incorrect free-energy basin. In all but the simplest cases, free-energy simulations still require experts in the field to prepare the system, run the calculations, and analyze the results in order to obtain robust predictions that can be confidently used to make decisions in drug discovery campaigns.

Recent concurrent advances in methodology development, computer hardware and simulation software has transformed our ability to make practical, quantitative predictions of relative ligand binding affinities to guide rational drug design. In the past, these calculations have been hampered by the lack of affordable software with highly efficient implementations of state-of-the-art methods on specialized hardware such as graphical processing units, combined with the paucity of available workflows to streamline throughput for real-world industry applications. Herein we discuss recent methodology development, GPU-accelerated implementation, and workflow creation for alchemical free energy simulation methods in the AMBER Drug Discovery Boost (AMBER-DD Boost) package available as a patch to AMBER20. Among the methodological advances are 1) new methods for the treatment of softcore potentials that overcome long standing end-point catastrophe and softcore imbalance problems and enable single-step alchemical transformations between ligands, 2) new adaptive enhanced sampling methods in the ”alchemical” (or ” λ”) dimension to accelerate convergence and obtain high precision ligand binding affinity predictions, 3) robust network-wide analysis methods that include cycle closure and reference constraints and restraints, and 4) practical workflows that enable streamlined calculations on large datasets to be performed. Benchmark calculations on various systems demonstrate that these tools deliver an outstanding combination of accuracy and performance, resulting in reliable high-throughput binding affinity predictions at affordable cost.

We develop a new Deep Potential - Range Correction (DPRc) machine learning potential for combined quantum mechanical/molecular mechanical (QM/MM) simulations of chemical reactions in the condensed phase. The new range correction enables short-ranged QM/MM interactions to be tuned for higher accuracy, and the correction smoothly vanishes within a specified cutoff. We further develop an active learning procedure for robust neural network training. We test the DPRc model and training procedure against a series of 6 non-enzymatic phosphoryl transfer reactions in solution that are important in mechanistic studies of RNA-cleaving enzymes. Specifically, we apply DPRc corrections to a base QM model and test its ability to reproduce free energy profiles generated from a target QM model. We perform comparisons using the MNDO/d and DFTB2 semiempirical models because they produce free energy profiles which differ significantly from each other, thereby providing us a rigorous stress test for the DPRc model and training procedure. The comparisons show that accurate reproduction of the free energy profiles requires correction of the QM/MM interactions out to 6 Å. We further find that the model's initial training benefits from generating data from temperature replica exchange simulations and including high-temperature configurations into the fitting procedure so the resulting models are trained to properly avoid high-energy regions. A single DPRc model was trained to reproduce 4 different reactions and yielded good agreement with the free energy profiles made from the target QM/MM simulations. The DPRc model was further demonstrated to be transferable to 2D free energy surfaces and 1D free energy profiles that were not explicitly considered in the training. Examination of the computational performance of the DPRc model showed that it was fairly slow when run on CPUs, but was sped up almost 100-fold when using an NVIDIA V100 GPUs, resulting in almost negligible overhead. The new DPRc model and training procedure provide a potentially powerful new tool for the creation of next-generation QM/MM potentials for a wide spectrum of free energy applications ranging from drug discovery to enzyme design.

Alchemical free energy methods, such as free energy perturbation (FEP) and thermodynamic integration (TI), become increasingly popular and crucial for drug design and discovery. However, the system preparation of alchemical free energy simulation is an error-prone, time-consuming, and tedious process for a large number of ligands. To address this issue, we have recently presented CHARMM-GUI Free Energy Calculator that can provide input and postprocessing scripts for NAMD and GENESIS FEP molecular dynamics systems. In this work, we extended three submodules of Free Energy Calculator to work with the full suite of GPU-accelerated alchemical free energy methods and tools in AMBER, including input and postprocessing scripts. The BACE1 (β-secretase 1) benchmark set was used to validate the AMBER-TI simulation systems and scripts generated by Free Energy Calculator. The overall results of relatively large and diverse systems are almost equivalent with different protocols (unified and split) and with different timesteps (1, 2, and 4 fs), with R2 > 0.9. More importantly, the average free energy differences between two protocols are small and reliable with four independent runs, with a mean unsigned error (MUE) below 0.4 kcal/mol. Running at least four independent runs for each pair with AMBER20 (and FF19SB/GAFF2.1/OPC force fields), we obtained a MUE of 0.99 kcal/mol and root-mean-square error of 1.31 kcal/mol for 58 alchemical transformations in comparison with experimental data. In addition, a set of ligands for T4-lysozyme was used to further validate our free energy calculation protocol whose results are close to experimental data (within 1 kcal/mol). In summary, Free Energy Calculator provides a user-friendly web-based tool to generate the AMBER-TI system and input files for high-throughput binding free energy calculations with access to the full set of GPU-accelerated alchemical free energy, enhanced sampling, and analysis methods in AMBER.

We redevelop the variational free energy profile (vFEP) method using a cardinal B-spline basis to extend the method for analyzing free energy surfaces (FESs) involving three or more reaction coordinates. We also implemented software for evaluating high-dimensional profiles based on the multistate Bennett acceptance ratio (MBAR) method which constructs an unbiased probability density from global reweighting of the observed samples. The MBAR method takes advantage of a fast algorithm for solving the unbinned weighted histogram (UWHAM)/MBAR equations which replaces the solution of simultaneous equations with a nonlinear optimization of a convex function. We make use of cardinal B-splines and multiquadric radial basis functions to obtain smooth, differentiable MBAR profiles in arbitrary high dimensions. The cardinal B-spline vFEP and MBAR methods are compared using three example systems that examine 1D, 2D, and 3D profiles. Both methods are found to be useful and produce nearly indistinguishable results. The vFEP method is found to be 150 times faster than MBAR when applied to periodic 2D profiles, but the MBAR method is 4.5 times faster than vFEP when evaluating unbounded 3D profiles. In agreement with previous comparisons, we find the vFEP method produces superior FESs when the overlap between umbrella window simulations decreases. Finally, the associative reaction mechanism of hammerhead ribozyme is characterized using 3D, 4D, and 6D profiles, and the higher-dimensional profiles are found to have smaller reaction barriers by as much as 1.5 kcal/mol. The methods presented here have been implemented into the FE-ToolKit software package along with new methods for network-wide free energy analysis in drug discovery.

Proton-coupled electron transfer reactions play critical roles in many aspects of sensory phototransduction. In the case of flavoprotein light sensors, reductive quenching of flavin excited states initiates chemical and conformational changes that ultimately transmit light signals to downstream targets. These reactions generally require neighboring aromatic residues and proton-donating side chains for rapid and coordinated electron and proton transfer to flavin. Although photoreduction of flavoproteins can produce either the anionic (ASQ) or neutral semiquinone (NSQ), the factors that favor one over the other are not well understood. Here we employ a biologically active variant of the light-oxygen-voltage (LOV) domain protein VVD devoid of the adduct-forming Cys residue (VVD-III) to probe the mechanism of flavin photoreduction and protonation. A series of isosteric and conservative residue replacements studied by rate measurements, fluorescence quantum yields, FTIR difference spectroscopy, and molecular dynamics simulations indicate that tyrosine residues facilitate charge recombination reactions that limit sustained flavin reduction, whereas methionine residues facilitate radical propagation and quenching and also gate solvent access for flavin protonation. Replacement of a single surface Met residue with Leu favors formation of the ASQ over the NSQ and desensitizes photoreduction to oxidants. In contrast, increasing site hydrophilicity by Gln substitution promotes rapid NSQ formation and weakens the influence of the redox environment. Overall, the photoreactivity of VVD-III can be understood in terms of redundant electron donors, internal hole quenching, and coupled proton transfer reactions that all depend upon protein conformation, dynamics, and solvent penetration.

We describe an efficient method for the simultaneous solution of all free energies within a relative binding free-energy (RBFE) network with cycle closure and experimental/reference constraint conditions using Bennett Acceptance Ratio (BAR) and Multistate BAR (MBAR) analysis. Rather than solving the BAR or MBAR equations for each transformation independently, the simultaneous solution of all transformations are obtained by performing a constrained minimization of a global objective function. The nonlinear optimization of the objective function is subjected to affine linear constraints that couple the free energies between the network edges. The constraints are used to enforce the closure of thermodynamic cycles within the RBFE network, and to enforce an additional set of linear constraint conditions demonstrated here to be subsets of (1 or 2) experimental values. We describe details of the practical implementation of the network BAR/MBAR procedure, including use of generalized coordinates in the minimization of the free-energy objective function, propagation of bootstrap errors from those coordinates, and performance and memory optimization. In some cases it is found that use of restraints in the optimization is more practical than use of generalized coordinates for enforcing constraint conditions. The fast BARnet and MBARnet methods are used to analyze the RBFEs of six prototypical protein–ligand systems, and it is shown that enforcement of cycle closure conditions reduces the error in the predictions only modestly, and further reduction in errors can be achieved when one or two experimental RBFEs are included in the optimization procedure. These methods have been implemented into FE-ToolKit, a new free-energy analysis toolkit. The BARnet/MBARnet framework presented here opens the door to new, more efficient and robust free-energy analysis with enhanced predictive capability for drug discovery applications.

As the conversation in higher education shifts from diversity to inclusion, the attrition rates of students in the STEM fields continue to be a point of discussion. Combined with the demand for expansion in the STEM workforce, various retention reforms have been proposed, implemented, and in some cases integrated into policy following evidence of success. Still, new findings, technological advances, and socio-cultural shifts inevitably necessitate an ongoing investigation as to how students approach learning. Among other factors, students who enter college without effective study skills are at much greater risk of being unsuccessful in their coursework. In order to construct an equitable learning environment, a mechanism must be developed to provide underprepared students with access to resources or interventions designed to refine the skills they need to be successful in the course. Early, reliable assessments can provide predictions of individual student outcomes in order to guide the development and implementation of such targeted interventions. In the present study, a model is developed to predict students’ odds of success based on their study approaches, as measured by their responses to twelve survey items from an existing instrument used in the Chemistry Education Research literature designed to measure students’ deep and surface learning approaches. The model’s prediction specificity ranges from 66.5% to 86.9% by semester. Two distinct sets of lower-performing students are identified in the data: those who align predominantly with surface approaches to learning versus those who indicate using both deep and surface approaches to learning. This supports the idea of a tailored approach to interventions, rather than a one-size-fits-all solution. Results from this instrument were correlated to students’ reported study methods and beliefs.

Predicting protein-ligand binding affinities and the associated thermodynamics of biomolecular recognition is a primary objective of structure-based drug design. Alchemical free energy simulations offer a highly accurate and computationally efficient route to achieving this goal. While the AMBER molecular dynamics package has successfully been used for alchemical free energy simulations in academic research groups for decades, widespread impact in industrial drug discovery settings has been minimal due to previous limitations within the AMBER alchemical code, coupled with challenges in system setup and postprocessing workflows. Through a close academia-industry collaboration we have addressed many of the previous limitations with an aim to improve accuracy, efficiency and robustness of alchemical binding free energy simulations in industrial drug discovery applications. Here, we highlight some of the recent advances in AMBER20 with a focus on alchemical binding free energy (BFE) calculations, which are less computationally intensive than alternative binding free energy methods where full binding/unbinding paths are explored. In addition to scientific and technical advances in AMBER20, we also describe the essential practical aspects associated with running relative alchemical BFE calculations along with recommendations for best practices, highlighting the importance not only of the alchemical simulation code, but also the auxiliary functionalities and expertise required to obtain accurate and reliable results. This work is intended to provide a contemporary overview of the scientific, technical, and practical issues associated with running relative BFE simulations in AMBER20, with a focus on real-world drug discovery applications.

The Teaching Internship is a credit-bearing program composed of undergraduate near peer instructors (teaching interns, or TIs) that offers supplemental assistance for students in the General Chemistry courses. With fellow undergraduates serving as a role model and student–faculty liaison, the benefits of near peer instruction have been well-documented. Because TIs develop a dual role of student and instructor over time, they afford a unique opportunity to explore the middle area of the expert/novice spectrum. Identifying the most influential components of the TI role may allow practitioners to implement these components in other ways for different groups of students. The present work provides a description of the TI model and uses a mixed-methods approach to analyze how the peer leadership role impacted the TIs’ attitudes about learning chemistry. Quantitative results show that TIs do hold predominantly expert-like learning attitudes compared to the General Chemistry population from which they are selected; however, evidence of novice thinking is still observed in some areas. This survey data was then used to inform a qualitative approach. Further analysis indicated that TIs’ responses on survey items were context-dependent, and that peer leadership experiences were associated with expert learning attitudes and appear to be influential in the development of these attitudes. These findings suggest that these factors should be taken into account when drawing general conclusions from survey results.

Progress in the development of GPU-accelerated free energy simulation software has enabled practical applications on complex biological systems and fueled efforts to develop more accurate and robust predictive methods. In particular, this work re-examines concerted (a.k.a., one-step or unified) alchemical transformations commonly used in the prediction of hydration and relative binding free energies (RBFEs). We first classify several known challenges in these calculations into three categories: endpoint catastrophes, particle collapse, and large gradient-jumps. While endpoint catastrophes have long been addressed using softcore potentials, the remaining two problems occur much more sporadically and can result in either numerical instability (i.e. complete failure of a simulation) or inconsistent estimation (i.e. stochastic convergence to an incorrect result). The particle collapse problem stems from an imbalance in short-range electrostatic and repulsive interactions and can, in principle, be solved by appropriately balancing the respective softcore parameters. However, the large gradient-jump problem itself arises from the sensitivity of the free energy to large values of the softcore parameters, as might be used in trying to solve the particle collapse issue. Often no satisfactory compromise exists with the existing softcore potential form. As a framework for solving these problems, we developed a new family of smoothstep softcore (SSC) potentials motivated by an analysis of the derivatives along the alchemical path. The smoothstep polynomials generalize the monomial functions that are used in most implementations and provide an additional path-dependent smoothing parameter. The effectiveness of this approach is demonstrated on simple, yet pathological cases that illustrate the three problems outlined. With appropriate parameter selection we find that a second-order SSC(2) potential does at least as well as the conventional approach and provides a vast improvement in terms of consistency across all cases. Lastly, we compare the concerted SSC(2) approach against the gold-standard stepwise (a.k.a., decoupled or multi-step) scheme over a large set of RBFE calculations as might be encountered in drug discovery.

With advancements in GPU-accelerated free energy methods, it is now possible to obtain sufficiently high precision in free energy calculations to rigorously stress test implementations for consistency, reproducibility and reliability. Herein we  rovide high precision validation tests that examine alchemical transformations of a small molecule data set that has been used elsewhere to examine the reproducibility of free energy calculations across different molecular simulation software  packages. We demonstrate that the most recent, updated AMBER18 provides consistent free energy results in both the gas phase and in solution. We first show, in the context of thermodynamic integration (TI), that results are invariant with  respect to “split” (e.g., stepwise decharge-vdW-recharge) versus “unified” protocols. This brought to light a subtle inconsistency in previous versions of AMBER that was traced to the improper treatment of 1-4 vdW and electrostatic interactions involving atoms across the softcore boundary. We illustrate that, under the assumption that the ensembles produced by different legs of the alchemical transformation between molecules “A” and “B” in the gas phase and aqueous phase are very small, the  inconsistency on the relative hydration free energy is minimal. However, for general cases where the ensembles are shown to be substantially different, these errors can be large. Finally, we demonstrate that results for relative hydration free energy simulations are independent of TI or multistate Bennett’s acceptance ratio (MBAR) analysis, invariant to the specific choice of the softcore region, and agree with results derived from absolute hydration free energy values.

Accurate prediction of the absolute or relative protein−ligand binding affinity is one of the major tasks in computer-aided drug design projects, especially in the stage of lead optimization. In principle, the alchemical free energy (AFE) methods such as thermodynamic integration (TI) or free-energy perturbation (FEP) can fulfill this task, but in practice, a lot of hurdles prevent them from being routinely applied in daily drug design projects, such as the demanding computing resources, slow computing processes, unavailable or inaccurate force field parameters, and difficult and unfriendly setting up and post-analysis procedures. In this study, we have exploited practical protocols of applying the CPU (central processing unit)-TI and newly developed GPU (graphic processing unit)-TI modules and other tools in the AMBER software package, combined with ff14SB/GAFF1.8 force fields, to conduct efficient and accurate AFE calculations on protein−ligand binding free energies. We have tested 134 protein−ligand complexes in total for four target proteins (BACE, CDK2, MCL1, and PTP1B) and obtained overall comparable performance with the commercial Schrodinger FEP+ program (Wang et al. J. Am. Chem. Soc. 2015, 137, 2695−2703). The achieved accuracy fits within the requirements for computations to generate effective guidance for experimental work in drug lead optimization, and the needed wall time is short enough for practical application. Our verified protocol provides a practical solution for routine AFE calculations in real drug design projects.

The Varkud satellite ribozyme catalyses site-specific RNA cleavage and ligation, and serves as an important model system to understand RNA catalysis. Here, we combine stereospecific phosphorothioate substitution, precision nucleobase mutation and linear free-energy relationship measurements with molecular dynamics, molecular solvation theory and ab initio quantum mechanical/molecular mechanical free-energy simulations to gain insight into the catalysis. Through this confluence of theory and experiment, we unify the existing body of structural and functional data to unveil the catalytic mechanism in unprecedented detail, including the degree of proton transfer in the transition state. Further, we provide evidence for a critical Mg2+ in the active site that interacts with the scissile phosphate and anchors the general base guanine in position for nucleophile activation. This novel role for Mg2+ adds to the diversity of known catalytic RNA strategies and unifies functional features observed in the Varkud satellite, hairpin and hammerhead ribozyme classes.

We develop an L-platform/L-scaffold framework we hypothesize may serve as a blueprint to facilitate site-specific RNA-cleaving nucleic acid enzyme design. Building on the L-platform motif originally described by Suslov and coworkers, we identify new critical scaffolding elements required to anchor a conserved general base guanine ("L-anchor") and bind functionally important metal ions at the active site ("L-pocket"). Molecular simulations, together with a broad range of experimental structural and functional data, connect the L-platform/L-scaffold elements to necessary and sufficient conditions for catalytic activity. We demonstrate that the L-platform/L-scaffold framework is common to 5 of the 9 currently known naturally occurring ribozyme classes (Twr, HPr, VSr, HHr, Psr), and intriguingly from a design perspective, the framework also appears in an artificially engineered DNAzyme (8-17dz). The flexibility of the L-platform/L-scaffold framework is illustrated on these systems, highlighting modularity and trends in the variety of known general acid moieties that are supported. These trends give rise to two distinct catalytic paradigms, building on the classifications proposed by Wilson and coworkers and named for the implicated general base and acid. The "G+A" paradigm (Twr, HPr, VSr) exclusively utilizes nucleobase residues for chemistry, and the "G+M⁺" paradigm (HHr, 8-17dz, Psr) involves structuring of the "L-pocket" metal ion binding site for recruitment of a divalent metal ion that plays an active role in the chemical steps of the reaction. Finally, the modularity of the L-platform/L-scaffold framework is illustrated in the VS ribozyme where the "L-pocket" assumes the functional role of the "L-anchor" element, highlighting a distinct mechanism, but one that is functionally linked with the hammerhead ribozyme.

An unique catalytic strategy was recently reported for the glmS ribozyme [Bingaman et al., Nat. Chem. Biol. 2017, 13, 439−445] that involves promotion of productive hydrogen bonding of the O2′ nucleophile to facilitate its activation. We provide broad evidence of this strategy in the hammerhead, pistol, and VS ribozymes and 8-17 DNAzyme, enabled by a functionally important divalent metal ion that interacts with the scissile phosphate and disrupts nonproductive competitive hydrogen bonding with the O2′ nucleophile. This strategy, designated tertiary gamma (3°γ) catalysis, illustrates an additional role for divalent ions in ribozyme catalysis.

We use the PBE0/6-31G* density functional method to perform ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations under periodic boundary conditions with rigorous electrostatics using the ambient potential composite Ewald method in order to test the convergence of MM → QM/MM free energy corrections for the prediction of 17 small-molecule solvation free energies and eight ligand binding free energies to T4 lysozyme. The “indirect” thermodynamic cycle for calculating free energies is used to explore whether a series of reference potentials improve the statistical quality of the predictions. Specifically, we construct a series of reference potentials that optimize a molecular mechanical (MM) force field’s parameters to reproduce the ab initio QM/MM forces from a QM/MM simulation. The optimizations form a systematic progression of successively expanded parameters that include bond, angle, dihedral, and charge parameters. For each reference potential, we calculate benchmark quality reference values for the MM → QM/MM correction by performing the mixed MM and QM/MM Hamiltonians at 11 intermediate states, each for 200 ps. We then compare forward and reverse application of Zwanzig’s relation, thermodynamic integration (TI), and Bennett’s acceptance ratio (BAR) methods as a function of reference potential, simulation time, and the number of simulated intermediate states. We find that Zwanzig’s equation is inadequate unless a large number of intermediate states are explicitly simulated. The TI and BAR mean signed errors are very small even when only the end-state simulations are considered, and the standard deviations of the TI and BAR errors are decreased by choosing a reference potential that optimizes the bond and angle parameters. We find a robust approach for the data sets of fairly rigid molecules considered here is to use bond + angle reference potential together with the end-state-only BAR analysis. This requires QM/MM simulations to be performed in order to generate reference data to parametrize the bond + angle reference potential, and then this same simulation serves a dual purpose as the full QM/MM end state. The convergence of the results with respect to time suggests that computational resources may be used more efficiently by running multiple simulations for no more than 50 ps, rather than running one long simulation.

We perform molecular dynamics simulations, based on recent crystallographic data, on the 8–17 DNAzyme at four states along the reaction pathway to determine the dynamical ensemble for the active state and transition state mimic in solution. A striking finding is the diverse roles played by Na+ and Pb2+ ions in the electrostatically strained active site that impact all four fundamental catalytic strategies, and share commonality with some features recently inferred for naturally occurring hammerhead and pistol ribozymes. The active site Pb2+ ion helps to stabilize in-line nucleophilic attack, provides direct electrostatic transition state stabilization, and facilitates leaving group departure. A conserved guanine residue is positioned to act as the general base, and is assisted by a bridging Na+ ion that tunes the pKa and facilitates in-line fitness. The present work provides insight into how DNA molecules are able to solve the RNA-cleavage problem, and establishes functional relationships between the mechanism of these engineered DNA enzymes with their naturally evolved RNA counterparts. This adds valuable information to our growing body of knowledge on general mechanisms of phosphoryl transfer reactions catalyzed by RNA, proteins and DNA.

The pistol ribozyme (Psr) is among the most recently discovered RNA enzymes, and has been the subject of experiments aimed at elucidating mechanism. Recent biochemical studies have revealed exciting clues about catalytic interactions in the active site not apparent from available crystallographic data. The present work unifies the interpretation of the existing body of structural and functional data on Psr by providing a dynamical model for the catalytically active state in solution from molecular simulation. Our results suggest that a catalytic Mg2+ ion makes inner-sphere contact with G33:N7 and outer-sphere coordination to the pro-Rp of the scissile phosphate, promoting electrostatic stabilization of the dianionic transition state and neutralization of the developing charge of the leaving group through a metal-coordinated water molecule that is made more acidic by a hydrogen bond donated from the 2’OH of P32. This model is consistent with experimental activity-pH and mutagenesis data, including sensitivity to G33(7cG) and phosphorothioate substitution/metal ion rescue. The model suggests several experimentally testable predictions including stereospecific thio substitutions and G42X (X=xanthine) mutations, some of which have appeared and been validated during the review of this work. Further, the model identifies striking similarities of Psr to the hammerhead ribozyme (HHr), including: similar global fold, organization of secondary structure around an active site 3-way junction, catalytic metal ion binding mode, and guanine general base. However, the specific binding mode and role of the Mg2+ ion, as well as a conserved 2’-OH in the active site, are inter-related but subtly different between the ribozymes.

A predictive understanding of the mechanisms of RNA cleavage is important for the design of emerging technology built from biological and synthetic molecules that have promise for new biochemical and medicinal applications. Over the past 15 years, RNA cleavage reactions involving 2′-O-transphosphorylation have been discussed using a simplified framework introduced by Breaker that consists of four fundamental catalytic strategies (designated α, β, γ, and δ) that contribute to rate enhancement. As more detailed mechanistic data emerge, there is need for the framework to evolve and keep pace. We develop an ontology for discussion of strategies of enzymes that catalyze RNA cleavage via 2′-O-transphosphorylation that stratifies Breaker’s framework into primary (1°), secondary (2°), and tertiary (3°) contributions to enable more precise interpretation of mechanism in the context of structure and bonding. Further, we point out instances where atomic-level changes give rise to changes in more than one catalytic contribution, a phenomenon we refer to as “functional blurring”. We hope that this ontology will help clarify our conversations and pave the path forward toward a consensus view of these fundamental and fascinating mechanisms. The insight gained will deepen our understanding of RNA cleavage reactions catalyzed by natural protein and RNA enzymes, as well as aid in the design of new engineered DNA and synthetic enzymes.

The catalytic properties of RNA have been a subject of fascination and intense research since their discovery over 30 years ago. Very recently, several classes of nucleolytic ribozymes have emerged and been characterized structurally. Among these, the twister ribozyme has been center-stage and a topic of debate about its architecture and mechanism owing to conflicting interpretations of different crystal structures and in some cases conflicting interpretations of the same functional data. In the present work, we attempt to clean up the mechanistic "debris" generated by twister ribozymes using a comprehensive computational RNA enzymology approach aimed to provide a unified interpretation of existing structural and functional data. Simulations in the crystalline environment and in solution provide insight into the origins of observed differences in crystal structures and coalesce on a common active site architecture and dynamical ensemble in solution. We use GPU-accelerated free energy methods with enhanced sampling to ascertain microscopic nucleobase pKa values of the implicated general acid and base, from which predicted activity–pH profiles can be compared directly with experiments. Next, ab initio quantum mechanical/molecular mechanical (QM/MM) simulations with full dynamic solvation under periodic boundary conditions are used to determine mechanistic pathways through multidimensional free energy landscapes for the reaction. We then characterize the rate-controlling transition state and make predictions about kinetic isotope effects and linear free energy relations. Computational mutagenesis is performed to explain the origin of rate effects caused by chemical modifications and to make experimentally testable predictions. Finally, we provide evidence that helps to resolve conflicting issues related to the role of metal ions in catalysis. Throughout each stage, we highlight how a conserved L-platform structural motif, together with a key L-anchor residue, forms the characteristic active site scaffold enabling each of the catalytic strategies to come together for not only the twister ribozyme but also the majority of the known small nucleolytic ribozyme classes.

Crystal simulations provide useful tools, along with solution simulations, to test nucleic acid force fields, but should be interpreted with care owing to the difficulty of establishing the environmental conditions needed to reproduce experimental crystal packing. These challenges underscore the need to construct proper protocols for carrying out crystal simulations and analyzing results to identify the origin of deviations from crystallographic data. Toward this end, we introduce a novel framework for B-factor decomposition into additive intramolecular, rotational, and translational atomic fluctuation components and partitioning of each of these components into individual asymmetric unit and lattice contributions. We apply the framework to a benchmark set of A-DNA, Z-DNA, and B-DNA double helix systems of various chain lengths. Overall, the intramolecular deviations from the crystal were quite small (≤1.0 Å), suggesting high accuracy of the force field, whereas crystal packing was not well reproduced. The present work establishes a framework to conduct and analyze crystal simulations that ultimately take on issues of crystal packing and can provide insight into nucleic acid force fields.

It is generally accepted that hydrogen tunneling enhances both primary and secondary H/D kinetic isotope effects (KIEs) over what would be expected under the assumptions of classical barrier transition. Previous studies have exclusively shown that the effects of tunneling upon primary H/D KIEs have been much larger than those observed for secondary H/D KIEs. Here we present a series of experimental H/D KIE results associated with the Chugaev elimination of methyl xanthate derived from ß-phenylethanol over the temperature range of 180 to 290 °C. Intramolecular H/D KIEs computed according to classical transition state theory (TST) are markedly overestimated relative to experimentally measured values. Experimental intermolecular H/D KIEs and direct dynamic calculations based on canonical variational transition state theory (CVT) with small-curvature tunneling correction (SCT) reveal that this result is largely the consequence of extraordinary tunneling enhancement of the secondary H/D KIE. This unexpected behavior is examined in the context of other similar hydrogen transfer reactions.

Abstract Site binding of ions and water shapes nucleic acids folding, dynamics and biological function, complementing the more diffuse territorial binding. Unlike territorial binding, prediction of site-specific binding to nucleic acids remains an unsolved challenge in computational biophysics. This work presents a new toolset based on the 3D-RISM molecular solvation theory and topological analysis that predicts cation and water site-binding to nucleic acids. 3D-RISM is shown to accurately capture alkali cations and water binding to the central channel, transversal loops and grooves of the Oxytricha nova's telomeres G-quadruplex (Oxy-GQ), in agreement with high resolution crystallographic data. To improve the computed cations occupancy along the Oxy-GQ central channel, it was necessary to refine and validate new cation-oxygen parameters using structural and thermodynamic data available for crown ethers and ion channels. This single set of parameters that describes both localized and delocalized binding to various biological systems is used to gain insight into cation occupancy along the Oxy-GQ channel under various salt conditions. The paper concludes with prospects for extending the method to predict divalent metals binding to nucleic acids. This work advances the forefront of theoretical methods able to provide predictive insight into ion atmosphere effects on nucleic acids function.

Maintaining a proper balance between specific intermolecular interactions and non-specific solvent interactions is of critical importance in molecular simulations, especially when predicting binding affinities or reaction rates in the condensed phase. The most rigorous metric for characterizing solvent affinity are solvation free energies, which correspond to a transfer from the gas phase into solution. Due to the drastic change of the electrostatic environment during this process, it is also a stringent test of polarization response in the model. Here, we employ both the CHARMM fixed charge and polarizable force fields to predict hydration free energies of twelve simple solutes. The resulting classical ensembles are then reweighted to obtain QM/MM hydration free energies using a variety of QM methods, including MP2, Hartree–Fock, density functional methods (BLYP, B3LYP, M06-2X) and semi-empirical methods (OM2 and AM1 ). Our simulations test the compatibility of quantum-mechanical methods with molecular-mechanical water models and solute Lennard–Jones parameters. In all cases, the resulting QM/MM hydration free energies were inferior to purely classical results, with the QM/MM Drude force field predictions being only marginally better than the QM/MM fixed charge results. In addition, the QM/MM results for different quantum methods are highly divergent, with almost inverted trends for polarizable and fixed charge water models. While this does not necessarily imply deficiencies in the QM models themselves, it underscores the need to develop consistent and balanced QM/MM interactions. Both the QM and the MM component of a QM/MM simulation have to match, in order to avoid artifacts due to biased solute–solvent interactions. Finally, we discuss strategies to improve the convergence and efficiency of multi-scale free energy simulations by automatically adapting the molecular-mechanics force field to the target quantum method.

We report progress in graphics processing unit (GPU)-accelerated molecular dynamics and free energy methods in Amber18. Of particular interest is the development of alchemical free energy algorithms, including free energy perturbation and thermodynamic integration methods with support for nonlinear soft-core potential and parameter interpolation transformation pathways. These methods can be used in conjunction with enhanced sampling techniques such as replica exchange, constant-pH molecular dynamics, and new 12–6–4 potentials for metal ions. Additional performance enhancements have been made that enable appreciable speed-up on GPUs relative to the previous software release.

In this work, we employ simple model systems to evaluate the relative performance of two of the most important free energy methods: The Zwanzig equation (also known as "Free energy perturbation") and Bennett’s acceptance ratio method (BAR). Although our examples should be transferable to other kinds of free energy simulations, we focus on applications of multi-scale free energy simulations. Such calculations are especially complex, since they connect two different levels of theory with very different requirements in terms of speed, accuracy, sampling and parallelisability. We try to reconcile all those different factors by developing some simple criteria to guide the early stages of the development of a free energy protocol. This is accomplished by quantifying how many intermediate steps and how many potential energy evaluations are necessary in order to reach a certain level of convergence.

The ability of RNA molecules to selectively and efficiently catalyze complex chemical transformation has vast implications. Our understanding of the mechanisms of RNA catalysis have been greatly advanced by the study of small nucleolytic RNA enzymes, or ribozymes, that have evolved naturally in viruses and living organisms, or artificially through high-throughput in vitro selection techniques [1]. Experimental structural and mechanistic work, along with computational simulations have provided deep insight into the mechanism of these model ribozyme systems. Very recently, there has been a surge in progress in the determination of crystallographic structures of ribozymes [2–6] that have provided a departure point for theoretical investigations that aim to bridge the gap between structural and mechanistic measurements, and provide a detailed dynamical picture of mechanism at atomic-level resolution. In this way, molecular simulations have the potential to unify the interpretations of a broad range of experimental data and establish a consensus view of mechanism [7]. Ultimately, multiscale simulations, together with experiments, afford the tools needed to gain predictive insight into catalysis, including control factors that regulate selectivity and reactivity, that may guide rational design efforts. In the present work, results from multiscale molecular simulations are presented for a series of ribozymes for which crystallographic data has recently become available, including the twister [2], Varkud satallite virus (VS) [3] and pistol [4] ribozymes. These results uncover recurring themes, as well as new twists, in the catalytic strategies taken by ribozymes that are apparent only when broadly analyzing their structure, biochemical characterization and detailed mechanisms predicted by molecular simulations. The interpretations of experimental data afforded by new multiscale molecular simulation results uncover general principles and provide predictive insight into the catalytic mechanisms of nucleolytic ribozymes that may guide rational design efforts.

Recently, 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 C7, or the N3 position of residue C7 itself.

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.

Molecular simulations are powerful tools for providing atomic-level details into complex chemical and physical processes that occur in the condensed phase. For strongly interacting systems where quantum many-body effects are known to play an important role, density-functional methods are often used to provide the model with the potential energy used to drive dynamics. These methods, however, suffer from two major drawbacks. First, they are often too computationally intensive to practically apply to large systems over long time scales, limiting their scope of application. Second, there remain challenges for these models to obtain the necessary level of accuracy for weak non-bonded interactions to obtain quantitative accuracy for a wide range of condensed phase properties. Quantum mechanical force fields (QMFFs) provide a potential solution to both of these limitations. In this review, we address recent advances in the development of QMFFs for condensed phase simulations. In particular, we examine the development of QMFF models using both approximate and ab initio density-functional models, the treatment of short-ranged non-bonded and long-ranged electrostatic interactions, and stability issues in molecular dynamics calculations. Example calculations are provided for crystalline systems, liquid water, and ionic liquids. We conclude with a perspective for emerging challenges and future research directions.

We report the implementation of the thermodynamic integration method on the pmemd module of the AMBER 16 package on GPUs (pmemdGTI). The pmemdGTI code typically delivers over 2 orders of magnitude of speed-up relative to a single CPU core for the calculation of ligand-protein binding affinities with no statistically significant numerical differences and thus provides a powerful new tool for drug discovery applications.

The hammerhead ribozyme is a well-studied nucleolytic ribozyme that catalyzes the self-cleavage of the RNA phosphodiester backbone. Despite experimental and theoretical efforts, key questions remain about details of the mechanism with regard to the activation of the nucleophile by the putative general base guanine (G12). Straightforward interpretation of the measured activity–pH data implies the pKa value of the N1 position in the G12 nucleobase is significantly shifted by the ribozyme environment. Recent crystallographic and biochemical work has identified pH-dependent divalent metal ion binding at the N7/O6 position of G12, leading to the hypothesis that this binding mode could induce a pKa shift of G12 toward neutrality. We present computational results that support this hypothesis and provide a model that unifies the interpretation of available structural and biochemical data and paints a detailed mechanistic picture of the general base step of the reaction. Experimentally testable predictions are made for mutational and rescue effects on G12, which will give further insights into the catalytic mechanism. These results contribute to our growing knowledge of the potential roles of divalent metal ions in RNA catalysis.

Modern crystal structure refinement programs rely on geometry restraints to overcome the challenge of a low data-to-parameter ratio. While the classical Engh and Huber restraints work well for standard amino-acid residues, the chemical complexity of small-molecule ligands presents a particular challenge. Most current approaches either limit ligand restraints to those that can be readily described in the Crystallographic Information File (CIF) format, thus sacrificing chemical flexibility and energetic accuracy, or they employ protocols that substantially lengthen the refinement time, potentially hindering rapid automated refinement workflows. PHENIX-AFITT refinement uses a full molecular-mechanics force field for user-selected small-molecule ligands during refinement, eliminating the potentially difficult problem of finding or generating high-quality geometry restraints. It is fully integrated with a standard refinement protocol and requires practically no additional steps from the user, making it ideal for high-throughput workflows. PHENIX-AFITT refinements also handle multiple ligands in a single model, alternate conformations and covalently bound ligands. Here, the results of combining AFITT and the PHENIX software suite on a data set of 189 protein-ligand PDB structures are presented. Refinements using PHENIX-AFITT significantly reduce ligand conformational energy and lead to improved geometries without detriment to the fit to the experimental data. For the data presented, PHENIX-AFITT refinements result in more chemically accurate models for small-molecule ligands.

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.

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.

Recent advances in biomolecular modeling have emphasized the importance of inclusion of explicit electronic polarizabilty, and a description of electrostatic interactions that includes atomic multipoles; however, these additional levels of treatment necessarily increase a model’s computational cost. Ultimately, the decision as to whether inclusion of these more rigorous levels are justified rests on the degree to which they impact the specific application areas of interest, balanced with the overhead of their computational cost. The purpose of this book is to stimulate the exchange of effective Many-Body Effects and Electrostatics in Biomolecules strategies used to describe many-body effects and electrostatics across the quantum, classical, and coarse-grained modeling regimes

We present results from molecular dynamics simulations and free energy calculations of the twister ribozyme at different stages along the reaction path in order to gain insight into its mechanism. 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, there currently lacks 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 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.

Solvent D2O and 18O kinetic isotope effects on RNA 2'-O-transphosphorylation catalyzed by Zn2+ demonstrate an altered transition state relative to specific base catalysis. A recent model from DFT calculations involving inner sphere coordination to the non-bridging and leaving group oxygens is consistent with the data.

Biological catalysis hinges on the precise structural integrity of an active site that binds and transforms its substrates and meeting this requirement presents a unique challenge for RNA enzymes. Functional RNAs, including ribozymes, fold into their active conformations within rugged energy landscapes that often contain misfolded conformers. Here we uncover and characterize one such “off-pathway” species within an active site after overall folding of the ribozyme is complete. The Tetrahymena group I ribozyme (E) catalyzes cleavage of an oligonucleotide substrate (S) by an exogenous guanosine (G) cofactor. We tested whether specific catalytic interactions with G are present in the preceding E•S•G and E•G ground-state complexes. We monitored interactions with G via the effects of 2'- and 3'-deoxy (–H) and -amino (–NH2) substitutions on G binding. These and prior results reveal that G is bound in an inactive configuration within E•G, with the nucleophilic 3'-OH making a nonproductive interaction with an active site metal ion termed MA and with the adjacent 2'-OH making no interaction. Upon S binding, a rearrangement occurs that allows both –OH groups to contact a different active site metal ion, termed MC, to make what are likely to be their catalytic interactions. The reactive phosphoryl group on S promotes this change, presumably by repositioning the metal ions with respect to G. This conformational transition demonstrates local rearrangements within an otherwise folded RNA, underscoring RNA's difficulty in specifying a unique conformation and highlighting Nature's potential to use local transitions of RNA in complex function

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.

A fully quantum mechanical force field (QMFF) based on a modified “divide-and-conquer” (mDC) framework is applied to a series of molecular simulation applications, using a generalized Particle Mesh Ewald method extended to multipolar charge densities. Simulation results are presented for three example applications: liquid water, p-nitrophenylphosphate reactivity in solution, and crystalline N,N-dimethylglycine. Simulations of liquid water using a parametrized mDC model are compared to TIP3P and TIP4P/Ew water models and experiment. The mDC model is shown to be superior for cluster binding energies and generally comparable for bulk properties. Examination of the dissociative pathway for dephosphorylation of p-nitrophenylphosphate shows that the mDC method evaluated with the DFTB3/3OB and DFTB3/OPhyd semiempirical models bracket the experimental barrier, whereas DFTB2 and AM1/d-PhoT QM/MM simulations exhibit deficiencies in the barriers, the latter for which is related, in part, to the anomalous underestimation of the p-nitrophenylate leaving group pKa. Simulations of crystalline N,N-dimethylglycine are performed and the overall structure and atomic fluctuations are compared with the experiment and the general AMBER force field (GAFF). The QMFF, which was not parametrized for this application, was shown to be in better agreement with crystallographic data than GAFF. Our simulations highlight some of the application areas that may benefit from using new QMFFs, and they demonstrate progress toward the development of accurate QMFFs using the recently developed mDC framework.

The well-studied mechanism of ribonuclease A is believed to involve concerted general acid-base catalysis by two histidine residues, His12 and His119. The basic features of this mechanism are often cited to explain rate enhancement by both protein and RNA enzymes that catalyze RNA 2'-O-transphosphorylation. Recent kinetic isotope effect analyses and computational studies are providing a more chemically detailed description of the mechanism of RNase A and the rate limiting transition state. Overall, the results support an asynchronous mechanism for both solution and ribonuclease catalyzed reactions in which breakdown of a transient dianoinic phosphorane intermediate by 5'OP bond cleavage is rate limiting. Relative to non-enzymatic reactions catalyzed by specific base, a smaller KIE on the 5'O leaving group and a less negative ßLG are observed for RNase A catalysis. Quantum mechanical calculations consistent with these data support a model in which electrostatic and H-bonding interactions with the non-bridging oxygens and proton transfer from His119 render departure of the 5'O less advanced and stabilize charge buildup in the transition state. Both experiment and computation indicate advanced 2'OP bond formation in the rate limiting transition state. However, this feature makes it difficult to resolve the chemical steps involved in 2'O activation. Thus, modeling the transition state for RNase A catalysis underscores those elements of its chemical mechanism that are well resolved, as well as highlighting those where ambiguity remains. This article is part of a Special Issue entitled: Enzyme Transition States from Theory and Experiment.

Understanding how enzymes catalyze specific reactions has yielded important insights into biology, provided the necessary context for the design of novel catalysts, and has facilitated the development of enzyme inhibitors for therapeutic purposes. The contributions of Professor W. Wallace Cleland helped to forge the modern framework for steady state enzyme kinetics, and together with his colleagues he made pioneering contributions to establishing heavy atom isotope effect analyses as a powerful approach for investigating the chemical mechanisms of enzyme-catalyzed reactions. Professor Cleland collaborated with and trained numerous talented enzymologists and generations of scientists. Several facets of W.W. Cleland's contributions to science are illustrated and reflected in this special issue that focuses on using kinetic isotope effects to understand chemical and enzymatic reaction mechanisms.

Experimental analysis of kinetic isotope effects represents an extremely powerful approach for gaining information about the transition state structure of complex reactions not available through other methodologies. The implementation of this approach to the study of nucleic acid chemistry requires the synthesis of nucleobases and nucleotides enriched for heavy isotopes at specific positions. In this review, we highlight current approaches to the synthesis of nucleic acids enriched site specifically for heavy oxygen and nitrogen and their application in heavy atom isotope effect studies. This article is part of a special issue titled: Enzyme Transition States from Theory and Experiment.

Divalent metal ions, due to their ability to stabilize high concentrations of negative charge, are important for RNA folding and catalysis. Detailed models derived from the structures and kinetics of enzymes and from computational simulations have been developed. However, in most cases the specific catalytic modes involving metal ions and their mechanistic roles and effects on transition state structures remain controversial. Valuable information about the nature of the transition state is provided by measurement of kinetic isotope effects (KIEs). However, KIEs reflect changes in all bond vibrational modes that differ between the ground state and transition state. QM calculations are therefore essential for developing structural models of the transition state and evaluating mechanistic alternatives. Herein, we present computational models for Zn2+ binding to RNA 2'O-transphosphorylation reaction models that aid in the interpretation of KIE experiments. Different Zn2+ binding modes produce distinct KIE signatures, and one binding mode involving two zinc ions is in close agreement with KIEs measured for non-enzymatic catalysis by Zn2+ aquo ions alone. Interestingly, the KIE signatures in this specific model are also very close to those in RNase A catalysis. These results allow a quantitative connection to be made between experimental KIE measurements and transition state structure and bonding, and provide insight into RNA 2'O-ransphosphorylation reactions catalyzed by metal ions and enzymes. This article is part of a Special Issue entitled: Enzyme Transition States from Theory and Experiment.

The composition of the ion atmosphere surrounding nucleic acids affects their folding, condensation and binding to other molecules. It is thus of fundamental importance to gain predictive insight into the formation of the ion atmosphere and thermodynamic consequences when varying ionic conditions. An early step toward this goal is to benchmark computational models against quantitative experimental measurements. Herein, we test the ability of the three dimensional reference interaction site model (3D-RISM) to reproduce preferential interaction parameters determined from ion counting (IC) experiments for mixed alkali chlorides and dsDNA. Calculations agree well with experiment with slight deviations for salt concentrations >200 mM and capture the observed trend where the extent of cation accumulation around the DNA varies inversely with its ionic size. Ion distributions indicate that the smaller, more competitive cations accumulate to a greater extent near the phosphoryl groups, penetrating deeper into the grooves. In accord with experiment, calculated IC profiles do not vary with sequence, although the predicted ion distributions in the grooves are sequence and ion size dependent. Calculations on other nucleic acid conformations predict that the variation in linear charge density has a minor effect on the extent of cation competition.

The hepatitis delta virus ribozyme is an efficient catalyst of RNA 2'-O-transphosphorylation and has emerged as a key experimental system for identifying and characterizing fundamental features of RNA catalysis. Recent structural and biochemical data have led to a proposed mechanistic model whereby an active site Mg2+ ion facilitates deprotonation of the O2' nucleophile, and a protonated cytosine residue (C75) acts as an acid to donate a proton to the O5' leaving group as noted in a previous study. This model assumes that the active site Mg2+ ion forms an inner-sphere coordination with the O2' nucleophile and a nonbridging oxygen of the scissile phosphate. These contacts, however, are not fully resolved in the crystal structure, and biochemical data are not able to unambiguously exclude other mechanistic models. In order to explore the feasibility of this model, we exhaustively mapped the free energy surfaces with different active site ion occupancies via quantum mechanical/molecular mechanical (QM/MM) simulations. We further incorporate a three-dimensional reference interaction site model for the solvated ion atmosphere that allows these calculations to consider not only the rate associated with the chemical steps, but also the probability of observing the system in the presumed active state with the Mg2+ ion bound. The QM/MM results predict that a pathway involving metal-assisted nucleophile activation is feasible based on the rate-controlling transition state barrier departing from the presumed metal-bound active state. However, QM/MM results for a similar pathway in the absence of Mg2+ are not consistent with experimental data, suggesting that a structural model in which the crystallographically determined Mg2+ is simply replaced with Na+ is likely incorrect. It should be emphasized, however, that these results hinge upon the assumption of the validity of the presumed Mg2+-bound starting state, which has not yet been definitively verified experimentally, nor explored in depth computationally. Thus, further experimental and theoretical study is needed such that a consensus view of the catalytic mechanism emerges.

L1 ligase (L1L) molecular switch is an in vitro optimized synthetic allosteric ribozyme that catalyzes the regioselective formation of a 5'-to-3' phosphodiester bond, a reaction for which there is no known naturally occurring RNA catalyst. L1L serves as a proof of principle that RNA can catalyze a critical reaction for prebiotic RNA self-replication according to the RNA world hypothesis. L1L crystal structure captures two distinct conformations that differ by a reorientation of one of the stems by around 80 Å and are presumed to correspond to the active and inactive state, respectively. It is of great interest to understand the nature of these two states in solution and the pathway for their interconversion. In this study, we use explicit solvent molecular simulation together with a novel enhanced sampling method that utilizes concepts from network theory to map out the conformational transition between active and inactive states of L1L. We find that the overall switching mechanism can be described as a three-state/two-step process. The first step involves a large-amplitude swing that reorients stem C. The second step involves the allosteric activation of the catalytic site through distant contacts with stem C. Using a conformational space network representation of the L1L switch transition, it is shown that the connection between the three states follows different topographical patterns: the stem C swing step passes through a narrow region of the conformational space network, whereas the allosteric activation step covers a much wider region and a more diverse set of pathways through the network.

1,2,3,4-Diepoxybutane (DEB) is considered the ultimate carcinogenic metabolite of 1,3-butadiene, an important industrial chemical and environmental pollutant present in urban air. Although it preferentially modifies guanine within DNA, DEB induces a large number of A ? T transversions, suggesting that it forms strongly mispairing lesions at adenine nucleobases. We now report the discovery of three potentially mispairing exocyclic adenine lesions of DEB: N6,N6-(2,3-dihydroxybutan-1,4-diyl)-2'-deoxyadenosine (compound 2), 1,N6-(2-hydroxy-3-hydroxymethylpropan-1,3-diyl)-2'-deoxyadenosine (compound 3), and 1,N6-(1-hydroxymethyl-2-hydroxypropan-1,3-diyl)-2'-deoxyadenosine (compound 4). The structures and stereochemistry of the novel DEB-dA adducts were determined by a combination of UV and NMR spectroscopy, tandem mass spectrometry, and independent synthesis. We found that synthetic N6-(2-hydroxy-3,4-epoxybut-1-yl)-2'-deoxyadenosine (compound 1) representing the product of N6-adenine alkylation by DEB spontaneously cyclizes to form 3 under aqueous conditions or 2 under anhydrous conditions in the presence of an organic base. Compound 3 can be interconverted with 4 by a reversible unimolecular pericyclic reaction favoring 4 as a more thermodynamically stable product. Both 3 and 4 are present in double stranded DNA treated with DEB in vitro and in liver DNA of laboratory mice exposed to 1,3-butadiene by inhalation. We propose that in DNA under physiological conditions, DEB alkylates the N-1 position of adenine in DNA to form N1-(2-hydroxy-3,4-epoxybut-1-yl)-adenine adducts, which undergo an SN2-type intramolecular nucleophilic substitution and rearrangement to give 3 (minor) and 4 (major). Formation of exocyclic DEB-adenine lesions following exposure to 1,3-butadiene provides a possible mechanism of mutagenesis at the A:T base pairs.

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983. © 2009 Wiley Periodicals, Inc.J Comput Chem, 2009.

The present study employs density-functional electronic structure methods to investigate the effect of chemical modification at the C5 position of cytosine. A series of experimentally motivated chemical modifications are considered, including alkyl, halogen, aromatic, fused ring, and strong s and p withdrawing functional groups. The effect of these modifications on cytosine geometry, electronic structure, proton affinities, gas phase basicities, cytosine–guanine base pair hydrogen bond network and corresponding nucleophilicity at guanine are examined. Ultimately, these results play a part in dissecting the effect of endogenous cytosine methylation on the reactivity of neighboring guanine toward carcinogens and DNA alkylating agents.