Free Energy Density of a Fluid and Its Role in Solvation and Binding.

The concept that a fluid has a position-dependent free energy density appears in the literature but has not been fully developed or accepted. We set this concept on an unambiguous theoretical footing via the following strategy. First, we set forth four desiderata that should be satisfied by any definition of the position-dependent free energy density, f(R), in a system comprising only a fluid and a rigid solute: its volume integral, plus the fixed internal energy of the solute, should be the system free energy; it deviates from its bulk value, fbulk, near a solute but should asymptotically approach fbulk with increasing distance from the solute; it should go to zero where the solvent density goes to zero; and it should be well-defined in the most general case of a fluid made up of flexible molecules with an arbitrary interaction potential. Second, we use statistical thermodynamics to formulate a definition of the free energy density that satisfies these desiderata. Third, we show how any free energy density satisfying the desiderata may be used to analyze molecular processes in solution. In particular, because the spatial integral of f(R) equals the free energy of the system, it can be used to compute free energy changes that result from the rearrangement of solutes as well as the forces exerted on the solutes by the solvent. This enables the use of a thermodynamic analysis of water in protein binding sites to inform ligand design. Finally, we discuss related literature and address published concerns regarding the thermodynamic plausibility of a position-dependent free energy density. The theory presented here has applications in theoretical and computational chemistry and may be further generalizable beyond fluids, such as to solids and macromolecules.

AmberTools.

AmberTools is a free and open-source collection of programs used to set up, run, and analyze molecular simulations. The newer features contained within AmberTools23 are briefly described in this Application note.

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.

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.

Thermodynamic Decomposition of Solvation Free Energies with Particle Mesh Ewald and Long-Range Lennard-Jones Interactions in Grid Inhomogeneous Solvation Theory.

Grid Inhomogeneous Solvation Theory (GIST) maps out solvation thermodynamic properties on a fine meshed grid and provides a statistical mechanical formalism for thermodynamic end-state calculations. However, differences in how long-range nonbonded interactions are calculated in molecular dynamics engines and in the current implementation of GIST have prevented precise comparisons between free energies estimated using GIST and those from other free energy methods such as thermodynamic integration (TI). Here, we address this by presenting PME-GIST, a formalism by which particle mesh Ewald (PME)-based electrostatic energies and long-range Lennard-Jones (LJ) energies are decomposed and assigned to individual atoms and the corresponding voxels they occupy in a manner consistent with the GIST approach. PME-GIST yields potential energy calculations that are precisely consistent with modern simulation engines and performs these calculations at a dramatically faster speed than prior implementations. Here, we apply PME-GIST end-state analyses to 32 small molecules whose solvation free energies are close to evenly distributed from 2 kcal/mol to -17 kcal/mol and obtain solvation energies consistent with TI calculations (R2 = 0.99, mean unsigned difference 0.8 kcal/mol). We also estimate the entropy contribution from the second and higher order entropy terms that are truncated in GIST by the differences between entropies calculated in TI and GIST. With a simple correction for the high order entropy terms, PME-GIST obtains solvation free energies that are highly consistent with TI calculations (R2 = 0.99, mean unsigned difference = 0.4 kcal/mol) and experimental results (R2 = 0.88, mean unsigned difference = 1.4 kcal/mol). The precision of PME-GIST also enables us to show that the solvation free energy of small hydrophobic and hydrophilic molecules can be largely understood based on perturbations of the solvent in a region extending a few solvation shells from the solute. We have integrated PME-GIST into the open-source molecular dynamics analysis software CPPTRAJ.

An online repository of solvation thermodynamic and structural maps of SARS-CoV-2 targets.

SARS-CoV-2 recently jumped species and rapidly spread via human-to-human transmission to cause a global outbreak of COVID-19. The lack of effective vaccine combined with the severity of the disease necessitates attempts to develop small molecule drugs to combat the virus. COVID19_GIST_HSA is a freely available online repository to provide solvation thermodynamic maps of COVID-19-related protein small molecule drug targets. Grid inhomogeneous solvation theory maps were generated using AmberTools cpptraj-GIST, 3D reference interaction site model maps were created with AmberTools rism3d.snglpnt and hydration site analysis maps were created using SSTMap code. The resultant data can be applied to drug design efforts: scoring solvent displacement for docking, rational lead modification, prioritization of ligand- and protein- based pharmacophore elements, and creation of water-based pharmacophores. Herein, we demonstrate the use of the solvation thermodynamic mapping data. It is hoped that this freely provided data will aid in small molecule drug discovery efforts to defeat SARS-CoV-2.

