Linear-scaling soft-core scheme for alchemical free energy calculations

Alchemical free energy calculations involving the removal or insertion of atoms into condensed phase systems generally make use of soft-core scaling of nonbonded interactions, designed to circumvent numerical instabilities that arise from weakly interacting "hard" atoms in close proximity. Current methods model soft-core atoms by introducing a nonlinear dependence between the shape of the interaction potential and the strength of the interaction. In this article, we propose a soft-core method that avoids introducing such a nonlinear dependence, through the application of a smooth flattening of the potential energy only in a region that is energetically accessible under normal conditions. We discuss the benefits that this entails and explore a selection of applications, including enhanced methods for the estimation of free energy differences and for the automated optimization of the placement of intermediate states in multistage alchemical calculations.     

Using free energy of binding calculations to improve the accuracy of virtual screening predictions

Virtual screening of small molecule databases against macromolecular targets was used to identify binding ligands and predict their lowest energy bound conformation (i.e., pose). AutoDock4-generated poses were rescored using mean-field pathway decoupling free energy of binding calculations and evaluated if these calculations improved virtual screening discrimination between bound and nonbound ligands. Two small molecule databases were used to evaluate the effectiveness of the rescoring algorithm in correctly identifying binders of L99A T4 lysozyme. Self-dock calculations of a database containing compounds with known binding free energies and cocrystal structures largely reproduced experimental measurements, although the mean difference between calculated and experimental binding free energies increased as the predicted bound poses diverged from the experimental poses. In addition, free energy rescoring was more accurate than AutoDock4 scores in discriminating between known binders and nonbinders, suggesting free energy rescoring could be a useful approach to reduce false positive predictions in virtual screening experiments.     

On the use of orientational restraints and symmetry corrections in alchemical free energy calculations

Alchemical free energy calculations are becoming a useful tool for calculating absolute binding free energies of small molecule ligands to proteins. Here, we find that the presence of multiple metastable ligand orientations can cause convergence problems when distance restraints alone are used. We demonstrate that the use of orientational restraints can greatly accelerate the convergence of these calculations. However, even with this acceleration, we find that sufficient sampling requires substantially longer simulations than are used in many published protocols. To further accelerate convergence, we introduce a new method of configuration space decomposition by orientation which reduces required simulation lengths by at least a factor of 5 in the cases examined. Our method is easily parallelizable, well suited for cases where a ligand cocrystal structure is not available, and can utilize initial orientations generated by docking packages.     

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.     

Optimal pairwise and non-pairwise alchemical pathways for free energy calculations of molecular transformation in solution phase

We estimate the global minimum variance path for computing the free energy insertion into or deletion of small molecules from a dense fluid. We perform this optimization over all pair potentials, irrespective of functional form, using functional optimization with a two-body approximation for the radial distribution function. Surprisingly, the optimal pairwise path obtained via this method is almost identical to the path obtained using a optimized generalized "soft core" potential reported by Pham and Shirts [J. Chem. Phys. 135, 034114 (2011)]. We also derive the lowest variance non-pairwise potential path for molecular insertion or deletion and compare its efficiency to the pairwise path. Under certain conditions, non-pairwise pathways can reduce the total variance by up to 60% compared to optimal pairwise pathways. However, optimal non-pairwise pathways do not appear generally feasible for practical free energy calculations because an accurate estimate of the free energy, the parameter that is itself is desired, is required for constructing this non-pairwise path. Additionally, simulations at most intermediate states of these non-pairwise paths have significantly longer correlation times, often exceeding standard simulation lengths for solvation of bulky molecules. The findings suggest that the previously obtained soft core pathway is the lowest variance pathway for molecular insertion or deletion in practice. The findings also demonstrate the utility of functional optimization for determining the efficiency of thermodynamic processes performed with molecular simulation.     

Quantitative calculations of antibody--antigen binding: steroid--DB3 binding energies by the linear interaction energy method

Linear interaction energy/molecular dynamics calculations have been used to compute steroid/antibody binding energies. The absolute binding affinities of 10 steroids to antibody DB3 and of a hapten to catalytic antibody 1E9 are computed and compared to experiment. A detailed analysis of the molecular origins of the observed binding patterns is provided. The binding energy of an untested steroid is predicted.     

Blinded predictions of host-guest standard free energies of binding in the SAMPL5 challenge

