Validation of Free Energy Methods in AMBER

Calculated hydration free energies become less accurate with increases in molecular weight

In order for computer-aided drug design to fulfil its long held promise of delivering new medicines faster and cheaper, extensive development and validation work must be done first. This pertains particularly to molecular dynamics force fields where one important aspect-the hydration free energy (HFE) of small molecules-is often insufficiently analyzed. While most benchmarking studies report excellent accuracies of calculated hydration free energies-usually within 2 kcal/mol of experimental values-we find that deeper analysis reveals significant shortcomings. Herein, we report a dependence of HFE prediction errors on ligand molecular weight-the higher the weight, the bigger the prediction error and the higher the probability the calculated result is erroneous by a large amount. We show that in the drug-like molecular weight region, HFE predictions can easily be off by 5 kcal/mol or more. This is likely to be highly problematic in a drug discovery and development setting. We make our HFE results and molecular descriptors freely and fully available in order to encourage deeper analysis of future molecular dynamics results and facilitate development of the next generation of force fields.

Ligand Binding Affinity Prediction for Membrane Proteins with Alchemical Free Energy Calculation Methods

Alchemical relative binding free energy (ΔΔG) calculations have shown high accuracy in predicting ligand binding affinity and have been used as important tools in computer-aided drug discovery and design. However, there has been limited research on the application of ΔΔG methods to membrane proteins despite the fact that these proteins represent a significant proportion of drug targets, play crucial roles in biological processes, and are implicated in numerous diseases. In this study, to predict the binding affinity of ligands to G protein-coupled receptors (GPCRs), we employed two ΔΔG calculation methods: thermodynamic integration (TI) with AMBER and the alchemical transfer method (AToM) with OpenMM. We calculated ΔΔG values for 53 transformations involving four class A GPCRs and evaluated the performance of AMBER-TI and AToM-OpenMM. In addition, we conducted tests using different numbers of windows and varying simulation times to achieve reliable ΔΔG results and to optimize resource utilization. Overall, both AMBER-TI and AToM-OpenMM show good agreement with the experimental data. Our results validate the applicability of AMBER-TI and AToM-OpenMM for optimization of lead compounds targeting membrane proteins.

Software Infrastructure for Next-Generation QM/MM-ΔMLP Force Fields

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.

Binding of Cholesterol to the N-Terminal Domain of the NPC1L1 Transporter: Analysis of the Epimerization-Related Binding Selectivity and Loop Mutations

Cholesterol is a fat-like substance with a pivotal physiological relevance in humans, and its homeostasis is tightly regulated by various cellular processes, including the import in the small intestine and the reabsorption in the biliary ducts by the Niemann-Pick C1 Like 1 (NPC1L1) importer. NPC1L1 can mediate the absorption of a variety of sterols but strikingly exhibits a large sensitivity to cholesterol epimerization. This study examines the molecular basis of the epimerization-related selective binding of cholesterol by combining extended unbiased molecular dynamics simulations of the apo and holo species of the N-terminal domain of wild-type NPC1L1, in conjunction with relative binding free energy, umbrella sampling, and well-tempered metadynamics calculations. The analysis of the results discloses the existence of two distinct binding modes for cholesterol and epi-cholesterol. The former binds deeper in the cavity, forming key hydrogen-bond interactions with Q95, S56, and a water molecule. In contrast, epi-cholesterol is shifted ca. 3 Å to the mouth of the cavity and the transition to the Q95 site is prevented by an energetic barrier of 4.1 kcal·mol-1. Thus, the configuration of the hydroxyl group of cholesterol, together with the presence of a structural water molecule, is a key feature for effective absorption. Finally, whereas these findings may seemingly be challenged by single-point mutations that impair cholesterol transport but have a mild impact on the binding of cholesterol to the Q95 binding site, our results reveal that they have a drastic influence on the conformational landscape of the α8/β7 loop in the apo species compared to the wild-type protein. Overall, the results give support to the functional role played by the α8/β7 loop in regulating the access of ligands to NPC1L1, and hence to interpreting the impact of these mutations on diseases related to disruption of sterol absorption, paving the way to understanding certain physiological dysfunctions.

Modern Alchemical Free Energy Methods for Drug Discovery Explained

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.

ACES: Optimized Alchemically Enhanced Sampling

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.

AMBER Drug Discovery Boost Tools: Automated Workflow for Production Free-Energy Simulation Setup and Analysis (ProFESSA)

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.

In Silico Positional Analogue Scanning with Amber GPU-TI

