The effect of novel 3R,6R-bis (4-hydroxy benzyl) piperazine-2,5-dione (BHBPPD) derivatives on the expression of caspases in gastric cancer: A molecular docking and dynamics simulation

2.1 Cell culture and gene expression study

2.2 Geo-analysis and finding common factors that are related to the hsa04210 kegg pathway

2.3 Molecular docking

Protein Data Bank (PDB) files corresponding to TNFRSF10B (ID: 1du3), GZMB (ID: 1iau), Jun (ID: 1jun), SPTA1 (ID: 1owa), APAF1 (ID: 1z6t), BID (ID: 2m5i), LMNA (ID: 2xv5), ACTB (ID: 3byh), FADD (ID: 3ezq), CFLAR (ID: 3 h13), CASP6 (ID: 3v6m), CYCS (ID: 3zcf), CASP3 (ID: 5i9b), CASP8 (ID: 5 l08), and CASP9 (ID: 1jxq) were retrieved from the Protein Data Bank server (https://www.rcsb.org). Subsequently, water and HETATOM components were removed from the structures. The three-dimensional structure (PDB) of the BHBPPD compound (Fig. 2) was visualized using Marvin Sketch software and optimized using Gaussian software with the B3LYP/6-311G (d,p) method. Autodock software version 4.2 was employed to conduct molecular docking studies. The BHBPPD compound ligand was docked with the protein mentioned above structures to determine optimal binding sites and the most stable free energy state of the ligand-receptor complex. A grid box was created for each receptor-ligand interaction according to the specifications outlined in Table 1. Docking simulations were performed with the Genetic Algorithm and Lamarckian GA parameters, running 200 molecular docking iterations. The autodock4 version of Linux was used to generate results files (dlg), which were subsequently analyzed (Morris et al., 2009). Discovery Studio Visualizer 2021 was utilized to quantify hydrophobic interactions and hydrogen bonds between receptors and ligands. Additionally, the type and number of amino acids in near to the binding sites were determined (Wallace, Laskowski, & Thornton, 1995).

Close

Fig. 2

Structure of BHBPPD.

Export to PPT

2.4 Molecular dynamics simulations

Molecular dynamics is a computational technique that simulates the physical movements of atoms and molecules as they interact over time. This method allows the study of atomic motion within a specified time frame, often focusing on the system's behavior far from equilibrium. To achieve meaningful results, the initial simulation step must guide the system toward stability, where thermodynamic and structural properties stabilize. This process is based on Newton's second law of motion for atoms. In this study, the molecular dynamics simulation of the target complex was conducted using Gromacs 2021.1 software. The simulation was conducted under constant pressure and temperature conditions (NPT) using the g54a7 force field and the SPC216 water solvent model. The simulation extended for 50 ns (ns). Two simulation models were chosen, each involving specific proteins with the highest binding energy from the docking results:

1. TNFRSF10B and CYCS are in their free state (without ligand).

2. TNFRSF10B and CYCS in the presence of BHBPPD ligand.

The latter mode included adding ligand information to the n.top and gro files, with itp and gro ligand files obtained from the Prodrg database, consistent with the Gromacs force field. The simulation process encompassed energy minimization, system neutralization, and balance steps. The steepest descent mining algorithm integrated 50,000 steps for energy minimization. Subsequently, the leap-frog algorithm was utilized for 150,000 steps to achieve balance using NVT (canonical ensemble) and NPT (isobaric-isothermal ensemble) ensembles. A Nose-Hoover thermostat and isotropic Monte Carlo (MC) barostat maintained a pressure of 1 atm. Coordinates were saved every ten ps during production runs, using a time step of 1 fs. The simulation was performed with three replications and reported.

For trajectory analysis and visualization, Visual Molecular Dynamics 1.9.3 (Humphrey, Dalke, & Schulten, 1996) was employed. 2D plots were generated to illustrate the intrinsic dynamical stabilities of the complex, including root-mean-square fluctuations (RMSF), root-mean-square deviation (RMSD), hydrogen bond interactions, and radius of gyration (Rg).

3.1 Cell proliferation and cell cytotoxicity assays

Following treatment with varying concentrations of BHBPPD, a dose-dependent reduction in the viability of the cancer cell line was observed (Fig. 3), while no such effect was observed in the normal cell line (Fig. 4). Our analysis indicates that the AGS cell line exhibited lower toxicity to BHBPPD at IC50 concentrations (below 45 μM) during drug treatments, whereas no significant toxicity was observed in the normal cell line at these concentrations. Primary cytotoxicity experiments demonstrated that the drugs hindered the growth of cancer cells compared to normal cells after 24, 48, and 72 h of treatment. Additionally, the cytotoxicity notably increased following the introduction of BHBPPD at various concentrations.

Close

Fig. 3

Cytotoxic effects of BHBPPD. AGS Cells were cultured for 3 times with different doses of drug. Cell viability was evaluated by MTT. The values represent the means of three independent experiments (mean ± SE).

Export to PPT

Close

Fig. 4

Cytotoxic effects of BHBPPD. HDF Cells were cultured for 3 times with different doses of drug. Cell viability was evaluated by MTT. The values represent the means of three independent experiments (mean ± SE).

Export to PPT

3.2 Apoptosis-related genes transcripts level

Following established standard protocols, relative transcripts of caspase 3, 9, 8, and GR were quantified using real-time PCR. Notably, the expression of CASP3 in AGS cells treated with BHBPPD exhibited a significant increase compared to the control at IC90 concentrations (P = 0.03 and 0.001, respectively) (Fig. 5A). Additionally, the expression of CASP9 showed a notable increase compared to the control (p = 0.005) in AGS cells treated with BHBPPD at IC90 concentrations (Fig. 5B). Treatment of AGS cells with BHBPPD resulted in a concentration-dependent increase in CASP8 expression at both 72 h and IC90 concentrations (Fig. 5C). Conversely, the expression of GR in BHBPPD-treated cells experienced a decrease (Fig. 5D). The study introduced a novel piperazine derivative that demonstrated a robust anticancer effect through the evaluation of cancer cell proliferation and the induction of caspase-dependent apoptosis. Caspases were positioned at the forefront of apoptotic signaling and recognized as the primary agents responsible for cleaving cellular proteins, a critical process for the degradation of dying cells. Caspases act as proteolytic enzymes that cleave cytoskeletal proteins upon activation (Sharifi, Hosseinali, Saboury, Szegezdi, & Falahati, 2019). Furthermore, the pivotal role of GR was highlighted in the cellular shift from GSH-to-GSSG, which can determine the fate of a cell. A decrease in GR gene expression was be linked with mitochondrial GSH/GSSG imbalance, consequently inducing apoptotic signaling.

Close

Fig. 5

Relative quantifications of CASP3, CASP9, CASP8, and GR transcripts in the AGS cell line. Cells were incubated with IC₅₀ and IC₉₀ concentrations of BHBPPD for 72 h. β-actin was used as a favored internal control. The results are presented as mean ± SE of three independent experiments each in triplicate compared with untreated samples. Symbols (***), (**), and (*) represent the mean differences between treated and untreated cell transcriptions as P < 0.001, P < 0.01, and P < 0.05 using one-way ANOVA, respectively.

Export to PPT

3.3 Identification of DEGs

In the analysis of the GSE4651 dataset, 340 gastric cancer samples and 218 normal samples were examined, while the GSE10262 dataset consisted of 18 samples, comprising 12 gastric cancer samples and six normal samples. Each chip's series were individually investigated using R software, and subsequently, differentially expressed genes (DEGs) were identified using adjusted P-value < 0.01 and logFC ≥ 2 or logFC ≤  − 2 as the cutoff criteria. In the context of GSE4651 analysis, a total of 3582 DEGs were identified, of which 2844 were upregulated genes and, 738 were downregulated (Fig. 6B). The presentation of the DEGs' fold change is depicted in Fig. 6A. Similarly, the analysis of GSE10262 revealed a total of 5859 DEGs, among which 3016 were downregulated and, 2843 were upregulated genes (Fig. 7B). The fold change performance of the DEGs is illustrated in Fig. 7A. Fig. 8 visually represents all the protein factors related to apoptosis, sourced from the hs02410 pathway of the KEGG database. The association between the 9441 genes obtained from the apoptotic protein factors shown in Fig. 4 was assessed using the STRING database. Ultimately, the protein factors that exhibited associations and were shared between the GSE4651 and GSE10262 datasets were selected and highlighted in Fig. 9.

Close

Fig. 6

(A) Heat map of DEGs. (B) Volcano plot of genes detected in Gastric cancer. Red means upregulated and downregulated DEGs; black means no difference (GSE4651 dataset).

Export to PPT

Close

Fig. 7

(A) Heat map of DEGs. (B) Volcano plot of genes detected in Gastric cancer. Red means upregulated and downregulated DEGs; black means no difference (GSE4651 dataset).

Export to PPT

Close

Fig. 8

Protein factors participating in apoptosis obtained from Kegg server with hsa0410 pathway.

Export to PPT

Close

Fig. 9

The final protein–protein interaction (PPI) network of the Caspas3, 8, and 9 with common genes that in apoptosis between GSE4651 and GSE10262 datasets.

Export to PPT

3.4 Molecular docking

Virtual initial screening was conducted to evaluate the potential interaction of the BHBPPD compound with a selection of key apoptotic factors, including TNFRSF10B, GZMB, JUN, SPTA1, APAF1, BID, LMNA, ACTB, FADD, CFLAR, CASP6, CYCS, CASP3, CASP8, and CASP9, using Autodock. The docking results provided the docking score represented as Gibbs free energy ΔG (kcal/mol), which accounts for the sum of all interactions between the receptor and ligand energies. Detailed docking results are outlined in Tables 1 and Supplementary S1. Notably, the BHBPPD ligand exhibited a ΔG of −9.6 kcal/mol when binding to TNFRSF10B, in comparison to −8.22 kcal/mol for CYCS. In the case of TNFRSF10B, hydrogen bonding occurred with Ans(L)199, Asp(G)95, Gln(L)205, and Arg(G)101 amino acids, while hydrophobic interactions were established with Glu(G)124, Arg(G)104, Ile(G)123, Arg(G)104, Thr(L)200, Asn(L)262, Lys(L)204, Ser(G)96, Asp(L)203, Lys(L)201, Phe(G)91, Asp(G)12, and Cys(G)125. Similarly, BHBPPD formed hydrogen bonds with CYCS, involving amino acids Lys(A)27, Tyr(A)97, Glu(A)104, Ala(A)101, and Lys(A)7. Additionally, hydrophobic and anionic bonds were established with Thr(A)102, His(A)33, Lys(A)22, Ile(A)11, Ser(A)15, His(A)18, and Thr19, along with an anionic bond with Glu(A)21 (Table 1). Furthermore, the study delved into the role of TNFRSF10A and TNFRSF10B as cell surface receptors that interact with tumor necrosis factor-related apoptosis-inducing ligand (TRAIL), mediating the extrinsic apoptosis pathway. However, the mechanisms underlying the transcriptional regulation of these receptors remain relatively uncharacterized (Li et al., 2015).

The TNFRSF10B protein, encoded by the gene mentioned above, belongs to the TNF-receptor superfamily, encompassing an intracellular death domain. Activation of this receptor may occur via tumor necrosis factor-related apoptosis-inducing ligand (TNFSF10/TRAIL/APO-2L), transmitting an apoptosis signal. FADD, a death domain with an adaptor protein, has been identified as necessary for apoptosis mediated by this protein. Various transcript variants encoding different isoforms, as well as one non-coding transcript, have been identified for this gene. Furthermore, DDIT3′s interaction with phospho-JUN results in the formation of the DDIT3-phospho-JUN complex, which binds to the AP-1 binding site within the TNFRSF10A promoter region. It was demonstrated that KAT2A interacts effectively with the N-terminal area of DDIT3, thereby regulating TNFRSF10a and TNFRSF10b expression in endoplasmic reticulum stress-mediated apoptosis in human lung cancer cells (Li et al., 2015). Apoptosis, a highly regulated form of programmed cell death, is pivotal in maintaining multicellular organism homeostasis and development. Cytochrome c serves as a crucial element in the activation of the apoptotic intrinsic pathway, initiating the caspase cascade through interaction with Apaf-1. Recent studies by Elena-Real et al. identified 14–3-3ε as a cytosolic cytochrome c target and demonstrated its role in caspase cascade activation by inhibiting Apaf-1 (Elena-Real et al., 2018). Additionally, Parray et al. investigated the interaction of polyethylene glycol with cytochrome c and suggested that specific concentrations of crowders might hold potential as therapeutic agents for various diseases (Parray et al., 2021). The docking outcomes of compounds targeting significant apoptotic factors, as highlighted by geo-analysis, indicated strong potential interactions with CYCS and TNFRSF10B-mediated apoptotic factors, making them promising candidates for therapeutic targeting.

3.5 Molecular dynamic (MD) simulations

Molecular dynamics simulations were conducted on the TNFRSF10B and CYCS proteins complexed with BHBPPD. The systems were solvated in a water box and simulated for 50 ns. Table 2 presents the average and standard deviations of kinetic (EKCMT) and potential (EPTOT) energies, total (ETOT) energy, volume (nm3), density (kg/m3), and enthalpy, which were analyzed to assess the system's behavior. The molecular dynamics analysis was carried out using GROMACS software. Various parameters such as solvent accessible surface area (SASA), hydrogen bonding, root mean square deviation (RMSD), and radius of gyration (Rg), as well as root mean square fluctuations (RMSF) at 300 K, were obtained for analysis (Table 3). The RMSD of the TNFRSF10B and CYCS systems in the presence of BHBPPD over a 50 ns interval indicated a stable state, with a more pronounced reduction in RMSD observed for CYCS when bound to the ligand. The results confirmed equilibrium in the MD simulations for the studied systems (Fig. 10A) (Gharaghani, Khayamian, & Ebrahimi, 2013; Suo et al., 2018). Table 2 provides insight into the energy levels of the four molecular dynamics simulation systems. The data reveals that BHBPPD can decrease potential energy levels and contribute to the overall stability of the simulation system, particularly for TNFRSF10B. Root mean square fluctuations (RMSF) of residues were compared among free TNFRSF10B and CYCS and their complexes with BHBPPD to assess residue mobility. Results shows that in the presence of a BHBPPD, the amount of RMSF decreases, and as a result, the two proteins TNFRSF10B and CYCS decrease their flexibility and become more stable, which is a significant difference in the amount of RMSF in the absence of ligand (Fig. 10 C&D). The radius of gyration (Rg) was used to gauge enzyme compactness, with data confirming a decrease in Rg for both TNFRSF10B and CYCS in the presence of BHBPPD.Rg of TNFRSF10B (3.17±0.19) and CYCS (3.05±0.21) in 50 ns time interval and it verifies the gyration radius of TNFRSF10B (2.20±0.05) and CYCS (2.33±0.09) decrease when BHBPPD is present (Table 3 and Fig. 10B). Hydrogen bonding interactions between BHBPPD and the proteins led to alterations in accessible surface area (ASA) and contributed to their structural stability (Fig. 10F) (Gan et al., 2019; G. Wang, Wang, Li, & Lu, 2011). BHBPPD can interact with TNFRSF10B (0.85±0.78 hydrogen bonds) and CYCS (2.90±0.92 hydrogen bonds) and alter the hydrogen bonds between TNFRSF10B (266.80±16.83 to 269.40±21.37) and CYCS (284.36±9.47 to 274.40±9.60) and the surrounding solution. The solvent-accessible surface area (ASA) analysis revealed a reduction in surface area for both TNFRSF10B and CYCS in the presence of BHBPPD, aligning with the decrease in Rg. Structural changes of TNFRSF10B in the presence of the ligand were more substantial compared to CYCS. Notably, the YIPF2-TNFRSF10B axis was associated with the malignant progression of lung cancer, and YIPF2 enhanced chemotherapeutic agent-mediated apoptosis by augmenting TNFRSF10B recycling to the plasma membrane in NSCLC cells (Falschlehner, Emmerich, Gerlach, & Walczak, 2007; Y. Wang et al., 2020).

Close

Fig. 10

Stability of TNFRSF10B and CYSC and complex they in the presence of BHBPPD. Average root means square deviation (RMSD) for protein fit to Alpha carbon (A) for 50 ns, and average root mean square fluctuation of whole protein following stabilization. (B) Radiues of gyration of TNFRSF10B and CYSC and TNFRSF10B and CYSC complex in the presence of BHBPPD (C and D) shows the RMSF values the TNFRSF10B (C) and CYSC (D) form and ligand (BHBPPD) bound forms of the proteins. (F) Indicates the number of hydrogen bonds between TNFRSF10B and CYSC proteins and water solvent (spc216) in the presence of ligand and absence.

Export to PPT

Proteins TNFRSF10B and CYCS assume pivotal roles within the intricate framework of the apoptosis pathway. TNFRSF10B is a cell membrane receptor that interacts with the apoptosis-inducing ligand TRAIL (Tumor Necrosis Factor-Related Apoptosis-Inducing Ligand). Upon TRAIL binding to TNFRSF10B, an orchestrated series of molecular events is set in motion, triggering a cascade that culminates in activating the apoptotic pathway within the cellular milieu. Additionally, CYCS, an essential constituent of the cytochrome C complex, plays a decisive role in endogenous apoptosis. It is actively engaged in fundamental processes, encompassing the liberation of mitochondrially stored factors and orchestrating a multi-component structure, the apoptosome, formed through the assembly of apoptotic protease activating factor-1 (Apaf-1). These functions collectively contribute to the meticulous regulation of programmed cell death (Farooqi, Gadaleta, Ranieri, Fayyaz, & Marech, 2017; Kretz et al., 2018). Table 4 structural changes of proteins TNFRSF10B and CYCS in the presence and absence of ligand BHBPPD during 50 nanoseconds of the simulation time. The results showed that the structural changes of TNFRSF10B in the presence of ligand are greater than CYCS, so the value of the structure is equal to 147.58 for TNFRSF10B and 162.99 value when ligand is present. The largest change in TNFRSF10B structure was related to changes in the alpha-helix and coil levels (Fig. 11)

Close

Fig. 11

Shows probability of various secondary structure elements by residue of TNFRSF10B and CYSC proteins in the presence of ligand (BHBPPD) and absence.

Export to PPT

The combined analysis of molecular dynamics simulations, energy levels, RMSD, RMSF, Rg, ASA, and their influence on structural stability sheds light on the interactions between BHBPPD and TNFRSF10B and CYCS proteins, providing valuable insights for potential therapeutic targeting.