Edit · 2026 01 10 090741 Nature S41592 025 02724 0

Saving needs a token. Log in or register to get one. Edits are public and licensed CC-BY-SA-4.0.

# An experimental proxy of water displaceability for ligand discovery | Nature Methods

**URL:** https://www.nature.com/articles/s41592-025-02724-0

**Saved:** 2026-01-10T17:07:46.857Z

---

Article
Published: 27 June 2025
An experimental proxy of water displaceability for ligand discovery
Justin T. Seffernick & Marcus Fischer

Nature Methods 
volume
 22, 
pages
1476–1485 (2025)Cite this article

2660 Accesses

1 Citations

63 Altmetric

Metrics
details

Abstract

Understanding protein hydration at the atomic level is a hallmark challenge in deciphering ligand binding and guiding ligand discovery. The fact that 95% of deposited crystal structures are cryogenic currently limits the utility of crystallographic water features, because water networks change with temperature. Here, to increase our confidence in interpreting cryogenic water networks, we developed ColdBrew—a machine learning method that predicts the likelihood of cryogenic crystallographic water molecules appearing at room temperature. To link ColdBrew probabilities to water displaceability, we examined >1 million waters in ligand-bound cryogenic structures and found that ColdBrew sufficiently separated conserved waters from displaceable or absent ones. To link our metric to water energies, we applied inhomogeneous solvation theory to show that waters with high ColdBrew probability and low displaceability indeed have more favorable energies. Unlike expensive computational methods, ColdBrew provides a scalable approach to assessing the utility of water molecules for ligand discovery directly from experimental cryogenic crystal structures. We precalculated ColdBrew probabilities for >46 million waters in over 100,000 Protein Data Bank structures.

You have full access to this article via your institution.

Download PDF
Similar content being viewed by others
Complex water networks visualized by cryogenic electron microscopy of RNA
Article Open access
11 March 2025
Dehydration of a crystal hydrate at subglacial temperatures
Article Open access
12 April 2023
Implicit water model within the Zimm-Bragg approach to analyze experimental data for heat and cold denaturation of proteins
Article Open access
04 May 2021
Main

Water structure plays an important role in biological processes from catalysis1,2, folding and transport3 to stabilizing protein conformational states4 linked to ligand affinity5 and selectivity6. To resolve ordered water networks experimentally, ligand discoverers typically consult high-resolution datasets derived from X-ray crystallography using synchrotron radiation7 or X-ray free-electron lasers8 and, more recently, from cryo-electron microscopy9,10. The hope is that understanding water networks with atomic accuracy will instruct structure-based drug designers on which waters to spare, displace or engage during ligand discovery. However, the reality is that our confidence in water positions from X-ray crystallographic structures is thwarted by experimental caveats of resolving and modeling water molecules. One such caveat is the data collection temperature. Crystallographic datasets in the Protein Data Bank (PDB) are almost exclusively (~95%) collected at cryogenic (cryo) temperatures of ~100 K (ref. 11). We and others12,13 have shown that this common practice distorts protein conformational landscapes, including changes in water networks. For instance, waters may differ between temperature-paired, otherwise identical datasets and stabilize different loop and sidechain conformations2,14,15.

Looking beyond the model into the electron density maps may reveal that well-defined cryo waters are missing at room temperature (RT). For instance, a pair of structures of a ligand-bound heat shock protein 90 α (Hsp90α) N‐terminal domain has several cryo waters with clear electron density (PDB: 7S9F) that were not observed at RT (PDB: 7S9G), despite higher resolution at RT15 (Fig. 1a). At RT, water molecules within the crystal lattice have higher mobility and energy compared with cryo, which, in addition to the higher vapor pressure, make them more likely to escape from a RT crystal structure, leading to desolvation. For these reasons, cryo crystallography stabilizes additional waters, anecdotally about twice as many as found in RT structures2,16. In the absence of temperature-paired datasets, such differences go unnoticed and may mislead drug design efforts. Furthermore, even when analyses are conducted cautiously and systematically, pinpointing the contributions of individual waters remains difficult because of collateral water network effects5. Given the complexity of water dynamics and the uncertainties of modeling water networks, individual waters are often ignored. Among the options to account for water are more17,18,19 or less20,21,22,23,24 expensive calculations to complement or even serve as surrogates for experimental datasets. Although these methods have had notable successes demonstrating the benefit of including water in ligand discovery pipelines25, such calculations demand expertise, infrastructure and time. Yet, the resulting energies do not necessarily reflect the ground truth.

Fig. 1: ColdBrew predicts RT water networks from cryo X-ray structures.

a, An example of water discrepancies in a cryo–RT pair for BIIB021 bound to Hsp90α (PDB: 3QDD for cryo, 7S9G for RT) with electron density. Predictions for cryo waters in holo (3QDD, 7S9G) and apo (6GPW, 7S8Z) Hsp90α structures using the model shown in c. b, Plots for five quantitative water metrics show different distributions for the categories ‘present at RT’ and ‘not present at RT’. c, ROC curves (true positive rate (TPR) versus false positive rate (FPR)) for individual metrics predicting the category in order from best to worst. The ROC curves are similar for the random forest model (right) for the training set, test set and overall set and are best for the subset of waters within 7 Å of the ligand. d, Five examples of ColdBrew predictions for SiaP (2XA5, 6H75), GFP (6QQ8, 6QQA), trypsin (7BSA, 7WB9), CypD (3RDC, 4ZSC) and FABP3 (7WD6, 7FCX). The key to the color scheme in a, b and d is shown in the center box. e, Box plots showing the distribution of the distances by which waters present at both cryo and RT shifted upon temperature change across ColdBrew probability bins. Box plots were made using the following parameters: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; outliers not shown. For the nine bins, n = 944, 2,622, 3,398, 3,427, 4,042, 4,156, 4,240, 4,750 and 5,322, respectively. f, An example of ColdBrew probability use in ligand design: a KDN-bound low-probability water in the SiaP binding site is displaced by ligands Neu5Ac and Neu5Gc.

Full size image

By introducing ColdBrew, we aim to increase confidence in observing cryo waters at elevated temperatures at which biological processes and binding experiments take place and to garnish ligand discovery pipelines with information on the displaceability of individual water molecules. In this study, we demonstrate that water conservation across temperatures is related to five easy-to-calculate metrics from the cryo map and model: normalized B factor (Bnorm), real-space correlation coefficient (RSCC), electron density of individual atoms (EDIA)26,27, solvent-accessible surface area (SASA) and the number of hydrogen bonds to protein (HB). Some of these metrics have been previously suggested to relate to conservation in general28,29,30, although not as a function of temperature. Combining these metrics into a random forest model achieved good predictive accuracy of whether a cryo water is present at RT. We then link this predictive value to water displaceability by analyzing a diverse set of ligand-bound cryo crystal structures. We show that waters in cryo crystal structures that are probably present at RT are also more likely to be conserved as ligands change. Finally, we demonstrate the link between our RT predictive value and water energy from grid inhomogeneous solvation theory (GIST)17,31. Our metric separates conserved, absent and displaceable water distributions as effectively as expensive energy calculations. Therefore, ColdBrew provides a facile experimental proxy that leverages information on water networks from cryo crystal structures to facilitate pragmatic decisions regarding water displacement in ligand discovery.

Results
Machine learning predicts water likelihood to appear at RT

