Relative Binding Free Energies Research Articles

D3R grand challenge 2015: Evaluation of protein-ligand pose and affinity predictions.

The Drug Design Data Resource (D3R) ran Grand Challenge 2015 between September 2015 and February 2016. Two targets served as the framework to test community docking and scoring methods: (1) HSP90, donated by AbbVie and the Community Structure Activity Resource (CSAR), and (2) MAP4K4, donated by Genentech. The challenges for both target datasets were conducted in two stages, with the first stage testing pose predictions and the capacity to rank compounds by affinity with minimal structural data; and the second stage testing methods for ranking compounds with knowledge of at least a subset of the ligand-protein poses. An additional sub-challenge provided small groups of chemically similar HSP90 compounds amenable to alchemical calculations of relative binding free energy. Unlike previous blinded Challenges, we did not provide cognate receptors or receptors prepared with hydrogens and likewise did not require a specified crystal structure to be used for pose or affinity prediction in Stage 1. Given the freedom to select from over 200 crystal structures of HSP90 in the PDB, participants employed workflows that tested not only core docking and scoring technologies, but also methods for addressing water-mediated ligand-protein interactions, binding pocket flexibility, and the optimal selection of protein structures for use in docking calculations. Nearly 40 participating groups submitted over 350 prediction sets for Grand Challenge 2015. This overview describes the datasets and the organization of the challenge components, summarizes the results across all submitted predictions, and considers broad conclusions that may be drawn from this collaborative community endeavor.

Insight into drug resistance mechanisms and discovery of potential inhibitors against wild-type and L1196M mutant ALK from FDA-approved drugs.

Anaplastic lymphoma kinase (ALK) plays a crucial role in multiple malignant cancers. It is known as a well-established target for the treatment of ALK-dependent cancers. Even though substantial efforts have been made to develop ALK inhibitors, only crizotinib, ceritinib, and alectinib had been approved by the U.S. Food and Drug Administration for patients with ALK-positive non-small cell lung cancer (NSCLC). The secondary mutations with drug-resistance bring up difficulties to develop effective drugs for ALK-positive cancers. To give a comprehensive understanding of molecular mechanism underlying inhibitor response to ALK tyrosine kinase mutations, we established an accurate assessment for the extensive profile of drug against ALK mutations by means of computational approaches. The molecular mechanics-generalized Born surface area (MM-GBSA) method based on molecular dynamics (MD) simulation was carried out to calculate relative binding free energies for receptor-drug systems. In addition, the structure-based virtual screening was utilized to screen effective inhibitors targeting wild-type ALK and the gatekeeper mutation L1196M from 3180 approved drugs. Finally, the mechanism of drug resistance was discussed, several novel potential wild-type and L1196M mutant ALK inhibitors were successfully identified.

Predicting Binding Affinities for GPCR Ligands Using Free-Energy Perturbation

The rapid growth of structural information for G-protein-coupled receptors (GPCRs) has led to a greater understanding of their structure, function, selectivity, and ligand binding. Although novel ligands have been identified using methods such as virtual screening, computationally driven lead optimization has been possible only in isolated cases because of challenges associated with predicting binding free energies for related compounds. Here, we provide a systematic characterization of the performance of free-energy perturbation (FEP) calculations to predict relative binding free energies of congeneric ligands binding to GPCR targets using a consistent protocol and no adjustable parameters. Using the FEP+ package, first we validated the protocol, which includes a full lipid bilayer and explicit solvent, by predicting the binding affinity for a total of 45 different ligands across four different GPCRs (adenosine A2AAR, β1 adrenergic, CXCR4 chemokine, and δ opioid receptors). Comparison with experimental binding affinity measurements revealed a highly predictive ranking correlation (average spearman ρ = 0.55) and low root-mean-square error (0.80 kcal/mol). Next, we applied FEP+ in a prospective project, where we predicted the affinity of novel, potent adenosine A2A receptor (A2AR) antagonists. Four novel compounds were synthesized and tested in a radioligand displacement assay, yielding affinity values in the nanomolar range. The affinity of two out of the four novel ligands (plus three previously reported compounds) was correctly predicted (within 1 kcal/mol), including one compound with approximately a tenfold increase in affinity compared to the starting compound. Detailed analyses of the simulations underlying the predictions provided insights into the structural basis for the two cases where the affinity was overpredicted. Taken together, these results establish a protocol for systematically applying FEP+ to GPCRs and provide guidelines for identifying potent molecules in drug discovery lead optimization projects.

