Numerical errors in minimization based binding energy calculations

This work examines the effect of small input perturbations on binding energies computed from differences between energy minimized structures, such as the Prime MM-GBSA and MOE MM-GB/VI methods. The applied perturbations include translations of the cognate ligand in the binding site by a maximum of 0.1 Å along each coordinate or the permutation of the order of atoms of the cognate ligand without any changes to the atom coordinates. These seemingly inconsequential input changes can lead to as much as 17 kcal/mol differences in the computed binding energy. The calculated binding energies cluster around discrete values, which correspond to specific ligand poses. It appears that the largest variations are observed for target-ligand systems in which there is a possibility for multiple poses with strong hydrogen bonds. The barriers between different poses can appear fractal-like, making it difficult to predict which solution will be produced from a given input. Including protein flexibility in MM-GBSA calculations further increases numerical instability, and the protein strain terms seem to be the major factor contributing to this sensitivity. In such calculations it appears unwise to extend the flexible region beyond 6 Å.