Effects of Water Placement on Predictions of Binding Affinities for p38α MAP Kinase Inhibitors

WaterKit: Thermodynamic Profiling of Protein Hydration Sites

Water desolvation is one of the key components of the free energy of binding of small molecules to their receptors. Thus, understanding the energetic balance of solvation and desolvation resulting from individual water molecules can be crucial when estimating ligand binding, especially when evaluating different molecules and poses as done in High-Throughput Virtual Screening (HTVS). Over the most recent decades, several methods were developed to tackle this problem, ranging from fast approximate methods (usually empirical functions using either discrete atom-atom pairwise interactions or continuum solvent models) to more computationally expensive and accurate ones, mostly based on Molecular Dynamics (MD) simulations, such as Grid Inhomogeneous Solvation Theory (GIST) or Double Decoupling. On one hand, MD-based methods are prohibitive to use in HTVS to estimate the role of waters on the fly for each ligand. On the other hand, fast and approximate methods show an unsatisfactory level of accuracy, with low agreement with results obtained with the more expensive methods. Here we introduce WaterKit, a new grid-based sampling method with explicit water molecules to calculate thermodynamic properties using the GIST method. Our results show that the discrete placement of water molecules is successful in reproducing the position of crystallographic waters with very high accuracy, as well as providing thermodynamic estimates with accuracy comparable to more expensive MD simulations. Unlike these methods, WaterKit can be used to analyze specific regions on the protein surface, (such as the binding site of a receptor), without having to hydrate and simulate the whole receptor structure. The results show the feasibility of a general and fast method to compute thermodynamic properties of water molecules, making it well-suited to be integrated in high-throughput pipelines such as molecular docking.

Water: new aspect of hydrogen bonding in the solid state

All water-water contacts in the crystal structures from the Cambridge Structural Database with d _OO ≤ 4.0 Å have been found. These contacts were analysed on the basis of their geometries and interaction energies from CCSD(T)/CBS calculations. The results show 6729 attractive water-water contacts, of which 4717 are classical hydrogen bonds (d _OH ≤ 3.0 Å and α ≥ 120°) with most being stronger than -3.3 kcal mol^-1. Beyond the region of these hydrogen bonds, there is a large number of attractive interactions (2062). The majority are antiparallel dipolar interactions, where the O-H bonds of two water molecules lying in parallel planes are oriented antiparallel to each other. Developing geometric criteria for these antiparallel dipoles (β₁, β₂ ≥ 160°, 80 ≤ α ≤ 140° and T _HOHO > 40°) yielded 1282 attractive contacts. The interaction energies of these antiparallel oriented water molecules are up to -4.7 kcal mol^-1, while most of the contacts have interaction energies in the range -0.9 to -2.1 kcal mol^-1. This study suggests that the geometric criteria for defining attractive water-water interactions should be broader than the classical hydrogen-bonding criteria, a change that may reveal undiscovered and unappreciated interactions controlling molecular structure and chemistry.

Machine learning on ligand-residue interaction profiles to significantly improve binding affinity prediction

Structure-based virtual screenings (SBVSs) play an important role in drug discovery projects. However, it is still a challenge to accurately predict the binding affinity of an arbitrary molecule binds to a drug target and prioritize top ligands from an SBVS. In this study, we developed a novel method, using ligand-residue interaction profiles (IPs) to construct machine learning (ML)-based prediction models, to significantly improve the screening performance in SBVSs. Such a kind of the prediction model is called an IP scoring function (IP-SF). We systematically investigated how to improve the performance of IP-SFs from many perspectives, including the sampling methods before interaction energy calculation and different ML algorithms. Using six drug targets with each having hundreds of known ligands, we conducted a critical evaluation on the developed IP-SFs. The IP-SFs employing a gradient boosting decision tree (GBDT) algorithm in conjunction with the MIN + GB simulation protocol achieved the best overall performance. Its scoring power, ranking power and screening power significantly outperformed the Glide SF. First, compared with Glide, the average values of mean absolute error and root mean square error of GBDT/MIN + GB decreased about 38 and 36%, respectively. Second, the mean values of squared correlation coefficient and predictive index increased about 225 and 73%, respectively. Third, more encouragingly, the average value of the areas under the curve of receiver operating characteristic for six targets by GBDT, 0.87, is significantly better than that by Glide, which is only 0.71. Thus, we expected IP-SFs to have broad and promising applications in SBVSs.

Sensitivity of Binding Free Energy Calculations to Initial Protein Crystal Structure

Binding free energy calculations using alchemical free energy (AFE) methods are widely considered to be the most rigorous tool in the computational drug discovery arsenal. Despite this, the calculations suffer from accuracy, precision, and reproducibility issues. In this publication, we perform a high-throughput study of more than a thousand AFE calculations, utilizing over 220 μs of total sampling time, on three different protein systems to investigate the impact of the initial crystal structure on the resulting binding free energy values. We also consider the influence of equilibration time and discover that the initial crystal structure can have a significant effect on free energy values obtained at short timescales that can manifest itself as a free energy difference of more than 1 kcal/mol. At longer timescales, these differences are largely overtaken by important rare events, such as torsional ligand motions, typically resulting in a much higher uncertainty in the obtained values. This work emphasizes the importance of rare event sampling and long-timescale dynamics in free energy calculations even for routinely performed alchemical perturbations. We conclude that an optimal protocol should not only concentrate computational resources on achieving convergence in the alchemical coupling parameter (λ) space but also on longer simulations and multiple repeats.

Incorporating structural similarity into a scoring function to enhance the prediction of binding affinities

