Alchemical Free Energy Methods Research Articles

Exploring the Limits of the Generalized CHARMM and AMBER Force Fields through Predictions of Hydration Free Energy of Small Molecules.

Accurate force field parameters, potential energy functions, and receptor-ligand models are essential for modeling the solvation and binding of drug-like molecules to a receptor. A large and ever-growing chemical space of medicinally relevant scaffolds has also required these factors, especially force field parameters, to be highly transferable. Generalized force fields such as the CHARMM General Force Field (CGenFF) and the generalized AMBER force field (GAFF) have accomplished this feat along with other contemporaneous ones like OPLS. Here, we analyze the limits in the parametrization of drug-like small molecules by CGenFF and GAFF in terms of the various functional groups represented within them. Specifically, we link the presence of specific functional groups to the error in the absolute hydration free energy of over 600 small molecules, predicted by alchemical free energy methods implemented in the CHARMM program. Our investigation reveals that molecules with (i) a nitro group in CGenFF and GAFF are, respectively, over- or undersolubilized in aqueous medium, (ii) amine groups are undersolubilized more so in CGenFF than in GAFF, and (iii) carboxyl groups are more oversolubilized in GAFF than in CGenFF. We present our analyses of the potential factors underlying these trends. We also showcase the use of a machine-learning-based approach combined with the SHapley Additive exPlanations framework to attribute these trends to specific functional groups, which can be easily adopted to explore the limits of other general force fields.

Application of an alchemical free energy method for the prediction of thermostable DuraPETase variants

Non-equilibrium (NEQ) alchemical free energy calculations are an emerging tool for accurately predicting changes in protein folding free energy resulting from amino acid mutations. In this study, this method in combination with the Rosetta ddg monomer tool was applied to predict more thermostable variants of the polyethylene terephthalate (PET) degrading enzyme DuraPETase. The Rosetta ddg monomer tool efficiently enriched promising mutations prior to more accurate prediction by NEQ alchemical free energy calculations. The relative change in folding free energy of 96 single amino acid mutations was calculated by NEQ alchemical free energy calculation. Experimental validation of ten of the highest scoring variants identified two mutations (DuraPETaseS61M and DuraPETaseS223Y) that increased the melting temperature (Tm) of the enzyme by up to 1 °C. The calculated relative change in folding free energy showed an excellent correlation with experimentally determined Tm resulting in a Pearson’s correlation coefficient of r = − 0.84. Limitations in the prediction of strongly stabilizing mutations were, however, encountered and are discussed. Despite these challenges, this study demonstrates the practical applicability of NEQ alchemical free energy calculations in prospective enzyme engineering projects.Key points• Rosetta ddg monomer enriches stabilizing mutations in a library of DuraPETase variants• NEQ free energy calculations accurately predict changes in Tmof DuraPETase• The DuraPETase variants S223Y, S42M, and S61M have increased TmGraphical

Resolving coupled pH titrations using alchemical free energy calculations.

In a protein, nearby titratable sites can be coupled: the (de)protonation of one may affect the other. The degree of this interaction depends on several factors and can influence the measured . Here, we derive a formalism based on double free energy differences ( ) for quantifying the individual site values of coupled residues. As values can be obtained by means of alchemical free energy calculations, the presented approach allows for a convenient estimation of coupled residue s in practice. We demonstrate that our approach and a previously proposed microscopic formalism, can be combined with alchemical free energy calculations to resolve pH-dependent protein values. Toy models and both, regular and constant-pH molecular dynamics simulations, alongside experimental data, are used to validate this approach. Our results highlight the insights gleaned when coupling and microstate probabilities are analyzed and suggest extensions to more complex enzymatic contexts. Furthermore, we find that naïvely computed values that ignore coupling, can be significantly improved when coupling is accounted for, in some cases reducing the error by half. In short, alchemical free energy methods can resolve the values of both uncoupled and coupled residues.

Ligand-binding affinity prediction for membrane proteins with alchemical free-energy calculation methods

Ligand-binding affinity prediction for membrane proteins with alchemical free-energy calculation methods

Cation-π Interactions Greatly Influence Ion Transportability of the Light-Driven Sodium Pump KR2: A Molecular Dynamics Study.

