Design and Evaluation of New 6-Trifluoromethoxy-Isatin Derivatives as Potential CDK2 Inhibitors

TL;DR

This study designs and evaluates new isatin compounds with a trifluoromethoxy group as potential CDK2 inhibitors for cancer treatment.

Contribution

The paper introduces novel 6-trifluoromethoxy-isatin derivatives and evaluates their CDK2 inhibition potential using in silico methods.

Findings

The derivatives form stable complexes with CDK2 through hydrogen bonds and hydrophobic interactions.

The trifluoromethoxy group enhances ligand orientation and insertion into the CDK2 hydrophobic pocket.

ADMET analysis shows favorable absorption and toxicity profiles with moderate solubility issues.

Abstract

Cyclin-dependent kinase 2 (CDK2) plays a central role in cell cycle regulation and represents an important molecular target in anticancer drug development. In this study, a series of novel isatin derivatives substituted with a trifluoromethoxy group at the C6 position were designed and evaluated as potential CDK2 inhibitors using a comprehensive in silico approach. Density functional theory calculations were applied to analyze the electronic properties of the proposed compounds. Molecular docking and molecular dynamics simulations were used to investigate binding modes, conformational stability, and key interactions within the CDK2 active site. Binding free energies were estimated using the Molecular Mechanics Poisson–Boltzmann Surface Area (MMPBSA) method, while QSAR-based (Quantitative Structure–Activity Relationship) ADMET (Absorption, Distribution, Metabolism, Excretion, and Toxicity) analyses were performed to assess drug-likeness and pharmacokinetic profiles. The results indicate that the investigated derivatives form stable complexes with CDK2, supported by persistent hydrogen bonds in the hinge region and favorable hydrophobic interactions. The trifluoromethoxy substituent significantly affects ligand orientation and promotes deeper insertion into the hydrophobic pocket compared with previously studied isatin analogues. ADMET predictions suggest generally favorable absorption and toxicity profiles, with moderate solubility limitations. Overall, these findings support the potential of 6-trifluoromethoxy-isatin derivatives as promising CDK2 inhibitors and provide a basis for further experimental studies.

Figures

Full text

1. Introduction

Cancer continues to represent a substantial global health burden and remains one of the leading causes of death worldwide [1,2]. Despite significant advances in diagnostic techniques and the development of new therapeutic strategies, effective control and prevention of this group of diseases remain limited. A defining feature of cancer is uncontrolled cellular proliferation, which arises from disruptions in the tightly regulated mechanisms governing the cell cycle. Central to the regulation of these processes are enzymes belonging to the cyclin-dependent kinase family, which are involved in the control of cell cycle progression, transcriptional regulation, and DNA synthesis [3,4,5]. Among these enzymes, cyclin-dependent kinase 2 (CDK2) plays a pivotal role in regulating the transition from the G1 to the S phase, as well as progression through the S phase. Numerous studies have demonstrated that dysregulation, hyperactivation, or overexpression of CDK2 is frequently observed across a wide range of cancer types and is strongly linked to excessive proliferation, impaired apoptotic signaling, and genomic instability [6,7,8,9,10]. These features make CDK2 a compelling molecular target in anticancer drug development.

Although several CDK inhibitors have successfully entered clinical use, their application is often constrained by insufficient selectivity, off-target effects, and toxicity-related dose limitations [7,11,12,13]. The high degree of structural conservation within the ATP-binding pockets of CDK family members further complicates the development of selective inhibitors. As a result, there is a continuing need to identify new chemical frameworks capable of achieving stable, selective, and effective interactions with the CDK2 active site.

In this context, isatin (1H-indole-2,3-dione) has emerged as a valuable and pharmacologically relevant scaffold in medicinal chemistry. The oxindole core of isatin is particularly well suited for kinase inhibition due to its ability to establish key interactions within the ATP-binding region, especially in the hinge area of CDK2 defined by residues such as GLU81 and LEU83 [14]. Previous investigations have shown that strategic structural modifications of the isatin framework—most notably substitutions on the aromatic ring and functionalization at the C-3 position—can substantially improve binding affinity, inhibitory potency, and overall stability of ligand–kinase complexes [15,16,17,18].

Recent studies further underscore the significance of isatin-derived compounds in the development of cyclin-dependent kinase inhibitors, with particular emphasis on CDK2. Comprehensive structure–activity relationship analyses have demonstrated that isatin-based scaffolds, especially those functionalized at the C-3 position, can achieve high inhibitory potency and selectivity toward CDK2 through effective hinge-region interactions and stabilization within the ATP-binding pocket [19]. Moreover, isatin derivatives bearing hydrazide or hydrazone side chains—representing close structural analogues of the compounds discussed herein—have been reported to exhibit pronounced antiproliferative activity and kinase inhibition profiles [20,21,22,23,24]. The presence of hydrazide-linked moieties was shown to enhance hydrogen-bonding networks and conformational flexibility, facilitating improved accommodation within the CDK2 catalytic cleft. Collectively, these findings highlight isatin-based frameworks as a versatile and pharmacologically relevant platform for the rational design of CDK2 inhibitors and provide strong justification for further exploration of structurally related derivatives.