In the context of the SAMPL5 blinded challenge standard free energies of binding were predicted for a dataset of 22 small guest molecules and three different host molecules octa-acids (OAH and OAMe) and a cucurbituril (CBC). Three sets of predictions were submitted, each based on different variations of classical molecular dynamics alchemical free energy calculation protocols based on the double annihilation method. The first model (model A) yields a free energy of binding based on computed free energy changes in solvated and host-guest complex phases; the second (model B) adds long range dispersion corrections to the previous result; the third (model C) uses an additional standard state correction term to account for the use of distance restraints during the molecular dynamics simulations. Model C performs the best in terms of mean unsigned error for all guests (MUE \(3.2\,<\,3.4\,<\,3.6\,\text{kcal}\,\text{mol}^{-1}\)—95 % confidence interval) for the whole data set and in particular for the octa-acid systems (MUE \(1.7\,<\,1.9\,<\,2.1\,\text{kcal}\,\text{mol}^{-1}\)). The overall correlation with experimental data for all models is encouraging (\(R^2\, 0.65\,<\,0.70<0.75\)). The correlation between experimental and computational free energy of binding ranks as one of the highest with respect to other entries in the challenge. Nonetheless the large MUE for the best performing model highlights systematic errors, and submissions from other groups fared better with respect to this metric.     

Density functional calculations of surface free energies

We propose a general method of thermodynamic integration to find the free energy of a surface, where our integration parameter is taken to be the strain on the unit cell of the system (which in the example presented in this paper is simply the extension of the unit cell along the normal to the surface), and the integration is performed over the thermal average stress from a molecular dynamics run. In order to open up a vacuum gap in a continuous and reversible manner, an additional control interaction has been introduced. We also use temperature integration to find a linear relation for the temperature dependence of the free surface energy. These methods have been applied to the titanium dioxide (110) surface, using first principles density functional theory. A proof of principle calculation for zero temperature shows excellent agreement between the integral calculation and the difference in energy calculated by the DFT program. Calculations that have been performed at 295 and 1000 K give excellent agreement between the two integration methods.     

Using the fast fourier transform in binding free energy calculations

According to implicit ligand theory, the standard binding free energy is an exponential average of the binding potential of mean force (BPMF), an exponential average of the interaction energy between the unbound ligand ensemble and a rigid receptor. Here, we use the fast Fourier transform (FFT) to efficiently evaluate BPMFs by calculating interaction energies when rigid ligand configurations from the unbound ensemble are discretely translated across rigid receptor conformations. Results for standard binding free energies between T4 lysozyme and 141 small organic molecules are in good agreement with previous alchemical calculations based on (1) a flexible complex ( for 24 systems) and (2) flexible ligand with multiple rigid receptor configurations ( for 141 systems). While the FFT is routinely used for molecular docking, to our knowledge this is the first time that the algorithm has been used for rigorous binding free energy calculations. © 2017 Wiley Periodicals, Inc.     

Blinded predictions of standard binding free energies: lessons learned from the SAMPL6 challenge

In the context of the SAMPL6 challenges, series of blinded predictions of standard binding free energies were made with the SOMD software for a dataset of 27 host–guest systems featuring two octa-acids hosts (OA and TEMOA) and a cucurbituril ring (CB8) host. Three different models were used, ModelA computes the free energy of binding based on a double annihilation technique; ModelB additionally takes into account long-range dispersion and standard state corrections; ModelC additionally introduces an empirical correction term derived from a regression analysis of SAMPL5 predictions previously made with SOMD. The performance of each model was evaluated with two different setups; buffer explicitly matches the ionic strength from the binding assays, whereas no-buffer merely neutralizes the host–guest net charge with counter-ions. ModelC/no-buffer shows the lowest mean-unsigned error for the overall dataset (MUE 1.29?<?1.39?<?1.50 kcal mol^?1, 95% CI), while explicit modelling of the buffer improves significantly results for the CB8 host only. Correlation with experimental data ranges from excellent for the host TEMOA (R² 0.91?<?0.94?<?0.96), to poor for CB8 (R² 0.04?<?0.12?<?0.23). Further investigations indicate a pronounced dependence of the binding free energies on the modelled ionic strength, and variable reproducibility of the binding free energies between different simulation packages.     

Tinker‐OpenMM: Absolute and relative alchemical free energies using AMOEBA on GPUs

