Blind prediction of host-guest binding affinities: a new SAMPL3 challenge

Host-Guest Binding Free Energies à la Carte: An Automated OneOPES Protocol

Estimating absolute binding free energies from molecular simulations is a key step in computer-aided drug design pipelines, but the agreement between computational results and experiments is still very inconsistent. Both the accuracy of the computational model and the quality of the statistical sampling contribute to this discrepancy, yet disentangling the two remains a challenge. In this study, we present an automated protocol based on OneOPES, an enhanced sampling method that exploits replica exchange and can accelerate several collective variables to address the sampling problem. We apply this protocol to 37 host-guest systems. The simplicity of setting up the simulations and producing well-converged binding free energy estimates without the need to optimize simulation parameters provides a reliable solution to the sampling problem. This, in turn, allows for a systematic force field comparison and ranking according to the correlation between simulations and experiments, which can inform the selection of an appropriate model. The protocol can be readily adapted to test more force field combinations and study more complex protein-ligand systems, where the choice of an appropriate physical model is often based on heuristic considerations rather than systematic optimization.

Calculation of protein-ligand binding entropies using a rule-based molecular fingerprint

The use of fast in silico prediction methods for protein-ligand binding free energies holds significant promise for the initial phases of drug development. Numerous traditional physics-based models (e.g., implicit solvent models), however, tend to either neglect or heavily approximate entropic contributions to binding due to their computational complexity. Consequently, such methods often yield imprecise assessments of binding strength. Machine learning models provide accurate predictions and can often outperform physics-based models. They, however, are often prone to overfitting, and the interpretation of their results can be difficult. Physics-guided machine learning models combine the consistency of physics-based models with the accuracy of modern data-driven algorithms. This work integrates physics-based model conformational entropies into a graph convolutional network. We introduce a new neural network architecture (a rule-based graph convolutional network) that generates molecular fingerprints according to predefined rules specifically optimized for binding free energy calculations. Our results on 100 small host-guest systems demonstrate significant improvements in convergence and preventing overfitting. We additionally demonstrate the transferability of our proposed hybrid model by training it on the aforementioned host-guest systems and then testing it on six unrelated protein-ligand systems. Our new model shows little difference in training set accuracy compared to a previous model but an order-of-magnitude improvement in test set accuracy. Finally, we show how the results of our hybrid model can be interpreted in a straightforward fashion.

The Open Force Field Initiative: Open Software and Open Science for Molecular Modeling

Force fields are a key component of physics-based molecular modeling, describing the energies and forces in a molecular system as a function of the positions of the atoms and molecules involved. Here, we provide a review and scientific status report on the work of the Open Force Field (OpenFF) Initiative, which focuses on the science, infrastructure and data required to build the next generation of biomolecular force fields. We introduce the OpenFF Initiative and the related OpenFF Consortium, describe its approach to force field development and software, and discuss accomplishments to date as well as future plans. OpenFF releases both software and data under open and permissive licensing agreements to enable rapid application, validation, extension, and modification of its force fields and software tools. We discuss lessons learned to date in this new approach to force field development. We also highlight ways that other force field researchers can get involved, as well as some recent successes of outside researchers taking advantage of OpenFF tools and data.

The SAMPL9 host-guest blind challenge: an overview of binding free energy predictive accuracy

We report the results of the SAMPL9 host-guest blind challenge for predicting binding free energies. The challenge focused on macrocycles from pillar[n]-arene and cyclodextrin host families, including WP6, and bCD and HbCD. A variety of methods were used by participants to submit binding free energy predictions. A machine learning approach based on molecular descriptors achieved the highest accuracy (RMSE of 2.04 kcal mol-1) among the ranked methods in the WP6 dataset. Interestingly, predictions for WP6 obtained via docking tended to outperform all methods (RMSE of 1.70 kcal mol-1), most of which are MD based and computationally more expensive. In general, methods applying force fields achieved better correlation with experiments for WP6 opposed to the machine learning and docking models. In the cyclodextrin-phenothiazine challenge, the ATM approach emerged as the top performing method with RMSE less than 1.86 kcal mol-1. Correlation metrics of ranked methods in this dataset were relatively poor compared to WP6. We also highlight several lessons learned to guide future work and help improve studies on the systems discussed. For example, WP6 may be present in other microstates other than its -12 state in the presence of certain guests. Machine learning approaches can be used to fine tune or help train force fields for certain chemistry (i.e. WP6-G4). Certain phenothiazines occupy distinct primary and secondary orientations, some of which were considered individually for accurate binding free energies. The accuracy of predictions from certain methods while starting from a single binding pose/orientation demonstrates the sensitivity of calculated binding free energies to the orientation, and in some cases the likely dominant orientation for the system. Computational and experimental results suggest that guest phenothiazine core traverses both the secondary and primary faces of the cyclodextrin hosts, a bulky cationic side chain will primarily occupy the primary face, and the phenothiazine core substituent resides at the larger secondary face.

Tuning Potential Functions to Host-Guest Binding Data

Software to more rapidly and accurately predict protein-ligand binding affinities is of high interest for early-stage drug discovery, and physics-based methods are among the most widely used technologies for this purpose. The accuracy of these methods depends critically on the accuracy of the potential functions that they use. Potential functions are typically trained against a combination of quantum chemical and experimental data. However, although binding affinities are among the most important quantities to predict, experimental binding affinities have not to date been integrated into the experimental data set used to train potential functions. In recent years, the use of host-guest complexes as simple and tractable models of binding thermodynamics has gained popularity due to their small size and simplicity, relative to protein-ligand systems. Host-guest complexes can also avoid ambiguities that arise in protein-ligand systems such as uncertain protonation states. Thus, experimental host-guest binding data are an appealing additional data type to integrate into the experimental data set used to optimize potential functions. Here, we report the extension of the Open Force Field Evaluator framework to enable the systematic calculation of host-guest binding free energies and their gradients with respect to force field parameters, coupled with the curation of 126 host-guest complexes with available experimental binding free energies. As an initial application of this novel infrastructure, we optimized generalized Born (GB) cavity radii for the OBC2 GB implicit solvent model against experimental data for 36 host-guest systems. This refitting led to a dramatic improvement in accuracy for both the training set and a separate test set with 90 additional host-guest systems. The optimized radii also showed encouraging transferability from host-guest systems to 59 protein-ligand systems. However, the new radii are significantly smaller than the baseline radii and lead to excessively favorable hydration free energies (HFEs). Thus, users of the OBC2 GB model currently may choose between GB cavity radii that yield more accurate binding affinities and GB cavity radii that yield more accurate HFEs. We suspect that achieving good accuracy on both will require more far-reaching adjustments to the GB model. We note that binding free-energy calculations using the OBC2 model in OpenMM gain about a 10× speedup relative to corresponding explicit solvent calculations, suggesting a future role for implicit solvent absolute binding free-energy (ABFE) calculations in virtual compound screening. This study proves the principle of using host-guest systems to train potential functions that are transferrable to protein-ligand systems and provides an infrastructure that enables a range of applications.

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.