Our research group has previously focused on the rational design of isatin-based CDK2 inhibitors through a combination of computational and experimental methodologies. Approaches including quantum-mechanical calculations, molecular docking, molecular dynamics simulations, and ADMET prediction have been applied to optimize molecular properties relevant to both activity and pharmacokinetics. These studies highlighted the critical influence of substituent type and position on ligand orientation, hydrogen-bond formation, and penetration into the hydrophobic pocket of CDK2. In particular, compounds bearing electron-withdrawing groups or moieties capable of forming additional stabilizing interactions exhibited superior binding characteristics compared with known reference inhibitors [15,17,25].

Extending these findings, further diversification of the isatin scaffold appears to be a promising strategy toward improved CDK2 inhibition. Previous studies on isatin derivatives bearing nitro or methyl substituents at the 5-position [15,17] demonstrated that such modifications play a key role in modulating interactions of the isatin core by enforcing an alternative orientation relative to the active site and promoting additional intermolecular contacts. In this context, particular attention has been directed toward fluorinated substituents, as fluorine-containing groups are well recognized for their ability to modulate electronic distribution, lipophilicity, metabolic stability, and intermolecular interactions. Accordingly, incorporation of a trifluoromethoxy group at the 6-position of the isatin core may exert an even more pronounced effect, contributing more substantially to ligand orientation within the CDK2 binding pocket and to the formation of stabilizing interactions, thereby potentially enhancing inhibitory activity.

Based on this rationale, the present study focuses on the design and in silico evaluation of a new series of 6-trifluoromethoxy isatin derivatives as potential CDK2 inhibitors. An integrated computational strategy encompassing molecular docking, molecular dynamics simulations, MMPBSA binding free energy calculations, QSAR analysis, and ADMET profiling was employed to investigate binding affinity, structural stability, and drug-likeness of the proposed compounds. The obtained results provide deeper insight into the structure–activity relationships of isatin-based CDK2 inhibitors and further support the suitability of this scaffold for the development of novel anticancer therapeutics.

2. Results and Discussion

The group of compounds investigated in the present study belongs to the class of isatin derivatives, for which significant inhibitory activity against cyclin-dependent kinases (CDKs) has been repeatedly reported. The isatin molecular scaffold offers broad opportunities for the development of new classes of compounds through modifications introduced directly into the core structure, particularly at the fifth and sixth positions, as well as through reactions enabling the attachment of new moieties within the pyrrolidinone ring, for example via reactions with benzoylhydrazide derivatives. This approach, previously described in our earlier work [17], allows for the design of compounds with an improved fit within the active-site cavity and enhanced binding affinity, ultimately resulting in increased inhibitory potential. The model structure illustrating the geometry of the new group of inhibitors is shown in Figure 1, with labels R2, R3, and R4 indicating the potential positions of substituents on the benzoylhydrazide aromatic system. Based on commercially available substituted benzoylhydrazide derivatives, 34 new compounds were designed, and their inhibitory and pharmacological potential was evaluated using a wide range of computational chemistry methods.

2.1. Frontier Orbitals Analysis

The electronic properties of the designed group of compounds were investigated using DFT methods. The frontier orbital energies and derived parameters characterizing the potential reactivity of the studied inhibitors are presented in Table 1 and Table S2. The HOMO energies range from −6.5 eV (compounds 1, 2a, 2c, 2i, and the series of derivatives with substituents at position 3) to −5.9 eV (2f), while derivatives with substituents at position 4 exhibit intermediate values ranging from −6.0 eV to −6.2 eV. The higher HOMO energies observed for compounds 2f and 4f indicate an increased ability to donate electron density, which favors donor interactions with electrophilic active groups located in the CDK2 active site, such as LYS33 and LEU83. The LUMO energies fall within a narrow range of −2.4 to −2.6 eV, with the lowest values observed for structures 3h, 3k, 4b, and 4d, suggesting their potentially greater electron-accepting ability and slightly enhanced capacity to stabilize charge-transfer interactions with nucleophilic groups present in amino acid residues such as ASP145, GLU51, and GLU81.

The localization of the HOMO and LUMO orbitals of the studied compounds is presented in Figure S1 (Supplementary Materials). For most of the compounds, the localization of both types of orbitals exhibits a normalized and consistent distribution. Exceptions are derivatives containing strongly electron-donating substituents, such as –NH_2_ (2f, 4f) or mie–OH (2i, 4i), which significantly contribute to enhanced reactivity of the benzoylhydrazide moiety.

The energy gap (ΔEGAP) is a key descriptor reflecting the electronic polarizability and intramolecular charge-transfer capability of the studied compounds. Within the investigated ligand set, ΔE_GAP_ values range from 3.53 to 4.04 eV. Higher values, observed for structures 1, 2a, 2c, 2i, and 3a–3e, indicate slightly lower polarizability and greater kinetic stability, features that may favor more selective and structurally well-defined inhibitory binding. Conversely, derivatives 2f, 3k, and compounds from series 4 exhibit lower ΔE_GAP_ values, suggesting enhanced electronic flexibility and a greater propensity for polarization-driven interactions within the enzyme active site, including π–π contacts with aromatic residues such as PHE80 and PHE82.

The reported values of chemical hardness (η) fall within a relatively narrow range of 1.77–2.02 eV. Lower hardness values, accompanied by higher softness values (σ = 0.55, S ≈ 0.27), indicate that structures such as 2f, 3k, and compounds from series 4 may exhibit a slightly greater ability to adapt to the properties of the investigated active site.