In this study, we developed a novel algorithm to improve the screening performance of an arbitrary docking scoring function by recalibrating the docking score of a query compound based on its structure similarity with a set of training compounds, while the extra computational cost is neglectable. Two popular docking methods, Glide and AutoDock Vina were adopted as the original scoring functions to be processed with our new algorithm and similar improvement performance was achieved. Predicted binding affinities were compared against experimental data from ChEMBL and DUD-E databases. 11 representative drug receptors from diverse drug target categories were applied to evaluate the hybrid scoring function. The effects of four different fingerprints (FP2, FP3, FP4, and MACCS) and the four different compound similarity effect (CSE) functions were explored. Encouragingly, the screening performance was significantly improved for all 11 drug targets especially when CSE = S⁴ (S is the Tanimoto structural similarity) and FP2 fingerprint were applied. The average predictive index (PI) values increased from 0.34 to 0.66 and 0.39 to 0.71 for the Glide and AutoDock vina scoring functions, respectively. To evaluate the performance of the calibration algorithm in drug lead identification, we also imposed an upper limit on the structural similarity to mimic the real scenario of screening diverse libraries for which query ligands are general-purpose screening compounds and they are not necessarily structurally similar to reference ligands. Encouragingly, we found our hybrid scoring function still outperformed the original docking scoring function. The hybrid scoring function was further evaluated using external datasets for two systems and we found the PI values increased from 0.24 to 0.46 and 0.14 to 0.42 for A2AR and CFX systems, respectively. In a conclusion, our calibration algorithm can significantly improve the virtual screening performance in both drug lead optimization and identification phases with neglectable computational cost.

Thermodynamics and Kinetics of Drug-Target Binding by Molecular Simulation

Computational studies play an increasingly important role in chemistry and biophysics, mainly thanks to improvements in hardware and algorithms. In drug discovery and development, computational studies can reduce the costs and risks of bringing a new medicine to market. Computational simulations are mainly used to optimize promising new compounds by estimating their binding affinity to proteins. This is challenging due to the complexity of the simulated system. To assess the present and future value of simulation for drug discovery, we review key applications of advanced methods for sampling complex free-energy landscapes at near nonergodicity conditions and for estimating the rate coefficients of very slow processes of pharmacological interest. We outline the statistical mechanics and computational background behind this research, including methods such as steered molecular dynamics and metadynamics. We review recent applications to pharmacology and drug discovery and discuss possible guidelines for the practitioner. Recent trends in machine learning are also briefly discussed. Thanks to the rapid development of methods for characterizing and quantifying rare events, simulation's role in drug discovery is likely to expand, making it a valuable complement to experimental and clinical approaches.

How Well Does the Extended Linear Interaction Energy Method Perform in Accurate Binding Free Energy Calculations?

With continually increased computer power, molecular mechanics force field-based approaches, such as the endpoint methods of molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) and molecular mechanics generalized Born surface area (MM-GBSA), have been routinely applied in both drug lead identification and optimization. However, the MM-PB/GBSA method is not as accurate as the pathway-based alchemical free energy methods, such as thermodynamic integration (TI) or free energy perturbation (FEP). Although the pathway-based methods are more rigorous in theory, they suffer from slow convergence and computational cost. Moreover, choosing adequate perturbation routes is also crucial for the pathway-based methods. Recently, we proposed a new method, coined extended linear interaction energy (ELIE) method, to overcome some disadvantages of the MM-PB/GBSA method to improve the accuracy of binding free energy calculation. In this work, we have systematically assessed this approach using in total 229 protein-ligand complexes for eight protein targets. Our results showed that ELIE performed much better than the molecular docking and MM-PBSA method in terms of root-mean-square error (RMSE), correlation coefficient (R), predictive index (PI), and Kendall's τ. The mean values of PI, R, and τ are 0.62, 0.58, and 0.44 for ELIE calculations. We also explored the impact of the length of simulation, ranging from 1 to 100 ns, on the performance of binding free energy calculation. In general, extending simulation length up to 25 ns could significantly improve the performance of ELIE, while longer molecular dynamics (MD) simulation does not always perform better than short MD simulation. Considering both the computational efficiency and achieved accuracy, ELIE is adequate in filling the gap between the efficient docking methods and computationally demanding alchemical free energy methods. Therefore, ELIE provides a practical solution for the routine ranking of compounds in lead optimization.

The critical role of QM/MM X-ray refinement and accurate tautomer/protomer determination in structure-based drug design

Conventional protein:ligand crystallographic refinement uses stereochemistry restraints coupled with a rudimentary energy functional to ensure the correct geometry of the model of the macromolecule—along with any bound ligand(s)—within the context of the experimental, X-ray density. These methods generally lack explicit terms for electrostatics, polarization, dispersion, hydrogen bonds, and other key interactions, and instead they use pre-determined parameters (e.g. bond lengths, angles, and torsions) to drive structural refinement. In order to address this deficiency and obtain a more complete and ultimately more accurate structure, we have developed an automated approach for macromolecular refinement based on a two layer, QM/MM (ONIOM) scheme as implemented within our DivCon Discovery Suite and "plugged in" to two mainstream crystallographic packages: PHENIX and BUSTER. This implementation is able to use one or more region layer(s), which is(are) characterized using linear-scaling, semi-empirical quantum mechanics, followed by a system layer which includes the balance of the model and which is described using a molecular mechanics functional. In this work, we applied our Phenix/DivCon refinement method—coupled with our XModeScore method for experimental tautomer/protomer state determination—to the characterization of structure sets relevant to structure-based drug design (SBDD). We then use these newly refined structures to show the impact of QM/MM X-ray refined structure on our understanding of function by exploring the influence of these improved structures on protein:ligand binding affinity prediction (and we likewise show how we use post-refinement scoring outliers to inform subsequent X-ray crystallographic efforts). Through this endeavor, we demonstrate a computational chemistry ↔ structural biology (X-ray crystallography) "feedback loop" which has utility in industrial and academic pharmaceutical research as well as other allied fields.