An overview of the SAMPL8 host-guest binding challenge

The SAMPL series of challenges aim to focus the community on specific modeling challenges, while testing and hopefully driving progress of computational methods to help guide pharmaceutical drug discovery. In this study, we report on the results of the SAMPL8 host–guest blind challenge for predicting absolute binding affinities. SAMPL8 focused on two host–guest datasets, one involving the cucurbituril CB8 (with a series of common drugs of abuse) and another involving two different Gibb deep-cavity cavitands. The latter dataset involved a previously featured deep cavity cavitand (TEMOA) as well as a new variant (TEETOA), both binding to a series of relatively rigid fragment-like guests. Challenge participants employed a reasonably wide variety of methods, though many of these were based on molecular simulations, and predictive accuracy was mixed. As in some previous SAMPL iterations (SAMPL6 and SAMPL7), we found that one approach to achieve greater accuracy was to apply empirical corrections to the binding free energy predictions, taking advantage of prior data on binding to these hosts. Another approach which performed well was a hybrid MD-based approach with reweighting to a force matched QM potential. In the cavitand challenge, an alchemical method using the AMOEBA-polarizable force field achieved the best success with RMSE less than 1 kcal/mol, while another alchemical approach (ATM/GAFF2-AM1BCC/TIP3P/HREM) had RMSE less than 1.75 kcal/mol. The work discussed here also highlights several important lessons; for example, retrospective studies of reference calculations demonstrate the sensitivity of predicted binding free energies to ethyl group sampling and/or guest starting pose, providing guidance to help improve future studies on these systems.

Best practices for constructing, preparing, and evaluating protein-ligand binding affinity benchmarks [Article v0.1]

Free energy calculations are rapidly becoming indispensable in structure-enabled drug discovery programs. As new methods, force fields, and implementations are developed, assessing their expected accuracy on real-world systems (benchmarking) becomes critical to provide users with an assessment of the accuracy expected when these methods are applied within their domain of applicability, and developers with a way to assess the expected impact of new methodologies. These assessments require construction of a benchmark-a set of well-prepared, high quality systems with corresponding experimental measurements designed to ensure the resulting calculations provide a realistic assessment of expected performance when these methods are deployed within their domains of applicability. To date, the community has not yet adopted a common standardized benchmark, and existing benchmark reports suffer from a myriad of issues, including poor data quality, limited statistical power, and statistically deficient analyses, all of which can conspire to produce benchmarks that are poorly predictive of real-world performance. Here, we address these issues by presenting guidelines for (1) curating experimental data to develop meaningful benchmark sets, (2) preparing benchmark inputs according to best practices to facilitate widespread adoption, and (3) analysis of the resulting predictions to enable statistically meaningful comparisons among methods and force fields. We highlight challenges and open questions that remain to be solved in these areas, as well as recommendations for the collection of new datasets that might optimally serve to measure progress as methods become systematically more reliable. Finally, we provide a curated, versioned, open, standardized benchmark set adherent to these standards (PLBenchmarks) and an open source toolkit for implementing standardized best practices assessments (arsenic) for the community to use as a standardized assessment tool. While our main focus is free energy methods based on molecular simulations, these guidelines should prove useful for assessment of the rapidly growing field of machine learning methods for affinity prediction as well.

Nonbonded Force Field Parameters from Minimal Basis Iterative Stockholder Partitioning of the Molecular Electron Density Improve CB7 Host-Guest Affinity Predictions

Binding affinity prediction by means of computer simulation has been increasingly incorporated in drug discovery projects. Its wide application, however, is limited by the prediction accuracy of the free energy calculations. The main error sources are force fields used to describe molecular interactions and incomplete sampling of the configurational space. Organic host-guest systems have been used to address force field quality because they share similar interactions found in ligands and receptors, and their rigidity facilitates configurational sampling. Here, we test the binding free energy prediction accuracy for 14 guests with an aromatic or adamantane core and the CB7 host using molecular electron density derived nonbonded force field parameters. We developed a computational workflow written in Python to derive atomic charges and Lennard-Jones parameters with the Minimal Basis Iterative Stockholder method using the polarized electron density of several configurations of each guest in the bound and unbound states. The resulting nonbonded force field parameters improve binding affinity prediction, especially for guests with an adamantane core in which repulsive exchange and dispersion interactions to the host dominate.

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).

Thermodynamics of Pillararene•Guest Complexation: Blinded Dataset for the SAMPL9 Challenge

We report an investigation of the complexation between a water soluble pillararene host (WP6) and a panel of hydrophobic cationic guests (G1 - G20) by a combination of 1H NMR spectroscopy and isothermal titration calorimetry in phosphate buffered saline. We find that WP6 forms 1:1 complexes with Ka values in the 104 - 109 M-1 range driven by favorable enthalpic contributions. This thermodynamic dataset serves as blinded data for the SAMPL9 challenge.

Molecular Environment-Specific Atomic Charges Improve Binding Affinity Predictions of SAMPL5 Host-Guest Systems

Host-guest systems are widely used in benchmarks as model systems to improve computational methods for absolute binding free energy predictions. Recent advances in sampling algorithms for alchemical free energy calculations and the increase in computational power have made their binding affinity prediction primarily dependent on the quality of the force field. Here, we propose a new methodology to derive the atomic charges of host-guest systems based on quantum mechanics/molecular mechanics calculations and minimal basis iterative stockholder (MBIS) partitioning of the polarized electron density. A newly developed interface between the OpenMM and ORCA software packages provides D-MBIS charges that represent the guest's average electrostatic interactions in the hosts or the solvent. The simulation workflow also calculates the average energy required to polarize the guest in the bound and unbound state. Alchemical free energy calculations using the general Amber force field parameters with D-MBIS charges improve the binding affinity prediction of six guests bound to two octa acid hosts compared to the AM1-BCC charge set after correction with the average energetic polarization cost. This correction originates from the difference in potential energy that is required to polarize the guest in the bound and unbound state and contributes significantly to the binding affinity of anionic guests.

Energy-entropy method using multiscale cell correlation to calculate binding free energies in the SAMPL8 host-guest challenge

