Studies on the Configurational Stability of Tropolone-Ketone-, Ester-, and Aldehyde-Based Chiral Axes.

Recent studies have revealed that tropolone-amide aryl C-C(O) rotational barriers are dramatically higher than those of analogous benzamide-based systems, and as a result, they have an increased likelihood of displaying high configurational stability. Studies on other tropolone-based chiral axes are important to assess the generality of this phenomenon. Herein, we describe a series of studies on the rotational barriers of tropolone-ketone, tropolone-ester, and tropolone-aldehyde chiral axes. These studies are complemented with computational modeling of the dynamics of these and analogous benzenoid variants to illuminate the impact that tropolone may have on aryl-C(O) configurational stability.

Performance and Analysis of the Alchemical Transfer Method for Binding-Free-Energy Predictions of Diverse Ligands.

The Alchemical Transfer Method (ATM) is herein validated against the relative binding-free energies (RBFEs) of a diverse set of protein-ligand complexes. We employed a streamlined setup workflow, a bespoke force field, and AToM-OpenMM software to compute the RBFEs of the benchmark set prepared by Schindler and collaborators at Merck KGaA. This benchmark set includes examples of standard small R-group ligand modifications as well as more challenging scenarios, such as large R-group changes, scaffold hopping, formal charge changes, and charge-shifting transformations. The novel coordinate perturbation scheme and a dual-topology approach of ATM address some of the challenges of single-topology alchemical RBFE methods. Specifically, ATM eliminates the need for splitting electrostatic and Lennard-Jones interactions, atom mapping, defining ligand regions, and postcorrections for charge-changing perturbations. Thus, ATM is simpler and more broadly applicable than conventional alchemical methods, especially for scaffold-hopping and charge-changing transformations. Here, we performed well over 500 RBFE calculations for eight protein targets and found that ATM achieves accuracy comparable to that of existing state-of-the-art methods, albeit with larger statistical fluctuations. We discuss insights into the specific strengths and weaknesses of the ATM method that will inform future deployments. This study confirms that ATM can be applied as a production tool for RBFE predictions across a wide range of perturbation types within a unified, open-source framework.

Enhancing Protein-Ligand Binding Affinity Predictions Using Neural Network Potentials.

This letter gives results on improving protein-ligand binding affinity predictions based on molecular dynamics simulations using machine learning potentials with a hybrid neural network potential and molecular mechanics methodology (NNP/MM). We compute relative binding free energies with the Alchemical Transfer Method and validate its performance against established benchmarks and find significant enhancements compared with conventional MM force fields like GAFF2.

Validation of the Alchemical Transfer Method for the Estimation of Relative Binding Affinities of Molecular Series.

The accurate prediction of protein-ligand binding affinities is crucial for drug discovery. Alchemical free energy calculations have become a popular tool for this purpose. However, the accuracy and reliability of these methods can vary depending on the methodology. In this study, we evaluate the performance of a relative binding free energy protocol based on the alchemical transfer method (ATM), a novel approach based on a coordinate transformation that swaps the positions of two ligands. The results show that ATM matches the performance of more complex free energy perturbation (FEP) methods in terms of Pearson correlation but with marginally higher mean absolute errors. This study shows that the ATM method is competitive compared to more traditional methods in speed and accuracy and offers the advantage of being applicable with any potential energy function.

Taming multiple binding poses in alchemical binding free energy prediction: the ß-cyclodextrin host-guest SAMPL9 blinded challenge.

We apply the Alchemical Transfer Method (ATM) and a bespoke fixed partial charge force field to the SAMPL9 bCD host-guest binding free energy prediction challenge that comprises a combination of complexes formed between five phenothiazine guests and two cyclodextrin hosts. Multiple chemical forms, competing binding poses, and computational modeling challenges pose significant obstacles to obtaining reliable computational predictions for these systems. The phenothiazine guests exist in solution as racemic mixtures of enantiomers related by nitrogen inversions that bind the hosts in various binding poses, each requiring an individual free energy analysis. Due to the large size of the guests and the conformational reorganization of the hosts, which prevent a direct absolute binding free energy route, binding free energies are obtained by a series of absolute and relative binding alchemical steps for each chemical species in each binding pose. Metadynamics-accelerated conformational sampling was found to be necessary to address the poor convergence of some numerical estimates affected by conformational trapping. Despite these challenges, our blinded predictions quantitatively reproduced the experimental affinities for the ß-cyclodextrin host and, to a lesser extent, those with a methylated derivative. The work illustrates the challenges of obtaining reliable free energy data in in silico drug design for even seemingly simple systems and introduces some of the technologies available to tackle them.

