Non-equilibrium approach for binding free energies in cyclodextrins in SAMPL7: force fields and software

Fast free energy estimates from λ-dynamics with bias-updated Gibbs sampling

Relative binding free energy calculations have become an integral computational tool for lead optimization in structure-based drug design. Classical alchemical methods, including free energy perturbation or thermodynamic integration, compute relative free energy differences by transforming one molecule into another. However, these methods have high operational costs due to the need to perform many pairwise perturbations independently. To reduce costs and accelerate molecular design workflows, we present a method called λ-dynamics with bias-updated Gibbs sampling. This method uses dynamic biases to continuously sample between multiple ligand analogues collectively within a single simulation. We show that many relative binding free energies can be determined quickly with this approach without compromising accuracy. For five benchmark systems, agreement to experiment is high, with root mean square errors near or below 1.0 kcal mol-1. Free energy results are consistent with other computational approaches and within statistical noise of both methods (0.4 kcal mol-1 or less). Notably, large efficiency gains over thermodynamic integration of 18-66-fold for small perturbations and 100-200-fold for whole aromatic ring substitutions are observed. The rapid determination of relative binding free energies will enable larger chemical spaces to be more readily explored and structure-based drug design to be accelerated.

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.

Host Dynamics under General-Purpose Force Fields

Macrocyclic hosts as prototypical receptors to gaseous and drug-like guests are crucial components in pharmaceutical research. The external guests are often coordinated at the center of these macromolecular containers. The formation of host-guest coordination is accompanied by the broken of host-water and host-ion interactions and sometimes also involves some conformational rearrangements of the host. A balanced description of various components of interacting terms is indispensable. However, up to now, the modeling community still lacks a general yet detailed understanding of commonly employed general-purpose force fields and the host dynamics produced by these popular selections. To fill this critical gap, in this paper, we profile the energetics and dynamics of four types of popular macrocycles, including cucurbiturils, pillararenes, cyclodextrins, and octa acids. The presented investigations of force field definitions, refitting, and evaluations are unprecedently detailed. Based on the valuable observations and insightful explanations, we finally summarize some general guidelines on force field parametrization and selection in host-guest modeling.

Comprehensive Evaluation of End-Point Free Energy Techniques in Carboxylated-Pillar[6]arene Host–Guest Binding: III. Force-Field Comparison, Three-Trajectory Realization and Further Dielectric Augmentation

Host–guest binding, despite the relatively simple structural and chemical features of individual components, still poses a challenge in computational modelling. The extreme underperformance of standard end-point methods in host–guest binding makes them practically useless. In the current work, we explore a potentially promising modification of the three-trajectory realization. The alteration couples the binding-induced structural reorganization into free energy estimation and suffers from dramatic fluctuations in internal energies in protein–ligand situations. Fortunately, the relatively small size of host–guest systems minimizes the magnitude of internal fluctuations and makes the three-trajectory realization practically suitable. Due to the incorporation of intra-molecular interactions in free energy estimation, a strong dependence on the force field parameters could be incurred. Thus, a term-specific investigation of transferable GAFF derivatives is presented, and noticeable differences in many aspects are identified between commonly applied GAFF and GAFF2. These force-field differences lead to different dynamic behaviors of the macrocyclic host, which ultimately would influence the end-point sampling and binding thermodynamics. Therefore, the three-trajectory end-point free energy calculations are performed with both GAFF versions. Additionally, due to the noticeable differences between host dynamics under GAFF and GAFF2, we add additional benchmarks of the single-trajectory end-point calculations. When only the ranks of binding affinities are pursued, the three-trajectory realization performs very well, comparable to and even better than the regressed PBSA_E scoring function and the dielectric constant-variable regime. With the GAFF parameter set, the TIP3P water in explicit solvent sampling and either PB or GB implicit solvent model in free energy estimation, the predictive power of the three-trajectory realization in ranking calculations surpasses all existing end-point methods on this dataset. We further combine the three-trajectory realization with another promising modified end-point regime of varying the interior dielectric constant. The combined regime does not incur sizable improvements for ranks and deviations from experiment exhibit non-monotonic variations.

Comprehensive evaluation of end-point free energy techniques in carboxylated-pillar[6]arene host-guest binding: II. regression and dielectric constant