Resolving the problem of trapped water in binding cavities: prediction of host\u2013guest binding free energies in the SAMPL5 challenge by funnel metadynamics

The funnel metadynamics method enables rigorous calculation of the potential of mean force along an arbitrary binding path and thereby evaluation of the absolute binding free energy. A problem of such physical paths is that the mechanism characterizing the binding process is not always obvious. In particular, it might involve reorganization of the solvent in the binding site, which is not easily captured with a few geometrically defined collective variables that can be used for biasing. In this paper, we propose and test a simple method to resolve this trapped-water problem by dividing the process into an artificial host-desolvation step and an actual binding step. We show that, under certain circumstances, the contribution from the desolvation step can be calculated without introducing further statistical errors. We apply the method to the problem of predicting host–guest binding free energies in the SAMPL5 blind challenge, using two octa-acid hosts and six guest molecules. For one of the hosts, well-converged results are obtained and the prediction of relative binding free energies is the best among all the SAMPL5 submissions. For the other host, which has a narrower binding pocket, the statistical uncertainties are slightly higher; longer simulations would therefore be needed to obtain conclusive results.

Sensitivity in Binding Free Energies Due to Protein Reorganization

Tremendous recent improvements in computer hardware, coupled with advances in sampling techniques and force fields, are now allowing protein-ligand binding free energy calculations to be routinely used to aid pharmaceutical drug discovery projects. However, despite these recent innovations, there are still needs for further improvement in sampling algorithms to more adequately sample protein motion relevant to protein-ligand binding. Here, we report our work identifying and studying such clear and remaining needs in the apolar cavity of T4 lysozyme L99A. In this study, we model recent experimental results that show the progressive opening of the binding pocket in response to a series of homologous ligands.1 Even while using enhanced sampling techniques, we demonstrate that the predicted relative binding free energies (RBFE) are sensitive to the initial protein conformational state. Particularly, we highlight the importance of sufficient sampling of protein conformational changes and demonstrate how inclusion of three key protein residues in the "hot" region of the FEP/REST simulation improves the sampling and resolves this sensitivity, given enough simulation time.

Examining how protein sampling impacts the convergence of relative binding free energy calculations

Examining how protein sampling impacts the convergence of relative binding free energy calculations

Biomolecular Force Field Parameterization via Atoms-in-Molecule Electron Density Partitioning.

Molecular mechanics force fields, which are commonly used in biomolecular modeling and computer-aided drug design, typically treat nonbonded interactions using a limited library of empirical parameters that are developed for small molecules. This approach does not account for polarization in larger molecules or proteins, and the parametrization process is labor-intensive. Using linear-scaling density functional theory and atoms-in-molecule electron density partitioning, environment-specific charges and Lennard-Jones parameters are derived directly from quantum mechanical calculations for use in biomolecular modeling of organic and biomolecular systems. The proposed methods significantly reduce the number of empirical parameters needed to construct molecular mechanics force fields, naturally include polarization effects in charge and Lennard-Jones parameters, and scale well to systems comprised of thousands of atoms, including proteins. The feasibility and benefits of this approach are demonstrated by computing free energies of hydration, properties of pure liquids, and the relative binding free energies of indole and benzofuran to the L99A mutant of T4 lysozyme.

Calculation of Relative Binding Free Energy in the Water-Filled Active Site of Oligopeptide-Binding Protein A.