Relative Binding Free Energy Calculations for Ligands with Diverse Scaffolds with the Alchemical Transfer Method.

We present an extension of the alchemical transfer method (ATM) for the estimation of relative binding free energies of molecular complexes applicable to conventional, as well as scaffold-hopping, alchemical transformations. Named ATM-RBFE, the method is implemented in the free and open-source OpenMM molecular simulation package and aims to provide a simpler and more generally applicable route to the calculation of relative binding free energies than what is currently available. ATM-RBFE is based on sound statistical mechanics theory and a novel coordinate perturbation scheme designed to swap the positions of a pair of ligands such that one is transferred from the bulk solvent to the receptor binding site while the other moves simultaneously in the opposite direction. The calculation is conducted directly in a single solvent box with a system prepared with conventional setup tools, without splitting of electrostatic and nonelectrostatic transformations, and without pairwise soft-core potentials. ATM-RBFE is validated here against the absolute binding free energies of the SAMPL8 GDCC host-guest benchmark set and against protein-ligand benchmark sets that include complexes of the estrogen receptor ERα and those of the methyltransferase EZH2. In each case the method yields self-consistent and converged relative binding free energy estimates in agreement with absolute binding free energies and reference literature values, as well as experimental measurements.

Application of the alchemical transfer and potential of mean force methods to the SAMPL8 host-guest blinded challenge.

We report the results of our participation in the SAMPL8 GDCC Blind Challenge for host-guest binding affinity predictions. Absolute binding affinity prediction is of central importance to the biophysics of molecular association and pharmaceutical discovery. The blinded SAMPL series have provided an important forum for assessing the reliability of binding free energy methods in an objective way. In this challenge, we employed two binding free energy methods, the newly developed alchemical transfer method (ATM) and the well-established potential of mean force (PMF) physical pathway method, using the same setup and force field model. The calculated binding free energies from the two methods are in excellent quantitative agreement. Importantly, the results from the two methods were also found to agree well with the experimental binding affinities released subsequently, with R values of 0.89 (ATM) and 0.83 (PMF). These results were ranked among the best of the SAMPL8 GDCC challenge and second only to those obtained with the more accurate AMOEBA force field. Interestingly, the two host molecules included in the challenge (TEMOA and TEETOA) displayed distinct binding mechanisms, with TEMOA undergoing a dehydration transition whereas guest binding to TEETOA resulted in the opening of the binding cavity that remains essentially dry during the process. The coupled reorganization and hydration equilibria observed in these systems is a useful prototype for the study of these phenomena often observed in the formation of protein-ligand complexes. Given that the two free energy methods employed here are based on entirely different thermodynamic pathways, the close agreement between the two and their general agreement with the experimental binding free energies are a testament to the high quality and precision achieved by theory and methods. The study provides further validation of the novel ATM binding free energy estimation protocol and paves the way to further extensions of the method to more complex systems.

Developing end-point methods for absolute binding free energy calculation using the Boltzmann-quasiharmonic model.

