Binding affinity prediction with different force fields: examination of the linear interaction energy method

A systematic study of the linear interaction energy (LIE) method and the possible dependence of its parameterization on the force field and system (receptor binding site) is reported. We have calculated the binding free energy for nine different ligands in complex with P450cam using three different force fields (Amber95, Gromos87, and OPLS-AA). The results from these LIE calculations using our earlier parameterization give relative free energies of binding that agree remarkably well with the experimental data. However, the absolute energies are too positive for all three force fields, and it is clear that an additional constant term (gamma) is required in this case. Out of five examined LIE models, the same one emerges as the best for all three force fields, and this, in fact, corresponds to our earlier one apart from the addition of the constant gamma, which is almost identical for the three force fields. Thus, the present free energy calculations clearly indicate that the coefficients of the LIE method are independent of the force field used. Their relation to solvation free energies is also demonstrated. The only free parameter of the best model is gamma, which is found to depend on the hydrophobicity of the binding site. We also attempt to quantify the binding site hydrophobicity of four different proteins which shows that the ordering of gamma's for these sites reflects the fraction of hydrophobic surface area.     

Ligand binding affinity prediction by linear interaction energy methods

A recent method for estimating ligand binding affinities is extended. This method employs averages of interaction potential energy terms from molecular dynamics simulations or other thermal conformational sampling techniques. Incorporation of systematic deviations from electrostatic linear response, derived from free energy perturbation studies, into the absolute binding free energy expression significantly enhances the accuracy of the approach. This type of method may be useful for computational prediction of ligand binding strengths, e.g., in drug design applications.     

Calculation of ligand binding free energies from molecular dynamics simulations

A recently developed method for predicting binding affinities in ligand–receptor complexes, based on interaction energy averaging and conformational sampling by molecular dynamics simulation, is presented. Polar and nonpolar contributions to the binding free energy are approximated by a linear scaling of the corresponding terms in the average intermolecular interaction energy for the bound and free states of the ligand. While the method originally assumed the validity of electrostatic linear response, we show that incorporation of systematic deviations from linear response derived from free energy perturbation calculations enhances the accuracy of the approach. The method is applied to complexes of wild-type and mutant human dihydrofolate reductases with 2,4-diaminopteridine and 2,4-diaminoquinazoline inhibitors. It is shown that a binding energy accuracy of about 1 kcal/mol is attainable even for multiply ionized compounds, such as methotrexate, for which electrostatic interactions energies are very large. © 1998 John Wiley & Sons, Inc. Int J Quant Chem 69: 77–88, 1998     

The performance of ANI-ML potentials for ligand-n(H2O) interaction energies and estimation of hydration free energies from end-point MD simulations

Here, we investigate the performance of “Accurate NeurAl networK engINe for Molecular Energies” (ANI), trained on small organic compounds, on bulk systems including non-covalent interactions and applicability to estimate solvation (hydration) free energies using the interaction between the ligand and explicit solvent (water) from single-step MD simulations. The method is adopted from ANI using the Atomic Simulation Environment (ASE) and predicts the non-covalent interaction energies at the accuracy of wb97x/6-31G(d) level by a simple linear scaling for the conformations sampled by molecular dynamics (MD) simulations of ligand-n(H₂O) systems. For the first time, we test ANI potentials' abilities to reproduce solvation free energies using linear interaction energy (LIE) formulism by modifying the original LIE equation. Our results on ~250 different complexes show that the method can be accurate and have a correlation of R² = 0.88–0.89 (MAE <1.0 kcal/mol) to the experimental solvation free energies, outperforming current end-state methods. Moreover, it is competitive to other conventional free energy methods such as FEP and BAR with 15-20 × fold reduced computational cost.     

Predicting relative binding affinities of non-peptide HIV protease inhibitors with free energy perturbation calculations

The relative binding free energies in HIV protease of haloperidol thioketal (THK) and three of its derivatives were examined with free energy calculations. THK is a weak inhibitor (IC50 = 15 M) for which two cocrystal structures with HIV type 1 proteases have been solved [Rutenber, E. et al., J. Biol. Chem., 268 (1993) 15343]. A THK derivative with a phenyl group on C2 of the piperidine ring was expected to be a poor inhibitor based on experiments with haloperidol ketal and its 2- phenyl derivative (Caldera, P., personal communication). Our calculations predict that a 5-phenyl THK derivative, suggested based on examination of the crystal structure, will bind significantly better than THK. Although there are large error bars as estimated from hysteresis, the calculations predict that the 5-phenyl substituent is clearly favored over the 2-phenyl derivative as well as the parent compound. The unfavorable free energies of solvation of both phenyl THK derivatives relative to the parent compound contributed to their predicted binding free energies. In a third simulation, the change in binding free energy for 5-benzyl THK relative to THK was calculated. Although this derivative has a lower free energy in the protein, its decreased free energy of solvation increases the predicted G(bind) to the same range as that of the 2-phenyl derivative.     

