Roscovitine derivative optimization via molecular dynamics reveals potent and selective CDK5 binders

1 Introduction

Cyclin dependent kinases (CDKs) are a large family of serine/threonine kinases comprising 20 different members (25 if including CDK-like) that play a crucial role in the regulation of the cell cycle and other related cellular processes (Malumbres et al., 2009). Among them CDK1, CDK2, CDK3, CDK4, CDK6 and CDK7 are major regulators of cell cycle progression, while others are involved in a less direct way (Figure 1).

Cyclin dependent kinase 5 (CDK5) is an atypical member of the family, although it shares the basic structure and highly conserved sequences, especially in the catalytic site that binds to ATP (Echalier et al., 2010).

CDK5 does not directly participate in the regulatory progression of cell cycle control, nor is it activated by cyclins, but is activated by binding with non-cyclin activators named CDK5R1 (p35) and CDK5R2 (p39), or their respective truncations p25 and p29. Both the activators and the truncations are found abundantly in post-mitotic neurons, which led to defining CDK5 as a neuron-specific CDK (Cicero and Herrup, 2005).

Indeed, CDK5 shows multiple functions in neuronal tissues. It is involved in regulating neuronal survival, migration, neurite outgrowth, secretion, dopamine signaling, and cytoskeletal dynamics and participates in the pathological changes in neurodegenerative diseases (Dhavan and Tsai, 2001). However, there is growing evidence of its roles outside the nervous system. In several different cell types, CDK5 has been shown to regulate crucial cellular processes such as gene expression, cell migration, apoptosis and DNA damage response (DDR). Dysregulation of its expression and/or activity triggers pathological processes including cancer, senescence, diabetes, immune dysfunction and inflammation (Arif, 2012; Contreras-Vallejos et al., 2012). Coherently with these observations, CDK5 is generally dysregulated in various types of cancer and is linked with cancerous characteristics and prognosis (Contreras-Vallejos et al., 2012; Liu et al., 2017), making it a novel biomarker and promising therapeutic target in cancer treatment. From this point of view, one of the most striking CDK5 features is its activity in promoting cancer metastatic processes and chemo-resistance (Liu et al., 2017; Yuzhalin et al., 2024). This is in major part due to the protective effect that CDK5 exerts on PD-L1 which is a major actor in these processes (Kciuk et al., 2023; Liu et al., 2022). On one hand, CDK5 inhibits posttranslational modification of IRF2BP2, which increases IRF2 and IRF2BP2 abundance and sustains PD-L1 transcriptional repression after IFN-γ stimulation (Dorand et al., 2016); on the other hand, its activity leads to degradation of FBXO22, a member of the F-box family, subunits of the ubiquitin E3 ligase complex SCF (SKP1-cullin-F-box), which exhibits specificity for PD-L1 (De et al., 2021). When overexpressed, F-box proteins can function as tumor suppressors if their substrates are oncoproteins or as oncoproteins if their substrates are tumor suppressors (Johmura et al., 2020). FBOX22 expression and phosphorylation dependent activation are necessary to decrease the PD-L1 levels in cancer cells (De et al., 2021).