The periplasmic oligopeptide binding protein A (OppA) represents a well-known example of water-mediated protein-ligand interactions. Here, we perform free-energy calculations for three different ligands binding to OppA, using a thermodynamic integration approach. The tripeptide ligands share a high structural similarity (all have the sequence KXK), but their experimentally-determined binding free energies differ remarkably. Thermodynamic cycles were constructed for the ligands, and simulations conducted in the bound and (freely solvated) unbound states. In the unbound state, it was observed that the difference in conformational freedom between alanine and glycine leads to a surprisingly slow convergence, despite their chemical similarity. This could be overcome by increasing the softness parameter during alchemical transformations. Discrepancies remained in the bound state however, when comparing independent simulations of the three ligands. These difficulties could be traced to a slow relaxation of the water network within the active site. Fluctuations in the number of water molecules residing in the binding cavity occur mostly on a timescale larger than the simulation time along the alchemical path. After extensive simulations, relative binding free energies that were converged to within thermal noise could be obtained, which agree well with available experimental data.

How a Kinase Inhibitor Withstands Gatekeeper Residue Mutations.

Mutations in the gatekeeper residue of kinases have emerged as a key way through which cancer cells develop resistance to treatment. As such, the design of gatekeeper mutation resistant kinase inhibitors is a crucial way forward in increasing the efficacy of a broad range of anticancer drugs. In this work we use atomistic simulations to provide detailed thermodynamic and structural insight into how two inhibitors of cSrc kinase, namely, a commercial drug and type I kinase inhibitor Dasatinib and the type II inhibitor RL45, respectively fail and succeed in being effective against the T338M gatekeeper residue mutation in the kinase binding site. Given the well-known limitations of atomistic simulations in sampling biomolecular systems, we use an enhanced sampling technique called free energy perturbation with replica exchange solute tempering (FEP/REST). Our calculations find that the type I inhibitor Dasatinib binds favorably to the wild type but unfavorably to T338M mutated kinase, while RL45 binds favorably to both. The predicted relative binding free energies are well within 1 kcal/mol accuracy compared to experiments. We find that Dasatinib's impotency against gatekeeper residue mutations arises from a loss of ligand-kinase hydrogen bonding due to T338M mutation and from steric hindrance due to the presence of an inflexible phenyl ring close to the ligand. On the other hand, in the type II binding RL45, the central phenyl ring has very pronounced flexibility. This leads to the inhibitor overcoming effects of steric clashes on mutation and maintaining an electrostatically favorable "edge-to-face" orientation with a neighboring phenylalanine residue. Our work provides useful insight into the mechanisms of mutation resistant kinase inhibitors and demonstrates the usefulness of enhanced sampling techniques in computational drug design.

Virtual substitution scan via single-step free energy perturbation.

With the rapid expansion of our computing power, molecular dynamics (MD) simulations ranging from hundreds of nanoseconds to microseconds or even milliseconds have become increasingly common. The majority of these long trajectories are obtained from plain (vanilla) MD simulations, where no enhanced sampling or free energy calculation method is employed. To promote the "recycling" of these trajectories, we developed the Virtual Substitution Scan (VSS) toolkit as a plugin of the open-source visualization and analysis software VMD. Based on the single-step free energy perturbation (sFEP) method, VSS enables the user to post-process a vanilla MD trajectory for a fast free energy scan of substituting aryl hydrogens by small functional groups. Dihedrals of the functional groups are sampled explicitly in VSS, which improves the performance of the calculation and is found particularly important for certain groups. As a proof-of-concept demonstration, we employ VSS to compute the solvation free energy change upon substituting the hydrogen of a benzene molecule by 12 small functional groups frequently considered in lead optimization. Additionally, VSS is used to compute the relative binding free energy of four selected ligands of the T4 lysozyme. Overall, the computational cost of VSS is only a fraction of the corresponding multi-step FEP (mFEP) calculation, while its results agree reasonably well with those of mFEP, indicating that VSS offers a promising tool for rapid free energy scan of small functional group substitutions. © 2016 Wiley Periodicals, Inc. Biopolymers 105: 324-336, 2016.