The capabilities of the polarizable force fields for alchemical free energy calculations have been limited by the high computational cost and complexity of the underlying potential energy functions. In this work, we present a GPU‐based general alchemical free energy simulation platform for polarizable potential AMOEBA. Tinker‐OpenMM, the OpenMM implementation of the AMOEBA simulation engine has been modified to enable both absolute and relative alchemical simulations on GPUs, which leads to a ∼200‐fold improvement in simulation speed over a single CPU core. We show that free energy values calculated using this platform agree with the results of Tinker simulations for the hydration of organic compounds and binding of host–guest systems within the statistical errors. In addition to absolute binding, we designed a relative alchemical approach for computing relative binding affinities of ligands to the same host, where a special path was applied to avoid numerical instability due to polarization between the different ligands that bind to the same site. This scheme is general and does not require ligands to have similar scaffolds. We show that relative hydration and binding free energy calculated using this approach match those computed from the absolute free energy approach. © 2017 Wiley Periodicals, Inc.     

FEW: A workflow tool for free energy calculations of ligand binding

In the later stages of drug design projects, accurately predicting relative binding affinities of chemically similar compounds to a biomolecular target is of utmost importance for making decisions based on the ranking of such compounds. So far, the extensive application of binding free energy approaches has been hampered by the complex and time‐consuming setup of such calculations. We introduce the free energy workflow (FEW) tool that facilitates setup and execution of binding free energy calculations with the AMBER suite for multiple ligands. FEW allows performing free energy calculations according to the implicit solvent molecular mechanics (MM‐PB(GB)SA), the linear interaction energy, and the thermodynamic integration approaches. We describe the tool's architecture and functionality and demonstrate in a show case study on Factor Xa inhibitors that the time needed for the preparation and analysis of free energy calculations is considerably reduced with FEW compared to a fully manual procedure. © 2013 Wiley Periodicals, Inc.     

Efficiency of alchemical free energy simulations. II. Improvements for thermodynamic integration

We attempt to optimize the efficiency of thermodynamic integration, as defined by the minimal number of unphysical intermediate states required for the computation of accurate and precise free energy differences. The suitability of various numerical quadrature methods is tested. In particular, we compare the trapezoidal rule, Simpson's rule, Gauss-Legendre, Gauss-Kronrod-Patterson, and Clenshaw-Curtis integration, as well as integration based on a cubic spline approximation of the integrand. We find that Simpson's rule and spline integration are already significantly more efficient that the trapezoidal rule, i.e., correct free energy differences can be obtained using fewer λ-states. We demonstrate that Simpson's rule can be used advantageously with nonequidistant values of the abscissa, which increases the flexibility of the method. Efficiency is enhanced even further if higher order methods, such as Gauss-Legendre, Gauss-Kronrod-Patterson, or Clenshaw-Curtis integration, are used; no more than seven λ-states, which in the case of Clenshaw-Curtis integration include the physical end states, were required for accurate results in all test problems studied. Thus, the performance of thermodynamic integration can equal that of Bennett's acceptance ratio method. We also show, however, that the high efficiency found here relies on the particular functional form of the soft-core potential used; overall, thermodynamic integration is more susceptible to the details of the hybrid Hamiltonian used than Bennett's acceptance ratio method. Therefore, we recommend Bennett's acceptance ratio method as the most robust method to compute alchemical free energy differences; nevertheless, scenarios when thermodynamic integration may be preferable are discussed.     

Computational prediction of monosaccharide binding free energies to lectins with linear interaction energy models

The linear interaction energy (LIE) method to compute binding free energies is applied to lectin‐monosaccharide complexes. Here, we calculate the binding free energies of monosaccharides to the Ralstonia solanacearum lectin (RSL) and the Pseudomonas aeruginosa lectin‐II (PA‐IIL). The standard LIE model performs very well for RSL, whereas the PA‐IIL system, where ligand binding involves two calcium ions, presents a major challenge. To overcome this, we explore a new variant of the LIE model, where ligand–metal ion interactions are scaled separately. This model also predicts the saccharide binding preference of PA‐IIL on mutation of the receptor, which may be useful for protein engineering of lectins. © 2012 Wiley Periodicals, Inc.     

Fast and accurate predictions of binding free energies using MM‐PBSA and MM‐GBSA