The absolute electronegativity (χ) determined for the studied compounds ranges from 4.18 to 4.57 eV. Compounds with higher absolute electronegativity values exhibit a stronger electron-accepting character, which may represent a favorable factor contributing to stabilization of the inhibitor–active-site complex through stronger electrostatic interactions and charge transfer. The global electrophilicity index (ω) reaches the highest values for ligands 3h (5.30 eV), 3b–3e (5.11 eV), and 2c (5.01 eV). These values predispose these ligands to greater efficiency in interactions with nucleophilic amino acid residues within the CDK2 active site, potentially enhancing the stability of the formed complexes. Lower values of this parameter observed for ligands 2i (4.89 eV) and 4f (4.84 eV) indicate a more balanced distribution of donor–acceptor properties and potentially more diverse modes of interaction with the active site.

2.2. Molecular Docking

For the purposes of this study, 34 novel 6-trifluoromethoxyisatin derivatives were designed, including the native form and 33 derivatives bearing substituents within the benzoylhydrazide ring. Based on these compounds, a molecular docking procedure was performed to obtain the most representative conformations of complexes formed between the selected inhibitors and the active site of cyclin-dependent kinase 2 (CDK2). A representative complex formed by the native derivative is shown in Figure 2. The ligand core is tightly embedded within the hinge region of the protein, in close proximity to GLU81 and LEU83. Notably, the trifluoromethoxy group also participates in binding interactions, with fluorine atoms forming hydrogen bonds with ASP145 and LYS33. An important stabilizing factor of the resulting complex is provided by hydrophobic interactions, in which aromatic systems of amino acid residues located in this region (PHE80, PHE82), as well as extended aliphatic groups such as that of ILE10, play a significant role.

The chemical affinity and inhibition constants (IC) describing the binding capabilities of each considered ligand toward the CDK2 active site are presented in Table 2. The IC values were obtained based on the Van’t Hoff Equation (1).

$$[eqn]$$

The presented dataset allows for the characterization of changes in the binding capabilities of the studied ligands in relation to the presence and nature of substituents located at the R2, R3, and R4 positions. The binding affinities of the modified derivatives were evaluated relative to the reference value of −9.80 kcal/mol obtained for the native molecule lacking any substituents on the aromatic system.

For modifications located at the R2 position, the greatest diversity of observed effects is evident. The introduction of sterically bulky substituents, which may constitute significant steric hindrance—as observed for derivatives 2g, 2h, 2i, 2j, and 2d—results in a decrease in binding affinity, with the most pronounced reduction exceeding 1.5 kcal/mol. In contrast, derivatives 2a, 2c, 2f, and 2i exhibit improved binding capabilities; however, the observed increases remain within the range of −0.1 to −0.2 kcal/mol.

Within the group of derivatives bearing substituents at the R3 position, no structure was identified in which the presence of a modification negatively affected binding affinity relative to the native compound. An increase in affinity was observed for derivatives 3a, 3b, 3c, 3d, 3e, 3h, and 3k, with gains ranging from −0.1 to −0.2 kcal/mol.

The final group comprises derivatives modified at the R4 position. Within this population, three compounds exhibit reduced binding affinity compared to the native molecule, while seven derivatives show enhanced binding activity. Consistent with the trends observed for the other substitution patterns, the enhancements in binding affinity within this group remain modest and do not exceed −0.2 kcal/mol.

Analysis of the effects of individual substituents indicates that the introduction of groups such as –CH_3_, –CF_3_, and –F significantly enhances the binding capabilities of the investigated derivatives, regardless of their position within the benzoylhydrazide ring.

Within the investigated group, 18 derivatives can be identified for which the introduction of structural modifications resulted in enhanced binding affinity toward the active site of cyclin-dependent kinase 2 (CDK2). The properties of these complexes were analyzed in greater detail relative to those of the system formed by the native molecule. Table 3 presents parameters describing hydrogen-bond lengths as well as distances between aromatic systems involved in hydrophobic interactions.

The most critical interactions in the context of previous studies on isatin derivatives binding within the hinge region of CDK2 are hydrogen bonds formed with GLU81 and LEU83. For the investigated derivatives, the hydrogen bond involving the hydrogen atom of LEU83 falls within a narrow range of 1.73–1.90 Å, classifying these interactions as strong hydrogen bonds. Another key interaction involves the oxygen atom of GLU81; for most of the studied derivatives, the corresponding bond length is approximately 2.63 Å. An exception is derivative 2a, for which a slightly shorter distance of 2.36 Å was observed. These values indicate the presence of weak hydrogen bonds.

An additional interaction, not involving the isatin core, engages a hydrogen atom from the hydrazine moiety of the side chain and the oxygen atom of LEU83, with observed distances ranging from 2.25 to 2.59 Å. Another group of hydrogen bonds was identified for fluorine atoms of the trifluoromethoxy group located at the 6-position of the isatin core. This group participates in the formation of two hydrogen bonds: the first with ASP145, which for most of the analyzed systems exhibits distances between 1.76 and 1.78 Å, with exceptions observed for derivatives 2a (1.69 Å) and 3b (2.23 Å). The second interaction involves hydrogen atoms of the amino group of LYS33, with bond lengths predominantly ranging from 2.38 to 2.52 Å; an exception is structure 1, for which a distance of 3.01 Å was measured, representing a borderline value for hydrogen bonding.