Efficient free energy calculations on small molecule host-guest systems - a combined linear interaction energy/one-step perturbation approach

Two efficient methods to calculate binding affinities of ligands with proteins have been critically evaluated by using sixteen small ligand host-guest complexes. It is shown that both the one-step (OS) perturbation method and the linear interaction energy (LIE) method have complementing strengths and weaknesses and can be optimally combined in a new manner. The OS method has a sound theoretical basis to address the free energy of cavity formation, whereas the LIE approach is more versatile and efficient to calculate the free energy of adding charges to such cavities. The off-term, which is neglected in the original LIE equation, can be calculated without additional costs from the OS, offering a powerful synergy between the two methods. The LIE/OS approach presented here combines the best of two worlds and for the model systems studied here, is more accurate than and as efficient as the original methods. It has a sound theoretical background and no longer requires any empirical parameters. The method appears very well suited for application in lead-optimization programmes in drug research, where the structure and dynamics of a series of molecules is of interest, together with an accurate calculation of the binding free energy.     

Assessment of free energy predictors for ligand binding to a methyllysine histone code reader

Methyllysine histone code readers constitute a new promising group of potential drug targets. For instance, L3MBTL1, a malignant brain tumor (MBT) protein, selectively binds mono- and di-methyllysine epigenetic marks (KMe, KMe(2) ) that eventually results in the negative regulation of multiple genes through the E2F/Rb oncogenic pathway. There is a pressing need in potent and selective small-molecule probes that would enable further target validation and might become therapeutic leads. Such an endeavor would require efficient tools to assess the free energy of protein-ligand binding. However, due to an unparalleled function of the MBT binding pocket (i.e., selective binding to KMe/KMe(2) ) and because of its distinctive structure representing a small aromatic "cage," an accurate assessment of its binding affinity to a ligand appears to be a challenging task. Here, we report a comparative analysis of computationally affordable affinity predictors applied to a set of seven small-molecule ligands interacting with L3MBTL1. The analysis deals with novel ligands and targets, but applies widespread computational approaches and intuitive comparison metrics that makes this study compatible with and incremental to earlier large scale accounts on the efficiency of affinity predictors. Ultimately, this study has revealed three top performers, far ahead of the other techniques, including two scoring functions, PMF04 and PLP, along with a simulation-based method MM-PB/SA. We discuss why some methods may perform better than others on this target class, the limits of their application, as well as how the efficiency of the most CPU-demanding techniques could be optimized.     

Ab initio quantum mechanics-based free energy perturbation method for calculating relative solvation free energies

A free energy perturbation (FEP) method was developed that uses ab initio quantum mechanics (QM) for treating the solute molecules and molecular mechanics (MM) for treating the surroundings. Like our earlier results using AM1 semi empirical QMs, the ab initio QM/MM-based FEP method was shown to accurately calculate relative solvation free energies for a diverse set of small molecules that differ significantly in structure, aromaticity, hydrogen bonding potential, and electron density. Accuracy was similar to or better than conventional FEP methods. The QM/MM-based methods eliminate the need for time-consuming development of MM force field parameters, which are frequently required for drug-like molecules containing structural motifs not adequately described by MM. Future automation of the method and parallelization of the code for Linux 128/256/512 clusters is expected to enhance the speed and increase its use for drug design and lead optimization.     

Correcting electrostatic artifacts due to net-charge changes in the calculation of ligand binding free energies

Alchemically derived free energies are artifacted when the perturbed moiety has a nonzero net charge. The source of the artifacts lies in the effective treatment of the electrostatic interactions within and between the perturbed atoms and remaining (partial) charges in the simulated system. To treat the electrostatic interactions effectively, lattice-summation (LS) methods or cutoff schemes in combination with a reaction-field contribution are usually employed. Both methods render the charging component of the calculated free energies sensitive to essential parameters of the system like the cutoff radius or the box side lengths. Here, we discuss the results of three previously published studies of ligand binding. These studies presented estimates of binding free energies that were artifacted due to the charged nature of the ligands. We show that the size of the artifacts can be efficiently calculated and raw simulation data can be corrected. We compare the corrected results with experimental estimates and nonartifacted estimates from path-sampling methods. Although the employed correction scheme involves computationally demanding continuum-electrostatics calculations, we show that the correction estimate can be deduced from a small sample of configurations rather than from the entire ensemble. This observation makes the calculations of correction terms feasible for complex biological systems. To show the general applicability of the proposed procedure, we also present results where the correction scheme was used to correct independent free energies obtained from simulations employing a cutoff scheme or LS electrostatics. In this work, we give practical guidelines on how to apply the appropriate corrections easily.     