Understanding the physical forces underlying receptor-ligand binding requires robust methods for analyzing the binding thermodynamics. In end-point binding free energy methods the binding free energy is naturally decomposable into physically intuitive contributions such as the solvation free energy and configurational entropy that can provide insights. Here we present a new end-point method called EE-BQH (Effective Energy-Boltzmann-Quasiharmonic) which combines the Boltzmann-Quasiharmonic model for configurational entropy with different solvation free energy methods, such as the continuum solvent PBSA model and the integral equation-based 3D-RISM, to estimate the absolute binding free energy. We compare EE-BQH with other treatments of configurational entropy such as Quasiharmonic models in internal coordinates (QHIC) and in Cartesian coordinates (QHCC), and Normal Mode analysis (NMA), by testing them on the octa acids host-guest complexes from the SAMPL8 blind challenge. The accuracies in the calculated absolute binding free energies strongly depend on the configurational entropy and solvation free energy methods used. QHIC and BQH yield the best agreements with the established potential of mean force (PMF) estimates, with R2 of ∼0.7 and mean unsigned error of ∼1.7 kcal mol-1. These results from the end-point calculations are also in similar agreement with experiments. While 3D-RISM in combination with QHIC or BQH lead to reasonable correlations with the PMF results and experiments, the calculated absolute binding free energies are underestimated by ∼5 kcal mol-1. While the binding is accompanied by a significant reduction in the ligand translational/rotational entropy, the change in the torsional entropy in these host-guest systems is slightly positive. Compared with BQH, QHIC underestimates the reduction of configurational entropy because of the non-Gaussian probability distributions in the ligand rotation and a small number of torsions. The study highlights the crucial role of configurational entropy in determining binding and demonstrates the potential of using the new end-point method to provide insights in more complex protein-ligand systems.

New tetrahydroisoquinoline-based D3R ligands with an o-xylenyl linker motif.

The effect of rigidification of the n-butyl linker region of tetrahydroisoquinoline-containing D3R ligands via inclusion of an o-xylenyl motif was examined in this study. Generally, rigidification with an o-xylenyl linker group reduces D3R affinity and negatively impacts selectivity versus D2R for compounds possessing a 6-methoxy-1,2,3,4,-tetrahydroisoquinolin-7-ol primary pharmacophore group. However, D3R affinity appears to be regulated by the primary pharmacophore group and high affinity D3R ligands with 6,7-dihydroxy-1,2,3,4-tetrahydroisoquinoline and 6,7-dimethoxy-1,2,3,4-tetrahydroisoquinoline primary pharmacophore groups were identified. The results of this study also indicate that D3R selectivity versus the σ2R is dictated by the benzamide secondary pharmacophore group, this being facilitated with 4-substituted benzamides. Compounds 5s and 5t were identified as high affinity (Ki < 4 nM) D3R ligands. Docking studies revealed that the added phenyl ring moiety interacts with the Cys181 in D3R which partially accounts for the strong D3R affinity of the ligands.

Alchemical transformations for concerted hydration free energy estimation with explicit solvation.

We present a family of alchemical perturbation potentials that enable the calculation of hydration free energies of small- to medium-sized molecules in a single concerted alchemical coupling step instead of the commonly used sequence of two distinct coupling steps for Lennard-Jones and electrostatic interactions. The perturbation potentials we employ are non-linear functions of the solute-solvent interaction energy designed to focus sampling near entropic bottlenecks along the alchemical pathway. We present a general framework to optimize the parameters of alchemical perturbation potentials of this kind. The optimization procedure is based on the λ-function formalism and the maximum-likelihood parameter estimation procedure we developed earlier to avoid the occurrence of multi-modal distributions of the coupling energy along the alchemical path. A novel soft-core function applied to the overall solute-solvent interaction energy rather than individual interatomic pair potentials critical for this result is also presented. Because it does not require modifications of core force and energy routines, the soft-core formulation can be easily deployed in molecular dynamics simulation codes. We illustrate the method by applying it to the estimation of the hydration free energy in water droplets of compounds of varying size and complexity. In each case, we show that convergence of the hydration free energy is achieved rapidly. This work paves the way for the ongoing development of more streamlined algorithms to estimate free energies of molecular binding with explicit solvation.

Exploring the Free-Energy Landscape and Thermodynamics of Protein-Protein Association.