Conversely, CDK5-mediated phosphorylation of FBXO22 itself, or of a protein that mediates its stability, decreases the levels of FBXO22 and therefore increases the stability of its substrate PD-L1. Indeed, CDK5 disruption or catalytic inhibition attenuates tumor PD-L1 expression and promotes antitumor immunity (De et al., 2021). A main problem with the use of CDK inhibitors is that, in contrast to CDK5 inhibition, which decreases PD-L1 levels, inhibiting other CDKs, including CDK4/6, increases PD-L1 levels by disrupting the cullin3-SPOP pathway that mediates PD-L1 ubiquitination and degradation (Zhang et al., 2018). This can be a major problem in the chemotherapeutic use of CDKs inhibitors, as shown by the frequent choice of a combination therapy of CDK inhibitors and PD1-PDL1 specific antibodies (Zhang et al., 2021). It is therefore very relevant to select inhibitors with higher affinity for CDK5 as compared with the other CDKs, in particular CDK4/CDK6. In the large library of CDK inhibitors recently developed, the compounds that have the best potential in this direction are Roscovitine and its derivatives. Roscovitine [CY-202, (R)-Roscovitine, (Seliciclib), is a small molecule that acts through direct competition at the ATP-binding site. It is a broad-range purine inhibitor, which inhibits CDK1, CDK2, CDK5 and CDK7, but is a poor inhibitor for CDK4 and CDK6 (Cicenas et al., 2015). In in vitro experiments, it showed the lowest IC50 for CDK5 as compared with the other CDK tested (Cicenas et al., 2015; Meijer et al., 1997). Roscovitine is widely used as a biological tool in cell cycle, cancer, apoptosis and neurobiology studies (Meijer et al., 1997). Moreover, it is currently evaluated as a potential drug to treat cancers, neurodegenerative diseases, inflammation, viral infections, polycystic kidney disease and glomerulonephritis (Cicenas et al., 2015). Several Roscovitine derivatives were developed (Kabadi et al., 2012; Oumata et al., 2008; Bettayeb et al., 2008; Sroka et al., 2010) that could be even more selective for CDK5 but they were never evaluated for this aspect. From a computational standpoint, molecular dynamics (MD) simulations of large, solvated macromolecules (Castrignanò et al., 2000; Chillemi et al., 2001; Castrignanò et al., 2002; Gabellone et al., 2022) provide a powerful means to investigate complex systems at atomic resolution. Combining molecular docking with MD simulations and free energy calculation of ligand-receptor complexes enables a detailed, time-resolved analysis of their interactions (Cirigliano et al., 2016; Keretsu et al., 2020; Bò et al., 2021; Madeddu et al., 2022). Recent studies further highlight the value of integrating docking with extended molecular dynamics and energetic/collective-motion analyses to refine binding hypotheses beyond static docking scores. For instance, Padhi et al. (2021) investigated arbidol binding at the SARS-CoV-2 RBD/ACE2 interface by combining cavity-guided docking with 300 ns MD simulations and subsequent trajectory analyses to characterize stability and binding mechanisms. Similarly, Bhagat et al. (2025) adopted an integrative multiscale strategy (structure modeling, docking, coarse-grained and all-atom MD, complex simulations) coupled with PCA and free-energy landscape analyses to connect structural/dynamic perturbations to changes in substrate binding in human thymidine phosphorylase variants. These examples support the use of long-timescale MD combined with post-processing analyses as a robust approach to investigate protein-ligand recognition and stability. Alongside these more traditional computational approaches, recent advances in machine learning have introduced powerful models for assessing tumor sensitivity to therapeutic compounds (Bolis et al., 2017).

In this study, we integrate molecular docking, explicit-solvent molecular-dynamics (MD) simulations, MM/GBSA binding-free-energy calculations, principal-component analysis (PCA), and time-resolved ligand-receptor interaction profiling to benchmark four Roscovitine derivatives, i.e. CR8, DRF053, N&N1, and LGR1406 (Figure 2), against the ATP-binding sites of CDK2, CDK4, CDK5, CDK6, and CDK12, using ATP and the Roscovitine as reference.

This analysis aims to identify compounds that exhibit high affinity for CDK5 (and possibly CDK2) together with optimal discrimination against CDK4/CDK6, and to determine the chemical determinants, such as key hydrogen bonds, hydrophobic contacts, and steric or electrostatic clashes, that can guide the creation of more selective and potent compounds. These insights will support the rational design of next-generation analogues with enhanced efficacy and specificity. As this work is purely computational, the proposed affinity/selectivity trends should be regarded as hypothesis-generating and will require experimental validation.

2 Materials and methods

All primary simulations presented in this study were performed on a high-performance computing (HPC) infrastructure, which allowed us to generate multiple complex trajectories and to efficiently handle the large data volumes produced. Specifically, computations were carried out on the TIER-0 Leonardo supercomputing system (https://leonardo-supercomputer.cineca.eu/). To promote reproducibility of the computational workflow (Supplementary Figure S1) and facilitate reuse of the outputs for downstream analyses by other groups, we created a dedicated figshare project and deposited the results generated at each stage of the analysis (see Table 1, Datadir 1–7).

TABLE 1. Overview of the data files generated in this study and their availability on Figshare.

2.1 Preparation of inhibitor structures

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.

2.2 Molecular docking

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.

2.3 Molecular dynamics simulations

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).

2.4 MMPBSA binding free energy and principal component analysis

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.

2.5 Interaction 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.

2.6 Umbrella sampling

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.

2.7 Drug-likeness and ADME filters

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.

3 Result

The bioinformatic workflow employed in this study, and described in Supplementary Figure S1, follows a hierarchical approach to investigate the interactions between various Cyclin-Dependent Kinases (CDKs) and specific ligands. It integrates molecular docking, molecular dynamics simulations, binding free energy calculations, principal component analysis (PCA), and detection of specific ligand-residue interactions. All the primary simulations described in this study were run on a high-performance computing (HPC) infrastructure, enabling the generation of multiple complex trajectories and the efficient management of the large volumes of data produced. To ensure repeatability of the computational experiment and reuse of the data for downstream analysis by other interested groups, we created a project on figshare where we deposited the results obtained at each step of the pipeline (see Table 1, Datadir 1–7).

3.1 Molecular docking results

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.

TABLE 2. Best binding stability (kcal/mol) from molecular docking simulations of CDK family proteins with the listed ligands.

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.

3.2 Analysis of RMSD data across CDKs and ligands

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.

TABLE 3. RMSD values (nm) for the c-∝ of the protein of the CDK-ligand complexes and the Apoprotein over 1 μs molecular dynamics simulations, with analysis starting at 0.2 μs. RMSD values were computed using the gmx rms module in GROMACS and are reported as mean ± standard deviation (in nm), based on statistical analysis of the trajectory.

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.

TABLE 4. Root mean square deviation (RMSD) values for ligands in CDK-ligand complexes over 1 μs molecular dynamics simulations, with analysis starting at 0.2 μs. RMSD values were calculated using the gmx rms module in GROMACS and are reported as mean ± standard deviation (in nm), based on statistical analysis of the trajectory.

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.

3.3 Residue-specific flexibility of CDK-ligand complexes

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.

3.4 Gibbs free binding energy via MM/PBSA

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.

TABLE 5. Gibbs free binding energy values (ΔG) obtained from MM/PBSA calculations.

Estimated binding free energies were computed using the Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) method on 8,000 frames extracted at 10 ps intervals from 1 μs molecular dynamics trajectories, starting at 0.2 μs. The table reports results for CDK-ligand complexes, expressed as mean ± standard deviation in kilocalories per mole (kcal/mol), based on trajectory-derived energy components.

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.

