The present paper is intended to be a comprehensive assessment and rationalization, from a statistical mechanics perspective, of existing alchemical theories for binding free energy calculations of ligand−receptor systems. In detail, the statistical mechanics foundation of noncovalent interactions in ligand−receptor systems is revisited, providing a unifying treatment that encompasses the most important variants in the alchemical approaches from the seminal double annihilation method [Jorgensen et al. J. Chem. Phys. 1988; 89, 3742]
■ INTRODUCTIONThe determination of the binding free energy in ligand− receptor systems is the cornerstone of drug discovery. In the last few decades, traditional molecular docking techniques in computer-assisted drug design have been modified, integrated, or superseded using methodologies relying on a more realistic description of the drug−receptor system. It has become increasingly clear that a microscopic description of the solvent is a crucial ingredient to reliably rank the affinity of putative ligands for a given target. 1,2 In the framework of atomistic molecular dynamics (MD) simulations with explicit solvent, several methods for rigorously determining the absolute binding free energy in noncovalently bonded systems have been devised. Most of these methodologies are based on the socalled alchemical route 3−5 (see refs 6−9 for recent reviews). In this approach, proposed for the first time by Jorgensen et al. 10,11 and named double annihilation method (DAM), the absolute binding free energy of the ligand−receptor system was obtained by setting up a thermodynamic cycle as indicated in Figure 1 in which the basic quantities to be estimated are the annihilation or decoupling free energies of the ligand in the bound state and in bulk water, indicated hereafter as ΔG b and ΔG u , respectively. These decoupling free energies correspond to the two closing branches of the cycle and are obtained by discretizing the alchemical path connecting the fully interacting and fully decoupled ligand in a number of intermediate nonphysical states and then running MD simulations for each of these states.Alchemical states are defined by a λ coupling parameter entering the Hamiltonian, varying between 1 and 0, such that at λ = 1 and λ = 0 one has the fully interacting and gas-phase ligand, respectively. ΔG b and ΔG u are usually recovered as a sum of the contributions from each of the λ windows by applying the free energy perturbation method 12 (FEP). Alternatively, and equivalently, one can compute the canonical average of the derivative of the Hamiltonian at the discrete λ points, obtaining the decoupling free energy via numerical thermodynamic integration (TI). 13 Finally, the thermodynamic cycle is closed by computing the difference between the two decoupling free energies along the alchemical path, ΔG b and ΔG u , obtaining the dissociation free energy in solution.Gilson et al. 14 criticized Jorgensen's theory 10 by pointing out that the resulting dissociation free energy does not de...