MovableType Software for Fast Free Energy-Based Virtual Screening: Protocol Development, Deployment, Validation, and Assessment

For decades, the complicated energy surfaces found in macromolecular protein:ligand structures, which require large amounts of computational time and resources for energy state sampling, have been an inherent obstacle to fast, routine free energy estimation in industrial drug discovery efforts. Beginning in 2013, the Merz research group addressed this cost with the introduction of a novel sampling methodology termed “Movable Type” (MT). Using numerical integration methods, the MT method reduces the computational expense for energy state sampling by independently calculating each atomic partition function from an initial molecular conformation in order to estimate the molecular free energy using ensembles of the atomic partition functions. In this work, we report a software package, the DivCon Discovery Suite with the MovableType module from QuantumBio Inc., that performs this MT free energy estimation protocol in a fast, fully encapsulated manner. We discuss the computational procedures and improvements to the original work, and we detail the corresponding settings for this software package. Finally, we introduce two validation benchmarks to evaluate the overall robustness of the method against a broad range of protein:ligand structural cases. With these publicly available benchmarks, we show that the method can use a variety of input types and parameters and exhibits comparable predictability whether the method is presented with “expensive” X-ray structures or “inexpensively docked” theoretical models. We also explore some next steps for the method. The MovableType software is available at http://www.quantumbioinc.com/

An assessment of water placement algorithms in quantum mechanics/molecular mechanics modeling: the case of rhodopsins' first spectral absorption band maxima

Quantum mechanics/molecular mechanics (QM/MM) models are a widely used tool to obtain detailed insight into the properties and functioning of proteins. The outcome of QM/MM studies heavily depends on the quality of the applied QM/MM model. Prediction and right placement of internal water molecules in protein cavities is one of the critical parts of any QM/MM model construction. Herein, we performed a systematic study of four protein hydration algorithms. We tested these algorithms for their ability to predict X-ray-resolved water molecules for a set of membrane photosensitive rhodopsin proteins, as well as the influence of the applied water placement algorithms on the QM/MM calculated absorption maxima (λ_max) of these proteins. We used 49 rhodopsins and their intermediates with available X-ray structures as the test set. We found that a proper choice of hydration algorithms and setups is needed to predict functionally important water molecules in the chromophore-binding cavity of rhodopsins, such as the water cluster in the N-H region of bacteriorhodopsin or two water molecules in the binding pocket of bovine visual rhodopsin. The QM/MM calculated λ_max of rhodopsins is also quite sensitive to the applied protein hydration protocols. The best methodology allows obtaining an 18.0 nm average value for the absolute deviation of the calculated λ_max from the experimental λ_max. Although the major effect of water molecules on λ_max originates from the water molecules located in the binding pocket, the water molecules outside the binding pocket also affect the calculated λ_max mainly by causing a reorganization of the protein structure. The results reported in this study can be used for the evaluation and further development of hydration methodologies, in general, and rhodopsin QM/MM models, in particular.

GemSpot: A Pipeline for Robust Modeling of Ligands into Cryo-EM Maps

Producing an accurate atomic model of biomolecule-ligand interactions from maps generated by cryoelectron microscopy (cryo-EM) often presents challenges inherent to the methodology and the dynamic nature of ligand binding. Here, we present GemSpot, an automated pipeline of computational chemistry methods that take into account EM map potentials, quantum mechanics energy calculations, and water molecule site prediction to generate candidate poses and provide a measure of the degree of confidence. The pipeline is validated through several published cryo-EM structures of complexes in different resolution ranges and various types of ligands. In all cases, at least one identified pose produced both excellent interactions with the target and agreement with the map. GemSpot will be valuable for the robust identification of ligand poses and drug discovery efforts through cryo-EM.

Advances in the Treatment of Explicit Water Molecules in Docking and Binding Free Energy Calculations

Background:

The inclusion of direct effects mediated by water during the ligandreceptor recognition is a hot-topic of modern computational chemistry applied to drug discovery and development. Docking or virtual screening with explicit hydration is still debatable, despite the successful cases that have been presented in the last years. Indeed, how to select the water molecules that will be included in the docking process or how the included waters should be treated remain open questions.

Objective:

In this review, we will discuss some of the most recent methods that can be used in computational drug discovery and drug development when the effect of a single water, or of a small network of interacting waters, needs to be explicitly considered.

Results:

Here, we analyse the software to aid the selection, or to predict the position, of water molecules that are going to be explicitly considered in later docking studies. We also present software and protocols able to efficiently treat flexible water molecules during docking, including examples of applications. Finally, we discuss methods based on molecular dynamics simulations that can be used to integrate docking studies or to reliably and efficiently compute binding energies of ligands in presence of interfacial or bridging water molecules.

Conclusions:

Software applications aiding the design of new drugs that exploit water molecules, either as displaceable residues or as bridges to the receptor, are constantly being developed. Although further validation is needed, workflows that explicitly consider water will probably become a standard for computational drug discovery soon.

Modelling flexible protein–ligand binding in p38α MAP kinase using the QUBE force field

The accuracy of quantum mechanical bespoke (QUBE) force fields for protein–ligand binding free energy calculations are benchmarked against experiment.

Dissecting the Energetics of Intrinsically Disordered Proteins via a Hybrid Experimental and Computational Approach