Another class of interactions that is particularly important for stabilizing ligand–active-site complexes comprises hydrophobic interactions. These were quantified by measuring distances between aromatic systems capable of forming so-called stacking interactions, characteristic of aromatic moieties. In most cases, the interaction between the isatin core and the aromatic system of the PHE80 residue is described by distances slightly exceeding 4 Å. A notably shorter distance was observed only for derivative 2a (3.74 Å), for which the hydrogen-bond lengths involving LEU83 and GLU81 also indicate deeper embedding of this part of the molecule within the hydrophobic pocket compared to the remaining structures.

The final set of parameters related to hydrophobic interactions includes distances defined between the aromatic system of PHE82 and the aromatic moiety of the ligand side chain. The observed values range from 3.69 to 4.23 Å and indicate greater conformational freedom of this ligand fragment, resulting from the architecture of the active site and the nature of substituents present on the aromatic system of the analyzed derivatives. Overall, all reported distances fall within ranges compatible with the occurrence of stacking interactions between the ligand and active-site amino acid residues, while the dispersion of values reflects differences in the strength of these interactions.

A proper evaluation of the inhibitory potential of the proposed compounds requires comparison of the obtained results with data for systems exhibiting confirmed and experimentally validated inhibitory activity. Therefore, 3-((2,6-dichlorobenzylidene)hydrazono)indolin-2-one was selected as the reference molecule (RM). In vitro studies have clearly demonstrated that this compound exhibits inhibitory activity against CDK2. The activity of the reference molecule toward the CDK2 active site was described in a previous study [15], and the structure of the reference compound as well as its complex with the CDK2 active site are presented in Figure 3.

The chemical structure of the reference molecule exhibits significant structural similarity to the studied group of inhibitors. In both cases, isatin served as the starting scaffold for compound design; however, the distinguishing features include a trifluoromethoxy group substituted at position 6 and differences in the nature of the substituent attached to the side chain. Docking studies performed for the reference molecule indicated a markedly lower affinity toward the investigated active site. The obtained binding energy value of −8.62 kcal/mol [15] suggested weaker binding not only in comparison with the most active derivatives (−10.1 kcal/mol) but also relative to the native model structure (−9.8 kcal/mol).

The localization of the isatin core within the active-site space was generally consistent with the conformations observed for the studied group of compounds; however, noticeable differences in its orientation were evident. These differences were reflected in variations in the interaction distances with GLU81 and LYS83, as well as in the distance between PHE80 and the aromatic system of the core. In contrast, the orientation of the aromatic moiety in the side chain differed substantially, with distances measured between this moiety and PHE82 being at least 1 Å longer. Based on the observed differences in binding affinity and the identified interactions, it could be concluded that the proposed compounds may exhibit significantly greater potential for stabilizing interactions within the active site of the biological target.

2.3. Molecular Dynamics

All ligand structures for which an increase in binding affinity toward the active site was observed compared to the native structure were selected for the molecular dynamics simulation stage. One of the most fundamental yet critically important evaluations of systems obtained from molecular dynamics simulations is the assessment of their conformational stability and its variation over the temporal evolution of the studied systems. The analysis of ligand mobility and conformational variability provides a general indication of the correctness of inhibitor fitting within the active site at the docking stage and allows for evaluation of the stability of the formed complexes. Table 4 presents a set of data illustrating the averaged RMSD values and their standard deviations for the analyzed populations, considering both the entire simulation time and the final 60 ns of the trajectories. This distinction was introduced because the analysis of the time-dependent distributions of the studied systems (Figure S2) indicated that the initial 20 ns of the simulations were, in most cases, associated with significant conformational changes manifested by pronounced RMSD fluctuations. A superficial comparison of mean values and standard deviations might suggest that this issue mainly concerns the conformational properties of the active site; however, statistical testing based on one-way analysis of variance (ANOVA) for paired datasets revealed statistically significant differences between the compared populations for eleven ligands. In contrast, for parameters describing the active site, statistical consistency between conformer populations collected over 80 ns and 60 ns was observed in only one case (compound 3a). Therefore, based on these observations, subsequent analyses focusing on interaction characterization and binding affinity were restricted to the final 60 ns of the molecular dynamics simulations.

The analysis of the time evolution of ligand RMSD values, presented in Figure S2, clearly indicates that several groups of ligands exhibiting distinct conformational stability characteristics can be distinguished within the studied population. Some structures display significant conformational changes throughout the entire molecular dynamics production phase, cyclically affecting their structural properties. This behavior is observed for compounds 2c and 3c bearing fluorine substituents, which exhibit the most pronounced and relatively chaotic conformational fluctuations.

More ordered conformational transitions were observed for derivative 3b containing a trifluoromethyl substituent; in this case, the recorded increase in RMSD was associated with a 180° rotation of the aromatic moiety in the side chain. In contrast to derivatives 2c and 3c, this conformational change was a single event, and the adopted conformation remained stable thereafter.

Another group consists of systems that undergo relatively dramatic conformational changes during the initial 20 ns of the simulation but subsequently achieve a stable conformational arrangement within the active-site space, as observed for derivatives 3a, 3h, and 4a.

The final and most numerous group comprises molecules that underwent conformational rearrangements during the initial optimization stage relative to the starting geometry and then exhibited only minor oscillations throughout the molecular dynamics simulations. Representative examples of this behavior include derivatives 1, 2a, 2i, 3e, 4b, 4g, and 4i. Overall, the RMSD values clearly indicate that, over the course of the molecular dynamics simulations, each ligand underwent subtle conformational adjustments aimed at optimizing its fit within the active site.