An online repository of solvation thermodynamic and structural maps of SARS-CoV-2 targets.

SARS-CoV-2 recently jumped species and rapidly spread via human-to-human transmission to cause a global outbreak of COVID-19. The lack of effective vaccine combined with the severity of the disease necessitates attempts to develop small molecule drugs to combat the virus. COVID19_GIST_HSA is a freely available online repository to provide solvation thermodynamic maps of COVID-19-related protein small molecule drug targets. Grid Inhomogeneous Solvation Theory maps were generated using AmberTools cpptraj-GIST and Hydration Site Analysis maps were created using SSTmap code. The resultant data can be applied to drug design efforts: scoring solvent displacement for docking, rational lead modification, prioritization of ligand- and protein- based pharmacophore elements, and creation of water-based pharmacophores. Herein, we demonstrate the use of the solvation thermodynamic mapping data. It is hoped that this freely provided data will aid in small molecule drug discovery efforts to defeat SARS-CoV-2.

Inclusion of enclosed hydration effects in the binding free energy estimation of dopamine D3 receptor complexes.

Confined hydration and conformational flexibility are some of the challenges encountered for the rational design of selective antagonists of G-protein coupled receptors. We present a set of C3-substituted (-)-stepholidine derivatives as potent binders of the dopamine D3 receptor. The compounds are characterized biochemically, as well as by computer modeling using a novel molecular dynamics-based alchemical binding free energy approach which incorporates the effect of the displacement of enclosed water molecules from the binding site. The free energy of displacement of specific hydration sites is obtained using the Hydration Site Analysis method with explicit solvation. This work underscores the critical role of confined hydration and conformational reorganization in the molecular recognition mechanism of dopamine receptors and illustrates the potential of binding free energy models to represent these key phenomena.

Role of Displacing Confined Solvent in the Conformational Equilibrium of ß-Cyclodextrin.

This study investigates the role of hydration and its relationship to the conformational equilibrium of the host molecule ß-cyclodextrin. Molecular dynamics simulations indicate that the unbound ß-cyclodextrin exhibits two state behavior in explicit solvent due to the opening and closing of its cavity. In implicit solvent, these transitions are not observed, and there is one dominant conformation of ß-cyclodextrin with an open cavity. Based on these observations, we investigate the hypothesis that the expulsion of thermodynamically unfavorable water molecules into the bulk plays an important role in controlling the accessibility of the closed macrostate at room temperature. We compare the results of the molecular mechanics analytical generalized Born plus nonpolar solvation approach to those obtained through grid inhomogeneous solvation theory analysis with explicit solvation to elucidate the thermodynamic forces at play. The work illustrates the use of continuum solvent models to tease out solvation effects related to the inhomogeneity and the molecular nature of water and demonstrates the key role of the thermodynamics of enclosed hydration in driving the conformational equilibrium of molecules in solution.

Hidden bias in the DUD-E dataset leads to misleading performance of deep learning in structure-based virtual screening.

Recently much effort has been invested in using convolutional neural network (CNN) models trained on 3D structural images of protein-ligand complexes to distinguish binding from non-binding ligands for virtual screening. However, the dearth of reliable protein-ligand x-ray structures and binding affinity data has required the use of constructed datasets for the training and evaluation of CNN molecular recognition models. Here, we outline various sources of bias in one such widely-used dataset, the Directory of Useful Decoys: Enhanced (DUD-E). We have constructed and performed tests to investigate whether CNN models developed using DUD-E are properly learning the underlying physics of molecular recognition, as intended, or are instead learning biases inherent in the dataset itself. We find that superior enrichment efficiency in CNN models can be attributed to the analogue and decoy bias hidden in the DUD-E dataset rather than successful generalization of the pattern of protein-ligand interactions. Comparing additional deep learning models trained on PDBbind datasets, we found that their enrichment performances using DUD-E are not superior to the performance of the docking program AutoDock Vina. Together, these results suggest that biases that could be present in constructed datasets should be thoroughly evaluated before applying them to machine learning based methodology development.

A molecular reconstruction approach to site-based 3D-RISM and comparison to GIST hydration thermodynamic maps in an enzyme active site.