Binding Site Protonation and Uphill Occlusion Control the Na+/K+-Pump Selectivity

The Na+/K+-pump resides in the plasma membrane and maintains the physiological K+ and Na+ concentration gradient across the cell membrane. It operates via a ping-pong mechanism, making iterative transitions between inward-facing (E1) and outward-facing (E2) conformations. The E1 conformation binds three Na+ from the cytosol and exports them using the energy from ATP hydrolysis. The release of Na+ and the binding of K+ at the extracellular side trigger the structural transition to the occluded E2 state, which imports two K+, as the pump returns to the E1 conformation. Although the general features of the pumping cycle are known, the detailed mechanism governing the selectivity in different states of the pump is not well understood. Using free energy molecular dynamics simulations based on the crystal structures of the pump in its E1 and E2 states, the current study reveals a switch in the binding site protonation state upon the E1-E2 transitions, giving the ion binding site its specificity when the pump is in different states. Ion selectivity of the states along the pump cycle has also been calculated. The relative binding free energy of Na+ and K+ in the binding sites calculated along the conformational transition of the protein reveals a self-correcting mechanism that ensures the selectivity and the stoichiometry of the pumping cycle.

Interpreting Thermodynamic Profiles of Aminoadamantane Compounds Inhibiting the M2 Proton Channel of Influenza A by Free Energy Calculations.

The development of novel anti-influenza drugs is of great importance because of the capability of influenza viruses to occasionally cross interspecies barriers and to rapidly mutate. One class of anti-influenza agents, aminoadamantanes, including the drugs amantadine and rimantadine now widely abandoned due to virus resistance, bind to and block the pore of the transmembrane domain of the M2 proton channel (M2TM) of influenza A. Here, we present one of the still rare studies that interprets thermodynamic profiles from isothermal titration calorimetry (ITC) experiments in terms of individual energy contributions to binding, calculated by the computationally inexpensive implicit solvent/implicit membrane molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) approach, for aminoadamantane compounds binding to M2TM of the avian "Weybridge" strain. For all eight pairs of aminoadamantane compounds considered, the trend of the predicted relative binding free energies and their individual components, effective binding energies and changes in the configurational entropy, agrees with experimental measures (ΔΔG, ΔΔH, TΔΔS) in 88, 88, and 50% of the cases. In addition, information yielded by the MM-PBSA approach about determinants of binding goes beyond that available in component quantities (ΔH, ΔS) from ITC measurements. We demonstrate how one can make use of such information to link thermodynamic profiles from ITC with structural causes on the ligand side and, ultimately, to guide decision making in lead optimization in a prospective manner, which results in an aminoadamantane derivative with improved binding affinity against M2TM(Weybridge).

Accurate Binding Free Energy Predictions in Fragment Optimization.

Predicting protein-ligand binding free energies is a central aim of computational structure-based drug design (SBDD)--improved accuracy in binding free energy predictions could significantly reduce costs and accelerate project timelines in lead discovery and optimization. The recent development and validation of advanced free energy calculation methods represents a major step toward this goal. Accurately predicting the relative binding free energy changes of modifications to ligands is especially valuable in the field of fragment-based drug design, since fragment screens tend to deliver initial hits of low binding affinity that require multiple rounds of synthesis to gain the requisite potency for a project. In this study, we show that a free energy perturbation protocol, FEP+, which was previously validated on drug-like lead compounds, is suitable for the calculation of relative binding strengths of fragment-sized compounds as well. We study several pharmaceutically relevant targets with a total of more than 90 fragments and find that the FEP+ methodology, which uses explicit solvent molecular dynamics and physics-based scoring with no parameters adjusted, can accurately predict relative fragment binding affinities. The calculations afford R(2)-values on average greater than 0.5 compared to experimental data and RMS errors of ca. 1.1 kcal/mol overall, demonstrating significant improvements over the docking and MM-GBSA methods tested in this work and indicating that FEP+ has the requisite predictive power to impact fragment-based affinity optimization projects.