Krokinobacter eikastus rhodopsin 2 (KR2) is a typical light-driven sodium pump. Although wild-type KR2 exhibits high Na+ selectivity, mutagenesis performed on the residues constituting the entrance enables permeation of K+ and Cs+, while the underlying mechanism remains elusive. This study presents a comprehensive molecular dynamics investigation, including force field optimization, metadynamics, and alchemical free energy methods, to explore the N61L/G263F mutant of KR2, which exhibits transportability for K+ and Cs+. The introduced Phe263 residue can directly promote ion binding at the entrance through cation-π interactions, while the N61L mutation can enhance ion binding at Phe46 by relieving steric hindrance. These results suggest that cation-π interactions may significantly influence the ion transportability and selectivity of KR2, which can provide important insights for protein engineering and the design of artificial ion transporters.

Accurately Predicting Protein pKa Values Using Nonequilibrium Alchemy.

The stability, solubility, and function of a protein depend on both its net charge and the protonation states of its individual residues. pKa is a measure of the tendency for a given residue to (de)protonate at a specific pH. Although pKa values can be resolved experimentally, theory and computation provide a compelling alternative. To this end, we assess the applicability of a nonequilibrium (NEQ) alchemical free energy method to the problem of pKa prediction. On a data set of 144 residues that span 13 proteins, we report an average unsigned error of 0.77 ± 0.09, 0.69 ± 0.09, and 0.52 ± 0.04 pK for aspartate, glutamate, and lysine, respectively. This is comparable to current state-of-the-art predictors and the accuracy recently reached using free energy perturbation methods (e.g., FEP+). Moreover, we demonstrate that our open-source, pmx-based approach can accurately resolve the pKa values of coupled residues and observe a substantial performance disparity associated with the lysine partial charges in Amber14SB/Amber99SB*-ILDN, for which an underused fix already exists.

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.

Alchemical Free Energy Workflows for the Computation of Protein-Ligand Binding Affinities.

Alchemical free energy methods can be used for the efficient computation of relative binding free energies during preclinical drug discovery stages. In recent years, this has been facilitated further by the implementation of workflows that enable non-experts to quickly and consistently set up the required simulations. Given the correct input structures, workflows handle the difficult aspects of setting up perturbations, including consistently defining the perturbable molecule, its atom mapping and topology generation, perturbation network generation, running of the simulations via different sampling methods, and analysis of the results. Different academic and commercial workflows are discussed, including FEW, FESetup, FEPrepare, CHARMM-GUI, Transformato, PMX, QLigFEP, TIES, ProFESSA, PyAutoFEP, BioSimSpace, FEP+, Flare, and Orion. These workflows differ in various aspects, such as mapping algorithms or enhanced sampling methods. Some workflows can accommodate more than one molecular dynamics (MD) engine and use external libraries for tasks. Differences between workflows can present advantages for different use cases, however a lack of interoperability of the workflows' components hinders systematic comparisons.

Uni-GBSA: an open-source and web-based automatic workflow to perform MM/GB(PB)SA calculations for virtual screening.

Binding free energy calculation of a ligand to a protein receptor is a fundamental objective in drug discovery. Molecular mechanics/Generalized-Born (Poisson-Boltzmann) surface area (MM/GB(PB)SA) is one of the most popular methods for binding free energy calculations. It is more accurate than most scoring functions and more computationally efficient than alchemical free energy methods. Several open-source tools for performing MM/GB(PB)SA calculations have been developed, but they have limitations and high entry barriers to users. Here, we introduce Uni-GBSA, a user-friendly automatic workflow to perform MM/GB(PB)SA calculations, which can perform topology preparation, structure optimization, binding free energy calculation and parameter scanning for MM/GB(PB)SA calculations. It also offers a batch mode that evaluates thousands of molecules against one protein target in parallel for efficient application in virtual screening. The default parameters are selected after systematic testing on the PDBBind-2011 refined dataset. In our case studies, Uni-GBSA produced a satisfactory correlation with the experimental binding affinities and outperformed AutoDock Vina in molecular enrichment. Uni-GBSA is available as an open-source package at https://github.com/dptech-corp/Uni-GBSA. It can also be accessed for virtual screening from the Hermite web platform at https://hermite.dp.tech. A free Uni-GBSA web server of a lab version is available at https://labs.dp.tech/projects/uni-gbsa/. This increases user-friendliness because the web server frees users from package installations and provides users with validated workflows for input data and parameter settings, cloud computing resources for efficient job completions, a user-friendly interface and professional support and maintenance.

