Prediction of CB[8] host-guest binding free energies in SAMPL6 using the double-decoupling method

Expanded ensemble predictions of absolute binding free energies in the SAMPL9 host-guest challenge

As part of the SAMPL9 community-wide blind host-guest challenge, we implemented an expanded ensemble workflow to predict absolute binding free energies for 13 small molecules against pillar[6]arene. Notable features of our protocol include consideration of a variety of protonation and enantiomeric states for both host and guests, optimization of alchemical intermediates, and analysis of free energy estimates and their uncertainty using large numbers of simulation replicates performed using distributed computing. Our predictions of absolute binding free energies resulted in a mean absolute error of 2.29 kcal mol-1 and an R2 of 0.54. Overall, results show that expanded ensemble calculations using all-atom molecular dynamics simulations are a valuable and efficient computational tool in predicting absolute binding free energies.

Obtaining QM/MM binding free energies in the SAMPL8 drugs of abuse challenge: indirect approaches

Accurately predicting free energy differences is essential in realizing the full potential of rational drug design. Unfortunately, high levels of accuracy often require computationally expensive QM/MM Hamiltonians. Fortuitously, the cost of employing QM/MM approaches in rigorous free energy simulation can be reduced through the use of the so-called “indirect” approach to QM/MM free energies, in which the need for QM/MM simulations is avoided via a QM/MM “correction” at the classical endpoints of interest. Herein, we focus on the computation of QM/MM binding free energies in the context of the SAMPL8 Drugs of Abuse host–guest challenge. Of the 5 QM/MM correction coupled with force-matching submissions, PM6-D3H4/MM ranked submission proved the best overall QM/MM entry, with an RMSE from experimental results of 2.43 kcal/mol (best in ranked submissions), a Pearson’s correlation of 0.78 (second-best in ranked submissions), and a Kendall \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau$$\end{document}τ correlation of 0.52 (best in ranked submissions).

Recent Developments in Free Energy Calculations for Drug Discovery

The grand challenge in structure-based drug design is achieving accurate prediction of binding free energies. Molecular dynamics (MD) simulations enable modeling of conformational changes critical to the binding process, leading to calculation of thermodynamic quantities involved in estimation of binding affinities. With recent advancements in computing capability and predictive accuracy, MD based virtual screening has progressed from the domain of theoretical attempts to real application in drug development. Approaches including the Molecular Mechanics Poisson Boltzmann Surface Area (MM-PBSA), Linear Interaction Energy (LIE), and alchemical methods have been broadly applied to model molecular recognition for drug discovery and lead optimization. Here we review the varied methodology of these approaches, developments enhancing simulation efficiency and reliability, remaining challenges hindering predictive performance, and applications to problems in the fields of medicine and biochemistry.

Machine Learning and Enhanced Sampling Simulations for Computing the Potential of Mean Force and Standard Binding Free Energy

Computational capabilities are rapidly increasing, primarily because of the availability of GPU-based architectures. This creates unprecedented simulative possibilities for the systematic and robust computation of thermodynamic observables, including the free energy of a drug binding to a target. In contrast to calculations of relative binding free energy, which are nowadays widely exploited for drug discovery, we here push the boundary of computing the binding free energy and the potential of mean force. We introduce a novel protocol that leverages enhanced sampling, machine learning, and ad hoc algorithms to limit human intervention, computing time, and free parameters in free energy calculations. We first validate the method on a host-guest system, and then we apply the protocol to glycogen synthase kinase 3 beta, a protein kinase of pharmacological interest. Overall, we obtain a good correlation with experimental values in relative and absolute terms. While we focus on protein-ligand binding, the strategy is of broad applicability to any complex event that can be described with a path collective variable. We systematically discuss key details that influence the final result. The parameters and simulation settings are available at PLUMED-NEST to allow full reproducibility.

A replica exchange umbrella sampling (REUS) approach to predict host-guest binding free energies in SAMPL8 challenge