A Python tool to set up relative free energy calculations in GROMACS.

Free energy calculations based on molecular dynamics (MD) simulations have seen a tremendous growth in the last decade. However, it is still difficult and tedious to set them up in an automated manner, as the majority of the present-day MD simulation packages lack that functionality. Relative free energy calculations are a particular challenge for several reasons, including the problem of finding a common substructure and mapping the transformation to be applied. Here we present a tool, alchemical-setup.py, that automatically generates all the input files needed to perform relative solvation and binding free energy calculations with the MD package GROMACS. When combined with Lead Optimization Mapper (LOMAP;Liu et al. in J Comput Aided Mol Des 27(9):755-770, 2013), recently developed in our group, alchemical-setup.py allows fully automated setup of relative free energy calculations in GROMACS. Taking a graph of the planned calculations and a mapping, both computed by LOMAP, our tool generates the topology and coordinate files needed to perform relative free energy calculations for a given set of molecules, and provides a set of simulation input parameters. The tool was validated by performing relative hydration free energy calculations for a handful of molecules from the SAMPL4 challenge (Mobley et al. in J Comput Aided Mol Des 28(4):135-150, 2014). Good agreement with previously published results and the straightforward way in which free energy calculations can be conducted make alchemical-setup.py a promising tool for automated setup of relative solvation and binding free energy calculations.

Relative Binding Free Energies of Adenine and Guanine to Damaged and Undamaged DNA in Human DNA Polymeraseη: Clues for Fidelity and Overall Efficiency.

Human DNA polymerase η (Pol η) plays an essential protective role against skin cancer caused by cyclobutane thymine–thymine dimers (TTDs), a frequent form of DNA damage arising from exposure to the sun. This enzyme rescues stalled replication forks at the TTDs by inserting bases opposite these DNA defects. Herein we calculate binding free energies for a free deoxyribose nucleotide triphosphate, dATP or dGTP, to Pol η complexed with undamaged or damaged DNA. The calculations indicate that the binding of dATP to the enzyme–DNA complex is thermodynamically favored for TTD-containing DNA over undamaged DNA, most likely because of more extensive hydrogen-bonding interactions between the TTD and the enzyme that hold the TTD more rigidly in place. The calculations also illustrate that dATP binding is thermodynamically favored over dGTP binding at both thymine positions of the TTD, most likely due to more persistent and stable hydrogen-bonding interactions between the TTD and dATP than between the TTD and dGTP. This free energy difference is slightly greater for binding at the 5′ thymine position than at the 3′ thymine position, presumably because of stabilization arising from the A:T base pair formed at the 3′ position of the TTD in the previous step of Pol η function. All of these trends in binding free energies are consistent with experimental measurements of binding strength, fidelity, processivity, and overall efficiency. The insights gained from this analysis have implications for drug design efforts aimed at modifying the binding properties of this enzyme for improving cancer chemotherapy treatments.

Influence of Solvent in Controlling Peptide-Surface Interactions.

Protein binding to surfaces is an important phenomenon in biology and in modern technological applications. Extensive experimental and theoretical research has been focused in recent years on revealing the factors that govern binding affinity to surfaces. Theoretical studies mainly focus on examining the contribution of the individual amino acids or, alternatively, the binding potential energies of the full peptide, which are unable to capture entropic contributions and neglect the dynamic nature of the system. We present here a methodology that involves the combination of nonequilibrium dynamics simulations with strategic mutation of polar residues to reveal the different factors governing the binding free energy of a peptide to a surface. Using a gold-binding peptide as an example, we show that relative binding free energies are a consequence of the balance between strong interactions of the peptide with the surface and the ability for the bulk solvent to stabilize the peptide.

Imidazotriazines: Spleen Tyrosine Kinase (Syk) Inhibitors Identified by Free-Energy Perturbation (FEP).

