The study began with the structural preparation of CDK2, CDK5, CDK4, CDK6, and CDK12, based on the monomeric catalytic subunits. The starting coordinates were derived from crystal structures of the kinases in their biologically active complexes (PDB IDs: 2WIH, 1UNL, 7SJ3, 5L2S, and 7NXK, respectively), this ensured that the ATP-binding pockets were initiated in their active conformations, with the activation loops and C-helices properly oriented. The crystal structures of CDK4, CDK6, and CDK12 were of relatively low resolution (R-factor >3 Å), resulting in poorly defined regions within the protein scaffold, particularly in flexible loop domains. Although this limitation did not affect the molecular docking phase, since the active site resided in a well-resolved region of the protein, accurate reconstruction of the missing segments was essential for reliable molecular dynamics simulations. To address this, homology modeling was performed using RoseTTAFold via the Robetta server (Baek et al., 2021), based on the UniProt sequences of CDK4, CDK6, CDK12). This approach preserved the experimentally resolved regions from the X-ray structures while predicting plausible conformations for the missing loops, resulting in complete and structurally coherent starting models.

To ensure conformational stability prior to docking, all protein structures were subjected to energy minimization followed by short relaxation molecular dynamics simulations using GROMACS 2023.3 (https://www.gromacs.org). Molecular docking was then performed with AutoDock Vina 1.1.2 (Forli et al., 2016). A computational docking protocol was employed to evaluate the interaction between each protein (CDK2, CDK4, CDK5, CDK6, CDK12) and a diverse panel of six ligands: Roscovitine, CR8, DRF053, N&N1, LGR1406, and ATP. For each of the 30 protein–ligand complexes (5 proteins × 6 ligands), we performed 25 independent docking runs and retained the best-scoring pose (selected as the best among the 25 best poses obtained from the 25 docking launches). This step enabled the identification of the most favorable binding poses of 30 complexes and provided estimates of binding affinities.

To assess the stability of the docked complexes, the top 30 ranked poses for each ligand-protein system were selected for 1 μs (1,000 ns) molecular dynamics (MD) simulations. For comparison, five simulations of the apoprotein (ligand-free) were also performed under identical conditions. All simulations were performed using GROMACS 2024.2 with the AMBER99SB-ILDN force field (Loschwitz et al., 2021). Topologies for each compound were generated using Acpype v2023.10.27 (Sousa et al., 2012), which assigns AMBER atom types. This ensured the ATP-binding pocket was initiated in its biologically active conformation. The starting structures for the MD simulations (see Table 1 - Datadir 4) were sourced as follows: the 30 highest-scoring ligand-receptor complexes were taken directly from the docking results, while the five apoprotein starting structures were obtained from the final step of the Structure Relaxation protocol. Each complex was placed at the center of the simulation box, solvated with TIP3P water molecules (Mark and Lennart, 2001), and neutralized by adding appropriate counterions. Energy minimization was performed using the steepest descent algorithm, allowing up to 50,000 steps and terminating when the maximum force dropped below 1,000 kJ·mol⁻¹·nm⁻¹. The system was first equilibrated for 1 ns in the NVT ensemble at 300 K using the V-rescale thermostat, followed by 2 ns of NPT equilibration using the C-rescale barostat set to 1 bar. Long-range electrostatics were computed using the GPU-accelerated particle-mesh Ewald (PME) method for each protein-ligand complex. Each equilibrated system was then subjected to a 1 μs production run, integrated with a 2 fs timestep, with coordinates saved every 10 ps. The resulting trajectories were analyzed to characterize the time-dependent behavior of protein-ligand interactions. Two primary analyses, root-mean-square deviation (RMSD) and root-mean-square fluctuation (RMSF), were performed to assess the overall stability and residue-level flexibility of each complex. RMSD and RMSF curves are reported as time series from the individual trajectories for each CDK-ligand complex. To provide a quantitative estimate of variability, RMSD values were summarized as mean ± standard deviation computed over trajectory frames in the analysis window starting at 0.2 µs (Tables 3, 4).

The molecular mechanics Poisson-Boltzmann surface area (MM/PBSA) method (Kumari et al., 2014) was used to estimate the binding free energies of the 30 complexes analyzed in this study. Calculations were performed on 8,000 frames extracted from each 1 μs molecular dynamics trajectory, starting at 0.200 μs, to capture the most representative conformations of each protein-ligand complex. Principal component analysis (PCA) was employed to reduce the dimensionality of the simulation data, thereby allowing for the identification and analysis of large-scale, collective motions. This widely used technique projects high-dimensional datasets onto a smaller set of orthogonal axes (principal components) that capture the majority of the variance. The analysis began with the construction of the covariance matrix, whose eigen-decomposition produced eigenvalues and corresponding eigenvectors, using the built-in modules in GROMACS. Principal components associated with the largest eigenvalues (those accounting for the greatest proportion of variance) were retained for further analysis.

Ligand-receptor interactions were identified using the MD-ligand-receptor pipeline (Pieroni et al., 2023), focusing on those maintained throughout the majority of the simulation time (last 800 ns). Interactions were classified by bond type and by their persistence over the course of the trajectory. This analysis enabled a comparative assessment of interaction stability across the different CDK-ligand complexes.

To characterize the unbinding free energy landscapes of the LGR1406, DRF053, and N&N1 ligands from the CDK5 receptor, Umbrella Sampling (US) simulations were performed using the GROMACS 2024.2 suite. The initial structures were derived from molecular docking complexes; to define the reaction coordinate, the primary exit tunnel was identified using CAVER 3.0.3 PyMOL plugin (Pavelka et al., 2016), and the systems were structurally aligned to the calculated tunnel vector to facilitate pulling along the z-axis. Ligand topologies were generated using Acpype v2023.10.27 with the General Amber Force Field (GAFF), while the protein was parameterized with the AMBER99SB-ILDN force field. Subsequently, the complexes were solvated in a rectangular box with a 7 nm buffer along the pulling direction and neutralized with ions. The equilibration protocol was executed in the same way described in the MD section to relax the solvent interface. Steered Molecular Dynamics (SMD) was subsequently utilized to dissociate the ligand, applying a directional harmonic potential with a pulling rate of 0.005 nm/ps and a force constant of 600 kJ mol⁻¹ nm⁻² (Ngo et al., 2022). Each window underwent independent NPT equilibration followed by production US simulations with harmonic position restraints on the ligand. Throughout the SMD trajectory, as well as during the subsequent NPT equilibration of the extracted windows and the final production US runs, position restraints of 1,000 kJ mol⁻¹ nm⁻² were applied to the protein C-α atoms to prevent structural drift. The final Potential of Mean Force (PMF) profiles were reconstructed using the Weighted Histogram Analysis Method (WHAM) via the “gmx wham” module included in the GROMACS suite, with iterative checking to ensure sufficient histogram overlap between adjacent windows.

Drug-likeness and early ADME properties were computed with the SwissADME webserver (www.swissadme.ch) (Daina et al., 2017) for the six ligands. We recorded rule compliance and violations for the Lipinski, Veber, Ghose, Egan, and Muegge filters; screened for PAINS and Brenk structural alerts; and extracted MW, TPSA, and the consensus log Po/w.

Molecular docking simulations provided insight into the binding stability of various ligands to CDKs, revealing differences in interaction strength among the enzyme-ligand complexes. Docking scores, expressed in kcal/mol, represent predicted binding energies, with more negative values indicating stronger stability. Table 2 summarizes the best ligand-receptor binding stability scores (in kcal/mol) obtained from molecular docking simulations for each CDK-ligand complex. The docking results are reported in Table 1, Datadir 1.

Among the tested ligands, LGR1406 showed the strongest binding stability across all CDKs, with values ranging from −8.5 kcal/mol (CDK12) to −9.9 kcal/mol (CDK2). CR8 and DRF053 also displayed consistently strong interactions, particularly with CDK2 (−9.4 kcal/mol) and CDK5 (−9.2 kcal/mol for CR8; −9.0 kcal/mol for DRF053). In contrast, ATP, the physiological substrate of CDKs, showed the weakest binding stability, with docking scores ranging from −7.1 kcal/mol (CDK5) to −8.7 kcal/mol (CDK6). These results emphasize the enhanced binding capability of synthetic inhibitors within the ATP-binding pocket, reinforcing their potential as effective CDK inhibitors.

Comparative analysis among CDKs identified CDK2 as the most susceptible to inhibition, consistently displaying the most favorable binding energies, particularly with LGR1406 (−9.9 kcal/mol), DRF053 (−9.0 kcal/mol), and CR8 (−9.4 kcal/mol). CDK5 exhibited a similar pattern, with high stability for the same ligands. In contrast, CDK6 and CDK12 showed comparatively weaker interactions overall, especially with ATP (−8.7 kcal/mol and −7.7 kcal/mol, respectively) and Roscovitine (−7.6 kcal/mol and −7.8 kcal/mol, respectively).

Taken together, these findings indicate that, although all tested ligands exhibit varying degrees of binding stability across the CDK family, LGR1406, CR8, and DRF053 emerge as the most promising candidates for selective inhibition, particularly against CDK2 and CDK5. These results provide a solid basis for further validation through molecular dynamics simulations and experimental enzymatic assays.

Analysis of root mean square deviation (RMSD) provided insights into the structural stability of CDK2, CDK4, CDK5, CDK6, and CDK12 complexes with different ligands over the 1,000 ns simulation period.

The backbone conformational dynamics of five CDK isoforms were evaluated through C-α RMSD analysis in both apo and ligand-bound states (Table 3). The CDK isoforms exhibited distinct structural stability profiles, with CDK2 and CDK5 demonstrating the lowest backbone flexibility (RMSD values predominantly in the range of 0.19–0.30 nm), while CDK4, CDK6, and CDK12 showed higher conformational variability. CDK2 maintained remarkable structural rigidity across all conditions, with RMSD values ranging from 0.206 ± 0.023 nm (Roscovitine complex) to 0.297 ± 0.021 nm (N&N1 complex), values closely comparable to the apoprotein (0.239 ± 0.031 nm). Similarly, CDK5 exhibited exceptional stability with RMSD values between 0.192 ± 0.026 nm (ATP complex) and 0.267 ± 0.026 nm (apoprotein), suggesting minimal ligand-induced conformational changes in these isoforms.

In contrast, CDK6 displayed substantially higher backbone flexibility, particularly in the ATP-bound state (0.785 ± 0.042 nm), which represented nearly double the RMSD of the apoprotein (0.395 ± 0.033 nm) and significantly exceeded values observed with synthetic inhibitors (0.390–0.580 nm). This pronounced ATP-induced conformational variability suggests that CDK6 undergoes substantial structural rearrangements upon natural substrate binding. The synthetic inhibitors generally stabilized CDK6 relative to ATP, with CR8 achieving RMSD values (0.390 ± 0.029 nm) nearly identical to the apoprotein, while DRF053, LGR1406, and N&N1 showed intermediate flexibility (0.577–0.580 nm). CDK4 maintained consistently high RMSD values (0.379–0.464 nm) across all ligand-bound states, comparable to or exceeding the apoprotein (0.417 ± 0.072 nm), indicating inherent backbone flexibility in this isoform regardless of ligand occupancy. CDK12 demonstrated moderate structural stability with most complexes yielding RMSD values between 0.292 and 0.385 nm, though the N&N1 complex exhibited slightly elevated flexibility (0.435 ± 0.032 nm). Notably, Roscovitine and ATP consistently induced lower backbone deviations in CDK2, CDK5, and CDK12, suggesting these ligands may preferentially stabilize certain conformational states.

Among the ligands (mean RMSD values shown in Table 4), all ligands demonstrated stable binding configurations, with RMSD values ranging from 0.118 ± 0.016 nm to 0.288 ± 0.066 nm across the 0.2–1.0 μs analysis window. The natural substrate ATP exhibited the highest structural variability, particularly in the CDK12 complex (0.288 ± 0.066 nm) and CDK6 complex (0.267 ± 0.039 nm), while showing notably lower deviations in CDK4 (0.154 ± 0.015 nm) and CDK5 (0.129 ± 0.053 nm). In contrast, the synthetic inhibitor CR8 displayed remarkably low RMSD values in CDK12 (0.118 ± 0.016 nm), CDK6 (0.149 ± 0.012 nm), and CDK2 (0.164 ± 0.030 nm), suggesting particularly stable binding modes in these isoforms, though it showed increased flexibility in the CDK4 complex (0.286 ± 0.017 nm). The inhibitors Roscovitine, DRF053, LGR1406, and N&N1 exhibited intermediate and relatively uniform RMSD values across all CDK isoforms, typically ranging between 0.13 and 0.23 nm, indicating consistent binding behavior. Notably, LGR1406 and N&N1 demonstrated the most uniform stability profiles across all five CDK isoforms, with standard deviations generally below 0.035 nm, suggesting these compounds maintain similar binding conformations regardless of the CDK isoform.

From this point onward, the complete RMSD profiles of each ligand-receptor complex are available in the Figshare archive (Table 1 - Datadir 2). In the main article, we report, for each CDK, the RMSD trajectories of the complexes formed with the most stable ligands.

CDK2 exhibited the lowest overall RMSD values among the kinases, stabilizing around 0.25–0.30 nm for the protein c-∝. The ligands showed moderate fluctuations, with ATP and Roscovitine displaying the highest deviations, suggesting increased conformational flexibility within the binding pocket. In contrast, CR8 and DRF053 remained relatively stable throughout the simulation, indicating stronger and more persistent interactions with the active site (Figure 3). Also the Apoprotein remained close to the RMSD values of the complex ligand-protein.

CDK4 exhibited an RMSD profile comparable to that of CDK2, but with slightly higher fluctuations, stabilizing around 0.30–0.35 nm for the protein backbone. Ligand stability showed more pronounced variability: ATP displayed the highest fluctuations, while CR8 and DRF053 maintained lower RMSD values, indicating more stable binding interactions (Figure 4).

Interestingly, N&N1 and Roscovitine exhibited relatively stable behavior in this system, suggesting favorable binding. CDK5 showed greater structural variability, with protein RMSD values fluctuating between 0.5 and 0.7 nm, indicative of increased conformational dynamics. Ligand stability was variable: CR8 and N&N1 maintained lower RMSD values, consistent with a well-retained binding mode (Figure 5), whereas ATP exhibited marked fluctuations, potentially reflecting transient interactions or reduced binding affinity.

CDK6 followed a similar trend, with the protein stabilizing around 0.5–0.6 nm. Ligand RMSD values varied markedly depending on the compound. As observed for other CDKs, ATP exhibited relatively high fluctuations, supporting the notion that nucleotide ligands may induce greater conformational plasticity in this kinase. In contrast, DRF053 and CR8 showed lower RMSD values (Figure 6), indicative of stronger retention within the binding site.

CDK12 exhibited a slightly lower RMSD profile compared to CDK5 and CDK6, with protein RMSD values fluctuating around 0.3–0.4 nm. Ligand stability was more variable: ATP showed the highest fluctuations, whereas CR8 and LGR1406 displayed more restrained RMSD values, suggesting stable and persistent interactions with the kinase (Figure 7).

Comparative analysis of the RMSD data indicates that CDK2 is the most structurally stable kinase, followed by CDK4, whereas CDK5 and CDK6 exhibit the highest degrees of conformational flexibility. CDK12 falls within an intermediate stability range. Among the ligands, ATP consistently shows the highest fluctuations across all systems, suggesting a more dynamic and potentially less stable binding mode. In contrast, CR8 and DRF053 display the lowest RMSD values across multiple kinases, emerging as the most stable ligands. These findings underscore the role of the kinase structural framework in modulating ligand-induced stability, with certain ligands promoting conformational rigidity and others contributing to increased flexibility.

Root-mean-square fluctuation (RMSF) analysis enabled a comparative assessment of residue-specific flexibility across cyclin-dependent kinases (CDK2, CDK4, CDK5, CDK6, and CDK12) in complex with the ligands LGR1406, CR8, DRF053, N&N1, ATP, and Roscovitine. Overall, RMSF profiles revealed similar patterns across the kinases, with conserved core regions exhibiting low flexibility and terminal regions showing higher mobility. However, notable differences in residue-specific dynamics were observed depending on the specific CDK-ligand combination. Detailed results are stored in Table 1 - Datadir 3.

For CDK2, ATP binding results in a relatively low RMSF profile, indicating overall structural stabilization. In contrast, ligands such as CR8 and DRF053 induce increased fluctuations in specific regulatory regions, suggesting localized flexibility upon binding (Figure 8).

CDK4 exhibits a distinct behavior, characterized by enhanced flexibility in active-site loop regions regardless of the bound ligand. This pattern suggests a higher intrinsic propensity for conformational modulation within these regions (Figure 9).

In the case of CDK5, Roscovitine binding induces a pronounced rigidification of the protein structure compared to the other ligands. In contrast, DRF053 maintains a more evenly distributed fluctuation profile along the sequence, suggesting a different mode of interaction and dynamic modulation (Figure 10).

CDK6 and CDK12 exhibit RMSF patterns marked by pronounced flexibility in terminal regions and key catalytic loops. In particular, CDK6 complexes with ATP and CR8 show reduced mobility, suggesting stabilizing interactions, whereas binding of LGR1406 induces localized fluctuations within the catalytic region (Figure 11).

For CDK12, the most pronounced fluctuations are observed in the N- and C-terminal regions. Residue-specific mobility increases notably in the presence of DRF053 and CR8, while ATP binding appears to exert a stabilizing effect on the overall protein structure (Figure 12).

Overall, the RMSF analysis indicates that ligand binding differentially modulates the structural flexibility of CDKs, with effects dependent on both the specific kinase isoform and the chemical nature of the ligand. ATP-induced stabilization is particularly evident in CDK2, CDK5, and CDK6, whereas experimental compounds such as LGR1406 and DRF053 tend to increase flexibility in key regulatory regions, potentially affecting inhibition mechanisms and molecular recognition.

The analysis of Gibbs free binding energy (ΔG), performed using the MM/PBSA method, offers a quantitative evaluation of the thermodynamic stability of CDK-ligand complexes, enabling the identification of compounds with the highest binding affinities. Calculations were conducted on a dataset of 8,000 frames extracted from 1 μs molecular dynamics trajectories, capturing representative conformations of each complex. All free energy results are summarized in Table 5 and stored in Table 1 - Datadir 5.

Overall, the results indicate that CR8 and Roscovitine exhibit the most favorable binding profiles, with ΔG values consistently below −25 kcal/mol across multiple CDKs. Notably, CR8 shows the highest affinity for CDK12 (−35.75 ± 4.73 kcal/mol) and CDK6 (−28.53 ± 4.20 kcal/mol), while Roscovitine demonstrates particularly strong binding to CDK2 (−28.73 ± 5.26 kcal/mol) and CDK4 (−30.50 ± 4.95 kcal/mol). DRF053 also exhibits high binding affinity, especially towards CDK2 (−27.10 ± 4.28 kcal/mol) and CDK5 (−26.83 ± 5.25 kcal/mol).

By contrast, ATP consistently shows the least favorable binding energies (i.e., the least negative ΔG values), suggesting weaker interactions relative to synthetic inhibitors. This trend is particularly evident in CDK5 (−17.61 ± 13.31 kcal/mol) and CDK2 (−18.05 ± 11.84 kcal/mol). Interestingly, in CDK4, ATP achieves a more favorable ΔG (−31.07 ± 11.88 kcal/mol), although the high standard deviation indicates significant variability in its binding behavior.

Among the kinases, CDK12 shows the highest affinity for CR8 (−35.75 ± 4.73 kcal/mol), while CDK6 displays a strong preference for N&N1 (−34.55 ± 2.96 kcal/mol), suggesting that these ligands may serve as promising leads for selective inhibition. Additionally, LGR1406 shows strong affinity for CDK5 (−29.63 ± 3.62 kcal/mol), further supporting its potential as a CDK5-specific inhibitor.

Collectively, these findings confirm the high inhibitory potential of CR8, Roscovitine, and DRF053, which emerge as the most promising ligands across multiple CDKs. In combination with structural and molecular dynamics analyses, these results provide a compelling basis for further investigation of these compounds as selective CDK inhibitors.

To evaluate the dynamic behavior of ligands within the binding pockets of various CDKs, we performed Principal Component Analysis (PCA) on the atomic fluctuations of each ligand throughout the equilibrated portions of the MD trajectories. The projections onto the first two principal components (PC1 and PC2) reveal notable differences in conformational sampling and spatial confinement depending on the CDK partner (Supplementary Figure S2; Table 1 - Datadir 6).

Roscovitine consistently exhibits a compact and localized distribution across all CDKs, typically oscillating between two or three discrete conformational states, particularly in CDK2, CDK4, and CDK6. This confined sampling suggests highly restrained motion within the binding site, supporting its potential as a strong and stable inhibitor.

CR8 displays variable behavior depending on the CDK context. In complexes with CDK2 and CDK5, it shows intermediate flexibility with broader conformational sampling. In contrast, when bound to CDK4, CDK12, and especially CDK6, CR8 adopts a more compact distribution, indicative of enhanced stability and tight binding.

ATP, in contrast, exhibits the broadest conformational spread in all CDKs, consistent with greater internal flexibility and reduced conformational restraint within the binding pocket. This dynamic behavior aligns with its role as a natural substrate rather than a tight-binding inhibitor.

N&N1 generally show intermediate behavior, with moderate conformational confinement and a relatively compact shape in CDK6, but increased dispersion in CDK12. DRF053 also displays context-dependent variability, ranging from moderate spatial restriction (e.g., CDK2) to broader conformational freedom (e.g., CDK5), suggesting ligand-specific adaptability in different binding environments.

Overall, ligand-based PCA provides qualitative insights into conformational restraint across CDK-ligand systems. Roscovitine emerges as the most consistently stable ligand, followed by CR8 and N&N1. In contrast, ATP and DRF053 exhibit more extensive conformational sampling, potentially reflecting lower binding stability. These trends, when integrated with global and protein-based PCA results, reinforce the interpretation of Roscovitine and CR8 as ligands with higher inhibitory potential.

The following analysis, derived from molecular dynamics (MD) simulations of various ligand-protein complexes using the MD-Ligand-Receptor software, highlights key similarities and differences in binding modes across several Cyclin-Dependent Kinases (CDK2, CDK4, CDK5, CDK6, and CDK12) and a panel of ligands, including the natural substrate ATP and synthetic inhibitors such as Roscovitine, CR8, DRF053, LGR1406, and N&N1. Table 6 summarizes the temporal profiles of ligand-receptor interactions observed throughout the MD simulations. Interaction data results are stored in Table 1 - Datadir 7.

As summarized in the “Interaction Analysis” table, all Roscovitine-derived inhibitors anchor their heterobicyclic core to the canonical hinge doublet Ile10-Phe80 in CDK2, with additional reinforcement from Lys33. In contrast, ATP engages a more distal Arg157/199/200 cluster, a binding mode that contributes to the enhanced selectivity of synthetic inhibitors. In CDK5, the hydrophobic triad Ile10-Leu133-Phe80 serves as a conserved anchoring site, with LGR1406 further extending the interaction network to include Asp86, Gln130, and Asn144, correlating with its particularly favourable ΔG value.

In CDK4 and CDK6, inhibitors preferentially interact with the Phe93-Leu147 or Tyr24-Lys43 pairs, while explicitly avoiding the acidic Asp158/163 residues engaged by ATP. This selective binding behaviour may underlie their reduced interaction with kinases implicated in PD-L1 stabilization. In the case of CDK12, all ligands converge on the polar residue Asp105, yet only CR8 additionally engages Glu111 and His104, offering a structural explanation for its exceptionally high binding affinity.

Collectively, these interaction patterns demonstrate that LGR1406 and CR8 simultaneously optimize hydrophobic and polar contacts with the conserved catalytic motifs of CDK2 and CDK5, while minimizing engagement with the lateral acidic pockets typical of CDK4 and CDK6. This dual mode of binding provides a clear structural rationale for their pronounced selectivity.

Beyond direct protein-ligand contacts, the analysis of solvent-mediated interactions (Figure 13) reveals distinct solvation profiles for the natural substrate compared to the synthetic inhibitors. As indicated by the interaction persistence data, ATP consistently exhibits the highest reliance on water bridges, maintaining solvent-mediated contacts with values of 71.51% in CDK2 and 79.23% in CDK5. This high solvation correlates with the charged nature of the substrate, which also engages in extensive salt bridges (29.90% in CDK2 and 61.22% in CDK5). In contrast, the Roscovitine derivatives are primarily driven by hydrophobic forces, showing consistent hydrophobic interaction values of ∼80.00% across all complexes. However, they display notable isoform-specific variations in their water networks. In CDK2, N&N1 displays a water bridge occupancy (68.55%) comparable to the natural substrate, whereas LGR1406 exhibits the lowest solvation frequency (41.70%), suggesting a tighter, more solvent-excluded fit. Interestingly, this trend shifts in CDK5, where N&N1 becomes the least solvated ligand (45.05%), while LGR1406 maintains a more robust solvent network (57.22%). These variations suggest that while hydrophobic enclosure provides the basal affinity for the inhibitors, the modulation of water-mediated contacts contributes to the specific recognition profiles observed between CDK2 and CDK5.

To rigorously quantify the thermodynamic stability of the most relevant CDK5-ligand complexes, we performed Umbrella Sampling (US) simulations to reconstruct the Potential Mean Force (PMF) profiles along the unbinding reaction coordinate, ξ. The resulting free energy landscapes, illustrated in Figure 14 (summary free energy data are stored in Table 1 - Datafile 8, the sampling windows are addressed in Supplementary Figure S3), reveal a clear hierarchy in binding affinity among the studied inhibitors. LGR1406 exhibits the most robust interaction, characterized by a deep free energy well and a total binding free energy ( Δ G b i n d of 12.27 kcal/mol. The PMF profile for LGR1406 is notably rugged between ξ = 0.3 and 1.0 nm, featuring multiple local minima and high-energy barriers; this suggests that the unbinding process involves the sequential breaking of strong, specific non-covalent interactions and significant desolvation penalties.

In contrast, the N&N1 and DRF053 ligands display significantly lower dissociation energy costs. N&N1 demonstrates intermediate stability with a calculated binding estimate of 8.24 kcal/mol, showing a steep initial energy penalty upon exit from the binding pocket but a smoother transition to the bulk solvent compared to LGR1406. DRF053 emerged as the weakest binder in this series, with a limiting PMF value of 6.75 kcal/mol. The dissociation profile for DRF053 is characterized by a distinct barrier at around ξ = 0.4 nm, followed by a transient minimum, before plateauing at a lower energy than its counterparts.

Collectively, these simulations establish a binding affinity ranking of LGR1406, followed by N&N1 and finally DRF053, correlating the depth of the PMF wells with the structural stability of the respective ligand-protein interfaces.

The ADME properties calculated for Roscovitine and its derivatives (CR8, DRF053, N&N1, LGR1406) outline an overall favorable drug-likeness profile. All compounds comply with the Lipinski, Ghose, Veber, Egan, and Muegge filters and show no PAINS or Brenk alerts, reducing the risk of non-specific assay interference. Molecular weights fall within 354.45–431.53 g/mol; TPSA spans 87.37–100.78 Å²; and the consensus log Po/w ranges from 2.51 to 3.27. These values are consistent with moderately sized molecules of limited polarity and intermediate lipophilicity, compatible with good passive permeability and, overall, potential oral dosing (Veber/Egan compliance). Differences among analogues are small: CR8 and DRF053 display the highest TPSA (100.78 Ų) and slightly higher log Po/w (∼3.2–3.3), whereas Roscovitine and N&N1 are slightly less polar and less lipophilic (TPSA ∼87–88 Å²; log Po/w ∼2.5–2.8). Such variations may primarily impact permeability, distribution, and metabolic clearance, without indicating clear liabilities. Overall, the profile supports progression to experimental ADME verification (solubility, pKa, metabolic stability in microsomes/S9) and selectivity assessments on CDK targets.

However, extended profiling reveals distinct pharmacokinetic behaviors. In terms of metabolism, the series exhibits a high propensity for Cytochrome P450 interactions; all ligands are predicted inhibitors of the major isoforms CYP1A2 and CYP3A4, suggesting a potential risk for drug-drug interactions. Notably, N&N1 displays a slightly cleaner toxicity profile by sparing CYP2C19 and CYP2D6, whereas CR8 and DRF053 exhibit broader inhibitory activity across multiple isoenzymes. Furthermore, despite favorable lipophilicity, the specific permeation model predicts limited Blood-Brain Barrier (BBB) crossing for these derivatives in their current form. Consequently, while the structural properties support progression to experimental ADME verification (solubility, pKa, metabolic stability), the predicted CYP inhibition and CNS distribution profiles indicate that future development may require specific formulation strategies or structural refinement to optimize brain penetrance and minimize metabolic interference.