We compared our RBFEs for the ligand pairs with literature values. The systems considered here have previously been studied by other alchemical free energy approaches. Ben-Shalom et al. studied the Thrombin, Scytalone Dehydratase and BTK systems using four different methods. They found that two of the four methods, where the trapped water was placed in the binding site using a hybrid Monte Carlo Molecular Dynamics (MC/MD), along with alchemical transformations of the ligands, resulted in more accurate RBFE calculations, compared to methods where water placement was not carried out along with the ligand transformations. Our RBFEs for the Thrombin and the BTK systems (Table ) are in excellent agreement with the RBFEs reported in their study (-0.99 and -0.86 kcal mol -1 for the Thrombin and -3.35 and -3.99 kcal mol -1 for the BTK systems). Our RBFE value between the Thrombin ligand pair is also within the standard deviation from the experimentally reported RBFE of -0.61 kcal mol -1 . Though we found disagreement between our RBFE value between the Scytalone Dehydratase ligand pair (-1.59 ± 0.5 kcal mol -1 ) and the values from Ben-Shalom et al. (-3.24 and -2.59 kcal mol -1 ), our value is in a better agreement with the experimental value of -1.98 kcal mol -1 39 than the values in the other study. The Factor Xa system has previously been studied by Abel et al., using molecular dynamics simulations and a solvent analysis tech-nique based on inhomogeneous solvation theory. The previous study reported RBFE values of -1.73 and -1.95 kcal mol -1 using two different solvent functionals. Our RBFE value of -2.46 ± 0.5 kcal mol -1 between the Factor Xa ligand pair remains in line with the literature values, as well as with the experimental value of -1.98 kcal mol -1 . Another study by Wahl et al. has reported RBFEs for three of the five systems studied here: the Scytalone Dehydratase, the BACE1, and the BTK systems. In their study, they calculated the RBFEs between ligands using Free Energy Perturbation (FEP) calculations and Replica Exchange Solute Tempering (REST) starting from different solvation states of the binding site, which resulted in large hysteresis. The hysteresis was reduced by introducing a Grand Canonical Monte Carlo (GCMC) step before the FEP calculations. Between the C3d-C5d ligand pair of Scytalone Dehydratase (Figure ), they reported RBFEs of -0.91 and -1.33 kcal mol -1 for two different initial solvation conditions, which are both in agreement with the RBFE we calculated. Ross et al. presented an advancement of the FEP method and reported the Scytalone Dehydratase ligand RBFEs of 0.37 and 1.04 kcal mol -1 for two different initial solvation states. Our RBFE value (-1.59 ± 0.5 kcal mol -1 ) disagrees with the values reported by Ross et al., however, unlike the RBFEs reported by Ross et al., our RBFE agrees with the experimental RBFE of -1.98 kcal mol -1 . Our RBFE value of 0.32 ± 0.6 kcal mol -1 between the BACE1 ligand pair disagrees with the values reported by Wahl et al. (-2.15 and -1.78 kcal mol -1 ) and Ross et al. (-1.69 and -2.35 kcal mol -1 ). Neither we nor the previous studies could reproduce the experimental RBFE of -0.61 kcal mol -1 . Between the BTK ligand pair, our RBFE value of -4.08 ± 0.7 kcal mol -1 is in disagreement with the experimental RBFE of -1.98 kcal mol -1 , 44 as well as with the RBFE values reported by Wahl et al. (-1.76 and -3.32 kcal mol -1 ), however it is in agreement with the values reported by Ross et al. Overall, we obtained RBFE values for our systems largely in agreement with the RBFE values reported by previous computational studies that implemented the MC/MD and GCMC methods, as well as with the experimental RBFE values. However, for the BACE1 system, we observe disagreements with the literature values. For the BTK system, even though our RBFE value agrees with a previous computational study, it remains in disagreement with the experimental estimate by around 2 kcal mol -1 . This prompted us to investigate the BACE1 and the BTK systems in greater detail, and our findings are described in the following sections. 3.2.3. Ligand rotation and lack of conformational sampling of the protein affected our RBFE estimate between the BACE1 ligand pair. A sampling challenge we observed in our BACE1 simulations was the flexibility of the ligand. The BACE1 binding site can accommodate a rotation of the pyridine ring of the ligand C4j (Figure ). In our simulations, upon decoupling of the trapped water, the hydrogen bond between the trapped water and the nitrogen of the pyridine ring vanished and this led to rotations of the pyridine ring. The rotations of the pyridine ring caused several of the NES switches for ligand transformation to start from the rotated conformations, which resulted in different degrees of overlap in the forward and reverse NES work distributions of the ligand transformation across our replicates. This resulted in a large error contribution (of 0.4 kcal mol -1 ) of the ligand transformation edge in the thermodynamic cycle (edge K, Figure ). The previous study by Wahl et al. has also reported the rotations of the pyridine ring of the ligand C4j in their simulations.