Calculating Binding Free Energies with BEDAM

The BEDAM method^{29 (link)} computes the binding free energy
$\mathrm{\Delta }{F}_{AB}^{\circ }$ for the monovalent association of a receptor A and a ligand B by means of the expression
$$\mathrm{\Delta }{F}_{AB}^{\circ }=-\mathit{kT}ln\left[C^{\circ }{V}_{\text{site}}\int du{p}_0(u)e^{-\beta u}\right]=-\mathit{kT}ln{C}^{\circ }{V}_{\text{site}}+\mathrm{\Delta }{F}_{AB},$$ which follows, without approximations, from a well-established statistical mechanics theory of molecular association,^{11 (link)} where β = 1/kT, C∘ is the standard concentration of ligand molecules (set to C∘ = 1 M, or equivalently 1,668 Å⁻³), V_site is the volume of the binding site, and p₀(u) is the probability distribution of binding energies collected in an appropriate decoupled ensemble of conformations in which the ligand is confined in the binding site while the receptor and the ligand are both interacting only with the solvent continuum and not with each other. The binding energy
$$u({\mathbf{r}}_B,{\mathbf{r}}_A)=V({\mathbf{r}}_B,{\mathbf{r}}_A)-V({\mathbf{r}}_B)-V({\mathbf{r}}_A)$$ is defined for each conformation r = (r_B, r_A) of the complex as the difference between the effective potential energies V (r) of the associated and separated conformations of the complex without conformational rearrangements. In our implementation BEDAM employs an effective potential in which the solvent is represented implicitly by means of the AGBNP2 implicit solvent model^{37 (link)} together with the OPLS-AA^{51 ,52} force field for covalent and non-bonded interatomic interactions.
Eq. (1) explicitly indicates that the standard binding free energy is the sum of two terms. The first term, −kT lnC∘V_site, represents the entropic work to transfer the ligand from the solution environment at concentration C∘ into the binding site of the complex. This term depends only on the standard state and the definition of the complex macrostate. The second term, ΔF_AB, involving the Boltzmann-weighted integral of p₀(u), corresponds to the work for turning on the interactions between the receptor and the ligand when the ligand is confined in the binding site region.^{29 (link)}Eq. (1) also naturally leads to the the definition of a binding affinity density function k(u) = C∘V_sitep₀(u) exp(−βu) in terms of which the binding constant is written as
Based on Eq. (3) the binding affinity density k(u) can be interpreted as a measure of the contribution of the conformations of the complex with the binding energy, u, to the binding constant.^{29 (link)} We have shown that k(u) is proportional to p₁(u), the binding affinity density in the coupled ensemble of the complex, with a proportionality constant related to the binding free energy.^{29 (link)}The larger the value of the integral in Eq. (1), the more favorable is the binding free energy. The magnitude of the p₀(u) distribution at positive, unfavorable, values of the binding energy u reflects the entropic thermodynamic driving force which opposes binding, whereas the tail at negative, favorable, binding energies measures the energetic gain for binding due to the formation of ligand-receptor interactions. The interplay between these two opposing forces ultimately determines the strength of binding. We found that the ability of BEDAM to explicitly include both favorable energetic gains and unfavorable entropic losses to be essential to properly reproduce experimental binding affinities in a challenging set of candidate ligands to T4 lysozyme receptors whose estimates of binding affinity failed by simplified docking & scoring approaches.^{53 (link)}The accurate calculation of the important low energy tail of p₀(u) can not be accomplished by a brute-force collection of binding energy values from a simulation of the complex in the decoupled state because these are rarely sampled when the ligand is not guided by the interactions with the receptor. Instead we use biased sampling and parallel Hamiltonian Replica Exchange (HREM) in which swarms of coupled replicas of the system, differing in the value of an interaction parameter 0 ≤ λ ≤ 1 controlling the strength of ligand-receptor interactions, are simulated simultaneously. The replicas collectively sample a wide range of unfavorable, intermediate and favorable binding energies which are then unbiased and combined together by means of reweighting techniques.^{34 (link),35 (link)}BEDAM is based on biasing potentials of the form λ u(r) yielding a family of λ -dependent hybrid potentials of the form
$$V_{\lambda }(\mathbf{r})=V_0(\mathbf{r})+\lambda u(\mathbf{r}),$$ where
$$V_0(\mathbf{r})=V({\mathbf{r}}_B)+V({\mathbf{r}}_A),$$ is the potential energy function of the decoupled state. It is easy to see from Eqs. (2), (4), and (5) that V_λ=1 corresponds to the effective potential energy of the coupled complex and V_λ=0 corresponds to the state in which the receptor and ligand are not interacting (decoupled state). Intermediate values of λ trace an alchemical thermodynamic path connecting these two states. The binding free energy ΔF_AB is by definition the difference in free energy between these two states.
For later use we introduce here the reorganization free energy for binding
$\mathrm{\Delta }{G}_{\mathit{reorg}}^{\circ }$ defined by the expression¹⁰ $$\mathrm{\Delta }{F}_{AB}^{\circ }={\langle u\rangle }_1+\mathrm{\Delta }{G}_{\mathit{reorg}}^{\circ }$$ where 〈u〉₁ is the average binding energy at λ = 1 and
$\mathrm{\Delta }{F}_{AB}^{\circ }$ is the standard binding free energy.