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.

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.

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.

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.

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 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.

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.

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.

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.

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 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 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.

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.

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.

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.

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.

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

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.