Free energy drives a wide range of molecular processes such as solvation, binding, chemical reactions and conformational change. Given the central importance of binding, a wide range of methods exist to calculate it, whether based on scoring functions, machine-learning, classical or electronic structure methods, alchemy, or explicit evaluation of energy and entropy. Here we present a new energy-entropy (EE) method to calculate the host-guest binding free energy directly from molecular dynamics (MD) simulation. Entropy is evaluated using Multiscale Cell Correlation (MCC) which uses force and torque covariance and contacts at two different length scales. The method is tested on a series of seven host-guest complexes in the SAMPL8 (Statistical Assessment of the Modeling of Proteins and Ligands) "Drugs of Abuse" Blind Challenge. The EE-MCC binding free energies are found to agree with experiment with an average error of 0.9 kcal mol-1. MCC makes clear the origin of the entropy changes, showing that the large loss of positional, orientational, and to a lesser extent conformational entropy of each binding guest is compensated for by a gain in orientational entropy of water released to bulk, combined with smaller decreases in vibrational entropy of the host, guest and contacting water.

SAMPL7 blind challenge: quantum-mechanical prediction of partition coefficients and acid dissociation constants for small drug-like molecules

The physicochemical properties of a drug molecule determine the therapeutic effectiveness of the drug. Thus, the development of fast and accurate theoretical approaches for the prediction of such properties is inevitable. The participation to the SAMPL7 challenge is based on the estimation of logP coefficients and pK_a values of small drug-like sulfonamide derivatives. Thereby, quantum mechanical calculations were carried out in order to calculate the free energy of solvation and the transfer energy of 22 drug-like compounds in different environments (water and n-octanol) by employing the SMD solvation model. For logP calculations, we studied eleven different methodologies to calculate the transfer free energies, the lowest RMSE value was obtained for the M06L/def2-TZVP//M06L/def2-SVP level of theory. On the other hand, we employed an isodesmic reaction scheme within the macro pK_a framework; this was based on selecting reference molecules similar to the SAMPL7 challenge molecules. Consequently, highly well correlated pK_a values were obtained with the M062X/6-311+G(2df,2p)//M052X/6-31+G(d,p) level of theory.

Evaluation of log P, pKa, and log D predictions from the SAMPL7 blind challenge

The Statistical Assessment of Modeling of Proteins and Ligands (SAMPL) challenges focuses the computational modeling community on areas in need of improvement for rational drug design. The SAMPL7 physical property challenge dealt with prediction of octanol-water partition coefficients and pK_a for 22 compounds. The dataset was composed of a series of N-acylsulfonamides and related bioisosteres. 17 research groups participated in the log P challenge, submitting 33 blind submissions total. For the pK_a challenge, 7 different groups participated, submitting 9 blind submissions in total. Overall, the accuracy of octanol-water log P predictions in the SAMPL7 challenge was lower than octanol-water log P predictions in SAMPL6, likely due to a more diverse dataset. Compared to the SAMPL6 pK_a challenge, accuracy remains unchanged in SAMPL7. Interestingly, here, though macroscopic pK_a values were often predicted with reasonable accuracy, there was dramatically more disagreement among participants as to which microscopic transitions produced these values (with methods often disagreeing even as to the sign of the free energy change associated with certain transitions), indicating far more work needs to be done on pK_a prediction methods.

Binding free energy predictions in host-guest systems using Autodock4. A retrospective analysis on SAMPL6, SAMPL7 and SAMPL8 challenges

We systematically tested the Autodock4 docking program for absolute binding free energy predictions using the host-guest systems from the recent SAMPL6, SAMPL7 and SAMPL8 challenges. We found that Autodock4 behaves surprisingly well, outperforming in many instances expensive molecular dynamics or quantum chemistry techniques, with an extremely favorable benefit-cost ratio. Some interesting features of Autodock4 predictions are revealed, yielding valuable hints on the overall reliability of docking screening campaigns in drug discovery projects.

Deciphering the Supramolecular Organization of Multiple Guests Inside a Microporous MOF to Understand their Release Profile

Metal-organic frameworks (MOFs) give the opportunity of confining guest molecules into their pores even by a post-synthetic protocol. PUM168 is a Zn-based MOF characterized by microporous cavities that allows the encapsulation of a significant number of guest molecules. The pores engineered with different binding sites show a remarkable guest affinity towards a series of natural essential oils components, such as eugenol, thymol and carvacrol, relevant for environmental applications. Exploiting single crystal X-ray diffraction, it was possible to step-wisely monitor the rather complex three-components guest exchange process involving dimethylformamide (DMF, the pristine solvent) and binary mixtures of the flavoring agents. A picture of the structural evolution of the DMF-to-guest replacement occurring inside the MOF crystal was reached by a detailed single-crystal-to-single-crystal monitoring. The relation of the supramolecular arrangement in the pores with selective guests release was then investigated as a function of time and temperature by static headspace GC-MS analysis.

Quantum-mechanical property prediction of solvated drug molecules: what have we learned from a decade of SAMPL blind prediction challenges?

Joint academic-industrial projects supporting drug discovery are frequently pursued to deploy and benchmark cutting-edge methodical developments from academia in a real-world industrial environment at different scales. The dimensionality of tasks ranges from small molecule physicochemical property assessment over protein-ligand interaction up to statistical analyses of biological data. This way, method development and usability both benefit from insights gained at both ends, when predictiveness and readiness of novel approaches are confirmed, but the pharmaceutical drug makers get early access to novel tools for the quality of drug products and benefit of patients. Quantum-mechanical and simulation methods particularly fall into this group of methods, as they require skills and expense in their development but also significant resources in their application, thus are comparatively slowly dripping into the realm of industrial use. Nevertheless, these physics-based methods are becoming more and more useful. Starting with a general overview of these and in particular quantum-mechanical methods for drug discovery we review a decade-long and ongoing collaboration between Sanofi and the Kast group focused on the application of the embedded cluster reference interaction site model (EC-RISM), a solvation model for quantum chemistry, to study small molecule chemistry in the context of joint participation in several SAMPL (Statistical Assessment of Modeling of Proteins and Ligands) blind prediction challenges. Starting with early application to tautomer equilibria in water (SAMPL2) the methodology was further developed to allow for challenge contributions related to predictions of distribution coefficients (SAMPL5) and acidity constants (SAMPL6) over the years. Particular emphasis is put on a frequently overlooked aspect of measuring the quality of models, namely the retrospective analysis of earlier datasets and predictions in light of more recent and advanced developments. We therefore demonstrate the performance of the current methodical state of the art as developed and optimized for the SAMPL6 pKa and octanol-water log P challenges when re-applied to the earlier SAMPL5 cyclohexane-water log D and SAMPL2 tautomer equilibria datasets. Systematic improvement is not consistently found throughout despite the similarity of the problem class, i.e. protonation reactions and phase distribution. Hence, it is possible to learn about hidden bias in model assessment, as results derived from more elaborate methods do not necessarily improve quantitative agreement. This indicates the role of chance or coincidence for model development on the one hand which allows for the identification of systematic error and opportunities toward improvement and reveals possible sources of experimental uncertainty on the other. These insights are particularly useful for further academia-industry collaborations, as both partners are then enabled to optimize both the computational and experimental settings for data generation.