As approximately 95% of all crystal structures in the PDB are collected under cryo conditions11, water analysis typically starts with a cryo structure. Because water networks differ with temperature, we wondered whether we could predict RT water networks from cryo structures. To predict water conservation across temperatures, we curated a dataset of 242 matched cryo–RT PDB crystal structure pairs covering a diverse range of protein classes (Supplementary Fig. 1). Within the dataset, the median resolution was 1.73 Å at RT and 1.70 Å at cryo with a median resolution difference for pairs of 0.15 Å and a median pairwise backbone Cα root mean square deviation (r.m.s.d.) of 0.31 Å. In total, there were 0.97 and 0.59 waters per amino acid for cryo and RT structures, respectively. We first classified whether cryo waters were also present in the paired RT structure using clustering, defining these waters as ‘present at RT’ or ‘not present at RT’. Of the 71,145 cryo waters, fewer than half (32,901) were also present at RT. Figure 1a shows an example of a temperature pair for which only 10 of the 13 ligand-proximal cryo water molecules appeared at RT.

To link water conservation to temperature, we used five metrics, Bnorm, RSCC, EDIA26,27, SASA and HB, some of which had previously been used for temperature-independent assessment of water conservation28,29,30. Figure 1b shows the distribution of these metrics for the two categories (present at RT or not present at RT), each of which followed the expected trends: lower Bnorm values suggest less water mobility; higher RSCC reports on agreement between the map and modeled water (range 0–1); higher EDIA quantifies the electron density fit of a water (range 0–1.2); lower SASA tracks water buriedness (range 0–98.5 Å2); and higher HB count captures some of the enthalpic contributions of water–protein interactions (range 0–4). This provides statistically compelling evidence that the median of each metric differs between cryo waters present and not present at RT (P < 2.2 × 10−16).

To examine the ability of each metric to rank water molecules from most likely to be present at RT to least likely, we performed a receiver operating characteristic (ROC) analysis. ROC curves for each metric show predictions of the binary classification of present at RT (positive) or not present at RT (negative) (Fig. 1c). Each metric predicted the classification better than random, with an area under the ROC curve (AUC) >0.5 (AUC values of 0.747 for EDIA, 0.734 for Bnorm, 0.725 for RSCC, 0.720 for SASA and 0.672 for HB).

Finally, to predict the probability between 0 and 1 of finding each water at RT, we combined all five metrics into a unified prediction by training a random forest machine learning model. Henceforth, we refer to this predictive value as the ColdBrew probability. ROC curves for the training, test and overall sets showed similar AUC values of 0.810, 0.799 and 0.804, respectively (Fig. 1c), with an accuracy of 73% for the training set and 72% for the test set. Binary predictions showed good, consistent predictions for the training, test and overall sets (Supplementary Table 1). There was a good balance between the classes (present at RT versus not present at RT) and a good balance of precision and recall for each set; a cutoff of 0.5 yielded the highest overall accuracy and the best balance between precision and recall (Supplementary Fig. 2). The model importance scores of each metric from the trained model largely followed the same order as for the individual AUCs (Fig. 1c and Supplementary Fig. 3).

ColdBrew best predicts ligand-proximal waters

For the subset of liganded holo proteins, our predictions were best for the 2256 water molecules within 7 Å of ligand atoms, with an AUC of 0.838 (Fig. 1c). This is encouraging, as these sites are most interesting to drug designers. For the waters near ligand, the results were essentially the same between train and test as well, as shown in Supplementary Fig. 4. Figure 1d shows examples that demonstrate the utility of our method near bound ligands. Cases in which all cryo waters are retained at RT should yield high ColdBrew probabilities. For instance, SiaP has a strong, well-connected water network, with all 18 cryo water molecules appearing at RT (cryo: 2XA5, RT: 6H75). Indeed, our method predicted that 18/18 waters would have high ColdBrew probabilities ≥0.63: 72% had probabilities >0.8, and 44% had probabilities >0.9. Another example in which all binding-site cryo waters are retained at RT is ligand-bound GFP (cryo: 6QQ8, RT: 6QQA); all six cryo water molecules near the ligand have a ColdBrew probability ≥0.82. The same six water molecules have been described as forming a tight, ‘ice-like’ H-bond network that opposes isomerization of the chromophore32.

However, water networks are typically more heterogeneous between temperature classes. For instance, among ten binding-site cryo waters in holo trypsin, only seven appeared at RT (cryo: 7BSA, RT: 7WB9). ColdBrew correctly predicted the seven highest-probability waters to be present and the three lowest-probability waters not to appear at RT. Previous analyses of water roles in trypsin have identified additional waters flanking the binding pocket33 (Supplementary Fig. 5, waters 1–7). Expanding our analysis to these waters confirmed that the apo waters known to be more mobile and suboptimally bonded33 did indeed have lower ColdBrew probabilities. Another example is holo CypD, in which four out of six cryo waters present at RT had the highest ColdBrew probabilities, whereas the two lowest-probability waters were absent at RT (cryo: 3RDC, RT: 4ZSC). Although our ranking was ideal in these two cases (AUC 1), we typically expect our predictions of high- and low-probability waters to be more informative than those for the middle range. Revisiting the Hsp90α example (Fig. 1a), we find that 3/13 cryo waters for holo Hsp90α (cryo: 3QDD, RT: 7S9G) and 6/16 for apo Hsp90α (cryo: 6GPW, RT: 7S8Z) do not appear at RT. For both examples, although all predictions with ColdBrew probabilities ≤0.35 (2/2) or ≥0.65 (16/16) were correct, the respective AUCs of 0.87 and 0.91 indicated uncertainty in the middle range. Another example showing that ColdBrew is most predictive in these high and low ranges is holo FABP3, which has 14/17 cryo waters conserved across temperatures (cryo: 7WD6, RT: 7FCX); of the 4 lowest-probability waters, 3 were correctly not present at RT (AUC 0.94). Across all these examples, predictions were best in the ColdBrew probability ranges of 0–0.35 and 0.65–1, with accuracies of 83% (5/6) and 98% (54/55), respectively.

To test the generality of this observation and provide practical guidance on confidence ranges for high- and low-probability waters, we expanded our approach to the entire dataset of 242 diverse temperature pairs. We found that 84% of waters present at RT had ColdBrew probabilities >0.35 and that 91% of waters not present at RT had ColdBrew probabilities <0.65. Unsurprisingly, the binary accuracy was high in the ranges 0–0.35 and 0.65–1, at 80% and 82%, respectively, but it was only 58% in the 0.35–0.65 range.

Low-probability waters have greater positional uncertainty

Another trend we observed anecdotally in several examples, including CypD, FABP3 and trypsin, was that high-probability waters were less mobile than others (Fig. 1d). For instance, four temperature-conserved waters in CypD with ColdBrew probabilities of 0.94, 0.78, 0.62 and 0.25 moved by 0.1, 0.1, 0.3 and 0.6 Å, respectively. In trypsin, two waters with lower probabilities were previously noted to shift between temperatures33 (Supplementary Fig. 5). Again, we sought to confirm whether lower positional variation generally coincided with higher ColdBrew probability across the entire set of diverse proteins. Therefore, for all waters present at both cryo and RT, we calculated the distance shifted upon temperature change. Figure 1e shows box plots for these distances across various ColdBrew probability ranges. As hypothesized, high-probability waters shifted less (a median of 0.23 Å for the 0.9–1.0 range) than low-probability waters (a median of 0.89 Å for the 0–0.2 range) and had less variance (standard deviations of 0.21 Å for 0.9–1.0 and 0.67 Å for 0–0.2, respectively) (P < 2 × 10−16). This confirms that the positional uncertainty of a water generally increases as its ColdBrew probability decreases.

High-probability waters are less displaceable by ligands

Next, we hypothesized that waters present at both temperatures were generally more stable and, thus, less likely to be displaced by ligands. Figure 1f shows an example comparing a set of three congeneric ligands bound to SiaP, with the smallest ligand, 2-keto-3-deoxynononic acid (KDN), having a water network with relatively strong waters. Assuming that the water least likely to be present at RT (ColdBrew probability of 0.3) was the most displaceable, this may indicate a plausible direction for ligand growth. Indeed, larger substrates, such as Neu5Ac and Neu5Gc, displace that water. Functionally, this enables secondary transporters, such as tripartite ATP-independent periplasmic transporters, to expand their range of host-derived substrates for more effective immune evasion and host colonization by pathogenic bacteria5,34. Our method even detected subtle differences in the trypsin example (Supplementary Fig. 5). Among seven waters present at both temperatures, two with lower relative probabilities were indeed ligand displaceable, whereas the surrounding waters remained untouched.