Calculation of binding free energies of inhibitors to plasmepsin II

An understanding at the atomic level of the driving forces of inhibitor binding to the protein plasmepsin (PM) II would be of interest to the development of drugs against malaria. To this end, three state of the art computational techniques to compute relative free energies-thermodynamic integration (TI), Hamiltonian replica-exchange (H-RE) TI, and comparison of bound versus unbound ligand energy and entropy-were applied to a protein-ligand system of PM II and several exo-3-amino-7-azabicyclo[2.2.1]heptanes and the resulting relative free energies were compared with values derived from experimental IC(50) values. For this large and flexible protein-ligand system, the simulations could not properly sample the relevant parts of the conformational space of the bound ligand, resulting in failure to reproduce the experimental data. Yet, the use of Hamiltonian replica exchange in conjunction with thermodynamic integration resulted in enhanced convergence and computational efficiency compared to standard thermodynamic integration calculations. The more approximate method of calculating only energetic and entropic contributions of the ligand in its bound and unbound states from conventional molecular dynamics (MD) simulations reproduced the major trends in the experimental binding free energies, which could be rationalized in terms of energetic and entropic characteristics of the different structural and physico-chemical properties of the protein and ligands.     

New approach to free energy of solvation applying continuum models to molecular dynamics simulation

A new approach to the calculation of the free energy of solvation from trajectories obtained by molecular dynamics simulation is presented. The free energy of solvation is computed as the sum of three contributions originated at the cavitation of the solute by the solvent, the solute-solvent nonpolar (repulsion and dispersion) interactions, and the electrostatic solvation of the solute. The electrostatic term is calculated based on ideas developed for the broadly used continuum models, the cavitational contribution from the excluded volume by the Claverie-Pierotti model, and the Van der Waals term directly from the molecular dynamics simulation. The proposed model is tested for diluted aqueous solutions of simple molecules containing a variety of chemically important functions: methanol, methylamine, water, methanethiol, and dichloromethane. These solutions were treated by molecular dynamics simulations using SPC/E water and the OPLS force field for the organic molecules. Obtained free energies of solvation are in very good agreement with experimental data.     

Minimum MD simulation length required to achieve reliable results in free energy perturbation calculations: Case study of relative binding free energies of fructose‐1,6‐bisphosphatase inhibitors

In an attempt to establish the criteria for the length of simulation to achieve the desired convergence of free energy calculations, two studies were carried out on chosen complexes of FBPase‐AMP mimics. Calculations were performed for varied length of simulations and for different starting configurations using both conventional‐ and QM/MM‐FEP methods. The results demonstrate that for small perturbations, 1248 ps simulation time could be regarded a reasonable yardstick to achieve convergence of the results. As the simulation time is extended, the errors associated with free energy calculations also gradually tapers off. Moreover, when starting the simulation from different initial configurations of the systems, the results are not changed significantly, when performed for 1248 ps. This study carried on FBPase‐AMP mimics corroborates well with our previous successful demonstration of requirement of simulation time for solvation studies, both by conventional and ab initio FEP. The establishment of aforementioned criteria of simulation length serves a useful benchmark in drug design efforts using FEP methodologies, to draw a meaningful and unequivocal conclusion. © 2011 Wiley Periodicals, Inc. J Comput Chem, 2011     

Net charge changes in the calculation of relative ligand‐binding free energies via classical atomistic molecular dynamics simulation

