Tackling Halogenated Species with PBSA: Effect of Emulating the σ-Hole

To model halogen bond phenomena using classical force fields, an extra-point (EP) of charge is frequently introduced at a given distance from the halogen (X) to emulate the σ-hole. The resulting molecular dynamics (MD) trajectories can be used in subsequent molecular mechanics (MM) combined with Poisson–Boltzmann and surface area calculations (MM PBSA) to estimate protein–ligand binding free energies (∆G_bind). While EP addition improves the MM/MD description of halogen-containing systems, its effect on the calculation of solvation free energies (∆G_solv) using the PBSA approach is yet to be assessed. As the PBSA calculations depend, among other parameters, on the empirical assignment of radii (PB radii), a problematic issue arises since standard halogen radii are smaller than the typical X· · · EP distances (usually corresponding to R_min), thus placing the EP within the solvent dielectric. Herein, we performed a comprehensive study on the performance of PBSA (using three different setups) in the calculation of ∆Gsolv values for 142 halogenated compounds (bearing Cl, Br, or I) for which the experimental values are known. By conducting an optimization (minimizing the error against experimental values), we provide a new optimized set of halogen PB radii, for each PBSA setup, that should be used when the EP is located at R min in the context of GAFF. A simultaneous optimization of PB radii and X· · · EP distances shows that a wide range of distance/radius pairs can be used without significant loss of accuracy, therefore laying the basis for expanding this halogen radii optimization strategy to other force fields and EP implementations. As ligand ∆G_solv estimation is an important term in the determination of protein–ligand ∆G_bind , this work is particularly relevant in the framework of structure-based virtual screening and related computer-aided drug design routines.