To test the generality of this hypothesis, we assembled a dataset of 162 groups totaling 5,604 cryo structures, with each group composed of at least 10 structures with a ligand bound in the same pocket, as determined by clustering ligand center of mass. For example, PTP1B had four different ligand-binding locations, corresponding to four different groups. We then clustered all waters within each group. To quantify water conservation and displaceability within the dataset, we calculated three metrics: the fraction occupied by water, the fraction missing and the fraction displaced (by ligand or protein); fractions are the percentage of PDB structures within a group that had a water with the respective feature in that cluster. For 101 of these groups, there was also at least 1 matching apo structure available. Here, we assumed that, for instance, all 457 deposited apo structures of endothiapepsin were truly ligand-free, although negative hits from a fragment screen35 may indeed be overlooked holo structures36. Despite this caveat, the utility of these structures within the scope of our method is unaffected, as it enables any modeled water to be analyzed. For each water in this dataset, we also calculated our ColdBrew probability metric.

We used the three cluster-based metrics to classify waters near the ligand-binding pocket as conserved, displaceable or absent. This yielded 11,981 holo and 8,409 apo waters that were ‘conserved’, that is, that were never displaced and were present in more than 90% of structures within the group; 23,299 holo and 2,834 apo waters that were considered ‘absent’, that is, missing from more than 60% of structures; and 20,395 holo and 6,435 apo total waters that were ‘displaceable’, that is, that were displaced by ligand or protein in more than 10% of grouped structures. Across the dataset, our ColdBrew probability metric separated the distributions of displaceability categories well, as shown in Fig. 2a for all apo and ligand-bound groups (P < 2 × 10−16). This confirmed that conserved waters often have high ColdBrew probability, whereas displaceable or absent ones are typically predicted to be less likely to be present at RT. Said differently, waters with low ColdBrew probability are generally easier to displace by ligands than are waters with high ColdBrew probability.

Fig. 2: ColdBrew probability links to ligand displaceability and energy among cryo X-ray structures.

a, The distribution (box and kernel density estimate plots) of ColdBrew probability for apo and ligand-bound water molecules shows different distributions for three displaceability categories: conserved, displaceable and absent (n = 8,409, 6,435, 2,834, 11,981, 20,395 and 23,299, respectively). b, An example showing waters in Hsp90α colored by ColdBrew probability for apo and ligand-bound structures, with box plot distributions for both (the axes are the same as in a) (n = 10, 36, 27, 680, 1,170 and 1,395, respectively). A total of 234 known Hsp90α ligands are overlaid in black lines for orientation, and conserved waters are labeled as A–D. c, Correlation between GIST energy (in kcal mol−1) and apo ColdBrew probability combined for all systems (Hsp90α, endothiapepsin, β-lactamase and PTP1B). d, The distribution (box plots and kernel density estimate plots) of GIST energy for apo and ligand-bound water molecules for three displaceability categories: conserved, displaceable and absent (n = 1,430, 745 and 100, respectively). For a, b and d, box plots were made using the following parameters: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, outliers.

Full size image

Figure 2b shows a side view of the Hsp90α binding pocket with all bound ligands superimposed on the apo and holo waters. Box plots show similar separation between the categories for this system. Because there are many ligand-bound waters, it is easy to see how the highly conserved waters correlate with our metric. Notably, those same waters are also conserved in the apo structures. Waters A–D, in particular, are conserved across nearly all Hsp90α structures, including across all four human isoforms7.

ColdBrew probability and displaceability link to GIST energy

With the relation between ColdBrew probability, water conservation and ligand displaceability established, we next hypothesized that our metric would provide a simple, intuitive proxy for water energetics in ligand discovery. Specifically, to establish a link between ColdBrew probabilities and water energetics, we ran molecular dynamics (MD) simulations and generated GIST grids covering the entire simulation box for four distinct systems: endothiapepsin, Hsp90α, PTP1B and β-lactamase. For each, we ran independent calculations for two to four available apo starting structures. We then identified hydration sites from the MD + GIST grids based on areas of high grid density and calculated Gaussian-smoothed energy for these sites by using WaterKit23. Through clustering, we matched these hydration sites to the apo waters, when possible, and observed a moderate correlation (r = −0.60, with r = −0.62 for waters near ligand) between the MD + GIST energy and the ColdBrew probability of the apo crystal structure (Fig. 2c).

Next, we combined our displaceability analysis with the GIST energies for all hydration sites near ligands. The GIST energies separated the distributions of the displaceability categories well (Fig. 2d); P values indicated that the medians for all pairs of categories were different, with pairwise P values <2.1 × 10−7.

The fact that our easy-to-calculate metric (Fig. 2a) separates waters in these categories similarly to computationally expensive MD simulations (Fig. 2d) suggests that ColdBrew is a useful proxy for advising ligand discovery efforts that consider water displaceability.

Ligand discovery empirically avoids high-probability waters

Our ColdBrew metric generally separates waters that are conserved, displaceable and absent (Fig. 2). We will now illustrate the utility of the approach, using specific examples. For Hsp90α, box plots show a clear separation of categories (conserved, displaceable or absent) by ColdBrew probability for ligand-bound or apo waters or by GIST energy (Fig. 3a,b). To visualize the correlation of energy, ColdBrew probability and displaceability, we scale waters within the binding site by energy and color them by probability (Fig. 3c) or category (Fig. 3d). These colors match those in the box plots; each water corresponds to a black point on the scatterplot of energy versus ColdBrew probability in Fig. 3b. For positive correlations, we would expect to see two main scenarios. First, waters with favorable energy would have high probability and be conserved, that is, the large blue waters in Fig. 3c would match the large purple waters in Fig. 3d. Second, unfavorable waters would have low probability and would be displaced or absent, that is, the small red waters in Fig. 3c would match the small yellow/green waters in Fig. 3d. We see many instances of this correlation, for example, agreement between ColdBrew probability and energy as a large blue water in grid-cell J11, a medium white water in E10 and a small red water in D5. These waters then can be compared with the categorical coloring (J11 is conserved, E10 is displaceable and D5 is absent).

Fig. 3: Hsp90α ligands are empirically designed to avoid high-probability waters.

a, Results for individual systems: box plots show the distributions of ligand-bound ColdBrew probability, apo ColdBrew probability and GIST energy for the three displaceability classifications (n = 680, 1,170, 1,395, 10, 36, 27, 10, 50 and 90, respectively). Box plots were made using the following parameters: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, outliers. b, A scatterplot showing the relation between GIST energy (in kcal mol−1) and apo ColdBrew probability. c, apo waters that had corresponding GIST energy (that is, the same waters as the black points in b) colored by ColdBrew probability and sized by GIST energy (see key for a–d). A dark-gray surface representation of all ligands is shown along with a representative cartoon structure. d, The same waters colored by category and sized by GIST energy. e, Hsp90α apo water network for a cryo–RT pair in two views (side and front), with RT waters shown in pink and cryo waters colored on the basis of RT identity and labeled by predicted ColdBrew probability. f, Coloring apo waters by ColdBrew probability shows how the designed ligand PU3 extends the adenine (ADE) building block in two directions along paths lined by low-probability waters. A high-probability water (circled) is displaced if a polar group overcomes the displacement penalty by gaining a favorable interaction. g, These two paths are recapitulated by waters from two different apo structures near any ligand; see i for an analysis of boxed loop-responsive water. h, The overlay of all 234 observed ligands, shown in black, that empirically explore the two paths of low-probability waters. i, Protein conformational changes modulate the ColdBrew probability and, consequently, the GIST energy of the boxed water (white square in g) in the loop-in structure and the helical structure, as compared with the apo state.

