Iron (III) triflate catalyzed, one-pot three-component synthesis of novel α-aminophosphonate of 5-hydrazineyl-2-(phenylsulfonyl)pyridine as potent antiproliferative agents against MCF-7 and HCT-116: In vitro and in silico investigations on molecular docking and dynamics simulation

2.2.1. Synthesis of 5-hydrazineyl-2-(phenylsulfonyl)pyridine (1)

5-Chloro-2-(phenylsulfonyl)pyridine (9.09 mmol) reacted with hydrazine 99% (36.4 mmol) in ethanol (40 mL) at reflux temperature for 5h; the reaction completion was monitored using TLC (Scheme 1). The reaction was cooled, evaporated and triturated with diethyl ether and the crude was purified using column chromatography with ethyl acetate/petroleum ether (30:70) to give compound 1 in 53% yield (1.2g), which is used in the next reaction directly.

Close

Scheme 1.

Synthesis of 5-hydrazineyl-2-(phenylsulfonyl)pyridine (1).

Export to PPT

2.2.2. Synthesis of α-aminophosphonates (4a-f)

A solution of compound 1 and diethyl phosphite (3) in 1,2-Dichloroethane (DCE) was added to the corresponding benzaldehyde compounds (2a-f). The reaction mixture was refluxed using Iron triflate Fe(OTf)3 for 2 h to produce α-aminophosphonates (4a-f).

2.3.1. Ligands and protein preparation