Intrinsically disordered proteins (IDPs) play important roles in biology, but little is known about the energetics of their inter-residue interactions. Methods that have been successfully applied to analyze the energetics of globular proteins are not applicable to the fluctuating partially ordered ensembles populated by IDPs. A combined computational experimental strategy is introduced for analyzing the energetic role of individual residues in the free state of IDPs. The approach combines experimental measurements of the binding of wild-type and mutant IDPs to their partners with alchemical free energy calculations of the structured complexes. These data allow quantitative information to be deduced about the free state via a thermodynamic cycle. The approach is validated by the analysis of the effects of mutations upon the binding free energy of the ovomucoid inhibitor third binding domain to its partners and is applied to the C-terminal domain of the measles virus nucleoprotein, a 125-residue IDP involved in the RNA transcription and replication of measles virus. The analysis reveals significant inter-residue interactions in the unbound IDP and suggests a biological role for them. The work demonstrates that advances in force fields and computational hardware have now led to the point where it is possible to develop methods, which integrate experimental and computational techniques to reveal insights that cannot be studied using either technique alone.

The role of hydration effects in 5-fluorouridine binding to SOD1: insight from a new 3D-RISM-KH based protocol for including structural water in docking simulations

Misfolded Cu/Zn superoxide dismutase enzyme (SOD1) shows prion-like propagation in neuronal cells leading to neurotoxic aggregates that are implicated in amyotrophic lateral sclerosis (ALS). Tryptophan-32 (W32) in SOD1 is part of a potential site for templated conversion of wild type SOD1. This W32 binding site is located on a convex, solvent exposed surface of the SOD1 suggesting that hydration effects can play an important role in ligand recognition and binding. A recent X-ray crystal structure has revealed that 5-Fluorouridine (5-FUrd) binds at the W32 binding site and can act as a pharmacophore scaffold for the development of anti-ALS drugs. In this study, a new protocol is developed to account for structural (non-displaceable) water molecules in docking simulations and successfully applied to predict the correct docked conformation binding modes of 5-FUrd at the W32 binding site. The docked configuration is within 0.58 Å (RMSD) of the observed configuration. The docking protocol involved calculating a hydration structure around SOD1 using molecular theory of solvation (3D-RISM-KH, 3D-Reference Interaction Site Model-Kovalenko-Hirata) whereby, non-displaceable water molecules are identified for docking simulations. This protocol was also used to analyze the hydrated structure of the W32 binding site and to explain the role of solvation in ligand recognition and binding to SOD1. Structural water molecules mediate hydrogen bonds between 5-FUrd and the receptor, and create an environment favoring optimal placement of 5-FUrd in the W32 binding site.

Using AMBER18 for Relative Free Energy Calculations

With renewed interest in free energy methods in contemporary structure-based drug design, there is a pressing need to validate against multiple targets and force fields to assess the overall ability of these methods to accurately predict relative binding free energies. We computed relative binding free energies using graphics processing unit accelerated thermodynamic integration (GPU-TI) on a data set originally assembled by Schrödinger, Inc. Using their GPU free energy code (FEP+) and the OPLS2.1 force field combined with the REST2 enhanced sampling approach, these authors obtained an overall MUE of 0.9 kcal/mol and an overall RMSD of 1.14 kcal/mol. In our study using GPU-TI from AMBER with the AMBER14SB/GAFF1.8 force field but without enhanced sampling, we obtained an overall MUE of 1.17 kcal/mol and an overall RMSD of 1.50 kcal/mol for the 330 perturbations contained in this data set. A more detailed analysis of our results suggested that the observed differences between the two studies arise from differences in sampling protocols along with differences in the force fields employed. Future work should address the problem of establishing benchmark quality results with robust statistical error bars obtained through multiple independent runs and enhanced sampling, which is possible with the GPU-accelerated features in AMBER.

Effect of set up protocols on the accuracy of alchemical free energy calculation over a set of ACK1 inhibitors

Hit-to-lead virtual screening frequently relies on a cascade of computational methods that starts with rapid calculations applied to a large number of compounds and ends with more expensive computations restricted to a subset of compounds that passed initial filters. This work focuses on set up protocols for alchemical free energy (AFE) scoring in the context of a Docking–MM/PBSA–AFE cascade. A dataset of 15 congeneric inhibitors of the ACK1 protein was used to evaluate the performance of AFE set up protocols that varied in the steps taken to prepare input files (using previously docked and best scored poses, manual selection of poses, manual placement of binding site water molecules). The main finding is that use of knowledge derived from X-ray structures to model binding modes, together with the manual placement of a bridging water molecule, improves the R² from 0.45 ± 0.06 to 0.76 ± 0.02 and decreases the mean unsigned error from 2.11 ± 0.08 to 1.24 ± 0.04 kcal mol^-1. By contrast a brute force automated protocol that increased the sampling time ten-fold lead to little improvements in accuracy. Besides, it is shown that for the present dataset hysteresis can be used to flag poses that need further attention even without prior knowledge of experimental binding affinities.

Hydrogen-Bond-Dependent Conformational Switching: A Computational Challenge from Experimental Thermochemistry

We have compiled an experimental data set (SWITCH10) of equilibrium constants for a series of hydrogen-bond-dependent conformational switches. These organic molecules possess common functionalities and are representative in terms of size and composition of systems routinely studied computationally. They exist as two well-defined conformations which serve as a useful tool to benchmark computational estimates of experimental Gibbs energy differences. We examine the performance of HF theory and a variety of density functionals (B3LYP, B3LYP-D3, CAM-B3LYP, ωB97X-D, M06-2X) against these experimental benchmarks. Surprisingly, despite a strong similarity between the two switch conformations, the average errors (0.4-1.7 kcal·mol^-1) obtained across the data set for all methods are larger than obtained with HF calculations. B3LYP was found to outperform implicitly and explicitly dispersion-corrected functionals, with an average error smaller by 1 kcal·mol^-1. Unsystematic errors in the optimized structures were found to contribute to the relatively poor performance obtained, while quasi-rigid rotor harmonic oscillator thermal contributions are important in improving the accuracy of computed Gibbs energy differences. These results emphasize the challenge of quantitative accuracy in computing solution-phase thermochemistry for flexible systems and caution against the often used (but unstated) assumption of favorable error cancellation in comparing conformers or stereoisomers.