The calculation of binding free energies of charged species to a target molecule is a frequently encountered problem in molecular dynamics studies of (bio‐)chemical thermodynamics. Many important endogenous receptor‐binding molecules, enzyme substrates, or drug molecules have a nonzero net charge. Absolute binding free energies, as well as binding free energies relative to another molecule with a different net charge will be affected by artifacts due to the used effective electrostatic interaction function and associated parameters (e.g., size of the computational box). In the present study, charging contributions to binding free energies of small oligoatomic ions to a series of model host cavities functionalized with different chemical groups are calculated with classical atomistic molecular dynamics simulation. Electrostatic interactions are treated using a lattice‐summation scheme or a cutoff‐truncation scheme with Barker–Watts reaction‐field correction, and the simulations are conducted in boxes of different edge lengths. It is illustrated that the charging free energies of the guest molecules in water and in the host strongly depend on the applied methodology and that neglect of correction terms for the artifacts introduced by the finite size of the simulated system and the use of an effective electrostatic interaction function considerably impairs the thermodynamic interpretation of guest‐host interactions. Application of correction terms for the various artifacts yields consistent results for the charging contribution to binding free energies and is thus a prerequisite for the valid interpretation or prediction of experimental data via molecular dynamics simulation. Analysis and correction of electrostatic artifacts according to the scheme proposed in the present study should therefore be considered an integral part of careful free‐energy calculation studies if changes in the net charge are involved. © 2013 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Efficient calculation of relative binding free energies by umbrella sampling perturbation

An important task of biomolecular simulation is the calculation of relative binding free energies upon chemical modification of partner molecules in a biomolecular complex. The potential of mean force (PMF) along a reaction coordinate for association or dissociation of the complex can be used to estimate binding affinities. A free energy perturbation approach, termed umbrella sampling (US) perturbation, has been designed that allows an efficient calculation of the change of the PMF upon modification of a binding partner based on the trajectories obtained for the wild type reference complex. The approach was tested on the interaction of modified water molecules in aqueous solution and applied to in silico alanine scanning of a peptide‐protein complex. For the water interaction test case, excellent agreement with an explicit PMF calculation for each modification was obtained as long as no long range electrostatic perturbations were considered. For the alanine scanning, the experimentally determined ranking and binding affinity changes upon alanine substitutions could be reproduced within 0.1–2.0 kcal/mol. In addition, good agreement with explicitly calculated PMFs was obtained mostly within the sampling uncertainty. The combined US and perturbation approach yields, under the condition of sufficiently small system modifications, rigorously derived changes in free energy and is applicable to any PMF calculation. © 2014 Wiley Periodicals, Inc.     

Absolute free energies of binding of peptide analogs to the HIV-1 protease from molecular dynamics simulations

The constants of binding of five peptide analogs to the active site of the HIV-1 aspartic-protease are calculated based on a novel sampling scheme that is efficient and does not introduce any approximations in addition to the energy function used to describe the system. The results agree with experiments. The squared correlation coefficient of the calculated vs. the measured values is 0.79. The sampling scheme consists of a series of molecular dynamics integrations with biases. The biases are selected based on an estimate of the probability density function of the system in a way to explore the conformational space and to reduce the statistical error in the calculated binding constants. The molecular dynamics integrations are done with a vacuum potential using a short cutoff scheme. To estimate the probability density of the simulated system, the results of the molecular dynamics integrations are combined using an extension of the weighted histogram analysis method (C. Bartels, Chem. Phys. Letters 331 (2000) 446-454). The probability density of the solvated ligand-protein system is obtained by applying a correction for the use of the short cutoffs in the simulations and by taking into account solvation with an electrostatic term and a hydrophobic term. The electrostatic part of the solvation is determined by finite difference Poisson-Boltzmann calculations; the hydrophobic part of the solvation is set proportional to the solvent accessible surface area. Setting the hydrophobic surface tension parameter equal to 8 mol(-1) K(-1) A(-2), absolute binding constants are in the muM to nM range. This is in agreement with experiments. The standard errors determined from eight repeated binding constant determinations are a factor of 14 to 411. A single determination of a binding constant is done with 499700 steps of molecular dynamics integration and 4500 finite difference Poisson-Boltzmann calculations. The simulations can be analyzed with respect to conformational changes of the active site of the HIV-1 protease or the ligands upon binding and provide information that complements experiments and can be used in the drug development process.     

Effect of sampling on BACE‐1 ligands binding free energy predictions via MM‐PBSA calculations

