grand: A Python Module for Grand Canonical Water Sampling in OpenMM

Efficient Sampling of Cavity Hydration in Proteins with Nonequilibrium Grand Canonical Monte Carlo and Polarizable Force Fields

Prediction of the hydration levels of protein cavities and active sites is important to both mechanistic analysis and ligand design. Due to the unique microscopic environment of these buried water molecules, a polarizable model is expected to be crucial for an accurate treatment of protein internal hydration in simulations. Here we adapt a nonequilibrium candidate Monte Carlo approach for conducting grand canonical Monte Carlo simulations with the Drude polarizable force field. The GPU implementation enables the efficient sampling of internal cavity hydration levels in biomolecular systems. We also develop an enhanced sampling approach referred to as B-walking, which satisfies detailed balance and readily combines with grand canonical integration to efficiently calculate quantitative binding free energies of water to protein cavities. Applications of these developments are illustrated in a solvent box and the polar ligand binding site in trypsin. Our simulation results show that including electronic polarization leads to a modest but clear improvement in the description of water position and occupancy compared to the crystal structure. The B-walking approach enhances the range of water sampling in different chemical potential windows and thus improves the accuracy of water binding free energy calculations.

Computational models of amorphous ice for accurate simulation of cryo-EM images of biological samples

Simulations of cryo-electron microscopy (cryo-EM) images of biological samples can be used to produce test datasets to support the development of instrumentation, methods, and software, as well as to assess data acquisition and analysis strategies. To be useful, these simulations need to be based on physically realistic models which include large volumes of amorphous ice. The gold standard model for EM image simulation is a physical atom-based ice model produced using molecular dynamics simulations. Although practical for small sample volumes; for simulation of cryo-EM data from large sample volumes, this can be too computationally expensive. We have evaluated a Gaussian Random Field (GRF) ice model which is shown to be more computationally efficient for large sample volumes. The simulated EM images are compared with the gold standard atom-based ice model approach and shown to be directly comparable. Comparison with experimentally acquired data shows the Gaussian random field ice model produces realistic simulations. The software required has been implemented in the Parakeet software package and the underlying atomic models are available online for use by the wider community.

Structure and Dynamics of Confined Water Inside Diphenylalanine Peptide Nanotubes

Diphenylalanine (FF) peptides exhibit a unique ability to self-assemble into nanotubes with confined water molecules playing pivotal roles in their structure and function. This study investigates the structure and dynamics of diphenylalanine peptide nanotubes (FFPNTs) using all-atom molecular dynamics (MD) and grand canonical Monte Carlo combined with MD (GCMC/MD) simulations with both the CHARMM additive and Drude polarizable force fields. The occupancy and dynamics of confined water molecules were also examined. It was found that less than 2 confined water molecules per FF help stabilize the FFPNTs on the x-y plane. Analyses of the kinetics of confined water molecules revealed distinctive transport behaviors for bound and free water, and their respective diffusion coefficients were compared. Our results validate the importance of polarizable force field models in studying peptide nanotubes and provide insights into our understanding of nanoconfined water.

MD simulations of macromolecular crystals: Implications for the analysis of Bragg and diffuse scattering

Some of our most detailed information about structure and dynamics of macromolecules comes from X-ray-diffraction studies in crystalline environments. More than 170,000 atomic models have been deposited in the Protein Data Bank, and the number of observations (typically of intensities of Bragg diffraction peaks) is generally quite large, when compared to other experimental methods. Nevertheless, the general agreement between calculated and observed intensities is far outside the experimental precision, and the majority of scattered photons fall between the sharp Bragg peaks, and are rarely taken into account. This chapter considers how molecular dynamics simulations can be used to explore the connections between microscopic behavior in a crystalline lattice and observed scattering intensities, and point the way to new atomic models that could more faithfully recapitulate Bragg intensities and extract useful information from the diffuse scattering that lies between those peaks.

WaterKit: Thermodynamic Profiling of Protein Hydration Sites