Computed, high-resolution, spatial distributions of solvation energy and entropy can provide detailed information about the role of water in molecular recognition. While grid inhomogeneous solvation theory (GIST) provides rigorous, detailed thermodynamic information from explicit solvent molecular dynamics simulations, recent developments in the 3D reference interaction site model (3D-RISM) theory allow many of the same quantities to be calculated in a fraction of the time. However, 3D-RISM produces atomic-site, rather than molecular, density distributions, which are difficult to extract physical meaning from. To overcome this difficulty, we introduce a method to reconstruct molecular density distributions from atomic-site density distributions. Furthermore, we assess the quality of the resulting solvation thermodynamics density distributions by analyzing the binding site of coagulation Factor Xa with both GIST and 3D-RISM. We find good qualitative agreement between the methods for oxygen and hydrogen densities as well as direct solute-solvent energetic interactions. However, 3D-RISM predicts lower energetic and entropic penalties for moving water from the bulk to the binding site.

Simulating Water Exchange to Buried Binding Sites.

Traditional molecular dynamics (MD) simulations of proteins, which relies on integration of Newton's equations of motion, cannot efficiently equilibrate water occupancy for buried cavities in proteins. This leads to slow convergence of thermodynamic averages for such systems. We have addressed this challenge by efficiently integrating standard Metropolis Monte Carlo (MC) translational water moves with MD in the AMBER simulation package. The translational moves allow water to easily enter or exit buried sites in a thermodynamically correct way during a simulation. To maximize efficiency, the algorithm avoids moves that only interchange waters within the bulk around the protein instead focusing on moves that can transfer water between bulk and the protein interior. In addition, a steric grid allows avoidance of moves that would lead to obvious steric clashes, and a fast grid-based energy evaluation is used to reduce the number of expensive full energy calculations. The potential energy distribution produced using MC/MD was found to be statistically indistinguishable from that of control simulations using only MD, and the algorithm effectively equilibrated water across steric barriers and into binding pockets that are not accessible with pure MD. The MC/MD method introduced here should be of increasing utility for applications spanning protein folding, the elucidation of protein mechanisms, and free energy calculations for computer-aided drug design. It is available in version 18 release of the widely disseminated AMBER simulation package.

New Dopamine D3-Selective Receptor Ligands Containing a 6-Methoxy-1,2,3,4-tetrahydroisoquinolin-7-ol Motif.

A series of analogues featuring a 6-methoxy-1,2,3,4-tetrahydroisoquinolin-7-ol unit as the arylamine "head" group of a classical D3 antagonist core structure were synthesized and evaluated for affinity at dopamine D1, D2, and D3 receptors (D1R, D2R, D3R). The compounds generally displayed strong affinity for D3R with very good D3R selectivity. Docking studies at D2R and D3R crystal structures revealed that the molecules are oriented such that their arylamine units are positioned in the orthosteric binding pocket of D3R, with the arylamide "tail" units residing in the secondary binding pocket. Hydrogen bonding between Ser 182 and Tyr 365 at D3R stabilize extracellular loop 2 (ECL2), which in turn contributes to ligand binding by interacting with the "tail" units of the ligands in the secondary binding pocket. Similar interactions between ECL2 and the "tail" units were absent at D2R due to different positioning of the D2R loop region. The presence of multiple H-bonds with the phenol moiety of the headgroup of 7 and Ser192 accounts for its stronger D3R affinity as compared to the 6,7-dimethoxy-1,2,3,4-tetrahydroisoquinoline-containing analogue 8.

Water Pharmacophore: Designing Ligands using Molecular Dynamics Simulations with Water.

In this study, we demonstrate a method to construct a water-based pharmacophore model which can be utilized in the absence of known ligands. This method utilizes waters found in the binding pocket, sampled through molecular dynamics. Screening of compound databases against this water-based pharmacophore model reveals that this approach can successfully identify known binders to a target protein. The method was tested by enrichment studies of 7 therapeutically important targets and compared favourably to screening-by-docking with Glide. Our results suggest that even without experimentally known binders, pharmacophore models can be generated using molecular dynamics with waters and used for virtual screening.

Solvation Structure and Thermodynamic Mapping (SSTMap): An Open-Source, Flexible Package for the Analysis of Water in Molecular Dynamics Trajectories.

