Multidimensional Global Optimization and Robustness Analysis in the Context of Protein–Ligand Binding

Abstract

Accuracy of protein–ligand binding free energy calculations utilizing implicit solvent models is critically affected by parameters of the underlying dielectric boundary, specifically, the atomic and water probe radii. Here, a global multidimensional optimization pipeline is developed to find optimal atomic radii specifically for protein–ligand binding calculations in implicit solvent. The computational pipeline has these three key components: (1) a massively parallel implementation of a deterministic global optimization algorithm (VTDIRECT95), (2) an accurate yet reasonably fast generalized Born implicit solvent model (GBNSR6), and (3) a novel robustness metric that helps distinguish between nearly degenerate local minima via a postprocessing step of the optimization. A graph-based “kT-connectivity” approach to explore and visualize the multidimensional energy landscape is proposed: local minima that can be reached from the global minimum without exceeding a given energy threshold (kT) are considered to be connected. As an illustration of the capabilities of the optimization pipeline, we apply it to find a global optimum in the space of just five radii: four atomic (O, H, N, and C) radii and water probe radius. The optimized radii, ρW = 1.37 Å, ρC = 1.40 Å, ρH = 1.55 Å, ρN = 2.35 Å, and ρO = 1.28 Å, lead to a closer agreement of electrostatic binding free energies with the explicit solvent reference than two commonly used sets of radii previously optimized for small molecules. At the same time, the ability of the optimizer to find the global optimum reveals fundamental limits of the common two-dielectric implicit solvation model: the computed electrostatic binding free energies are still almost 4 kcal/mol away from the explicit solvent reference. The proposed computational approach opens the possibility to further improve the accuracy of practical computational protocols for binding free energy calculations.

An Effective MM/GBSA Protocol for Absolute Binding Free Energy Calculations: A Case Study on SARS-CoV-2 Spike Protein and the Human ACE2 Receptor

“…Therefore, the right choice of atomic radii is crucial to the accuracy of binding free energy estimation [64,65]. Two sets of atomic radii are used here: global optimal radii for ∆G bind calculations (OPT1) [43,44] and Bondi [66], see Table 1. The van der Waals radii determined by Bondi from molecular crystals and noble gas crystals are commonly known as "all-purpose" sets of intrinsic atomic radii in a wide range of molecular modeling applications [67] including continuum solvent calculations [68].…”

“…Through the MM/GBSA approach, the absolute binding free energy of a complex is calculated as the sum of gas-phase energy, solvation free energy, and entropic contributions averaged over several snapshots extracted from the main MD trajectory. A grid-based surface GB model is used for estimating the polar component of solvation-free energy, coupled with a new water model and atomic radii introduced earlier in [43,44]. Human H-Ras and the Ras-binding domain of C-Raf1, the so-called Ras-Raf complex [8,28], is chosen as the reference for the initial evaluation of the MM/GBSA model.…”

MMGB/SA Consensus Estimate of the Binding Free Energy Between the Novel Coronavirus Spike Protein to the Human ACE2 Receptor

“…Therefore, the right choice of atomic radii is crucial to the accuracy of binding free energy estimation. Three sets of atomic radii are used here: OPT1, 29 , 46 bondi, and mbondi2. The first two are listed in Tab.…”

“…A grid-based surface GB model is used for estimating the polar component of solvation free energy, coupled with water and atomic radii introduced earlier. 29 Human H-Ras and the Ras-binding domain of C-Raf1, so-called Ras-Raf complex, 6 , 30 is chosen as the reference for the initial evaluation of the MMGB/SA model. Final results are compared with those from a few available relevant studies.…”

“…Also for consistency with the same work, here we use γ = 0.0072 kcal/mol/A 2 . Atomic radii that form SASA not only play an important role in the non-polar component, but also enter the polar component OPT1,29,46 bondi, and mbondi2. The first two are listed in Tab.…”