Water desolvation is one of the key components of the free energy of binding of small molecules to their receptors. Thus, understanding the energetic balance of solvation and desolvation resulting from individual water molecules can be crucial when estimating ligand binding, especially when evaluating different molecules and poses as done in High-Throughput Virtual Screening (HTVS). Over the most recent decades, several methods were developed to tackle this problem, ranging from fast approximate methods (usually empirical functions using either discrete atom-atom pairwise interactions or continuum solvent models) to more computationally expensive and accurate ones, mostly based on Molecular Dynamics (MD) simulations, such as Grid Inhomogeneous Solvation Theory (GIST) or Double Decoupling. On one hand, MD-based methods are prohibitive to use in HTVS to estimate the role of waters on the fly for each ligand. On the other hand, fast and approximate methods show an unsatisfactory level of accuracy, with low agreement with results obtained with the more expensive methods. Here we introduce WaterKit, a new grid-based sampling method with explicit water molecules to calculate thermodynamic properties using the GIST method. Our results show that the discrete placement of water molecules is successful in reproducing the position of crystallographic waters with very high accuracy, as well as providing thermodynamic estimates with accuracy comparable to more expensive MD simulations. Unlike these methods, WaterKit can be used to analyze specific regions on the protein surface, (such as the binding site of a receptor), without having to hydrate and simulate the whole receptor structure. The results show the feasibility of a general and fast method to compute thermodynamic properties of water molecules, making it well-suited to be integrated in high-throughput pipelines such as molecular docking.

Accounting for Solvation Correlation Effects on the Thermodynamics of Water Networks in Protein Cavities

Macromolecular recognition and ligand binding are at the core of biological function and drug discovery efforts. Water molecules play a significant role in mediating the protein-ligand interaction, acting as more than just the surrounding medium by affecting the thermodynamics and thus the outcome of the binding process. As individual water contributions are impossible to measure experimentally, a range of computational methods have emerged to identify hydration sites in protein pockets and characterize their energetic contributions for drug discovery applications. Even though several methods model solvation effects explicitly, they focus on determining the stability of specific water sites independently and neglect solvation correlation effects upon replacement of clusters of water molecules, which typically happens in hit-to-lead optimization. In this work, we rigorously determine the conjoint effects of replacing all combinations of water molecules in protein binding pockets through the use of the RE-EDS multistate free-energy method, which combines Hamiltonian replica exchange (RE) and enveloping distribution sampling (EDS). Applications on the small bovine pancreatic trypsin inhibitor and four proteins of the bromodomain family illustrate the extent of solvation correlation effects on water thermodynamics, with the favorability of replacement of the water sites by pharmacophore probes highly dependent on the composition of the water network and the pocket environment. Given the ubiquity of water networks in biologically relevant protein targets, we believe our approach can be helpful for computer-aided drug discovery by providing a pocket-specific and a priori systematic consideration of solvation effects on ligand binding and selectivity.

Enhanced Grand Canonical Sampling of Occluded Water Sites Using Nonequilibrium Candidate Monte Carlo

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.

Water Networks in Complexes between Proteins and FDA-Approved Drugs

Water molecules at protein-ligand interfaces are often of significant pharmaceutical interest, owing in part to the entropy which can be released upon the displacement of an ordered water by a therapeutic compound. Protein structures may not, however, completely resolve all critical bound water molecules, or there may be no experimental data available. As such, predicting the location of water molecules in the absence of a crystal structure is important in the context of rational drug design. Grand canonical Monte Carlo (GCMC) is a computational technique that is gaining popularity for the simulation of buried water sites. In this work, we assess the ability of GCMC to accurately predict water binding locations, using a dataset that we have curated, containing 108 unique structures of complexes between proteins and Food and Drug Administration (FDA)-approved small-molecule drugs. We show that GCMC correctly predicts 81.4% of nonbulk crystallographic water sites to within 1.4 Å. However, our analysis demonstrates that the reported performance of water prediction methods is highly sensitive to the way in which the performance is measured. We also find that crystallographic water sites with more protein/ligand hydrogen bonds and stronger electron density are more reliably predicted by GCMC. An analysis of water networks revealed that more than half of the structures contain at least one ligand-contacting water network. In these cases, displacement of a water site by a ligand modification might yield unexpected results if the larger network is destabilized. Cooperative effects between waters should therefore be explicitly considered in structure-based drug design.

Open Binding Pose Metadynamics: An Effective Approach for the Ranking of Protein-Ligand Binding Poses

