Computations of standard binding free energies with molecular dynamics simulations

Abstract

An increasing number of studies have reported computations of the standard (absolute) binding free energy of small ligands to proteins using molecular dynamics (MD) simulations and explicit solvent molecules that are in good agreement with experiments. This encouraging progress suggests that physics-based approaches hold the promise of making important contributions to the process of drug discovery and optimization in the near future. Two types of approaches are principally used to compute binding free energies with MD simulations. The most widely known is the alchemical double decoupling method, in which the interaction of the ligand with its surroundings are progressively switched off. It is also possible to use a potential of mean force (PMF) method, in which the ligand is physically separated from the protein receptor. For both of these computational approaches, restraining potentials may be activated and released during the simulation for sampling efficiently the changes in translational, rotational, and conformational freedom of the ligand and protein upon binding. Because such restraining potentials add bias to the simulations, it is important that their effects be rigorously removed to yield a binding free energy that is properly unbiased with respect to the standard state. A review of recent results is presented, and differences in computational methods are discussed. Examples of computations with T4-lysozyme mutants, FKBP12, SH2 domain, and cytochrome P450 are discussed and compared. Remaining difficulties and challenges are highlighted.

Figure 1

Shown in the figure is path for the free energy computation accordinng to the alchemical DDM with restraining potentials. $L_{\mathrm{aq}}^{\mathrm{c}}$ denotes the configuration restricted ligand in bulk solvent. $L_{\mathrm{vac}}^{\mathrm{c}}$ denotes the non-interacting ligand with configurational restraint. $L_{\mathrm{vac}}^{\mathrm{c},\mathrm{r}}$ denotes the non-interacting ligand with both configurational and rotational restraints. The spring represents the translational restraint potential. See reference 34 for a complete theoretical formulation.

Figure 2

Shown in the figure is path for the free energy computation using a PMF-based method with restraining potentials. $L_{\mathrm{aq}}^{\mathrm{c}}$ denotes the configuration restricted ligand in bulk solvent. $L_{\mathrm{aq}}^{\mathrm{c},\mathrm{o}}$ denotes the ligand with both configuration and orientation restraints. $L_{\mathrm{aq}}^{\mathrm{c},\mathrm{o},\mathrm{a}}$ denotes the ligand in the bulk with all configurational, orientational and axial restraint potentials. $\Delta G_{\mathrm{a}}^{\text{bulk}}$ is related to S* in Eq (9). ΔG^bulk→site is related to I* in Eq (9). See reference 49 for a complete theoretical formulation.

Figure 3

T4 lysozyme L99A mutant with benzene bound in the cavity. The grey parts are treated as a mean field approximation with generalized solvent boundary potential. See reference 34 for computational details.

Figure 4

FKBP12 bound with ligand #8 studied previously.^, The grey parts are treated as a mean field approximation with generalized solvent boundary potential. See reference 40 for computational details.

Figure 5

SH2 domain with bound peptide pYEEI. The system was simulated with PBC and water is not shown for clarity. See reference 49 for computational details.

Figure 6

The binding site of cytochrome p450 binding site. In A is an overview of the simulation system with the burried binding site. In B, the binding site is not visible with a space-filling representation. In C, the site with no camphor isoccupied by water, in D camphor is bound. See reference 78 for theoretical formulation and computational details.

Figure 7

The jnk kinase bound with ligand 19 from the OpenEye statistical assessment of the modeling of proteins and ligands (SAMPL) challenge. The grey parts are treated as a mean field approximation with generalized solvent boundary potential. The computations were carried out according to the same protocol presented in previous studies.^,,

Figure 8

The cumulative experimental rank (vertical axis) is plotted as a function of the calculated (horizontal axis). The top five ligands from the computations are: 40, 27, 10, 47, 8 ; the experimental top five are: 10, 1, 27, 11, 38. The worse ten ligands from the computations are: 9, 43, 51, 6, 5, 21, 23, 56, 31, 17; the experimental worst ten are: 51, 4, 5, 6, 17, 23, 29, 30, 36, 43.