Collinearity of Decomposed Energy Terms in MM-P(G)BSA Binding Free Energy Calculations of Protein-Protein Complexes

Abstract

The molecular mechanics (MM) Poisson-Boltzmann (PB) or generalized Born (GB) surface area methods (MM-P(G)BSA) are among the most commonly used end state approaches for the calculation of the binding free energies (BFEs) in computational drug design and screening studies. Their thermodynamic cycle is based on the decomposition of the free energy change into several energy terms, including molecular mechanics (MM) electrostatic and van der Waals terms representing the gas phase binding energy change, as well as implicit solvation energies comprising polar solvation, calculated by solving PB equation or its reduced form, GB, and nonpolar solvation surface area (SA) terms. Although these terms are additive, and thus the sum of these components should yield the total free energy change, most of these terms are represented with coefficients that are physically meaningful but empirical. To date, a great deal of effort has been made to improve the accuracy of these methods in the prediction of experimental binding free energies, at least relatively, i.e., relative binding free energy (RBFE), resulting in totally empirical coefficients on energy terms and breaking the underlying physics. Furthermore, although these methods originated from protein-ligand RBFE estimation, there have been examples or even tutorials using these methods for protein-peptide or protein-protein interactions (PPIs), due to the same simple physics of binding at the bound state, without sufficient hesitation regarding their accuracy limitations. Here, we thoroughly evaluate MMP(G)BSA methods in all variants available in the literature for a diverse protein-protein complex set and reveal their true RBFE accuracy by means of Pearson correlation with experimental values, with at most R = 0.42. We also demonstrate that the assumption of independent fitting coefficients for decomposed energy terms could not only break the physics but also statistically represent overfitting. Through analytic derivation and large-scale molecular dynamics simulations, we show that (i) the protein-ligand (PL), protein-peptide or protein-protein (PP) Coulomb interaction energy and the GB solvation correction are almost perfectly collinear (R2[≥]0.99) reflecting their designed role as vacuum electrostatics plus solvent screening, and (ii) the van der Waals interaction and SA term also exhibit strong correlation. Interaction entropy (IE) and C2 entropy corrections, which are also found to be strongly dependent on each other, worsen the overall accuracy of the methods due to the large energetic fluctuations of these extended protein-protein systems. The findings of collinearity hold both at the level of instantaneous trajectory fluctuations and when averaged across a diverse set of 138 PP complexes and persist in both single-trajectory and three-trajectory MM-GBSA protocols. Our results show the limitation of the methods for protein-protein complexes and caution against using decomposed MM-GBSA terms as independent predictors in regression models and suggest instead combining correlated terms into effective polar, nonpolar, and entropic contributions, while noting the adverse impact of these entropy terms.