Novel structural features of CDK inhibition revealed by an ab initio computational method combined with dynamic simulations

The rational development of specific inhibitors for the ~500 protein kinases encoded in the human genome is impeded by a poor understanding of the structural basis for the activity and selectivity of small molecules that compete for ATP binding. Combining classical dynamic simulations with a novel ab initio computational approach linear-scalable to molecular interactions involving thousands of atoms, we have investigated the binding of five distinct inhibitors to the cyclin-dependent kinase CDK2. We report here that polarization and dynamic hydrogen bonding effects, so far undetected by crystallography, affect both their activity and selectivity. The effects arise from the specific solvation patterns of water molecules in the ATP binding pocket or the intermittent formation of hydrogen bonds during the dynamics of CDK/inhibitor interactions and explain the unexpectedly high potency of certain inhibitors such as 3-(3H-imidazol-4-ylmethylene)-5-methoxy-1,3-dihydro-indol-2-one (SU9516). The Lys89 residue in the ATP-binding pocket of CDK2 is observed to form temporary hydrogen bonds with the three most potent inhibitors. This residue is replaced in CDK4 by Thr89, whose shorter side-chain cannot form similar bonds, explaining the relative selectivity of the inhibitors for CDK2. Our results provide a generally applicable computational method for the analysis of biomolecular structures and reveal hitherto unrecognized features of the interaction between protein kinases and their inhibitors

💡 Research Summary

The paper addresses the longstanding challenge of rationally designing selective inhibitors for the cyclin‑dependent kinase CDK2, a key target in cancer therapy, by integrating large‑scale ab initio density‑functional theory (DFT) with classical molecular dynamics (MD). Five chemically diverse inhibitors—NU2058, NU6027, NU6102, the 9d‑variant of NU6027, and SU9516—were examined. The authors first defined a binding‑pocket model that includes all amino‑acid residues within 7 Å of the ligand, a size shown to converge binding energies within 0.15 kcal mol⁻¹ using a convergence test with the large inhibitor staurosporine. Linear‑scaling O(N) DFT codes (ONETEP and Siesta) were then employed to calculate gas‑phase binding enthalpies (ΔE_g) for each complex. Solvation contributions (ΔG_solv) were added using a Generalized Born/Surface Area (GBSA) continuum model, while entropic differences between inhibitors were assumed negligible for relative comparisons. This yields a simple expression for the relative free‑energy change: ΔΔG ≈ ΔΔE_g + ΔΔG_solv.

The key findings go beyond static crystal structures. MD simulations revealed that water molecules occupying the ATP‑binding pocket form transient hydrogen‑bond networks with the ligands, significantly modulating the electronic environment and enhancing binding affinity. These water‑mediated interactions, invisible in crystallography, were quantified by the DFT calculations and contributed up to several kcal mol⁻¹ to ΔΔG. Moreover, the side‑chain of Lys89 in CDK2 was observed to establish intermittent hydrogen bonds with the three most potent inhibitors (SU9516, NU6102, 9d‑NU6027) for 30–45 % of the simulation time. This dynamic interaction stabilizes the complex and explains the high potency of these compounds. In CDK4, the analogous position is occupied by Thr89, whose shorter side‑chain cannot form comparable bonds, providing a structural rationale for the observed CDK2‑selectivity.

When only static binding energies were considered, the predicted potency ranking disagreed with experimental IC₅₀/K_i values. Inclusion of the dynamic water network and Lys89 hydrogen‑bonding restored agreement, reproducing the experimental relative K_i values within ~1 kcal mol⁻¹. The authors also confirmed that dispersion forces, although important for absolute binding energies, contribute similarly across all inhibitors and thus do not affect relative potency.

Overall, the study demonstrates that (i) linear‑scaling ab initio DFT can be feasibly applied to protein‑ligand systems involving thousands of atoms, (ii) explicit treatment of solvent dynamics is essential for accurate prediction of inhibitor potency and selectivity, and (iii) the Lys89‑mediated hydrogen‑bond network is a decisive factor for CDK2 selectivity over CDK4. The combined computational framework offers a generalizable tool for the design of kinase inhibitors and highlights the importance of dynamic, water‑mediated interactions that are missed by static structural analyses alone.