Estimating the Roles of Protonation and Electronic Polarization in Absolute Binding Affinity Simulations

Abstract

Accurate prediction of binding free energies is critical to streamlining the drug development and protein design process. With the advent of GPU acceleration, absolute alchemical methods, which simulate the removal of ligand electrostatics and van der Waals interactions with the protein, have become routinely accessible and provide a physically rigorous approach that enables full consideration of flexibility and solvent interaction. However, standard explicit solvent simulations are unable to model protonation or electronic polarization changes upon ligand transfer from water to the protein interior, leading to inaccurate prediction of binding affinities for charged molecules. Here, we perform extensive simulation totaling ∼540 μs to benchmark the impact of modeling conditions on predictive accuracy for absolute alchemical simulations. Binding to urokinase plasminogen activator (UPA), a protein frequently overexpressed in metastatic tumors, is evaluated for a set of 10 inhibitors with extended flexibility, highly charged character, and titratable properties. We demonstrate that the alchemical simulations can be adapted to utilize the MBAR/PBSA method to improve the accuracy upon incorporating electronic polarization, highlighting the importance of polarization in alchemical simulations of binding affinities. Comparison of binding energy prediction at various protonation states indicates that proper electrostatic setup is also crucial in binding affinity prediction of charged systems, prompting us to propose an alternative binding mode with protonated ligand phenol and Hid-46 at the binding site, a testable hypothesis for future experimental validation.

Figure 1.

Chemical structures of the 10 evaluated UPA inhibitors. The molecules share a benzamidine-like scaffold with characteristic amidine group carrying positive charge, and extended tails comprised of a phenol group and other functional modifiers. The hydroxyl on the phenol is proposed to be titratable and samples deprotonated and protonated states during binding, altering the hydrogen bonding capability of the ligands. The inhibitors are categorized as small (those without the phenol group)—1C5X, 1C5Y, 1C5Z, and 1GI7—and big (for those with potentially charged phenols)—1GJ7, 1GJ8, 1GJA, 1GJB, 1GJD, 1O3P.

Figure 2.

Example inhibitor binding poses. (A) The protein and ligand form a network of polar interactions at two locations, at the base of the active site between the negatively charged Asp-192 and the positively charged amidine, and near the phenol hydroxyl with Ser-198 and His-46. (B) Electron density supports the positioning of the ligand hydroxyl unusually close to Ser-198. An exceptionally short hydrogen bond is formed between the phenol hydroxyl and Ser-198 hydroxyl with a distance of ~2.2 Å; this interaction may not be captured accordingly with typical force fields due to van der Waals repulsion.

Figure 3.

Relieving steric clash between the ligand phenol and Ser-198. (A) The distance between the ligand phenol oxygen and Ser-198 hydroxyl oxygen is recorded over the last 10 ns of equilibration to analyze sampled conformations and compared to the distance observed in the crystal structures. The trend observed is identical for both GAFF and GAFF2 force fields, all ligands except 1GJD twist away due to repulsive steric interactions but remain in hydrogen bonding range. 1GJD samples broad distances, indicating the initial hydrogen bond is detached. (B) Sample frame from the 1GJD simulation illustrates that the phenol hydroxyl rotates outward away from the protein, and the starting hydrogen bond is replaced with one between the peptide bond-like carbonyl and Ser-198. The inhibitor is colored green and labeled DRG.

Figure 4.

Baseline absolute alchemical binding predictions for UPA inhibitors. Evaluating the effects of simulation with 150 mM salt alone, deprotonated ligands alone, and with 150 mM salt and deprotonated ligands combined (1DOF All-HIP) on the baseline condition (fully protonated ligands, counterions added only up to neutralize system charge, and 1DOF restraints). The highest performance is observed with the 1DOF All-HIP condition with RMSE 2.50 kcal/mol and Pearson correlation 0.51.

Figure 5.

MBAR/PBSA binding affinity calculations. The All-HID condition shows the best agreement to experiment with consideration of polarization effects through solute dielectric scaling. In comparison to values from the standard alchemical transformation, RMSEs are reduced and Pearson correlations are improved for all conditions.