3.5 PCA analysis

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.

3.6 Interaction analysis

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.

TABLE 6. Persistent protein-ligand contacts identified through molecular dynamics simulations.

For each cyclin-dependent kinase (CDK2, CDK4, CDK5, CDK6, and CDK12), the listed residues form stable side-chain contacts (interaction occupancy ≥30% over the final 800 ns of 1 µs trajectories) with six ligands: the natural substrate ATP, Roscovitine and four inhibitor derivatives off the latter (CR8, DRF053, LGR1406, N&N1). Data were extracted using the MD-Ligand-Receptor workflow from thirty independent replica simulations. Conserved hinge anchors (Ile10, Phe80, and Lys33) are highlighted in bold, isoform-specific selectivity hotspots (e.g., Asp86/Gln130 in CDK5 or Asp158 in CDK4) are italicized, and within each cell, hydrophobic contacts are listed before polar/charged interactions. This layout underscores how CR8 and LGR1406 engage both the canonical hinge region and additional polar residues in CDK2 and CDK5, while avoiding the acidic Asp158/Asp163 pocket targeted by ATP in CDK4/6, thus providing structural insights into the observed affinity and selectivity trends.

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.

3.7 Evaluation of binding free energies via umbrella sampling

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.

3.8 ADME and drug-likeness assessment of roscovitine and analogues

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.

4 Discussion

CDKs emerged as one of the most relevant targets for tumor chemotherapy (Zhang et al., 2021). CDK5 is an atypical member of the CDKs family that is not directly involved in cell cycle regulation but has a very relevant role in promoting DNA damage response and immune response evasion in cancer. This is partly due to the protective effect that CDK5 exerts on PD-L1 which is a major actor in these processes. In this regard, CDK5 conflicts with other CDKs, namely CDK4 and CDK6, which instead promote PDL-1 degradation.

Targeting CDK5 with small molecules which could efficiently inhibit other CDKs, such as CDK2, but poorly affect CDK4 and CDK6, could be an extremely convenient strategy since it would block cell cycle progression and, at the same time, would repress PDL-1-induced immune-evasion and chemo-resistance (15). Indeed, CDK5 targeting has also been recently considered for killing of breast cancer cells via perturbation of redox homeostasis (Navaneetha Krishnan et al., 2018). Moreover, CDK5 seems involved in brain lesions in models of global and focal cerebral ischemia (Le Roy et al., 2021) and CDK5-selective inhibitors could be very useful in treating related pathologies without triggering severe collateral effects related to the inhibition of other CDKs.