Full size image

Next, we map ColdBrew probabilities onto Hsp90α apo waters (cryo: 6GPW, RT: 7S8Z) and color them by classification (Fig. 3e). This highlights that waters that are absent at RT (purple) are often displaced upon ligand binding. For example, a building block of known Hsp90α binders (adenine, ADE) superimposes mainly on non-RT waters. The exception is one circled high-probability water (Fig. 3f; cf. the water in grid-cell K10 in Fig. 3c,d) that results from the ligand paying the energetic cost of displacing this water by forming favorable interactions with the polar nitrogen. Our ColdBrew metric lays out two viable paths to build off this ligand by replacing low-probability waters. Remarkably, a designed ligand PU3 indeed proceeds down both of these paths. Likewise, low-probability waters line these 2 paths in 2 Hsp90α apo structures (Fig. 3g) and when all 234 ligand-bound crystal structures of Hsp90α are overlaid (Fig. 3h)37. Remarkably, all Hsp90α PDB ligands that were mostly designed in the absence of water information in tedious trial-and-error medchem cycles empirically follow these two paths elucidated by our method.

Protein conformational changes modify water displaceability

Ligand binding can cause conformational changes affecting water displaceability. For instance, Hsp90α binds ligands mainly via three discrete lid conformations: loop-in, loop-out and helical37. In >200 available structures of Hsp90α ligands, a particular water molecule (boxed in Fig. 3g) was conserved in 88% of loop-in structures and was never displaced by ligand. The same water was never present in loop-out structures. However, in the helical conformation, this water was displaced by a ligand in 18% of structures and was conserved in 71%. Given that this water was strong in the apo (7S8Y, loop-in, ColdBrew probability of 0.96, energy of −8.33 kcal mol−1), we would have expected it to remain untouched by ligands. This demonstrates how changes in the protein energy landscape, such as conformational changes, can lower the water’s mean ColdBrew probability, in this case from 0.87 to 0.74 because of the transition from the loop-in to the helical conformation, thereby facilitating its displacement (Fig. 3i).

Polar ligand atoms can offset penalties of stable waters

Another example of a correlation between ColdBrew probability distributions based on category and energy is β-lactamase (Fig. 4a–d). apo waters from 119 β-lactamase structures are superimposed with ligands and a catalytic serine, the covalent attachment point for some ligands (see Fig. 4e for H7M and Fig. 4f for all). Exemplified by the nanomolar noncovalent inhibitor H7M, ligands typically avoid high-probability waters, except for the water coordinated by the covalent serine (circled). The ligand replaces the stable serine-proximal water with a polar oxygen atom. Expanding this analysis to all 119 known ligands, both covalent and noncovalent ligands bind mainly in low-probability water areas (Supplementary Fig. 6). However, as for Hsp90α, some high-probability waters are displaced by groups for select ligands (circles i–v). The one most frequently displaced is the serine-adjacent water (ii). Our analysis shows that these high-probability waters are replaced exclusively by polar ligand groups (Fig. 4g).

Fig. 4: Ligands empirically avoid high-probability waters for β-lactamase.

a, β-Lactamase box plots of the distributions of ligand-bound ColdBrew probability, apo ColdBrew probability and GIST energy over the three displaceability classifications (n = 206, 250, 179, 14, 45, 19, 14, 110 and 167, respectively). Box plots were made using the following parameters: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, outliers. b, A scatterplot showing the relation between GIST energy and apo ColdBrew probability. c, apo waters (cf. black points in the scatterplot) colored by ColdBrew probability and sized by GIST energy. A dark-gray surface representation of all ligands is shown, along with a representative cartoon structure. d, The same waters colored by displaceability category and sized by GIST energy. e, β-Lactamase examples with apo water networks overlaid with H7M, a nanomolar-affinity noncovalent ligand, colored by ColdBrew probability at two angles and with the catalytic serine shown in green. f, Ligands avoid high-probability waters unless they can displace them with a polar group; see the circled waters. Overlay of 119 ligands. g, The circled waters in e and f are displaced by polar groups in the ligands. h, Context-dependent distribution of ligand-bound waters that are displaced by other ligands in more than half of the structures. The separation of the red and gray distributions shows that polar atoms are better able to displace high-probability waters than are nonpolar ones (P < 2 × 10−16).

Full size image

Finally, to determine whether this observation generally held true, we expanded our analysis to the entire ligand-grouped set of 5,604 structures (cf. data in Fig. 2). First, we tracked waters that were displaced by a ligand atom in >50% of the structures (Fig. 4h). We then split these waters into subsets in which the water was displaced by a polar versus nonpolar atom more than half of the time. Notably, both distributions are clearly separate and are skewed toward high ColdBrew probability (median 0.71) for the polar-displaced waters, whereas the nonpolar-displaced waters are skewed toward lower probability (median 0.38) (P < 2 × 10−16). This highlights the fact that ligands empirically follow low-probability waters and, for the exceptions, ligands generally displace a high-probability water with an equivalent polar ligand group to pay the energetic displacement cost.

Discussion

In this work, we have explored a simple experimental proxy of water displaceability from crystal structures. Three main implications for exploring structural water networks in ligand discovery have emerged. First, ColdBrew can predict the probability of cryo waters occurring at RT. Second, water molecules with high ColdBrew probability are typically less displaceable by ligands than are low-probability ones. Third, water displaceability and ColdBrew probability correlate with the GIST energy of a given water. These results encourage the utility of ColdBrew in leveraging water networks observed in cryo crystal structures for ligand discovery in lieu of expensive computational methods.

Our previous findings that water networks differ with temperature15 prompted the question as to whether we could predict which cryo waters would also be present at RT. We found that five metrics captured key aspects of what makes a cryo water (dis)appear at RT. Notably, information hiding in plain sight in the electron density map (EDIA and RSCC) or in the model (Bnorm) carried more information than the typically used measures of burial and hydrogen bonds. This should encourage users of structural data to focus on the primary data (the maps) instead of relying solely on map interpretations (the models). Predictions were best for water molecules near ligands, which is good news as binding sites are generally of greatest interest in ligand discovery and for understanding protein function38. It is tempting to speculate that this bias is introduced as a result of crystallographers paying most attention to carefully modeling binding-site waters. If our hypothesis is correct, our results should motivate improvements in water modeling in the field. In turn, improving our ability to evaluate solvation patterns may provide new hints of assessing the ligandability of protein sites and the potential tractability of targets for drug discovery39. Furthermore, as machine learning methods continue to increase in power and ubiquity, relying on input data quality, it is key to train these methods on correctly modeled dynamic features such as waters and other transient features. Also, as not all parts of the protein receive equal attention in current water modeling approaches, other areas of potential interest, for example, hidden allosteric sites14, rely mostly on automatic water-picking methods40,41 and may escape detection by downstream computational methods42. Although these problems appear to be present in current PDB datasets, the observed correlation between GIST-calculated energy and ColdBrew probability holds surprisingly well across the entire protein surface (r = −0.6) (Fig. 2c).

Our results support high confidence in the ColdBrew probability ranges of 0–0.35 and 0.65–1, for which approximately 80% of waters were predicted correctly. We recommend that, instead of interpreting the results exclusively as binary predictions, users should consider relative rankings of waters. The gray zone of waters with probabilities near 0.5 will include waters from both categories, present or absent, depending on the context of the binding site.