The stability of complexes formed by ligands with the active sites of the investigated biological targets is governed by the presence of a broad spectrum of interactions, including electrostatic and van der Waals interactions, hydrogen bonds, halogen bonds, as well as a wide range of hydrophobic interactions, including, among others, π–π stacking interactions. In order to unambiguously verify the stability of the studied complexes, a detailed analysis of hydrogen bonds identified at the docking stage as well as those formed during molecular dynamics simulations was performed. A summary of the cumulative values for individual hydrogen bonds is presented in Table 5 and Table S3.

For all analyzed compounds, several general trends can be observed with respect to both the number and quality of the formed hydrogen bonds. Among all interactions identified during the docking phase, only two hydrogen bonds retained full stability and persistence in the conformers sampled during molecular dynamics simulations: Ligand (H1)···(O) GLU81 and Ligand (O1)···(HN) LEU83. These interactions were present in 100% of the conformers obtained during the molecular dynamics simulations; however, significant changes in their bond lengths were observed compared to those identified at the docking stage, indicating a reorientation of the isatin core relative to the hinge region of the protein comprising LEU83 and GLU81.

The observed discrepancies are manifested by a pronounced shortening of the hydrogen bond formed with GLU81, which initially oscillated around 2.6 Å in the docked complexes, whereas in the vast majority of conformers (≥80%) its length ranged from 1.6 to 2.0 Å. The second key hydrogen bond involving LEU83 also exhibited slight changes in its properties. The originally identified distance of approximately 1.8 Å was observed, depending on the ligand, in only 26–40% of the analyzed conformers, while a considerably larger population of structures (~60%) exhibited distances in the range of 2.0–2.2 Å.

Much more pronounced changes were observed for the Ligand (H5)···(O) LEU83 hydrogen bond, whose significance nearly vanished in complexes formed by ligands 2a and 2c. For many derivatives, the presence of this interaction was detected in fewer than 50% of the conformers, and the recorded bond lengths indicated a dominant contribution of weak interactions or interactions merely fulfilling the boundary criteria for hydrogen bonding (compounds 1, 2f, 2i, 3c, 4c, 4d, 4e, and 4f). Conversely, complexes were also identified in which this hydrogen bond occurred with higher frequency (70–80%), with a substantial fraction of the population corresponding to hydrogen bonds of intermediate strength (compounds 3a, 3b, 3d, 3e, 3h, 4a, and 4b).

The aforementioned reorientation of the isatin core simultaneously influenced the ability of fluorine atoms from the trifluoromethoxy group located at position 6 of the central ring to participate in hydrogen bond formation. The strong hydrogen bond involving the hydrogen atom of ASP145 exhibited more diverse characteristics in the systems sampled during molecular dynamics simulations. No complex was identified in which this interaction was present in 100% of the conformers; nevertheless, several complexes demonstrated its significant contribution to complex stabilization. Notably, derivatives 3h and 4e exhibited this interaction in at least 95% of the analyzed conformers. Although the bond lengths increased noticeably compared to the initial docked structures, the dominant fraction of the population could still be classified as hydrogen bonds of moderate strength. A significant stabilizing role of the hydrogen bond formed with ASP145 was also observed in complexes formed by derivatives 1, 2i, 3b, 3c, 3d, 3k, 4a, 4d, and 4f, where the fraction of conformers confirming this interaction ranged from 70 to 85%.

The most labile hydrogen bond identified during the docking stage is the interaction formed with LYS33. Analysis of the collected conformational populations indicates that its occurrence ranges from less than 20% to approximately 46%, with the dominant contribution arising from bond lengths characteristic of weak hydrogen bonds. For the majority of the analyzed populations, the fraction of conformers confirming the presence of this interaction oscillates around 35%. Although this hydrogen bond may contribute to the overall spectrum of interactions stabilizing the complexes, it cannot be regarded as a dominant factor governing their formation.

In the case of the reference molecule, it can be observed that the hydrogen bond formed between GLU81 and the hydrogen atom of the isatin core exhibits characteristics similar to those identified for the newly proposed compounds. In contrast, significant differences are evident in the properties of the hydrogen bonds formed with the hydrogen atom of LEU83. For this interaction, only a negligible contribution of strong hydrogen bonds is observed (bond length range 1.6–1.8 Å, approximately 6%), whereas hydrogen bonds of intermediate strength constitute the dominant group, accounting for approximately 71% of the conformers.

Another important factor determining the stability of the studied complexes is the broad range of hydrophobic interactions. A key aspect of this type of interaction is π–π stacking between aromatic systems, the strength of which is primarily determined by the distance between the interacting aromatic rings. Figure 4 presents a cumulative analysis illustrating the distribution of distances between the aromatic system of PHE80 and the isatin core (Figure 4a), as well as between PHE82 and the aromatic ring located in the side chain of the studied ligands (Figure 4b).

In the case of the first interaction, it can be observed that for the majority of the analyzed systems, the dominant fraction (70–80% of conformers) corresponds to distances oscillating around 4.0 Å. This value correlates well with the distances identified at the docking stage, and a similar interaction profile is also observed for the reference molecule. A subset of derivatives exhibiting a significantly lower tendency for deep insertion of the molecular core into the hydrophobic pocket can also be identified, namely derivatives 1, 2c, 3a, 4d, and 4f.