Predicting the correct pose of a ligand binding to a protein and its associated binding affinity is of great importance in computer-aided drug discovery. A number of approaches have been developed to these ends, ranging from the widely used fast molecular docking to the computationally expensive enhanced sampling molecular simulations. In this context, methods such as coarse-grained metadynamics and binding pose metadynamics (BPMD) use simulations with metadynamics biasing to probe the binding affinity without trying to fully converge the binding free energy landscape in order to decrease the computational cost. In BPMD, the metadynamics bias perturbs the ligand away from the initial pose. The resistance of the ligand to this bias is used to calculate a stability score. The method has been shown to be useful in reranking predicted binding poses from docking. Here, we present OpenBPMD, an open-source Python reimplementation and reinterpretation of BPMD. OpenBPMD is powered by the OpenMM simulation engine and uses a revised scoring function. The algorithm was validated by testing it on a wide range of targets and showing that it matches or exceeds the performance of the original BPMD. We also investigated the role of accurate water positioning on the performance of the algorithm and showed how the combination with a grand-canonical Monte Carlo algorithm improves the accuracy of the predictions.

An open-source molecular builder and free energy preparation workflow

Automated free energy calculations for the prediction of binding free energies of congeneric series of ligands to a protein target are growing in popularity, but building reliable initial binding poses for the ligands is challenging. Here, we introduce the open-source FEgrow workflow for building user-defined congeneric series of ligands in protein binding pockets for input to free energy calculations. For a given ligand core and receptor structure, FEgrow enumerates and optimises the bioactive conformations of the grown functional group(s), making use of hybrid machine learning/molecular mechanics potential energy functions where possible. Low energy structures are optionally scored using the gnina convolutional neural network scoring function, and output for more rigorous protein-ligand binding free energy predictions. We illustrate use of the workflow by building and scoring binding poses for ten congeneric series of ligands bound to targets from a standard, high quality dataset of protein-ligand complexes. Furthermore, we build a set of 13 inhibitors of the SARS-CoV-2 main protease from the literature, and use free energy calculations to retrospectively compute their relative binding free energies. FEgrow is freely available at https://github.com/cole-group/FEgrow, along with a tutorial.

Enhancing sampling of water rehydration upon ligand binding using variants of grand canonical Monte Carlo

Water plays an important role in mediating protein-ligand interactions. Water rearrangement upon a ligand binding or modification can be very slow and beyond typical timescales used in molecular dynamics (MD) simulations. Thus, inadequate sampling of slow water motions in MD simulations often impairs the accuracy of the accuracy of ligand binding free energy calculations. Previous studies suggest grand canonical Monte Carlo (GCMC) outperforms normal MD simulations for water sampling, thus GCMC has been applied to help improve the accuracy of ligand binding free energy calculations. However, in prior work we observed protein and/or ligand motions impaired how well GCMC performs at water rehydration, suggesting more work is needed to improve this method to handle water sampling. In this work, we applied GCMC in 21 protein-ligand systems to assess the performance of GCMC for rehydrating buried water sites. While our results show that GCMC can rapidly rehydrate all selected water sites for most systems, it fails in five systems. In most failed systems, we observe protein/ligand motions, which occur in the absence of water, combine to close water sites and block instantaneous GCMC water insertion moves. For these five failed systems, we both extended our GCMC simulations and tested a new technique named grand canonical nonequilibrium candidate Monte Carlo (GCNCMC). GCNCMC combines GCMC with the nonequilibrium candidate Monte Carlo (NCMC) sampling technique to improve the probability of a successful water insertion/deletion. Our results show that GCNCMC and extended GCMC can rehydrate all target water sites for three of the five problematic systems and GCNCMC is more efficient than GCMC in two out of the three systems. In one system, only GCNCMC can rehydrate all target water sites, while GCMC fails. Both GCNCMC and GCMC fail in one system. This work suggests this new GCNCMC method is promising for water rehydration especially when protein/ligand motions may block water insertion/removal.

py-MCMD: Python Software for Performing Hybrid Monte Carlo/Molecular Dynamics Simulations with GOMC and NAMD