Positional analogue scanning (PAS) is an accepted strategy for multiparameter lead optimization (MPO) in drug discovery. Small structural changes as introduced by PAS can lead to 10-fold changes in binding potency in ∼10-20% of cases, a significant parameter shift irrespective of other MPO objectives. Sometimes performing a complete PAS is challenging due to resource and time constraints, building block availability, or difficulty in synthesis. Calculating relative binding free energies (RBFEs) for all positions can contribute to prioritizing the most promising analogues for synthesis. We tested a well-established RBFE calculation method, Amber GPU-TI, for 20 positional analogue scans in 14 test systems (cyclin-dependent kinase 8 (CDK8), hepatitis C virus nonstructural protein 5B (HCV NS5B), tankyrase, RAC-α serine/threonine-protein kinase (Akt), phosphodiesterase 1B (PDE1B), orexin/hypocretin receptor type 1 (OX1R), orexin/hypocretin receptor type 2 (OX2R), histone acetyltransferase K (lysine) acetyltransferase 6A (KAT6A), peroxisome proliferator-activated receptor γ (PPARγ), extracellular signal-regulated kinases (ERK1/2), coactivator-associated arginine methyltransferase 1 (PRMT4), α_vβ₆, bromodomain 1 (BD1), human immunodeficiency virus-1 (HIV-1) entry) involving nitrogen, methyl, halogen, methoxy, and hydroxyl scans with at least four analogues per set. Among the 66 analogue positions explored, we found that in 18 cases Amber GPU-TI calculations predicted a more than 10-fold change in potency. In all of these cases, the experimentally observed direction of potency changes agreed with the predictions. In 16 cases, more than 10-fold changes in experimental potency were observed. Again, in all of these cases, Amber GPU-TI predicted the direction of the potency changes correctly. In none of these cases would a decision made for or against synthesis based on a 10-fold change in potency have resulted in missing an important analogue. Therefore, in silico RBFE calculations using Amber GPU-TI can meaningfully contribute to the prioritization of positional analogues before synthesis.

Best practices for constructing, preparing, and evaluating protein-ligand binding affinity benchmarks [Article v0.1]

Free energy calculations are rapidly becoming indispensable in structure-enabled drug discovery programs. As new methods, force fields, and implementations are developed, assessing their expected accuracy on real-world systems (benchmarking) becomes critical to provide users with an assessment of the accuracy expected when these methods are applied within their domain of applicability, and developers with a way to assess the expected impact of new methodologies. These assessments require construction of a benchmark-a set of well-prepared, high quality systems with corresponding experimental measurements designed to ensure the resulting calculations provide a realistic assessment of expected performance when these methods are deployed within their domains of applicability. To date, the community has not yet adopted a common standardized benchmark, and existing benchmark reports suffer from a myriad of issues, including poor data quality, limited statistical power, and statistically deficient analyses, all of which can conspire to produce benchmarks that are poorly predictive of real-world performance. Here, we address these issues by presenting guidelines for (1) curating experimental data to develop meaningful benchmark sets, (2) preparing benchmark inputs according to best practices to facilitate widespread adoption, and (3) analysis of the resulting predictions to enable statistically meaningful comparisons among methods and force fields. We highlight challenges and open questions that remain to be solved in these areas, as well as recommendations for the collection of new datasets that might optimally serve to measure progress as methods become systematically more reliable. Finally, we provide a curated, versioned, open, standardized benchmark set adherent to these standards (PLBenchmarks) and an open source toolkit for implementing standardized best practices assessments (arsenic) for the community to use as a standardized assessment tool. While our main focus is free energy methods based on molecular simulations, these guidelines should prove useful for assessment of the rapidly growing field of machine learning methods for affinity prediction as well.

Relative Binding Free Energy between Chemically Distant Compounds Using a Bidirectional Nonequilibrium Approach

In the context of advanced hit-to-lead drug design based on atomistic molecular dynamics simulations, we propose a dual topology alchemical approach for calculating the relative binding free energy (RBFE) between two chemically distant compounds. The method (termed NE-RBFE) relies on the enhanced sampling of the end-states in bulk and in the bound state via Hamiltonian Replica Exchange, alchemically connected by a series of independent and fast nonequilibrium (NE) simulations. The technique has been implemented in a bidirectional fashion, applying the Crooks theorem to the NE work distributions for RBFE predictions. The dissipation of the NE process, negatively affecting accuracy, has been minimized by introducing a smooth regularization based on shifted electrostatic and Lennard-Jones non bonded potentials. As a challenging testbed, we have applied our method to the calculation of the RBFEs in the recent host-guest SAMPL international contest, featuring a macrocyclic host with guests varying in the net charge, volume, and chemical fingerprints. Closure validation has been successfully verified in cycles involving compounds with disparate Tanimoto coefficients, volume, and net charge. NE-RBFE is specifically tailored for massively parallel facilities and can be used with little or no code modification on most of the popular software packages supporting nonequilibrium alchemical simulations, such as Gromacs, Amber, NAMD, or OpenMM. The proposed methodology bypasses most of the entanglements and limitations of the standard single topology RBFE approach for strictly congeneric series based on free-energy perturbation, such as slowly relaxing cavity water, sampling issues along the alchemical stratification, and the need for highly overlapping molecular fingerprints.

Force Field Parameterization of Actinyl Molecular Cations Using the 12-6-4 Model