An exceptional case within the analyzed group is derivative 2a, for which nearly 85% of the conformers exhibit distances between the aromatic systems oscillating around 3.25 Å. Such pronounced penetration of the molecular core into the hydrophobic pocket is further supported by the high fraction of very strong hydrogen bonds formed with GLU81 (approximately 63%).

For the second set of analyzed aromatic interactions, a relatively uniform profile is observed for most of the studied systems, characterized by a dominant population of conformers in which the distance between the aromatic systems fluctuates around 4.0 Å. A slightly different behavior is observed for derivatives 2f and 3h, where approximately 40% of the conformers display significantly shorter distances, close to 3.25 Å. The only derivative for which a pronounced increase in the distance between the aromatic systems relative to the initial docked structure is observed is compound 3e, for which more than half of the analyzed population exhibits distances oscillating around 4.75 Å.

A noteworthy and expected exception is the reference structure (RM), which, due to the presence of symmetrically distributed chloro substituents, exhibits limited ability to interact within the active-site region containing the aromatic system of PHE82.

The final stage of the molecular dynamics–based analyses, aimed at extending the characterization of the studied systems beyond structural stability and conformational variability of the active sites, was based on the Molecular Mechanics Poisson–Boltzmann Surface Area (MMPBSA) method. The calculations were performed using systems extracted from the final and most stable phase of the molecular dynamics simulations, covering the last 60 ns of the trajectories. The collected data, presenting averaged enthalpic contributions to binding affinity along with standard deviations, are shown in Figure 5. It can be observed that although a considerable population of conformers exhibited higher binding affinity values than the native ligand (i.e., 2a, 3d, 3h, 3k, 4a, 4e, and 4f), statistical analysis using one-way analysis of variance (ANOVA) with the pairwise multiple comparison Holm–Šidák method did not indicate the presence of statistically significant differences between the compared datasets.

Within the studied group, derivatives that showed statistically significantly lower binding affinities relative to the native system included compounds 2c, 3a, and the reference molecule (RM). In contrast, a substantially broader group of derivatives demonstrated statistically significant improvements in binding affinity relative to the reference molecule, including compounds 1, 2a, 3d, 3h, 3k, 4a, 4e, and 4f.

Ligand–binding site interactions are inherently bidirectional; not only does the ligand undergo conformational adaptation to the active site, but the active site itself may also experience ligand-induced conformational changes. To assess this phenomenon, extended molecular dynamics simulations of 160 ns were performed for representative systems, including the apo protein as well as complexes with the native ligand and selected derivatives bearing substituents at the R2, R3, and R4 positions (2a, 3h, and 4e), which also exhibited the highest affinity toward the active site. Based on the collected trajectories, root mean square fluctuation (RMSF) values were calculated for individual amino acid residues comprising CDK2. A comparative summary of the RMSF profiles is presented in Figure 6a,b. The largest structural fluctuations of the analysed protein are observed in two residue ranges, namely 35–55 and 145–165. The former corresponds to a loop region and a fragment of an α-helix located in close proximity to the active site, whereas the latter represents the most flexible segment of the protein, namely the T-loop. More specific and less pronounced changes are observed for other residue ranges forming β-sheet elements that delimit the active-site cavity in the upper (residues 6–15, 20–30, and 70–85) and lower planes (residues 130–140).

The observed fluctuations clearly indicate that the presence of a bound ligand significantly influences the conformational properties of the active site and the structural dynamics of the amino acid residues forming this region.

2.4. ADMET Analysis

Assessment of the interaction potential of a candidate drug with its biological target constitutes a fundamental aspect of drug design; however, it is accompanied by other essential activities aimed at evaluating the overall molecular properties and in vivo behavior of the investigated substances. In this context, methods based on quantitative structure–activity relationship (QSAR) analysis are particularly useful. Such tools enable the evaluation of a broad spectrum of compound properties, including absorption, distribution, metabolism, excretion, and toxicity (ADMET) [26,27,28,29]. The selected set of descriptors characterizing the properties of the investigated group of compounds is presented in Table 6.

To facilitate a better interpretation of the ADMET parameters describing the studied group of compounds, a broader set of reference compounds was included in the analysis. This reference group comprises substances not only investigated in clinical studies but also widely used in the development of new and effective anticancer therapies, namely roscovitine (RM2) [30,31] and dinaciclib (RM3) [32], which are anticancer drugs acting as cyclin-dependent kinases (CDK’s) inhibitors.

One of the most basic parameters considered in this dataset is molecular weight (MW), which for the analysed population of compounds ranges from 363.08 to 433.05 u. All identified values fall below the commonly accepted threshold (MW < 500 u) [33]. The subsequent parameters describe the numbers of hydrogen bond donors (nHD) and hydrogen bond acceptors (nHA). In both cases, the values identified for the investigated compounds comply with the criteria defined for molecules with pharmacological potential [33].

Comparison of the aforementioned parameters for the reference molecule 1 (RM) and the investigated group of compounds reveals that the reference compound possesses a lower molecular weight and a lower degree of saturation with functional groups acting as hydrogen bond donors and acceptors. This trend is also reflected in another analyzed parameter, namely the topological polar surface area (TPSA) [34], for which the values recorded for the investigated compounds are significantly higher and, in several cases, nearly twice as large. The observed disparities in TPSA values may indicate differences in the ability of the studied compounds to permeate environments where lipophilicity is a critical factor. Nevertheless, this observation does not imply an inability of the investigated compounds in this regard, but rather suggests a relative advantage of the reference compound 1. When the parameter values obtained for the studied group of compounds are compared with those calculated for roscovitine (RM2) and dinaciclib (RM3), it can be observed that they are characterized by very similar numbers of hydrogen bond donors and acceptors, as well as comparable TPSA values, indicating a high degree of similarity in the overall physicochemical profiles of the compounds under comparison.