In this study, we report binding free energy calculations of various drugs-of-abuse to Cucurbit-[8]-uril as part of the SAMPL8 blind challenge. Force-field parameters were obtained from force-matching with different quantum mechanical levels of theory. The Replica Exchange Umbrella Sampling (REUS) approach was used with a cylindrical restraint to enhance the sampling of host–guest binding. Binding free energy was calculated by pulling the guest molecule from one side of the symmetric and cylindrical host, then into and through the host, and out the other side (bidirectional) as compared to pulling only to the bound pose inside the cylindrical host (unidirectional). The initial results with force-matched MP2 parameter set led to RMSE of 4.68 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text{kcal}}/{\text{mol}}$$\end{document}kcal/mol from experimental values. However, the follow-up study with CHARMM generalized force field parameters and force-matched PM6-D3H4 parameters resulted in RMSEs from experiment of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.65$$\end{document}2.65 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.72 {\text{kcal}}/{\text{mol}}$$\end{document}1.72kcal/mol, respectively, which demonstrates the potential of REUS for accurate binding free energy calculation given a more suitable description of energetics. Moreover, we compared the free energies for the so called bidirectional and unidirectional free energy approach and found that the binding free energies were highly similar. However, one issue in the bidirectional approach is the asymmetry of profile on the two sides of the host. This is mainly due to the insufficient sampling for these larger systems and can be avoided by longer sampling simulations. Overall REUS shows great promise for binding free energy calculations.

Evaluating the Performance of Water Models with Host-Guest Force Fields in Binding Enthalpy Calculations for Cucurbit[7]uril-Guest Systems

Computational prediction of thermodynamic components with computational methods has become increasingly routine in computer-aided drug design. Although there has been significant recent effort and improvements in the calculation of free energy, the prediction of enthalpy (and entropy) remains underexplored. Furthermore, there has been relatively little work reported so far that attempts to comparatively assess how well different force fields and water models perform in conjunction with each other. Here, we report a comprehensive assessment of force fields and water models using host-guest systems that mimic many features of protein-ligand systems. These systems are computationally inexpensive, possibly because of their small size compared to protein-ligand systems. We present absolute enthalpy calculations using the multibox approach on a set of 25 cucurbit[7]uril-guest pairs. Eight water models were considered (TIP3P, TIP4P, TIP4P-Ew, SPC, SPC/E, OPC, TIP5P, Bind3P), along with five force fields commonly used in the literature (GAFFv1, GAFFv2, CGenFF, Parsley, and SwissParam). We observe that host-guest binding enthalpies are strongly sensitive to the selection of force field and water model. In terms of water models, we find that TIP3P and its derivative Bind3P are the best performing models for this particular host-guest system. The performance is generally better for aliphatic compounds than for aromatic ones, suggesting that aromaticity remains a difficult property to include accurately in these simple force fields.

Binding Thermodynamics of Host-Guest Systems with SMIRNOFF99Frosst 1.0.5 from the Open Force Field Initiative

Designing ligands that bind their target biomolecules with high affinity and specificity is a key step in small-molecule drug discovery, but accurately predicting protein-ligand binding free energies remains challenging. Key sources of errors in the calculations include inadequate sampling of conformational space, ambiguous protonation states, and errors in force fields. Noncovalent complexes between a host molecule with a binding cavity and a druglike guest molecule have emerged as powerful model systems. As model systems, host-guest complexes reduce many of the errors in more complex protein-ligand binding systems, as their small size greatly facilitates conformational sampling, and one can choose systems that avoid ambiguities in protonation states. These features, combined with their ease of experimental characterization, make host-guest systems ideal model systems to test and ultimately optimize force fields in the context of binding thermodynamics calculations. The Open Force Field Initiative aims to create a modern, open software infrastructure for automatically generating and assessing force fields using data sets. The first force field to arise out of this effort, named SMIRNOFF99Frosst, has approximately one tenth the number of parameters, in version 1.0.5, compared to typical general small molecule force fields, such as GAFF. Here, we evaluate the accuracy of this initial force field, using free energy calculations of 43 α and β-cyclodextrin host-guest pairs for which experimental thermodynamic data are available, and compare with matched calculations using two versions of GAFF. For all three force fields, we used TIP3P water and AM1-BCC charges. The calculations are performed using the attach-pull-release (APR) method as implemented in the open source package, pAPRika. For binding free energies, the root-mean-square error of the SMIRNOFF99Frosst calculations relative to experiment is 0.9 [0.7, 1.1] kcal/mol, while the corresponding results for GAFF 1.7 and GAFF 2.1 are 0.9 [0.7, 1.1] kcal/mol and 1.7 [1.5, 1.9] kcal/mol, respectively, with 95% confidence ranges in brackets. These results suggest that SMIRNOFF99Frosst performs competitively with existing small molecule force fields and is a parsimonious starting point for optimization.