Understanding the sources of error in MBAR through asymptotic analysis.

Many sampling strategies commonly used in molecular dynamics, such as umbrella sampling and alchemical free energy methods, involve sampling from multiple states. The Multistate Bennett Acceptance Ratio (MBAR) formalism is a widely used way of recombining the resulting data. However, the error of the MBAR estimator is not well-understood: previous error analyses of MBAR assumed independent samples. In this work, we derive a central limit theorem for MBAR estimates in the presence of correlated data, further justifying the use of MBAR in practical applications. Moreover, our central limit theorem yields an estimate of the error that can be decomposed into contributions from the individual Markov chains used to sample the states. This gives additional insight into how sampling in each state affects the overall error. We demonstrate our error estimator on an umbrella sampling calculation of the free energy of isomerization of the alanine dipeptide and an alchemical calculation of the hydration free energy of methane. Our numerical results demonstrate that the time required for the Markov chain to decorrelate in individual states can contribute considerably to the total MBAR error, highlighting the importance of accurately addressing the effect of sample correlation.

High-Throughput Molecular Dynamics-Based Alchemical Free Energy Calculations for Predicting the Binding Free Energy Change Associated with the Selected Omicron Mutations in the Spike Receptor-Binding Domain of SARS-CoV-2.

The ongoing pandemic caused by SARS-CoV-2 has gone through various phases. Since the initial outbreak, the virus has mutated several times, with some lineages showing even stronger infectivity and faster spread than the original virus. Among all the variants, omicron is currently classified as a variant of concern (VOC) by the World Health Organization, as the previously circulating variants have been replaced by it. In this work, we have focused on the mutations observed in omicron sub lineages BA.1, BA.2, BA.4 and BA.5, particularly at the receptor-binding domain (RBD) of the spike protein that is responsible for the interactions with the host ACE2 receptor and binding of antibodies. Studying such mutations is particularly important for understanding the viral infectivity, spread of the disease and for tracking the escape routes of this virus from antibodies. Molecular dynamics (MD) based alchemical free energy calculations have been shown to be very accurate in predicting the free energy change, due to a mutation that could have a deleterious or a stabilizing effect on either the protein itself or its binding affinity to another protein. Here, we investigated the significance of five spike RBD mutations on the stability of the spike protein binding to ACE2 by free energy calculations using high throughput MD simulations. For comparison, we also used conventional MD simulations combined with a Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) based approach, and compared our results with the available experimental data. Overall, the alchemical free energy calculations performed far better than the MM-GBSA approach in predicting the individual impact of the mutations. When considering the experimental variation, the alchemical free energy method was able to produce a relatively accurate prediction for N501Y, the mutant that has previously been reported to increase the binding affinity to hACE2. On the other hand, the other individual mutations seem not to have a significant effect on the spike RBD binding affinity towards hACE2.

Molecular basis of Toxoplasma gondii oryzalin resistance from a novel α-tubulin binding site model

Oryzalin (ORY) is a dinitroaniline derivative that inhibits the microtubule polymerization in plants and parasitic protozoa by selectively binding to the α-tubulin subunit. This herbicidal agent exhibits good antiprotozoal activity against major human parasites, such as Toxoplasma gondii (toxoplasmosis), Leishmania mexicana (leishmaniasis), and Plasmodium falciparum (malaria). Previous chemical mutagenesis assays on T. gondii α-tubulin (TgAT) have identified key mutations that lead to ORY resistance. Herein, we employed alchemical free energy methods and molecular dynamics simulations to determine if the ORY resistance mutations either decrease the TgAT's affinity of the compound or increase the protein stability. Our results here suggest that L136F and V202F mutations significantly decrease the affinity of ORY to TgAT, while T239I and V252L mutations diminish TgAT's flexibility. On the other hand, protein stability predictors determined that R243S mutation reduces TgAT stability due to the loss of its salt bridge interaction with E27. Interestingly, molecular dynamics simulations confirm that the loss of this key interaction leads to ORY binding site closure. Our study provides a better insight into the TgAT-ORY interaction, further supporting our recently proposed ORY-binding site.

