HIV-1 protease (hivpr) is a 99 amino acids viral aspartic protease that for role to cleave viral polyproteins to mature structural proteins and viral enzymes . hivpr is a homodimeric enzyme whose active site exhibits perfect twofold symmetry in the absence of bound inhibitor or substrate . The active-site region is covered by two symmetry-related β hairpins (a glycine-rich loop), known as the "flaps". It is part of substrate binding site and plays an important role in the substrate binding through a bridging water molecule, W301 . This water molecule forms two hydrogen bonds with the NH groups of Ile50 and Ile50' in hivpr, and at the same time engages hydrogen bond acceptor groups at the P1' and P2 sites like the carbonyls of the bound inhibitor KNI-272 . This bridging water molecule plays a critical role in the interaction between hivpr and its inhibitor and is in fact observed in nearly all hivpr complexes, except when specifically targeted for its displacement . For example, the position of this water molecule has been exploited to design cyclic urea inhibitors able to bind by taking advantage of its displacement entropy . The GIST analysis on hivpr showed convergence problems with the MD trajectories, requiring the MD simulations to be extended from 30 to 100 ns (Fig. ). The difficulty reaching convergence might be responsible for the overall inferior performance of WaterKit on this system, with relatively lower agreement with MD results with respect to other systems (Fig. , and Fig. ). A number of previous studies used different sampling protocols and methods to estimate desolvation free energy of W301 in complex with inhibitors. ∆G energies reported using MD-based methods range from -15.12 kcal/mol (Li and Lazaridis 122 ), to -3.1 ± 0.6 and -3.2 ± 0.4 kcal/mol (Hamelberg and McCammon 141 ; Lu et al ). In another study, Barillari et al. used multiple crystallographic structures to estimate ∆G energies ranging from -10 ± 0.5 to -7.1 ± 0.5 kcal/mol. WaterKit predicted the position of W301 in the PDB structure 2zye within 0.3 Å from the crystallographic coordinates, with an estimated ∆G of -5.7 kcal/mol (E sw = -3.7, E ww = -3.0, TS trans = -0.6 and TS orient = -0.5 kcal/mol) (Fig. ). Despite its well-defined receptor interaction pattern, MD simulations could not identify a well-defined hydration site, but predicted a network of multiple hydrogen bond-connected water molecules occupying positions in between the experimental coordinates of the conserved water (RMSD from W301 between 1.4 and 2.1 Å, Fig. ) and the atoms from the bound ligand KNI-272. We hypothesized that this result might be due to a possible bias in the crystallographic coordinates used. In absence of the inhibitor, the positions of Ile50 backbone amide nitrogens in the flaps appear to favor dynamic configurations with two water molecules interacting with each other and with each nitrogen independently. This, in turn, could lead to an unbalance between proteinwater and water-water interactions, where protein-water interactions at position W301 are not strong enough to identify a stable energy minimum with a single water coordination. Such bias would be likely reinforced by the protein mobility restrains used in the MD simulations. This hypothesis was tested by repeating MD simulations using a different PDB (1hpx) bound to the same ligand (KNI-272) using an identical protocol. The new MD simulations were more successful in identifying a hydration site closer to W301 position in 2 out of 3 of the replicates, while WaterKit showed results consistent with the previous calculations (Fig. ). Other stable water molecules surrounding bound ligands are considered responsible for stabilizing their binding , either by potentially bridging ligands, or by being displaced by them. The geometric analysis of the water network reconstructed by WaterKit shows that a nearly-tetrahedral water cluster recovers the positions of the terminal methyls of the N-(ter-butyl amide group, suggesting that this group is optimally structured for their displacement. Similarly, another water cluster occupies the space engaged by the thiomethyl group of the ligand, while other more flattened clusters mimic the hydrophobic rings of the ligand. WaterKit also predicted the position of several water molecules in correspondence of key hydrogen bond acceptor and donoramide atoms on the inhibitor used to mimic the peptidic substrate of the enzyme. Of these amide water molecules, only one was detected by MD simulations. Furthermore, both methods identified the positions of several structurally conserved water molecules in the vicinity of the ligand, such as HOH566, HOH608, HOH349, HOH422 (2zye residue numbering), which can be potentially exploited to further increase ligand affinity. Only the first two water molecules were identified by MD simulations, while Wa-terKit recovered all of them. Finally, both WaterKit and MD simulations identify the water molecule bound in between side chains of the catalytic Asp diad, and which is known to be actively involved in the catalytic mechanism of the enzyme . The data on hivpr highlights how relatively minor structural variations, such as those in the catalytic site flap region, can dramatically influence and alter water molecule positions calculated with constrained MD simulations. Conversely, the ensemble sampling used in WaterKit, appears to be less sensitive to the initial conformation used and able to systematically find the position of key structural water molecules, including the crystallographic water 301. Ultimately, while MD is capable Fatty acid binding protein adipocyte (fabp4). The Fatty acid binding protein adipocyte protein (fabp4) is composed of 132 amino acids arranged in 10-stranded anti-parallel β-sheet forming a β-barrel. Within this β-barrel is a large cavity with a volume of approximately 950 Å, formed by the β-barrel and two α helices that are acting as a cavity "lid" . Despite the hydrophobic nature of the endogenous fatty acids ligands the cavity is lined with both hydrophobic and polar amino acids. It was observed that several water molecules are conserved regardless of the presence or absence of ligand . Those bound water molecules serve as side chain extensions for hydrophilic residues, forming a surface against which ligands are laying 144 . It was also hypothesizing that those internalized polar residues, and herein presence of conserved water molecules, are important to avoid the structure to collapse into an hydrophobic core, eliminating the cavity 144 . Visual inspection of the crystallographic structure (PDB: 2nnq) shows only two possible accesses for water molecules to exchange between the cavity and the bulk solvent (near residues Thr60 and Arg126), making fabp4 a good candidate for comparing water sampling methods in deeply buried cavity. Results shows little to no correlation between MD simulations and WaterKit. For a more in-depth analysis and due to convergence issues, MD simulation triplicates were extended from 30 to 200 ns. While WaterKit is not dependent on the solvent accessibility of a given site, the difficulty of exchange with the bulk can have a dramatic effect on MD simulations, unless specific methods are applied to facilitate bulk/cavity exchanges . In order to further characterize the discrepancies observed between WaterKit and the MD simulations, we analyzed the water density in fabp4 cavity. The water density was estimated based on the volume of the cavity using PyVol 146 and the number of water molecules using a sphere of 8.0 Å radius centered on the ligand. The reference bulk density used was 0.0334 molecule per Å . The average water density in MD simulations is about 10 % lower than the bulk density, with 0.90 ± 0.07 times the average bulk density (Fig. ; minimum and maximum density of 0.57 and 1.28 bulk density, respectively). WaterKit results appear to be more stable across the packing trajectory, with a calculated water density average of about 1.03 ± 0.07 bulk density, (Fig. ; minimum and maximum density of 0.70 and also 1.28 bulk density, respectively). The average density calculated with the MD simulations is similar to what was reported to the buried cavity of cytochromes P450, which was found to have an average density up to ≈ 20 % lower than the bulk solvent 147 . However, significant transient deviations from the mean can be observed during MD simulations of fabp4, in which the cavity goes through a temporary de-wetting phase (water density 40 % lower than the bulk density, around 100 ns in at least one of the triplicates) and a high-density phase (water density 1.1 times the bulk density, around 80 ns in another replicate). Previous studies showed that drying transition could naturally occur in confined regions . However, due to the restraints on heavy atoms, the MD simulations show more difficulty to compensate for transient density fluctuations, and therefore reaching convergence, likely generate artifacts in the solvation energy. Conversely, the anal-ysis suggests that WaterKit generate more consistent water configurations with reduced density fluctuations. The packing method prevents also de-wetting phenomena since water positions are sampled independently for each frame, which could also helps reaching convergence more rapidly. The difference in average density between the two methods translates approximately to 4 additional water molecules in the enclosed cavity with respect to MD results. These factors ultimately can be responsible of the reduced overall agreement with MD simulations.