Based on preliminary simulations, the equilibration time for the vanish stage simulations was increased to 3 ns per window to remove systematic error from slow movement of water into the binding site. However, as shown by Figure , occasionally water failed to enter the binding site during the entire simulation at values of λ where several waters entered the binding site during replicate runs. Therefore, some systematic error likely persisted in some Figure : Waters (shown as spheres) in the binding site at λ = 0.4 for B1 at 6.48 ns. Only binding site waters (as defined as the overlap of two spheres of radius 8 Å centred on the N atom in Pro1A and CG2 in Val106A) are shown. The run 1 (red ligand/ waters) trajectory was superimposed on that from run 3 (blue ligand/ no waters present) by aligning MIF. Surface generated based on the ligands alone to approximately show the binding pocket.The ligand must be pushed to the side of the binding pocket to provide space for the water. Rendered with PyMOL. cases. This appeared to be true for the B1 simulations. Irrespective of whether the force constants were fit to simulation (B1), or set to 10 (B1-10) or 20 (B1-20) kcal mol -1 Å-2 or kcal mol -1 rad -2 , the results were more negative than for B2 and B3. In most cases the differences were significant at 95 % confidence. In addition, these simulations generally showed the greatest uncertainties between replicates, demonstrating greater random error. This was despite the B1 restraints being selected as optimum based on the minimum variance algorithm, and the protein anchor points being based on a stable Asn forming part of a βsheet, highlighting the difficulty of selecting optimal restraints. It seemed unlikely that the offset was due to the restraint of the phenol group, because multiple distance restraint schemes based on this group did not produce such negative free energies of vanishing (see Section 2.4). To ensure that this was not due to flexibility of the Asn side-chain, the calculation was repeated using ligand anchors in the phenol but protein anchors only in the backbone and Cα of this Asn. The result was very similar (-7.88 ± 0.95 kcal mol -1 ), showing that the issue was not side-chain flexibility. Finally, the calculation was repeated using ligand atoms in the phenol but protein anchors in a different residue (Ile64A) -see B1-P. This resulted in a ∆G o Bound of -6.69 ± 0.28 kcal mol -1 which had a much smaller uncertainty and was closer to the results for B2 and B3. The difference between B1 and B1-P was significantly different based on a Student's t-test, although it is acknowledged that the assumption of normally distributed free energy differences is likely to be incorrect. In this system especially, where the dominant source of error appears to be the slow hopping of water between energy minima, the central limit theorem is unlikely to be applicable. Regardless, the B1-P results were substantially less negative than the B1-10 and B1-20 results, suggesting that basing the protein anchors on the Asn for the B1 restraints may have made the simulations more susceptible to water sampling issues. It is possible that the restraints between the phenol and this Asn, which quickly forms a H-bond to water upon its entry to the binding site, created particularly high barriers to the entry of water, giving rise to systematic error. Regardless, the difference observed suggests that the approach of running independent replicate simulations with different restraints, as taken by Alibay et al., is sensible.