The BACE‐1 enzyme is a prime target to find a cure to Alzheimer's disease. In this article, we used the MM‐PBSA approach to compute the binding free energies of 46 reported ligands to this enzyme. After showing that the most probable protonation state of the catalytic dyad is mono‐protonated (on ASP32), we performed a thorough analysis of the parameters influencing the sampling of the conformational space (in total, more than 35 μs of simulations were performed). We show that ten simulations of 2 ns gives better results than one of 50 ns. We also investigated the influence of the protein force field, the water model, the periodic boundary conditions artifacts (box size), as well as the ionic strength. Amber03 with TIP3P, a minimal distance of 1.0 nm between the protein and the box edges and a ionic strength of I = 0.2 M provides the optimal correlation with experiments. Overall, when using these parameters, a Pearson correlation coefficient of R = 0.84 (R ² = 0.71) is obtained for the 46 ligands, spanning eight orders of magnitude of K _d (from 0.017 nm to 2000 μM, i.e., from −14.7 to −3.7 kcal/mol), with a ligand size from 22 to 136 atoms (from 138 to 937 g/mol). After a two‐parameter fit of the binding affinities for 12 of the ligands, an error of RMSD = 1.7 kcal/mol was obtained for the remaining ligands. © 2017 Wiley Periodicals, Inc.     

Overcoming dissipation in the calculation of standard binding free energies by ligand extraction

This article addresses calculations of the standard free energy of binding from molecular simulations in which a bound ligand is extracted from its binding site by steered molecular dynamics (MD) simulations or equilibrium umbrella sampling (US). Host–guest systems are used as test beds to examine the requirements for obtaining the reversible work of ligand extraction. We find that, for both steered MD and US, marked irreversibilities can occur when the guest molecule crosses an energy barrier and suddenly jumps to a new position, causing dissipation of energy stored in the stretched molecule(s). For flexible molecules, this occurs even when a stiff pulling spring is used, and it is difficult to suppress in calculations where the spring is attached to the molecules by single, fixed attachment points. We, therefore, introduce and test a method, fluctuation‐guided pulling, which adaptively adjusts the spring's attachment points based on the guest's atomic fluctuations relative to the host. This adaptive approach is found to substantially improve the reversibility of both steered MD and US calculations for the present systems. The results are then used to estimate standard binding free energies within a comprehensive framework, termed attach‐pull‐release, which recognizes that the standard free energy of binding must include not only the pulling work itself, but also the work of attaching and then releasing the spring, where the release work includes an accounting of the standard concentration to which the ligand is discharged. © 2013 Wiley Periodicals, Inc.     

Calculation of relative binding affinities of fructose 1,6-bisphosphatase mutants with adenosine monophosphate using free energy perturbation method

The free energy perturbation (FEP) methodology is the most accurate means of estimating relative binding affinities between inhibitors and protein variants. In this article, the importance of hydrophobic and hydrophilic residues to the binding of adenosine monophosphate (AMP) to the fructose 1,6-bisphosphatase (FBPase), a target enzyme for type-II diabetes, was examined by FEP method. Five mutations were made to the FBPase enzyme with AMP inhibitor bound: 113Tyr --> 113Phe, 31Thr --> 31Ala, 31Thr --> 31Ser, 177Met --> 177Ala, and 30Leu --> 30Phe. These mutations test the strength of hydrogen bonds and van der Waals interactions between the ligand and enzyme. The calculated relative free energies indicated that: 113Tyr and 31Thr play an important role, each via two hydrogen bonds affecting the binding affinity of inhibitor AMP to FBPase, and any changes in these hydrogen bonds due to mutations on the protein will have significant effect on the binding affinity of AMP to FBPase, consistent to experimental results. Also, the free energy calculations clearly show that the hydrophilic interactions are more important than the hydrophobic interactions of the binding pocket of FBPase.     

On the accurate estimation of free energies using the jarzynski equality

The Jarzynski equality is one of the most widely celebrated and scrutinized nonequilibrium work theorems, relating free energy to the external work performed in nonequilibrium transitions. In practice, the required ensemble average of the Boltzmann weights of infinite nonequilibrium transitions is estimated as a finite sample average, resulting in the so-called Jarzynski estimator, . Alternatively, the second-order approximation of the Jarzynski equality, though seldom invoked, is exact for Gaussian distributions and gives rise to the Fluctuation-Dissipation estimator . Here we derive the parametric maximum-likelihood estimator (MLE) of the free energy considering unidirectional work distributions belonging to Gaussian or Gamma families, and compare this estimator to . We further consider bidirectional work distributions belonging to the same families, and compare the corresponding bidirectional to the Bennett acceptance ratio () estimator. We show that, for Gaussian unidirectional work distributions, is in fact the parametric MLE of the free energy, and as such, the most efficient estimator for this statistical family. We observe that and perform better than and , for unidirectional and bidirectional distributions, respectively. These results illustrate that the characterization of the underlying work distribution permits an optimal use of the Jarzynski equality. © 2018 Wiley Periodicals, Inc.