WaterKit: Thermodynamic Profiling of Protein Hydration Sites

Abstract

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.

Fig. 1.

Schematic representation of the four main steps of the WaterKit sampling protocol: 1) identification of anchor points, and affinity grid maps preparation 2) placement order (a) and sampling (b and c) methods used to generate an ensemble of discrete water molecule conformations; 3) minimization of the fully hydrated system; 4) the estimation of the thermodynamic properties; of water molecules using GIST.

Fig. 2.

Difference between spherical maps using (A) AutoDock-Vina non-directional donor-acceptor oxygen probe (O_DA) and (B) TIP3P water model for the quaternary amine group (R-NH3⁺) in lysine residues. The spherical maps were contoured at 85 % of the most favorable energy value, corresponding to −0.22 and −11.56 kcal/mol for the AutoDock-Vina O_DA probe and TIP3P model, respectively.

Fig. 3.

Results for each WaterKit model tested and averaged over all the studied systems (except PDB 1hpx for hivpr protein), varying the hydrogen bond anchor points definition (${\mathrm{A}\mathrm{P}}_{\text{polar}}$: polar HB donor/acceptors only; ${\mathrm{A}\mathrm{P}}_{\text{all}}$: polar and non-polar hydrogen atoms), the spherical water model (TIP3P water, Vina oxygen donor-acceptor) and the number of minimization steps. All the hydration sites identified in the ligand binding pocket were compared to MD simulations. The water placement with a distance cutoff of 1.0 Å was evaluated using the True Positive Rate (TPR) (sensitivity), the Positive Predictive Value (PPV) (precision). The different GIST energy components, ${\mathrm{E}}_{\text{sw}},{\mathrm{E}}_{\mathrm{w}\mathrm{w}},{\mathrm{T}\mathrm{S}}_{\text{trans}},{\mathrm{T}\mathrm{S}}_{\text{orient}}$ as well as $\mathrm{\Delta }\mathrm{G}$, were compared using the Root Mean Square Deviation (RMSD) and the coefficient of determination (r²).

Fig. 4.

(A) Results for each system studied using the best WaterKit model: polar Hydrogen Bond Anchors (HBA); positions identified using the Vina O_DA spherical water model interactions relaxed with 100 steps of minimization. All hydration sites identified in the ligand binding pocket were compared to MD simulations. The number of True Positive (TP) corresponds to the total number of water molecules correctly identified by WaterKit, using a distance cutoff of 1.0 Å. The different GIST energy components, ${\mathrm{E}}_{\mathrm{s}\mathrm{w}},{\mathrm{E}}_{\mathrm{w}\mathrm{w}},{\mathrm{T}\mathrm{S}}_{\text{trans}.},{\mathrm{T}\mathrm{S}}_{\text{orient}}.$ as well as $\mathrm{\Delta }\mathrm{G}$, were compared using the Root Mean Square Deviation (RMSD) and the coefficient of determination (${\mathrm{r}}^2$). RMSD and ${\mathrm{r}}^2$ calculated for the MD triplicates is reported as a reference. (B) The water placement was evaluated using the True Positive Rate (TPR) (sensitivity) and the Positive Predictive Value (PPV) (precision) metrics. The sensitivity and precision were calculated on all the identified hydration sites, as well on different set of hydration sites based on their ${\mathrm{E}}_{\mathrm{s}\mathrm{w}}$ energy from −5 to −2 kcal/mol.

Fig. 5.

Hydration sites found in ital ligand binding site (PDB: 2ica) using MD simulations (triplicates) and WaterKit. Hydration sites found in MD simulation are represented as blue spheres. Hydration sites found with WaterKit are colored in green (found an equivalent in MD simulations) and red (no equivalent found in MD simulations) and represented as spheres. The protein is shown as cartoon (white) with side-chains in sticks. The co-crystallized ligand is colored in orange and represented in stick.

Fig. 6.