In this work, a set of 12-6-4 force fields (FFs) parameters were developed for the actinyl molecular cations, AnO2n+ (n = 1, 2), from uranium to plutonium for classical molecular dynamics (MD) for four water models: TIP3P, SPC/E, OPC3, and TIP4Pew. Such a non-bonded potential model taking into account the induced dipole between the metallic center and the surrounding molecules has shown better performances for various cations than the classic 12-6 non-bonded potentials. The parametrization method proposed elsewhere for metallic cations has been extended to these molecular cations. In contrast to the actinyl 12-6 FFs from the literature, the new models reproduce correctly both solvation and thermodynamic properties, thanks to the inclusion of the induced dipole term (C4). The transferability of such force fields was assessed by performing MD simulations of carbonato actinyl species, which are highly implicated in actinide migration or actinide extraction from seawater. A highly satisfying agreement was found when comparing the EXAFS signals computed from our MD simulation to the experimental ones. The set of FFs developed here opens new possibilities for the study of actinide chemistry.

Alchemical free energy simulations without speed limits. A generic framework to calculate free energy differences independent of the underlying molecular dynamics program

We describe the theory of the so‐called common‐core/serial‐atom‐insertion (CC/SAI) approach to compute alchemical free energy differences and its practical implementation in a Python package called Transformato. CC/SAI is not tied to a specific biomolecular simulation program and does not rely on special purpose code for alchemical transformations. To calculate the alchemical free energy difference between several small molecules, the physical end‐states are mutated into a suitable common core. Since this only requires turning off interactions, the setup of intermediate states is straightforward to automate. Transformato currently supports CHARMM and OpenMM as back ends to carry out the necessary molecular dynamics simulations, as well as post‐processing calculations. We validate the method by computing a series of relative solvation free energy differences.

Precise Binding Free Energy Calculations for Multiple Molecules Using an Optimal Measurement Network of Pairwise Differences

Alchemical binding free energy (BFE) calculations offer an efficient and thermodynamically rigorous approach to in silico binding affinity predictions. As a result of decades of methodological improvements and recent advances in computer technology, alchemical BFE calculations are now widely used in drug discovery research. They help guide the prioritization of candidate drug molecules by predicting their binding affinities for a biomolecular target of interest (and potentially selectivity against undesirable antitargets). Statistical variance associated with such calculations, however, may undermine the reliability of their predictions, introducing uncertainty both in ranking candidate molecules and in benchmarking their predictive accuracy. Here, we present a computational method that substantially improves the statistical precision in BFE calculations for a set of ligands binding to a common receptor by dynamically allocating computational resources to different BFE calculations according to an optimality objective established in a previous work from our group and extended in this work. Our method, termed Network Binding Free Energy (NetBFE), performs adaptive BFE calculations in iterations, re-optimizing the allocations in each iteration based on the statistical variances estimated from previous iterations. Using examples of NetBFE calculations for protein binding of congeneric ligand series, we demonstrate that NetBFE approaches the optimal allocation in a small number (≤5) of iterations and that NetBFE reduces the statistical variance in the BFE estimates by approximately a factor of 2 when compared to a previously published and widely used allocation method at the same total computational cost.

Mechanistic Understanding from Molecular Dynamics in Pharmaceutical Research 2: Lipid Membrane in Drug Design

We review the use of molecular dynamics (MD) simulation as a drug design tool in the context of the role that the lipid membrane can play in drug action, i.e., the interaction between candidate drug molecules and lipid membranes. In the standard "lock and key" paradigm, only the interaction between the drug and a specific active site of a specific protein is considered; the environment in which the drug acts is, from a biophysical perspective, far more complex than this. The possible mechanisms though which a drug can be designed to tinker with physiological processes are significantly broader than merely fitting to a single active site of a single protein. In this paper, we focus on the role of the lipid membrane, arguably the most important element outside the proteins themselves, as a case study. We discuss work that has been carried out, using MD simulation, concerning the transfection of drugs through membranes that act as biological barriers in the path of the drugs, the behavior of drug molecules within membranes, how their collective behavior can affect the structure and properties of the membrane and, finally, the role lipid membranes, to which the vast majority of drug target proteins are associated, can play in mediating the interaction between drug and target protein. This review paper is the second in a two-part series covering MD simulation as a tool in pharmaceutical research; both are designed as pedagogical review papers aimed at both pharmaceutical scientists interested in exploring how the tool of MD simulation can be applied to their research and computational scientists interested in exploring the possibility of a pharmaceutical context for their research.

CHARMM-GUI Free Energy Calculator for Practical Ligand Binding Free Energy Simulations with AMBER

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

Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery

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 because of the 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.

Improved Alchemical Free Energy Calculations with Optimized Smoothstep Softcore Potentials

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 vast improvement in terms of consistency across all cases. Last, we compare the concerted SSC(2) approach against the gold-standard stepwise (a.k.a., decoupled or multistep) scheme over a large set of RBFE calculations as might be encountered in drug discovery.