End-point free energy calculations as a powerful tool have been widely applied in protein-ligand and protein-protein interactions. It is often recognized that these end-point techniques serve as an option of intermediate accuracy and computational cost compared with more rigorous statistical mechanic models (e.g., alchemical transformation) and coarser molecular docking. However, it is observed that this intermediate level of accuracy does not hold in relatively simple and prototypical host-guest systems. Specifically, in our previous work investigating a set of carboxylated-pillar[6]arene host-guest complexes, end-point methods provide free energy estimates deviating significantly from the experimental reference, and the rank of binding affinities is also incorrectly computed. These observations suggest the unsuitability and inapplicability of standard end-point free energy techniques in host-guest systems, and alteration and development are required to make them practically usable. In this work, we consider two ways to improve the performance of end-point techniques. The first one is the PBSA_E regression that varies the weights of different free energy terms in the end-point calculation procedure, while the second one is considering the interior dielectric constant as an additional variable in the end-point equation. By detailed investigation of the calculation procedure and the simulation outcome, we prove that these two treatments (i.e., regression and dielectric constant) are manipulating the end-point equation in a somehow similar way, i.e., weakening the electrostatic contribution and strengthening the non-polar terms, although there are still many detailed differences between these two methods. With the trained end-point scheme, the RMSE of the computed affinities is improved from the standard ~ 12 kcal/mol to ~ 2.4 kcal/mol, which is comparable to another altered end-point method (ELIE) trained with system-specific data. By tuning PBSA_E weighting factors with the host-specific data, it is possible to further decrease the prediction error to ~ 2.1 kcal/mol. These observations along with the extremely efficient optimized-structure computation procedure suggest the regression (i.e., PBSA_E as well as its GBSA_E extension) as a practically applicable solution that brings end-point methods back into the library of usable tools for host-guest binding. However, the dielectric-constant-variable scheme cannot effectively minimize the experiment-calculation discrepancy for absolute binding affinities, but is able to improve the calculation of affinity ranks. This phenomenon is somehow different from the protein-ligand case and suggests the difference between host-guest and biomacromolecular (protein-ligand and protein-protein) systems. Therefore, the spectrum of tools usable for protein-ligand complexes could be unsuitable for host-guest binding, and numerical validations are necessary to screen out really workable solutions in these 'prototypical' situations.

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.

Evaluating the use of absolute binding free energy in the fragment optimisation process

Key to the fragment optimisation process within drug design is the need to accurately capture the changes in affinity that are associated with a given set of chemical modifications. Due to the weakly binding nature of fragments, this has proven to be a challenging task, despite recent advancements in leveraging experimental and computational methods. In this work, we evaluate the use of Absolute Binding Free Energy (ABFE) calculations in guiding fragment optimisation decisions, retrospectively calculating binding free energies for 59 ligands across 4 fragment elaboration campaigns. We first demonstrate that ABFEs can be used to accurately rank fragment-sized binders with an overall Spearman's r of 0.89 and a Kendall τ of 0.67, although often deviating from experiment in absolute free energy values with an RMSE of 2.75 kcal/mol. We then also show that in several cases, retrospective fragment optimisation decisions can be supported by the ABFE calculations. Comparing against cheaper endpoint methods, namely N_wat-MM/GBSA, we find that ABFEs offer better ranking power and correlation metrics. Our results indicate that ABFE calculations can usefully guide fragment elaborations to maximise affinity.

Alchemical absolute protein-ligand binding free energies for drug design

The recent advances in relative protein-ligand binding free energy calculations have shown the value of alchemical methods in drug discovery. Accurately assessing absolute binding free energies, although highly desired, remains a challenging endeavour, mostly limited to small model cases. Here, we demonstrate accurate first principles based absolute binding free energy estimates for 128 pharmaceutically relevant targets. We use a novel rigorous method to generate protein-ligand ensembles for the ligand in its decoupled state. Not only do the calculations deliver accurate protein-ligand binding affinity estimates, but they also provide detailed physical insight into the structural determinants of binding. We identify subtle rotamer rearrangements between apo and holo states of a protein that are crucial for binding. When compared to relative binding free energy calculations, obtaining absolute binding free energies is considerably more challenging in large part due to the need to explicitly account for the protein in its apo state. In this work we present several approaches to obtain apo state ensembles for accurate absolute ΔG calculations, thus outlining protocols for prospective application of the methods for drug discovery.

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.

Binding of Inhibitors to BACE1 Affected by pH-Dependent Protonation: An Exploration from Multiple Replica Gaussian Accelerated Molecular Dynamics and MM-GBSA Calculations

To date, inhibiting the activity of β-amyloid cleaving enzyme 1 (BACE1) has been considered an efficient approach for treating Alzheimer's disease (AD). In the current work, multiple replica Gaussian accelerated molecular dynamics (MR-GaMD) simulations and the molecular mechanics general Born surface area (MM-GBSA) method were combined to investigate the effect of pH-dependent protonation on the binding of the inhibitors CS9, C6U, and 6WE to BACE1. Dynamic analyses based on the MR-GaMD trajectory show that pH-dependent protonation strongly affects the structural flexibility, correlated motions, and dynamic behavior of inhibitor-bound BACE1. According to the constructed free energy profiles, in the protonated state at low pH, inhibitor-bound BACE1 tends to populate at more conformations than in high pH. The binding free energies calculated by MM-GBSA suggest that inhibitors possess stronger binding abilities under the protonation conditions at high pH than under the protonation conditions at low pH. Moreover, pH-dependent protonation exerts a significant effect on the hydrogen bonding interactions of CS9, C6U, and 6WE to BACE1, which correspondingly alters the binding abilities of the three inhibitors to BACE1. Furthermore, in different protonated environments, three inhibitors share common interaction clusters and similar binding sites in BACE1, which are reliably used as efficient targets for the design of potent inhibitors of BACE1.

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.