Predicting antibiotic resistance in complex protein targets using alchemical free energy methods.

Drug resistant Mycobacterium tuberculosis, which mostly results from single nucleotide polymorphisms in antibiotic target genes, poses a major threat to tuberculosis treatment outcomes. Relative binding free energy (RBFE) calculations can rapidly predict the effects of mutations, but this approach has not been tested on large, complex proteins. We use RBFE calculations to predict the effects of M. tuberculosis RNA polymerase and DNA gyrase mutations on rifampicin and moxifloxacin susceptibility respectively. These mutations encompass a range of amino acid substitutions with known effects and include large steric perturbations and charged moieties. We find that moderate numbers (n = 3–15) of short RBFE calculations can predict resistance in cases where the mutation results in a large change in the binding free energy. We show that the method lacks discrimination in cases with either a small change in energy or that involve charged amino acids, and we investigate how these calculation errors may be decreased.

Large Scale Study of Ligand-Protein Relative Binding Free Energy Calculations: Actionable Predictions from Statistically Robust Protocols.

The accurate and reliable prediction of protein–ligand binding affinities can play a central role in the drug discovery process as well as in personalized medicine. Of considerable importance during lead optimization are the alchemical free energy methods that furnish an estimation of relative binding free energies (RBFE) of similar molecules. Recent advances in these methods have increased their speed, accuracy, and precision. This is evident from the increasing number of retrospective as well as prospective studies employing them. However, such methods still have limited applicability in real-world scenarios due to a number of important yet unresolved issues. Here, we report the findings from a large data set comprising over 500 ligand transformations spanning over 300 ligands binding to a diverse set of 14 different protein targets which furnish statistically robust results on the accuracy, precision, and reproducibility of RBFE calculations. We use ensemble-based methods which are the only way to provide reliable uncertainty quantification given that the underlying molecular dynamics is chaotic. These are implemented using TIES (Thermodynamic Integration with Enhanced Sampling). Results achieve chemical accuracy in all cases. Ensemble simulations also furnish information on the statistical distributions of the free energy calculations which exhibit non-normal behavior. We find that the “enhanced sampling” method known as replica exchange with solute tempering degrades RBFE predictions. We also report definitively on numerous associated alchemical factors including the choice of ligand charge method, flexibility in ligand structure, and the size of the alchemical region including the number of atoms involved in transforming one ligand into another. Our findings provide a key set of recommendations that should be adopted for the reliable application of RBFE methods.

Optimizing Multisite λ-Dynamics Throughput with Charge Renormalization.

With the ability to sample combinations of alchemical perturbations at multiple sites off a small molecule core, multisite λ-dynamics (MSλD) has become an attractive alternative to conventional alchemical free energy methods for exploring large combinatorial chemical spaces. However, current software implementations dictate that combinatorial sampling with MSλD must be performed with a multiple topology model (MTM), which is nontrivial to create by hand, especially for a series of ligand analogues which may have diverse functional groups attached. This work introduces an automated workflow, referred to as msld_py_prep, to assist in the creation of a MTM for use with MSλD. One approach for partitioning partial atomic charges between ligands to create a MTM, called charge renormalization, is also presented and rigorously evaluated. We find that msld_py_prep greatly accelerates the preparation of MSλD ready-to-use files and that charge renormalization can provide a successful approach for MTM generation, as long as bookending calculations are applied to correct small differences introduced by charge renormalization. Charge renormalization also facilitates the use of many different force field parameters with MSλD, broadening the applicability of MSλD for computer-aided drug design.

Alchemical Free Energy Methods Applied to Complexes of the First Bromodomain of BRD4.