Comparison of free-energy methods using a tripeptide-water model system

We investigate the ability of several free-energy calculation methods to combine two alchemical changes. We use Bennett acceptance ratio (BAR), thermodynamic integration (TI), extended TI (X-TI), and enveloping distribution sampling (EDS) to perturb a water molecule, which is restrained to an amino acid that is also being perturbed. In addition to these pairwise methods, we present two two-dimensional approaches, EDS-TI and two-dimensional TI (2D-TI). We compare feasibility, efficiency and usability of these methods in regard to our simple model system, which mimics the displacement of a water molecule in the active site of a protein on residue mutation. The correct treatment of structural water has been shown to greatly aid binding affinity calculations in some cases that remained elusive otherwise. This is of broad interest in, for example, drug design, and we conclude that thus far, the pairwise method BAR and also the newer X-TI remain the most suitable methods to treat this problem as long as few end states are involved. © 2018 Wiley Periodicals, Inc.

Calculate protein-ligand binding affinities with the extended linear interaction energy method: application on the Cathepsin S set in the D3R Grand Challenge 3

We participated in the Cathepsin S (CatS) sub-challenge of the Drug Design Data Resource (D3R) Grand Challenge 3 (GC3) in 2017 to blindly predict the binding poses of 24 CatS-bound ligands, the binding affinity ranking of 136 ligands, and the binding free energies of a subset of 33 ligands in Stage 1A and Stage 2. Our submitted predictions ranked relatively well compared to the submissions from other participants. Here we present our methodologies used in the challenge. For the binding pose prediction, we employed the Glide module in the Schrodinger Suite 2017 and AutoDock Vina. For the binding affinity/free energy prediction, we carried out molecular dynamics simulations of the complexes in explicit water solvent with counter ions, and then estimated the binding free energies with our newly developed model of extended linear interaction energy (ELIE), which is inspired by two other popular end-point approaches: the linear interaction energy (LIE) method, and the molecular mechanics with Poisson-Boltzmann surface area solvation method (MM/PBSA). Our studies suggest that ELIE is a good trade-off between efficiency and accuracy, and it is appropriate for filling the gap between the high-throughput docking and scoring methods and the rigorous but much more computationally demanding methods like free energy perturbation (FEP) or thermodynamics integration (TI) in computer-aided drug design (CADD) projects.

Detailed potential of mean force studies on host-guest systems from the SAMPL6 challenge

Accurately predicting receptor-ligand binding free energies is one of the holy grails of computational chemistry with many applications in chemistry and biology. Many successes have been reported, but issues relating to sampling and force field accuracy remain significant issues affecting our ability to reliably calculate binding free energies. In order to explore these issues in more detail we have examined a series of small host-guest complexes from the SAMPL6 blind challenge, namely octa-acids (OAs)-guest complexes and Curcurbit[8]uril (CB8)-guest complexes. Specifically, potential of mean force studies using umbrella sampling combined with the weighted histogram method were carried out on both systems with both known and unknown binding affinities. We find that using standard force fields and straightforward simulation protocols we are able to obtain satisfactory results, but that simply scaling our results allows us to significantly improve our predictive ability for the unknown test sets: the overall RMSD of the binding free energy versus experiment is reduced from 5.59 to 2.36 kcal/mol; for the CB8 test system, the RMSD goes from 8.04 to 3.51 kcal/mol, while for the OAs test system, the RSMD goes from 2.89 to 0.95 kcal/mol. The scaling approach was inspired by studies on structurally related known benchmark sets: by simply scaling, the RMSD was reduced from 6.23 to 1.19 kcal/mol and from 2.96 to 0.62 kcal/mol for the CB8 benchmark system and the OA benchmark system, respectively. We find this scaling procedure to correct absolute binding affinities to be highly effective especially when working across a "congeneric" series with similar charge states. It is less successful when applied to mixed ligands with varied charges and chemical characteristics, but improvement is still realized in the present case. This approach suggests that there are large systematic errors in absolute binding free energy calculations that can be straightforwardly accounted for using a scaling procedure. Random errors are still an issue, but near chemical accuracy can be obtained using the present strategy in select cases.

Predicting resistance of clinical Abl mutations to targeted kinase inhibitors using alchemical free-energy calculations

The therapeutic effect of targeted kinase inhibitors can be significantly reduced by intrinsic or acquired resistance mutations that modulate the affinity of the drug for the kinase. In cancer, the majority of missense mutations are rare, making it difficult to predict their impact on inhibitor affinity. This complicates the practice of precision medicine, pairing of patients with clinical trials, and development of next-generation inhibitors. Here, we examine the potential for alchemical free-energy calculations to predict how kinase mutations modulate inhibitor affinities to Abl, a major target in chronic myelogenous leukemia (CML). We find these calculations can achieve useful accuracy in predicting resistance for a set of eight FDA-approved kinase inhibitors across 144 clinically-identified point mutations, achieving a root mean square error in binding free energy changes of 1.1 0.9 1.3 kcal/mol (95% confidence interval) and correctly classifying mutations as resistant or susceptible with 88 82 93 % accuracy. Since these calculations are fast on modern GPUs, this benchmark establishes the potential for physical modeling to collaboratively support the rapid assessment and anticipation of the potential for patient mutations to affect drug potency in clinical applications.