Accurate Receptor-Ligand Binding Free Energies from Fast QM Conformational Chemical Space Sampling

Small molecule receptor-binding is dominated by weak, non-covalent interactions such as van-der-Waals hydrogen bonding or electrostatics. Calculating these non-covalent ligand-receptor interactions is a challenge to computational means in terms of accuracy and efficacy since the ligand may bind in a number of thermally accessible conformations. The conformational rotamer ensemble sampling tool (CREST) uses an iterative scheme to efficiently sample the conformational space and calculates energies using the semi-empirical 'Geometry, Frequency, Noncovalent, eXtended Tight Binding' (GFN2-xTB) method. This combined approach is applied to blind predictions of the modes and free energies of binding for a set of 10 drug molecule ligands to the cucurbit[n]urils CB[8] receptor from the recent 'Statistical Assessment of the Modeling of Proteins and Ligands' (SAMPL) challenge including morphine, hydromorphine, cocaine, fentanyl, and ketamine. For each system, the conformational space was sufficiently sampled for the free ligand and the ligand-receptor complexes using the quantum chemical Hamiltonian. A multitude of structures makes up the final conformer-rotamer ensemble, for which then free energies of binding are calculated. For those large and complex molecules, the results are in good agreement with experimental values with a mean error of 3 kcal/mol. The GFN2-xTB energies of binding are validated by advanced density functional theory calculations and found to be in good agreement. The efficacy of the automated QM sampling workflow allows the extension towards other complex molecular interaction scenarios.

Computation of host-guest binding free energies with a new quantum mechanics based mining minima algorithm

A new method called QM-VM2 is presented that efficiently combines statistical mechanics with quantum mechanical (QM) energy potentials in order to calculate noncovalent binding free energies of host-guest systems. QM-VM2 efficiently couples the use of semi-empirical QM (SEQM) energies and geometry optimizations with an underlying molecular mechanics (MM) based conformational search, to find low SEQM energy minima, and allows for processing of these minima at higher levels of ab initio QM theory. A progressive geometry optimization scheme is introduced as a means to increase conformational sampling efficiency. The newly implemented QM-VM2 is used to compute the binding free energies of the host molecule cucurbit[7]uril and a set of 15 guest molecules. The results are presented along with comparisons to experimentally determined binding affinities. For the full set of 15 host-guest complexes, which have a range of formal charges from +1 to +3, SEQM-VM2 based binding free energies show poor correlation with experiment, whereas for the ten +1 complexes only, a significant correlation (R2 = 0.8) is achieved. SEQM-VM2 generation of conformers followed by single-point ab initio QM calculations at the dispersion corrected restricted Hartree-Fock-D3(BJ) and TPSS-D3(BJ) levels of theory, as post-processing corrections, yields a reasonable correlation with experiment for the full set of host-guest complexes (R2 = 0.6 and R2 = 0.7, respectively) and an excellent correlation for the +1 formal charge set (R2 = 1.0 and R2 = 0.9, respectively), as long as a sufficiently large basis set (triple-zeta quality) is employed. The importance of the inclusion of configurational entropy, even at the MM level, for the achievement of good correlation with experiment was demonstrated by comparing the calculated ΔE values with experiment and finding a considerably poorer correlation with experiment than for the calculated free energy ΔE - TΔS. For the complete set of host-guest systems with the range of formal charges, it was observed that the deviation of the predicted binding free energy from experiment correlates somewhat with the net charge of the systems. This observation leads to a simple empirical interpolation scheme to improve the linear regression of the full set.

Testing automatic methods to predict free binding energy of host-guest complexes in SAMPL7 challenge

The design of new host–guest complexes represents a fundamental challenge in supramolecular chemistry. At the same time, it opens new opportunities in material sciences or biotechnological applications. A computational tool capable of automatically predicting the binding free energy of any host–guest complex would be a great aid in the design of new host systems, or to identify new guest molecules for a given host. We aim to build such a platform and have used the SAMPL7 challenge to test several methods and design a specific computational pipeline. Predictions will be based on machine learning (when previous knowledge is available) or a physics-based method (otherwise). The formerly delivered predictions with an RMSE of 1.67 kcal/mol but will require further work to identify when a specific system is outside of the scope of the model. The latter is combines the semiempirical GFN2B functional, with docking, molecular mechanics, and molecular dynamics. Correct predictions (RMSE of 1.45 kcal/mol) are contingent on the identification of the correct binding mode, which can be very challenging for host–guest systems with a large number of degrees of freedom. Participation in the blind SAMPL7 challenge provided fundamental direction to the project. More advanced versions of the pipeline will be tested against future SAMPL challenges.

SAMPL7 Host-Guest Challenge Overview: assessing the reliability of polarizable and non-polarizable methods for binding free energy calculations

The SAMPL challenges focus on testing and driving progress of computational methods to help guide pharmaceutical drug discovery. However, assessment of methods for predicting binding affinities is often hampered by computational challenges such as conformational sampling, protonation state uncertainties, variation in test sets selected, and even lack of high quality experimental data. SAMPL blind challenges have thus frequently included a component focusing on host-guest binding, which removes some of these challenges while still focusing on molecular recognition. Here, we report on the results of the SAMPL7 blind prediction challenge for host-guest affinity prediction. In this study, we focused on three different host-guest categories-a familiar deep cavity cavitand series which has been featured in several prior challenges (where we examine binding of a series of guests to two hosts), a new series of cyclodextrin derivatives which are monofunctionalized around the rim to add amino acid-like functionality (where we examine binding of two guests to a series of hosts), and binding of a series of guests to a new acyclic TrimerTrip host which is related to previous cucurbituril hosts. Many predictions used methods based on molecular simulations, and overall success was mixed, though several methods stood out. As in SAMPL6, we find that one strategy for achieving reasonable accuracy here was to make empirical corrections to binding predictions based on previous data for host categories which have been studied well before, though this can be of limited value when new systems are included. Additionally, we found that alchemical free energy methods using the AMOEBA polarizable force field had considerable success for the two host categories in which they participated. The new TrimerTrip system was also found to introduce some sampling problems, because multiple conformations may be relevant to binding and interconvert only slowly. Overall, results in this challenge tentatively suggest that further investigation of polarizable force fields for these challenges may be warranted.

AMOEBA binding free energies for the SAMPL7 TrimerTrip host-guest challenge

