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.