Networks of water molecules can play a critical role at the protein-ligand interface, and can directly influence drug-target interactions. Grand canonical methods aid in the sampling of these water molecules, where conventional molecular dynamics equilibration times are often long, by allowing waters to be inserted and deleted from the system, according to the chemical potential. Here, we present our open source Python module, grand (https://github.com/essex-lab/grand), which allows molecular dynamics simulations to be performed in conjunction with grand canonical Monte Carlo sampling, using the OpenMM simulation engine. We demonstrate the accuracy of this module by reproducing the density of bulk water observed from constant pressure simulations.Application of this code to the bovine pancreatic trypsin inhibitor protein reproduces three buried crystallographic water sites that are poorly sampled using conventional molecular dynamics.

In this review we examine computational approaches to explore the structure and thermodynamics of water binding in protein–drug complexes

Water often plays a key role in protein structure, molecular recognition, and mediating protein−ligand interactions. Thus, free energy calculations must adequately sample water motions, which often proves challenging in typical MD simulation time scales. Thus, the accuracy of methods relying on MD simulations ends up limited by slow water sampling. Particularly, as a ligand is removed or modified, bulk water may not have time to fill or rearrange in the binding site. In this work, we focus on several molecular dynamics (MD) simulation-based methods attempting to help rehydrate buried water sites: BLUES, using nonequilibrium candidate Monte Carlo (NCMC); grand, using grand canonical Monte Carlo (GCMC); and normal MD. We assess the accuracy and efficiency of these methods in rehydrating target water sites. We selected a range of systems with varying numbers of waters in the binding site, as well as those where water occupancy is coupled to the identity or binding mode of the ligand. We analyzed the rehydration of buried water sites in binding pockets using both clustering of trajectories and direct analysis of electron density maps. Our results suggest both BLUES and grand enhance water sampling relative to normal MD and grand is more robust than BLUES, but also that water sampling remains a major challenge for all of the methods tested. The lessons we learned for these methods and systems are discussed.

Water molecules play a key role in many biomolecular systems, particularly when bound at protein–ligand interfaces. However, molecular simulation studies on such systems are hampered by the relatively long time scales over which water exchange between a protein and solvent takes place. Grand canonical Monte Carlo (GCMC) is a simulation technique that avoids this issue by attempting the insertion and deletion of water molecules within a given structure. The approach is constrained by low acceptance probabilities for insertions in congested systems, however. To address this issue, here, we combine GCMC with nonequilibium candidate Monte Carlo (NCMC) to yield a method that we refer to as grand canonical nonequilibrium candidate Monte Carlo (GCNCMC), in which the water insertions and deletions are carried out in a gradual, nonequilibrium fashion. We validate this new approach by comparing GCNCMC and GCMC simulations of bulk water and three protein binding sites. We find that not only is the efficiency of the water sampling improved by GCNCMC but that it also results in increased sampling of ligand conformations in a protein binding site, revealing new water-mediated ligand-binding geometries that are not observed using alternative enhanced sampling techniques.

The influenza A M2 wild-type (WT) proton channel is the target of the anti-influenza drug rimantadine. Rimantadine has two enantiomers, though most investigations into drug binding and inhibition have used a racemic mixture. Solid-state NMR experiments using the full length-M2 WT have shown significant spectral differences that were interpreted to indicate tighter binding for (R)- vs (S)-rimantadine. However, it was unclear if this correlates with a functional difference in drug binding and inhibition. Using X-ray crystallography, we have determined that both (R)- and (S)-rimantadine bind to the M2 WT pore with slight differences in the hydration of each enantiomer. However, this does not result in a difference in potency or binding kinetics, as shown by similar values for k on, k off, and K d in electrophysiological assays and for EC50 values in cellular assays. We concluded that the slight differences in hydration for the (R)- and (S)-rimantadine enantiomers are not relevant to drug binding or channel inhibition. To further explore the effect of the hydration of the M2 pore on binding affinity, the water structure was evaluated by grand canonical ensemble molecular dynamics simulations as a function of the chemical potential of the water. Initially, the two layers of ordered water molecules between the bound drug and the channel’s gating His37 residues mask the drug’s chirality. As the chemical potential becomes more unfavorable, the drug translocates down to the lower water layer, and the interaction becomes more sensitive to chirality. These studies suggest the feasibility of displacing the upper water layer and specifically recognizing the lower water layers in novel drugs.