Hydration sites found in hsp90α ligand binding site using MD simulations (triplicates) and WaterKit (PDB: 1uyg). (A) comparison of MD and WaterKit predictions; (B) WaterKit hydration sites near the purine ring; (C) MD hydration sites near the purine ring. Hydration sites found with MD are represented by blue spheres. Hydration sites found with WaterKit shown as spheres colored in green (if found also in MD simulations) and red (no equivalent found in MD simulations) and represented by spheres. Spheres are labelled by their calculated enthalpy solute-water energy (${\mathrm{E}}_{\mathrm{s}\mathrm{w}}$) in kcal/mol. The protein is shown in white with cartoons for the secondary structure and sticks for side chains and the bound ligand. Hydrogen bonds are shown as yellow dotted lines.

Fig. 7.

Hydration sites found around the crystallographic water W301 in hivpr (PDB: 2zye) using (A) WaterKit and (B) MD simulations (triplicates). Hydration sites found are represented by spheres colored in blue (found in MD simulation), green (found by WaterKit and present in MD simulations) and red (found by WaterKit but not present in MD simulations). Spheres are labelled by their calculated free energy ($\mathrm{\Delta }\mathrm{G}$) in kcal/mol. The protein is colored in white and represented in cartoon secondary structure and side-chains in sticks. The co-crystallized inhibitor KNI-272 is shown as orange sticks. Hydrogen bonds are represented by yellow dotted lines. Oxygen density (gO) map from WaterKit and MD simulation averages, contoured at 5.0 bulk density, is colored in green cyan and deep purple, respectively. (C) Network of hydration sites predicted by WaterKit within 1.5 Å from the heavy atoms of ligand KNI-272 bound in HIV-1 PR. Hydration sites found by Waterkit are represented by spheres and colored in green (if also found in MD simulations) and red (no equivalent found in MD simulations). KNI-272 ligand is shown as gray sticks. Hydration sites overlapping with key hydrogen bond features are surrounded by an additional transparent sphere.

Fig. 8.

Density expressed in bulk density unit in the ligand binding pocket of fabp4 (A) during MD simulations (triplicates) of 200 ns long each and (B) in the ensemble of 10,000 frames generated with WaterKit. The average density in WaterKit and MD simulations are represented by dotted red lines. The histograms represent the observed probability distribution of the density.

Fig. 9.

Comparison between (A) conserved crystallographic water molecules found in apo structures and (B-F) hydration sites predicted with WaterKit in the neuraminidase active site. Crystallographic water molecules are shown in red; all water molecules predicted by WaterKit are in gray, and water molecules are shown and WaterKit hydration sites are shown as red and green spheres, respectively. (A) Sialic acid bound to neuraminidase (PDB id:2bat ) with crystallographic water molecules from apo neuraminidase structures (PDB entries 3nn9, 4nn9, 5nn9, 6nn9, 6crd, 6d3b, 6mcx) within 1.5 of ligand heavy atoms shown in in B-F. (B) Sialic acid (PDB id: 2bat) superposed to WaterKit results. (C) Oseltamivir (PDB id: 2ht7, ). (D) Zanamivir (PDB id: 3ckz ). (E) FeqGuDFSA (PDB id: 3w09 ). (F) BANA206 (PDB id: 1b9v )

Fig. 10.

Hydration sites found (A) forming a five-membered ring-like structure in the streptavidin ligand binding site and (B) overlapping with biotin in the ligand binding site using MD simulations (triplicates) and WaterKit (WK). Hydration sites found in MD simulation are represented by blue spheres. Hydration sites found with WaterKit are colored in green (found an equivalent in MD simulations) and red (no equivalent found in MD simulations) and represented by spheres. The spheres are labelled by their enthalpy solute-water energy (${\mathrm{E}}_{\mathrm{s}\mathrm{w}}$) in kcal/mol. Crystallographic water molecules found in apo structures of the streptavidin (pdb ids: 7knk, 7ek8 and 7ek9) are represented by crosses colored in red. The protein is colored in white and represented in cartoon with side-chains in sticks. The oxygen density (gO) map from WaterKit and MD simulations, contoured at 5.0 bulk density, is colored in green cyan and deep purple, respectively. The gO map for MD simulations was obtained by average the gO map from each MD replicates. Hydrogen bonds between hydration sites in MD are represented by yellow dotted lines.