==== Front J Chem Theory Comput J Chem Theory Comput ct jctcce Journal of Chemical Theory and Computation 1549-9618 1549-9626 American Chemical Society 37285579 10.1021/acs.jctc.3c00139 Article Comparison of Receptor–Ligand Restraint Schemes for Alchemical Absolute Binding Free Energy Calculations https://orcid.org/0000-0003-0474-5475 Clark Finlay † https://orcid.org/0000-0002-4531-4375 Robb Graeme ‡ Cole Daniel J. ¶ https://orcid.org/0000-0003-0360-1760 Michel Julien *† † EaStCHEM School of Chemistry, University of Edinburgh, David Brewster Road, Edinburgh EH9 3FJ, United Kingdom ‡ Oncology R&D, AstraZeneca, Cambridge CB4 0WG, United Kingdom ¶ School of Natural and Environmental Sciences, Newcastle University, Newcastle upon Tyne NE1 7RU, United Kingdom * E-mail: julien.michel@ed.ac.uk. 07 06 2023 27 06 2023 19 12 36863704 02 02 2023 © 2023 The Authors. Published by American Chemical Society 2023 The Authors https://creativecommons.org/licenses/by/4.0/ Permits the broadest form of re-use including for commercial purposes, provided that author attribution and integrity are maintained (https://creativecommons.org/licenses/by/4.0/). Alchemical absolute binding free energy calculations are of increasing interest in drug discovery. These calculations require restraints between the receptor and ligand to restrict their relative positions and, optionally, orientations. Boresch restraints are commonly used, but they must be carefully selected in order to sufficiently restrain the ligand and to avoid inherent instabilities. Applying multiple distance restraints between anchor points in the receptor and ligand provides an alternative framework without inherent instabilities which may provide convergence benefits by more strongly restricting the relative movements of the receptor and ligand. However, there is no simple method to calculate the free energy of releasing these restraints due to the coupling of the internal and external degrees of freedom of the receptor and ligand. Here, a method to rigorously calculate free energies of binding with multiple distance restraints by imposing intramolecular restraints on the anchor points is proposed. Absolute binding free energies for the human macrophage migration inhibitory factor/MIF180, system obtained using a variety of Boresch restraints and rigorous and nonrigorous implementations of multiple distance restraints are compared. It is shown that several multiple distance restraint schemes produce estimates in good agreement with Boresch restraints. In contrast, calculations without orientational restraints produce erroneously favorable free energies of binding by up to approximately 4 kcal mol–1. These approaches offer new options for the deployment of alchemical absolute binding free energy calculations. AstraZeneca 10.13039/100004325 NA Engineering and Physical Sciences Research Council 10.13039/501100000266 NA UK Research and Innovation 10.13039/100014013 MR/T019654/1 document-id-old-9ct3c00139 document-id-new-14ct3c00139 ccc-price ==== Body pmc1 Introduction The in silico prediction of protein–ligand binding affinities is an important problem in drug discovery. The ability to rapidly and accurately calculate affinities for arbitrary protein–ligand systems would allow the efficient prioritization of compounds for synthesis and testing, accelerating the hit-to-lead and lead optimization stages of drug discovery.1 Recent improvements in computing power and automation have brought this vision closer to realization.2−6 In particular, alchemical methods are ideally suited for application during the hit-to-lead and lead optimization stages of drug discovery,1,7 as well as in the later stages of virtual screening.8 Along with path-based methods,9 alchemical simulations form a class of exact (in the limit of complete sampling and a perfect description of the potential energy) methods based on molecular dynamics or Monte Carlo sampling which provide binding affinity predictions of greater accuracy than alternatives.10 Modern computing resources now support routine use of alchemical relative binding free energy (RBFE) calculations to add value during drug discovery campaigns.1 RBFE calculations avoid the computationally intractable challenge of converging unbiased simulations of ligand binding and unbinding by instead gradually interconverting two structurally similar ligands. Interconversion proceeds through unphysical “alchemical” intermediates and is done in both the bound and unbound states.11 Based on a thermodynamic cycle, the free energy differences for each step are summed to yield the difference in the free energy of binding between the ligands. RBFE calculations are used routinely, and protocols for their robust deployment have been researched in detail.12−15 However, as a result of the requirement for a common ligand core, binding pose, and binding site, the following valuable problems typically lie outside the scope of RBFE calculations:8,16(1) Calculating the RBFEs of structurally dissimilar ligands to a common target. (2) Calculating RBFEs of the same ligand to the same protein with different binding poses. (3) Calculating the RBFEs of the same ligand to different targets. (4) Calculating the absolute binding free energy of a given ligand to a given target. Alchemical absolute binding free energy (ABFE) calculations escape these limitations by following a more general thermodynamic cycle in which the ligand’s intermolecular interactions are completely turned off.8 In principle, these calculations can be used to calculate the binding free energies of structurally diverse molecules to varied targets, making them attractive for drug discovery. However, the alchemical ABFE framework presents challenges not encountered during an RBFE calculation. Restraints must be applied between the protein and ligand to avoid convergence issues such as those associated with the ligand “wandering” out of the binding site as its intermolecular interactions are removed.17 Restraints may also be required to avoid errors in the calculated binding free energies when the bound state is implicitly defined to include configurations where the ligand is anywhere in the entire simulation box relative to the receptor, as is the case when restraints are not used. However, these errors only affect weak binders (Section S1). In general, it is nontrivial to select the optimum receptor–ligand restraints. Furthermore, ABFE calculations can be challenging to converge and therefore computationally costly because the ligand is completely removed.18,19 As a result, application studies still combine RBFE and ABFE, with ABFE applied more successfully to low molecular-weight compounds.20 Thus, there are barriers to the routine application of ABFE calculations. The performance and accessibility of ABFE calculations would be improved if it was trivial to select receptor–ligand restraints which resulted in stable simulations and produced optimal convergence. While progress has been made in this direction with tools for automated or partially automated restraint selection,3−5,18,21 there is still no restraint type or selection method which completely solves this issue. Receptor–ligand restraints of a variety of forms have been proposed. Following early work utilizing restraints on a single ligand atom,22,23 the first theoretically rigorous approach involving restraints on all of the external degrees of freedom (DoF) of the ligand was the Body Restraint Algorithm of Hermans and Wang.24 Later, the Virtual Bond Algorithm (VBA) of Boresch et al. was introduced,17 which involves restraining one distance, two bond angles, and three dihedral angles between six anchor points defined by the receptor and ligand. This provided a more convenient method to restrain the relative external DoF of the receptor and ligand, along with a simple analytical correction for releasing the restraints. The VBA has found widespread use and is often referred to as “Boresch restraints”. However, despite their popularity, Boresch restraints suffer from a number of limitations and must be carefully applied to avoid numerical instabilities and sampling issues.25 For instance, if the anchor points are tied to the positions of highly flexible portions of the ligand or protein which do not strongly interact, then the restraints will be unable to maintain a binding pose similar to the restrained and interacting system, potentially leading to slow convergence of free energy estimates. Since only six relative external degrees of freedom can be restrained within this framework, there are limits on the extent to which ligand motions can be restricted. Thus, additional restraints on the intramolecular degrees of freedom of the ligand may be required to improve convergence for flexible ligands.26 Furthermore, if the restraints are poorly chosen, small changes in the Cartesian coordinates of the anchor points can result in large jumps in the six DoF defined in the VBA framework, resulting in the application of large forces, which can cause simulations to crash. This frequently occurs when sets of three contiguous anchor points approach collinearity,3 which can result in the application of large forces through the dihedral restraints. Alternative restraint schemes have recently been proposed to address these issues. Fu et al. proposed a method in which the restrained six external DoF are derived by finding the optimal rotation of the ligand which minimizes its root-mean-square deviation (RMSD) with respect to a protein–ligand complex reference structure (after correcting for rotation and translation of the protein).27 By moving away from the six anchor points of the VBA, this was intended to simplify the selection of stable and efficient restraints. The “distance to bound configuration” (DBC) restraint is also intended to simplify restraint selection and to minimize the variance of free energy estimates for removing the ligand’s intermolecular interactions.28−30 This is achieved by directly restraining the RMSD of a subset of the ligand coordinates within the frame of reference of the binding site in order to optimally restrict the accessible configurational volume as the ligand intermolecular interactions are removed. However, because this scheme couples the internal and external degrees of freedom of the protein and ligand, there is no simple way to calculate the free energy of releasing the noninteracting ligand to the standard state. This necessitates a final stage to release the DBC restraints to a single harmonic restraint, for which the free energy of release is simple to calculate. Another alternative to Boresch restraints is to restrain the distance between multiple receptor–ligand atom pairs. These restraints offer several advantages; for example, they can be intuitively selected to match native receptor–ligand interactions such as hydrogen bonds, thus closely mimicking the interacting state. This may accelerate convergence by tightly restricting ligand motion while intermolecular interactions are removed. Furthermore, these restraints do not suffer from the numerical instabilities inherent to the Boresch restraints scheme. Indeed, multiple distance restraints were used in an early study of the binding of biotin to streptavidin.31 However, the naive application of multiple distance restraints is theoretically incorrect, because they introduce coupling between the internal and external degrees of freedom of the protein and ligand, preventing the rigorous calculation of the free energy of releasing the noninteracting ligand.32 Despite this, a recent implementation has been described which relied on the assumption that the restraints were sufficiently weak that such coupling was negligible and that the free energy of turning on the restraints was close to zero.33 However, this was not verified, and the scheme has not been systematically compared to Boresch restraints. To address this, this study compares the absolute binding free energies obtained using Boresch restraints and different implementations of multiple distance restraints for a single ligand (MIF-180) binding to a single protein (human macrophage migration inhibitory factor, or MIF). This was suggested as a good model system by Qian et al.34 because MIF is a pharmaceutically relevant protein, but of moderate size (342 residues), and no major conformational changes occur in the protein upon ligand binding (Figure 1). In this work, the standard binding free energy of MIF-180 was calculated using multiple sets of Boresch restraint parameters. The results were compared to those produced by nonrigorous implementations of multiple distance restraints similar to that of Mendoza-Martinez et al.33 and two rigorous multiple distance restraint schemes: one inspired by Salari et al.28 and one newly developed. Figure 1 (A, B) MIF with MIF180 bound, rendered with PyMOL.35 (C) Hydrogen bonding interactions between MIF and MIF180 in the tautomerase active site. Redrawn from Qian et al.34 2 Theory 2.1 Alchemical Absolute Binding Free Energy Calculations ABFEs can be computed using an alchemical cycle (Figure 2). The standard free energy of binding is calculated by adding up the terms around the cycle:1 where “Free” and “Bound” indicate that the ligand is in solution or in the receptor binding site. “Discharge” means removal of ligand Coulombic interactions. “Vanish” indicates removal of the ligand Lennard-Jones (LJ) interactions, and “Restrain” means introduction of intermolecular restraints between the receptor and ligand. ΔGReleaseo is the correction for releasing the receptor–ligand restraints when the ligand has no intermolecular interactions, and ΔGSym.Corr. accounts for symmetries broken by the introduction of restraints.36 The bound leg refers to all calculations where the ligand is in the binding site and includes the symmetry correction, and the free leg describes all calculations where the ligand is in solution. Figure 2 An alchemical thermodynamic cycle for the calculation of ABFEs. The red dashed lines indicate protein–ligand restraints. It is generally computationally intractable to obtain ΔGBindo by direct simulation, but the alchemical cycle allows ΔGBindo to be obtained through a series of states which are less challenging to sample at equilibrium. ΔGReleaseo is calculated without simulation. When only the intermolecular components of the ligand nonbonded interactions (Coulombic or Lennard-Jones (LJ)) are removed, this is often termed “decoupling”, while “annihilation” may refer to removal of both the inter- and intramolecular components.37 However, the terminology of Gilson et al. is used here:38 “decoupling” denotes removal of of the ligand intermolecular interactions while enforcing receptor–ligand restraints, irrespective of how the intramolecular interactions are treated. 2.2 Receptor–Ligand Restraints The free energy of releasing the restraint on the decoupled ligand is given by the ratio of configurational integrals2 where kB is the Boltzmann constant, T is the temperature, and ZState3 and ZState4 are the configurational integrals for States 3 and 4, as defined in Figure 2. However, this result is dependent on the size of the water box in State 3, VBox. A more useful quantity is the standard free energy of binding3 which is independent of VBox. Vo = 1660 Å3 is the standard state volume. This ratio could be evaluated using a simulation in which the restrained decoupled ligand is completely released into the simulation box, but this would be slow to converge. Instead, this ratio must be simplified so that it can be evaluated without simulation. It can be shown (Section S2) that eq 3 can be written4 where Wr(xExt) is the potential of mean force (PMF) of the receptor–ligand restraint energy with respect to the six relative receptor–ligand external DoF, xExt. The form of the Jacobian determinant, |J|, depends on the coordinate transformation used to extract the relative external degrees of freedom from the internal degrees of freedom of the complex. ΔGPreorg. accounts for straining of the receptor and decoupled ligand when Wr(xExt) is at its minimum, while ΔGDistort. accounts for further distortion of the receptor and ligand when Wr(xExt) is not at its minimum. For an arbitrary set of receptor–ligand restraints, there is no straightforward way to evaluate this expression. The standard solution to obtain ΔGReleaseo is to select a set of receptor–ligand restraints for which the restraint energy, Ur(xExt), is a function of only the receptor–ligand relative external degrees of freedom. In this case, ΔGPreorg. = ΔGDistort. = 0, because the relative external coordinates of the decoupled ligand and receptor are uncorrelated with the internal DoF and Wr(xExt) = Ur(xExt). Restraints of this form are described as not coupling the internal and external degrees of freedom of the receptor and ligand. Intuitively, such restraints do no “squeeze” or “stretch” the receptor or decoupled ligand. The free energy of releasing these restraints is5 which can be integrated directly. Thus, the ideal receptor–ligand restraints would not couple the internal and external degrees of freedom of the receptor and ligand, so that eq 5 is valid. The ideal restraints would also ensure optimal convergence of the bound stages. This might be achieved by mimicking the native receptor–ligand interactions as closely as possible,16 thus minimally perturbing the fully interacting complex while maximally restricting the accessible configurational volume during decoupling.29 However, the extent to which this can be achieved is limited when restraining only six DoF. Hence, it may be desirable to use restraints which do couple the these internal and relative external degrees of freedom in order to accelerate convergence of the decoupling calculations. In this case, the restraints which do couple the these degrees of freedom should be released through simulation to those which do not (and the associated free energy change accounted for), or some degree of error must be tolerated in the calculation of ΔGReleaseo. Finally, the ideal restraints would lack instabilities, and would be simple to select in an easily automatable manner.4,27,39 2.3 Boresch Restraints Boresch restraints (Figure 3) only affect the six relative external degrees of freedom of the ligand with respect to the receptor, and do not couple the receptor and ligand internal and external degrees of freedom.17 As a result, eq 5 can be easily evaluated by numerical integration of the expression6 where u(x) is the restraining potential applied to the degree of freedom x. These degrees of freedom are defined in Figure 3. The second term in eq 6 integrates over the position of anchor atom A with respect to the receptor coordinates, and does not depend on anchor atoms B and C. The third term integrates over the orientation of the ligand with respect to the receptor. Hence, Boresch restraints can be used to restrain only the position of the anchor point A with respect to the receptor, by setting u(θB), u(ϕB), and u(ϕC) to 0. Figure 3 General form of Boresch restraints.17 Three anchor points (a, b, and c) are selected based on the receptor (green) coordinates and three (A, B, and C) are selected based on the ligand (blue) coordinates. The restrained six external degrees of freedom are defined as one distance (r), two bond angles (θA and θB), and three dihedral angles (ϕA, ϕB, and ϕC). In this diagram, the contiguous anchor points b, a, and A are close to collinear and therefore these restraints would be expected to be unstable. Although it is common to use harmonic restraining potentials and to select the anchor points as atomic positions, these are not constraints of the framework. For example, periodic dihedral restraints can be used,5 and anchor points may be derived from multiple atomic positions or centers of mass.34 Harmonic restraints are used in this work. In this case, eq 6 can be evaluated analytically as7 where K denotes a force constant and 0 denotes an equilibrium value.17 This assumes that r, sin θA, and sin θB can be taken out of the integrals in eq 6 and replaced by their equilibrium values. As mentioned previously, when the anchor points are arranged so that large changes in the six DoF defined in the VBA framework result from small changes to the Cartesian coordinates of atoms, this can result in large forces and simulation crashes. To avoid such instabilities, anchor points must be carefully selected to avoid the collinearity of any three contiguous anchor points. Algorithms have been proposed for the selection of Boresch restraints.3,16,18,40−42 At a minimum, these aim to select stable restraints based on the geometry of the complex, while more sophisticated methods aim to enhance convergence by directly (e.g., based on H-bonds) or indirectly (based on minimum total variance of the distance, angles, and dihedrals) mimicking strong receptor–ligand interactions based on a short unrestrained simulation. However, there is no obviously superior method which has been shown to guarantee selection of numerically stable restraints with optimal convergence properties. 2.4 Multiple Distance Restraints Restraints schemes based on multiple distance restraints do not suffer from the inherent instabilities of Boresch restraints, and allow the ligand to be restrained to a greater extent (Figure 4). In this work, harmonic or flat-bottomed distance restraints were used. Using the former, Ur is given by8 and with the latter the total restraint energy is given by9 where rn is the distance between atoms in a restrained pair, rn,0 is the equilibrium distance, rn,fb is the flat-bottomed radius, and Kn is the force constant for restrained pair n. N is the total number of restrained pairs. Figure 4 A ligand restrained using multiple distance restraints. Four distance restraints (r1–r4) are applied between pairs of atoms in the ligand (blue) and the receptor (green). Applying distance restraints between several receptor–ligand atom pairs can allow greater restriction of ligand movement than a typical six DoF restraint scheme. However, multiple distance restraints can couple the internal and external degrees of freedom of the receptor and ligand, as shown by the distortion of the ligand, preventing rigorous calculation of the free energy of releasing the decoupled ligand. The main flaw of multiple distance restraints is that, applied naively, they may couple the internal and external degrees of freedom of the receptor and ligand, preventing the simplification of eq 4 to 5 and hence the calculation of ΔGReleaseo. Here, we investigate three approaches to circumvent this difficulty (Figure 5). The first is to implement multiple distance restraints in such a way that the error introduced is negligible (the “naive” approach). This is the basis for recent implementations in the Michel group,33,43,44 in which sets of relatively permissive flat-bottomed restraints were used. An approximate value for ΔGReleaseo is calculated by numerical integration of10 where xBox, yBox, and zBox are the side lengths of the simulation box, and x, y, z, ψ, θ, and ϕ are the Cartesian coordinates of the center of mass and the Euler angles of the ligand in the frame of reference of the receptor. This is evaluated by taking the average intramolecular coordinates of the anchor points from a simulation of State 4. The receptor and ligand are then assumed to be rigid, allowing Ur(x, y, z, ψ, θ, ϕ) to be calculated by translating and rotating the ligand anchor points with respect to the receptor anchor points and evaluating Ur, Rigid(r1, ..., rN) = Ur, Rigid(x, y, z, ψ, θ, ϕ) at each point, where “Rigid” shows that the anchor points have been fixed to their average intramolecular positions. Figure 5 Multiple distance restraints schemes. From top to bottom, all schemes start in State 3 and progress to State 7 (Figure 2). For the example shown, the “naive” (nonrigorous) multiple distance restraint scheme may suffer from large systematic error due to distortion of the flexible ligand, which is not accounted for in the −ΔGReleaseo stage. This is avoided in the other schemes by either applying strong intramolecular restraints, or releasing all but one distance restraint. This assumes that ΔGPreorg. = 0 (eq 4), and produces a bias toward more negative free energies of binding, because ΔGPreorg. ≥ 0. The error will be substantial when the average intramolecular positions of the anchor points for the decoupled complex are very different to those in the free ligand and receptor, or the restraints substantially restrict the conformational freedom of the receptor or ligand, thus enforcing substantial “preorganization”. The magnitude of the error would be expected to increase with the number and strength of restraints, and decreasing volume of the flat-bottom region. This is likely to be particularly problematic for systems where there are substantial changes in the conformations of the ligand and the binding site upon binding. In addition, the use of average positions assumes that ΔGDistort. = 0 and Ur, Rigid(x, y, z, ψ, θ, ϕ) = Wr(x, y, z, ψ, θ, ϕ). This produces a slight bias toward more positive free energies of binding, because a flexible system will distort to minimize the sum Wr(x, y, z, ψ, θ, ϕ) + ΔGDistort., Point(x, y, z, ψ, θ, ϕ) ≤ Ur, Rigid(x, y, z, ψ, θ, ϕ) for a given relative position and orientation. ΔGDistort., Point(x, y, z, ψ, θ, ϕ) is related to ΔGDistort. by integration over the six external DoF. Overall, neglecting flexibility is expected to result in an erroneously negative ΔGReleaseo due to neglect of ΔGPreorg.. Obtaining an accurate free energy of binding under these approximations depends upon the restraints being sufficiently permissive and there being no large changes in the average intramolecular coordinates of the anchor points during binding. Here, we investigate an alternative scheme to eliminate these sources of error: applying intramolecular restraints to rigidify the anchor points. The free energy of applying and releasing these in the free and bound stages can be explicitly calculated with simulation. The restraints should be sufficiently strong to “preorganize” the system, and substantially stronger than the intermolecular restraints. This guarantees that ΔGPreorg. = ΔGDistort. = 0 and that Wr(x, y, z, ψ, θ, ϕ) = Ur, Rigid(x, y, z, ψ, θ, ϕ), rendering eq 10 rigorous. An alternative rigorous approach to the implementation of multiple distance restraints is to evaluate the free energy change, ΔGTo Dist. Rest., for releasing them to a single harmonic or flat-bottomed restraint after decoupling. Once only a single distance restraint is active, the free energy of release can be calculated exactly, as is done for the DBC restraint.28,29 The free energy of releasing a single distance restraint can be calculated using11 where r is the distance between two anchor atoms, and Ur(r) takes the form of eq 8 or 9. 3 Methods 3.1 System Preparation The MIF/MIF180 systems were set up approximately following the methodology of Qian et al.34 The structures of MIF and MIF180 were obtained from the crystal structures of complexes of MIF with both MIF180 (relatively poor resolution −2.6 Å, PDB ID 4WR8) and the structurally similar MIF190 (relatively high resolution −1.8 Å, PDB ID 4WRB).45 Five extra copies of the biological assembly were removed from 4WR8 by deleting all chains from D onward. The higher resolution structure was superimposed on the lower resolution structure by aligning chain B from 4WR8 with chain A from 4WRB using PyMOL.35 Crystallographic waters (from 4WRB) were retained, and all other nonprotein and nonligand atoms were discarded. For atoms with alternative locations, the location with the greatest occupancy was selected. Thirteen missing atoms were added using pdb4amber.46 H++ (version 3.2) and PROPKA (version 3.4.0) were used to suggest the protonation state of all residues.47−49 Consistent with the literature for the apo protein,50 a large pKa shift was predicted in both cases for the N-terminal proline residue which is present in the binding site, which produces the neutral form at a pH of 7. The suggested protonation sites for all histidines were taken from H++. Hydrogens were added to MIF using H++ and a hydrogen ion was removed from the N of each of the protonated terminal prolines so as to maintain the hydrogen bond to Tyr36. Based on the results of Qian et al., the neutrality of the terminal proline was maintained in the complex. There are no parameters available for neutral N-terminal proline in the AMBER ff14SB force field.51 Therefore, antechamber 21.0 was used to parametrize the neutral proline with an NME-capped C-terminus, using AM1-BCC partial charges and AMBER atom types.52,53 The remainder of the MIF protein was parametrized using the AMBER ff14SB force field. Hydrogens were added to MIF180 using Open Babel (version 3.0.0), and BioSimSpace (version 2020.1.0) was used to parametrize the ligand in the syn conformation (Figure S2) with the GAFF2.11 force field and AM1-BCC charges using antechamber 21.0.54,55 3.1.1 Absolute Binding Free Energy Calculations Alchemical ABFE calculations were performed using the double decoupling method,38 according to the cycle shown in Figure 2 and the multiple distance restraint schemes shown in Figure 5:12 13 14 where ΔGRest. Scheme collects any additional terms which are specific to the restraints scheme used. The relevant restraints scheme for each form of ΔGRest. Scheme is shown on the right of eq 14: “Boresch” includes both Boresch restraints and the “non-rigorous” implementation of multiple distance restraints shown in Figure 5, “intramol. rigid.” refers to multiple distance restraints with intramolecular rigidification, and “release” denotes multiple distance restraints with release to a single distance restraint. ΔGRigid. Recept. and ΔGRigid. Lig. are the free energy changes for intramolecular rigidification of receptor and ligand anchor points, −ΔGRigid. Complex accounts for the release of intermolecular restraints in the receptor-interacting ligand complex, and ΔGTo Dist. Rest. is the free energy of releasing several distance restraints to a single distance restraint. For the intramolecular rigidification scheme, the intramolecular restraints are applied in State 3, and released in State 7. ΔGSym. Corr. is at least – kBT ln 3 for this system, because MIF is a homotrimer with three equivalent binding sites, and all restraints restrict the ligand to a single binding site. ΔGSym. Corr. may be larger if the restraints also break symmetries arising from the structure of the ligand (Section S6). ΔGFree and ΔGBoundo are the overall free and bound leg contributions to ΔGBindo, where ΔGBoundo includes all terms other than ΔGFree, Discharge and ΔGFree, Vanish. No correction is included to account for the truncation of the tails of the LJ potentials because initial simulations yielded negligible corrections (≈ 0.1 kcal mol–1).56 Protein–ligand restraints were introduced, and charges and LJ interactions were incrementally removed by scaling the coupling parameter, λ, from 0 to 1. For calculations with Boresch restraints and multiple distance restraints without intramolecular rigidification or release to a single distance restraint, the force constants of the protein–ligand restraints and the magnitude of the charges were scaled linearly with λ. The soft-core potential implemented in Sire (with a LJ soft-core parameter set to 2.0), which is based on the potentials of Zacharias et al.57 and Michel et al.,58 was used to scale the LJ interactions. Eight evenly spaced windows and 18 nonevenly spaced λ windows (0.000, 0.028, 0.056, 0.111, 0.167, 0.222, 0.278, 0.333, 0.389, 0.444, 0.500, 0.556, 0.611, 0.667, 0.722, 0.778, 0.889, and 1.000) were used for the free discharging and vanishing stages, respectively. Six nonevenly spaced windows (0.000, 0.125, 0.250, 0.375, 0.500, and 1.000), eight evenly spaced windows, and 36 nonevenly spaced windows (31 evenly spaced windows from λ = 0 to 0.750, then five evenly spaced windows from λ = 0.750 to 1.000) were used for the bound restraining, discharging, and vanishing stages, respectively. The window spacings were selected to yield sufficient overlap without excessive numbers of windows based on initial test simulations. The intramolecular components of both the Coulombic and LJ interactions between ligand atoms were completely removed. For calculations where multiple distance restraints were released to a single distance restraint, ΔGReleaseo was calculated by numerical integration of eq 11, and the force constants were scaled with λ5 over 21 evenly spaced λ windows.29 The same protocol was used to introduce and remove all restraints in the multiple distance restraints simulations with intramolecular restraints. Restraints were selected to optimally mimic native protein–ligand interactions by postprocessing a 6 ns simulation of the fully interacting complex.23,26 From the first frame, all heavy atoms in the protein within 10 Å of the ligand, and all heavy atoms in the ligand were selected. To avoid anchor points with poor correlation in position, the distances between all possible protein–ligand atom pairs from this selection were tracked over the trajectory, and the 200 pairs with the lowest standard deviation were selected. For multiple distance restraints, only the lowest variance pair for any anchor point was retained, provided that neither of the anchor points had already been selected for use in another restraint. For Boresch restraints, all pairs were taken as candidate anchor points a and A (Figure 3). For each pair, adjacent heavy atoms were selected to complete the sets of Boresch anchor points. These sets were ordered by increasing total variance of the Boresch DoF, as done by Alibay,40,41 and sets of anchor points were discarded if the average values of θA or θB were below 30 or above 150 degrees. The equilibrium values for all restraints were taken to be their average values during the unrestrained simulation. Force constants were selected so that in the decoupled state, the harmonic restraints would generate the same distributions as observed in the coupled state.23 Gaussian distributions in the coupled state were assumed and the variances of the Boresch DoF were used to calculate the force constants (Section S7). For the Boresch restraints, ΔGReleaseo was calculated by numerical integration of eq 6. This was used in preference to the analytical correction to avoid potential errors introduced by the approximations required to derive eq 7. For multiple distance restraints, numerical integration of eq 10 was performed using the “standardstatecorrection” utility available within Sire,59 using all frames of the trajectory, a buffer of 5 Å, a translational volume element of 0.25 Å, and 30 orientations per [0, 2π] Euler angle interval. 3.1.2 Molecular Dynamics Protocol Solvation and equilibration were performed using BioSimSpace.55 The protein–ligand complex was placed in a periodic cube of side 84 Å (determined by the longest edge of the axis-aligned bounding box plus 15 Å, padding on each side) and solvated with TIP3P water molecules.60 Then, 150 mM NaCl was added. The system was energy minimized using PMEMD (50,000 steps).46 Equilibration in the NVT ensemble was performed using PMEMD (5 ps with all nonsolvent atoms restrained and heating from 0 to 298 K, followed by 50 ps with only backbone atoms restrained, then 50 ps with no restraints), followed by equilibration in the NPT ensemble at 1 atm and 298 K using PMEMD.CUDA (400 ps with all nonsolvent heavy atoms restrained, followed by 2 ns with no restraints). The free ligand was solvated with TIP3P water and 150 mM NaCl in a periodic box of side length 40 Å. Then, 50,000 steps of minimization were performed using PMEMD. Equilibration in the NVT ensemble was performed using PMEMD (5 ps with all nonsolvent atoms restrained and heating from 0 to 298 K, followed by 50 ps with no restraints), followed by equilibration in the NPT ensemble with PMEMD.CUDA (1 atm and 298 K with restraints on nonsolvent heavy atoms for 200 ps followed by 2 ns with no restraints). A Langevin thermostat and Berendsen barostat were used for the relevant equilibration steps.61 All alchemical simulations were performed using the software SOMD,62 available within Sire (version 2022.2.0).59 SOMD was modified to allow the use of Boresch restraints, the scaling of restraints with λ5, and the simultaneous use of different restraints. The code implementing Boresch restraints has been integrated into the main branch of Sire. An Andersen thermostat (collision frequency 10 ps–1) and Monte Carlo barostat (25 time steps between isotropic box scaling attempts) were used to maintain a temperature and pressure of 298 K and 1 atm.63,64 A time step of 4 fs was used in combination with the leapfrog Verlet integrator and hydrogen mass repartitioning (using a repartitioning factor of 4).65 All bond lengths were constrained. The reaction field method was used with a dielectric constant of 78.3,66 and a cutoff of 12 Å was used for all nonbonded interactions. Energy minimization was performed prior to each simulation with a maximum of 1000 iterations. The bound stage vanish λ windows were run for 8 ns, and all others for 6 ns. Free energy differences for each stage were estimated using the Multistate Bennett-Acceptance Ratio (MBAR) for the final 5 ns of all simulations.67,68 Coordinates were saved every 20 ps. All simulations were repeated five times with independent starting velocities. Because the MBAR uncertainties estimated from single runs were small compared to the variation between repeat runs, errors are reported as 95% confidence intervals based on the deviation between independent replicates, assuming Gaussian distributions and using t-values for 4 degrees of freedom. Student’s t-test was used to assess evidence for a significant difference at 95% confidence. 4 Results and Discussion 4.1 Boresch Restraints 4.1.1 Restraint Selection From an initial set of restraining simulations, two binding poses were identified (Figure 6) which interconverted slowly on the time scale of the simulations (6 ns). To allow comparison of different restraints for a single binding pose, all restraints were fit to binding pose A, other than a single set of Boresch restraints which was fit to pose B. The calculation for pose B was carried out to allow comparison to the experimental free energy of binding. Figure 6 Alternative binding poses A (panel A) and B (panel B). Interconversion occurred rarely on the time scale of the simulation (6 ns). The Asn on which restraints B1 and B1 poseB are based is shown. Rendered with PyMOL.35 Three sets of Boresch restraints were selected initially for pose A (Figure 7). The first was the best-scoring set (B1) based on the minimum-variance algorithm, and mimics the phenol-Asn97C hydrogen bond (Figure 7). To test varied anchor point positions, the second set (B2) was selected as the top-scoring restraints with anchor points outwith the phenol moiety; these were based on the triazole ring and were seventh best-scoring overall. Finally, selection was constrained to the quinoline moiety (B3). The parameters of restraint sets B1, B2, and B3 are given in Table S1. The protein anchor points were located in residues with low root-mean-square fluctuations (RMSFs) of the α-carbon positions although this was not directly targeted by the algorithm (Figure S3). B1-poseB was the highest scoring set of anchor points based on a simulation including only pose B. Similar to the B1 restraints, these were based on the phenol group of MIF-180 and Asn97C. Figure 7 Sets of Boresch anchor points B1 (red), B2 (blue), and B3 (green). Anchor points are circled or shown as spheres. Windows rendered with PyMOL.35 A discussion of the challenges faced during restraint selection, symmetry corrections, and the strengths and limitations of the restraint selection algorithm is given in Section S6. Improvements to the restraints selection algorithm are also proposed; in particular, we recommend scoring possible restraints using a metric calculated from the variances of the DoF of the prospective restraints, rather than using the total variances directly. 4.1.2 Results with Force Constants Fit to Simulation Calculations were performed with Boresch restraints with force constants fit to simulation. Force constants were fit based on the variance of the Boresch DoF in State 7, as discussed in Section S7. Several calculations were performed for pose A (B1, B2, and B3) to allow comparison between restraints (Table 1). A single calculation was performed for pose B (B1-poseB) to allow overall comparison with experiment by combining the results for both poses. Table 1 Bound Leg Contributions to ΔGBindo with Boresch Restraintsa   Restraints Contribution B1 B2 B3 B1-P B3-Pb B1-10 B2-10 B3-10c B1-20 –ΔGReleaseo 9.76 9.90 9.10 10.02 8.91 6.69 6.39 6.53 7.90 –ΔGBound, Vanish –2.06 ± 1.08 –2.20 ± 0.94 –0.73 ± 0.25 –2.22 ± 0.19 –0.68 ± 0.06 –1.04 ± 1.21 0.09 ± 1.03 0.58 ± 1.96 –1.98 ± 1.18 –ΔGBound, Discharge –12.99 ± 0.63 –11.90 ± 0.48 –11.86 ± 0.22 –11.76 ± 0.19 –11.71 ± 0.77 –12.46 ± 0.59 –11.62 ± 0.47 –12.04 ± 0.70 –12.64 ± 0.19 –ΔGBound, Restrain –1.48 ± 0.04 –1.74 ± 0.03 –1.70 ± 0.11 –1.66 ± 0.08 –1.77 ± 0.11 –0.40 ± 0.01 –0.59 ± 0.08 –0.94 ± 0.21 –0.77 ± 0.20 ΔGSym. Corr. –1.06 –0.65 –0.65 –1.06 –0.65 –1.06 –0.65 –0.65 –1.06 ΔGBoundo –7.83 ± 1.25 –6.60 ± 1.05 –5.84 ± 0.36 –6.69 ± 0.28 –5.91 ± 0.78 –8.27 ± 1.35 –6.37 ± 1.13 –6.53 ± 2.09 –8.55 ± 1.21 a Results for binding pose A, in kcal mol–1. Uncertainties stated as 95% confidence intervals based on the variance of five replicate runs, assuming Gaussian distributions. Restraint parameters are given in Table S1; -10 and -20 denote that all force constants were set to 10 or 20 kcal mol–1 Å–2 [rad–2], respectively. b Run 1 excluded from average due to undersampling of water in binding site during the vanishing stage. c Due to simulations crashing, results based on four independent replicates with partial completion of many lambda windows. The results shown for −ΔGReleaseo were calculated by numerical integration of eq 6. While use of the Boresch analytical correction introduces large errors in certain regimes (force constants very weak, r very short, θA or θB close to 0 or π rad), only small deviations between the analytical and numerical corrections, no greater than 0.04 kcal mol–1, were found for any of the Boresch restraints. This is in accordance with previous studies.69 The average of repeat runs for B1-3 and B1-poseB generally showed good convergence as assessed by lack of drift with increasing sampling time (Section S8). However, there were substantial differences between replicate runs, most notably during the vanishing stage, which generally contributed the greatest uncertainty to ΔGBoundo. In all cases, this uncertainty could be traced back to a few windows around λ = 0.4 (Figure 8). Figure 8 (A) PMF along λ for the bound vanish stage for B1. (B) Average number of waters in the binding site (defined as the overlap of two spheres of radius 8 Å centered on the N atom in Pro1A, and CG2 in Val106A) against λ during the vanish stage for B1. There is a strong correlation between the water occupancy of the binding sites and the divergence of the PMFs around λ = 0.4. The shaded area shows the 95% confidence interval, and the solid blue line shows the mean. The majority of this uncertainty can be attributed to the entry of water to the binding site. Binding site water was defined as being simultaneously within 8 Å of the N atom in Pro1A, and CG2 in Val106A, which are on opposite sides of the binding site. At the start of the vanishing stage there were no waters in the binding site, which increased to an average of approximately 4.5 after decoupling, in good agreement with the crystal structure of free MIF (PDB ID 1GD0).70 A sudden jump in water occupancy of the binding site occurred around λ = 0.4, the region of divergence of the PMFs for B1. Here, it was found that water may only enter the binding site by “forcing” the partially vanished ligand to the side of the binding site (Figure 9). Figure 9 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 Å centered 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.35 In cases where the binding site showed high water occupancy (runs 1, 2, and 4), the resultant strain favored vanishing of the ligand, producing divergence of the vanish stage PMF toward more favorable free energies of vanishing. The opposite was true for runs 3 and 5, which had very low average water occupancies at λ = 0.4. This is in accord with the results of Rogers et al.,71 who found negative correlation between water occupancy and the gradient of the potential energy with respect to λ at intermediate stages of the vanishing leg. Although dramatic fluctuations in the number of water molecules in the binding sites were observed at higher values of λ, this did not translate into additional uncertainty, because at this stage the LJ terms are mostly removed and the ligand is able to pass through atoms relatively smoothly. Some dips in the water occupancy were found to be due to the rotation of the side chain of Met2A to obstruct the end of the binding site. 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 8, 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 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 ΔGBoundo 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.72 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.,40 is sensible. The B3 simulations were also rerun with protein anchor atoms on a different residue (Ala38A), see B3-P. There was very good agreement between runs 2–5 for the bound vanish stage, but ΔGBoundo was around 3 kcal mol–1 more negative for run 1 due to comparatively low number of waters in the binding site over just 3 λ windows (see Section S9). This highlights the importance of correctly sampling rehydration of the binding site upon decoupling.73 It has been demonstrated that hybrid sampling approaches combining molecular dynamics with Monte Carlo water moves, in both the μVT and NPT ensembles, can improve the performance of relative binding free energy calculations.74−76 In this system, it is possible that these approaches may perform poorly as a result of low acceptance probabilities, because water has to strain the partially decoupled ligand to enter the binding site around λ = 0.4, where proper sampling of rehydration is most critical. Nonequilibirum candidate Monte Carlo may overcome this by allowing for relaxation of the ligand position.73 The strong dependence of the free energy on the correct sampling of binding-site water may make this a good system for testing methods for enhancing sampling of rehydration. Regardless, the ΔGBoundo obtained after discarding the data from run 1 (−5.91 ± 0.78 kcal mol–1) was very similar to that obtained for B3 (−5.84 ± 0.36 kcal mol–1). Ignoring the results for B1, which seemed to be especially affected by water sampling issues, the remaining simulations with force constants fit to simulation showed generally good agreement (within 1 kcal mol–1), demonstrating reasonable reproducibility with this restraints scheme. Averaging B2, B3, B1-P, and B3-P yielded a ΔGBoundo of −6.26 ± 0.72 kcal mol–1 for binding pose A. Although the aim of this study was to compare the values of ΔGBoundo obtained with different restraints, a single set of five replicate calculations were carried out for the free leg to allow comparison with experiment. Comparison with the experiment cannot be used for the comparison of restraint schemes as there are likely to be systematic errors from the force field, but it is used here as a crude check on the overall results. There was excellent agreement between replicates of the free leg simulations, and convergence was achieved quickly (Section S10), yielding ΔGFree = −3.08 ± 0.14 kcal mol–1 (ΔGFree, Discharge = 9.03 ± 0.07 kcal mol–1, and ΔGFree, Vanish = −12.11 ± 0.11 kcal mol–1). This was in spite of syn–anti interconversion, which was observed during at least one run for every lambda window, but which occurred slowly on the time scale of the simulations. ΔGBoundo for binding pose B (with the B1-poseB restraints) was −6.63 ± 1.11 kcal mol–1 (−ΔGReleaseo = 9.76 kcal mol–1, – ΔGBound, Discharge = −2.66 ± 1.04 kcal mol–1, – ΔGBound, Vanish = −10.77 ± 0.39 kcal mol–1, – ΔGBound, Restrain = −1.90 ± 0.1 kcal mol–1, and ΔGSym. Corr. = −1.06 kcal mol–1). Combining the average result for pose A (−6.28 ± 0.49 kcal mol–1, using all results from Table 1 excluding all B1 calculations and B3-10 due to the issues discussed) with the result for pose B according to ΔGBoundo = −kBT ln (exp(−βΔGBound, 1o) + exp(−βΔGBound, 2o)) yielded an overall ΔGBoundo of −6.89 ± 0.74 kcal mol–1 and ΔGBindo = −9.97 ± 0.76 kcal mol–1.37 The overall result was in good agreement with the experimental binding free energy of −8.98 ± 0.28 kcal mol–1,77 but more negative than the value calculated by Qian et al. using molecular dynamics and the AMBER ff14SB and GAFF force fields (−7.47 ± 0.99 kcal mol–1).34 However, a clear comparison with the results of Qian et al. is prevented by a number of methodological differences. For example, they only observe the anti and syn conformers of MIF180 during the free and bound legs, respectively, and they apply a penalty of 1.60 kcal mol–1 to account for this. Here, interconversion was observed in both the free and bound states, and no correction was applied. Furthermore, Qian et al. do not apply a symmetry correction to account for the 3-fold symmetry of MIF, and do not perform calculations for an alternative binding pose. 4.1.3 Results without Orientational Restraints It was found that orientational restraints were essential for achieving reliable free energy estimates. The requirement for orientational restraints was investigated by repeating the B1 and B2 calculations without the orientational component of the restraint (KθB, KϕB, and KϕC were set to 0), see B1-o and B2-o in Table 2. B1 was also repeated setting all force constants other than Kr to 0, thus retaining only a single distance restraint, see B1-d. This resulted in large and significant shifts to more negative free energies of binding by 1.72, 3.84, and 2.05 kcal mol–1, for B1-o, B2-o, and B1-d, respectively. Despite this, there was no obvious drift of the free energies with increasing sampling time (Section S11). Table 2 Bound Leg Contributions to ΔGBindo without Orientational Restraintsa   Restraints   B1-o B2-o B1-d –ΔGReleaseo 4.94 4.51 1.08 –ΔGBound, Vanish –0.60 ± 0.71 –1.47 ± 0.99 1.67 ± 1.11 –ΔGBound, Discharge –12.50 ± 0.40 –11.93 ± 0.20 –11.89 ± 0.44 –ΔGBound, Restrain –0.74 ± 0.01 –0.91 ± 0.14 –0.09 ± 0.00 ΔGSym. Corr. –0.65 –0.65 –0.65 ΔGBoundo –9.55 ± 0.82 –10.44 ± 1.02 –9.88 ± 1.19 a All values in kcal mol–1. Uncertainties stated as 95% confidence intervals based on the variance of five replicate runs, assuming Gaussian distributions; -o and -d mean that all force constants other than kr, kθA, and kϕA, or kr were set to 0, respectively. The lack of orientational restraints allows the mixing of binding poses A and B, but this is also not a plausible source of the error introduced. In the limit of perfect sampling, the free energy of binding can be no more negative than that calculated by combining the free energies of binding of the two poses as was done previously. The exception to this would be if there were other binding poses which were numerous or more favorable, which seems unlikely. Instead, the negative offset is very likely due to the failure to sample all relevant orientations at intermediate values of λ during vanishing. The close agreement between B1-o and B1-d suggests that the offset is due to the removal of the orientational component of the restraint. As the LJ interactions are removed, the sampling of orientations different to that of the binding pose will become favorable. The gradient of the free energy change as the LJ interactions are removed will likely be less positive in these alternative orientations, because they were high in energy when the LJ interactions were at full strength. Therefore, failure to sample these orientations due to large barriers should give excessively positive free energies of vanishing, resulting in erroneously negative free energies of binding, as is observed. This is illustrated by the divergence of the PMF for the bound vanish stage for B2-o (Figure 10). The divergence shows a strong correlation with the orientational sampling at λ = 0.325, and not with the presence of water in the binding site (Figure S16). For run 4, the plots of ϕC show that sampling was largely restricted to the orientation of the original binding mode. This resulted in the most positive gradient of the PMF and the most favorable free energy of binding. During runs 1, 2, and 5, the ligand rotated lengthwise in the binding site by around 90° with respect to the initial pose, resulting in a less positive PMF gradient. During run 3, the ligand rotated around 180° lengthwise in the binding site with respect to the initial pose. The gradient here was evidently even less positive, resulting in a substantially less negative free energy of binding. The slow interconversion between orientations explains the lack of drift of the results with increasing simulation time; there are large barriers between orientations which prevent equilibrium sampling on the time scale of the simulations. To support this explanation, the free energy of releasing B1 to B1-o in the decoupled state was calculated by scaling the strength of the orientational force constants with λ, using the same set of λ windows as for the bound vanish stage. The free energy difference for releasing the restraints calculated by simulation was 4.92 ± 0.02 kcal mol–1 (MBAR 95% C.I. estimate for a single simulation), in close agreement with the difference calculated by numerical integration by subtracting the two ΔGReleaseo terms (4.82 kcal mol–1). This confirms that sampling a large increase in configurational space is not problematic when there are no barriers, at least when the growth is sufficiently slow; the issue is sampling the “rugged” configurational space at intermediate stages of decoupling. Comparison to prior work on the use of orientational restraints is given in Section S13. Figure 10 (A) The PMFs for the bound vanish stage for B2-o. These diverge around λ = 0.325. (B) The unrestrained Boresch DoF ϕC for B2-o at λ = 0.325, showing the presence of multiple slowly interconverting orientations. Snapshots taken from run 3 (0.06 ns), run 5 (5.06 ns), and run 3 (7.22 ns), from left to right. The anchor points used to define ϕC are shown as spheres. Windows rendered with PyMOL.35 The offset introduced by the removal of orientational restraints may be reduced or removed through the use of the Hamiltonian-replica exchange (HREX) method. This is because in the fully decoupled state, there are no barriers to orientational rearrangements, and HREX has been found to improve sampling when the barriers in configurational space are low in at least one state.78 This allows free sampling of varied orientations, and mixing of these into the intermediate λ states using HREX may improve orientational sampling. However, it may not remove the bias; Lapelosa et al. found that convergence of their HREX ABFE calculations for a large and flexible ligand could only be achieved with orientational restraints.79 4.1.4 Performance of Common Default Force Constants The performance of the Boresch restraints with the force constants fit to simulation were compared to those with all force constants set to the common defaults of 10 or 20 kcal mol–1 Å–2 [rad–2], denoted by -10 and -20 (Table 1).37,80,81 No significant differences in ΔGBoundo were found compared to when the force constants were fit to simulation, which was unsurprising given that the theoretical independence of the binding free energy with respect to the strength of restraints has been previously confirmed.17,69 If there were any improvements in precision or increases in the rate of convergence with the force constants fit to simulation, these were not observed above the noise generated by other sources of error. The only difference when default force constants were used was that several simulations crashed, very likely due to the collinearity of contiguous anchor points. While all simulations completed successfully when the force constants were fit to simulation, 1 λ window failed for B2-10 (this was rerun) and 12 failed for B3-10. For B3-10, the minimum energy penalty arising from the restraints for setting a θA or θB to 0 or 180° was approximately 5 kBT, meaning that collinearity was relatively likely and crashes may have been anticipated. However, for B2-10, the minimum penalty for collinearity was approximately 10 kBT. This makes crashes in the decoupled state highly unlikely, but when the ligand is still interacting with the protein it may become trapped in unusual orientations which distort the anchor points toward collinearity. This was the cause of the crash for B3-10 run 1: during the vanish λ = 0.475 window, the simulation failed after θA approaches 0 (Figure 11). This occurred because the ligand became trapped underneath the terminal proline residue, resulting in θA approaching collinearity. A better restraint selection algorithm would have discounted restraints if the energy penalty from the restraints for collinearity was below some threshold, rather than only checking the equilibrium angles. Regardless, this highlights that instabilities with Boresch restraints can be an issue even when sensible restraint selections are made. It also illustrates that fitting the force constants to simulation (or at least using higher force constants) can produce more stable restraints. Figure 11 Restrained angle θA for B3–10 at λ = 0.475 during the vanish stage. For run 1, θA tended toward 0 as the ligand became trapped under the terminal proline and the simulation crashed. Anchor points used in the definition of θA are shown as spheres. Snapshots taken from run 1 at 1.12 ns (upper image) and 5.30 ns (lower image). Windows rendered with PyMOL.35 4.2 Multiple Distance Restraints 4.2.1 With Intramolecular Rigidification Multiple distance restraints provide a framework which is free from the inherent instabilities of Boresch restraints, and which allows more complete restraint of the ligand than Boresch restraints alone. However, their naive application renders the ABFE framework theoretically inexact. To illustrate this, harmonic distance restraints were applied to every heavy atom in the ligand and their lowest variance unique partner heavy atom in the protein (protocol M-All, 22 receptor–ligand distance restraints). As expected, this produced an excessively negative ΔGBoundo estimate, in excess of 3 kcal mol–1 more negative than most of the Boresch restraints, indicating that ΔGPreorg. was large (Table 3). The distance restraint dictionaries used for all protocols are given in Section S14. Table 3 Bound Leg Contributions to ΔGBindo with Multiple Distance Restraintsa   Restraints Contribution M-All M-Rig M-Rig-Nb M-All-R M-Hand-R M-Hand M-Hand-1 ΔGRigid. Lig. – 0.50 ± 0.00 – – – – – ΔGRigid. Recept. – 10.36 ± 0.09 – – – – – –ΔGReleaseo 15.68 ± 0.37 10.05 ± 0.17 9.97 ± 0.03 2.40 1.44 4.35 ± 0.27 5.82 ± 0.26 –ΔGTo Dist. Rest. – – – 16.66 ± 0.29 3.07 ± 0.24 – – –ΔGBound, Vanish –7.71 ± 1.05 –2.48 ± 1.25 –2.83 ± 0.90 –7.71 ± 1.05 1.58 ± 1.05 1.58 ± 1.05 1.81 ± 1.06 –ΔGBound, Discharge –13.54 ± 0.45 –13.23 ± 0.46 –12.97 ± 0.10 –13.54 ± 0.45 –11.89 ± 0.43 –11.89 ± 0.43 –12.99 ± 0.64 –ΔGBound, Restrain –3.67 ± 0.10 –1.33 ± 0.03 –1.35 ± 0.03 –3.67 ± 0.10 –0.02 ± 0.02 –0.02 ± 0.02 –0.28 ± 0.12 –ΔGRigid. Complex – –8.31 ± 0.08 – – – – – ΔGSym. Corr. –1.06 –1.06 –1.06 –0.65 –0.65 –0.65 –0.65 ΔGBoundo –10.31 ± 1.20 –5.50 ± 1.35 –8.24 ± 0.91 –6.51 ± 1.18 –6.47 ± 1.16 –6.64 ± 1.16 –6.30 ± 1.27 a Results for binding pose A, in kcal mol–1. Uncertainties stated as 95% confidence intervals based on the variance of five replicate runs assuming Gaussian distributions. The distance restraint dictionaries used for all protocols are given in Section S14; -N indicates that the protocol was repeated with no intramolecular rigidification, -R repetition with release to the single strongest distance restraint, and -1 repetition with the flat-bottomed diameter set to 1 Å for all restraints. ΔGTo Dist. Rest. is the free energy of releasing the multiple distance restraints to a single distance restraint. b 10 replicate runs were used. The convergence of −ΔGReleaseo with respect to the parameters of the Sire standard state correction utility was confirmed (Section S15). This was contrasted with the rigorous multiple distance restraints scheme with intramolecular rigidification of anchor points (Figure 12), referred to as M-Rig. Anchor points were selected as for M-All, except that only anchor points in the phenol moiety of the ligand were used in order to avoid the tight restraints restricting rotatable bonds in the ligand; this would likely require an enhanced sampling approach, such as umbrella sampling,82 to restrain bond rotation before applying the restrictive intermolecular restraints. This resulted in seven receptor–ligand distance restraints. The intramolecular restraints were implemented as harmonic distance restraints between all pairs of anchor atoms within the given molecule, with force constants of 75 kcal mol–1 Å–2. For simulations where −ΔGReleaseo was calculated using eq 10, it was confirmed that the estimate had converged with respect to the number of points used for numerical integration (Section S15). This scheme gave ΔGBoundo = −5.50 ± 1.35 kcal mol –1, in good agreement with the previous Boresch calculations, providing proof-of-concept of a rigorous implementation of multiple distance restraints. Figure 12 Anchor points (red) used for the M-Rig restraints. Intermolecular restraints are shown as red dashed lines, while intramolecular restraints not shown. Rendered with PyMOL.35 The magnitude of ΔGPreorg. can be roughly estimated as ΔGRigid. Recept. + ΔGRigid. Lig. – ΔGRigid. Complex = 2.55 ± 0.17 kcal mol–1 (see Section S18). This is in excellent agreement with the difference observed when the procedure was repeated without intramolecular restraints (M-Rig-N), which was 2.74 ± 1.63 kcal mol–1, providing evidence that this is the main source of the error observed when multiple distance restraints are applied in a nonrigorous manner. The intramolecular restraints are required to prevent distortion of the intramolecular degrees of freedom by the intermolecular restraints, which would otherwise introduce error into the ΔGReleaseo calculated using eq 10. The strength of the intramolecular restraints required to eliminate this error is dependent on how permissive the intermolecular restraints are: highly restrictive intermolecular restraints necessitate aggressive rigidification. For the M-Rig calculations, strong intermolecular restraints were used, and therefore, strong intramolecular restraints were also required. The large value of ΔGPreorg. which resulted from the strong intermolecular restraints allowed us to show that that intramolecular rigidification eliminates this error, by comparison of M-Rig and M-Rig-N. However, the strong intermolecular restraints resulted in a significant additional computational cost during the rigidification stage. The λ protocols for the rigidification stages could likely have been optimized to reduce the total number of windows to approximately 25 (Section S17), but this would still result in an additional computational cost comparable to the vanishing stage. Hence, we would not recommend multiple distance restraints with intramolecular rigidification when using strong intermolecular restraints. Instead, this scheme is expected to perform best with more permissive intermolecular restraints. An efficient implementation may be as follows:(1) Select flat-bottomed, rather than harmonic, intermolecular restraints and select the flat-bottomed regions so to include almost all distances sampled during the unrestrained simulation of the fully interacting complex. (2) Perform the rigidification simulations simultaneously starting from States 3 and 4. Monitor the convergence of ΔGPreorg. as the intramolecular restraint strength is increased and stop the calculations once convergence is achieved. (3) Calculate ΔGReleaseo with eq 10 based on the trajectory for the most strongly rigidified version of State 4. In cases where such flat-bottomed restraints are used, and there is little rearrangement of the ligand or binding site upon decoupling, ΔGPreorg. would be expected to converge immediately and the above scheme would reduce to the naive multiple distance restraints scheme. When this is not the case, the naive scheme would be expected to yield incorrect results, while the above scheme should remain correct. 4.2.2 With Release to Single Distance Restraint Based on the approach taken with DBC restraints,28,29 the free energy of releasing all but the strongest distance restraint after decoupling was calculated for M-All (M-All-R). Because a single distance restraint does not couple the internal and external degrees of freedom of the protein and ligand, this allows the rigorous calculation of −ΔGReleaseo using eq 11, removing the requirement for intramolecular restraints. In contrast to M-All, M-All-R produced a ΔGBoundo of −6.51 ± 1.18 kcal mol–1, in good agreement with the Boresch calculations. There appeared to be a slight drift in ΔGTo Dist. Rest. with time (Figure S20) toward more negative values, which may be due to the requirement for the ligand to sample all points on the surface of a sphere surrounding the protein anchor point upon decoupling. This may be improved by releasing to a center-of-mass restraint centered on the binding site, reducing the volume which must be sampled.29 In addition, scaling the restraint potential differently between the end points may improve convergence; instead of the λ5 scaling used here, a soft bond stretch potential would likely perform well for the removal of these harmonic restraints.83 This process was repeated using a significantly more permissive multiple distance restraint scheme, where four flat-bottomed distance restraints were selected to mimic the four protein–ligand hydrogen bonds shown in Figure 1. This scheme is denoted M-Hand-R as the anchor points were selected by hand, although automated selection to match hydrogen bonds would be straightforward. The radius of the flat-bottomed region was selected to be as small as possible without the restraints engaging at any point during the restraint fitting simulation, ensuring that ΔGBound, Restrain ≈ 0. The force constants for the half-harmonic potentials were 40 kcal mol–1. This again produced a ΔGBoundo in good agreement with the Boresch results and M-Rig (−6.47 ± 1.16 kcal mol–1), indicating that the relatively permissive restraints sufficiently restricted orientational sampling. Furthermore, no substantial drift in the estimate with simulation time was observed (Figure S21). Compared to the naive and Boresch schemes, this scheme required a single additional release stage, which was relatively computationally affordable. Furthermore, computational cost could be further reduced with optimization of the λ schedule for M-Hand-R, because excellent overlap was achieved between many nonconsecutive λ windows (Section S17). The overlap matrices for the vanishing stages were very similar regardless of the restraint scheme, and the number of windows required for the vanishing stages could not have been reduced below approximately 30 in any case (see Figure S24 as a representative example). In contrast, for both M-Hand-R and M-Rig-R, it appears that 10 windows would be sufficient, or even fewer in the case of M-Hand-R. For M-Hand-R, convergence appears to be achieved after around 1 ns sampling per window (not including the 1 ns equilibration), which seems broadly similar to the bound vanish stage results (not including the 3 ns equilibration). Therefore, when the intermolecular restraints are not extremely numerous and strong, the additional cost associated with the release stage with an optimized λ schedule is expected to be substantially less than a third of the vanish stage. Although there was no evidence for improved convergence of the decoupling simulations with M-All-R over the Boresch schemes, this might be observed in other systems with highly flexible ligands, where ligand conformational sampling is the dominant source of uncertainty. However, the use of Boresch restraints in combination with RMSD-based restraints on the conformation of the ligand may prove similarly effective.82 This scheme is in some ways similar to the DBC scheme, in that a complex restraint (set of restraints) involving many degrees of freedom is released to a single harmonic restraint to allow the calculation of ΔGReleaseo. While the DBC restraint is attractively simple—it consists of a single flat-bottomed restraint on the RMSD of a subset of ligand coordinates in the frame of reference of the binding site—multiple distance restraints schemes offer finer control over the strength of restraints applied to different subsections of the system. This may be beneficial, for example, in the case of a large and flexible ligand where only part of the ligand interacts strongly with the receptor. If the coordinates of all ligand heavy atoms were included in the DBC restraint, then the DBC coordinate would show very wide fluctuations during simulations of the bound state. Fitting the DBC restraint to encompass the 95th percentile of sampled DBC coordinates would result in a very weak restraint on all sections of the ligand, which may result in sampling issues. Flat-bottomed multiple distance restraints could be fit in a similar way, such that the flat-bottomed regions encompassed almost all distances measured during a simulation of the fully interacting complex. In contrast to the DBC restraints fit to all ligand heavy atoms, the multiple distance restraints fit to all heavy atoms would closely restrict the portion of the ligand which interacts strongly with the receptor, while allowing large fluctuations in the flexible portion which does not. While multiple distance restraints require many more parameters than the DBC restraint, they can be automatically selected from a simulation of the fully interacting complex using algorithms described in this work. Furthermore, if the user opted not to run such a simulation for restraint selection, distance restraints could be intuitively selected to match receptor–ligand interactions, and reasonable default parameters could be chosen. It may be challenging to select a reasonable flat-bottomed region for a DBC restraint without a trajectory. Averaging over the rigorous multiple distance restraint schemes (M-Rig, M-All-R, and M-Hand-R), the mean ΔGBoundo result for pose A was −6.16 ± 1.43 kcal mol–1. This was in good agreement with the mean result for pose A calculated with Boresch restraints (−6.28 ± 0.49 kcal mol–1, average of B2, B3, B1-P, and B3-P, and B2-10). Ignoring the contribution to the free energy of binding from pose B, for which no calculations were performed with multiple distance restraints, the free energies of binding would have been −9.24 ± 1.44 and −9.36 ± 0.37 kcal mol–1 for multiple distance restraints and Boresch restraints, respectively. 4.2.3 With a Large Flat-Bottomed Region Finally, a nonrigorous implementation of multiple distance restraints was tested, based on the assumption of no coupling between internal and external degrees of freedom in the limit of weak restraints. The schemes tested were M-Hand, a repetition of M-Hand-R without releasing to a single distance restraint, and M-Hand-1, a repetition of M-Hand with all flat-bottomed diameters set to 1 Å, a reduction from the average diameter of 2.7 Å for M-Hand-R. A ΔGBoundo of −6.64 ± 1.16 kcal mol –1 was obtained for M-Hand. This was not significantly different to M-Hand-R or M-Rig, suggesting that the approximations made by the scheme were minor in this case. The results for M-Hand and M-Hand-1 (−6.30 ± 1.27 kcal mol–1) were very similar, despite the slightly more restrictive restraints. This shows that it is possible to obtain equivalent free energies using the naive distance restraints scheme and rigorous schemes, so long as the coupling between the internal and relative external DoF of the protein and ligand is weak when the ligand is decoupled. However, increasingly negative binding free energies would be expected as more restrictive restraints are used and with increasing differences between the apo and holo conformations of the binding site, and the magnitude of the error may be be difficult to predict. 5 Conclusion The free energies of binding for MIF-180 to MIF calculated with varied sets of Boresch restraints were fairly self-consistent and in good agreement with experiment. However, removal of the orientational restraints produced estimates which were up to approximately 4 kcal mol–1 more negative, likely because the ligand failed to sample all relevant orientations as the LJ interactions were removed. It was found that the calculations were highly sensitive to the sampling of water in the binding site at intermediate stages of vanishing and that underhydration of the binding site during these stages over as few as 3 λ windows could shift the binding free energy by over 3 kcal mol–1 toward more favorable binding. Instabilities inherent to the Boresch restraint scheme were highlighted by the failure of several simulations, even with sensible restraint parameters which imposed a minimum energy penalty of around 10 kBT for collinear anchor points. The use of multiple distance restraints offers an alternative restraint scheme which lacks the instabilities of Boresch restraints and may improve convergence during decoupling by allowing greater restriction of ligand movements. The theory of multiple distance restraints was discussed, and a rigorous implementation of the multiple distance restraints scheme was proposed. This utilized intramolecular rigidification of anchor points to prevent coupling between the internal and external DoF. This was shown to produce free energy estimates in good agreement with the Boresch restraints, at least within the large uncertainties encountered with this system (Figure 13). This scheme incurred a substantial additional computational cost over Boresch restraints because aggressive rigidification was required to counter the strong intermolecular restraints, but the scheme may offer benefits where more permissive intermolecular restraints are used. Figure 13 Summary of results for ΔGBoundo obtained for binding pose A using a variety of restraints schemes. Uncertainties are 95% confidence intervals based on the variance of five replicate runs, assuming Gaussian distributions. Results for B1 have been omitted as they appeared to be more susceptible to water sampling issues. For B3-P, run 1 was excluded from the average due to undersampling of water in binding site during the vanishing stage. Another rigorous implementation of the multiple distance restraints scheme was tested, which involved releasing the multiple restraints to a single distance restraint after decoupling. In contrast to calculations performed entirely without orientational restraints, this scheme produced free energy estimates in good agreement with the Boresch restraints scheme, at a reduced computational cost compared to the scheme employing intramolecular restraints. The additional computational cost compared to Boresch restraints is expected to be less than a third of the vanish stage unless very many strong intermolecular restraints are used. Additional costs associated with rigorous multiple distance restraints schemes may be compensated for by convergence benefits in some systems, although that was not demonstrated in this work. The mean ΔGBoundo calculated based on pose A with rigorous implementations of multiple distance restraints (−6.16 ± 1.43 kcal mol–1) was in close agreement with the mean result calculated with Boresch restraints (−6.28 ± 0.49 kcal mol–1). Finally, a nonrigorous implementation of the multiple distance restraints scheme was tested, which relied on the assumption of negligible coupling between the internal and external DoF. With strong restraints, this assumption was violated and excessively negative free energies of binding were calculated, but quantities close to the rigorous estimates were obtained with sufficiently permissive restraints. However, it may be difficult to predict the magnitude of the error introduced by this implementation. The dominant source of uncertainty in these calculations appeared to be the sampling of water, and no convergence benefits were demonstrated with multiple distance restraints over Boresch restraints. Future work may investigate whether convergence benefits are observed in systems where ligand conformational sampling is the dominant source of uncertainty. Further comparison against a wider range of restraint schemes over a variety of systems is also the subject of future work. In summary, this work demonstrates that absolute binding free energies equivalent to those obtained with Boresch restraints can be calculated using multiple distance restraints. This framework is in principle more stable and may offer convergence benefits during decoupling, although this must be balanced against the additional computational cost incurred by the extra stages required for the rigorous schemes. However, a multiple distance restraints scheme utilizing many flat-bottomed potentials to closely restrain the ligand with minimal disturbance of the interacting system may allow the restraining stage to be neglected, as with the DBC restraint,29 while improving the convergence of the decoupling stages when the ligand is flexible. This work discussed the theory and demonstrated proof-of-principle of rigorous multiple distance restraint schemes; future work may investigate whether the scheme can offer performance benefits. Data Availability Statement Example input files, restraint selection code, and modified version of Sire required for the rigorous multiple distance restraints calculations are available on GitHub at michellab/Multiple-Distance-Restraints.. Supporting Information Available The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jctc.3c00139.Discussion of the theoretical requirement for receptor–ligand restraints, derivation of eq 4, diagram of the syn and anti conformations of MIF-180, parameters for all restraint sets, RMSFs of residues containing anchor points for Boresch restraints, discussion of challenges and potential improvements related to the restraint selection procedure, discussion of the force constant fitting procedure, convergence of free energies with simulation time, bound vanish stage PMF and plot of number of waters in the binding site for B3-P, plot of the number of waters in the binding site during the bound vaish stage for B2-o, comparison of the results without orientational restraints to the literature, results showing convergence of restraint correction for multiple distance restraints schemes, overlap matrices for selected multiple distance restraint protocols, and convergence of the free energy of preorganization for M-Rig. (PDF) Supplementary Material ct3c00139_si_001.pdf The authors declare the following competing financial interest(s): J.M. is a member of the Scientific Advisory Board of Cresset. Acknowledgments F.C. thanks AstraZeneca and the Engineering and Physical Sciences Research Council for Ph.D. studentship funding. D.J.C. acknowledges support from a UKRI Future Leaders Fellowship (Grant MR/T019654/1). Gratitude is expressed to Lester Hedges for support implementing ABFE restraints in SOMD and BioSimSpace. ==== Refs References Cournia Z. ; Allen B. ; Sherman W. Relative Binding Free Energy Calculations in Drug Discovery: Recent Advances and Practical Considerations. J. Chem. Inf. Model. 2017, 57 , 2911–2937. 10.1021/acs.jcim.7b00564.29243483 Salomon-Ferrer R. ; Götz A. W. ; Poole D. ; Le Grand S. ; Walker R. C. Routine Microsecond Molecular Dynamics Simulations with AMBER on GPUs. 2. Explicit Solvent Particle Mesh Ewald. J. Chem. Theory Comput. 2013, 9 , 3878–3888. 10.1021/ct400314y.26592383 Heinzelmann G. ; Gilson M. K. Automation of Absolute Protein-Ligand Binding Free Energy Calculations for Docking Refinement and Compound Evaluation. Sci. Rep. 2021, 11 , 1116 10.1038/s41598-020-80769-1.33441879 Fu H. ; Chen H. ; Cai W. ; Shao X. ; Chipot C. BFEE2: Automated, Streamlined, and Accurate Absolute Binding Free-Energy Calculations. J. Chem. Inf. Model. 2021, 61 , 2116–2123. 10.1021/acs.jcim.1c00269.33906354 Rizzi A. ; Chodera J. ; Naden L. ; Beauchamp K. ; Albanese S. ; Grinaway P. ; Prada-Gracia D. ; Rustenburg B. ; Ajsilveira ; Saladi S. ; Boehm K. ; Gmach J. ; Rodríguez-Guerra J. Choderalab/Yank: 0.25.2 - Bugfix Release. Zenodo, 2019. 10.5281/zenodo.3534289. Huggins D. J. ; Biggin P. C. ; Dämgen M. A. ; Essex J. W. ; Harris S. A. ; Henchman R. H. ; Khalid S. ; Kuzmanic A. ; Laughton C. A. ; Michel J. ; et al. Biomolecular simulations: From dynamics and mechanisms to computational assays of biological activity. WIREs Comput. Mol. Sci. 2019, 9 , e1393 10.1002/wcms.1393. De Simone A. ; Georgiou C. ; Ioannidis H. ; Gupta A. A. ; Juárez-Jiménez J. ; Doughty-Shenton D. ; Blackburn E. A. ; Wear M. A. ; Richards J. P. ; Barlow P. N. ; et al. A computationally designed binding mode flip leads to a novel class of potent tri-vector cyclophilin inhibitors. Chem. Sci. 2019, 10 , 542–547. 10.1039/C8SC03831G.30746096 Cournia Z. ; Allen B. K. ; Beuming T. ; Pearlman D. A. ; Radak B. K. ; Sherman W. Rigorous Free Energy Simulations in Virtual Screening. J. Chem. Inf. Model. 2020, 60 , 4153–4169. 10.1021/acs.jcim.0c00116.32539386 Gumbart J. C. ; Roux B. ; Chipot C. Standard Binding Free Energies from Computer Simulations: What Is the Best Strategy?. J. Chem. Theory Comput. 2013, 9 , 794–802. 10.1021/ct3008099.23794960 Breznik M. ; Ge Y. ; Bluck J. P. ; Briem H. ; Hahn D. F. ; Christ C. D. ; Mortier J. ; Mobley D. L. ; Meier K. Prioritizing Small Sets of Molecules for Synthesis through In-silico Tools: A Comparison of Common Ranking Methods. ChemMedChem. 2023, 18 , e202200425 10.1002/cmdc.202200425.36240514 Mey A. S. ; Allen B. K. ; Bruce Macdonald H. E. ; Chodera J. D. ; Hahn D. F. ; Kuhn M. ; Michel J. ; Mobley D. L. ; Naden L. N. ; Prasad S. ; Rizzi A. ; Scheen J. ; Shirts M. R. ; Tresadern G. ; Xu H. Best Practices for Alchemical Free Energy Calculations [Article v1.0]. LiveCoMS 2020, 2 , 18378 10.33011/livecoms.2.1.18378.34458687 Mey A. S. ; Juárez-Jiménez J. ; Hennessy A. ; Michel J. Blinded predictions of binding modes and energies of HSP90-α ligands for the 2015 D3R grand challenge. Bioorg. Med. Chem. 2016, 24 , 4890–4899. 10.1016/j.bmc.2016.07.044.27485604 Mey A. S. ; Jiménez J. J. ; Michel J. Impact of domain knowledge on blinded predictions of binding energies by alchemical free energy calculations. J. Comput. Aided Mol. Des. 2018, 32 , 199–210. 10.1007/s10822-017-0083-9.29134431 Loeffler H. H. ; Bosisio S. ; Duarte Ramos Matos G. ; Suh D. ; Roux B. ; Mobley D. L. ; Michel J. Reproducibility of free energy calculations across different molecular simulation software packages. J. Chem. Theory Comput. 2018, 14 , 5567–5582. 10.1021/acs.jctc.8b00544.30289712 Granadino-Roldan J. M. ; Mey A. S. ; Pérez González J. J. ; Bosisio S. ; Rubio-Martinez J. ; Michel J. Effect of set up protocols on the accuracy of alchemical free energy calculation over a set of ACK1 inhibitors. PLoS One 2019, 14 , e0213217 10.1371/journal.pone.0213217.30861030 Lee T.-S. ; Allen B. K. ; Giese T. J. ; Guo Z. ; Li P. ; Lin C. ; McGee T. D. ; Pearlman D. A. ; Radak B. K. ; Tao Y. ; Tsai H.-C. ; Xu H. ; Sherman W. ; York D. M. Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery. J. Chem. Inf. Model. 2020, 60 , 5595–5623. 10.1021/acs.jcim.0c00613.32936637 Boresch S. ; Tettinger F. ; Leitgeb M. ; Karplus M. Absolute Binding Free Energies: A Quantitative Approach for Their Calculation. J. Phys. Chem. B 2003, 107 , 9535–9551. 10.1021/jp0217839. Chen W. ; Cui D. ; Jerome S. V. ; Michino M. ; Lenselink E. B. ; Huggins D. ; Beautrait A. ; Vendome J. ; Abel R. ; Friesner R. A. ; Wang L. Enhancing Hit Discovery in Virtual Screening through Absolute Protein–Ligand Binding Free-Energy Calculations. J. Chem. Inf. Model. 2023, 63 , 3171 10.1021/acs.jcim.3c00013.37167486 Rizzi A. ; Jensen T. ; Slochower D. R. ; Aldeghi M. ; Gapsys V. ; Ntekoumes D. ; Bosisio S. ; Papadourakis M. ; Henriksen N. M. ; de Groot B. L. ; Cournia Z. ; Dickson A. ; Michel J. ; Gilson M. K. ; Shirts M. R. ; Mobley D. L. ; Chodera J. D. The SAMPL6 SAMPLing Challenge: Assessing the Reliability and Efficiency of Binding Free Energy Calculations. J. Comput. Aided Mol. Des. 2020, 34 , 601–633. 10.1007/s10822-020-00290-5.31984465 Georgiou C. ; McNae I. ; Wear M. ; Ioannidis H. ; Michel J. ; Walkinshaw M. Pushing the Limits of Detection of Weak Binding Using Fragment-Based Drug Discovery: Identification of New Cyclophilin Binders. J. Mol. Biol. 2017, 429 , 2556–2570. 10.1016/j.jmb.2017.06.016.28673552 Baumann H. M. ; Dybeck E. ; McClendon C. L. ; Pickard IV F. C. ; Gapsys V. ; Pérez-Benito L. ; Hahn D. F. ; Tresadern G. ; Mathiowetz A. M. ; Mobley D. L. Broadening the Scope of Binding Free Energy Calculations Using a Separated Topologies Approach. ChemRxiv Preprint, 2023. 10.26434/chemrxiv-2023-9678k Hermans J. ; Shankar S. The Free Energy of Xenon Binding to Myoglobin from Molecular Dynamics Simulation. Isr. J. Chem. 1986, 27 , 225–227. 10.1002/ijch.198600032. Roux B. ; Nina M. ; Pomès R. ; Smith J. Thermodynamic Stability of Water Molecules in the Bacteriorhodopsin Proton Channel: A Molecular Dynamics Free Energy Perturbation Study. Biophys. J. 1996, 71 , 670–681. 10.1016/S0006-3495(96)79267-6.8842206 Hermans J. ; Wang L. Inclusion of Loss of Translational and Rotational Freedom in Theoretical Estimates of Free Energies of Binding. Application to a Complex of Benzene and Mutant T4 Lysozyme. J. Am. Chem. Soc. 1997, 119 , 2707–2714. 10.1021/ja963568+. Procacci P. ; Macchiagodena M. On the NS-DSSB Unidirectional Estimates in the SAMPL6 SAMPLing Challenge. J. Comput. Aided Mol. Des. 2021, 35 , 1055–1065. 10.1007/s10822-021-00419-0.34625885 Wang J. ; Deng Y. ; Roux B. Absolute Binding Free Energy Calculations Using Molecular Dynamics Simulations with Restraining Potentials. Biophys. J. 2006, 91 , 2798–2814. 10.1529/biophysj.106.084301.16844742 Fu H. ; Cai W. ; Hénin J. ; Roux B. ; Chipot C. New Coarse Variables for the Accurate Determination of Standard Binding Free Energies. J. Chem. Theory Comput. 2017, 13 , 5173–5178. 10.1021/acs.jctc.7b00791.28965398 Salari R. ; Joseph T. ; Lohia R. ; Hénin J. ; Brannigan G. A Streamlined, General Approach for Computing Ligand Binding Free Energies and Its Application to GPCR-Bound Cholesterol. J. Chem. Theory Comput. 2018, 14 , 6560–6573. 10.1021/acs.jctc.8b00447.30358394 Ebrahimi M. ; Hénin J. Symmetry-Adapted Restraints for Binding Free Energy Calculations. J. Chem. Theory Comput. 2022, 18 , 2494–2502. 10.1021/acs.jctc.1c01235.35230113 Santiago-McRae E. ; Ebrahimi M. ; Sandberg J. W. ; Brannigan G. ; Hénin J. Computing absolute binding affinities by Streamlined Alchemical Free Energy Perturbation. bioRxiv Preprint, 2022. 10.1101/2022.12.09.519809. Miyamoto S. ; Kollman P. A. Absolute and Relative Binding Free Energy Calculations of the Interaction of Biotin and Its Analogs with Streptavidin Using Molecular Dynamics/Free Energy Perturbation Approaches. Proteins 1993, 16 , 226–245. 10.1002/prot.340160303.8346190 Gallicchio E. ; Levy R. M. Advances in Protein Chemistry and Structural Biology; Elsevier, 2011; Vol. 85 ; pp 27–80.21920321 Mendoza-Martinez C. ; Papadourakis M. ; Llabrés S. ; Gupta A. A. ; Barlow P. N. ; Michel J. Energetics of a Protein Disorder–Order Transition in Small Molecule Recognition. Chem. Sci. 2022, 13 , 5220–5229. 10.1039/D2SC00028H.35655546 Qian Y. ; Cabeza de Vaca I. ; Vilseck J. Z. ; Cole D. J. ; Tirado-Rives J. ; Jorgensen W. L. Absolute Free Energy of Binding Calculations for Macrophage Migration Inhibitory Factor in Complex with a Druglike Inhibitor. J. Phys. Chem. B 2019, 123 , 8675–8685. 10.1021/acs.jpcb.9b07588.31553604 The PyMOL Molecular Graphics System Version 2.5.4; Schrödinger, LLC, 2022. Duboué-Dijon E. ; Hénin J. Building Intuition for Binding Free Energy Calculations: Bound State Definition, Restraints, and Symmetry. J. Chem. Phys. 2021, 154 , 204101 10.1063/5.0046853.34241173 Mobley D. L. ; Chodera J. D. ; Dill K. A. On the Use of Orientational Restraints and Symmetry Corrections in Alchemical Free Energy Calculations. J. Chem. Phys. 2006, 125 , 084902 10.1063/1.2221683.16965052 Gilson M. ; Given J. ; Bush B. ; McCammon J. The Statistical-Thermodynamic Basis for Computation of Binding Affinities: A Critical Review. Biophys. J. 1997, 72 , 1047–1069. 10.1016/S0006-3495(97)78756-3.9138555 Fu H. ; Chen H. ; Blazhynska M. ; Goulard Coderc de Lacam E. ; Szczepaniak F. ; Pavlova A. ; Shao X. ; Gumbart J. C. ; Dehez F. ; Roux B. ; Cai W. ; Chipot C. Accurate Determination of Protein:Ligand Standard Binding Free Energies from Molecular Dynamics Simulations. Nat. Protoc. 2022, 17 , 1114–1141. 10.1038/s41596-021-00676-1.35277695 Alibay I. ; Magarkar A. ; Seeliger D. ; Biggin P. C. Evaluating the Use of Absolute Binding Free Energy in the Fragment Optimisation Process. Commun. Chem. 2022, 5 , 105 10.1038/s42004-022-00721-4.36697714 Alibay I. IAlibay/MDRestraintsGenerator: MDRestraintsGenerator 0.1.0. Zenodo, 2021. 10.5281/zenodo.4570556. Huggins D. J. Comparing the Performance of Different AMBER Protein Forcefields, Partial Charge Assignments, and Water Models for Absolute Binding Free Energy Calculations. J. Chem. Theory Comput. 2022, 18 , 2616–2630. 10.1021/acs.jctc.1c01208.35266690 Bosisio S. ; Mey A. S. J. S. ; Michel J. Blinded Predictions of Host-Guest Standard Free Energies of Binding in the SAMPL5 Challenge. J. Comput. Aided Mol. Des. 2017, 31 , 61–70. 10.1007/s10822-016-9933-0.27503495 Papadourakis M. ; Bosisio S. ; Michel J. Blinded predictions of standard binding free energies: lessons learned from the SAMPL6 challenge. J. Comput. Aided Mol. Des. 2018, 32 , 1047–1058. 10.1007/s10822-018-0154-6.30159717 Dziedzic P. ; Cisneros J. A. ; Robertson M. J. ; Hare A. A. ; Danford N. E. ; Baxter R. H. G. ; Jorgensen W. L. Design, Synthesis, and Protein Crystallography of Biaryltriazoles as Potent Tautomerase Inhibitors of Macrophage Migration Inhibitory Factor. J. Am. Chem. Soc. 2015, 137 , 2996–3003. 10.1021/ja512112j.25697265 Case D. A. ; Aktulga H. M. ; Belfon K. ; Ben-Shalom I. Y. ; Brozell S. R. ; Cerutti D. S. ; Cheatham T. E. III ; Cruzeiro V. W. D. ; Darden T. A. ; Duke R. E. ; Giambasu G. ; Gilson M. K. ; Gohlke H. ; Goetz A. W. ; Harris R. ; Izadi S. ; Izmailov S. A. ; Jin C. ; Kasavajhala K. ; Kaymak M. C. ; King E. ; Kovalenko A. ; Kurtzman T. ; Lee T. S. ; LeGrand S. ; Li P. ; Lin C. ; Liu J. ; Luchko T. ; Luo R. ; Machado M. ; Man V. ; Manathunga M. ; Merz K. M. ; Miao Y. ; Mikhailovskii O. ; Monard G. ; Nguyen H. ; O’Hearn K. A. ; Onufriev A. ; Pan F. ; Pantano S. ; Qi R. ; Rahnamoun A. ; Roe D. R. ; Roitberg A. ; Sagui C. ; Schott-Verdugo S. ; Shen J. ; Simmerling C. L. ; Skrynnikov N. R. ; Smith J. ; Swails J. ; Walker R. C. ; Wang J. ; Wei H. ; Wolf R. M. ; Wu X. ; Xue Y. ; York D. M. ; Zhao S. ; Kollman P. A. Amber 2021; University of California: San Francisco,2021. Olsson M. H. M. ; Søndergaard C. R. ; Rostkowski M. ; Jensen J. H. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical p Ka Predictions. J. Chem. Theory Comput. 2011, 7 , 525–537. 10.1021/ct100578z.26596171 Søndergaard C. R. ; Olsson M. H. M. ; Rostkowski M. ; Jensen J. H. Improved Treatment of Ligands and Coupling Effects in Empirical Calculation and Rationalization of p Ka Values. J. Chem. Theory Comput. 2011, 7 , 2284–2295. 10.1021/ct200133y.26606496 Anandakrishnan R. ; Aguilar B. ; Onufriev A. V. H++ 3.0: Automating pK Prediction and the Preparation of Biomolecular Structures for Atomistic Molecular Modeling and Simulations. Nucleic Acids Res. 2012, 40 , W537–W541. 10.1093/nar/gks375.22570416 Swope M. Direct Link between Cytokine Activity and a Catalytic Site for Macrophage Migration Inhibitory Factor. EMBO J. 1998, 17 , 3534–3541. 10.1093/emboj/17.13.3534.9649424 Maier J. A. ; Martinez C. ; Kasavajhala K. ; Wickstrom L. ; Hauser K. E. ; Simmerling C. ff14SB: Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB. J. Chem. Theory Comput. 2015, 11 , 3696–3713. 10.1021/acs.jctc.5b00255.26574453 Jakalian A. ; Bush B. L. ; Jack D. B. ; Bayly C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: I. Method. J. Comput. Chem. 2000, 21 , 132–146. 10.1002/(SICI)1096-987X(20000130)21:2<132::AID-JCC5>3.0.CO;2-P. Jakalian A. ; Jack D. B. ; Bayly C. I. Fast, Efficient Generation of High-Quality Atomic Charges. AM1-BCC Model: II. Parameterization and Validation. J. Comput. Chem. 2002, 23 , 1623–1641. 10.1002/jcc.10128.12395429 O’Boyle N. M. ; Banck M. ; James C. A. ; Morley C. ; Vandermeersch T. ; Hutchison G. R. Open Babel: An Open Chemical Toolbox. J. Cheminform. 2011, 3 , 33 10.1186/1758-2946-3-33.21982300 Hedges L. ; Mey A. ; Laughton C. ; Gervasio F. ; Mulholland A. ; Woods C. ; Michel J. BioSimSpace: An Interoperable Python Framework for Biomolecular Simulation. J. Open Source Softw. 2019, 4 , 1831 10.21105/joss.01831. Shirts M. R. ; Mobley D. L. ; Chodera J. D. ; Pande V. S. Accurate and efficient corrections for missing dispersion interactions in molecular simulations. J. Phys. Chem. B 2007, 111 , 13052–13063. 10.1021/jp0735987.17949030 Zacharias M. ; Straatsma T. P. ; McCammon J. A. Separation-shifted Scaling, a New Scaling Method for Lennard-Jones Interactions in Thermodynamic Integration. J. Chem. Phys. 1994, 100 , 9025–9031. 10.1063/1.466707. Michel J. ; Verdonk M. L. ; Essex J. W. Protein-Ligand Complexes: Computation of the Relative Free Energy of Different Scaffolds and Binding Modes. J. Chem. Theory Comput. 2007, 3 , 1645–1655. 10.1021/ct700081t.26627610 Woods C. J. ; Mey A. S. J. S. ; Calabrò G. ; Michel J. Sire Molecular Simulation Framework, 2019. https://github.com/OpenBioSim/sire. Jorgensen W. L. ; Chandrasekhar J. ; Madura J. D. ; Impey R. W. ; Klein M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79 , 926–935. 10.1063/1.445869. Berendsen H. J. C. ; Postma J. P. M. ; van Gunsteren W. F. ; DiNola A. ; Haak J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81 , 3684–3690. 10.1063/1.448118. Calabrò G. ; Woods C. J. ; Powlesland F. ; Mey A. S. J. S. ; Mulholland A. J. ; Michel J. Elucidation of Nonadditive Effects in Protein–Ligand Binding Energies: Thrombin as a Case Study. J. Phys. Chem. B 2016, 120 , 5340–5350. 10.1021/acs.jpcb.6b03296.27248478 Andersen H. C. Molecular Dynamics Simulations at Constant Pressure and/or Temperature. J. Chem. Phys. 1980, 72 , 2384–2393. 10.1063/1.439486. Åqvist J. ; Wennerström P. ; Nervall M. ; Bjelic S. ; Brandsdal B. O. Molecular Dynamics Simulations of Water and Biomolecules with a Monte Carlo Constant Pressure Algorithm. Chem. Phys. Lett. 2004, 384 , 288–294. 10.1016/j.cplett.2003.12.039. Hockney R. ; Goel S. ; Eastwood J. Quiet high-resolution computer models of a plasma. J. of Comput. Phys. 1974, 14 , 148–158. 10.1016/0021-9991(74)90010-2. Tironi I. G. ; Sperb R. ; Smith P. E. ; van Gunsteren W. F. A Generalized Reaction Field Method for Molecular Dynamics Simulations. J. Chem. Phys. 1995, 102 , 5451–5459. 10.1063/1.469273. Shirts M. R. ; Chodera J. D. Statistically Optimal Analysis of Samples from Multiple Equilibrium States. J. Chem. Phys. 2008, 129 , 124105 10.1063/1.2978177.19045004 Bennett C. H. Efficient Estimation of Free Energy Differences from Monte Carlo Data. J. Comput. Phys. 1976, 22 , 245–268. 10.1016/0021-9991(76)90078-4. Deng Y. ; Roux B. Calculation of Standard Binding Free Energies: Aromatic Molecules in the T4 Lysozyme L99A Mutant. J. Chem. Theory Comput. 2006, 2 , 1255–1273. 10.1021/ct060037v.26626834 Orita M. ; Yamamoto S. ; Katayama N. ; Aoki M. ; Takayama K. ; Yamagiwa Y. ; Seki N. ; Suzuki H. ; Kurihara H. ; Sakashita H. ; Takeuchi M. ; Fujita S. ; Yamada T. ; Tanaka A. Coumarin and Chromen-4-one Analogues as Tautomerase Inhibitors of Macrophage Migration Inhibitory Factor: Discovery and X-ray Crystallography. J. Med. Chem. 2001, 44 , 540–547. 10.1021/jm000386o.11170644 Rogers K. E. ; Ortiz-Sánchez J. M. ; Baron R. ; Fajer M. ; de Oliveira C. A. F. ; McCammon J. A. On the Role of Dewetting Transitions in Host–Guest Binding Free Energy Calculations. J. Chem. Theory Comput. 2013, 9 , 46–53. 10.1021/ct300515n.23316123 Bhati A. P. ; Coveney P. V. Large Scale Study of Ligand–Protein Relative Binding Free Energy Calculations: Actionable Predictions from Statistically Robust Protocols. J. Chem. Theory Comput. 2022, 18 , 2687–2702. 10.1021/acs.jctc.1c01288.35293737 Ge Y. ; Wych D. C. ; Samways M. L. ; Wall M. E. ; Essex J. W. ; Mobley D. L. Enhancing Sampling of Water Rehydration on Ligand Binding: A Comparison of Techniques. J. Chem. Theory Comput. 2022, 18 , 1359–1381. 10.1021/acs.jctc.1c00590.35148093 Ben-Shalom I. Y. ; Lin Z. ; Radak B. K. ; Lin C. ; Sherman W. ; Gilson M. K. Accounting for the central role of interfacial water in protein–ligand binding free energy calculations. J. Chem. Theory Comput. 2020, 16 , 7883–7894. 10.1021/acs.jctc.0c00785.33206520 Ben-Shalom I. Y. ; Lin C. ; Kurtzman T. ; Walker R. C. ; Gilson M. K. Simulating water exchange to buried binding sites. J. Chem. Theory Comput. 2019, 15 , 2684–2691. 10.1021/acs.jctc.8b01284.30835999 Ross G. A. ; Russell E. ; Deng Y. ; Lu C. ; Harder E. D. ; Abel R. ; Wang L. Enhancing water sampling in free energy calculations with grand canonical Monte Carlo. J. Chem. Theory Comput. 2020, 16 , 6061–6076. 10.1021/acs.jctc.0c00660.32955877 Cisneros J. A. ; Robertson M. J. ; Valhondo M. ; Jorgensen W. L. A Fluorescence Polarization Assay for Binding to Macrophage Migration Inhibitory Factor and Crystal Structures for Complexes of Two Potent Inhibitors. J. Am. Chem. Soc. 2016, 138 , 8630–8638. 10.1021/jacs.6b04910.27299179 Hahn D. F. ; König G. ; Hünenberger P. H. Overcoming Orthogonal Barriers in Alchemical Free Energy Calculations: On the Relative Merits of łambda-Variations, łambda-Extrapolations, and Biasing. J. Chem. Theory Comput. 2020, 16 , 1630–1645. 10.1021/acs.jctc.9b00853.31995374 Lapelosa M. ; Gallicchio E. ; Levy R. M. Conformational transitions and convergence of absolute binding free energy calculations. J. Chem. Theory Comput. 2012, 8 , 47–60. 10.1021/ct200684b.22368530 Khalak Y. ; Tresadern G. ; Aldeghi M. ; Baumann H. M. ; Mobley D. L. ; de Groot B. L. ; Gapsys V. Alchemical Absolute Protein–Ligand Binding Free Energies for Drug Design. Chem. Sci. 2021, 12 , 13958–13971. 10.1039/D1SC03472C.34760182 Aguayo-Ortiz R. ; Dominguez L. Unveiling the Possible Oryzalin-Binding Site in the α-Tubulin of ToxoplasmaGondii. ACS Omega 2022, 7 , 18434–18442. 10.1021/acsomega.2c00729.35694483 Woo H.-J. ; Roux B. Calculation of Absolute Protein–Ligand Binding Free Energy from Computer Simulations. Proc. Natl. Acad. Sci. U.S.A. 2005, 102 , 6825–6830. 10.1073/pnas.0409005102.15867154 Wang L. ; Deng Y. ; Wu Y. ; Kim B. ; LeBard D. N. ; Wandschneider D. ; Beachy M. ; Friesner R. A. ; Abel R. Accurate Modeling of Scaffold Hopping Transformations in Drug Discovery. J. Chem. Theory Comput. 2017, 13 , 42–54. 10.1021/acs.jctc.6b00991.27933808