Computational discovery of novel anticancer agents: A comprehensive study of phenothiazine derivatives

2.1 Computational analysis

The chemical structures of the phenothiazine derivatives (1-14) were designed using ChemDraw, while the reference drug structures were retrieved in SDF format from the PubChem database (Figs. 1 and S1 in the Supplementary Information). Geometry optimization of all ligands and reference compounds was performed with Gaussian16 (Frisch et al., 2016) at the B3LYP/6-31G(d,p) level of theory, which has been extensively validated for similar molecular systems (Uddin et al., 2025; Uddin et al., 2023; Uddin et al., 2022). Frequency calculations were conducted to confirm the absence of imaginary frequencies, ensuring that all optimized structures represented true energy minima.

Molecular electrostatic potential (MEP) maps and frontier molecular orbital (FMO) distributions, including the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO), were generated using GaussView 6 (https://gaussian.com/gaussview6/). Representative HOMO-LUMO plots, MEP maps, and optimized structures are provided in the Supplementary Information (Tables S1-S18, Figs. S2–S7). The FMO energy gap (E_gap), an indicator of molecular stability, was calculated as the difference between EHOMO and ELUMO. These orbital energies were further used to derive a series of global chemical reactivity descriptors: ionization potential (IP), electron affinity (EA), electronegativity (χ), chemical potential (µ), global hardness (η), global softness (σ), electrophilicity index (ω), maximum charge acceptance (ΔNmax), and energy change (ΔE) (Uddin et al., 2025; Uddin et al., 2023; Uddin et al., 2022). The descriptors were computed using the following standard relations:

$\begin{array}{c}\text{IP}=-{\text{E}}_{\text{HOMO}};\text{EA}=-{\text{E}}_{\text{LUMO}};\text{µ}=(\text{IP}+\text{EA})/2;\\ \text{H}=(\text{IP}-\text{EA});\chi =-\eta ;\omega =\hat{\mu }2/\left(2\eta \right);\sigma =1/\eta ;\end{array}$

$\Delta {\text{N}}_{\text{max}}=–(\mu /\eta );\Delta \text{E }=–\omega ;{\text{E}}_{\text{gap}}={\text{E}}_{\text{LUMO}}–{\text{E}}_{\text{HOMO}}$

This computational workflow provided insights into the stability, electronic properties, and reactivity of the synthesized phenothiazine derivatives, forming the basis for subsequent docking and MD simulations.

2.2 Evaluation of physicochemical and pharmacokinetic properties

We evaluated the drug-likeness, medicinal chemistry, lipophilicity, physicochemical properties, pharmacokinetics, and aqueous solubility of the compounds using publicly available web-based platforms, including AdmetSAR (http://lmmd.ecust.edu.cn/admetsar2/) and SwissADME (www.swissadme.ch) (Daina et al., 2017; Yang et al., 2019). For this purpose, the Simplified Molecular Input Line Entry System (SMILES) format was employed. The chemical structures of compounds 1-14 were initially drawn using ChemBioDraw Ultra 14.0 and subsequently converted into SMILES format (Table S19, SI) to generate data in MDL Molfile format.

2.3 Pharmacological attributes

The PASS web server (http://www.way2drug.com/PASSOnline/predict.php) was used to predict the pharmacological activities of compounds 1-14 and the reference drugs. The analysis suggested multiple potential effects, including anaphylatoxin receptor antagonism, antineoplastic activity, NADPH oxidase inhibition, and GPI-PLD inhibition (Druzhilovskiy et al., 2017). The PASS database incorporates information on biologically active substances approved in Russia and the United States (Zardecki et al., 2016).

2.4 Molecular docking

2.5 Molecular dynamics

MD simulations were performed using GROMACS 2021.6 (https://manual.gromacs.org/2021.6/index.html) to examine interactions between protein-compound complexes (3, 9 with 4KQ8 and 9, 12 with 1CKP) as well as the unbound 4KQ8 structure. The docked protein-ligand complexes were first separated into individual protein and ligand files. The protein structure was prepared by removing hydrogens and generating the topology using GROMACS tools. The ligand structure was converted to PDB format, assigned a protonation state corresponding to pH 7.4, and processed with ACPYPE to generate GROMACS-compatible topologies. Subsequently, the protein and ligand topologies were combined, and the resulting complex was placed in a triclinic simulation box with a 1 nm buffer distance from the edges. The system was then solvated, and counterions were added to neutralize the overall charge. The final set of files included the solvated structure (GRO), topology (TOP), and positional restraints (ITP), which were used in subsequent energy minimization and MD simulations. Atomic charges were derived using the AM1-BCC method with Antechamber, and GROMACS topologies were generated using the ACPYPE model. The AMBER99SB force field was applied for protein topology generation, and the Gaff force field was used for ligand parameterization. Topologies were generated via the Galaxy servers (https://usegalaxy.eu/login/start and https://usegalaxy.org.au/). The systems were solvated in a triclinic box with the TIP3P model, neutralized with Na⁺/Cl⁻ ions, and equilibrated under NVT conditions at 300 K for 3000 ps using the leapfrog algorithm (Uddin et al., 2025; Uddin et al., 2023; Uddin et al., 2022). Production runs were conducted for 100 ns at constant temperature and pressure.

Trajectory analyses were performed with GROMACS utilities (gmx rmsd, gmx gyrate, gmx rmsf, gmx hbond) to monitor structural stability, compactness, and hydrogen bonding. Principal component analysis (PCA) was carried out using the Bio3D package (https://bio.tools/bio3d) to assess conformational dynamics of compound-protein complexes at 300 K, with additional PCA performed for compound 10 at 300, 305, 310, and 320 K. Cosine content analysis was included to evaluate the convergence of motion. These simulations provided key insights into the stability and dynamic behavior of the investigated protein-ligand systems.

3.1 Validation of protein structure

The structural integrity of CDK2 and aromatase was assessed using multiple validation tools. Ramachandran plot analysis via PROCHECK (https://www.ebi.ac.uk/thornton-srv/software/PROCHECK/), indicated that 88.8% of CDK2 residues and 84.2% of aromatase residues were in favored regions (Fig. 2a). ERRAT quality (https://www.doe-mbi.ucla.edu/errat/) scores were 92.51% and 92.08%, respectively (Fig. 2b), while Verify3D confirmed that 82.80% of CDK2 and 74.56% of aromatase residues achieved acceptable 3D–1D compatibility (Fig. 2c). ProSA-web analysis yielded Z-scores of –6.05 (CDK2) and –8.96 (aromatase), consistent with experimentally determined structures (Fig 2d). Local energy profiles further demonstrated stable residue distribution across both proteins (Fig 2e). Ligand structures were energy-minimized using OpenBabel in PyRx (https://pyrx.sourceforge.io/). Protein–protein interaction (PPI) networks were constructed using STRING (http://string-db.org/), considering interactions with a confidence score >0.4. CDK2 exhibited the highest confidence score (0.999), while aromatase also showed strong connectivity (0.986) (Figs. 2f and 2g), confirming their biological significance within cellular pathways.

Close

Fig. 2.

Structural validation of CDK2 and aromatase using (a) Ramachandran plots, (b) ERRAT analysis, (c) Verify3D profiles, (d) ProSA-web Z-scores, (e) ProSA-web local quality plots, (f) STRINGS output and (g) Interaction partner: (I) CDK2 and (II) aromatase.

Export to PPT

3.2 Frontier molecular orbital analysis

FMO analysis was performed to evaluate the stability and reactivity of compounds 1-14. The HOMO-LUMO energy gap (E_gap) is a key descriptor of chemical reactivity, softness, hardness, electrophilicity, and chemical potential (Cohen & Benson, 1993). Smaller E_gap values indicate higher reactivity and lower stability, while larger gaps correspond to greater stability and reduced reactivity. Among the tested derivatives, compound 9 showed the highest E_gap (4.096 eV), suggesting enhanced stability, whereas compound 10 exhibited the lowest E_gap (0.393 eV), indicating high reactivity and low stability (Table 1). Overall, E_gap values ranged from 0.393 to 4.096 eV, with the order: 9 (4.096) > 1 (3.806) > 13 (3.672) > 14 (3.671) > 2 (3.654) > 6 (3.036) > 12 (3.021) > 5 (3.008) > 4 (2.823) > 11 (2.788) > 3 (2.731) > 8 (2.716) > 7 (2.100) > 10 (0.393). For comparison, reference drugs ranked as follows: letrozole (5.335 eV) > tamoxifen (4.934 eV) > flavopiridol (4.146 eV) > doxorubicin (3.568 eV).

The anticancer compounds generally displayed higher reactivity than reference drugs, supported by lower hardness, increased electronegativity, electrophilicity, and dipole moments. Notably, compounds 3 and 12 showed favorable electrophilicity indices and strong dipole moments, suggesting enhanced binding potential (Figs. S2‒S7 in the SI). MEP maps revealed electron-rich (red) and electron-deficient (blue) regions, highlighting probable nucleophilic and electrophilic sites, respectively. These findings identify compounds 3 and 12 as promising candidates for further docking and biological evaluation.

3.3 In silico molecular docking

Computational docking simulations were performed to evaluate the interactions of the natural anticancer compounds (1-14) and reference drugs with six target proteins: CDK2 (PDB: 1CKP), SKP1/SKP2 (PDB: 1FQV), transthyretin (PDB: 4D7B), aromatase (PDB: 4KQ8), adenylyl cyclase (PDB: 5D0R), and GRP78 (PDB: 5F1X). Docking results (Tables 2-3, S20‒S25; Figs. 3, S8–S11) showed that most compounds exhibited stronger binding affinities than reference drugs. Compound 12 displayed the highest affinity toward CDK2 (−10.9 kcal∙mol⁻¹) and GRP78 (−8.9 kcal∙mol⁻¹), while compound 9 bound strongly to SKP1/SKP2 (−10.0 kcal∙mol⁻¹). Compound 13 showed dual targeting with TTR and AC (both −7.3 kcal∙mol⁻¹), and compound 3 exhibited high affinity for aromatase (−10.0 kcal∙mol⁻¹).

Table 3. Ligand-protein interacting amino acid residues of ligand 3 against aromatase (PDB: 4KQ8) and ligand 12 against CDK2 (PDB: 1CKP).

Complex	Interacting residues	Distance (Å)	Type of interaction	2D diagram of interaction
3-Aromatase	A:LYS396 A:TYR397 A:TYR317 A:TYR380 A:PRO385 A:PRO385 A:TYR317 A:TYR380 A:TYR380 A:PHE383 A:PHE386 A:PHE388	3.06 2.55 4.00 5.07 4.78 4.13 5.07 4.64 3.84 4.43 4.88 4.52	Conven. HB Conven. HB Pi-Sigma Pi-Sulfur Alkyl Alkyl Pi-Alkyl Pi-Alkyl Pi-Alkyl Pi-Alkyl Pi-Alkyl Pi-Alkyl
12-CDK2	A:ASP86 A:VAL18 A:LEU134 A:VAL18 A:LEU134 A:VAL18 A:ALA31 A:LYS33 A:LYS88	2.09 3.99 4.32 4.12 5.46 5.40 4.83 5.27 5.09	Conven.HB Pi-Sigma Alkyl Pi-Alkyl Pi-Alkyl Pi-Alkyl Pi-Alkyl Pi-Alkyl Pi-Alkyl

Close

Fig. 3.

Molecular docking results: (a) compound orientation within the binding pocket; (b) visualization of the active site; (c) hydrogen bond interactions; and (d) 2D protein–ligand interaction diagrams for ligand 12 with CDK2 (1CKP) and (e) compound orientation within the binding pocket; (f) visualization of the active site; (g) hydrogen bond interactions; and (h) 2D protein–ligand interaction diagrams for ligand 3 with aromatase (4KQ8)..

Export to PPT

Detailed interaction analyses (Tables 3; S20–S25) revealed multiple hydrogen bonds, pi–pi stacking, pi–alkyl, and other noncovalent contacts stabilizing the compounds in the protein active sites. For example, compound 12 formed 14 distinct interactions in CDK2, including hydrogen bonds with ASP145 and multiple pi–alkyl contacts with VAL164, LEU134, and ALA31. Compound 3 established 10 interactions in aromatase, including hydrogen bonding and pi–pi stacking with TYR424. Comparative docking with tamoxifen highlighted additional stabilizing contacts in both CDK2 and aromatase.

Bond distribution analysis showed that CDK2 interactions were dominated by pi–alkyl (35%), hydrogen bonds (11%), and pi–sigma (15%), with key residues ALA31, ILE10, LEU134, and VAL18 frequently involved. Aromatase interactions were enriched in hydrogen bonds (18%), pi–alkyl (26%), and pi–stacking contacts (13%), with TYR361, TYR424, LYS440, and ARG403 identified as critical residues. Overall, compounds 3 and 12 emerged as lead candidates, showing higher binding affinities than reference drugs, favorable pharmacokinetic profiles, and stable contacts during MD simulations. These findings support their potential as inhibitors of CDK2 and aromatase for anticancer therapy.

Beyond basic docking, this study examined the interactions of 14 candidate anticancer compounds with CDK2 and aromatase, alongside four reference drugs (Figs. 4 and S12–S13 in the SI). Binding analyses revealed multiple interaction types, including hydrogen bonds, π–alkyl, and π–π interactions, which stabilize ligands within protein binding pockets. For CDK2, the dominant interactions included 20 conventional hydrogen bonds (11%), 32 alkyl (17%), 28 π–σ (15%), and 67 π–alkyl (35%), with ALA31, ILE10, LEU134, and VAL18 as key residues. Aromatase interactions were characterized by 25 hydrogen bonds (18%), 37 π–alkyl (26%), 19 π–π stacking (13%), and 18 alkyl (13%), with notable contributions from TYR361, TYR424, LYS440, and ARG403. Compounds 3 and 12 consistently showed strong binding to aromatase and CDK2, respectively, outperforming reference drugs. Their favorable physicochemical and pharmacokinetic profiles, together with stable contacts during MD simulations, underscore their potential as lead inhibitors.

Close

Fig. 4.

(a) Interacting residues, (b) distribution of noncovalent interactions, and (c) interaction maps of compounds 3, 9, 12, and the reference drug tamoxifen with (I) CDK2 and (II) aromatase.

Export to PPT

3.4 Analysis of physicochemical and pharmacokinetic properties

The physicochemical profiles of the anticancer compounds were evaluated against Lipinski’s “rule of five” and Veber’s rule to assess oral bioavailability (Lipinski et al., 2001; Veber et al., 2002). According to Lipinski, suitable oral drugs should have MW ≤ 500 g/mol, log P ≤ 5, ≤ 5 H-bond donors, ≤ 10 H-bond acceptors, and TPSA ≤ 140 Ų, while Veber’s rule requires TPSA ≤ 140 Ų and < 10 rotatable bonds. SwissADME analysis (Table S26) showed that most compounds met these criteria, with bioavailability scores of 1, supporting their potential as orally active drugs. All compounds and reference drugs, except doxorubicin, complied with Lipinski’s criteria, though none satisfied the mLogP threshold. Despite this, all compounds and reference drugs met Veber’s criteria, indicating favorable molecular flexibility. Drug-likeness scores, calculated using Molinspiration (Table S26), revealed that compounds 9, 11, and 12 exhibited positive ion channel modulator (ICM) values, while several (9-12) also showed positive GPCR, nuclear receptor ligand, enzyme inhibitor, and kinase inhibitor activities. In particular, compound 11 recorded the highest NRL value (0.20) and EI value (0.33). In contrast, all 14 compounds displayed negative protease inhibitor values, whereas most reference drugs showed positive PI scores.

Pharmacokinetic properties were further assessed using SwissADME and admetSAR (Tables 5 and S27 in the SI) and include human intestinal absorption (HIA), blood-brain barrier (BBB) hERG inhibition, cytochrome P450 enzyme inhibition (CYP3A4 and CYP2C19), synthetic accessibility (SA) score, and plasma protein binding (PPB) (Fukunishi & Nakamura, 2011). ADMET characteristics were assessed using established prediction measures. High HIA (>85%) and BBB penetration were marked as + and low as -, with probabilities in parentheses. PPB (0-1) values represent plasma protein binding strength. Inhibition of CYP3A4 and CYP2C19 is shown by +/- probability. Subcellular localization predicts likely sites (M = mitochondria, N = nucleus) using confidence scores. The synthetic accessibility (SA) scale runs from 1 (easy) to 10 (difficult), representing the ease of synthesis. These scales provide insight into the substances’ pharmacokinetic and drug-likeness characteristics. All compounds exhibited excellent intestinal absorption (HIA > 0.95), with compounds 3, 4, 7, 8, and 9 reaching complete absorption (HIA = 1.0). BBB penetration values ranged from 3.52 (compound 13) to 4.77 (compound 5), while reference drugs varied between 1.48 and 4.98. CYP2C19 inhibition was higher in the candidate compounds compared to references, though tamoxifen showed the strongest inhibition overall. CYP3A4 inhibition was also more pronounced in the test compounds, with compound 10 showing the highest inhibition (0.9876). Subcellular localization analysis indicated that most compounds were mitochondrial, with compounds 5 and 6 exhibiting lysosomal activity. Among reference drugs, doxorubicin and flavopiridol localized to the nucleus, while letrozole and tamoxifen were mitochondrial. Overall, these results suggest that compounds, particularly 9, 11, and 12, display favorable physicochemical and pharmacokinetic profiles. While deviations in lipophilicity may affect absorption and distribution, their compliance with Lipinski’s and Veber’s rules supports their potential as orally bioavailable anticancer agents.

3.5 Analysis of pharmacological activities

The pharmacological activities of 14 naturally derived anticancer compounds were predicted using the PASS (Prediction of Activity Spectra for Substances) platform, which applies Multilevel Neighborhoods of Atoms (MNA) descriptors to assess biological activity with high accuracy (Druzhilovskiy et al., 2017). PASS generates probability values for activity (Pa) and inactivity (Pi) within the range 0–1. A Pa > 0.7 indicates a strong likelihood of experimental confirmation, values between 0.5 and 0.7 suggest moderate potential (possibly reflecting novel activity), and values < 0.5 represent low probability.

All tested compounds showed strong activity as anaphylatoxin receptor antagonists (Pa = 0.761–0.891), with compound 12 displaying the highest activity. Among the reference drugs, flavopiridol exhibited strong antagonistic potential (Pa = 0.859), while letrozole was inactive. Predicted antineoplastic activity was observed across all compounds (Pa = 0.221–0.524), though with lower probability compared to doxorubicin (Pa = 0.960). NADPH oxidase inhibition was absent in compounds 1, 2, 4, 7, and 8, as well as in all reference drugs (Table S28, SI), but other compounds demonstrated modest to strong inhibition (Pa = 0.103–0.759). For GPI-PLD inhibition, compounds 1, 2, 4, 7, and 8 were inactive, whereas compound 10 showed the highest predicted activity (Pa = 0.815). Among reference drugs, tamoxifen (Pa = 0.786) was active, while doxorubicin and flavopiridol were inactive. These results suggest that compounds 3, 10, 12, and 14 exhibit the most promising pharmacological profiles and warrant further experimental validation.

3.6 In Silico molecular dynamics

The structural stability and binding interactions of compounds 9 and 12 with CDK2 (PDB ID: 1CKP) and compounds 3 and 9 with aromatase (PDB ID: 4KQ8) were examined through 100 ns MD simulations. Docking results had previously indicated strong binding, with compound 12 showing a binding affinity of –10.9 kcal∙mol⁻¹ to CDK2 and compound 3 showing –10.0 kcal∙mol⁻¹ to aromatase. To validate these findings, MD simulations were performed for both protein-ligand complexes and their corresponding apo forms. Stability and conformational changes were assessed using RMSD, RMSF, radius of gyration (Rg), hydrogen bonding, and PCA. Selected complexes were also simulated at 300, 305, 310, and 320 K to investigate temperature-dependent dynamics. MD simulations were performed at 300 K, 305 K, 310 K, and 320 K to evaluate the thermal stability and resilience of the complexes under both physiological and mildly elevated temperature conditions, simulating cellular stress or fever. This approach allows for the assessment of conformational robustness beyond standard physiological settings.

RMSD analysis demonstrated increased stability of ligands upon binding. For CDK2, the compound 9-protein complex exhibited the narrowest RMSD fluctuations (0.15-0.61 nm), while the compound 12-CDK2 complex ranged from 0.10-2.30 nm (Figs. 5 and S14−S20 in the SI). For aromatase, compound 3 stabilized within 0.05-0.59 nm, and compound 9 showed even lower deviation (0.10-0.37 nm). In contrast, the free ligands exhibited markedly higher RMSD values, typically 4.2-6.0 nm, underscoring the stabilizing influence of protein binding. Slightly higher RMSD values observed in some systems reflect flexible loop regions accommodating ligand motion, rather than indicating structural instability.

Close

Fig. 5.

RMSD profiles: Compounds 9 (black curve) and 12 (red curve) at 300 K and docked complex of the protein CDK2 (PDB: 1CKP) with compounds 9 (green curve) and 12 (blue curve) at 300 K; (b) Compounds 3 (black curve) and 9 (red curve) at 300 K and docked complex of the protein aromatase (PDB: 4KQ8) with compounds 3 (green curve) and 9 (blue curve) at 300 K; (c) Compound 12 at 300 (black curve), 305 (red curve), 310 (green curve), and 320 K (blue curve); and target protein CDK2 with compound 12 at 300 (yellow curve), 305 (brown curve), 310 (grey curve), and 320 K (purple curve); (d) Compound 3 at 300 (black curve), 305 (red curve), 310 (green curve), and 320 K (blue curve); and target protein aromatase with compound 3 at (yellow curve), 305 (brown curve), 310 (grey curve), and 320 K (purple curve)

Export to PPT

Temperature-dependent simulations further suggest stability across all systems. For compound 12, free RMSD values varied between 2.6-6.3 nm, whereas in complex with CDK2, deviations were restricted to 0.05-2.3 nm at 300 K and narrowed further (0.19-0.48 nm) at 320 K. Compound 3 displayed free RMSD fluctuations of 4.4-5.9 nm, while its aromatase complex remained within 0.04-1.1 nm across all temperatures. Minimal scattering in RMSD profiles indicates robust, long-lasting interactions with the proteins, particularly at 310 K, suggesting thermally stable binding. Overall, compounds 3, 9, and 12 formed highly stable complexes with CDK2 and aromatase, highlighting their potential as strong inhibitors and supporting further preclinical evaluation.

Root-Mean-Square Fluctuation (RMSF) analysis was performed to examine residue-level flexibility in the studied complexes (Figs. S21–S23, SI). RMSF values reflect atomic mobility, with higher values indicating greater flexibility and lower values corresponding to structural rigidity. For the free compounds (Fig. S21A), compound 9 showed lower fluctuations (0.025–0.25 nm), whereas compound 12 exhibited slightly broader flexibility (0.025–0.30 nm), particularly around residues 4860 and 4890. In CDK2 complexes (Fig. S21b), compound 9 displayed fluctuations between 0.075 and 0.65 nm, while compound 12 ranged from 0.10 to 0.80 nm, with notable differences at residues 3400 and 7500. These results suggest that both ligands form stable complexes, though compound 12 showed localized flexibility that may affect binding dynamics.

For aromatase interactions, compound 3 displayed RMSF values of 0.01–0.16 nm, while compound 9 ranged from 0.02–0.18 nm (Fig. S21C). In the aromatase–ligand complexes (Fig. S21D), both compounds showed flexibility between 0.05 and 0.60 nm, with the largest fluctuations observed in the 2000–4000 residue region, while the 1000–2000 region remained more rigid. Overall, RMSF analysis indicated stable binding of compounds 3, 9, and 12 to CDK2 and aromatase, with localized flexibility revealing potential differences in interaction dynamics that may be relevant for drug design.

The radius of gyration (Rg) is a key parameter in MD simulations used to evaluate molecular compactness and conformational stability. Rg profiles for compounds 9 and 12 with CDK2, and compounds 3 and 9 with aromatase, were monitored over 100 ns (Supporting Information Figs. S24–S32). For the free ligands, compound 9 exhibited Rg values between 0.431–0.481 nm, while compound 12 ranged from 0.48–0.535 nm (Fig. S24A). Upon binding to CDK2, the protein maintained stable conformations, with Rg fluctuations of 1.975–2.06 nm (compound 9) and 1.999–2.06 nm (compound 12). Similarly, free compounds 3 and 9 displayed Rg values of 0.397–0.436 nm and 0.379–0.449 nm, respectively, while their complexes with aromatase showed stable compactness, ranging from 2.25–2.35 nm (compound 3) and 2.25–2.31 nm (compound 9) (Fig. S24B). These results suggest increased molecular dispersion upon binding but overall stable conformations in the protein-ligand complexes.

Thermal stability was further assessed. Compound 12 exhibited consistent Rg values across 300–320 K (0.48–0.535 nm), while its CDK2 complex remained stable with minimal fluctuations (1.95–2.09 nm) (Fig. S24C). Compound 3 showed similar stability across temperatures (0.39–0.435 nm), with its aromatase complex maintaining compact conformations between 2.24–2.33 nm, except for a transient deviation at 310 K (Fig. S24D). Overall, both CDK2–compound 12 and aromatase–compound 3 complexes retained compact and stable structures across simulations, indicating robust protein–ligand interactions under physiological and elevated thermal conditions.

Hydrogen bonding is a key determinant of protein–ligand binding strength and specificity. As illustrated in Fig. 6 (Supporting Information Figs. S33–S35), the number of intermolecular hydrogen bonds varied between 0 and 2 throughout the 100 ns MD simulations for both CDK2 and aromatase in complexes with compounds 9 and 12 (CDK2) and compounds 3 and 9 (aromatase). These fluctuations reflect conformational adjustments, ligand flexibility, and protein-specific structural dynamics. CDK2 complexes with compounds 9 and 12 exhibited dynamic yet consistent hydrogen bond formation (Fig. 6a), whereas aromatase complexes with compounds 3 and 9 showed comparable trends, indicating transient but stable interactions (Fig. 6b). Temperature-dependent analyses further suggested the persistence of hydrogen bonding across 300–320 K for both CDK2–ligand (Fig. 6c) and aromatase–ligand complexes (Fig. 6d), highlighting their structural resilience under thermal variation. Supplementary Figs. S33–S35 provide additional insights into the temporal evolution of hydrogen bonding. Compounds 9 and 12 with CDK2 and compounds 3 and 9 with aromatase both showed 100% hydrogen bond occupancies. Collectively, these results underscore the dynamic nature of hydrogen bonding and its critical role in stabilizing protein–ligand complexes, with important implications for binding affinity, specificity, and potential therapeutic efficacy.

Close

Fig. 6.

Time-dependent hydrogen bond profiles (ps) for (a) CDK2 complexes with compounds 9 (black curve) and 12 (red curve); (b) aromatase complexes with compounds 3 (black curve) and 9 (red curve); (c) CDK2–compound 12 complex at 300 (black curve), 305 (red curve), 310 (green curve), and 320 K (blue curve); and (d) aromatase–compound 3 complex at 300 (black curve), 305 (red curve), 310 (green curve), and 320 K (blue curve)

Export to PPT

To better characterize ligand–protein interactions over the course of the simulations, MD-derived interaction occupancies and contact heatmaps were generated. Hydrogen bond and hydrophobic contact occupancies were computed over the 100 ns MD trajectories to quantify the persistence of key interactions. Hydrogen bond occupancies were calculated using a distance cutoff of 0.35 nm and an angle cutoff of 30°, and expressed as the percentage of frames in which the interaction was maintained.

In addition, contact heatmaps were generated to visualize the frequency and spatial distribution of contacts between ligand atoms and protein residues throughout the simulations. These heatmaps highlight key residues involved in stable interactions and provide a clear overview of interaction hotspots and their temporal stability. As shown in Fig. 6, panels (A–D) illustrate hydrogen bond and hydrophobic contact occupancies for selected protein–ligand complexes over the entire 100 ns trajectory. Consistent interaction patterns, represented by dense regions in the heatmaps and sustained occupancy profiles, confirm the stability of the ligand binding modes within the active site.

Monitoring temperature stability is essential in MD simulations to ensure physiologically relevant conditions. Throughout the 100 ns simulations, all protein–ligand systems maintained a stable temperature, confirming well-controlled simulations and reliable structural and energetic analyses. As shown in Figs. S36‒S40 in the SI, potential energy fluctuations reflected dynamic atom-to-atom interactions within the complexes. For example, in CDK2 complexes (Fig. S36a), compound 9 exhibited potential energy values between –5.078 × 10⁵ and –5.013 × 10⁵ kJ∙mol⁻¹, while compound 12 ranged from –5.079 × 10⁵ to –5.015 × 10⁵ kJ∙mol⁻¹. In aromatase complexes (Fig. S36b), compound 3 fluctuated between –6.935 × 10⁵ and –6.875 × 10⁵ kJ∙mol⁻¹, and compound 9 varied from –6.935 × 10⁵ to –6.865 × 10⁵ kJ∙mol⁻¹. These observations underscore the inherent dynamism of protein–ligand interactions.

Electrostatic contributions were further assessed through Coulomb (short-range) energy analysis. For compounds 9 and 12, energy values ranged from –6.235 × 10⁵ to –6.15 × 10⁴ kJ∙mol⁻¹, while compounds 3 and 9 showed a narrower and more stable distribution between –8.57 × 10⁵ and –8.46 × 10⁵ kJ∙mol⁻¹ (Figs. S36c–d, S37‒S38, SI). The predominance of negative values indicates stable electrostatic interactions, supporting strong and persistent binding within the protein pockets. While these data confirm robust stabilization, more detailed energetic evaluations, such as binding free energy calculations, will be required to fully dissect interaction specificity and strength.

The Lennard–Jones short-range (LJ-SR) energy, reflecting van der Waals interactions, remained stable across the 100 ns simulations, indicating a consistent balance between attractive and repulsive forces within the complexes. As shown in Supporting Information Figs. S39‒S40, compounds 9 and 12 displayed LJ-SR energy values between 56,500 and 62,000 kJ∙mol⁻¹, while compounds 3 and 9 exhibited slightly higher values ranging from 74,000 to 81,000 kJ∙mol⁻¹. The stability of these energy profiles highlights persistent steric interactions at the protein-ligand interface, ensuring proper ligand positioning and orientation within the binding pocket. These findings emphasize the combined role of electrostatic and van der Waals forces in maintaining complex stability during MD simulations.

3.7 Principal component analysis

PCA was performed on the MD trajectories of the CDK2-compound 12 and aromatase–compound 3 complexes at four temperatures (300, 305, 310, and 320 K), focusing on Cα atoms. Using the Bio3D tool, PCA captured collective motions and conformational changes, with the first three eigenvectors (PC1–PC3) accounting for the most significant dynamic variations (Fig. S41, SI). Sample distributions were visualized as color gradients from blue (initial) to red (final), reflecting temporal progression. To determine the dominant motions during the simulation, PCA of the protein and protein–ligand complexes was carried out. The contributions of PC1, PC2, and PC3 to the overall motion were examined in order to evaluate the stability of the system. PC1, which stands for the direction of maximum variance, captures the most significant collective motions in each system; it should be noted that the PC1 of a protein–ligand complex and the PC1 of the free protein are derived independently, so they may correspond to different dominant motions. The variations in the PCA profiles show how ligand binding changes the conformational dynamics and stability of the protein. As summarized in Table 5, the CDK2–compound 12 complex showed the greatest PC1 variance at 305 K (18.51%), whereas the aromatase–compound 3 complex displayed the highest PC1 variance at 320 K (16.51%). Across all conditions, PC1 consistently accounted for the largest proportion of variance, highlighting its role in describing major conformational transitions.

Convergence values (H) further assessed dynamic stability. For CDK2–compound 12, H values were 0.75, 0.69, 0.88, and 0.02 (Table 4) at 300, 305, 310, and 320 K, respectively; for aromatase–compound 3, values were 0.84, 0.09, 0.88, and 0.81 at the same temperatures. These values, within the convergence threshold (0 < H < 0.5) (Gubbi et al., 2006), indicate adequate sampling and structural stability. Notably, the strongest associations were observed at 310 K, suggesting favorable binding of compound 12 with CDK2 and compound 3 with aromatase under physiological conditions. Collectively, these results reinforce the potential of compound 12 as a CDK2 inhibitor and compound 3 as an aromatase inhibitor, warranting further in vivo validation.

The MTT assay was used to assess the viability of HeLa carcinoma cells in the presence of different amounts of chemicals (Venkatesan et al., 2022). The HeLa cells were treated with substances at dosages of 1 mg, 10 mg, and 1000 mg following a day of cell incubation. As Table 6 illustrates, cell viability declines as compound concentration rises. Every derivative exhibited strong inhibition. HeLa cells were used in a lactate dehydrogenase (LDH) release assay to assess the cytotoxicity of the produced substances (Venkatesan et al., 2022). By measuring the amount of LDH released into the supernatant once membrane reliability is lost, this assay measures cell lysis. The substances at concentrations of 1 mg, 10 mg, and 1000 mg increased the release of LDH in a dose-dependent way.