Among CDKs inhibitors, Roscovitine has been shown to have some selectivity toward CDK5, as compared with CDK4 and CDK6 in vitro. We analyzed the structural interaction of Roscovitine and some of its derivatives, namely CR8, DFR053, N&N1, and LGR1406, with different CDKs by molecular docking and dynamics in order to identify those showing a predicted high affinity for CDK5 (and eventually for CDK2), accompanied by optimal discrimination against CDK4/CDK6 and to derive chemical clues to eventually produce other, more selective and efficient compounds.

From molecular docking data, LGR1406 and CR8 exhibited the greatest binding affinity for CDK5 and LGR1406 showed the best discrimination between CDK5 and CDK4/CDK6 (Table 2).

The analysis of root mean square deviation (RMSD) provides insights into the structural stability of CDK2, CDK4, CDK5, CDK6, and CDK12 complexes with different ligands over the simulation period (1,000 ns). The data (Tables 3, 4) suggested an optimal stability of Roscovitine, LGR1406 and N&N1 complexed with CDK5. Once again, the stability of LGR1406 complexed with CDK4 and CDK6 was lower.

According to the Gibbs Free energy estimates of molecular dynamics, all four derivatives outperform Roscovitine (ΔG = −19.37 ± 4.19) in binding to CDK5. The main differences observed with the docking based estimates are represented by the case of N&N1 which in the docking estimate showed the same affinity as the precursor and by the larger energy difference between Roscovitine and all four derivatives observed in the molecular dynamics analysis. Once again, the best affinity was observed in the case of LGR1406 complexed with CDK5 (ΔG = −29.63 Kcal/mol), followed by DFR053 (ΔG = −26.83 Kcal/mol) and N&N1 (ΔG = −26.69 Kcal/mol).

Notably, Umbrella Sampling analysis provided an independent estimate of the unbinding thermodynamics for the three lead compounds on CDK5, confirming the same affinity hierarchy observed across the other computational metrics. In particular, the PMF profiles indicate a markedly higher dissociation free energy for LGR1406 (∼12.27 kcal/mol) compared with N&N1 (∼8.24 kcal/mol) and DRF053 (∼6.75 kcal/mol), consistent with a more strongly anchored binding mode and higher desolvation penalties during unbinding.

As depicted in Figure 15, increase in the affinity of LGR1406, N&N1 and DRF056, can be due to stronger interactions with the key amino acid residues of CDK5. In particular, heterobicyclic pyrazolopyrimidine (LGR1406), pyrroletriazine (N&N1) and purine (DRF056) scaffolds of these three ligands are oriented differently in the enzyme binding site than the purine core of Roscovitine.

Comparing Roscovitine with LGR1406, it can be seen that the pyrazole ring of the latter is on the opposite side with respect to the imidazole portion of the former (Figure 15a). Therefore, nitrogen atoms of pyrazole make two hydrogen bonds with Cys83 of CDK5. In N&N1 one nitrogen of the five-membered heterocycle and the exocyclic amine group interact with Cys 83 (Figure 15b) as well as in DRF056 imidazole and external amine make hydrogen bonds with the Cys 83 residue. On the contrary, Roscovitine does not interact with the above mentioned cystine.

Additionally, LGR1406 also exhibited clear discrimination versus CDK4 (ΔΔG = 5.54 Kcal/mol) and CDK6 (ΔΔG = 5.94 Kcal/mol) (Table 1). Similar results were obtained for N&N1 and DFR053. The promising in silico ADME and drug-likeness profiles (Table 7) further bolster the potential of these compounds, particularly LGR1406 and N&N1, to progress as viable drug candidates, as their physicochemical properties are conducive to good bioavailability and reduced off-target liabilities. From this analysis we can conclude that LGR1406 could be the best choice to selectively inhibit CDK5 protective action on PDL-1, followed by N&N1. This could be a good starting point for an experimental analysis of their effects in living cells.

TABLE 7. In silico drug-likeness and early ADME profiling of CDK ligands, Roscovitine and analogues (CR8, DRF053, N&N1, LGR1406).

The table reports compliance/violations for Lipinski, Ghose, Veber, Egan, and Muegge filters; predicted Blood-Brain Barrier (BBB) permeability, inhibition profiles for major Cytochrome P450 (CYP) isoforms (CYP1A2, CYP2C19, CYP2C9, CYP2D6, and CYP3A4); PAINS and Brenk structural alerts; molecular weight (MW), topological polar surface area (TPSA), and consensus log Po/w.