Use of Interaction Energies in QM/MM Free Energy Simulations

The use of the most accurate (i.e., QM or QM/MM) levels of theory for free energy simulations (FES) is typically not possible. Primarily, this is because the computational cost associated with the extensive configurational sampling needed for converging FES is prohibitive. To ensure the feasibility of QM-based FES, the "indirect" approach is generally taken, necessitating a free energy calculation between the MM and QM/MM potential energy surfaces. Ideally, this step is performed with standard free energy perturbation (Zwanzig's equation) as it only requires simulations be carried out at the low level of theory; however, work from several groups over the past few years has conclusively shown that Zwanzig's equation is ill-suited to this task. As such, many approximations have arisen to mitigate difficulties with Zwanzig's equation. One particularly popular notion is that the convergence of Zwanzig's equation can be improved by using interaction energy differences instead of total energy differences. Although problematic numerical fluctuations (a major problem when using Zwanzig's equation) are indeed reduced, our results and analysis demonstrate that this "interaction energy approximation" (IEA) is theoretically incorrect, and the implicit approximation invoked is spurious at best. Herein, we demonstrate this via solvation free energy calculations using IEA from two different low levels of theory to the same target high level. Results from this proof-of-concept consistently yield the wrong results, deviating by ∼1.5 kcal/mol from the rigorously obtained value.

Accelerating QM/MM Free Energy Computations via Intramolecular Force Matching

The calculation of free energy differences between levels of theory has numerous potential pitfalls. Chief among them is the lack of overlap, i.e., ensembles generated at one level of theory (e.g., "low") not being good approximations of ensembles at the other (e.g., "high"). Numerous strategies have been devised to mitigate this issue. However, the most straightforward approach is to ensure that the "low" level ensemble more closely resembles that of the "high". Ideally, this is done without increasing computational cost. Herein, we demonstrate that by reparametrizing classical intramolecular potentials to reproduce high level forces (i.e., force matching) configurational overlap between a "low" (i.e., classical) and "high" (i.e., quantum) level can be significantly improved. This procedure is validated on two test cases and results in vastly improved convergence of free energy simulations.

Comparison of the umbrella sampling and the double decoupling method in binding free energy predictions for SAMPL6 octa-acid host-guest challenges

We calculate the absolute binding free energies of tetra-methylated octa-acids host-guest systems as a part of the SAMPL6 blind challenge (receipt ID vq30p). We employed two different free energy simulation methods, i.e., the umbrella sampling (US) and double decoupling method (DDM). The US method was used with the weighted histogram analysis method (WHAM) (US-WHAM scheme). In the DDM scheme, Hamiltonian replica-exchange method (HREM) was combined with the Bennett acceptance ratio (BAR) (HREM-BAR scheme). We obtained initial binding poses via molecular docking using GalaxyDock-HG program, which is developed for the SAMPL challenge. The root mean square deviation (RMSD) and the mean absolute deviations (MAD) using US-WHAM scheme were 1.33 and 1.02 kcal/mol, respectively. The MAD was the top among all submissions, however the correlation with respect to experiment was unexceptional. While the RMSD and MAD via HREM-BAR scheme were greater than US-WHAM scheme, (i.e., 2.09 and 1.76 kcal/mol), their correlations were slightly better than US-WHAM. The correlation between the two methods was high. Further discussion on the DDM method can be found in a companion paper by Han et al. (receipt ID 3z83m) in the same issue.

Force matching as a stepping stone to QM/MM CB[8] host/guest binding free energies: a SAMPL6 cautionary tale

Use of quantum mechanical/molecular mechanical (QM/MM) methods in binding free energy calculations, particularly in the SAMPL challenge, often fail to achieve improvement over standard additive (MM) force fields. Frequently, the implementation is through use of reference potentials, or the so-called "indirect approach", and inherently relies on sufficient overlap existing between MM and QM/MM configurational spaces. This overlap is generally poor, particularly for the use of free energy perturbation to perform the MM to QM/MM free energy correction at the end states of interest (e.g., bound and unbound states). However, by utilizing MM parameters that best reproduce forces obtained at the desired QM level of theory, it is possible to lessen the configurational disparity between MM and QM/MM. To this end, we sought to use force matching to generate MM parameters for the SAMPL6 CB[8] host-guest binding challenge, classically compute binding free energies, and apply energetic end state corrections to obtain QM/MM binding free energy differences. For the standard set of 11 molecules and the bonus set (including three additional challenge molecules), error statistics, such as the root mean square deviation (RMSE) were moderately poor (5.5 and 5.4 kcal/mol). Correlation statistics, however, were in the top two for both standard and bonus set submissions ([Formula: see text] of 0.42 and 0.26, [Formula: see text] of 0.64 and 0.47 respectively). High RMSE and moderate correlation strongly indicated the presence of systematic error. Identifiable issues were ameliorated for two of the guest molecules, resulting in a reduction of error and pointing to strong prospects for the future use of this methodology.

Overview of the SAMPL6 host-guest binding affinity prediction challenge

Accurately predicting the binding affinities of small organic molecules to biological macromolecules can greatly accelerate drug discovery by reducing the number of compounds that must be synthesized to realize desired potency and selectivity goals. Unfortunately, the process of assessing the accuracy of current computational approaches to affinity prediction against binding data to biological macromolecules is frustrated by several challenges, such as slow conformational dynamics, multiple titratable groups, and the lack of high-quality blinded datasets. Over the last several SAMPL blind challenge exercises, host-guest systems have emerged as a practical and effective way to circumvent these challenges in assessing the predictive performance of current-generation quantitative modeling tools, while still providing systems capable of possessing tight binding affinities. Here, we present an overview of the SAMPL6 host-guest binding affinity prediction challenge, which featured three supramolecular hosts: octa-acid (OA), the closely related tetra-endo-methyl-octa-acid (TEMOA), and cucurbit[8]uril (CB8), along with 21 small organic guest molecules. A total of 119 entries were received from ten participating groups employing a variety of methods that spanned from electronic structure and movable type calculations in implicit solvent to alchemical and potential of mean force strategies using empirical force fields with explicit solvent models. While empirical models tended to obtain better performance than first-principle methods, it was not possible to identify a single approach that consistently provided superior results across all host-guest systems and statistical metrics. Moreover, the accuracy of the methodologies generally displayed a substantial dependence on the system considered, emphasizing the need for host diversity in blind evaluations. Several entries exploited previous experimental measurements of similar host-guest systems in an effort to improve their physical-based predictions via some manner of rudimentary machine learning; while this strategy succeeded in reducing systematic errors, it did not correspond to an improvement in statistical correlation. Comparison to previous rounds of the host-guest binding free energy challenge highlights an overall improvement in the correlation obtained by the affinity predictions for OA and TEMOA systems, but a surprising lack of improvement regarding root mean square error over the past several challenge rounds. The data suggests that further refinement of force field parameters, as well as improved treatment of chemical effects (e.g., buffer salt conditions, protonation states), may be required to further enhance predictive accuracy.