OpenBioSim is pleased to announce the initial public release of loch, our new CUDA accelerated Grand Canonical Monte Carlo (GCMC) water sampling engine, built on top of Sire, BioSimSpace, and OpenMM.
Water molecules in protein-ligand interactions
In drug design it is well appreciated that desolvation of protein binding sites and water-mediated protein-ligand interactions play a crucial role in determining the binding energetics of a ligand to a protein. Atomistic, explicit solvent molecular dynamics (MD) simulations are widely used to tackle this problem. However because such MD simulations are typically carried out on relatively short (sub microsecond) timescales it can be challenging to equilibrate the water content of certain protein binding sites. Further, kinetic trapping of water molecules in unrealistically high energy states is a concern for practitioners of FEP calculations.
Grand-Canonical Monte Carlo to the Rescue
The field of molecular simulation has developed numerous methodologies to address this problem. One of the most successful approach leverages the Grand Canonical Monte Carlo (GCMC) ensemble to enable insertion and deletion of water molecules over the course of a molecular dynamics simulation. Foundational GCMC work was done by Mezei in 1987. Starting from the late 90s significant extensions of the framework to handle the treatment of water molecules in biomolecular systems were pursued by several research groups, notably the group of Professor Jonathan Essex at the University of Southampton (see this review article for a summary of the body of work). Nowadays several implementations of GCMC water sampling are available in different MD packages (see references below for the implementations that inspired the design of loch).
While the GCMC methodology is robust, deployment in production settings requires careful engineering for efficient integration with Molecular Dynamics packages. The chief issue is that modern MD simulation packages leverage efficiently GPUs to accelerate the computation of a timestep in a MD trajectory. Naive deployment of Monte Carlo algorithms on GPUs doesn’t scale as well because of the number of operations per cycle in a Monte Carlo step is a lot smaller than in a Molecular dynamics step. Other issues include overheads associated with pausing, updating and restarting a MD simulation when interleaving chunks of MD sampling time with MC cycles.
Loch is an efficient GCMC solution to sample the water molecules that matter
To enable parallelisation of insertion and deletion attepts, loch leverages GPU capabilities using a custom CUDA kernel for nonbonded interactions. This allows thousands of GCMC trials to be attempted in parallel, significantly enhancing sampling efficiency compared to traditional CPU-based implementations that perform sequential attempts via the OpenMM Python API. Additionally, electrostatics for GCMC attempts are computed using the reaction field (RF) method, with accepted candidates being re-evaluated with a correction step based on the difference between the reaction field and Particle Mesh Ewald (PME) potential energies. The use of an approximate potential for trials moves leads to a substantial speed-up in GCMC move evaluation. loch has been designed to be modular, allowing standalone GCMC sampling, or integration with OpenMM-based molecular dynamics simulation code.
GCMC bulk density for TIP3P water at a temperature of 298K and pressure of 1bar compared to constant pressure (NPT) simulations using both RF and PME electrostatics. In both cases, the agreement with NPT is excellent.
GCMC water clusters (overlaid on the crystal structure waters) obtained from simulation of bovine pancreatic trypsin inhibitor (BPTI). Blue spheres indicate the positions of the crystallographic waters, while red spheres indicate the position of the water oxygen atom from simulation that were closest to the centre of each cluster. Darker red spheres indicate clusters with higher occupancy. The results show good agreement with the crystallographic waters, as well as identifying several additional lower occupancy water sites.
loch has been integrated in the somd2 package (that is still under development), enabling protein-ligand FEP MD/GCMC cycles with modest overheads.
loch is based on grand, M. L. Samways, et. al., J. Chem. Inf. Model., 2020, 60, 4436-4441 (DOI), O. J. Melling et. al., J. Chem. Theory Comput., 2023 (DOI).The parallelisation strategy comes from Ross et. al, (JCTC 2020), and we use the approximate potential method introduced by Gelb (J. Chem. Phys., 2003).
Written by Lester Hedges.