A key limitation of this work is the absence of experimental corroboration. Although these in silico metrics are promising, metabolic stability assessments (e.g., microsomal intrinsic clearance) are not covered by this computational workflow and will be required in the validation phase. Furthermore microsecond-scale MD, MM/PBSA estimates, and persistent interaction fingerprints provide a mechanistic rationale to prioritize compounds (notably LGR1406), these approaches cannot substitute for direct measurements of enzymatic inhibition, binding thermodynamics/kinetics, or cellular target engagement. Accordingly, our conclusions are framed as a computational prioritization of candidate CDK5 inhibitors to guide subsequent wet-lab testing.

5 Conclusion

This computational study demonstrates that rational optimization of Roscovitine derivatives can yield highly selective CDK5 inhibitors with favorable drug-like properties. Because this is a purely computational study, the reported affinity/selectivity rankings remain predictions that require experimental confirmation. Through integrated molecular docking, microsecond-scale molecular dynamics simulations, MM/PBSA free energy calculations, and comprehensive interaction profiling, we identified LGR1406 as the most promising candidate, exhibiting superior binding affinity for CDK5 (ΔG = −29.63 kcal/mol) coupled with significant discrimination against CDK4 (ΔΔG = 5.54 kcal/mol) and CDK6 (ΔΔG = 5.94 kcal/mol). Consistently, the Umbrella Sampling/PMF analysis of ligand dissociation from CDK5 identified LGR1406 as the most stable complex, showing the highest dissociation penalty relative to N&N1 and DRF053, thereby further supporting its prioritization as the lead CDK5 inhibitor in this series. The enhanced selectivity profile of LGR1406, N&N1, and DRF053 stems from their ability to engage both the conserved Ile10-Phe80 hinge anchor and the CDK5-specific Cys83-Asp86-Gln130 polar network, while avoiding the acidic Asp158/163 pocket preferentially occupied by ATP in CDK4/6. Notably, the pyrazolopyrimidine scaffold of LGR1406 adopts an inverted orientation relative to Roscovitine’s purine core, enabling dual hydrogen bonding with Cys83 that is absent in the parent compound. These structural insights, combined with favorable ADME predictions and compliance with all major drug-likeness filters, establish LGR1406 and N&N1 as lead compounds warranting experimental validation. Given that selective CDK5 inhibition can simultaneously suppress tumor PD-L1 expression while avoiding the PD-L1-stabilizing effects associated with CDK4/6 inhibition, these compounds represent a promising therapeutic strategy to overcome immune evasion and chemoresistance in cancer. Future work should focus on in vitro enzymatic assays, cellular efficacy studies, and structure-based optimization to further enhance potency and selectivity, potentially leading to next-generation CDK5-targeted therapies for oncology and neurodegenerative disease.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchbi.2026.1741100/full#supplementary-material

Schematic of the study workflow. Structures of CDK2, CDK5, CDK4, CDK6, and CDK12 (PDB IDs: 2WIH, 1UNL, 7SJ3, 5L2S, 7NXK) were used, with Robetta-based homology modeling when needed. All structures were relaxed in GROMACS. Docking (AutoDock Vina, 25 runs/system) was performed with Roscovitine, CR8, DRF053, N&N1, LGR1406, and ATP. Complexes and Apoprotein underwent 1000 ns GROMACS MD simulations. Analyses included ligand-receptor interactions, PCA, and MM/PBSA. Supplementary material at the source

Principal Component Analysis (PCA) of ligand dynamics in complex with CDKs. The plots show the projection of the conformational ensembles of each ligand onto the first two principal components (PC1 and PC2), derived from molecular dynamics simulations. Each panel corresponds to a specific CDK: (a) CDK2, (b) CDK4, (c) CDK5, (d) CDK6, and (e) CDK12. The PCA highlights the dominant modes of motion of the ligands within the binding pocket, allowing comparative evaluation of conformational variability across different CDK-ligand complexes. Supplementary material at the source

The plots display the frequency counts for the 50 sampling windows used to construct the Potential of Mean Force (PMF) profiles. The significant overlap between adjacent windows along the entire reaction coordinate confirms adequate sampling and continuous coverage for the Weighted Histogram Analysis Method (WHAM). The panels correspond to the three investigated systems: (a) LGR1406-CDK5, (b) DRF053-CDK5, and (c) N&N1-CDK5. The color gradient (blue to red) represents the progression of the reaction coordinate from the bound state to the fully dissociated state. Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Footnotes

Edited by: Vikash Kumar Dubey, Indian Institute of Technology (BHU), India

Reviewed by: Aditya Kumar Padhi, Indian Institute of Technology (BHU), India

Narra Siva Krishna, VNR Vignana Jyothi Institute of Engineering and Technology, India