Accurate absolute free energies for ligand-protein binding based on non-equilibrium approaches

The accurate calculation of the binding free energy for arbitrary ligand-protein pairs is a considerable challenge in computer-aided drug discovery. Recently, it has been demonstrated that current state-of-the-art molecular dynamics (MD) based methods are capable of making highly accurate predictions. Conventional MD-based approaches rely on the first principles of statistical mechanics and assume equilibrium sampling of the phase space. In the current work we demonstrate that accurate absolute binding free energies (ABFE) can also be obtained via theoretically rigorous non-equilibrium approaches. Our investigation of ligands binding to bromodomains and T4 lysozyme reveals that both equilibrium and non-equilibrium approaches converge to the same results. The non-equilibrium approach achieves the same level of accuracy and convergence as an equilibrium free energy perturbation (FEP) method enhanced by Hamiltonian replica exchange. We also compare uni- and bi-directional non-equilibrium approaches and demonstrate that considering the work distributions from both forward and reverse directions provides substantial accuracy gains. In summary, non-equilibrium ABFE calculations are shown to yield reliable and well-converged estimates of protein-ligand binding affinity.

Challenges Encountered Applying Equilibrium and Nonequilibrium Binding Free Energy Calculations

Binding free energy calculations have become increasingly valuable to drive decision making in drug discovery projects. However, among other issues, inadequate sampling can reduce accuracy, limiting the value of the technique. In this paper, we apply absolute binding free energy calculations to ligands binding to T4 lysozyme L99A and HSP90 using equilibrium and nonequilibrium approaches. We highlight sampling problems encountered in these systems, such as slow side chain rearrangements and slow changes of water placement upon ligand binding. These same types of challenges are also likely to show up in other protein-ligand systems, and we propose some strategies to diagnose and test for such problems in alchemical free energy calculations. We also explore similarities and differences in how the equilibrium and the nonequilibrium approaches handle these problems. Our results show the large amount of work still to be done to make free energy calculations robust and reliable and provide insight for future research in this area.

Doxorubicin Encapsulation in Carbon Nanotubes Having Haeckelite or Stone-Wales Defects as Drug Carriers: A Molecular Dynamics Approach

Doxorubicin (DOX), a recognized anticancer drug, forms stable associations with carbon nanotubes (CNTs). CNTs when properly functionalized have the ability to anchor directly in cancerous tumors where the release of the drug occurs thanks to the tumor slightly acidic pH. Herein, we study the armchair and zigzag CNTs with Stone-Wales (SW) defects to rank their ability to encapsulate DOX by determining the DOX-CNT binding free energies using the MM/PBSA and MM/GBSA methods implemented in AMBER16. We investigate also the chiral CNTs with haeckelite defects. Each haeckelite defect consists of a pair of square and octagonal rings. The armchair and zigzag CNT with SW defects and chiral nanotubes with haeckelite defects predict DOX-CNT interactions that depend on the length of the nanotube, the number of present defects and nitrogen doping. Chiral nanotubes having two haeckelite defects reveal a clear dependence on the nitrogen content with DOX-CNT interaction forces decreasing in the order 0N > 4N > 8N. These results contribute to a further understanding of drug-nanotube interactions and to the design of new drug delivery systems based on CNTs.

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.

SAMPL7 blind predictions using nonequilibrium alchemical approaches

In the context of the SAMPL7 challenge, we computed, employing a non-equilibrium (NE) alchemical technique, the standard binding free energy of two series of host-guest systems, involving as a host the Isaac's TrimerTrip, a Cucurbituril-like open cavitand, and the Gilson's Cyclodextrin derivatives. The adopted NE alchemy combines enhanced sampling molecular dynamics simulations with driven fast out-of-equilibrium alchemical trajectories to recover the free energy via the Jarzynski and Crooks NE theorems. The GAFF2 non-polarizable force field was used for the parametrization. Performances were acceptable and similar in accuracy to those we submitted for Gibb's Deep Cavity Cavitands in the previous SAMPL6 host-guest challenge, confirming the reliability of the computational approach and exposing, in some cases, some important deficiencies of the GAFF2 non-polarizable force field.

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.