There has been significant interest in spleen tyrosine kinase (Syk) owing to its role in a number of disease states, including autoimmunity, inflammation, and cancer. Ongoing therapeutic programs have resulted in several compounds that are now in clinical use. Herein we report our optimization of the imidazopyrazine core scaffold of Syk inhibitors through the use of empirical and computational approaches. Free-energy perturbation (FEP) methods with MCPRO+ were undertaken to calculate the relative binding free energies for several alternate scaffolds. FEP was first applied retrospectively to determine if there is any predictive value; this resulted in 12 of 13 transformations being predicted in a directionally correct manner. FEP was then applied in a prospective manner to evaluate 17 potential targets, resulting in the realization of imidazotriazine 17 (3-(4-(3,4-dimethoxyphenylamino)imidazo[1,2-f][1,2,4]triazin-2-yl)benzamide), which shows a tenfold improvement in activity relative to the parent compound and no increase in atom count. Optimization of 17 led to compounds with nanomolar cellular activity.

Estimation of relative binding free energies of agonists for the α7 Nicotinic Acetylcholine Receptor: trying to improve clustering analysis in MM-GBSA

Estimation of relative binding free energies of agonists for the α7 Nicotinic Acetylcholine Receptor: trying to improve clustering analysis in MM-GBSA

Binding affinities by alchemical perturbation using QM/MM with a large QM system and polarizable MM model.

The most general way to improve the accuracy of binding-affinity calculations for protein-ligand systems is to use quantum-mechanical (QM) methods together with rigorous alchemical-perturbation (AP) methods. We explore this approach by calculating the relative binding free energy of two synthetic disaccharides binding to galectin-3 at a reasonably high QM level (dispersion-corrected density functional theory with a triple-zeta basis set) and with a sufficiently large QM system to include all short-range interactions with the ligand (744-748 atoms). The rest of the protein is treated as a collection of atomic multipoles (up to quadrupoles) and polarizabilities. Several methods for evaluating the binding free energy from the 3600 QM calculations are investigated in terms of stability and accuracy. In particular, methods using QM calculations only at the endpoints of the transformation are compared with the recently proposed non-Boltzmann Bennett acceptance ratio (NBB) method that uses QM calculations at several stages of the transformation. Unfortunately, none of the rigorous approaches give sufficient statistical precision. However, a novel approximate method, involving the direct use of QM energies in the Bennett acceptance ratio method, gives similar results as NBB but with better precision, ∼3 kJ/mol. The statistical error can be further reduced by performing a greater number of QM calculations.

Quantum Chemical Studies on Detail Mechanism of Nitrosylation of NAMI-A-HSA Adduct.

Hydrolysis of NAMI-A in NAMI-A-HSA (HSA = human serum albumin) and nitrosylation of hydrolyzed NAMI-A-HSA adduct have been studied in detail using density functional theory method. It has been observed that the chloride exchange reaction with water in the NAMI-A-HSA adduct follows an interchange dissociative mechanism passing through an unstable heptacoordinated activated complex. The computed free energy of activation (ΔG) and rate constant (k) for the hydrolysis process in aqueous medium are observed to be 24.85 kcal mol(-1) and 3.81 × 10(-6) s(-1), respectively. Nitrosylation of hydrolyzed NAMI-A-HSA adduct with nitric oxide is found to be thermodynamically more favorable with the incorporation of solvent effect and provides a detailed understanding related to the antimetastatic activity of the NAMI-A drug. This investigation shows that nitric oxide coordinates linearly to NAMI-A-HSA adduct leading to the reduction of ruthenium(III) to more active ruthenium(II), with the reduction potential of -2.32 V. Negative relative solvation and relative binding free energies suggest that the hydrolysis and nitrosylation reactions are found to be thermodynamically favorable and faster. Our computed results provide a detailed thermodynamics and kinetics which may be highly beneficial for understanding antimetastatic activity as well as the nitric oxide scavenging ability of NAMI-A.