Accurate and rapid predictions of the binding affinity of a compound to a target are one of the ultimate goals of computer aided drug design. Alchemical approaches to free energy estimations follow the path from an initial state of the system to the final state through alchemical changes of the energy function during a molecular dynamics simulation. Herein, we explore the accuracy and efficiency of two such techniques: relative free energy perturbation (FEP) and multisite lambda dynamics (MSλD). These are applied to a series of inhibitors for the bromodomain-containing protein 4 (BRD4). We demonstrate a procedure for obtaining accurate relative binding free energies using MSλD when dealing with a change in the net charge of the ligand. This resulted in an impressive comparison with experiment, with an average difference of 0.4 ± 0.4 kcal mol–1. In a benchmarking study for the relative FEP calculations, we found that using 20 lambda windows with 0.5 ns of equilibration and 1 ns of data collection for each window gave the optimal compromise between accuracy and speed. Overall, relative FEP and MSλD predicted binding free energies with comparable accuracy, an average of 0.6 kcal mol–1 for each method. However, MSλD makes predictions for a larger molecular space over a much shorter time scale than relative FEP, with MSλD requiring a factor of 18 times less simulation time for the entire molecule space.

Addressing Intersite Coupling Unlocks Large Combinatorial Chemical Spaces for Alchemical Free Energy Methods.

Alchemical free energy methods are playing a growing role in molecular design, both for computer-aided drug design of small molecules and for computational protein design. Multisite λ dynamics (MSλD) is a uniquely scalable alchemical free energy method that enables more efficient exploration of combinatorial alchemical spaces encountered in molecular design, but simulations have typically been limited to a few hundred ligands or sequences. Here, we focus on coupling between sites to enable scaling to larger alchemical spaces. We first discuss updates to the biasing potentials that facilitate MSλD sampling to include coupling terms and show that this can provide more thorough sampling of alchemical states. We then harness coupling between sites by developing a new free energy estimator based on the Potts models underlying direct coupling analysis, a method for predicting contacts from sequence coevolution, and find it yields more accurate free energies than previous estimators. The sampling requirements of the Potts model estimator scale with the square of the number of sites, a substantial improvement over the exponential scaling of the standard estimator. This opens up exploration of much larger alchemical spaces with MSλD for molecular design.

A combined maximum-caliber and FEP approach to predict how mutations perturb protein stability and folding kinetics

Over the last decade, computational protein design has mainly focused on estimating changes in protein stability using molecular dynamics-based alchemical free energy methods such as free energy perturbation (FEP) and thermodynamic integration (TI). Much less attention has been paid to predicting folding kinetics, due to the complexity of the reaction and difficulty in conformational sampling on slow timescales. Yet, the ability to predict changes to folding kinetics would be highly valuable, both for determining biological function and for computational design.

BLaDE: A Basic Lambda Dynamics Engine for GPU-Accelerated Molecular Dynamics Free Energy Calculations.

There is an accelerating interest in practical applications of alchemical free energy methods to problems in protein design, constant pH simulations, and especially computer-aided drug design. In the present paper, we describe a basic lambda dynamics engine (BLaDE) that enables alchemical free energy simulations, including multisite λ dynamics (MSλD) simulations, on graphical processor units (GPUs). We find that BLaDE is 5 to 8 times faster than the current GPU implementation of MSλD-based free energy calculations in CHARMM. We also demonstrate that BLaDE running standard molecular dynamics attains a performance competitive with and sometimes exceeding that of the highly optimized OpenMM GPU code. BLaDE is available as a standalone program and through an API in CHARMM.

Computing the Liquidus of Binary Monatomic Salt Mixtures with Direct Simulation and Alchemical Free Energy Methods.

We describe and validate a free-energy-based method for computing the liquidus for binary solid-liquid phase diagrams in molecular simulations of monatomic salts. The method is demonstrated by calculating the liquidus for LiCl-KCl and MgCl2-KCl salt mixtures with the polarizable ion model (PIM). The free-energy-based method is cross-validated with direct coexistence simulations. Both techniques show excellent agreement with one another. Though the predictions of the PIM disagree with experiments, we use our free-energy-based approach to decouple the contributions of liquid mixture nonidealities and pure component solid-liquid equilibrium to the phase diagram. In both mixtures, the PIM accurately reproduces the liquid phase nonidealities but fails to predict the liquidus because it does not accurately predict the pure component melting temperature of LiCl or MgCl2.