Impact of domain knowledge on blinded predictions of binding energies by alchemical free energy calculations

The Drug Design Data Resource (D3R) consortium organises blinded challenges to address the latest advances in computational methods for ligand pose prediction, affinity ranking, and free energy calculations. Within the context of the second D3R Grand Challenge several blinded binding free energies predictions were made for two congeneric series of Farsenoid X Receptor (FXR) inhibitors with a semi-automated alchemical free energy calculation workflow featuring FESetup and SOMD software tools. Reasonable performance was observed in retrospective analyses of literature datasets. Nevertheless, blinded predictions on the full D3R datasets were poor due to difficulties encountered with the ranking of compounds that vary in their net-charge. Performance increased for predictions that were restricted to subsets of compounds carrying the same net-charge. Disclosure of X-ray crystallography derived binding modes maintained or improved the correlation with experiment in a subsequent rounds of predictions. The best performing protocols on D3R set1 and set2 were comparable or superior to predictions made on the basis of analysis of literature structure activity relationships (SAR)s only, and comparable or slightly inferior, to the best submissions from other groups.

Relative Binding Free Energy Calculations in Drug Discovery: Recent Advances and Practical Considerations

Accurate in silico prediction of protein-ligand binding affinities has been a primary objective of structure-based drug design for decades due to the putative value it would bring to the drug discovery process. However, computational methods have historically failed to deliver value in real-world drug discovery applications due to a variety of scientific, technical, and practical challenges. Recently, a family of approaches commonly referred to as relative binding free energy (RBFE) calculations, which rely on physics-based molecular simulations and statistical mechanics, have shown promise in reliably generating accurate predictions in the context of drug discovery projects. This advance arises from accumulating developments in the underlying scientific methods (decades of research on force fields and sampling algorithms) coupled with vast increases in computational resources (graphics processing units and cloud infrastructures). Mounting evidence from retrospective validation studies, blind challenge predictions, and prospective applications suggests that RBFE simulations can now predict the affinity differences for congeneric ligands with sufficient accuracy and throughput to deliver considerable value in hit-to-lead and lead optimization efforts. Here, we present an overview of current RBFE implementations, highlighting recent advances and remaining challenges, along with examples that emphasize practical considerations for obtaining reliable RBFE results. We focus specifically on relative binding free energies because the calculations are less computationally intensive than absolute binding free energy (ABFE) calculations and map directly onto the hit-to-lead and lead optimization processes, where the prediction of relative binding energies between a reference molecule and new ideas (virtual molecules) can be used to prioritize molecules for synthesis. We describe the critical aspects of running RBFE calculations, from both theoretical and applied perspectives, using a combination of retrospective literature examples and prospective studies from drug discovery projects. This work is intended to provide a contemporary overview of the scientific, technical, and practical issues associated with running relative binding free energy simulations, with a focus on real-world drug discovery applications. We offer guidelines for improving the accuracy of RBFE simulations, especially for challenging cases, and emphasize unresolved issues that could be improved by further research in the field.

Validating the Water Flooding Approach by Comparing It to Grand Canonical Monte Carlo Simulations

The study of the function of proteins on a quantitative level requires consideration of the water molecules in and around the protein. This requirement presents a major computational challenge due to the fact that the insertion of water molecules can have a very high activation barrier and would require a long simulation time. Recently, we developed a water flooding (WF) approach which is based on a postprocessing Monte Carlo ranking of possible water configurations. This approach appears to provide a very effective way for assessing the insertion free energies and determining the most likely configurations of the internal water molecules. Although the WF approach was used effectively in modeling challenging systems that have not been addressed reliably by other microscopic approaches, it was not validated by a comparison to the more rigorous grand canonical Monte Carlo (GCMC) method. Here we validate the WF approach by comparing its performance to that of the GCMC method. It is found that the WF approach reproduces the GCMC results in well-defined test cases but does so much faster. This established the WF approach as a useful strategy for finding correct water configurations in proteins and thus to provide a powerful way for studies of the functions of proteins.

Computationally-guided optimization of small-molecule inhibitors of the Aurora A kinase-TPX2 protein-protein interaction

Free energy perturbation theory, in combination with enhanced sampling of protein-ligand binding modes, is evaluated in the context of fragment-based drug design, and used to design two new small-molecule inhibitors of the Aurora A kinase-TPX2 protein-protein interaction.

Understanding ligand-receptor non-covalent binding kinetics using molecular modeling

Kinetic properties may serve as critical differentiators and predictors of drug efficacy and safety, in addition to the traditionally focused binding affinity. However the quantitative structure-kinetics relationship (QSKR) for modeling and ligand design is still poorly understood. This review provides an introduction to the kinetics of drug binding from a fundamental chemistry perspective. We focus on recent developments of computational tools and their applications to non-covalent binding kinetics.

Resolving the problem of trapped water in binding cavities: prediction of host-guest 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.

Hydration of proteins and nucleic acids: Advances in experiment and theory. A review

BACKGROUND

Most biological processes involve water, and the interactions of biomolecules with water affect their structure, function and dynamics.

SCOPE OF REVIEW

This review summarizes the current knowledge of protein and nucleic acid interactions with water, with a special focus on the biomolecular hydration layer. Recent developments in both experimental and computational methods that can be applied to the study of hydration structure and dynamics are reviewed, including software tools for the prediction and characterization of hydration layer properties.

MAJOR CONCLUSIONS

In the last decade, important advances have been made in our understanding of the factors that determine how biomolecules and their aqueous environment influence each other. Both experimental and computational methods contributed to the gradually emerging consensus picture of biomolecular hydration.

GENERAL SIGNIFICANCE

