A Hierarchical Free Energy Framework for Ranking CDK Inhibitors: From PMF Barriers to Alchemical Precision

A Hierarchical Free Energy Framework for Ranking CDK Inhibitors: From PMF Barriers to Alchemical Precision

 

Ganga Graziela Christina¹, Chen qu²

 

chrisnalove23@gmail.com, Zhejiang University of Science and Technology, Hangzhou 310000, Zhejiang, China.

chrisnalove23@gmail.com, Zhejiang University of Science and Technology, Hangzhou 310000, Zhejiang, China.

 

Abstract

   Accurate prediction of binding free energies for structurally homologous kinase targets remains one of the most challenging tasks in computational drug discovery. CDK1 and CDK2 provide an exceptionally stringent benchmark: their active sites differ by only 0.72 Å backbone RMSD, yet experimental ITC data reveal up to 170-fold selectivity differences for certain inhibitors, corresponding to a ΔΔG of only 2.87 kcal/mol, which is smaller than the typical error of most computational methods. Here we establish and validate a three-tier hierarchical free energy protocol applied to five CDK–inhibitor complexes (Dinaciclib, AZD5438, and CGP74514A bound to CDK1 and CDK2). Tier 1 (MM/PBSA endpoint calculations) achieves a qualitative affinity ranking but fails quantitatively, yielding R ≈ 0.45 and a critical rank inversion for CDK2-Dinaciclib. Tier 2 (umbrella sampling with 56 harmonically restrained windows and WHAM reconstruction) correctly identifies the primary free energy barrier to dissociation in the 6-12 Å center-of-mass distance regime and correctly rank-orders all complexes, though absolute values are systematically underestimated. Tier 3 (alchemical double-decoupling via 62 λ-states with thermodynamic integration and MBAR cross-validation converged to within 0.2 kcal/mol) achieves near-chemical accuracy with a mean absolute error of ~1.11 kcal/mol and R ≈ 1.0.

Keywords: free energy calculation; alchemical methods; umbrella sampling; MM/PBSA; CDK inhibitors; thermodynamic integration; MBAR; kinase selectivity

1. Introduction

   The prediction of protein–ligand binding free energies with sufficient accuracy to guide medicinal chemistry decisions is a long-standing goal of computational drug discovery. Modern alchemical FEP protocols can achieve chemical accuracy (~1 kcal/mol) for congeneric ligand series, and prospective FEP campaigns have demonstrated genuine utility in reducing synthesis cycles [1] . Yet highly homologous kinase pairs where the goal is to distinguish selectivity between nearly identical binding sites remain a particularly demanding frontier [2]  .

   A fundamental relationship quantifies the problem's sensitivity. At 300 K, an error of 1.4 kcal/mol in the calculated ΔGbind corresponds to an order-of-magnitude error in the predicted Kd. For CDK1 and CDK2, the experimental selectivity window for the best inhibitors spans only a ΔΔG of 2.87 kcal/mol, less than twice this error threshold. A computational method must achieve sub-kcal/mol accuracy to discriminate CDK2-selective from non-selective compounds reliably [3]  .

   CDK1 and CDK2 share a global backbone RMSD of ~0.72 Å and near-superimposable ATP-binding pockets. Yet, ITC measurements reveal striking isoform selectivity for inhibitors, including Dinaciclib (~23-fold CDK2 preference) and AZD5438 (~170-fold CDK2 preference; Kd = 26 ± 3 nM for CDK2 versus 4,400 ± 3,200 nM for CDK1). Previous studies applying MM/PBSA to these systems correctly reproduced qualitative trends but demonstrated rank inversions when entropic contributions are poorly estimated. Higher-level methods have not previously been systematically applied and benchmarked across this full panel.

   Pathway-based methods, such as steered by molecular dynamics and umbrella sampling, offer a physically interpretable complement to alchemical approaches: they construct the potential of mean force (PMF) along the actual ligand dissociation coordinate, revealing not only the thermodynamic binding free energy but also the kinetic barrier structure[4] . This mechanistic information is inaccessible to alchemical methods but highly valuable for understanding drug residence time  [5] .

   In this work, we report the systematic application and benchmarking of three tiers of free energy calculation  MM/PBSA endpoint methods, umbrella-sampling-based PMF reconstruction, and alchemical TI/MBAR double-decoupling across five CDK–inhibitor complexes [6]. We benchmark each tier against experimental ITC reference values (Table 1), identify the failure modes of lower-level methods, and demonstrate the conditions under which near-chemical accuracy is achievable.

Table 1. Experimental ITC reference data for CDK–inhibitor complexes.