Next, we showed that waters not present at RT were generally more displaceable. Although this is indicated by the separation of the interquartile range and median of the distributions of the box plots in Fig. 2a, the data range ‘whiskers’ for displaceable waters span the entire spectrum of ColdBrew probabilities. This highlights the fact that likely no water is irreplaceable, which makes functional sense considering the high context dependence of water energetics that drives biology. For instance, if ligands or protein sidechains can pay the energetic price, they will displace conserved waters despite their high ColdBrew probability, especially if a strong interaction is gained upon displacement. So, more good news for ligand discoverers, although, in practice, this may be difficult to achieve and anticipate.

These findings indicate three ways by which ligand discoverers can approach a high-probability water. Our guidelines are to (1) leave the high-probability water untouched and focus on building ligands along paths of low-probability waters; (2) displace the high-probability water with an appropriate polar ligand group that ‘pays’ the energetic displacement cost; or (3) alter the protein conformation to reshape the energetic landscape of the water network. The latter may lower the water displacement cost and enable ligands to bind in unexpected regions with high-probability apo water networks. While one cannot expect to easily predict protein changes prospectively, the Hsp90α example demonstrates how our method can track the consequences of such changes retrospectively (Fig. 3i).

The correlation when linking water conservation and ColdBrew probability to GIST energy is striking, considering the caveats regarding both current computational uncertainties and experimental modeling mishaps. For computational methods, it is worth remembering that GIST calculations do not represent the ground truth. Experimentally, an intrinsic caveat of a water-focused approach is that water modeling remains arduous. First, water placement is resolution dependent. Water visibility increases with resolution; for instance, at 2-Å resolution, we expect approximately 1 water per amino acid, whereas at 1.0-Å resolution, that ratio increases to approximately 1.6 (ref. 43). Resolutions worse than 2.5 Å often make it impossible to model water into electron density maps with confidence. To minimize the risk of confusing genuinely missing waters with those that are not visible because of lower resolution, we kept the median resolution differences tight for temperature pairs (0.16 Å). Another pitfall is the presence of mismodeled water density blobs that are ignored or mistaken for crystallographic additives such as buffer solute, ions and cryo protectants, or vice versa. Generally, careful modeling will distinguish a metal ion from a water, for example, on the basis of B factors, but, at the PDB scale, we expect there to be many such mishaps, compromising correlative analyses such as our ColdBrew approach42,44. Water networks are vulnerable to experimental conditions. Preventing desolvation during crystal handling and data collection, for instance by using oils45 to coat crystals or by controlling humidity levels46, is important to ensure the observed difference are due to temperature, not evaporation.

Finally, sidechain and ligand flexibility with alternate conformations and, therefore, occupancies <1 interfere with the interpretation of water positions47. Although there are tools for modeling alternate conformations of sidechains into electron density maps48,49, water molecules with correlated alternate conformations still require cautious manual annotation. For these reasons, even crystal structures of the same system may be unexpectedly heterogeneous, causing seemingly random differences from the perspective of water networks. Considering these notable caveats and the fact that all our metrics rely on correct water modeling, it is remarkable that our approach worked so well. It also suggests that better predictions may be obtained for more carefully modeled waters in the future, either by improving the water modeling toolbox or by cross-pollinating13 information on water networks across well-studied systems.

Overall, these results illustrate the ability of ColdBrew predictions to serve as an intuitive proxy for water energy straight from experimental cryo structures, in lieu of computationally expensive methods. This provides guidance on which waters should be displaced or left untouched during ligand design and discovery. Understanding water conservation at physiological temperatures is also important to understanding a protein’s function, for instance to identify highly conserved waters that participate in hydrolysis of catalytic intermediates38. To make the predictions widely accessible, we have provided ColdBrew probability predictions for more than 46 million waters in more than 107,000 PDB cryo structures in addition to making the code freely available50,51.

Methods
Temperature-matched dataset

Assembly of our dataset of cryo–RT pairs began with a snapshot of the PDB from October 2023. The basic criteria for considering PDB structures were as follows: protein as the only polymer, refinement resolution ≤2.5 Å, data collected via X-ray diffraction, experimental data available and r-free (an established cross-validation residual calculated from a subset of reflections excluded from refinement) ≤0.3. Only monomers and homomers were considered. Some other structures were excluded because of problems with the metadata deposited in the PDB, such as multiple listed temperatures. Structures without water or with D2O present were also discarded. The diffraction source temperature criteria for RT and cryo structures were ≥270 K and <150 K, respectively. Cryo PDB structures were matched to all available RT PDB structures based on (1) an exact sequence match, but allowing tails to differ in length and ignoring the first and last 5 residues in the chain; (2) a space group match; (3) resolution within 0.55 Å; (4) Cα r.m.s.d. ≤1.5 Å; (5) a chain identity match within the asymmetric unit; and (6) a ligand identity match. The entire asymmetric unit was used for these analyses. Rather than forcing all ligands to match, because there may be crystallographic additives of no interest, this study focused on a subset enriched in more biologically relevant ligands. As this type of approach has been used previously52, ligands (for the sole purpose of matching of pairs) from the AlphaFold ligand exclusion list, ions and standard residues were ignored53. Based on a manual inspection, the following ligands were also ignored: PG5, XE, AR, NO, KR, 2PE, SAC, XYP, OXY, 1PG, 2FU, CEM, P2K, TLA, KCX, XYS, P4K, SF4, NH3, CSO, FRU, 4OS, BMA, TYS, GAL, GLC, 12P, CME, PTR, FUC, MAN, NGA, NAG and BGC. To reduce redundancy, these pairs were filtered down to one pair per combination of sequence, space group and ligand(s) by selecting the RT structure with the best resolution and the cryo structure with the resolution closest to that of the RT structure. A flowchart illustrating the temperature-matched dataset construction is provided in Supplementary Fig. 7a. This procedure resulted in 242 cryo–RT matched pairs (listed in Supplementary Data 1). Some systems had multiple pairs (matched with a different space group or ligand(s)), but no single system predominated. PTP1B contributed the most cryo waters, representing <7% of the total, with most systems (130/155) contributing <1% of the total (Supplementary Fig. 1d).

To compare the water molecules in the cryo–RT pairs, a clustering approach was used, similar to previous methods54. The oxygen atoms of water molecules were clustered using complete-linkage hierarchical clustering with a radius of 2.4 Å, in an in-house script using the Python package scipy 1.7.3 and sklearn 1.0.2 (ref. 55).

Across all pairs, there were a total of 72,318 cryo and 43,693 RT waters before any filtering. In this study, cryo waters were classified as ‘present at RT’ if present in a cluster also occupied by an RT water and as ‘not present at RT’ if in a cluster not occupied by an RT water. Cryo waters that overlapped (within 1 Å) with any nonwater atom in the RT structure were not considered in this analysis as they could not be reasonably considered missing (908 overlapped with an atom in the RT structure). Furthermore, some RT waters not present at cryo temperature (10,742 in total, 684 of which overlapped with another atom in the cryo structure) were not considered in this study.

Visualization was created using seaborn 0.11.2 and matplotlib 3.5.1. Box plots were made using the following parameters: center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; points, outliers.

Calculation of water metrics in crystal structures

Bnorm was calculated as previously described5, with all protein atoms being used to calculate the mean and standard deviation for normalization. EDIA was calculated using the command-line tool EDIAscorer26,27. The number of hydrogen bonds was calculated using hbplus after adding hydrogens in PyMOL56,57. The HB count used here for each water was the sum of the backbone and sidechain hydrogen bonds to a given water molecule, which was the total number of hydrogen bonds between the water of interest and protein atoms. SASA (absolute) was calculated using the naccess program58. RSCC was calculated using phenix.real_space_correlation41. EDIA could not be calculated for 268 of the cryo waters, so the final counts were 32,901 present at RT and 38,244 not present at RT. Unless otherwise specified, default parameters were used.

Significance testing