The subsequent parameters presented relate to the solubility of the investigated compounds, namely the octanol/water partition coefficient (logP) and water solubility (logD). The reported data indicate that, in most cases, the studied substances exhibit values that slightly exceed the optimal range recommended for pharmacologically active compounds; however, they remain well below the established threshold limits [33]. In the analysis of the solubility-related parameters of the studied compounds, it can be observed that although the investigated substances do not exhibit optimal solubility properties, in most cases they demonstrate more favorable characteristics than roscovitine (RM2). This suggests that the potentially lower solubility of the studied compounds should not significantly limit their applicability or therapeutic potential in anticancer treatments.

Issues associated with limited solubility in common solvents can be substantially mitigated through the use of binary solvent systems based on trace amounts of organic solvents, which are capable of significantly improving the solubility of the investigated substances, as demonstrated in studies conducted on analogous classes of isatin derivatives [25].

Based on the descriptors presented thus far, the properties of the investigated compounds can be evaluated using the most widely recognized drug-likeness criterion, namely Lipinski’s Rule of Five. This rule directly links the pharmacological potential of candidate compounds to compliance with the following guidelines: molecular weight (MW ≤ 500), logP (≤5), number of hydrogen bond acceptors (nHA ≤ 10), and number of hydrogen bond donors (nHD ≤ 5) [33]. The values presented unequivocally indicate that the proposed group of compounds satisfies all the requirements of this rule.

Another important set of parameters describes blood–brain barrier (BBB) permeability and clearly indicates that none of the proposed compounds is predicted to exhibit difficulties in crossing the barrier between the bloodstream and the central nervous system. An additional crucial aspect of drug absorption is its ability to be absorbed and permeate following oral administration. The reported values of the Caco-2 [35] and human intestinal absorption (HIA) parameters indicate that all proposed compounds exhibit favorable absorption properties upon oral delivery.

The hERG parameter represents a critical safety indicator used to assess a compound’s potential to inhibit hERG potassium (K^+^) channels [29], an effect that may contribute to the development of long QT syndrome (LQTS), cardiac arrhythmias, and Torsade de Pointes. The obtained results indicate that all analyzed derivatives exhibit a moderate risk of this type of inhibitory activity, comparable to that observed for the reference compounds RM2 and RM3, while the risk associated with RM remains negligible.

The predicted Ames toxicity values suggest that none of the proposed compounds is expected to display mutagenic potential. In contrast, the reference molecules show a markedly different profile, with a substantial probability of mutagenicity ranging from 39% (RM3) to 83% (RM2).

The final evaluated parameter, CARC, reflecting carcinogenic potential, shows considerable variability among the analyzed compounds. Certain derivatives, such as 3b, 4b, and 4e, exhibit negligible carcinogenic risk, whereas others, including 2f and 4f, demonstrate an intermediate carcinogenic potential. Nevertheless, for all investigated compounds, the predicted carcinogenic probability remains lower than that observed for the reference molecules RM and RM2.

3. Materials and Methods

3.1. Ab-Initio Calculations

The geometric and electronic properties of the studied compounds were evaluated using Density Functional Theory (DFT) [36,37]. The structures of the investigated inhibitors were prepared using GaussView 6.0.16, and all quantum-chemical calculations were performed with the Gaussian 16 Rev. C.01 package [38], while the analysis of the frontier orbitals was carried out using the Avogadro 1.2.0 application [39]. Geometry optimizations and molecular properties were calculated at the B3LYP/6-31G (d,p) level of theory [40,41]. The energies of the HOMO and LUMO orbitals were used as source data for the estimation of several chemical reactivity descriptors, providing an opportunity to characterize a molecule’s chemical nature and reactivity. The analyses were based on the values of the energy gap (ΔE_GAP_) (2), absolute hardness (η) (3), absolute electronegativity (χ) (4), chemical potential (μ) (5), absolute softness (σ) (6), global softness (S) (7), global electrophilicity index (ω) (8), and the maximum amount of additional electronic charge (ΔN_MAX_) (9) [42,43,44].

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

3.2. The Docking Procedure and In-Silico ADMET Analysis

The geometric structure of the enzyme cyclin-dependent kinase 2 (CDK2; PDB ID: 1E9H) was downloaded from the Protein Data Bank (PDB) [45]. During the docking procedure, two programs were used: initial preparation was performed using the AutoDockTools 1.5.6 package [46], while docking simulations were carried out using AutoDock Vina software [47]. During the initial preparatory steps, all non-polar hydrogen atoms were removed from all chemical structures, including the protein and the inhibitors; additionally, the number of rotatable bonds was defined for the considered ligands. The grid box dimensions of 16 × 14 × 22 Å were set, ensuring that the docking procedure fully covered the active-site space and enabled a reliable assessment of ligand affinity. The docking protocol was validated by re-docking the co-crystallized ligand into the CDK2 active site. The resulting geometry closely matched the experimental conformation, with the RMSD between the predicted pose and the original crystallographic binding mode well below 2.0 Å, confirming the reliability and accuracy of the selected docking parameters. For each inhibitor candidate, the docking procedure was repeated five times using different random seed values to ensure proper execution and to obtain representative and reproducible results. In all considered complexes, the predicted conformations were highly consistent and corresponded to the poses with the lowest predicted binding affinity, further supporting the validity of the docking approach. During the simulations, the exhaustiveness parameter was set to 20, which ensured an appropriate balance between precision and computational time.

