Benchmarking Free Energy Computational Methods for Revealing the Interactions Driving PARP1 Selective Inhibition

Accurate prediction of inhibitor selectivity across protein paralogues remains a central challenge in computational drug discovery. Here, we systematically benchmark three computational methods--Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA), free energy perturbation (FEP) and potential of mean force (PMF) calculations--in their ability to recapitulate PARP1 versus PARP2 selectivity for eight clinically relevant PARP enzyme inhibitors used in ovarian, breast and prostate tumors among others. We demonstrate how MM/PBSA calculations offer rapid and qualitative insights, but show pronounced sensitivity to the chosen static conformational pose, being particularly challenging for ligands with subtle energetic differences between distinct protein paralogues. In contrast, both FEP and PMF calculations using atomistic models with explicit solvent result in substantially improved agreement with experimental binding affinities. The FEP method exhibits the strongest quantitative correlation with experimental binding free energy differences, remarkably reproducing selectivity trends even among nearly isoenergetic complexes. Notably, our structural contact analysis reveals how contact connectivity controls ligand selectivity, providing valuable mechanistic and molecular insight into the key residues that stabilize each inhibitor in both protein enzymes. Together, our multi-method computational study contributes to elucidate potential chemical modifications across the ligand chemical space to enhance potency and specificity, informing the future design and evaluation of selective inhibitors for precision oncology, including therapies targeting homologous recombination-deficient cancers.