Accuracy of continuum electrostatic calculations based on three common dielectric boundary definitions

Abstract

We investigate the influence of three common definitions of the solute/solvent dielectric boundary (DB) on the accuracy of the electrostatic solvation energy ΔGel computed within the Poisson Boltzmann and the generalized Born models of implicit solvation. The test structures include small molecules, peptides and small proteins; explicit solvent ΔGel are used as accuracy reference. For common atomic radii sets BONDI, PARSE (and ZAP9 for small molecules) the use of van der Waals (vdW) DB results, on average, in considerably larger errors in ΔGel than the molecular surface (MS) DB. The optimal probe radius ρw for which the MS DB yields the most accurate ΔGel varies considerably between structure types. The solvent accessible surface (SAS) DB becomes optimal at ρw ~ 0.2 Å (exact value is sensitive to the structure and atomic radii), at which point the average accuracy of ΔGel is comparable to that of the MS-based boundary. The geometric equivalence of SAS to vdW surface based on the same atomic radii uniformly increased by ρw gives the corresponding optimal vdW DB. For small molecules, the optimal vdW DB based on BONDI + 0.2 Å radii can yield ΔGel estimates at least as accurate as those based on the optimal MS DB. Also, in small molecules, pairwise charge-charge interactions computed with the optimal vdW DB are virtually equal to those computed with the MS DB, suggesting that in this case the two boundaries are practically equivalent by the electrostatic energy criteria. In structures other than small molecules, the optimal vdW and MS dielectric boundaries are not equivalent: the respective pairwise electrostatic interactions in the presence of solvent can differ by up to 5 kcal/mol for individual atomic pairs in small proteins, even when the total ΔGel are equal. For small proteins, the average decrease in pairwise electrostatic interactions resulting from the switch from optimal MS to optimal vdW DB definition can be mimicked within the MS DB definition by doubling of the solute dielectric constant. However, the use of the higher interior dielectric does not eliminate the large individual deviations between pairwise interactions computed within the two DB definitions. It is argued that while the MS based definition of the dielectric boundary is more physically correct in some types of practical calculations, the choice is not so clear in some other common scenarios.

Quantum chemistry in arbitrary dielectric environments: Theory and implementation of nonequilibrium Poisson boundary conditions and application to compute vertical ionization energies at the air/water interface

“…The usual PCM prescription using atom-centered vdW spheres, however, proves to be problematic when explicit water molecules are included in the QM region, leading to unphysical high-dielectric regions between these explicit solvent molecules. Oddly, this problem is not often discussed in the quantum chemistry literature although a similar problem in biomolecular simulation has been widely discussed, [132][133][134][135][136] where in the context of Poisson-Boltzmann electrostatics calculations the vdW cavity construction may leave highdielectric regions in the hydrophobic interior of a protein. In this work, we introduced a "hybrid" cavity model that avoids this problem.…”

Accuracy comparison of several common implicit solvent models and their implementations in the context of protein-ligand binding

“…This medium has predetermined electrostatic properties including a specified dielectric constant; the solute is separated from the solvent by a dielectric boundary (DB), the results of practical calculations are very sensitive to the choice of DB [47], with several possible, generally non-equivalent, choices [20]. In this study the PCM, S-GB and COSMO models implemented in the DISOLV program and the PCM model in the MCBHSOLV program employ the same DB (the Solvent Excluded Surface or SES) constructed as follows [10, 26, 27, 30].…”

“…The experimental hydration energies were also known for these molecules and the values are given in [38] and in the respective Supplementary materials. The second set consists of 19 small charge neutral proteins with no more than 500 atoms per structure: [20] (PDB IDs: 1az6, 1bh4, 1bku, 1brv, 1byy, 1cmr, 1dfs, 1dmc, 1eds, 1fct, 1fmh, 1fwo, 1g26, 1ha9, 1hzn, 1paa, 1qfd, 1qk7, 1scy), which were selected from the larger set [21]. The third set used here for the desolvation energy calculations, consists of 15 protein-ligand complexes (PDB IDs: 1b11, 1bkf, 1f40, 1fb7, 1fkb, 1fkf, 1fkg, 1fkh, 1fkj, 1fkl, 1pbk, 1zp8, 2fke, 2hah, 3kfp) [22].…”

“…The choice of the ZAP9 [60, 61] radii set for small molecules is based on the results presented in [14], where it has been shown that radii set ZAP9 more precisely reproduce experimental values of solvation energies for small molecules. In case of proteins and complexes we have taken the radii set mbondi2 because it has also provided good results [20]. …”

“…good correspondence of the calculated hydration energies to experimentally measured ones for a set of more than 400 small molecules [26]. Taking into account that DISOLV and MCBHSOLV programs are used successfully for new inhibitors development [32–34] at the post-processing stage [10] and in the gridless docking procedure [35] in the frame of the MMFF94 force field [31] there is a need to compare results of DISOLV and MCBHSOLV calculations with ones performed with other implicit solvent programs and parameterization methods [14, 20–22, 25] for same sets of molecules. Here we will also examine a relatively recent addition to the family of GB models, GBNSR6, which has already demonstrated high accuracy in estimates of hydration free energies of small molecules [14].…”