The molecular structure of compounds 4b/4e was accurately sketched using ChemDraw Ultra software. The initial structural representation, provided in Simplified Molecular Input Line Entry System (SMILES) notation, was subsequently converted to the Protein Data Bank (PDB) format. This conversion was facilitated by a custom Python script employing the RDKit cheminformatics toolkit, an open-source library for molecular manipulation and analysis (https://www.rdkit.org/).

The structures of MD2 and CDK2 was acquired from the RCSB, a curated archive of experimentally determined three-dimensional structures of biological macromolecules (https://www.rcsb.org/ligand/DM2/https://www.rcsb.org/structure/8FP0). The retrieved CDK2 structure (PDB ID: 8FP0) was refined using UCSF Chimera, a molecular visualization and analysis system. This refinement process involved the removal of extraneous polypeptide chains and bound ligands to isolate the target protein. Subsequently, any unresolved regions within the CDK2 structure were modeled using the loop modeling functionality within MODELLER, a comparative protein structure modeling program [38]. Finally, to optimize the refined CDK2 structure, energy minimization was performed using the CHARMM force field, implemented within the OpenMM molecular dynamics simulation toolkit. This minimization protocol consisted of 100 steps.

2.3.2. Molecular docking

The interaction between 4b/4e and DM2 chemical structures and the CDK2 protein was gauged using a computational docking approach. Molecular docking was performed with GNINA, a docking program that utilizes a convolutional neural network (CNN) ensemble to improve prediction accuracy [39]. A blind docking protocol was employed, where the search space, defined by a grid box, encompassed the entire CDK2 protein structure. Ten ligand poses in PDBQT format were generated, and the pose exhibiting the highest binding affinity was selected for further analysis.

To characterize the CDK2-ligand interactions, PyMOL was used to generate complex files representing the top-ranked pose. The Protein-Ligand Interaction Profiler (PLIP) web server (https://plip-tool.biotec.tu-dresden.de) was then used to analyze these complex files and identify key binding residues within the protein-ligand complexes. Complementary 2D interaction diagrams were generated using PoseView (https://proteins.plus/) to visualize the specific interactions. Furthermore, functional domains of CDK2 were identified using the InterPro database (https://www.ebi.ac.uk/interpro/), providing context for the observed interactions based on known protein family and domain information.

2.3.3. Molecular dynamic simulations

Molecular dynamics (MD) simulations were conducted using Maestro version 12.0.012 (Schrödinger, LLC, New York, NY). The CDK2-ligand complex was prepared using the Protein Preparation Wizard, which included preprocessing steps, energy minimization, and removal of crystallographic water molecules. The system was solvated using the simple point charge (SPC) water model, and counterions were added to neutralize the charge. The MD simulations were performed for 100 ns under constant temperature conditions (300 K).

Post-simulation trajectory analysis was performed to evaluate the stability and dynamics of the complex. Root mean square deviation (RMSD) and root mean square fluctuation (RMSF) were calculated using the interaction diagram module. Additionally, the radius of gyration (rGyr), polar surface area (PSA), and solvent-accessible surface area (SASA) were determined. A trajectory-based 3D structure analysis was conducted over the entire simulation period using the trj2mae.py utility script within Maestro.

Binding free energies were calculated using the Prime module. The calculations were performed on every 100th trajectory frame, providing a decomposition of the binding free energy into van der Waals, electrostatic, polar solvation, non-polar solvation, and entropic contributions. This decomposition allowed for a detailed understanding of the energetic factors driving ligand binding.

Principal Component Analysis (PCA) was carried out using the bio3d package in R to analyze the dominant modes of motion within the simulated complexes. The trajectory files were first converted from .xtc to .dcd format. The analysis focused on the C-alpha atoms, and their coordinates were extracted from the trajectory. The pca.xyz function, employing the singular value decomposition (SVD) method, was applied to these coordinates. The first two principal components (PCs), capturing the largest variance in the data, were used to represent the essential dynamics and conformational changes, providing insights into the complexes’ stability and flexibility.

Dynamic Cross-Correlation Matrix (DCCM) analysis was employed to investigate the correlated motions of Cα atoms within the protein-ligand complex using the bio3d package in R. The analysis, performed on the 0-100 ns trajectory, calculated the cross-correlation coefficient (Mxy) between the displacement vectors of all Cα atom pairs. The resulting cross-correlation values, ranging from -1 (perfect anti-correlation) to +1 (perfect correlation), with 0 indicating no correlation, were visualized as a heatmap. This allowed for identifying regions exhibiting correlated or anti-correlated motions, providing insights into the internal dynamics of the complex.

3.2.1. In-vitro cytotoxic activity

The resulted α-aminophosphonate phenylsulfonyl derivatives 4a-f were evaluated as antiproliferative agents to study their effects on tumor cells using the 3-[4,5-dimethyl-2-thiazolyl)-2,5-diphenyl-2H-tetrazolium bromide (MTT) assay against two types of human cancer cell lines breast (MCF-7) and colon (HCT-116), which is based on the cleavage of the tetrazolium salt by mitochondrial dehydrogenases in viable cells [18,41]. The percentages of intact cells were calculated and compared to those of the control activities of these compounds against the two carcinoma cell lines in comparison with DM2 as a standard drug. While DM2 primarily exerts its anticancer effects through topoisomerase inhibition, it was selected as a reference due to its well-established efficacy in comparison to the newly designed hybrid agents. Although the chemical structure of DM2 differs from our hybrid compounds, its broad-spectrum anticancer activity serves as a useful benchmark for evaluating the antiproliferative potential of our compounds [42,43].

It is apparent from Figure 3 and Table 2 that all the investigated compounds showed no significant cytotoxic activity against MCF-7 human breast cancer cells when compared to DM2.

Close

Figure 3.

Dose-dependent cytotoxic activities of compounds (4a-4f) against MCF-7 cancer cells according to the MTT assay.

Export to PPT

It is interesting to note that the compounds (4e) and (4b) were the most potent cytotoxic compounds against HCT-116 cell lines comparable to DM2 with IC₅₀ 9.2 and 9.5 µg/ml, respectively, as illustrated in Figure 4 and Table 2.

Close

Figure 4.

Dose-dependent cytotoxic activities of compounds (4a-4f) against HCT-116 cancer cells according to the MTT assay.

Export to PPT

Structure-Activity Relationship (SAR) for In-Vitro Cytotoxic Activity:

The SAR analysis of the α-aminophosphonate phenylsulfonyl derivatives (4a-f) provides insights into the structural features responsible for the observed cytotoxic activity against human cancer cell lines, particularly the breast cancer (MCF-7) and colon cancer (HCT-116) cell lines.

The compounds 4a-f displayed varying cytotoxicity, but none were as effective against MCF-7 breast cancer cells when compared to DM2. This suggests that the presence of para substituted phenyl group in the hybrid compounds may not significantly influence activity against this particular cell line.

Compounds 4e and 4b emerged as the most potent cytotoxic agents against the HCT-116 colon cancer cell line, showing IC₅₀ values of 9.2 µg/ml and 9.5 µg/ml, respectively. These values were comparable to the cytotoxicity of DM2.

The differences in cytotoxicity may be contributed to the presence of one of the following:

The presence of a phenylsulfonyl group in the structure of the compounds could enhance solubility and facilitate interactions with the cell membrane, improving cellular uptake [44,45].

The differences in cytotoxicity between compounds could be attributed to the substitutions on the aromatic rings. For example, Br in compound 4e or methyl in 4b could influence interactions with cellular targets or modify the compound’s stability.

The ability of these compounds to engage in hydrophobic interactions within the cancer cell environment may also be critical for their efficacy, particularly in the context of the breast cancer cells, where membrane penetration and cellular interactions are important for the drug’s effectiveness [46].

3.2.2.1. Molecular docking

The prospective candidates 4e, 4b and the reference drug DM2 were evaluated through in silico docking studies to gain further insight into the binding patterns of the target ligands with cyclin-dependent kinase 2 (CDK2), whose dysregulation is linked to various cancers, making it a significant target for therapeutic interventions [47,48]. CDK2 regulates the cell cycle, particularly the transition from the G1 to the S phase [49]. CDK2’s involvement in cancer is multifaceted, influencing cell proliferation, drug resistance, and the development of targeted therapies [50,51]. Nemours research reported that alterations in CDK2 signaling pathways are associated with cancers such as colon, breast, ovarian, and prostate cancer [52-54].

The experimentally obtained CDK2 structure was docked against 4b /4e and the reference drug DM2 using the GNINA, yielding binding affinity scores of -8.6 kcal/mol, -8.4 kcal/mol and -9.3 kcal/mol. Based on the molecular docking of DM2 with CDK2, we observed that the chosen drug binds to the same region as 4e and 4b within minimal differences in binding residues as depicted in Figures 5 and 6.

Close

Figure 5.

The shared binding pocket of the protein CDK2 (PDB ID: 8FP0) is visualized with the docked poses of the three ligands: 4e (yellow), 4b (green), and DM2 (pink).

Export to PPT

Close

Figure 6.

Interaction analysis of 4b/4e and DM2 with CDK2. (a) Visualization of CDK2 domains and structure. The 2D interaction diagrams of (b) 4e, (c) 4b and (d) DM2 with CDK2 residues highlights hydrophobic interactions (green) and hydrogen bonds (purple arrows).

Export to PPT

Upon docking the two compounds, 4b and 4e, it was revealed that both compounds engage in hydrophobic interactions and hydrogen bonding within the protein structure. 4e interacts hydrophobically with ILE (10), VAL (18), ALA (31), PHE (82), LYS (88), LYS (89), GLN (131), LEU (134), and GLU (162). Additionally, it forms hydrogen bonds with GLU (12), GLN (131), and ASN (132). In contrast, 4b interacts hydrophobically with ILE (10), VAL (18), ALA (31), GLN (131), and LEU (134), thereby showing less hydrophobic interactions, while forming hydrogen bonds with GLU (12), GLN (131), and ASN (132), similar to 4e.

Critically, both compounds interact with GLN (131), a residue within the active site (residues 123-135). Additionally, ILE (10), GLU (12) and VAL (18), located in the ATP-binding site, participate in hydrophobic interactions with both compounds. However, 4e establishes additional hydrophobic contacts with PHE (82), LYS (88), LYS (89), and GLU (162), interactions absent in 4b binding. Furthermore, both compounds form hydrogen bonds with GLU (12) in the ATP-binding region and ASN (132) near the C-terminal of the active site.

Both 4b and 4e interact with CDK2’s ATP-binding and active sites, implying they could affect CDK2 function. However, 4e’s more extensive hydrophobic interactions suggest a more substantial binding efficacy than 4b. The 2D interactions between CDK2- 4e, CDK2- 4b and CDK2- DM2 are shown in Figures 6(b-d).

3.2.3. Molecular dynamics simulation

The conformational stability of the CDK2- 4e complex was assessed via molecular dynamics simulations. The complex’s root-mean-square deviation (RMSD) remained within a narrow range of 2.0 Å to 3.5 Å throughout the simulation period, indicating a generally stable conformation. No abrupt conformational changes were observed. Following an initial equilibration period of approximately 4 ns, during which the RMSD stabilized around 2.5 Å, the complex maintained high stability, exhibiting only minor fluctuations. A gradual decrease in RMSD from 3.5 Å to 3.1 Å was noted between 60 ns and the end of the simulation. In comparison, the CDK2- 4b complex exhibited similar conformational RMSD values to CDK2- 4e until approximately 40 ns. However, after this point, the RMSD of CDK2- 4b gradually increased, diverging significantly from that of CDK2- 4e. By the end of the simulation, the RMSD of CDK2- 4b was approximately 1 Å higher than that of CDK2- 4e. This difference suggests that, while both complexes show initial stability, CDK2- 4b exhibits signs of minor instability and deviates from its initial conformation over approximately half of the simulation period. RMSD visualizations are shown in Figure 7(a).

Close

Figure 7.

Conformational and interaction analysis of CDK2 bound to 4b and 4e. (a) RMSD trajectories over time for CDK2-4b and CDK2-4e complexes. (b) RMSF analysis showing fluctuations at key residues. (c) SASA of the ligands throughout the simulation. (d) Ligand RMSD (e-f) Ligand stability as observed from the Rg and PSA plots. (g) Time-series analysis of protein-ligand contact frequencies with key interacting residues identified. (h) Residual interaction plots of 4b and 4e with CDK2.

Export to PPT

The RMSF analysis revealed that both complexes have similar fluctuations at key residue regions of the protein, such as around the 45th amino acid and the C-terminal. However, CDK2- 4b was observed to have higher fluctuations on the N-terminal when compared to CDK2- 4e. Similarly, several minor differential fluctuations were observed between 150-250 amino acid regions, Figure 7(b).

In a contrasting observation, it was revealed that both ligands show a different accessibility to the solvent through SASA analysis, especially after the 25th ns CDK2- 4b, whereby gradually increased from 350 Å to ∼800.0 Å whereas within the same timeframe, CDK2- 4e gradually decreased till 80th ns to 250 Å before going to around 400 Å at the end of the simulation, as shown in Figure 7(c,d). Interestingly, both ligands showed a similar radius of gyration (∼5.0 Å) and polar surface area (∼145-155.0 Å) were nearly identical, Figure 7(e,f).

Interestingly, the ligand 4e demonstrated highly significant stability while showing a sudden peak at the end of the simulation period, around 90-100 ns. The observed high energy elevation between 90-100 ns appears to stem from a sudden alteration in the protein residues interacting with the ligand. This involves a transition to using VAL-18, HIS-84, and LEU-83 more prominently, while interactions with GLY-11 and LYS-9, consistent throughout the simulation, become less significant. Interestingly, 3D trajectory analysis indicates the ligand remains securely in its binding pocket despite this energetic fluctuation. Despite this peak, the ligand remained relatively stable overall. In contrast, 4b exhibited marked instability, reaching peaks of 600.0 Å before settling back to ∼35 Å by the end of the simulation. Even then, the ligand continued to fluctuate substantially throughout the simulation, suggesting a high instability of the ligand when bound to CDK2, indicating a potential imbalance in interactions, as illustrated in Figure 7(d).

Protein-ligand contact time-series analysis revealed consistent interactions between the ligands and several key residues. The most prominent interacting residues for 4b throughout the simulation included 15-TYR, 33-LYS, 132-ASN, 131-GLN, 165-THR, 14-THR, 6-LYS, 9-LYS, 10-ILE and 16-GLY, which are comparatively more when compared to 4e which included 12-GLU, 20-LYS, 10-ILE, 131-GLN and 33-LYS. Interestingly, 12-GLU showed the highest interaction factor of 56%, followed by 33-LYS, both of which are in the protein’s active site. In CDK2- 4b, residue 15-TYR had the most prominent interaction at 27.83%, followed by 33-LYS (16.87%), 132-ASN(13.55%), 131-GLN (9.98%), 16-GLY (9.85%), 10-ILE (8.99%), 9-LYS(4.68%), 6-LYS(3.57%) 14-THR (2.96%), and 165-THR (1.72%). In CDK2- 4e, other notable interactions included 131-GLN (11.12%), 10-ILE (5.61%), and 20-LYS (2.13%), shown in Figure 7(g-h).

Additional DCCM analysis revealed the correlated motions of residues in CDK2 when bound to 4b and 4e. The overall similarity in the correlation patterns for both ligands suggests that the global dynamics of CDK2 remain relatively consistent regardless of the ligand bound. However, subtle differences in the strength and extent of these correlations may indicate ligand-specific effects on the local flexibility of certain protein regions, potentially influencing CDK2’s functional state, Figure 8(a) and (b).

Close

Figure 8.

Global dynamics and binding analysis of CDK2 with 4b and 4e. (a) DCCM analysis shows correlated motions of residues for both complexes. (b) PCA scatter plots illustrate the conformational changes and flexibility of CDK2 upon binding each ligand. (c-e) MMGBSA analysis comparing the binding free energy profiles of the two complexes, highlighting stability shifts and binding affinities over time, indicating that 4b does not stay bound to the CDK2 protein over the simulation period.

Export to PPT

PCA provides insight into the major conformational changes experienced by CDK2 in the presence of each ligand. The scatter plots show the distribution of simulation snapshots along the first two principal components (PC1 and PC2), representing the system’s dominant motions. The clustering and spread of these points indicate how the conformational landscape of CDK2 is modulated by each ligand. For CDK2 bound to 4b, the data points show a more defined clustering, suggesting that the protein explores a narrower range of conformations, implying a relatively stable binding mode. In contrast, for 4e, a wider spread of points is observed, which suggests greater conformational flexibility and possibly less stable binding. Greater conformational diversity may indicate that 4e allows CDK2 to sample a broader range of states, potentially influencing its enzymatic activity. If excessive flexibility is induced, it could lead to decreased binding affinity or altered functional dynamics of CDK2, Figure 8(c) and (d).

The MMGBSA analysis estimates the binding free energy for each ligand throughout the simulation. Lower (more negative) binding free energy values correspond to more potent and more favorable interactions between the ligand and CDK2. The CDK2- 4b complex exhibits a relatively stable binding energy profile, fluctuating around a more favorable range (∼ -35 to -25 kcal/mol); however, after the 70th ns, the ligand completely detaches from the protein, indicating a completely unstable binding and potentially meaningless use case of this ligand as an inhibitor for CDK2. In contrast, CDK2- 4e shows larger fluctuations and a less favorable average binding energy, as it starts at around -57 kcal/mol and then dips to -25 kcal/mol. It then gradually increases for the remainder of the simulation and reaches back to -45 kcal/mol. A differential shift around the same timeframe is observed where 4b loses all its interactions while 4e gains more free energy, indicating a stronger bound complex, Figure 8(e). Based on these observations, 4e serves as a potential inhibitor against CDK2 due to strong binding affinity and bound conformation at the protein’s active site.