An improved knowledge of the structural and thermodynamic properties of the hydration layer will enable a detailed understanding of the various biological processes in which it is involved, with implications for a wide range of applications, including protein-structure prediction and structure-based drug design.

Prediction of Water Binding to Protein Hydration Sites with a Discrete, Semiexplicit Solvent Model

Buried water molecules are ubiquitous in protein structures and are found at the interface of most protein-ligand complexes. Determining their distribution and thermodynamic effect is a challenging yet important task, of great of practical value for the modeling of biomolecular structures and their interactions. In this study, we present a novel method aimed at the prediction of buried water molecules in protein structures and estimation of their binding free energies. It is based on a semiexplicit, discrete solvation model, which we previously introduced in the context of small molecule hydration. The method is applicable to all macromolecular structures described by a standard all-atom force field, and predicts complete solvent distribution within a single run with modest computational cost. We demonstrate that it indicates positions of buried hydration sites, including those filled by more than one water molecule, and accurately differentiates them from sterically accessible to water but void regions. The obtained estimates of water binding free energies are in fair agreement with reference results determined with the double decoupling method.

Evaluation of Selected Classical Force Fields for Alchemical Binding Free Energy Calculations of Protein-Carbohydrate Complexes

Protein-carbohydrate recognition is crucial in many vital biological processes including host-pathogen recognition, cell-signaling, and catalysis. Accordingly, computational prediction of protein-carbohydrate binding free energies is of enormous interest for drug design. However, the accuracy of current force fields (FFs) for predicting binding free energies of protein-carbohydrate complexes is not well understood owing to technical challenges such as the highly polar nature of the complexes, anomerization, and conformational flexibility of carbohydrates. The present study evaluated the performance of alchemical predictions of binding free energies with the GAFF1.7/AM1-BCC and GLYCAM06j force fields for modeling protein-carbohydrate complexes. Mean unsigned errors of 1.1 ± 0.06 (GLYCAM06j) and 2.6 ± 0.08 (GAFF1.7/AM1-BCC) kcal·mol(-1) are achieved for a large data set of monosaccharide ligands for Ralstonia solanacearum lectin (RSL). The level of accuracy provided by GLYCAM06j is sufficient to discriminate potent, moderate, and weak binders, a goal that has been difficult to achieve through other scoring approaches. Accordingly, the protocols presented here could find useful applications in carbohydrate-based drug and vaccine developments.

Spatial analysis and quantification of the thermodynamic driving forces in protein-ligand binding: binding site variability

The thermodynamic driving forces behind small molecule-protein binding are still not well-understood, including the variability of those forces associated with different types of ligands in different binding pockets. To better understand these phenomena we calculate spatially resolved thermodynamic contributions of the different molecular degrees of freedom for the binding of propane and methanol to multiple pockets on the proteins Factor Xa and p38 MAP kinase. Binding thermodynamics are computed using a statistical thermodynamics based end-point method applied on a canonical ensemble comprising the protein-ligand complexes and the corresponding free states in an explicit solvent environment. Energetic and entropic contributions of water and ligand degrees of freedom computed from the configurational ensemble provide an unprecedented level of detail into the mechanisms of binding. Direct protein-ligand interaction energies play a significant role in both nonpolar and polar binding, which is comparable to water reorganization energy. Loss of interactions with water upon binding strongly compensates these contributions leading to relatively small binding enthalpies. For both solutes, the entropy of water reorganization is found to favor binding in agreement with the classical view of the "hydrophobic effect". Depending on the specifics of the binding pocket, both energy-entropy compensation and reinforcement mechanisms are observed. It is notable to have the ability to visualize the spatial distribution of the thermodynamic contributions to binding at atomic resolution showing significant differences in the thermodynamic contributions of water to the binding of propane versus methanol.

A large-scale test of free-energy simulation estimates of protein-ligand binding affinities

We have performed a large-scale test of alchemical perturbation calculations with the Bennett acceptance-ratio (BAR) approach to estimate relative affinities for the binding of 107 ligands to 10 different proteins. Employing 20-Å truncated spherical systems and only one intermediate state in the perturbations, we obtain an error of less than 4 kJ/mol for 54% of the studied relative affinities and a precision of 0.5 kJ/mol on average. However, only four of the proteins gave acceptable errors, correlations, and rankings. The results could be improved by using nine intermediate states in the simulations or including the entire protein in the simulations using periodic boundary conditions. However, 27 of the calculated affinities still gave errors of more than 4 kJ/mol, and for three of the proteins the results were not satisfactory. This shows that the performance of BAR calculations depends on the target protein and that several transformations gave poor results owing to limitations in the molecular-mechanics force field or the restricted sampling possible within a reasonable simulation time. Still, the BAR results are better than docking calculations for most of the proteins.

The hydrophobic effect revisited--studies with supramolecular complexes imply high-energy water as a noncovalent driving force

Traditional descriptions of the hydrophobic effect on the basis of entropic arguments or the calculation of solvent-occupied surfaces must be questioned in view of new results obtained with supramolecular complexes. In these studies, it was possible to separate hydrophobic from dispersive interactions, which are strongest in aqueous systems. Even very hydrophobic alkanes associate significantly only in cavities containing water molecules with an insufficient number of possible hydrogen bonds. The replacement of high-energy water in cavities by guest molecules is the essential enthalpic driving force for complexation, as borne out by data for complexes of cyclodextrins, cyclophanes, and cucurbiturils, for which complexation enthalpies of up to -100 kJ mol(-1) were reached for encapsulated alkyl residues. Water-box simulations were used to characterize the different contributions from high-energy water and enabled the calculation of the association free enthalpies for selected cucurbituril complexes to within a 10% deviation from experimental values. Cavities in artificial receptors are more apt to show the enthalpic effect of high-energy water than those in proteins or nucleic acids, because they bear fewer or no functional groups in the inner cavity to stabilize interior water molecules.