We report the free-energy landscape and thermodynamics of the protein-protein association responsible for the drug-induced multimerization of HIV-1 integrase (IN). Allosteric HIV-1 integrase inhibitors promote aberrant IN multimerization by bridging IN-IN intermolecular interactions. However, the thermodynamic driving forces and kinetics of the multimerization remain largely unknown. Here, we explore the early steps in the IN multimerization by using umbrella sampling and unbiased molecular dynamics simulations in explicit solvent. In direct simulations, the two initially separated dimers spontaneously associate to form near-native complexes that resemble the crystal structure of the aberrant tetramer. Most strikingly, the effective interaction of the protein-protein association is very short-ranged: the two dimers associate rapidly within tens of nanoseconds when their binding surfaces are separated by d ≤ 4.3 Å (less than two water diameters). Beyond this distance, the oligomerization kinetics appears to be diffusion controlled with a much longer association time. The free-energy profile also captured the crucial role of allosteric IN inhibitors in promoting multimerization and explained why several C-terminal domain mutations are remarkably resistant to the drug-induced multimerization. The results also show that at small separation, the protein-protein binding process contains two consecutive phases with distinct thermodynamic signatures. First, interprotein water molecules are expelled to the bulk, resulting in a small increase in entropy, as the solvent entropy gain from the water release is nearly cancelled by the loss of side-chain entropies as the two proteins approach each other. At shorter distances, the two dry binding surfaces adapt to each other to optimize their interaction energy at the expense of further protein configurational entropy loss. Although the binding interfaces feature clusters of hydrophobic residues, overall, the protein-protein association in this system is driven by enthalpy and opposed by entropy.

3,7-Dihydroxytropolones Inhibit Initiation of Hepatitis B Virus Minus-Strand DNA Synthesis.

Initiation of protein-primed (-) strand DNA synthesis in hepatitis B virus (HBV) requires interaction of the viral reverse transcriptase with epsilon (ε), a cis-acting regulatory signal located at the 5' terminus of pre-genomic RNA (pgRNA), and several host-encoded chaperone proteins. Binding of the viral polymerase (P protein) to ε is necessary for pgRNA encapsidation and synthesis of a short primer covalently attached to its terminal domain. Although we identified small molecules that recognize HBV ε RNA, these failed to inhibit protein-primed DNA synthesis. However, since initiation of HBV (-) strand DNA synthesis occurs within a complex of viral and host components (e.g., Hsp90, DDX3 and APOBEC3G), we considered an alternative therapeutic strategy of allosteric inhibition by disrupting the initiation complex or modifying its topology. To this end, we show here that 3,7-dihydroxytropolones (3,7-dHTs) can inhibit HBV protein-primed DNA synthesis. Since DNA polymerase activity of a ribonuclease (RNase H)-deficient HBV reverse transcriptase that otherwise retains DNA polymerase function is also abrogated, this eliminates direct involvement of RNase (ribonuclease) H activity of HBV reverse transcriptase and supports the notion that the HBV initiation complex might be therapeutically targeted. Modeling studies also provide a rationale for preferential activity of 3,7-dHTs over structurally related α-hydroxytropolones (α-HTs).

A grid-based algorithm in conjunction with a gaussian-based model of atoms for describing molecular geometry.

A novel grid-based method is presented, which in conjunction with a smooth Gaussian-based model of atoms, is used to compute molecular volume (MV) and surface area (MSA). The MV and MSA are essential for computing nonpolar component of free energies. The objective of our grid-based approach is to identify solute atom pairs that share overlapping volumes in space. Once completed, this information is used to construct a rooted tree using depth-first method to yield the final volume and SA by using the formulations of the Gaussian model described by Grant and Pickup (J. Phys Chem, 1995, 99, 3503). The method is designed to function uninterruptedly with the grid-based finite-difference method implemented in Delphi, a popular and open-source package used for solving the Poisson-Boltzmann equation (PBE). We demonstrate the time efficacy of the method while also validating its performance in terms of the effect of grid-resolution, positioning of the solute within the grid-map and accuracy in identification of overlapping atom pairs. We also explore and discuss different aspects of the Gaussian model with key emphasis on its physical meaningfulness. This development and its future release with the Delphi package are intended to provide a physically meaningful, fast, robust and comprehensive tool for MM/PBSA based free energy calculations. © 2019 Wiley Periodicals, Inc.

Massive-Scale Binding Free Energy Simulations of HIV Integrase Complexes Using Asynchronous Replica Exchange Framework Implemented on the IBM WCG Distributed Network.