To calculate P values for the continuous variables (Bnorm, EDIA, RSCC and SASA), the R package ‘clusrank’ was used to fit pairwise Wilcoxon rank-based tests that account for clustering, with clusters defined by a PDB structure or a pair of PDB structures, followed by Holm adjustment for multiple testing59,60. To calculate P values for the discrete variable (HB), the R package ‘htestClust’ was used to fit a cluster-weighted chi-squared test of independence, with clusters being defined by the pair of PDB structures61,62.

Model training and analysis

Starting with our dataset of 71,145 classified cryo waters, machine learning models were trained to predict whether cryo waters would appear at RT, using the Python package sklearn55. Because there was a slight class imbalance between ‘present at RT’ (32,901) and ‘not present at RT’ (38,244), Synthetic Minority Over-sampling Technique was first performed, using the Python package imblearn63. Next, the dataset was split into training and test datasets, randomly using a 75:25 train-to-test ratio. A random forest classification model was trained with a maximum depth of 7. We define the result of the model as the ColdBrew probability, that is, the probability of a cryo water being present at RT. The code to calculate the metrics and ColdBrew probability is available via GitHub at https://github.com/TheFischerLab/ColdBrew. The exact values of the result may vary depending on compiler, software versions or choices made in experimental data acquisition. For the results here (which are also available via Zenodo at ref. 50), we used the following software versions: PyMol 2.5.0, Phenix 1.20.1, Naccess 2.1.1, HBPLUS 3.2 (compiled on Mac M1) and EDIAScorer 1.1.0. Furthermore, there were a total of 2,256 waters within 7 Å of ligands (of those, 1,691 were in the training set and 565 were in the test set).

Ligand-grouped dataset

The dataset used to quantify displaceability for cryo structures again began with the following criteria: cryo temperature (diffraction source temperature <150 K), protein as the only polymer, refinement resolution ≤2.5 Å, data collected performed using X-ray diffraction, experimental data available, and r-free ≤0.3. Again, only monomers and homomers were considered. The ‘A chain’ in the .pdb file was considered, if available. Some other structures were excluded because of problems with the metadata deposited in the PDB, such as multiple listed temperatures. Structures without water or with D2O present were also discarded. The same ligand criteria were used as for the temperature-matched dataset. Only structures with biologically relevant ligands were considered, and structures were removed if they contained any other nonprotein or nonwater atoms within 7.5 Å of the ligand so as not to interfere with the analysis. Furthermore, only structures with one copy of a ligand were allowed. Groups were assembled on the basis of UniProt ID alone (that is, without considering the exact sequence match or space group). The Cα r.m.s.d. of each group member was required to be ≤1.5 Å to an arbitrary group representative. Groups were further segmented by clustering the ligand center of mass (complete-linkage hierarchical clustering with a radius of 10 Å) to compare only ligands bound in similar locations. Only groups with at least ten members were analyzed. A flowchart illustrating the ligand-grouped dataset construction is provided in Supplementary Fig. 7b. This procedure resulted in 162 groups, accounting for 5,604 cryo structures. The PDB IDs for this dataset are listed in Supplementary Data 2.

A search was also performed for apo structures, defined as meeting all match criteria listed above except the ligand match (and having no other ligands based on our definition). However, each possible apo PDB structure was explicitly checked to ensure there were no nonwater or nonprotein atoms within 7.5 Å of the ligand in the holo structure. From our definition, 101 groups had at least 1 apo structure, for a total of 1,277 apo structures. The PDB IDs for this dataset are listed in Supplementary Data 3. Note that our apo analysis is dominated by endothiapepsin, with a total of 457 apo PDB structures and 3 endothiapepsin groups.

To quantify the conservation and displaceability of waters, all water molecules in the ligand-bound structures within each group were first clustered using the same method as described above. Then, three metrics were calculated for each cluster: the fraction of PDB structures with a water in that cluster (the fraction occupied by water), the fraction of PDB structures without a water in that cluster (the fraction missing) and the fraction of PDB structures in which the cluster center overlapped a ligand or protein atom (the fraction displaced). From these metrics, waters near the ligand binding pocket (within 4 Å of a ligand) were classified as ‘conserved’ (fraction occupied by water >0.9, fraction displaced 0), ‘absent’ (fraction missing >0.6) or ‘displaceable’ (fraction displaced >0.1). We note that these categories are not mutually exclusive, but the fractions do add up to 1 for each water. It is possible to have a water that is both absent and displaceable (for example, if fraction occupied by water is 0.1, fraction displaced is 0.2 and fraction missing is 0.7) and some waters are unclassified (for example, if fraction occupied by water is 0.5, fraction displaced is 0.0 and fraction missing is 0.5).

MD and GIST energy calculations

MD simulations for four systems (endothiapepsin, Hsp90α, PTP1B and β-lactamase) were performed using various apo crystal structures as input. The PDB IDs used as starting points can be found in Supplementary Data 3. All calculations were performed using the CHARMM36m force field64. To prepare each system, all nonprotein atoms were removed from the .pdb file. Visual molecular dynamics (VMD)65 was used to solvate each apo structure in a 20-Å padded TIP3P water box, via an in-house tcl script. The systems were not ionized, so as not to interfere with water networks, as is standard for this type of analysis31. MD simulations were performed using NAMD66 (version 3.0 with graphics processing unit acceleration). The systems were minimized and equilibrated in two initial steps each. The first minimization stage was 20,000 steps, with all protein atoms being restrained using a force constant of 10 kcal mol−1 Å−2, that is, minimizing solvent only. The second minimization stage was also 20,000 steps, but with a restraint force constant of 5 kcal mol−1 Å−2 for protein heavy atoms. In the first equilibration stage, the system was heated to a temperature of 300 K by systematically weakening the restraints, for a total of 380 ps under npt (constant number of particles, pressure and temperature) conditions. All subsequent calculations were performed with a restraint force constant of 5 kcal mol−1 Å−2 for protein heavy atoms. The second equilibration stage was a 20-ps npt simulation. Before the production run, the system was further equilibrated with 5-ns npt simulations and 5-ns nvt (constant number of particles, volume and temperature) simulations. Finally, the production run was a 20-ns nvt simulation. GIST grids were then calculated from the production trajectories, using cpptraj V5.1.0 (ref. 67) with grid dimensions large enough to cover the entire protein, with padding. Grid spacing was 0.5 Å. Waterkit23 was used to analyze the GIST grids from the MD production run. Hydration sites and Gaussian-smoothed energies for those hydration sites were calculated using default parameters.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability

Dataset tables containing PDB IDs for cryo–RT pairs are provided in Supplementary Data 1. The cryo PDB IDs for the ligand displacement dataset (holo and apo) are also included in Supplementary Data 2 and 3, respectively. Precalculated ColdBrew probabilities for over 107,000 PDB structures are available via Zenodo at https://doi.org/10.5281/zenodo.13909324 (ref. 50).

Code availability

The code for calculating ColdBrew probabilities is available via GitHub at https://github.com/TheFischerLab/ColdBrew.

References

Lee, Y.-H. et al. Crystal structure of a trapped phosphoenzyme during a catalytic reaction. Nat. Struct. Biol. 4, 615–618 (1997).

Article
 
CAS
 
PubMed
 
Google Scholar

Keedy, D. A. et al. Crystal cryocooling distorts conformational heterogeneity in a model Michaelis complex of DHFR. Structure 22, 899–910 (2014).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Ho, J. D. et al. Crystal structure of human aquaporin 4 at 1.8 Å and its mechanism of conductance. Proc. Natl Acad. Sci. USA 106, 7437–7442 (2009).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Bronstein, A. & Marx, A. Water stabilizes an alternate turn conformation in horse heart myoglobin. Sci. Rep. 13, 6094 (2023).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Darby, J. F. et al. Water networks can determine the affinity of ligand binding to proteins. J. Am. Chem. Soc. 141, 15818–15826 (2019).

Article
 
CAS
 