The ADMET (Absorption, Distribution, Metabolism, Excretion, and Toxicity) properties of the considered molecules were evaluated using the ADMETlab 3.0 platform [48], which enables quantitative structure–activity relationship (QSAR) analysis based on two-dimensional molecular structures described by SMILES notation. Among the numerous available descriptors, those defining molecular properties (molecular weight, number of hydrogen bond donors (nHD), number of hydrogen bond acceptors (nHA), topological polar surface area (TPSA), logP, logD, and logS, absorption (Caco-2 permeability and human intestinal absorption (HIA)), distribution (blood–brain barrier penetration), and toxicity (hERG blockers, AMES toxicity, and carcinogenicity (CARC)) were selected to obtain a comparative characterization of the properties of the candidate molecules exhibiting the highest binding affinity to the CDK2 active site.

3.3. The Molecular Dynamics Simulations

During all molecular dynamics simulations, the protein component of the studied complexes (CDK2) was described using the ff14SB force field [49], while the ligand structures were parameterized using the General Amber Force Field (GAFF) [50]. Partial atomic charges for the inhibitors were derived using the Merz–Kollman scheme followed by the RESP fitting procedure at the HF/6-31G* level of theory [51]. Each molecular system considered during the molecular dynamics stage was neutralized with chloride ions and immersed in a periodic box of TIP3P water molecules. The initial phase consisted of energy minimization and gradual heating of the system to 300 K; temperature regulation was carried out using a Langevin thermostat [52]. After these preparatory steps, production molecular dynamics simulations were performed for 80 ns under periodic boundary conditions, with the SHAKE algorithm applied to constrain bonds involving hydrogen atoms.

The analysis of the studied complexes, including the estimation of structural stability, flexibility, and mobility of individual components of the tested systems, was carried out using the VMD package [53]. During the structural analysis of interactions within the active-site region, hydrogen bonds were defined according to the following criteria: a hydrogen–acceptor distance ≤ 3.0 Å, a donor–acceptor distance ≤ 3.5 Å, and a donor–hydrogen–acceptor angle greater than 90°. For selected systems, namely the apo CDK2 protein and its complexes with compounds 1, 2a, 3h, and 4e, extended molecular dynamics (MD) simulations of 160 ns were performed to investigate the potential effects of ligand binding on the active site. The simulations were analyzed in terms of root mean square fluctuations (RMSF) to characterize the flexibility and conformational dynamics of individual amino acid residues within the active site and its immediate surroundings. The molecular mechanics Poisson–Boltzmann surface area (MMPBSA) method [54] was employed to estimate the enthalpic contributions to binding affinity values characterizing the complexes formed by ligands with the CDK2 active site. The calculations were performed based on 100 randomly selected conformers of the analyzed systems from the last 60 ns of the molecular dynamics simulations. The molecular dynamics simulations were performed using the AMBER 14 package [55]. One-way analysis of variance (ANOVA) with the pairwise multiple-comparison Holm–Šidák method [56,57] was employed in the statistical analysis of mean values and standard deviations for RMSD and enthalpic contributions to binding affinity. All statistical tests were conducted with an overall significance level of α = 0.05.

4. Conclusions

The conducted studies, including quantum mechanical calculations, molecular docking, molecular dynamics simulations, and QSAR-based ADMET analysis, enabled a comprehensive evaluation of the investigated compounds as potential cyclin-dependent kinase 2 (CDK2) inhibitors. Structural diversity arising from electron-donating and electron-withdrawing substituents on the aromatic system allowed assessment of binding capabilities and identification of derivatives exhibiting the highest conformational stability within the enzyme active site. Docking results revealed relatively low variability in affinity across the compound set and emphasized the dominant role of the isatin core bearing a novel substituent that reorients the molecular scaffold and actively participates in intermolecular interactions.

The trifluoromethoxy substituent at the C6 position was identified as a key determinant of active-site orientation through hydrogen-bond formation with ASP145 and LYS33, producing a binding geometry distinct from previously studied C5-substituted isatin derivatives and altering interaction distances with GLU81 and LYS83. Hydrogen-bond and distance distribution analyses indicated deeper embedding of the hinge region within the hydrophobic pocket, particularly for representative derivatives 2a, 3h, and 4a. Molecular dynamics simulations further confirmed favorable accommodation and conformational stability of selected compounds—1, 2a, 3e, 3h, 3k, 4f, and 4i—as evidenced by stable RMSD profiles and consistently low enthalpy values.

ADMET evaluation demonstrated compliance with fundamental drug-likeness criteria, with only moderate solubility limitations that remain below pharmacological exclusion thresholds and may be mitigated using binary solubility enhancers, as reported for analogous isatin derivatives. Additional descriptors indicated favorable absorption, distribution, and low toxicity or mutagenicity risk.

Overall, the strong affinity toward CDK2 and favorable physicochemical and pharmacokinetic properties—especially for the representative derivatives highlighted above—identify the proposed compounds as promising candidates for further investigation, including chemical synthesis and cellular evaluation.