To perform massive-scale replica exchange molecular dynamics (REMD) simulations for calculating binding free energies of protein-ligand complexes, we implemented the asynchronous replica exchange (AsyncRE) framework of the binding energy distribution analysis method (BEDAM) in implicit solvent on the IBM World Community Grid (WCG) and optimized the simulation parameters to reduce the overhead and improve the prediction power of the WCG AsyncRE simulations. We also performed the first massive-scale binding free energy calculations using the WCG distributed computing grid and 301 ligands from the SAMPL4 challenge for large-scale binding free energy predictions of HIV-1 integrase complexes. In total there are ∼10000 simulated complexes, ∼1 million replicas, and ∼2000 µs of aggregated MD simulations. Running AsyncRE MD simulations on the WCG requires accepting a trade-off between the number of replicas that can be run (breadth) and the number of full RE cycles that can be completed per replica (depth). As compared with synchronous Replica Exchange (SyncRE) running on tightly coupled clusters like XSEDE, on the WCG many more replicas can be launched simultaneously on heterogeneous distributed hardware, but each full RE cycle requires more overhead. We compared the WCG results with that from AutoDock and more advanced RE simulations including the use of flattening potentials to accelerate sampling of selected degrees of freedom of ligands and/or receptors related to slow dynamics due to high energy barriers. We propose a suitable strategy of RE simulations to refine high throughput docking results which can be matched to corresponding computing resources: from HPC clusters, to small or medium-size distributed campus grids, and finally to massive-scale computing networks including millions of CPUs like the resources available on the WCG.

Perturbation potentials to overcome order/disorder transitions in alchemical binding free energy calculations.

We investigate the role of order/disorder transitions in alchemical simulations of protein-ligand absolute binding free energies. We show, in the context of a potential of mean force description, that for a benchmarking system (the complex of the L99A mutant of T4 lysozyme with 3-iodotoluene) and for a more challenging system relevant for medicinal applications (the complex of the farnesoid X receptor with inhibitor 26 from a recent D3R challenge) that order/disorder transitions can significantly hamper Hamiltonian replica exchange sampling efficiency and slow down the rate of equilibration of binding free energy estimates. We further show that our analytical model of alchemical binding combined with the formalism developed by Straub et al. for the treatment of order/disorder transitions of molecular systems can be successfully employed to analyze the transitions and help design alchemical schedules and soft-core functions that avoid or reduce the adverse effects of rare binding/unbinding transitions. The results of this work pave the way for the application of these techniques to the alchemical estimation with explicit solvation of hydration free energies and absolute binding free energies of systems undergoing order/disorder transitions.

The generalized Boltzmann distribution is the only distribution in which the Gibbs-Shannon entropy equals the thermodynamic entropy.

We show that the generalized Boltzmann distribution is the only distribution for which the Gibbs-Shannon entropy equals the thermodynamic entropy. This result means that the thermodynamic entropy and the Gibbs-Shannon entropy are not generally equal, but rather the equality holds only in the special case where a system is in equilibrium with a reservoir.

Efficient gaussian density formulation of volume and surface areas of macromolecules on graphical processing units.

We present an algorithm to efficiently compute accurate volumes and surface areas of macromolecules on graphical processing unit (GPU) devices using an analytic model which represents atomic volumes by continuous Gaussian densities. The volume of the molecule is expressed by means of the inclusion-exclusion formula, which is based on the summation of overlap integrals among multiple atomic densities. The surface area of the molecule is obtained by differentiation of the molecular volume with respect to atomic radii. The many-body nature of the model makes a port to GPU devices challenging. To our knowledge, this is the first reported full implementation of this model on GPU hardware. To accomplish this, we have used recursive strategies to construct the tree of overlaps and to accumulate volumes and their gradients on the tree data structures so as to minimize memory contention. The algorithm is used in the formulation of a surface area-based non-polar implicit solvent model implemented as an open source plug-in (named GaussVol) for the popular OpenMM library for molecular mechanics modeling. GaussVol is 50 to 100 times faster than our best optimized implementation for the CPUs, achieving speeds in excess of 100 ns/day with 1 fs time-step for protein-sized systems on commodity GPUs. © 2017 Wiley Periodicals, Inc.