As part of the SAMPL7 host-guest binding challenge, the AMOEBA force field was applied to calculate the absolute binding free energy for 16 charged organic ammonium guests to the TrimerTrip host, a recently reported acyclic cucurbituril-derived clip host structure with triptycene moieties at its termini. Here we report binding free energy calculations for this system using the AMOEBA polarizable atomic multipole force field and double annihilation free energy methodology. Conformational analysis of the host suggests three families of conformations that do not interconvert in solution on a time scale available to nanosecond molecular dynamics (MD) simulations. Two of these host conformers, referred to as the "indent" and "overlap" structures, are capable of binding guest molecules. As a result, the free energies of all 16 guests binding to both conformations were computed separately, and combined to produce values for comparison with experiment. Initial ranked results submitted as part of the SAMPL7 exercise had a mean unsigned error (MUE) from experimental binding data of 2.14 kcal/mol. Subsequently, a rigorous umbrella sampling reference calculation was used to better determine the free energy difference between unligated "indent" and "overlap" host conformations. Revised binding values for the 16 guests pegged to this umbrella sampling reference reduced the MUE to 1.41 kcal/mol, with a correlation coefficient (Pearson R) between calculated and experimental binding values of 0.832 and a rank correlation (Kendall τ) of 0.65. Overall, the AMOEBA results demonstrate no significant systematic error, suggesting the force field provides an accurate energetic description of the TrimerTrip host, and an appropriate balance of solvation and desolvation effects associated with guest binding.

Separating Enthalpic, Configurational, and Solvation Entropic Components in Host-Guest Binding: Application to Cucurbit[7]uril Complexes through a Full In Silico Approach via Water Nanodroplets

Cucurbiturils are a family of supramolecular hosts obtained by condensation of glycoluril and formaldehyde. Cucurbit[7]uril, CB[7], is the most prominent member of the family for its biomolecular interest, arising from its mild solubility in water and for its strong binding with a large variety of guests containing nonpolar fragments such as adamantanes and ferrocene. For instance, CB[7] encapsulates diamantane diammonium iodide with an attomolar dissociation constant, a value unmatched even in natural encapsulation processes. Computational chemistry has been extensively employed to describe the enthalpic-entropic compensation principle of the molecular recognition process of cucurbituril hosts, but the synergistic contribution of experimental data is required for accurate results to be obtained. This paper proposes the first fully theoretical model able to reconcile the calculated thermodynamics of the complexation process with the experimental data obtained by calorimetry (ITC) for cucurbit[7]uril. The model allows the isolation and estimation of all of the enthalpic and entropic contributions coming from solute and solvent alike to the whole host-guest binding event and enables the straightforward calculation of the contribution of the solvation entropy to the binding.

Toward Understanding CB[7]-Based Supramolecular Diels-Alder Catalysis

Cucurbiturils (CBs) are robust and versatile macrocyclic compounds, often used as molecular hosts in complex supramolecular systems. In previous work, remarkable catalytic activity has been observed for asymmetric cycloadditions under very mild conditions. Herein, we investigate the nature of supramolecular catalysis using DFT calculations and QM/MM techniques. We discuss induced conformational changes, electrostatic shielding effects from the highly polar aqueous environment and cooperativity in hydrogen bonding of the substrates in explicit water using QM/MM simulation techniques. Our results show little specificity for the chosen molecules, suggesting an excellent opportunity to expand the scope for catalytic use of these supramolecular macrocyclic containers.

SAMPL7: Host-guest binding prediction by molecular dynamics and quantum mechanics

Statistical Assessment of Modeling of Proteins and Ligands (SAMPL) challenges provide routes to compare chemical quantities determined using computational chemistry approaches to experimental measurements that are shared after the competition. For this effort, several computational methods have been used to calculate the binding energies of Octa Acid (OA) and exo-Octa Acid (exoOA) host-guest systems for SAMPL7. The initial poses for molecular dynamics (MD) were generated by molecular docking. Binding free energy calculations were performed using molecular mechanics combined with Poisson-Boltzmann or generalized Born surface area solvation (MMPBSA/MMGBSA) approaches. The factors that affect the utility of the MMPBSA/MMGBSA approaches including solvation, partial charge, and solute entropy models were also analyzed. In addition to MD calculations, quantum mechanics (QM) calculations were performed using several different density functional theory (DFT) approaches. From SAMPL6 results, B3PW91-D3 was found to overestimate binding energies though it was effective for geometry optimizations, so it was considered for the DFT geometry optimizations in the current study, with single-point energy calculations carried out with B2PLYP-D3 with double-, triple-, and quadruple-ζ level basis sets. Accounting for dispersion effects, and solvation models was deemed essential for the predictions. MMGBSA and MMPBSA correlated better to experiment when used in conjunction with an empirical/linear correction.

Ligand Gaussian Accelerated Molecular Dynamics (LiGaMD): Characterization of Ligand Binding Thermodynamics and Kinetics

Calculations of ligand binding free energies and kinetic rates are important for drug design. However, such tasks have proven challenging in computational chemistry and biophysics. To address this challenge, we have developed a new computational method, ligand Gaussian accelerated molecular dynamics (LiGaMD), which selectively boosts the ligand nonbonded interaction potential energy based on the Gaussian accelerated molecular dynamics (GaMD) enhanced sampling technique. Another boost potential could be applied to the remaining potential energy of the entire system in a dual-boost algorithm (LiGaMD_Dual) to facilitate ligand binding. LiGaMD has been demonstrated on host-guest and protein-ligand binding model systems. Repetitive guest binding and unbinding in the β-cyclodextrin host were observed in hundreds-of-nanosecond LiGaMD_Dual simulations. The calculated guest binding free energies agreed excellently with experimental data with <1.0 kcal/mol errors. Compared with converged microsecond-time scale conventional molecular dynamics simulations, the sampling errors of LiGaMD_Dual simulations were also <1.0 kcal/mol. Accelerations of ligand kinetic rate constants in LiGaMD simulations were properly estimated using Kramers' rate theory. Furthermore, LiGaMD allowed us to capture repetitive dissociation and binding of the benzamidine inhibitor in trypsin within 1 μs simulations. The calculated ligand binding free energy and kinetic rate constants compared well with the experimental data. In summary, LiGaMD provides a powerful enhanced sampling approach for characterizing ligand binding thermodynamics and kinetics simultaneously, which is expected to facilitate computer-aided drug design.

Multi-phase Boltzmann weighting: accounting for local inhomogeneity in molecular simulations of water-octanol partition coefficients in the SAMPL6 challenge

Accurately computing partition coefficients is a pivotal part of drug discovery. Specifically, octanol-water partition coefficients can provide information into hydrophobicity of drug-like molecules, as well as a de facto representation of membrane permeability. However, one challenge facing the computation of partition coefficients is the need to encapsulate various microscopic environments. These include areas of largely bulk solvent (i.e., either water or octanol) or regions where octanol is saturated with water or areas of higher salt concentration. Also, tautomeric effects require consideration. Thus, we present a Boltzmann weighting approach that incorporates transfer free energies across varying microscopic media, as well as varying tautomeric state, to compute partition coefficients in the SAMPL6 challenge.

