Absolute binding free energy calculations using molecular dynamics simulations with restraining potentials.
ABSTRACT: The absolute (standard) binding free energy of eight FK506-related ligands to FKBP12 is calculated using free energy perturbation molecular dynamics (FEP/MD) simulations with explicit solvent. A number of features are implemented to improve the accuracy and enhance the convergence of the calculations. First, the absolute binding free energy is decomposed into sequential steps during which the ligand-surrounding interactions as well as various biasing potentials restraining the translation, orientation, and conformation of the ligand are turned "on" and "off." Second, sampling of the ligand conformation is enforced by a restraining potential based on the root mean-square deviation relative to the bound state conformation. The effect of all the restraining potentials is rigorously unbiased, and it is shown explicitly that the final results are independent of all artificial restraints. Third, the repulsive and dispersive free energy contribution arising from the Lennard-Jones interactions of the ligand with its surrounding (protein and solvent) is calculated using the Weeks-Chandler-Andersen separation. This separation also improves convergence of the FEP/MD calculations. Fourth, to decrease the computational cost, only a small number of atoms in the vicinity of the binding site are simulated explicitly, while all the influence of the remaining atoms is incorporated implicitly using the generalized solvent boundary potential (GSBP) method. With GSBP, the size of the simulated FKBP12/ligand systems is significantly reduced, from approximately 25,000 to 2500. The computations are very efficient and the statistical error is small ( approximately 1 kcal/mol). The calculated binding free energies are generally in good agreement with available experimental data and previous calculations (within approximately 2 kcal/mol). The present results indicate that a strategy based on FEP/MD simulations of a reduced GSBP atomic model sampled with conformational, translational, and orientational restraining potentials can be computationally inexpensive and accurate.
Project description:Advanced free energy perturbation molecular dynamics (FEP/MD) simulation methods are available to accurately calculate absolute binding free energies of protein-ligand complexes. However, these methods rely on several sophisticated command scripts implementing various biasing energy restraints to enhance the convergence of the FEP/MD calculations, which must all be handled properly to yield correct results. Here, we present a user-friendly Web interface, CHARMM-GUI Ligand Binder ( http://www.charmm-gui.org/input/gbinding ), to provide standardized CHARMM input files for calculations of absolute binding free energies using the FEP/MD simulations. A number of features are implemented to conveniently set up the FEP/MD simulations in highly customizable manners, thereby permitting an accelerated throughput of this important class of computations while decreasing the possibility of human errors. The interface and a series of input files generated by the interface are tested with illustrative calculations of absolute binding free energies of three nonpolar aromatic ligands to the L99A mutant of T4 lysozyme and three FK506-related ligands to FKBP12. Statistical errors within individual calculations are found to be small (~1 kcal/mol), and the calculated binding free energies generally agree well with the experimental measurements and the previous computational studies (within ~2 kcal/mol). Therefore, CHARMM-GUI Ligand Binder provides a convenient and reliable way to set up the ligand binding free energy calculations and can be applicable to pharmaceutically important protein-ligand systems.
Project description:An increasing number of studies have reported computations of the standard (absolute) binding free energy of small ligands to proteins using molecular dynamics (MD) simulations and explicit solvent molecules that are in good agreement with experiments. This encouraging progress suggests that physics-based approaches hold the promise of making important contributions to the process of drug discovery and optimization in the near future. Two types of approaches are principally used to compute binding free energies with MD simulations. The most widely known is the alchemical double decoupling method, in which the interaction of the ligand with its surroundings are progressively switched off. It is also possible to use a potential of mean force (PMF) method, in which the ligand is physically separated from the protein receptor. For both of these computational approaches, restraining potentials may be activated and released during the simulation for sampling efficiently the changes in translational, rotational, and conformational freedom of the ligand and protein upon binding. Because such restraining potentials add bias to the simulations, it is important that their effects be rigorously removed to yield a binding free energy that is properly unbiased with respect to the standard state. A review of recent results is presented, and differences in computational methods are discussed. Examples of computations with T4-lysozyme mutants, FKBP12, SH2 domain, and cytochrome P450 are discussed and compared. Remaining difficulties and challenges are highlighted.
Project description:Gleevec, a well-known cancer therapeutic agent, is an effective inhibitor of several tyrosine kinases, including Abl and c-Kit, but displays less potency to inhibit closely homologous tyrosine kinases, such as Lck and c-Src. Because many structural features of the binding site are highly conserved in these homologous kinases, the molecular determinants responsible for the binding specificity of Gleevec remain poorly understood. To address this issue, free energy perturbation molecular dynamics (FEP/MD) simulations with explicit solvent was used to compute the binding affinity of Gleevec to Abl, c-Kit, Lck, and c-Src. The results of the FEP/MD calculations are in good agreement with experiments, enabling a detailed and quantitative dissection of the absolute binding free energy in terms of various thermodynamic contributions affecting the binding specificity of Gleevec to the kinases. Dominant binding free energy contributions arises from the van der Waals dispersive interaction, compensating about two-thirds of the unfavorable free energy penalty associated with the loss of translational, rotational, and conformational freedom of the ligand upon binding. In contrast, the contributions from electrostatic and repulsive interactions nearly cancel out due to solvent effects. Furthermore, the calculations show the importance of the conformation of the kinase activation loop. Among the kinases examined, Abl provides the most favorable binding environment for Gleevec via optimal protein-ligand interactions and a small free energy cost for loss of the translational, rotational, and conformational freedom upon ligand binding. The FEP/MD calculations additionally reveal that Lck and c-Src provide similar nonbinding interactions with the bound-Gleevec, but the former pays less entropic penalty for the ligand losing its translational, rotational, and conformational motions to bind, examining the empirically observed differential binding affinities of Gleevec between the two Src-family kinases.
Project description:The hypothetical scanning molecular dynamics (HSMD) method is used here for calculating the absolute free energy of binding, ?A(0) of the complex of the protein FKBP12 with the ligand SB2 (also denoted L8) - a system that has been studied previously for comparing the performance of different methods. Our preliminary study suggests that considering long-range electrostatics is imperative even for a hydrophobic ligand such as L8. Therefore the system is modeled by the AMBER force field using Particle Mesh Ewald (PME). HSMD consists of three stages applied to both the ligand-solvent and ligand-protein systems. (1) A small set of system configurations (frames) is extracted from an MD trajectory. (2) The entropy of the ligand in each frame is calculated by a reconstruction procedure. (3) The contribution of water and protein to ?A(0) is calculated for each frame by gradually increasing the ligand-environment interactions from zero to their full value using thermodynamic integration (TI). Unlike the conventional methods, the structure of the ligand is kept fixed during TI, and HSMD is thus free from the end-point problem encountered with the double annihilation method (DAM); therefore, the need for applying restraints is avoided. Furthermore, unlike the conventional methods, the entropy of the ligand and water is obtained directly as a byproduct of the simulation. In this paper, in addition to the difference in the internal entropies of the ligand in the two environments, we calculate for the first time the external entropy of the ligand, which provides a measure for the size of the active site. We obtain ?A(0) = -10.7 ±1.0 as compared to the experimental values -10.9 and -10.6 kcal/mol. However, a protein/water system treated by periodic boundary conditions grows significantly with increasing protein size and the computation of ?A(0) would become expensive by all methods. Therefore, we also apply HSMD to FKBP12-L8 described by the GSBP/SSBP model of Roux's group (implemented in the software CHARMM) where only part of the protein and water around the active site are considered and long-range electrostatic effects are taken into account. For comparison this model was also treated by the double decoupling method (DDM). The two methods have led to comparable results for ?A(0) which are somewhat lower than the experimental value. The ligand was found to be more confined in the active site described by GSBP/SSBP than by PME where its entropy in solvent is larger than in the active site by 1.7 and by 5.5 kcal/mol, respectively.
Project description:Molecular docking is widely used to obtain binding modes and binding affinities of a molecule to a given target protein. Despite considerable efforts, however, prediction of both properties by docking remains challenging mainly due to protein's structural flexibility and inaccuracy of scoring functions. Here, an integrated approach has been developed to improve the accuracy of binding mode and affinity prediction and tested for small molecule MDM2 and MDMX antagonists. In this approach, initial candidate models selected from docking are subjected to equilibration MD simulations to further filter the models. Free energy perturbation molecular dynamics (FEP/MD) simulations are then applied to the filtered ligand models to enhance the ability in predicting the near-native ligand conformation. The calculated binding free energies for MDM2 complexes are overestimated compared to experimental measurements mainly due to the difficulties in sampling highly flexible apo-MDM2. Nonetheless, the FEP/MD binding free energy calculations are more promising for discriminating binders from nonbinders than docking scores. In particular, the comparison between the MDM2 and MDMX results suggests that apo-MDMX has lower flexibility than apo-MDM2. In addition, the FEP/MD calculations provide detailed information on the different energetic contributions to ligand binding, leading to a better understanding of the sensitivity and specificity of protein-ligand interactions.
Project description:Large and flexible ligands gain increasing interest in the development of bioactive agents. They challenge the applicability of computational ligand optimization strategies originally developed for small molecules. Free energy perturbation (FEP) is often used for predicting binding affinities of small molecule ligands, however, its use for more complex ligands remains limited. Herein, we report the structure-based design of peptide macrocycles targeting the protein binding site of human adaptor protein 14-3-3. We observe a surprisingly strong dependency of binding affinities on relatively small variations in substituent size. FEP was performed to rationalize observed trends. To account for insufficient convergence of FEP, restrained calculations were performed and complemented with extensive REST MD simulations of the free ligands. These calculations revealed that changes in affinity originate both from altered direct interactions and conformational changes of the free ligand. In addition, MD simulations provided the basis to rationalize unexpected trends in ligand lipophilicity. We also verified the anticipated interaction site and binding mode for one of the high affinity ligands by X-ray crystallography. The introduced fully-atomistic simulation protocol can be used to rationalize the development of structurally complex ligands which will support future ligand maturation efforts.
Project description:Monte Carlo free energy perturbation (MC/FEP) calculations have been applied to compute the relative binding affinities of 17 congeneric pyridazo-pyrimidinone inhibitors of the protein p38? MAP kinase. Overall correlation with experiment was found to be modest when the complexes were hydrated using a traditional procedure with a stored solvent box. Significant improvements in accuracy were obtained when the MC/FEP calculations were repeated using initial solvent distributions optimized by the water placement algorithm JAWS. The results underscore the importance of accurate placement of water molecules in a ligand binding site for the reliable prediction of relative free energies of binding.
Project description:Alchemical free energy simulations have long been utilized to predict free energy changes for binding affinity and solubility of small molecules. However, while the theoretical foundation of these methods is well established, seamlessly handling many of the practical aspects regarding the preparation of the different thermodynamic end states of complex molecular systems and the numerous processing scripts often remains a burden for successful applications. In this work, we present CHARMM-GUI Free Energy Calculator (http://www.charmm-gui.org/input/fec) that provides various alchemical free energy perturbation molecular dynamics (FEP/MD) systems with input and post-processing scripts for NAMD and GENESIS. Four submodules are available: Absolute Ligand Binder (for absolute ligand binding FEP/MD), Relative Ligand Binder (for relative ligand binding FEP/MD), Absolute Ligand Solvator (for absolute ligand solvation FEP/MD), and Relative Ligand Solvator (for relative ligand solvation FEP/MD). Each module is designed to build multiple systems of a set of selected ligands at once for high-throughput FEP/MD simulations. The capability of Free Energy Calculator is illustrated by absolute and relative solvation FEP/MD of a set of ligands and absolute and relative binding FEP/MD of a set of ligands for T4-lysozyme in solution and the adenosine A2A receptor in a membrane. The calculated free energy values are overall consistent with the experimental and published free energy results (within ?1 kcal/mol). We hope that Free Energy Calculator is useful to carry out high-throughput FEP/MD simulations in the field of biomolecular sciences and drug discovery.
Project description:Estimating the correct binding modes of ligands in protein-ligand complexes is crucial not only in the drug discovery process but also for elucidating potential toxicity mechanisms. In the current paper, we propose a computational modeling workflow using the combination of docking, classical molecular dynamics (cMD), accelerated molecular dynamics (aMD) and free-energy perturbation (FEP+ protocol) for identification of possible ligand binding modes. It was applied for investigation of selected perfluorocarboxyl acids (PFCAs) in the PPAR? nuclear receptor. Although both regular and induced fit docking failed to reproduce the experimentally determined binding mode of the ligands when docked into a non-native X-ray structure, cMD and aMD simulations successfully identified the most probable binding conformations. Moreover, multiple binding modes were identified for all of these compounds and the shorter-chain PFCAs continuously moved between a few energetically favorable binding conformations. On the basis of MD predictions of binding conformations, we applied the default and also redesigned FEP+ sampling protocols, which accurately reproduced experimental differences in the binding energies. Thus, the preliminary MD simulations can also provide helpful information about correct setup of the FEP+ calculations. These results show that the PFCA binding modes were accurately predicted and that the FEP+ protocol can be used to estimate free energies of binding of flexible ligands that are not typical druglike compounds. Our in silico workflow revealed the specific ligand-residue interactions within the ligand binding domain and the main characteristics of the PFCAs, and it was concluded that these compounds are week PPAR? partial agonists. This work also suggests a common pipeline for identification of ligand binding modes, ligand-protein dynamics description, and relative free-energy calculations.
Project description:Recent improvements to the free energy perturbation (FEP) calculations, especially FEP+?, established their utility for pharmaceutical lead optimization. Herein, we propose a modified version of the FEP/REST (i.e., replica exchange with solute tempering) sampling protocol, based on detail studies on several targets by probing a large number of perturbations with different sampling schemes. Improved FEP+ binding affinity predictions for regular flexible-loop motions and considerable structural changes can be obtained by extending the prior to REST (pre-REST) sampling time from 0.24?ns/? to 5?ns/? and 2?×?10?ns/?, respectively. With this new protocol, much more precise ??G values of the individual perturbations, including the sign of the transformations and decreased error were obtained. We extended the REST simulations from 5?ns to 8?ns to achieve reasonable free energy convergence. Implementing REST to the entire ligand as opposed to solely the perturbed region, and also some important flexible protein residues (pREST region) in the ligand binding domain (LBD) has considerably improved the FEP+ results in most of the studied cases. Preliminary molecular dynamics (MD) runs were useful for establishing the correct binding mode of the compounds and thus precise alignment for FEP+?. Our improved protocol may further increase the FEP+ accuracy.