Strategies to calculate water binding free energies in protein-ligand complexes

Water molecules are commonplace in protein binding pockets, where they can typically form a complex between the protein and a ligand or become displaced upon ligand binding. As a result, it is often of great interest to establish both the binding free energy and location of such molecules. Several approaches to predicting the location and affinity of water molecules to proteins have been proposed and utilized in the literature, although it is often unclear which method should be used under what circumstances. We report here a comparison between three such methodologies, Just Add Water Molecules (JAWS), Grand Canonical Monte Carlo (GCMC), and double-decoupling, in the hope of understanding the advantages and limitations of each method when applied to enclosed binding sites. As a result, we have adapted the JAWS scoring procedure, allowing the binding free energies of strongly bound water molecules to be calculated to a high degree of accuracy, requiring significantly less computational effort than more rigorous approaches. The combination of JAWS and GCMC offers a route to a rapid scheme capable of both locating and scoring water molecules for rational drug design.

Enhanced Monte Carlo Sampling through Replica Exchange with Solute Tempering

With a view to improving the consistency of free energy perturbation calculations in Monte Carlo simulations of protein–ligand complexes, we have implemented the replica exchange with solute tempering (REST) method in the MCPRO software. By augmenting the standard REST approach with regular attempted jumps in selected dihedral angles, our combined method facilitates sampling of ligand binding modes that are separated by high free energy barriers and ensures that computed free energy changes are considerably less dependent on the starting conditions and the chosen mutation pathway than those calculated with standard Monte Carlo sampling. We have applied the enhanced sampling method to the calculation of the activities of seven non-nucleoside inhibitors of HIV-1 reverse transcriptase, and its Tyr181Cys variant, and have shown that a range of binding orientations is possible depending on the nature of the ligand and the presence of mutations at the binding site.

Lead optimization mapper: automating free energy calculations for lead optimization

Alchemical free energy calculations hold increasing promise as an aid to drug discovery efforts. However, applications of these techniques in discovery projects have been relatively few, partly because of the difficulty of planning and setting up calculations. Here, we introduce lead optimization mapper, LOMAP, an automated algorithm to plan efficient relative free energy calculations between potential ligands within a substantial library of perhaps hundreds of compounds. In this approach, ligands are first grouped by structural similarity primarily based on the size of a (loosely defined) maximal common substructure, and then calculations are planned within and between sets of structurally related compounds. An emphasis is placed on ensuring that relative free energies can be obtained between any pair of compounds without combining the results of too many different relative free energy calculations (to avoid accumulation of error) and by providing some redundancy to allow for the possibility of error and consistency checking and provide some insight into when results can be expected to be unreliable. The algorithm is discussed in detail and a Python implementation, based on both Schrödinger's and OpenEye's APIs, has been made available freely under the BSD license.

Roles for ordered and bulk solvent in ligand recognition and docking in two related cavities

A key challenge in structure-based discovery is accounting for modulation of protein-ligand interactions by ordered and bulk solvent. To investigate this, we compared ligand binding to a buried cavity in Cytochrome c Peroxidase (CcP), where affinity is dominated by a single ionic interaction, versus a cavity variant partly opened to solvent by loop deletion. This opening had unexpected effects on ligand orientation, affinity, and ordered water structure. Some ligands lost over ten-fold in affinity and reoriented in the cavity, while others retained their geometries, formed new interactions with water networks, and improved affinity. To test our ability to discover new ligands against this opened site prospectively, a 534,000 fragment library was docked against the open cavity using two models of ligand solvation. Using an older solvation model that prioritized many neutral molecules, three such uncharged docking hits were tested, none of which was observed to bind; these molecules were not highly ranked by the new, context-dependent solvation score. Using this new method, another 15 highly-ranked molecules were tested for binding. In contrast to the previous result, 14 of these bound detectably, with affinities ranging from 8 µM to 2 mM. In crystal structures, four of these new ligands superposed well with the docking predictions but two did not, reflecting unanticipated interactions with newly ordered waters molecules. Comparing recognition between this open cavity and its buried analog begins to isolate the roles of ordered solvent in a system that lends itself readily to prospective testing and that may be broadly useful to the community.

Large scale affinity calculations of cyclodextrin host-guest complexes: Understanding the role of reorganization in the molecular recognition process

Host-guest inclusion complexes are useful models for understanding the structural and energetic aspects of molecular recognition. Due to their small size relative to much larger protein-ligand complexes, converged results can be obtained rapidly for these systems thus offering the opportunity to more reliably study fundamental aspects of the thermodynamics of binding. In this work, we have performed a large scale binding affinity survey of 57 β-cyclodextrin (CD) host guest systems using the binding energy distribution analysis method (BEDAM) with implicit solvation (OPLS-AA/AGBNP2). Converged estimates of the standard binding free energies are obtained for these systems by employing techniques such as parallel Hamitionian replica exchange molecular dynamics, conformational reservoirs and multistate free energy estimators. Good agreement with experimental measurements is obtained in terms of both numerical accuracy and affinity rankings. Overall, average effective binding energies reproduce affinity rank ordering better than the calculated binding affinities, even though calculated binding free energies, which account for effects such as conformational strain and entropy loss upon binding, provide lower root mean square errors when compared to measurements. Interestingly, we find that binding free energies are superior rank order predictors for a large subset containing the most flexible guests. The results indicate that, while challenging, accurate modeling of reorganization effects can lead to ligand design models of superior predictive power for rank ordering relative to models based only on ligand-receptor interaction energies.