Data-Driven Mapping of Gas-Phase Quantum Calculations to General Force Field Lennard-Jones Parameters

Molecular dynamics simulations are helpful tools for a range of applications, ranging from drug discovery to protein structure determination. The successful use of this technology largely depends on the potential function, or force field, used to determine the potential energy at each configuration of the system. Most force fields encode all of the relevant parameters to be used in distinct atom types, each associated with parameters for all parts of the force field, typically bond stretches, angle bends, torsions, and nonbonded terms accounting for van der Waals and electrostatic interactions. Much attention has been paid to the nonbonded parameters and their derivation, which are important in particular due to their governance of noncovalent interactions, such as protein-ligand binding. Parametrization involves adjusting the nonbonded parameters to minimize the error between simulation results and experimental properties, such as heats of vaporization and densities of neat liquids. In this setting, determining the best set of atom types is far from trivial, and the large number of parameters to be fit for the atom types in a typical force field can make it difficult to approach a true optimum. Here, we utilize a previously described Minimal Basis Iterative Stockholder (MBIS) method to carry out an atoms-in-molecules partitioning of electron densities. Information from these atomic densities is then mapped to Lennard-Jones parameters using a set of mapping parameters much smaller than the typical number of atom types in a force field. This approach is advantageous in two ways: it eliminates atom types by allowing each atom to have unique Lennard-Jones parameters, and it greatly reduces the number of parameters to be optimized. We show that this approach yields results comparable to those obtained with the typed GAFF 1.7 force field, even when trained on a relatively small amount of experimental data.

Quantum chemical predictions of water-octanol partition coefficients applied to the SAMPL6 logP blind challenge

Theoretical approaches for predicting physicochemical properties are valuable tools for accelerating the drug discovery process. In this work, quantum chemical methods are used to predict water-octanol partition coefficients as a part of the SAMPL6 blind challenge. The SMD continuum solvent model was employed with MP2 and eight DFT functionals in conjunction with correlation consistent basis sets to determine the water-octanol transfer free energy. Several tactics towards improving the predictions of the partition coefficient were examined, including increasing the quality of basis sets, considering tautomerization, and accounting for inhomogeneities in the water and n-octanol phases. Evaluation of these various schemes highlights the impact of modeling approaches across different methods. With the inclusion of tautomers and adjustments to the permittivity constants, the best predictions were obtained with smaller basis sets and the O3LYP functional, which yielded an RMSE of 0.79 logP units. The results presented correspond to the SAMPL6 logP submission IDs: DYXBT, O7DJK, and AHMTF.

A blind SAMPL6 challenge: insight into the octanol-water partition coefficients of drug-like molecules via a DFT approach

In this study quantum mechanical methods were used to predict the solvation energies of a series of drug-like molecules both in water and in octanol, in the context of the SAMPL6 n-octanol/water partition coefficient challenge. In pharmaceutical design, n-octanol/water partition coefficient, LogP, describes the drug's hydrophobicity and membrane permeability, thus, a well-established theoretical method that rapidly determines the hydrophobicity of a drug, enables the progress of the drug design. In this study, the solvation free energies were obtained via six different methodologies (B3LYP, M06-2X and ωB97XD functionals with 6-311+G** and 6-31G* basis sets) by taking into account the environment implicitly; the methodology chosen (B3LYP/6-311+G**) was used later to evaluate ΔG_solv by using explicit water as solvent. We optimized each conformer in different solvents separately, our calculations have shown that the stability of the conformers is highly dependent on the solvent environment. We have compared implicitly and explicitly solvated systems, the interaction of one explicit water with drug-molecules at the proper location leads to the prediction of more accurate LogP values.

Solvation Free Energy of Self-Assembled Complexes: Using Molecular Dynamics to Understand the Separation of ssDNA-Wrapped Single-Walled Carbon Nanotubes

Molecular dynamics simulations were used to characterize the self-assembly of single-stranded DNA (ssDNA) on a (6,5) single-walled carbon nanotube (SWCNT) in aqueous solution for the purpose of gaining an improved theoretical understanding of separation strategies for SWCNTs using ssDNA as a dispersant. Four separate ssDNA sequences, ((TAT)4, TTA(TAT)2ATT, C12, (GTC)2GT), at various levels of loading, were chosen for study based on published experimental work showing selective extraction of particular SWCNT species based on the ssDNA dispersant sequence. We develop a unique workflow based on free energy perturbation (FEP) and use this to determine the relative solubility of these complexes due to the adsorption of the ssDNA on the SWCNT surface, and hence, rank the favorability of separations observed during experiments. Results qualitatively agree with experiments and indicate that the nucleobase sequence of the adsorbed ssDNA greatly affects the free energy of complex solvation which ultimately drives SWCNT separation. Further, to elucidate the underlying physics governing the ssDNA-SWCNT solubility rankings, we also present calculations for four structural characteristics of ssDNA adsorption. We demonstrate that a unique type of intra-strand hydrogen bonding is the most important factor contributing to the stability of the ssDNA-SWCNT complexes and show how these adsorption characteristics are coupled with the FEP results.

Octanol-water partition coefficient measurements for the SAMPL6 blind prediction challenge

Partition coefficients describe the equilibrium partitioning of a single, defined charge state of a solute between two liquid phases in contact, typically a neutral solute. Octanol-water partition coefficients ([Formula: see text]), or their logarithms (log P), are frequently used as a measure of lipophilicity in drug discovery. The partition coefficient is a physicochemical property that captures the thermodynamics of relative solvation between aqueous and nonpolar phases, and therefore provides an excellent test for physics-based computational models that predict properties of pharmaceutical relevance such as protein-ligand binding affinities or hydration/solvation free energies. The SAMPL6 Part II octanol-water partition coefficient prediction challenge used a subset of kinase inhibitor fragment-like compounds from the SAMPL6 [Formula: see text] prediction challenge in a blind experimental benchmark. Following experimental data collection, the partition coefficient dataset was kept blinded until all predictions were collected from participating computational chemistry groups. A total of 91 submissions were received from 27 participating research groups. This paper presents the octanol-water log P dataset for this SAMPL6 Part II partition coefficient challenge, which consisted of 11 compounds (six 4-aminoquinazolines, two benzimidazole, one pyrazolo[3,4-d]pyrimidine, one pyridine, one 2-oxoquinoline substructure containing compounds) with log P values in the range of 1.95-4.09. We describe the potentiometric log P measurement protocol used to collect this dataset using a Sirius T3, discuss the limitations of this experimental approach, and share suggestions for future log P data collection efforts for the evaluation of computational methods.

Triptycene Walled Glycoluril Trimer: Synthesis and Recognition Properties