py-MCMD, an open-source Python software, provides a robust workflow layer that manages communication of relevant system information between the simulation engines NAMD and GOMC and generates coherent thermodynamic properties and trajectories for analysis. To validate the workflow and highlight its capabilities, hybrid Monte Carlo/molecular dynamics (MC/MD) simulations are performed for SPC/E water in the isobaric-isothermal (NPT) and grand canonical (GC) ensembles as well as with Gibbs ensemble Monte Carlo (GEMC). The hybrid MC/MD approach shows close agreement with reference MC simulations and has a computational efficiency that is 2 to 136 times greater than traditional Monte Carlo simulations. MC/MD simulations performed for water in a graphene slit pore illustrate significant gains in sampling efficiency when the coupled-decoupled configurational-bias MC (CD-CBMC) algorithm is used compared with simulations using a single unbiased random trial position. Simulations using CD-CBMC reach equilibrium with 25 times fewer cycles than simulations using a single unbiased random trial position, with a small increase in computational cost. In a more challenging application, hybrid grand canonical Monte Carlo/molecular dynamics (GCMC/MD) simulations are used to hydrate a buried binding pocket in bovine pancreatic trypsin inhibitor. Water occupancies produced by GCMC/MD simulations are in close agreement with crystallographically identified positions, and GCMC/MD simulations have a computational efficiency that is 5 times better than MD simulations. py-MCMD is available on GitHub at https://github.com/GOMC-WSU/py-MCMD.

Comparison of Grand Canonical and Conventional Molecular Dynamics Simulation Methods for Protein-Bound Water Networks

Water molecules play important roles in all biochemical processes. Therefore, it is of key importance to obtain information of the structure, dynamics, and thermodynamics of water molecules around proteins. Numerous computational methods have been suggested with this aim. In this study, we compare the performance of conventional and grand-canonical Monte Carlo (GCMC) molecular dynamics (MD) simulations to sample the water structure, as well GCMC and grid-based inhomogeneous solvation theory (GIST) to describe the energetics of the water network. They are evaluated on two proteins: the buried ligand-binding site of a ferritin dimer and the solvent-exposed binding site of galectin-3. We show that GCMC/MD simulations significantly speed up the sampling and equilibration of water molecules in the buried binding site, thereby making the results more similar for simulations started from different states. Both GCMC/MD and conventional MD reproduce crystal-water molecules reasonably for the buried binding site. GIST analyses are normally based on restrained MD simulations. This improves the precision of the calculated energies, but the restraints also significantly affect both absolute and relative energies. Solvation free energies for individual water molecules calculated with and without restraints show a good correlation, but with large quantitative differences. Finally, we note that the solvation free energies calculated with GIST are ∼5 times larger than those estimated by GCMC owing to differences in the reference state.

Accurate determination of protein:ligand standard binding free energies from molecular dynamics simulations

Designing a reliable computational methodology to calculate protein:ligand standard binding free energies is extremely challenging. The large change in configurational enthalpy and entropy that accompanies the association of ligand and protein is notoriously difficult to capture in naive brute-force simulations. Addressing this issue, the present protocol rests upon a rigorous statistical mechanical framework for the determination of protein:ligand binding affinities together with the comprehensive Binding Free-Energy Estimator 2 (BFEE2) application software. With the knowledge of the bound state, available from experiments or docking, application of the BFEE2 protocol with a reliable force field supplies in a matter of days standard binding free energies within chemical accuracy, for a broad range of protein:ligand complexes. Limiting undesirable human intervention, BFEE2 assists the end user in preparing all the necessary input files and performing the post-treatment of the simulations towards the final estimate of the binding affinity.

Comparing the Performance of Different AMBER Protein Forcefields, Partial Charge Assignments, and Water Models for Absolute Binding Free Energy Calculations

Identifying chemical starting points is a vital first step in small molecule drug discovery and can take significant time and money. For this reason, computational approaches to virtual screening are of great interest as they can lower the cost and shorten timeframes. However, simple approaches such as molecular docking and pharmacophore screening are of limited accuracy and provide a low probability of success. Alchemical binding free energies represent a promising approach for virtual screening as they naturally incorporate the key effects of water molecules, protein flexibility, and binding entropy. However, the calculations are technically very challenging, with performance depending on the specific forcefield used. For this reason, it is important that the community has access to benchmark test sets to assess prediction accuracy. In this paper, we present an approach to alchemical binding free energies using OpenMM. We identify effective simulation parameters using an existing BRD4(1) test set and present two new benchmark sets (cMET and PDE2A) that can be used in the community for validation purposes. Our findings also highlight the effectiveness of some AMBER forcefields, in particular, AMBER ff15ipq.