No.	PDB	Complex	Kd (nM)	ΔGexp (kcal/mol)
1	6GU6	CDK1–Dinaciclib	955 ± 246	−8.15 ± 0.12
2	4KD1	CDK2–Dinaciclib	41 ± 14	−9.35 ± 0.16
3	6GU7	CDK1–AZD5438	4,400 ± 3,200	−6.79 ± 0.42
4	6GUH	CDK2–AZD5438	26 ± 3	−9.54 ± 0.07
5	6GUK	CDK2–CGP74514A	715 ± 15	−7.24 ± 0.01

 

2. Computational Methods

2.1 MM/PBSA End-Point Calculations (Tier 1)

   MM/PBSA binding free energies were computed from snapshots extracted from the final 10 ns (90–100 ns) of equilibrated 100 ns production MD trajectories, using the g_mmpbsa tool interfaced with GROMACS 4.6.5. The polar solvation contribution was computed by solving the linearised Poisson-Boltzmann equation using the APBS package; the non-polar contribution was estimated as a linear function of the solvent-accessible surface area [7] . The entropic correction was applied via the second-order cumulant expansion (C2) approximation, which estimates binding entropy from the variance of receptor–ligand interaction energies across trajectory snapshots  [8] [9] .

2.2 Steered Molecular Dynamics and Umbrella Sampling (Tier 2)

   Dissociation pathways for all five complexes were predicted using MoMA-LigPath, generating continuous conformational snapshots spanning a center-of-mass (COM) distance ξ from 3 Å to 25 Å  [10] . Steered molecular dynamics (SMD) simulations employed a harmonic spring (k = 1,000 kJ/mol/nm²) pulling the ligand COM along ξ at a constant velocity of 0.5 Å/ns, selected following systematic testing at four speeds. For each complex, 50 independent SMD trajectories were generated with Gaussianity of work distributions verified by the Shapiro-Wilk test (α = 0.05). Non-equilibrium work was analyzed using the Jarzynski equality and a second-order cumulant estimator [10]  [11].

   Umbrella sampling (US) windows were distributed along ξ from 3 Å to 25 Å: 44 base windows at 0.5 Å spacing, refined to 0.25 Å within the critical ξ = 6–12 Å regime, with 12 additional windows, yielding 56 windows per complex. A harmonic biasing potential of 3,000 kJ/mol/nm² was applied. Each window underwent 500 ps of equilibration followed by 10 ns of production sampling using PLUMED 2.8, interfaced with GROMACS. The unbiased PMF was reconstructed using WHAM, with 200 bootstrap samples to estimate statistical uncertainties.

2.3 Alchemical Double-Decoupling TI/MBAR (Tier 3)

   Absolute binding free energies were computed via double-decoupling thermodynamic cycles (Figure 1). Each leg comprised 31 λ-windows: 11 for electrostatic annihilation and 20 for van der Waals decoupling, for a total of 62 discrete alchemical states per system  [12]  . Within each window, 500 ps equilibration was followed by 5 ns production sampling. Free energies were extracted by Gaussian quadrature over the TI integrand. MBAR served as the secondary estimator; TI and MBAR results were required to agree within 0.2 kcal/mol before acceptance. Phase-space overlap between adjacent λ pairs was monitored; any pair below a 0.03 overlap threshold triggered insertion of additional windows  [13] .

Figure 1. Thermodynamic cycle for absolute binding free energy calculation via double-decoupling. ΔGbind = ΔG3 − ΔG1, where ΔG3 is alchemical decoupling in the bound environment, and ΔG1 is decoupling in bulk solvent.

3. Results

3.1 Tier 1: MM/PBSA  Performance and Failure Modes

   MM/PBSA calculations capture the correct qualitative trend but fail significantly in quantitative accuracy. The most critical failure is a rank inversion for CDK2–Dinaciclib: the method calculates ΔGbind of −16.87 kcal/mol, whereas −28.52 kcal/mol is calculated for CDK1–Dinaciclib, directly contradicting the experimental finding that CDK2 binds ~23-fold more tightly. The origin is a massive entropic penalty (−TΔS = 28.04 kcal/mol) that the C2 approximation captures but cannot correctly balance against enthalpic gains from the induced-fit reorganization. Across all five complexes, the Pearson correlation between calculated and experimental ΔG values is R ≈ = 0.45, which is insufficient for reliable compound ranking [14] [15]  .

3.2 Dissociation Pathway Geometry

   MoMA-LigPath successfully generated physically plausible dissociation pathways for all five complexes, producing continuous conformational snapshots spanning ξ = 3 Å to ξ = 25 Å (Figure 2). The predicted pathways route the ligand through the mouth of the ATP-binding cleft between the N- and C-terminal lobes. CDK2 complex pathways traverse a narrower exit channel, consistent with a tighter hydrophobic enclosure. In contrast, CDK1 pathways show a more open channel geometry, reflecting the greater rigidity and looser lobe packing of this isoform  [16] [17] .