We report the synthesis of a new acyclic CB[n]-type host (1) that features a central glycoluril trimer capped by triptycene sidewalls. Host 1 has good solubility in water (≈ 3 mM) and does not undergo strong self-association (K_s = 480 M^-1). We probed the geometry of the complexes by analyzing the complexation induced changes in the ¹H NMR spectra and measured the complexation thermodynamics by isothermal titration calorimetry. The conformation of 1 and its packing in the solid state was revealed by single crystal x-ray diffraction measurements.

Continuous Evaluation of Ligand Protein Predictions: A Weekly Community Challenge for Drug Docking

Docking calculations can accelerate drug discovery by predicting the bound poses of ligands for a targeted protein. However, it is not clear which docking methods work best. Furthermore, predicting poses requires steps outside the docking algorithm itself, such as preparation of the protein and ligand, and it is not known which components are most in need of improvement. The Continuous Evaluation of Ligand Protein Predictions (CELPP) is a blinded prediction challenge designed to address these issues. Participants create a workflow to predict protein-ligand binding poses, which is then tasked with predicting 10-100 new protein-ligand crystal structures each week. CELPP evaluates the accuracy of each workflow's predictions and posts the scores online. The results can be used to identify the strengths and weaknesses of current approaches, help map docking problems to the algorithms most likely to overcome them, and illuminate areas of unmet need in structure-guided drug design.

Ultrahigh binding affinity of a hydrocarbon guest inside cucurbit[7]uril enhanced by strong host–guest charge matching

Electrostatic interactions as a strong driving force for neutral apolar host–guest complexation is revealed via a large-scale computational approach.

Toward Achieving Efficient and Accurate Ligand-Protein Unbinding with Deep Learning and Molecular Dynamics through RAVE

In this work, we demonstrate how to leverage our recent iterative deep learning-all atom molecular dynamics (MD) technique "Reweighted autoencoded variational Bayes for enhanced sampling (RAVE)" (Ribeiro, Bravo, Wang, Tiwary, J. Chem. Phys. 2018, 149, 072301) for investigating ligand-protein unbinding mechanisms and calculating absolute binding free energies, Δ G_b, when plagued with difficult to sample rare events. In order to do so, we introduce a simple but powerful extension to RAVE that allows learning a reaction coordinate expressed as a piecewise function that is linear over all intervals. Such an approach allows us to retain the physical interpretation of a RAVE-derived reaction coordinate while making the method more applicable to a wider range of complex biophysical problems. As we will demonstrate, using as our test-case the slow dissociation of benzene from the L99A variant of lysozyme, the RAVE extension led to observing an unbinding event in 100% of the independent all-atom MD simulations, all within 3-50 ns for a process that takes on an average close to few hundred milliseconds, which reflects a 7 orders of magnitude acceleration relative to straightforward MD. Furthermore, we will show that without the use of time-dependent biasing, clear back-and-forth movement between metastable intermediates was achieved during the various simulations, demonstrating the caliber of the RAVE-derived piecewise reaction coordinate and bias potential, which together drive efficient and accurate sampling of the ligand-protein dissociation event. Last, we report the results for Δ G_b, which via very short MD simulations, can form a strict lower-bound that is ∼2-3 kcal/mol off from experiments. We believe that RAVE, together with its multidimensional extension that we introduce here, will be a useful tool for simulating the slow unbinding process of practical ligand-protein complexes in an automated manner with minimal use of human intuition.

Escaping Atom Types in Force Fields Using Direct Chemical Perception

Traditional approaches to specifying a molecular mechanics force field encode all the information needed to assign force field parameters to a given molecule into a discrete set of atom types. This is equivalent to a representation consisting of a molecular graph comprising a set of vertices, which represent atoms labeled by atom type, and unlabeled edges, which represent chemical bonds. Bond stretch, angle bend, and dihedral parameters are then assigned by looking up bonded pairs, triplets, and quartets of atom types in parameter tables to assign valence terms and using the atom types themselves to assign nonbonded parameters. This approach, which we call indirect chemical perception because it operates on the intermediate graph of atom-typed nodes, creates a number of technical problems. For example, atom types must be sufficiently complex to encode all necessary information about the molecular environment, making it difficult to extend force fields encoded this way. Atom typing also results in a proliferation of redundant parameters applied to chemically equivalent classes of valence terms, needlessly increasing force field complexity. Here, we describe a new approach to assigning force field parameters via direct chemical perception. Rather than working through the intermediary of the atom-typed graph, direct chemical perception operates directly on the unmodified chemical graph of the molecule to assign parameters. In particular, parameters are assigned to each type of force field term (e.g., bond stretch, angle bend, torsion, and Lennard-Jones) based on standard chemical substructure queries implemented via the industry-standard SMARTS chemical perception language, using SMIRKS extensions that permit labeling of specific atoms within a chemical pattern. We use this to implement a new force field format, called the SMIRKS Native Open Force Field (SMIRNOFF) format. We demonstrate the power and generality of this approach using examples of specific molecules that pose problems for indirect chemical perception and construct and validate a minimalist yet very general force field, SMIRNOFF99Frosst. We find that a parameter definition file only ∼300 lines long provides coverage of all but <0.02% of a 5 million molecule drug-like test set. Despite its simplicity, the accuracy of SMIRNOFF99Frosst for small molecule hydration free energies and selected properties of pure organic liquids is similar to that of the General Amber Force Field, whose specification requires thousands of parameters. This force field provides a starting point for further optimization and refitting work to follow.

pKa measurements for the SAMPL6 prediction challenge for a set of kinase inhibitor-like fragments

Determining the net charge and protonation states populated by a small molecule in an environment of interest or the cost of altering those protonation states upon transfer to another environment is a prerequisite for predicting its physicochemical and pharmaceutical properties. The environment of interest can be aqueous, an organic solvent, a protein binding site, or a lipid bilayer. Predicting the protonation state of a small molecule is essential to predicting its interactions with biological macromolecules using computational models. Incorrectly modeling the dominant protonation state, shifts in dominant protonation state, or the population of significant mixtures of protonation states can lead to large modeling errors that degrade the accuracy of physical modeling. Low accuracy hinders the use of physical modeling approaches for molecular design. For small molecules, the acid dissociation constant (pK_a) is the primary quantity needed to determine the ionic states populated by a molecule in an aqueous solution at a given pH. As a part of SAMPL6 community challenge, we organized a blind pK_a prediction component to assess the accuracy with which contemporary pK_a prediction methods can predict this quantity, with the ultimate aim of assessing the expected impact on modeling errors this would induce. While a multitude of approaches for predicting pK_a values currently exist, predicting the pK_as of drug-like molecules can be difficult due to challenging properties such as multiple titratable sites, heterocycles, and tautomerization. For this challenge, we focused on set of 24 small molecules selected to resemble selective kinase inhibitors-an important class of therapeutics replete with titratable moieties. Using a Sirius T3 instrument that performs automated acid-base titrations, we used UV absorbance-based pK_a measurements to construct a high-quality experimental reference dataset of macroscopic pK_as for the evaluation of computational pK_a prediction methodologies that was utilized in the SAMPL6 pK_a challenge. For several compounds in which the microscopic protonation states associated with macroscopic pK_as were ambiguous, we performed follow-up NMR experiments to disambiguate the microstates involved in the transition. This dataset provides a useful standard benchmark dataset for the evaluation of pK_a prediction methodologies on kinase inhibitor-like compounds.