We have developed SSTMap, a software package for mapping structural and thermodynamic water properties in molecular dynamics trajectories. The package introduces automated analysis and mapping of local measures of frustration and enhancement of water structure. The thermodynamic calculations are based on Inhomogeneous Fluid Solvation Theory (IST), which is implemented using both site-based and grid-based approaches. The package also extends the applicability of solvation analysis calculations to multiple molecular dynamics (MD) simulation programs by using existing cross-platform tools for parsing MD parameter and trajectory files. SSTMap is implemented in Python and contains both command-line tools and a Python module to facilitate flexibility in setting up calculations and for automated generation of large data sets involving analysis of multiple solutes. Output is generated in formats compatible with popular Python data science packages. This tool will be used by the molecular modeling community for computational analysis of water in problems of biophysical interest such as ligand binding and protein function.

Testing inhomogeneous solvation theory in structure-based ligand discovery.

Binding-site water is often displaced upon ligand recognition, but is commonly neglected in structure-based ligand discovery. Inhomogeneous solvation theory (IST) has become popular for treating this effect, but it has not been tested in controlled experiments at atomic resolution. To do so, we turned to a grid-based version of this method, GIST, readily implemented in molecular docking. Whereas the term only improves docking modestly in retrospective ligand enrichment, it could be added without disrupting performance. We thus turned to prospective docking of large libraries to investigate GIST's impact on ligand discovery, geometry, and water structure in a model cavity site well-suited to exploring these terms. Although top-ranked docked molecules with and without the GIST term often overlapped, many ligands were meaningfully prioritized or deprioritized; some of these were selected for testing. Experimentally, 13/14 molecules prioritized by GIST did bind, whereas none of the molecules that it deprioritized were observed to bind. Nine crystal complexes were determined. In six, the ligand geometry corresponded to that predicted by GIST, for one of these the pose without the GIST term was wrong, and three crystallographic poses differed from both predictions. Notably, in one structure, an ordered water molecule with a high GIST displacement penalty was observed to stay in place. Inclusion of this water-displacement term can substantially improve the hit rates and ligand geometries from docking screens, although the magnitude of its effects can be small and its impact in drug binding sites merits further controlled studies.

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.

Solvation thermodynamic mapping of molecular surfaces in AmberTools: GIST.

The expulsion of water from surfaces upon molecular recognition and nonspecific association makes a major contribution to the free energy changes of these processes. In order to facilitate the characterization of water structure and thermodynamics on surfaces, we have incorporated Grid Inhomogeneous Solvation Theory (GIST) into the CPPTRAJ toolset of AmberTools. GIST is a grid-based implementation of Inhomogeneous Fluid Solvation Theory, which analyzes the output from molecular dynamics simulations to map out solvation thermodynamic and structural properties on a high-resolution, three-dimensional grid. The CPPTRAJ implementation, called GIST-cpptraj, has a simple, easy-to-use command line interface, and is open source and freely distributed. We have also developed a set of open-source tools, called GISTPP, which facilitate the analysis of GIST output grids. Tutorials for both GIST-cpptraj and GISTPP can be found at ambermd.org. © 2016 Wiley Periodicals, Inc.

Enthalpic Breakdown of Water Structure on Protein Active-Site Surfaces.

The principles underlying water reorganization around simple nonpolar solutes are well understood and provide the framework for the classical hydrophobic effect, whereby water molecules structure themselves around solutes so that they maintain favorable energetic contacts with both the solute and the other water molecules. However, for certain solute surface topographies, water molecules, due to their geometry and size, are unable to simultaneously maintain favorable energetic contacts with both the surface and neighboring water molecules. In this study, we analyze the solvation of ligand-binding sites for six structurally diverse proteins using hydration site analysis and measures of local water structure, in order to identify surfaces at which water molecules are unable to structure themselves in a way that maintains favorable enthalpy relative to bulk water. These surfaces are characterized by a high degree of enclosure, weak solute-water interactions, and surface constraints that induce unfavorable pair interactions between neighboring water molecules. Additionally, we find that the solvation of charged side chains in an active site generally results in favorable enthalpy but can also lead to pair interactions between neighboring water molecules that are significantly unfavorable relative to bulk water. We find that frustrated local structure can occur not only in apolar and weakly polar pockets, where overall enthalpy tends to be unfavorable, but also in charged pockets, where overall water enthalpy tends to be favorable. The characterization of local water structure in these terms may prove useful for evaluating the displacement of water from diverse protein active-site environments.