Figure 2. MoMA-LigPath predicted dissociation pathways for CDK–inhibitor complexes. Each panel shows the trajectory from the bound state (left) to the fully solvated state (right). CDK1 in cyan ribbons; CDK2 in grey. Dinaciclib in red spheres; AZD5438/CGP74514A in green.

 

3.3 SMD Force Profiles and Rupture Mechanics

   SMD force-time profiles from 50 independent trajectories per complex reveal the mechanical signature of the primary rupture event (Figure 3). For CDK1–Dinaciclib, the pulling force reaches a peak of approximately 800 kJ/(mol·nm) at ~300 ps, corresponding to the rupture of dominant protein–ligand interactions. Following this peak, the force drops sharply, confirming the ligand has achieved a fully solvated state. Work distributions passed the Shapiro-Wilk normality test (α = 0.05) for all five complexes, validating the use of the second-order cumulant estimator  [16] .

Figure 3. Representative SMD force-time profile for CDK1–Dinaciclib. The peak at ~300 ps marks the rupture of primary intermolecular interactions; the subsequent force decay confirms the transition to the fully solvated state.

3.4 Umbrella Sampling: Validation and PMF Profiles

   Sampling histograms confirm continuous, substantial overlap across the entire reaction coordinate (Figure 4), validating the choice of biasing potential force constant and window spacing. The resulting PMF profile for CDK1–Dinaciclib (Figure 5) shows the steepest gradient between 0.6 and 1.2 nm, in the 6–12 Å regime, where simultaneous rupture of conserved hinge-region hydrogen bonds and reorganization of the interfacial hydration shell occur. Beyond ξ = 2.0 nm, the free energy reaches a stable plateau signifying complete solvation. US-derived binding free energies range from −3.98 to −5.13 kcal/mol (Table 2), systematically underestimating experimental values by ~2–4 kcal/mol without standard-state volume corrections, but the relative ordering of all five complexes is correct [18]  .

Figure 4. Representative umbrella sampling histograms along the reaction coordinate ξ, confirming adequate phase-space overlap between adjacent windows, a prerequisite for accurate PMF reconstruction by WHAM.

Figure 5. PMF profile for CDK1–Dinaciclib dissociation. Free energy (kcal/mol) as a function of COM distance ξ (nm). Error bars represent 200-bootstrap statistical uncertainty.

Table 2. Umbrella sampling ΔGbind values versus experimental ITC data.

No.	Complex	ΔGbind US (kcal/mol)	ΔGexp (kcal/mol)
1	CDK1–Dinaciclib	−4.58 ± 0.49	−8.15 ± 0.12
2	CDK2–Dinaciclib	−5.13 ± 0.51	−9.35 ± 0.16
3	CDK1–AZD5438	−4.10 ± 0.42	−6.79 ± 0.42
4	CDK2–AZD5438	−4.98 ± 0.53	−9.54 ± 0.07
5	CDK2–CGP74514A	−3.98 ± 0.38	−7.24 ± 0.01

 

3.5 Alchemical λ-Decoupling Profiles

   The alchemical decoupling of Dinaciclib illustrates the thermodynamic character of the two-stage transformation. During the electrostatic annihilation stage (Figure 6), the cumulative ΔG decreases monotonically to approximately −72 kT at full decoupling. The van der Waals decoupling stage (Figure 7) exhibits a characteristic non-monotonic ΔG profile, with an initial peak of approximately 8.7 kT per window as repulsive Lennard-Jones cores are softened, followed by a decrease in later windows  [19] . TI and MBAR estimates agreed within the 0.2 kcal/mol acceptance criterion across all five systems, confirming convergence.

Figure 6. Electrostatic free energy changes for Dinaciclib decoupling in bulk solvent. Cumulative ΔG (kT) for transitions between five discrete λele states at fixed λvdW = 0.00.

Figure 7. van der Waals free energy changes during Dinaciclib decoupling. Individual ΔG per window transition for λvdW States 6–20 at fixed λele = 1.00.

3.6 Method Benchmarking

   A direct comparison of predictive performance across all three tiers is summarised in Figure 8 and Table 3. The FEP/alchemical tier achieves R ≈ 1.0 with experimental data, far exceeding the performance of umbrella sampling and MM/PBSA. The FEP mean absolute error of ~1.11 kcal/mol closely approaches the 1.0 kcal/mol chemical accuracy threshold. A consistent systematic overestimation of approximately 1.1 kcal/mol is observed as a well-understood consequence of fixed-charge, non-polarisable force fields. Despite this absolute bias, relative ΔΔG values remain highly accurate, enabling reliable interpretation of structure-activity relationships   .