Enhancing Sampling of Water Rehydration on Ligand Binding: A Comparison of Techniques

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.

How the Choice of Force-Field Affects the Stability and Self-Assembly Process of Supramolecular CTA Fibers

In recent years, computational methods have become an essential element of studies focusing on the self-assembly process. Although they provide unique insights, they face challenges, from which two are the most often mentioned in the literature: the temporal and spatial scale of the self-assembly. A less often mentioned issue, but not less important, is the choice of the force-field. The repetitive nature of the supramolecular structure results in many similar interactions. Consequently, even a small deviation in these interactions can lead to significant energy differences in the whole structure. However, studies comparing different force-fields for self-assembling systems are scarce. In this article, we compare molecular dynamics simulations for trifold hydrogen-bonded fibers performed with different force-fields, namely GROMOS, CHARMM General Force Field (CGenFF), CHARMM Drude, General Amber Force-Field (GAFF), Martini, and polarized Martini. Briefly, we tested the force-fields by simulating: (i) spontaneous self-assembly (none form a fiber within 500 ns), (ii) stability of the fiber (observed for CHARMM Drude, GAFF, MartiniP), (iii) dimerization (observed for GROMOS, GAFF, and MartiniP), and (iv) oligomerization (observed for CHARMM Drude and MartiniP). This system shows that knowledge of the force-field behavior regarding interactions in oligomer and larger self-assembled structures is crucial for designing efficient simulation protocols for self-assembling systems.

Fast Equilibration of Water between Buried Sites and the Bulk by Molecular Dynamics with Parallel Monte Carlo Water Moves on Graphical Processing Units

Molecular dynamics (MD) simulations of proteins are commonly used to sample from the Boltzmann distribution of conformational states, with wide-ranging applications spanning chemistry, biophysics, and drug discovery. However, MD can be inefficient at equilibrating water occupancy for buried cavities in proteins that are inaccessible to the surrounding solvent. Indeed, the time needed for water molecules to equilibrate between the bulk solvent and the binding site can be well beyond what is practical with standard MD, which typically ranges from hundreds of nanoseconds to a few microseconds. We recently introduced a hybrid Monte Carlo/MD (MC/MD) method, which speeds up the equilibration of water between buried cavities and the surrounding solvent, while sampling from the thermodynamically correct distribution of states. While the initial implementation of the MC functionality led to considerable slowing of the overall simulations, here we address this problem with a parallel MC algorithm implemented on graphical processing units. This results in speed-ups of 10-fold to 1000-fold over the original MC/MD algorithm, depending on the system and simulation parameters. The present method is available for use in the AMBER simulation software.

Rimantadine Binds to and Inhibits the Influenza A M2 Proton Channel without Enantiomeric Specificity

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 kon, koff, and Kd 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.

A Cloud Computing Platform for Scalable Relative and Absolute Binding Free Energy Predictions: New Opportunities and Challenges for Drug Discovery

Free energy perturbation (FEP) has become widely used in drug discovery programs for binding affinity prediction between candidate compounds and their biological targets. However, limitations of FEP applications also exist, including, but not limited to, high cost, long waiting time, limited scalability, and breadth of application scenarios. To overcome these problems, we have developed XFEP, a scalable cloud computing platform for both relative and absolute free energy predictions using optimized simulation protocols. XFEP enables large-scale FEP calculations in a more efficient, scalable, and affordable way, for example, the evaluation of 5000 compounds can be performed in 1 week using 50-100 GPUs with a computing cost roughly equivalent to the cost for the synthesis of only one new compound. By combining these capabilities with artificial intelligence techniques for goal-directed molecule generation and evaluation, new opportunities can be explored for FEP applications in the drug discovery stages of hit identification, hit-to-lead, and lead optimization based not only on structure exploitation within the given chemical series but also including evaluation and comparison of completely unrelated molecules during structure exploration in a larger chemical space. XFEP provides the basis for scalable FEP applications to become more widely used in drug discovery projects and to speed up the drug discovery process from hit identification to preclinical candidate compound nomination.