In the drug discovery process, accurate methods of computing the affinity of small molecules with a biological target are strongly needed. This is particularly true for molecular docking and virtual screening methods, which use approximated scoring functions and struggle in estimating binding energies in correlation with experimental values. Among the various methods, MM‐PBSA and MM‐GBSA are emerging as useful and effective approaches. Although these methods are typically applied to large collections of equilibrated structures of protein‐ligand complexes sampled during molecular dynamics in water, the possibility to reliably estimate ligand affinity using a single energy‐minimized structure and implicit solvation models has not been explored in sufficient detail. Herein, we thoroughly investigate this hypothesis by comparing different methods for the generation of protein‐ligand complexes and diverse methods for free energy prediction for their ability to correlate with experimental values. The methods were tested on a series of structurally diverse inhibitors of Plasmodium falciparum DHFR with known binding mode and measured affinities. The results showed that correlations between MM‐PBSA or MM‐GBSA binding free energies with experimental affinities were in most cases excellent. Importantly, we found that correlations obtained with the use of a single protein‐ligand minimized structure and with implicit solvation models were similar to those obtained after averaging over multiple MD snapshots with explicit water molecules, with consequent save of computing time without loss of accuracy. When applied to a virtual screening experiment, such an approach proved to discriminate between true binders and decoy molecules and yielded significantly better enrichment curves. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Hydration free energy calculations by the acceptance ratio method

The hydration energy difference between the alanine and glycine zwitter ions was calculated by both the free energy perturbation method and the acceptance ratio method. The calculations were carried out by using different increments of the mutation parameter λ, δλ = ? 0.05, ?0.10, and ?0.20. The free energy difference calculated by the acceptance ratio method was found to be approximately the same as an average of the two free energy differences in the forward and the backward directions calculated by the perturbation method. The results by the perturbation method were significantly affected by large δλ as compared with that by the acceptance ratio method. The statistical error caused by decreasing the simulation time for sampling equilibrium configurations is discussed.     

Virtual screening of integrase inhibitors by large scale binding free energy calculations: the SAMPL4 challenge

As part of the SAMPL4 blind challenge, filtered AutoDock Vina ligand docking predictions and large scale binding energy distribution analysis method binding free energy calculations have been applied to the virtual screening of a focused library of candidate binders to the LEDGF site of the HIV integrase protein. The computational protocol leveraged docking and high level atomistic models to improve enrichment. The enrichment factor of our blind predictions ranked best among all of the computational submissions, and second best overall. This work represents to our knowledge the first example of the application of an all-atom physics-based binding free energy model to large scale virtual screening. A total of 285 parallel Hamiltonian replica exchange molecular dynamics absolute protein-ligand binding free energy simulations were conducted starting from docked poses. The setup of the simulations was fully automated, calculations were distributed on multiple computing resources and were completed in a 6-weeks period. The accuracy of the docked poses and the inclusion of intramolecular strain and entropic losses in the binding free energy estimates were the major factors behind the success of the method. Lack of sufficient time and computing resources to investigate additional protonation states of the ligands was a major cause of mispredictions. The experiment demonstrated the applicability of binding free energy modeling to improve hit rates in challenging virtual screening of focused ligand libraries during lead optimization.     

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.     

Using physics-based pose predictions and free energy perturbation calculations to predict binding poses and relative binding affinities for FXR ligands in the D3R Grand Challenge 2

Computer-aided drug design has become an integral part of drug discovery and development in the pharmaceutical and biotechnology industry, and is nowadays extensively used in the lead identification and lead optimization phases. The drug design data resource (D3R) organizes challenges against blinded experimental data to prospectively test computational methodologies as an opportunity for improved methods and algorithms to emerge. We participated in Grand Challenge 2 to predict the crystallographic poses of 36 Farnesoid X Receptor (FXR)-bound ligands and the relative binding affinities for two designated subsets of 18 and 15 FXR-bound ligands. Here, we present our methodology for pose and affinity predictions and its evaluation after the release of the experimental data. For predicting the crystallographic poses, we used docking and physics-based pose prediction methods guided by the binding poses of native ligands. For FXR ligands with known chemotypes in the PDB, we accurately predicted their binding modes, while for those with unknown chemotypes the predictions were more challenging. Our group ranked #1st (based on the median RMSD) out of 46 groups, which submitted complete entries for the binding pose prediction challenge. For the relative binding affinity prediction challenge, we performed free energy perturbation (FEP) calculations coupled with molecular dynamics (MD) simulations. FEP/MD calculations displayed a high success rate in identifying compounds with better or worse binding affinity than the reference (parent) compound. Our studies suggest that when ligands with chemical precedent are available in the literature, binding pose predictions using docking and physics-based methods are reliable; however, predictions are challenging for ligands with completely unknown chemotypes. We also show that FEP/MD calculations hold predictive value and can nowadays be used in a high throughput mode in a lead optimization project provided that crystal structures of sufficiently high quality are available.