Absolute binding free energies for the SAMPL6 cucurbit[8]uril host-guest challenge via the AMOEBA polarizable force field

As part of the SAMPL6 host-guest blind challenge, the AMOEBA force field was applied to calculate the absolute binding free energy for a cucurbit[8]uril host complexed with 14 diverse guests, ranging from small, rigid structures to drug molecules. The AMOEBA results from the initial submission prompted an investigation into aspects of the methodology and parameterization employed. Lessons learned from the blind challenge include: a double annihilation scheme (electrostatics and van der Waals) is needed to obtain proper sampling of guest conformations, annihilation of key torsion parameters of the guest are recommended for flexible guests, and a more thorough analysis of torsion parameters is warranted. When put in to practice with the AMOEBA model, the lessons learned improved the MUE from 2.63 to 1.20 kcal/mol and the RMSE from 3.62 to 1.68 kcal/mol, respectively. Overall, the AMOEBA protocol for determining absolute binding free energies benefitted from participation in the SAMPL6 host-guest blind challenge and the results suggest the implementation of the methodology in future host-guest calculations.

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.

SAMPL6 host-guest challenge: binding free energies via a multistep approach

In this effort in the SAMPL6 host-guest binding challenge, a combination of molecular dynamics and quantum mechanical methods were used to blindly predict the host-guest binding free energies of a series of cucurbit[8]uril (CB8), octa-acid (OA), and tetramethyl octa-acid (TEMOA) hosts bound to various guest molecules in aqueous solution. Poses for host-guest systems were generated via molecular dynamics (MD) simulations and clustering analyses. The binding free energies for the structures obtained via cluster analyses of MD trajectories were calculated using the MMPBSA method and density functional theory (DFT) with the inclusion of Grimme's dispersion correction, an implicit solvation model to model the aqueous solution, and the resolution-of-the-identity (RI) approximation (MMPBSA, RI-B3PW91-D3, and RI-B3PW91, respectively). Among these three methods tested, the results for OA and TEMOA systems showed MMPBSA and RI-B3PW91-D3 methods can be used to qualitatively rank binding energies of small molecules with an overbinding by 7 and 37 kcal/mol respectively, and RI-B3PW91 gave the poorest quality results, indicating the importance of dispersion correction for the binding free energy calculations. Due to the complexity of the CB8 systems, all of the methods tested show poor correlation with the experimental results. Other quantum mechanical approaches used for the calculation of binding free energies included DFT without the RI approximation, utilizing truncated basis sets to reduce the computational cost (memory, disk space, CPU time), and a corrected dielectric constant to account for ionic strength within the implicit solvation model.

Binding free energies in the SAMPL6 octa-acid host-guest challenge calculated with MM and QM methods

We have estimated free energies for the binding of eight carboxylate ligands to two variants of the octa-acid deep-cavity host in the SAMPL6 blind-test challenge (with or without endo methyl groups on the four upper-rim benzoate groups, OAM and OAH, respectively). We employed free-energy perturbation (FEP) for relative binding energies at the molecular mechanics (MM) and the combined quantum mechanical (QM) and MM (QM/MM) levels, the latter obtained with the reference-potential approach with QM/MM sampling for the MM → QM/MM FEP. The semiempirical QM method PM6-DH+ was employed for the ligand in the latter calculations. Moreover, binding free energies were also estimated from QM/MM optimised structures, combined with COSMO-RS estimates of the solvation energy and thermostatistical corrections from MM frequencies. They were performed at the PM6-DH+ level of theory with the full host and guest molecule in the QM system (and also four water molecules in the geometry optimisations) for 10-20 snapshots from molecular dynamics simulations of the complex. Finally, the structure with the lowest free energy was recalculated using the dispersion-corrected density-functional theory method TPSS-D3, for both the structure and the energy. The two FEP approaches gave similar results (PM6-DH+/MM slightly better for OAM), which were among the five submissions with the best performance in the challenge and gave the best results without any fit to data from the SAMPL5 challenge, with mean absolute deviations (MAD) of 2.4-5.2 kJ/mol and a correlation coefficient (R²) of 0.77-0.93. This is the first time QM/MM approaches give binding free energies that are competitive to those obtained with MM for the octa-acid host. The QM/MM-optimised structures gave somewhat worse performance (MAD = 3-8 kJ/mol and R² = 0.1-0.9), but the results were improved compared to previous studies of this system with similar methods.

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 ligand binding affinity using on- and off-rates for the SAMPL6 SAMPLing challenge

Interest in ligand binding kinetics has been growing rapidly, as it is being discovered in more and more systems that ligand residence time is the crucial factor governing drug efficacy. Many enhanced sampling methods have been developed with the goal of predicting ligand binding rates ([Formula: see text]) and/or ligand unbinding rates ([Formula: see text]) through explicit simulation of ligand binding pathways, and these methods work by very different mechanisms. Although there is not yet a blind challenge for ligand binding kinetics, here we take advantage of experimental measurements and rigorously computed benchmarks to compare estimates of [Formula: see text] calculated as the ratio of two rates: [Formula: see text]. These rates were determined using a new enhanced sampling method based on the weighted ensemble framework that we call "REVO": Reweighting of Ensembles by Variance Optimization. This is a further development of the WExplore enhanced sampling method, in which trajectory cloning and merging steps are guided not by the definition of sampling regions, but by maximizing trajectory variance. Here we obtain estimates of [Formula: see text] and [Formula: see text] that are consistent across multiple simulations, with an average log10-scale standard deviation of 0.28 for on-rates and 0.56 for off-rates, which is well within an order of magnitude and far better than previously observed for previous applications of the WExplore algorithm. Our rank ordering of the three host-guest pairs agrees with the reference calculations, however our predicted [Formula: see text] values were systematically lower than the reference by an average of 4.2 kcal/mol. Using tree network visualizations of the trajectories in the REVO algorithm, and conformation space networks for each system, we analyze the results of our sampling, and hypothesize sources of discrepancy between our [Formula: see text] values and the reference. We also motivate the direct inclusion of [Formula: see text] and [Formula: see text] challenges in future iterations of SAMPL, to further develop the field of ligand binding kinetics prediction and modeling.