By using different evaluation strategies, we systemically evaluated the performance of Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) and Molecular Mechanics/Poisson–Boltzmann Surface Area (MM/PBSA) methodologies based on more than 1800 protein–ligand crystal structures in the PDBbind database. The results can be summarized as follows: (1) for the one-protein-family/one-binding-ligand case which represents the unbiased protein–ligand complex sampling, both MM/GBSA and MM/PBSA methodologies achieve approximately equal accuracies at the interior dielectric constant of 4 (with r_p = 0.408 ± 0.006 of MM/GBSA and r_p = 0.388 ± 0.006 of MM/PBSA based on the minimized structures); while for the total dataset (1864 crystal structures), the overall best Pearson correlation coefficient (r_p = 0.579 ± 0.002) based on MM/GBSA is better than that of MM/PBSA (r_p = 0.491 ± 0.003), indicating that biased sampling may significantly affect the accuracy of the predicted result (some protein families contain too many instances and can bias the overall predicted accuracy). Therefore, family based classification is needed to evaluate the two methodologies; (2) the prediction accuracies of MM/GBSA and MM/PBSA for different protein families are quite different with r_p ranging from 0 to 0.9, whereas the correlation and ranking scores (an averaged r_p/r_s over a list of protein folds and also representing the unbiased sampling) given by MM/PBSA (r_p-score = 0.506 ± 0.050 and r_s-score = 0.481 ± 0.052) are comparable to those given by MM/GBSA (r_p-score = 0.516 ± 0.047 and r_s-score = 0.463 ± 0.047) at the fold family level; (3) for the overall prediction accuracies, molecular dynamics (MD) simulation may not be quite necessary for MM/GBSA (r_p-minimized = 0.579 ± 0.002 and r_p-1ns = 0.564 ± 0.002), but is needed for MM/PBSA (r_p-minimized = 0.412 ± 0.003 and r_p-1ns = 0.491 ± 0.003). However, for the individual systems, whether to use MD simulation is depended. (4) both MM/GBSA and MM/PBSA may be unable to give successful predictions for the ligands with high formal charges, with the Pearson correlation coefficient ranging from 0.621 ± 0.003 (neutral ligands) to 0.125 ± 0.142 (ligands with a formal charge of 5). Therefore, it can be summarized that, although MM/GBSA and MM/PBSA perform similarly in the unbiased dataset, for the currently available crystal structures in the PDBbind database, compared with MM/GBSA, which may be used in multi-target comparisons, MM/PBSA is more sensitive to the investigated systems, and may be more suitable for individual-target-level binding free energy ranking. This study may provide useful guidance for the post-processing of docking based studies.

With the rapid development of computational techniques and hardware, more rigorous and precise theoretical models have been used to predict the binding affinities of a large number of small molecules to biomolecules. By employing continuum solvation models, the MM/GBSA and MM/PBSA methodologies achieve a good balance between low computational cost and reasonable prediction accuracy. In this study, we have thoroughly investigated the effects of interior dielectric constant, molecular dynamics (MD) simulations, and the number of top-scored docking poses on the performance of the MM/GBSA and MM/PBSA rescoring of docking poses for three tyrosine kinases, including ABL, ALK, and BRAF. Overall, the MM/PBSA and MM/GBSA rescoring achieved comparative accuracies based on a relatively higher solute (or interior) dielectric constant (i.e. ε = 2, or 4), and could markedly improve the ‘screening power’ and ‘ranking power’ given by Autodock. Moreover, with a relatively higher solute dielectric constant, the MM/PBSA or MM/GBSA rescoring based on the best scored docking poses and the multiple top-scored docking poses gave similar predictions, implying that much computational cost can be saved by considering the best scored docking poses only. Besides, compared with the rescoring based on the minimized structures, the rescoring based on the MD simulations might not be completely necessary due to its negligible impact on the docking performance. Considering the much higher computational demand of MM/PBSA, MM/GBSA with a high solute dielectric constant (ε = 2 or 4) is recommended for the virtual screening of tyrosine kinases.

Entropy effects play an important role in drug–target interactions, but the entropic contribution to ligand-binding affinity is often neglected by end-point binding free energy calculation methods, such as MM/GBSA and MM/PBSA, due to the expensive computational cost of normal mode analysis (NMA). Here, we systematically investigated entropy effects on the prediction power of MM/GBSA and MM/PBSA using >1500 protein–ligand systems and six representative AMBER force fields. Two computationally efficient methods, including NMA based on truncated structures and the interaction entropy approach, were used to estimate the entropic contributions to ligand–target binding free energies. In terms of the overall accuracy, we found that, for the minimized structures, in most cases the inclusion of the conformational entropies predicted by truncated NMA (enthalpy_{nmode_min_9Å}) compromises the overall accuracy of MM/GBSA and MM/PBSA compared with the enthalpies calculated based on the minimized structures (enthalpy_min). However, for the MD trajectories, the binding free energies can be improved by the inclusion of the conformation entropies predicted by either truncated-NMA for a relatively high dielectric constant (ε_in = 4) or the interaction entropy method for ε_in = 1–4. In terms of reproducing the absolute binding free energies, the binding free energies estimated by including the truncated-NMA entropies based on the MD trajectories (ΔG_{nmode_md_9Å}) give the lowest average absolute deviations against the experimental data among all the tested strategies for both MM/GBSA and MM/PBSA. Although the inclusion of the truncated NMA based on the MD trajectories (ΔG_{nmode_md_9Å}) for a relatively high dielectric constant gave the overall best result and the lowest average absolute deviations against the experimental data (for the ff03 force field), it needs too much computational time. Alternatively, considering that the interaction entropy method does not incur any additional computational cost and can give comparable (at high dielectric constant, ε_in = 4) or even better (at low dielectric constant, ε_in = 1–2) results than the truncated-NMA entropy (ΔG_{nmode_md_9Å}), the interaction entropy approach is recommended to estimate the entropic component for MM/GBSA and MM/PBSA based on MD trajectories, especially for a diverse dataset. Furthermore, we compared the predictions of MM/GBSA with six different AMBER force fields. The results show that the ff03 force field (ff03 for proteins and gaff with AM1-BCC charges for ligands) performs the best, but the predictions given by the tested force fields are comparable, implying that the MM/GBSA predictions are not very sensitive to force fields.