Figure 8. Comparison of predictive accuracy across three computational tiers. Pearson R for MM/PBSA (black), Umbrella Sampling (red), and FEP/alchemical methods (green) against experimental ITC ΔG values.

Table 3. FEP/alchemical ΔGbind values versus experimental ITC data.

No.	Complex	ΔGbind FEP (kcal/mol)	ΔGexp (kcal/mol)
1	CDK1–Dinaciclib	−9.21 ± 0.92	−8.15 ± 0.12
2	CDK2–Dinaciclib	−10.42 ± 0.69	−9.35 ± 0.16
3	CDK1–AZD5438	−8.16 ± 0.75	−6.79 ± 0.42
4	CDK2–AZD5438	−10.37 ± 0.86	−9.54 ± 0.07
5	CDK2–CGP74514A	−8.48 ± 0.87	−7.24 ± 0.01

 

4. Discussion

   The three-tier hierarchy reflects a fundamental trade-off between computational cost and predictive accuracy. MM/PBSA's failure is instructive: the rank inversion for CDK2–Dinaciclib arises from a systematic inability to balance large entropic penalties with large enthalpic gains. This is a fundamental limitation of endpoint methods for systems in which binding involves substantial conformational adaptation — precisely the induced-fit scenario characterizing CDK2  [19] .

   The umbrella sampling tier resolves this by constructing the full PMF along the physical dissociation coordinate. The identification of the 6–12 Å regime as the primary barrier zone is mechanistically meaningful: this distance range corresponds precisely to the point at which the ligand breaks its final direct hydrogen bonds with the hinge region and has not yet shed the structured water molecules mediating indirect contacts. The systematic underestimation of absolute ΔGbind by the US tier (~2–4 kcal/mol) arises from the absence of standard-state corrections, which contribute approximately 2–3 kcal/mol at standard concentration.

   The FEP tier's near-chemical accuracy demonstrates that the alchemical framework, properly implemented with 62 λ-states and rigorous convergence criteria, can resolve subtle selectivity differences even when ΔΔG is only 2.87 kcal/mol. The residual ~1.1 kcal/mol systematic overestimation points to force-field polarisability as the next methodological frontier. Practical recommendations: MM/PBSA for rapid virtual screening; US for mechanistic insights into dissociation pathways; FEP for final compound ranking and SAR interpretation [20]  .

5. Conclusion

   We have established and validated a three-tier hierarchical free energy framework for CDK inhibitor ranking, benchmarked against experimental ITC data for five CDK1/CDK2–inhibitor complexes. MM/PBSA (Tier 1) achieves a qualitative ranking but fails quantitatively, with R ≈ 0.45 and a critical rank inversion driven by uncontrolled entropy estimation. Umbrella sampling (Tier 2) correctly ranks all complexes. It localizes the primary dissociation barrier to the 6–12 Å COM distance regime, but systematically underestimates absolute binding free energies without a standard-state correction. Alchemical TI/MBAR (Tier 3) achieves near-chemical accuracy (MAE ~1.11 kcal/mol, R ≈ 1.0)  and correctly reproduces CDK2 > CDK1 selectivity even for the most challenging comparison (ΔΔG = 2.21 kcal/mol). The established protocol provides a validated computational template for the design of next-generation CDK inhibitors.

References.

 [1]         A. V. Sadybekov and V. Katritch, “Computational approaches streamlining drug discovery,” Nature, vol. 616, no. 7958, pp. 673–685, Apr. 2023, doi: 10.1038/s41586-023-05905-z.

[2]          R. Abel, L. Wang, E. D. Harder, B. J. Berne, and R. A. Friesner, “Advancing Drug Discovery through Enhanced Free Energy Calculations,” Acc. Chem. Res., vol. 50, no. 7, pp. 1625–1632, Jul. 2017, doi: 10.1021/acs.accounts.7b00083.

[3]          D. J. Wood et al., “Differences in the Conformational Energy Landscape of CDK1 and CDK2 Suggest a Mechanism for Achieving Selective CDK Inhibition,” Cell Chemical Biology, vol. 26, no. 1, pp. 121-130.e5, Jan. 2019, doi: 10.1016/j.chembiol.2018.10.015.

[4]          C. Jarzynski, “Nonequilibrium Equality for Free Energy Differences,” Phys. Rev. Lett., vol. 78, no. 14, pp. 2690–2693, Apr. 1997, doi: 10.1103/PhysRevLett.78.2690.