A combined treatment of hydration and dynamical effects for the modeling of host-guest binding thermodynamics: the SAMPL5 blinded challenge.

As part of the SAMPL5 blinded experiment, we computed the absolute binding free energies of 22 host-guest complexes employing a novel approach based on the BEDAM single-decoupling alchemical free energy protocol with parallel replica exchange conformational sampling and the AGBNP2 implicit solvation model specifically customized to treat the effect of water displacement as modeled by the Hydration Site Analysis method with explicit solvation. Initial predictions were affected by the lack of treatment of ionic charge screening, which is very significant for these highly charged hosts, and resulted in poor relative ranking of negatively versus positively charged guests. Binding free energies obtained with Debye-Hückel treatment of salt effects were in good agreement with experimental measurements. Water displacement effects contributed favorably and very significantly to the observed binding affinities; without it, the modeling predictions would have grossly underestimated binding. The work validates the implicit/explicit solvation approach employed here and it shows that comprehensive physical models can be effective at predicting binding affinities of molecular complexes requiring accurate treatment of conformational dynamics and hydration.

Parameterization of an effective potential for protein-ligand binding from host-guest affinity data.

Force field accuracy is still one of the "stalemates" in biomolecular modeling. Model systems with high quality experimental data are valuable instruments for the validation and improvement of effective potentials. With respect to protein-ligand binding, organic host-guest complexes have long served as models for both experimental and computational studies because of the abundance of binding affinity data available for such systems. Binding affinity data collected for cyclodextrin (CD) inclusion complexes, a popular model for molecular recognition, is potentially a more reliable resource for tuning energy parameters than hydration free energy measurements. Convergence of binding free energy calculations on CD host-guest systems can also be obtained rapidly, thus offering the opportunity to assess the robustness of these parameters. In this work, we demonstrate how implicit solvent parameters can be developed using binding affinity experimental data and the binding energy distribution analysis method (BEDAM) and validated using the Grid Inhomogeneous Solvation Theory analysis. These new solvation parameters were used to study protein-ligand binding in two drug targets against the HIV-1 virus and improved the agreement between the calculated and the experimental binding affinities. This work illustrates how benchmark sets of high quality experimental binding affinity data and physics-based binding free energy models can be used to evaluate and optimize force fields for protein-ligand systems. Copyright © 2015 John Wiley & Sons, Ltd.

Large scale free energy calculations for blind predictions of protein-ligand binding: the D3R Grand Challenge 2015.

We describe binding free energy calculations in the D3R Grand Challenge 2015 for blind prediction of the binding affinities of 180 ligands to Hsp90. The present D3R challenge was built around experimental datasets involving Heat shock protein (Hsp) 90, an ATP-dependent molecular chaperone which is an important anticancer drug target. The Hsp90 ATP binding site is known to be a challenging target for accurate calculations of ligand binding affinities because of the ligand-dependent conformational changes in the binding site, the presence of ordered waters and the broad chemical diversity of ligands that can bind at this site. Our primary focus here is to distinguish binders from nonbinders. Large scale absolute binding free energy calculations that cover over 3000 protein-ligand complexes were performed using the BEDAM method starting from docked structures generated by Glide docking. Although the ligand dataset in this study resembles an intermediate to late stage lead optimization project while the BEDAM method is mainly developed for early stage virtual screening of hit molecules, the BEDAM binding free energy scoring has resulted in a moderate enrichment of ligand screening against this challenging drug target. Results show that, using a statistical mechanics based free energy method like BEDAM starting from docked poses offers better enrichment than classical docking scoring functions and rescoring methods like Prime MM-GBSA for the Hsp90 data set in this blind challenge. Importantly, among the three methods tested here, only the mean value of the BEDAM binding free energy scores is able to separate the large group of binders from the small group of nonbinders with a gap of 2.4 kcal/mol. None of the three methods that we have tested provided accurate ranking of the affinities of the 147 active compounds. We discuss the possible sources of errors in the binding free energy calculations. The study suggests that BEDAM can be used strategically to discriminate binders from nonbinders in virtual screening and to more accurately predict the ligand binding modes prior to the more computationally expensive FEP calculations of binding affinity.