Improving Prediction Accuracy of Binding Free Energies and Poses of HIV Integrase Complexes Using the Binding Energy Distribution Analysis Method with Flattening Potentials

Abstract

To accelerate conformation sampling of slow dynamics from receptor or ligand, we introduced flattening potentials on selected bonded and nonbonded intramolecular interactions to the binding energy distribution analysis method (BEDAM) for calculating absolute binding free energies of protein-ligand complexes using an implicit solvent model and implemented flattening BEDAM using the asynchronous replica exchange (AsyncRE) framework for performing large scale replica exchange molecular dynamics (REMD) simulations. The advantage of using the flattening feature to reduce high energy barriers was exhibited first by the p-xylene-T4 lysozyme complex, where the intramolecular interactions of a protein side chain on the binding site were flattened to accelerate the conformational transition of the side chain from the trans to the gauche state when the p-xylene ligand is present in the binding site. Much more extensive flattening BEDAM simulations were performed for 53 experimental binders and 248 nonbinders of HIV-1 integrase which formed the SAMPL4 challenge, with the total simulation time of 24.3 μs. We demonstrated that the flattening BEDAM simulations not only substantially increase the number of true positives (and reduce false negatives) but also improve the prediction accuracy of binding poses of experimental binders. Furthermore, the values of area under the curve (AUC) of receiver operating characteristic (ROC) and the enrichment factors at 20% cutoff calculated from the flattening BEDAM simulations were improved significantly in comparison with that of simulations without flattening as we previously reported for the whole SAMPL4 database. Detailed analysis found that the improved ability to discriminate the binding free energies between the binders and nonbinders is due to the fact that the flattening simulations reduce the reorganization free energy penalties of binders and decrease the overlap of binding free energy distributions of binders relative to that of nonbinders. This happens because the conformational ensemble distributions for both the ligand and protein in solution match those at the fully coupled (complex) state more closely when the systems are more fully sampled after the flattening potentials are applied to the intermediate states.

Figure 1:

Cartoon representation of the thermodynamic cycles of the BEDAM method.

Figure 2:

Cartoon representations of T4 Lysozyme complexes with p-xylene. Green: p-xylene ligand; Red: trans conformation (χ≈ 180° for the side-chain N-CA-CB-CG1 dihedral angle) of Val111; Blue: gauche conformation (χ≈−60°).

Figure 3:

The time series of sidechain dihedral angle (N-CA-CB-CG1) of Val111 of T4 Lysozyme in the complex with p-xylene (at λ_b = 1). a) normal BEDAM simulation without flattening started from a gauche state (χ₁ ≈ −60°); b) normal BEDAM simulation without flattening started from the trans state (χ₁ ≈ 180°); c) BEDAM with the flattening of all torsional (including 1–4) interactions on the sidechain related to the CA-CB bond of the Val111 sidechain (started from the trans state); d) In addition to the flattened interactions in c), nonbonded interactions such as L-J and Coulomb interactions from all atoms involved in c) are also flattened (started from the trans state). Note both conformations of χ₁ ≈ −60° and +60° belong to the gauche state due to the symmetry of Val111 sidechain. The binding free energies calculated from these four different BEDAM simulations are listed in Table 2.

Figure 4:

The binding energies at λ_b = 1 (fully coupled state) calculated from the BEDAM simulations with three different flattening strategies as a function of ligand RMSD from the crystal pose for two HIV-IN complexes (a) avx38747 and (c) avx38788 from SAMPL4. (b) and (d) show the corresponding 2D schematic graphs of ligands with the rotatable bonds shown as in red arrows for applying flattening potentials. The black squares show the starting conformations from the AutoDock program.

Figure 5:

The binding energies at λ_b = 1 (fully coupled state) obtained from three different types of BEDAM simulations as a function of ligand RMSD from the crystal pose for two HIV-IN complexes avx38747 (a and b) and avx38788 (c and d) from SAMPL4. For (a) and (c), the red points show the 1000 snapshots from the BEDAM simulations without flattening and started from the crystal poses; the green points denote the 1000 snapshots from the similar BEDAM simulations as the red points but started from the docked poses; the blue points display the 1000 snapshots from the flattening BEDAM simulations started from the same docked poses as that of the green points. In (b) and (d), the ligands in red show the crystal poses, those in green exhibit the docking poses predicted by AutoDock, and those in blue display the best poses from the flattening BEDAM simulations that consistent with the corresponding crystal poses.

Figure 6:

The boxplots of the BEDAM simulation results for the SAMPL4 53 experimental binders from different types of runs. a) binding free energies; b) minimum values of binding energies from the trajectory snapshots at λ_b = 1; c) minimum values of ligand RMSDs from the corresponding crystal poses calculated from the trajectory snapshots at λ_b = 1. The exact median and mean values are listed in Table 4.

Figure 7:

The binding free energies in Fig. 6a are decomposed into two parts. a) $\Delta G_{II}^o=<u{>}_{RL}$, average binding energies at λ_b = 1; b) $\Delta G_I^o=\Delta G_{reorg}^o$, reorganization free energies.

Figure 8:

The boxplots of the BEDAM simulation results for the whole set of SAMPL4 ligand from different types of runs. a) binding free energies $\left(\Delta G_b^o\right)$; b) average binding energies $\left(\Delta G_{II}^o=<u{>}_{RL}\right)$ at λ_b = 1; c) reorganization free energies $\left(\Delta G_I^o=\Delta G_{reorg}^o\right)$

Figure 9:

ROC curves of different types of BEDAM simulations for the whole database of SAMPL4 using two different predictions: a) binding free energies $\left(\Delta G_b^o\right)$ and b) average binding energies $\left(\Delta G_{II}^o=<u{>}_{RL}\right)$.