PubMed
 
Google Scholar

Sleigh, S., Seavers, P., Wilkinson, A., Ladbury, J. & Tame, J. Crystallographic and calorimetric analysis of peptide binding to OppA protein. J. Mol. Biol. 291, 393–415 (1999).

Article
 
CAS
 
PubMed
 
Google Scholar

Stachowski, T. R., Nithianantham, S., Vanarotti, M., Lopez, K. & Fischer, M. Pan-HSP90 ligand binding reveals isoform-specific differences in plasticity and water networks. Protein Sci 32, e4629 (2023).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Doyle, M. D. et al. Water networks in photosystem II using crystalline molecular dynamics simulations and room-temperature XFEL serial crystallography. J. Am. Chem. Soc. 145, 14621–14635 (2023).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Cushing, V. I. et al. High-resolution cryo-EM of the human CDK-activating kinase for structure-based drug design. Nat. Commun. 15, 2265 (2024).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Chari, A. & Stark, H. Prospects and limitations of high-resolution single-particle cryo-electron microscopy. Annu. Rev. Biophys. 52, 391–411 (2023).

Article
 
CAS
 
PubMed
 
Google Scholar

Fischer, M. Macromolecular room temperature crystallography. Q. Rev. Biophys. 54, e1 (2021).

Article
 
PubMed
 
Google Scholar

Fraser, J. S. et al. Accessing protein conformational ensembles using room-temperature X-ray crystallography. Proc. Natl Acad. Sci. USA 108, 16247–16252 (2011).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Keedy, D. A. Journey to the center of the protein: allostery from multitemperature multiconformer X-ray crystallography. Acta Crystallogr. D 75, 123–137 (2019).

Article
 
CAS
 
Google Scholar

Fischer, M., Shoichet, B. K. & Fraser, J. S. One crystal, two temperatures: cryocooling penalties alter ligand binding to transient protein sites. ChemBioChem 16, 1560–1564 (2015).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Stachowski, T. R., Vanarotti, M., Seetharaman, J., Lopez, K. & Fischer, M. Water networks repopulate protein–ligand interfaces with temperature. Angew. Chem. Int. Ed. 61, e202112919 (2022).

Article
 
CAS
 
Google Scholar

Dunge, A. et al. Exploring serial crystallography for drug discovery. IUCrJ 11, 831–842 (2024).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Lazaridis, T. Inhomogeneous fluid approach to solvation thermodynamics. 1. Theory. J. Phys. Chem. B 102, 3531–3541 (1998).

Article
 
CAS
 
Google Scholar

Michel, J., Tirado-Rives, J. & Jorgensen, W. L. Energetics of displacing water molecules from protein binding sites: consequences for ligand optimization. J. Am. Chem. Soc. 131, 15403–15411 (2009).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Li, Z. & Lazaridis, T. Computing the thermodynamic contributions of interfacial water. Comput. Drug Discov. Design 819, 393–404 (2012).

Article
 
CAS
 
Google Scholar

Bayden, A. S., Moustakas, D. T., Joseph-McCarthy, D. & Lamb, M. L. Evaluating free energies of binding and conservation of crystallographic waters using SZMAP. J. Chem. Inf. Model. 55, 1552–1565 (2015).

Article
 
CAS
 
PubMed
 
Google Scholar

Raymer, M. L. et al. Predicting conserved water-mediated and polar ligand interactions in proteins using a K-nearest-neighbors genetic algorithm. J. Mol. Biol. 265, 445–464 (1997).

Article
 
CAS
 
PubMed
 
Google Scholar

Sun, H., Zhao, L., Peng, S. & Huang, N. Incorporating replacement free energy of binding‐site waters in molecular docking. Proteins 82, 1765–1776 (2014).

Article
 
CAS
 
PubMed
 
Google Scholar

Eberhardt, J. & Forli, S. WaterKit: thermodynamic profiling of protein hydration sites. J. Chem. Theory Comput. 19, 2535–2556 (2023).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Abel, R., Young, T., Farid, R., Berne, B. J. & Friesner, R. A. Role of the active-site solvent in the thermodynamics of factor Xa ligand binding. J. Am. Chem. Soc. 130, 2817–2831 (2008).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Zsidó, B. Z. & Hetényi, C. The role of water in ligand binding. Curr. Opin. Struct. Biol. 67, 1–8 (2021).

Article
 
PubMed
 
Google Scholar

Meyder, A., Nittinger, E., Lange, G., Klein, R. & Rarey, M. Estimating electron density support for individual atoms and molecular fragments in X-ray structures. J. Chem. Inf. Model. 57, 2437–2447 (2017).

Article
 
CAS
 
PubMed
 
Google Scholar

Nittinger, E., Schneider, N., Lange, G. & Rarey, M. Evidence of water molecules—a statistical evaluation of water molecules based on electron density. J. Chem. Inf. Model. 55, 771–783 (2015).

Article
 
CAS
 
PubMed
 
Google Scholar

García-Sosa, A. T., Mancera, R. L. & Dean, P. M. WaterScore: a novel method for distinguishing between bound and displaceable water molecules in the crystal structure of the binding site of protein–ligand complexes. J. Mol. Model. 9, 172–182 (2003).

Article
 
PubMed
 
Google Scholar

Xiao, W. et al. Predicting conserved water molecules in binding sites of proteins using machine learning methods and combining features. Comput. Math. Methods Med 2022, 5104464 (2022).

Article
 
PubMed
 
PubMed Central
 
Google Scholar

Carugo, O. Correlation between occupancy and B factor of water molecules in protein crystal structures. Protein Eng 12, 1021–1024 (1999).

Article
 
CAS
 
PubMed
 
Google Scholar

Balius, T. E. et al. Testing inhomogeneous solvation theory in structure-based ligand discovery. Proc. Natl Acad. Sci. USA 114, E6839–E6846 (2017).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Helms, V., Straatsma, T. & McCammon, J. A. Internal dynamics of green fluorescent protein. J. Phys. Chem. B 103, 3263–3269 (1999).

Article
 
CAS
 
Google Scholar

Schiebel, J. et al. Intriguing role of water in protein–ligand binding studied by neutron crystallography on trypsin complexes. Nat. Commun. 9, 3559 (2018).

Article
 
PubMed
 
PubMed Central
 
Google Scholar

Mulligan, C., Fischer, M. & Thomas, G. H. Tripartite ATP-independent periplasmic (TRAP) transporters in bacteria and archaea. FEMS Microbiol. Rev. 35, 68–86 (2011).

Article
 
CAS
 
PubMed
 
Google Scholar

Schiebel, J. et al. High-throughput crystallography: reliable and efficient identification of fragment hits. Structure 24, 1398–1409 (2016).

Article
 
CAS
 
PubMed
 
Google Scholar

Mehlman, T., Ginn, H. M. & Keedy, D. A. An expanded trove of fragment-bound structures for the allosteric enzyme PTP1B from computational reanalysis of large-scale crystallographic data. Struct. 32, 1231–1238 (2024).

Article
 
CAS
 
Google Scholar

Stachowski, T. R. & Fischer, M. Large-scale ligand perturbations of the protein conformational landscape reveal state-specific interaction hotspots. J. Med. Chem. 65, 13692–13704 (2022).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Dasgupta, M. et al. Mix-and-inject XFEL crystallography reveals gated conformational dynamics during enzyme catalysis. Proc. Natl Acad. Sci. USA 116, 25634–25640 (2019).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Vukovic, S. & Huggins, D. J. Quantitative metrics for drug–target ligandability. Drug Discov. Today 23, 1258–1266 (2018).

Article
 
CAS
 
PubMed
 
Google Scholar

Murshudov, G. N. et al. REFMAC5 for the refinement of macromolecular crystal structures. Acta Crystallogr. D 67, 355–367 (2011).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Liebschner, D. et al. Macromolecular structure determination using X-rays, neutrons and electrons: recent developments in Phenix. Acta Crystallogr. D 75, 861–877 (2019).