[5]          G. M. Torrie and J. P. Valleau, “Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling,” Journal of Computational Physics, vol. 23, no. 2, pp. 187–199, Feb. 1977, doi: 10.1016/0021-9991(77)90121-8.

[6]          W. Liu, X. Xu, and C. Mu, “Experimental Study on Two-Phase Explosion Suppression of Gas/Pulverized Coal by Explosion Suppressant,” ACS Omega, vol. 7, no. 19, pp. 16644–16652, May 2022, doi: 10.1021/acsomega.2c00987.

[7]          D. Van Der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark, and H. J. C. Berendsen, “GROMACS: Fast, flexible, and free,” J Comput Chem, vol. 26, no. 16, pp. 1701–1718, Dec. 2005, doi: 10.1002/jcc.20291.

[8]          R. Kumari, R. Kumar, Open Source Drug Discovery Consortium, and A. Lynn, “g_mmpbsa —A GROMACS Tool for High-Throughput MM-PBSA Calculations,” J. Chem. Inf. Model., vol. 54, no. 7, pp. 1951–1962, Jul. 2014, doi: 10.1021/ci500020m.

[9]          M. H. Bikoumou Gambat, H. Li, and J. Zhang, “Breaking barriers to integrated digital Delivery (IDD) adoption in the construction Industry: Evidence from five Central African countries,” Ain Shams Engineering Journal, vol. 17, no. 7, p. 104216, Jul. 2026, doi: 10.1016/j.asej.2026.104216.

[10]        J. I. Baños-Sanz et al., “Crystal structure and functional insights into uracil-DNA glycosylase inhibition by phage ϕ29 DNA mimic protein p56,” Nucleic Acids Research, vol. 41, no. 13, pp. 6761–6773, Jul. 2013, doi: 10.1093/nar/gkt395.

[11]        S. A. Maxime, B. G. M. Horacio, and M. S. Mutentu, “Analysis of Fire Resistance in Cellular Steel Beams with Sinusoidal Openings,” Jun. 2024, doi: 10.5281/ZENODO.11473823.

[12]        S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen, and P. A. Kollman, “THE weighted histogram analysis method for free‐energy calculations on biomolecules. I. The method,” J Comput Chem, vol. 13, no. 8, pp. 1011–1021, Oct. 1992, doi: 10.1002/jcc.540130812.

[13]        M. J. Abraham et al., “GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers,” SoftwareX, vol. 1–2, pp. 19–25, Sep. 2015, doi: 10.1016/j.softx.2015.06.001.

[14]        P. V. Klimovich, M. R. Shirts, and D. L. Mobley, “Guidelines for the analysis of free energy calculations,” J Comput Aided Mol Des, vol. 29, no. 5, pp. 397–411, May 2015, doi: 10.1007/s10822-015-9840-9.

[15]        M. S. Mutentu, B. G. M. Horacio, and Y. Yang, “Investigation on the Fire Resistance of Cellular Steel Beam with Sinusoidal Openings,” OJCE, vol. 13, no. 04, pp. 637–663, 2023, doi: 10.4236/ojce.2023.134043.

[16]        D. J. Wood et al., “Differences in the Conformational Energy Landscape of CDK1 and CDK2 Suggest a Mechanism for Achieving Selective CDK Inhibition,” Cell Chemical Biology, vol. 26, no. 1, pp. 121-130.e5, Jan. 2019, doi: 10.1016/j.chembiol.2018.10.015.

[17]        Y. Wang, M. L. Diabakanga Batatana, and M. H. Bikoumou Gambat, “Public perceptions of government policies to COVID-19: A comparative study in six African countries,” Heliyon, vol. 10, no. 3, p. e24888, Feb. 2024, doi: 10.1016/j.heliyon.2024.e24888.

[18]        A. V. Sadybekov and V. Katritch, “Computational approaches streamlining drug discovery,” Nature, vol. 616, no. 7958, pp. 673–685, Apr. 2023, doi: 10.1038/s41586-023-05905-z.

[19]        W. Liu, X. Xu, and C. Mu, “Experimental Study on Two-Phase Explosion Suppression of Gas/Pulverized Coal by Explosion Suppressant,” ACS Omega, vol. 7, no. 19, pp. 16644–16652, May 2022, doi: 10.1021/acsomega.2c00987.

[20]        M. R. Shirts and J. D. Chodera, “Statistically optimal analysis of samples from multiple equilibrium states,” The Journal of Chemical Physics, vol. 129, no. 12, p. 124105, Sep. 2008, doi: 10.1063/1.2978177.