Article
 
CAS
 
Google Scholar

Gao, Y., Thorn, V. & Thorn, A. Errors in structural biology are not the exception. Acta Crystallogr. D 79, 206–211 (2023).

Article
 
CAS
 
Google Scholar

Carugo, O. & Bordo, D. How many water molecules can be detected by protein crystallography? Acta Crystallogr. D 55, 479–483 (1999).

Article
 
CAS
 
PubMed
 
Google Scholar

Zheng, H. et al. CheckMyMetal: a macromolecular metal-binding validation tool. Acta Crystallogr. D 73, 223–233 (2017).

Article
 
CAS
 
Google Scholar

Doukov, T., Herschlag, D. & Yabukarski, F. Instrumentation and experimental procedures for robust collection of X-ray diffraction data from protein crystals across physiological temperatures. J. Appl. Crystallogr. 53, 1493–1501 (2020).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Farley, C., Burks, G., Siegert, T. & Juers, D. H. Improved reproducibility of unit-cell parameters in macromolecular cryocrystallography by limiting dehydration during crystal mounting. Acta Crystallogr. D 70, 2111–2124 (2014).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Biela, A. et al. Impact of ligand and protein desolvation on ligand binding to the S1 pocket of thrombin. J. Mol. Biol. 418, 350–366 (2012).

Article
 
CAS
 
PubMed
 
Google Scholar

Stachowski, T. R. & Fischer, M. FLEXR: automated multi-conformer model building using electron-density map sampling. Acta Crystallogr. D 79, 354–367 (2023).

Article
 
CAS
 
Google Scholar

Riley, B. T. et al. qFit 3: protein and ligand multiconformer modeling for X-ray crystallographic and single-particle cryo-EM density maps. Protein Sci. 30, 270–285 (2021).

Article
 
CAS
 
PubMed
 
Google Scholar

Seffernick, J. ColdBrew probabilities. Zenodo https://doi.org/10.5281/zenodo.13909324 (2024).

Seffernick, J. T. & Fischer, M. ColdBrew. GitHub https://github.com/TheFischerLab/ColdBrew (2025).

Wankowicz, S. A., de Oliveira, S. H., Hogan, D. W., van den Bedem, H. & Fraser, J. S. Ligand binding remodels protein side-chain conformational heterogeneity. eLife 11, e74114 (2022).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Performance and Structural Coverage of the Latest, In-Development AlphaFold Model (Google DeepMind AlphaFold Team & Isomorphic Labs Team, 2023).

Sanschagrin, P. C. & Kuhn, L. A. Cluster analysis of consensus water sites in thrombin and trypsin shows conservation between serine proteases and contributions to ligand specificity. Protein Sci 7, 2054–2064 (1998).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Pedregosa, F. et al. Scikit-learn: machine learning in Python. J. Mach. Learn. Res. 12, 2825–2830 (2011).

Google Scholar

McDonald, I., Naylor, D., Jones, D. & Thornton, J. HBPLUS Computer Program (Department of Biochemistry and Molecular Biology, Univ. College London, 1993).

The PyMOL Molecular Graphics System, Version 3.0 (Schrödinger, LLC, 2024).

Hubbard, S. J. NACCESS-Computer Program (UK: University College London, 1993).

Holm, S. A simple sequentially rejective multiple test procedure. Scand. J. Stat. 6, 65–70 (1979).

Google Scholar

Jiang, Y., Lee, M.-L. T., He, X., Rosner, B. & Yan, J. Wilcoxon rank-based tests for clustered data with R package clusrank. J. Stat. Softw. 96, 1–26 (2020).

Article
 
CAS
 
Google Scholar

Gregg, M., Datta, S. & Lorenz, D. Variance estimation in tests of clustered categorical data with informative cluster size. Stat. Methods Med. Res. 29, 3396–3408 (2020).

Article
 
PubMed
 
PubMed Central
 
Google Scholar

Gregg, M., Datta, S. & Lorenz, D. htestClust: a package for marginal inference of clustered data under informative cluster size. R. J. 14, 54–66 (2022).

Article
 
Google Scholar

LemaÃŽtre, G., Nogueira, F. & Aridas, C. K. Imbalanced-learn: a Python toolbox to tackle the curse of imbalanced datasets in machine learning. J. Mach. Learn. Res. 18, 1–5 (2017).

Google Scholar

Huang, J. et al. CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nat. Methods 14, 71–73 (2017).

Article
 
CAS
 
PubMed
 
Google Scholar

Humphrey, W., Dalke, A. & Schulten, K. VMD: visual molecular dynamics. J. Mol. Graph. 14, 33–38 (1996).

Article
 
CAS
 
PubMed
 
Google Scholar

Phillips, J. C. et al. Scalable molecular dynamics with NAMD. J. Comput. Chem. 26, 1781–1802 (2005).

Article
 
CAS
 
PubMed
 
PubMed Central
 
Google Scholar

Roe, D. R. & Cheatham, T. E. III PTRAJ and CPPTRAJ: software for processing and analysis of molecular dynamics trajectory data. J. Chem. Theory Comput. 9, 3084–3095 (2013).

Article
 
CAS
 
PubMed
 
Google Scholar

Download references

Acknowledgements

We thank the St. Jude High-Performance Computing Facility for ongoing support; T. R. Stachowski for sharing code, advice and proofreading; A. E. Seffernick for statistical guidance; and K. A. Laycock for scientific editing of the paper. This work was in part supported by the American Lebanese Syrian Associated Charities (ALSAC) and by NIH grant R35GM142772 (to M.F.). The funders had no role in the study design, data collection and analysis, decision to publish or preparation of the manuscript. The content of this Article is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Author information
Authors and Affiliations

Department of Chemical Biology and Therapeutics, St. Jude Children’s Research Hospital, Memphis, TN, USA

Justin T. Seffernick & Marcus Fischer

Contributions

M.F. and J.T.S. jointly conceived the project, designed the study, interpreted the results and wrote the paper. J.T.S. collected the data, developed the methods and performed the analyses.

Corresponding author

Correspondence to Marcus Fischer.

Ethics declarations
Competing interests

The authors declare no competing interests.

Peer review
Peer review information

Nature Methods thanks Stefano Forli, Henry van den Bedem and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: Arunima Singh, in collaboration with the Nature Methods team.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information
Supplementary Information

Supplementary Figs. 1–7, Table 1 and pseudocode for the ColdBrew algorithm.

Reporting Summary
Supplementary Data 1

Set of 242 temperature-matched (cryo–RT) pairs used to develop and test ColdBrew predictions.

Supplementary Data 2

Groups (162 total) of holo PDBs at various binding sites to test water displaceability based on ColdBrew probability.

Supplementary Data 3

For those with available apo structures (101 groups with at least 1 apo), PDBs in groups corresponding to those in Supplementary Data 2.

Rights and permissions

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

Reprints and permissions

About this article
Cite this article

Seffernick, J.T., Fischer, M. An experimental proxy of water displaceability for ligand discovery. Nat Methods 22, 1476–1485 (2025). https://doi.org/10.1038/s41592-025-02724-0

Download citation

Received
08 November 2024

Accepted
09 May 2025

Published
27 June 2025

Version of record
27 June 2025

Issue date
July 2025

DOI
https://doi.org/10.1038/s41592-025-02724-0

Share this article

Anyone you share the following link with will be able to read this content:

Get shareable link

Provided by the Springer Nature SharedIt content-sharing initiative

Subjects
Computational biophysics
Drug discovery
Machine learning
This article is cited by
Unfreezing structural biology for drug discovery
Timothy R. StachowskiMarcus Fischer

Nature Chemical Biology (2025)

Cancel

Preview · live