Investigating chlorothiophene-based chalcones as anticancer agents: Synthesis, biological evaluation, and computational insight via the p53-targeted pathway

2.1. General procedure of the synthesis of chlorothiophene-based chalcones

A mixture of 2-acetyl-5-chlorothiophene or 2-acetyl-4,5-dichlorothiophene (0.01 mol) and the corresponding aromatic aldehyde (0.01 mol) in methanol (20 mL) was stirred for 24 h in the presence of 40% KOH (4 mL). The reaction was monitored by thin layer chromatography (TLC) with the eluent n-hexane/ethyl acetate (7:3). The workup procedure was performed by pouring the mixture into crushed ice and acidifying it with 5% HCl. The resulting solid was washed and purified using a recrystallization process with ethyl acetate and then dried in a desiccator. All chemicals and reagents were obtained from both Sigma-Aldrich Co. LLC and Shanghai Macklin Biochemical Co., Ltd. and utilized without further purification. Electrothermal’s IA9100 melting point apparatus was used to determine melting points (m.p.). Purity level analysis was conducted using a Dionex Ultimate 3000 system with a C₁₈ column and a UV detector set at 254 nm. The mobile phase consisted of 0.1% formic acid in water (v/v) (Phase A) and 0.1% formic acid in acetonitrile (v/v) (Phase B). The gradient program was structured as follows: starting with 32.5% B in A, transitioning to 40% B in A at 4 mins, maintaining 40% B in A until 8 mins, increasing to 45% B in A at 10 mins, reaching 95% B in A at 15 mins, holding at 95% B in A until 18 mins, and then returning to the initial 32.5% B in A by 19 mins, maintaining this condition until the run ended at 35 mins. The flow rate was set at 1 mL/min. The injection volume for the analysis was 10 µL, with samples prepared at a concentration of 1 mg/mL, dissolved in acetonitrile [22]. Molecular mass was obtained using gas chromatography with mass spectrometry (Shimadzu QP-2010S). ¹H NMR (500 MHz) and ¹³C NMR (125 MHz) spectroscopy were performed using the JEOL JNM-ECA 500 NMR spectrometer. The chemical shifts (δ) were expressed in parts per million (ppm) downfield from tetramethyl silane as an internal standard, and the coupling constant (J) was expressed in Hertz (Hz). Infrared spectra were obtained from PerkinElmer^® spectrum two attenuated total reflectance-infrared (ATR-IR) spectrometer. TLC on silica gel Al foil was visualized under a UV lamp.

2.2. Anticancer activity

Approximately 10,000 cells per 100 µL were seeded into 96-well microplates and incubated in a CO₂ incubator for 24 h. The media were removed, and the cells were first treated with a serial dilution solution of the compounds and then incubated for 24 h. Subsequently, 100 µL of an 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyl-2H-tetrazolium bromide (MTT) solution was added to each well. The plate was left to incubate for 4 h. To prevent crystal formation, 100 µL of sodium dodecyl sulfate (SDS) stopper solution was added to each well, and the plates were left overnight. Following the complete dissolution of all crystals, the absorbance (A) was measured at 595 nm using a microplate reader. Then, the percentage of dead cells was calculated using Eq. (1):

The obtained percentage of cell death is transformed into probit numbers, creating a curve of log concentration (x) vs. probit number (y) using the equation y = ax + b. When y equals 5, the corresponding x value represents the Log IC₅₀. The IC₅₀ is determined by calculating the antilog of the x value, representing the minimum compound concentration required to induce 50% cell death. The selectivity index was calculated using Eq. (2):

2.3. Computational methods

3.1. Structural design

Recent advancements in anticancer drug design have placed considerable emphasis on modifying the structure of molecules containing heteroaryl substitutions. Bio-isosteric replacement has emerged as a key strategy to retain comparable properties within heteroaromatic systems while enhancing biological activity [54]. Our research has been consistently dedicated to analyzing and developing heteroaryl chalcones and their derivatives, specifically those containing five-member rings. We observed a trend indicating that thiophene significantly enhances the toxicity of a compound to cancer cells compared to other five-member rings heterocycles such as pyrrole and furan, especially against HeLa cervical cancer cells [55-57]. Thiophene is commonly employed as a substitute for benzene in drug discovery because of its recognized similar aromaticity, rendering it an effective mimic of the phenyl moiety.

Moreover, our previous research considered the impact of two substituents, i.e., chloro and methoxy groups, as crucial functional groups influencing the toxicity and selectivity of the compound [58]. Several studies have investigated the introduction of substituents in thiophene and have observed alterations in activity [59-61]. Incorporating a chlorine atom has been noted to significantly improve the stability, potency, and specificity of the binding site. This enhancement arises from the capability of the chlorine atom to fit into tight and deep hydrophobic pockets within biological targets, thereby facilitating stronger and more selective interactions [62,63]. This attribute is particularly advantageous because our binding target consists of hydrophobic pockets. However, achieving selectivity is a significantly challenging task. Hence, we explore other substituents to optimize selectivity, aiming for compounds that remain inactive in normal cells. Although a comprehensive study of anticancer activity on this scaffold has not yet been conducted, we consider it the most promising starting point. Therefore, we focus further optimization efforts on this framework.

3.2. Chemistry

The chalcone scaffold is highly regarded as a versatile pharmacophore in medicinal chemistry for designing diverse therapeutic agents. In this study, we expanded upon our previous work [18] with the chalcone scaffold by successfully synthesizing a series of chlorothiophene-based chalcone derivatives. The synthetic route employed is straightforward, as illustrated in Figure 2.

Close

Figure 2.

Reaction scheme.

Export to PPT

The required starting material was the chloro-substituted acetyl thiophene (2-acetyl-5-chlorothiophene and 2-acetyl-4,5-dichlorothiophene) and appropriate aldehydes with methoxy and chloro substituents. Then, chalcone compounds C1–C9 were synthesized using base-catalyzed aldol with aqueous KOH (40%) in the presence of methanol as a solvent at room temperature by stirring. This condition is sufficient to produce high-yield products of chalcones, as shown in Table 1. The chemical features of the product analogs were determined using spectrophotometric techniques.

Both ¹H NMR and ¹³C NMR spectroscopy were used to analyze the synthesized chalcone compounds. These compounds have a conjugated enone group at the center of the molecule. The alkene C=C double bond positions are designated as the α and β positions, with the carbonyl C=O group referred to as the α position and the carbonyl C=O group referred to as the β position. In all compounds, the protons corresponding to these positions appeared as doublets because of their coupling with each other, with the coupling constant (J_Hα–Hβ) ranging from 15.4 to 15.7 ppm, indicating a trans configuration of alkene proton [64,18]. The carbon in the C=O group typically appeared in the most downfield region of all-carbon signals, with chemical shifts ranging from 180.38 to 181.88 ppm. The IR spectra also support the presence of this carbonyl system by showing strong absorbance at approximately 1,660 cm⁻¹, indicating a stretching vibration of the C=O carbonyl bond.

By analyzing the yield trend, we observed that benzaldehyde substitution at position 4 (para-position) resulted in higher yields compared with benzaldehyde substitution at position 2 (C4 vs. C7). Moreover, aldehydes with electron-withdrawing chloro substituents exhibited greater yield than those with electron-donating methoxy groups, attributed to an electron-withdrawing effect that enhances benzaldehyde reactivity. Conversely, the introduction of an additional chloro group to the thiophene system led to decreased yield (C4 vs. C8). The dichloro group in the thiophene system exerted a significant electron-withdrawing influence via an inductive effect, thereby reducing the reactivity of acetyl during the reaction.

3.3. Anticancer activity

Cytotoxicity tests were conducted using assays on various predominant types of cancer cells to assess the anticancer effects of chalcone compounds. Doxorubicin, cisplatin, and 5-fluorouracil (5-FU) were utilized as reference drugs. Table 2 shows that most compounds exhibit low IC₅₀ values against WiDr colorectal cancer cells, surpassing the efficacy of the positive control (i.e., 5-FU). Compounds C4 and C6 exhibited the most promising IC₅₀ value for this cell line at 0.77 and 0.45 µg/mL, respectively. However, C6 was observed to have moderate toxicity to normal cells, whereas C4 was non-toxic.

C1 exhibited low IC₅₀ values across all cancer cell lines, with moderate toxicity toward normal cells. The observed trend was the addition of one methoxy substituent (C2–C4) led to reduced toxicity toward normal cells, indicating a better selectivity index. However, this trend was not observed with the addition of the dimethoxy group (C5). Compounds with multiple chloro substituents did not yield favorable results, e.g., C6. Similarly, C8 and C9 exhibited more toxicity to normal cells. Thus, the optimal conditions lead to good activity toward WiDr colorectal cancer cells, as shown in Figure 3.

Close

Figure 3.

Optimized structures of compounds that potentially show good selectivity.

Export to PPT

3.4. Molecular docking

The docking protocol commenced with the validation process using the native ligand. The redocking procedure was conducted with a predefined grid box size, yielding RMSD values of 0.385 Å for MDM2 and 0.583 Å for Bcl-2 (Table S1). To validate the adequacy of the grid box size used in the docking procedure, our approach includes employing an RMSD threshold of <2 Å. Figure S1 illustrates the stacking condition of the native ligand before and after the validation process. For further analysis, we selected two samples, i.e., C4 and C6, known for their cytotoxicity against WiDr colorectal cancer cells.

In the current investigation, Table 3 shows the binding affinity of the compounds, and the protein-ligand interactions compared with each native ligand. The binding affinity of both compounds to both MDM2 and Bcl-2 proteins did not exhibit higher values than their respective native ligands. In complex to MDM2, the binding affinity values of −6.2 kcal/mol for C4 and −6.6 kcal/mol for C6 were observed. In comparison, the native ligand nutlin-3a has a value of −8.8 kcal/mol, which indicates that our compound has weaker binding than nutlin-3a. This trend is similarly reflected in the Bcl-2 protein, with C4 and C6 exhibiting values of −6.4 and −6.6 kcal/mol, respectively, whereas the native ligand has a value of −8.5 kcal/mol. C6 exhibits a slightly higher value than C4, indicating potentially stronger binding strength, as evidenced by its smaller IC₅₀ value.

In addition to assessing binding energies, we also evaluated the protein-ligand interactions that are essential for inhibiting the selected proteins. This comprehensive approach enabled us to consider the specific molecular interactions that contribute to the inhibition of target proteins, providing valuable insights into the efficacy of potential inhibitors. Before delving into further interactions, we need to comprehend the binding area of the p53 peptide to the MDM2 protein. The MDM2 domain consists of structural repetitions, comprising two portions, each spanning 45 and 38 amino acids, with similar structures (residues 26–70 and 71–108) which are linked by an approximate dyad axis of symmetry. The crystal structure of MDM2 bound to a p53 peptide elucidated the significance of Phe19, Trp23, and Leu26 as key residues, establishing a locally confined, hydrophobic interaction site within p53. Moreover, the MDM2 cleft comprises 14 hydrophobic and 11 aromatic amino acids, facilitating numerous van der Waals contacts with p53, including those with His96, Ile99, and Tyr100 [65]. The compact nature of the p53-MDM2 interface has facilitated the design of not only peptides but also small organic compounds capable of blocking this interaction. A compound such as nutlin-3a, along with other established cancer drugs, exert their inhibitory effects by targeting the hydrophobic cleft formed by specific amino acids, such as Leu54, Gly58, Ile61, Val93, His96, and Ile99 [66,67]. This targeted binding disrupts the interaction between MDM2 and p53, thereby restoring the tumor-suppressive functions of p53 and offering promising therapeutic strategies for cancer treatment.

Figure 4(a) and 4(b) present the 3D visualization of docked compounds within the binding pocket of MDM2. This docking protocol shows that both compounds lack hydrogen bonding interactions, like nutlin-3a. Instead, both compounds interact electrostatically (Pi-cation) with the amino acid His96, with distances of 3.974 and 4.188 Å for C4 and C6, respectively. Moreover, the thiophene moiety contributes to a Pi-sulfur interaction with Met62. Hydrophobic contacts are observed with Leu54, Ile61, Val93, and Ile99. These hydrophobic interactions play a crucial role in anchoring compounds to the binding site and contribute to the overall stability of the protein-ligand complex.

Close

Figure 4.

3D interaction visualization of murine double minute 2 (MDM2) complexed with (a) C4 and (b) C6. 3D interaction visualization of B-cell lymphoma 2 (Bcl-2) complexed with (c) C4 and (d) C6. Note: Ala: alanine, Arg: arginine, His: histidine, Leu: leucine, Ile: isoleucine, Met: methionine, Phe: phenylalanine, Tyr: tyrosine, Val: valine.

Export to PPT

The Bcl-2 proteins fold into a distinctive multiple α-helices bundle [68]. The hydrophobic groove within the anti-apoptotic protein Bcl-2, characterized by the presence of conserved Bcl-2 homology domains BH1 and BH2, serves as a critical binding site for regulating apoptotic pathways. This specialized region is pivotal in mediating protein-protein interactions and modulating cellular responses to apoptotic stimuli. This groove works like a hiding spot, trapping the BH3 domain of pro-apoptotic Bcl-2 proteins. This BH3 domain plays a critical role in correctly folding a hydrophobic pocket within the Bcl-2 protein. To disrupt the function of Bcl-2, an approach involving small molecules that can bind to this groove area is utilized. These small molecules target one or more P1–P4 sub-pockets within the BH3 groove of anti-apoptotic proteins [69]. By binding to these sub-pockets, these molecules induce the release of pro-apoptotic BH3-only proteins, which subsequently activate the apoptotic pathway, leading to cell death. Polar residues Tyr108 and Arg146 act as important anchor points, helping to stabilize the attached ligands [70,71].

Figure 4(c) and 4(d) show the 3D visualization of the docked compound inside the binding pocket of Bcl-2. Interestingly, neither the docked compound nor the native ligand forms hydrogen bonds in the complexes. However, the crucial amino acid Arg146 is observed to interact hydrophobically with the docked compound, supported by interactions with Phe112, Met115, and Val133. In addition, the thiophene moiety contributes to a Pi-sulfur interaction with Phe104. These interactions indicated that our compound had moderately bound to the Bcl-2 protein. Further exploration of this protein-ligand interaction will be discussed more comprehensively through MD simulation studies.

3.5. Molecular dynamics simulation

After the docking investigation, MD simulation was utilized to evaluate the dynamic stability of the compound in a solvent system. C4 was chosen as the most promising candidate for MD simulation because of its high toxicity and better selectivity. The simulation was conducted over 500 ns, with its trajectory result compared with their apoproteins. To validate the MD simulation and free binding energy calculation method, we attempted to replicate the system used by Kumar et al. [72] which employed AMBERFF. The summarized results are presented in Table S2.

Figure 5(a) depicts the stability of the C4-MDM2 complex, exhibiting a consistent RMSD profile within the range of 0.10–0.20 nm for the initial 300 ns, followed by an increase to 0.40 nm toward the end of the simulation. This deviation is noticeably larger than that of the stable RMSD of the apoprotein throughout the simulation. Interestingly, previous studies focusing on the development of potential ligands for MDM2 via MD simulation have shown a similar trend up to a maximum of 200 ns, which is in the same range of RMSD values [73,74]. However, the results obtained up to 500 ns provide new insights for analysis, offering information about events occurring after 300 ns of simulation. Figure 5(b) shows that the C4-Bcl2 complex exhibits a more consistent RMSD trend, like its apoprotein counterpart, and achieves equilibrium within the first 100 ns, with an RMSD increment from 0.10 nm to 0.40 nm. Compared with numerous studies of Bcl-2 complex MD simulations, typically spanning 100–200 ns, the range of RMSD values of our study was considered acceptable [75,76]. The stability persisted until the end of the simulation. This behavior was also noted in the apoprotein simulation, indicating that the introduction of the ligand did not significantly alter the protein backbone.

Close

Figure 5.

Root mean square deviation (RMSD) superimposition of the (a) backbones of murine double minute 2 (MDM2) (black) and C4-MDM2 (red), (b) the backbones of B-cell lymphoma 2 (Bcl-2) (blue) and C4-Bcl-2 (orange), (c) ligand C4 complexed with MDM2 (red), and (d) ligand C4 complexed with Bcl-2 (orange).

Export to PPT

The RMSD value of the ligand C4 complexed with MDM2 ranges from 0.05 to 0.15 nm (Figure 5c), with a significant change occurring at 100 ns. This finding is attributed to the formation of a hydrogen bond between the carbonyl group of the chalcone moiety and the amino acid residue Tyr67. Subsequently, the RMSD remains stable until the end of the simulation. Similarly, the RMSD value of ligand C4 complexed with Bcl-2 exhibits fluctuations within the same range (Figure 5d), i.e., it fluctuated along the simulation, with the highest deviation occurring at approximately 150 and 350 ns, reaching 0.18 nm.

The RMSF pattern of the protein backbone of the C4-MDM2 complex, depicted in Figure 6, exhibits fluctuations within the range of 0.10–0.20 nm. A significant deviation up to 0.23 nm is observed, which is attributed to the formation of a hydrogen bond between the carbonyl group of the chalcone and Tyr97, consistent with the RMSD profile findings. The RMSF pattern of the protein backbone of the C4-Bcl2 complex closely resembles the apoprotein. Minor discrepancies are noticed around amino acid residues 130–140 ns, particularly in the binding pocket area, which will be further explained by the DSSP result. In summary, the RMSF profile of the protein backbone during the simulation indicates that the introduction of the ligand has minimal impact on the protein structure.

Close

Figure 6.

Root mean square fluctuations (RMSF) superimposition of the (a) backbones of murine double minute 2 (MDM2) (black) and C4-MDM2 (red) and (b) backbones of Bcl-2 (blue) and C4-Bcl-2 (orange).

Export to PPT

The RoG serves as a vital metric for assessing the structural stability of proteins and reflecting their compactness level. Large RoG values indicate a protein molecule with a loosely packed structure, whereas small RoG values indicate a protein molecule with a tightly packed structure. In the case of the MDM2 protein, the RoG values ranged from 1.20 nm to 1.35 nm (Figure 7a). An increase in RoG values was observed after a simulation time of 300 ns, indicating significant structural changes consistent with the RMSD profile findings. This result indicates that the simulation of the C4-MDM2 complex may not be stable beyond 300 ns. The RoG values of the protein backbone of the C4-Bcl2 complex ranged from 1.40 to 1.55 nm and remained stable throughout the simulation (Figure 7b). This finding supports the notion that the C4-Bcl2 complex has greater stability and compactness than the C4-MDM2 complex.

Close

Figure 7.

Radius of gyration (RoG) superimposition of the (a) backbones of murine double minute 2 (MDM2) (black) and C4-MDM2 (red) and (b) backbones of B-cell lymphoma 2 (Bcl-2) (blue) and C4-Bcl-2 (orange). Solvent accessible surface area (SASA) superimposition of the (c) backbones of MDM2 (black) and C4-MDM2 (red) and (d) backbones of Bcl-2 (blue) and C4-Bcl-2 (orange).

Export to PPT

Determining the SASA parameters is pivotal for understanding the structural alterations occurring in proteins during the binding process. The SASA values, quantified in nm², provide insights into the flexibility of protein structures. An increase in SASA signifies the expansion of receptor macromolecules, whereas a decrease in SASA indicates greater compactness. In the case of the C4-MDM2 complex, the SASA values ranged from 45 to 54 nm² (Figure 7c). By contrast, in the case of the C4-Bcl2 complex, the SASA values ranged from 78 to 86 nm² (Figure 7d). Both were consistent during the simulation and reflected stability comparable to the apoprotein simulation.

Hydrogen bonds are pivotal in facilitating the formation of protein-ligand complexes, thereby influencing the stability and specificity of interactions. Figure 8(a) and 8(b) illustrate the bonding patterns observed in both complexes. The C4-Bcl-2 complex had a larger number of potential hydrogen bonds, i.e., a total of 4, than the C4-MDM2 complex, which exhibited only two potential hydrogen bonds. Specifically, in the C4-MDM2 complex, hydrogen bonds involving interactions with amino acids Tyr67 and Lys94 were captured at 167 ns (Figure 9a). The crucial amino acid His96 was weakly involved, primarily via van der Waals interactions. Conversely, in the C4-Bcl-2 complex, three potential hydrogen bonds were captured at 194 ns (Figure 9b), involving amino acids Leu137, Asn143, and Arg146, which are considered crucial for binding. The interactions involved the carbonyl moiety of the ligand and hydrogen atoms in the amino acid residues. Overall, the hydrogen bond analysis confirms that C4 exhibits a stronger binding affinity for Bcl-2 than MDM2 protein, highlighting its potential as a promising candidate for further investigation in anticancer drug development.

Close

Figure 8.

Number of hydrogens detected in the complexes of (a) C4-MDM2 and (b) C4-Bcl-2.

Export to PPT

Close

Figure 9.

Visualization of hydrogen bonds in (a) C4-MDM2 and (b) C4-Bcl-2 (green lines). Note: Leu: leucine, Tyr: tyrosine.

Export to PPT

The stability of the protein structures was evaluated using the RAM plot, which visually represents the statistical distribution of amino acids. For the MDM2 protein (Figure 10a), 88.9% (72 amino acids) of the 81 residues were in the favored region (red area), indicating a favorable conformational arrangement. Moreover, 11.1% (9 amino acids) was detected in the additional allowed region (yellow area), indicating minor deviations from the ideal conformation. Similarly, for the Bcl-2 protein post-MD simulation, the RAM plot revealed that 85.8% (115 amino acids) of the 134 residues were situated in the favored region (Figure 10b). Furthermore, 12.7% (17 amino acids) was observed in the additional allowed region, signifying acceptable structural deviations. Only 1.5% (2 amino acids) were in the generously allowed region, indicating minimal structural strain. None of the amino acid residues in each protein were in the disallowed area (white area) after the MD simulation, indicating the overall structural stability of the proteins under investigation.

Close

Figure 10.

Ramachandran plot (RAM) plots of (a) MDM2 and (b) Bcl-2, and DSSP plots of (c) MDM2 and (d) Bcl-2, and FEL plots of RMSD vs RoG after simulation of (e) MDM2 and (f) Bcl-2 complexed with C4 after MD simulation.

Export to PPT

Supporting the RMSD and RMSF profile findings, Figure 10(c) and 10(d) illustrate the DSSP plots of both MDM2 and Bcl-2 proteins after MD simulation using C4. The increase in RMSD value after 300 ns is confirmed to be influenced by the loss of helical structure around amino acid residues 18–25, detected by a change in color from blue to yellow in the DSSP spectrum. As shown in Figure S2a, the helical structure loosens into a coil structure. The crucial amino acid residue His96 also shows slight changes, as is observed at 343 ns, when the helical structure begins to loosen. These two changes are predicted to contribute to the fluctuation of the RMSF profile. In the DSSP profile of the C4-Bcl-2 complex, the essential binding area is dominated by α-helical loops. In the binding pocket of Arg146, the structure remains consistent throughout the simulation, indicating no structural change around this area. However, slight changes were observed in the amino acid Tyr108 area, and helical loop opening was observed in amino acid residues above 200 (Figure S2b).

Using the RoG and RMSD values extracted from the trajectories, we constructed a contour map of the FEL, enabling us to pinpoint a dominant cluster encompassing all complexes. In Figure 10(e) and 10(f), the blue area denotes the minimum free binding energy state, whereas the cyan, green, and yellow areas denote the metastable conformational states. From each complex, a single frame was selected from the cluster with the lowest energy level, ensuring acceptable structural deviations. The C4-MDM2 complex has the lowest free binding energy structures (∆G = 0) with confirmation located around RoG = 1.24519 nm and RMSD = 0.12813 nm. This confirmation was captured from the frames at 124.64 ns, where the RMSD is at a local minimum energy, indicating high consistency with its starting point. However, after 300 ns, the increasing RMSD indicates structural changes, including the opening of helices in the crucial amino acid residue His96 area and other residues in the tail side from 15 to 18. This finding indicates that, after 300 ns, the introduction of ligand C4 compromises its stability.

Meanwhile, the C4-Bcl2 complex exhibits the lowest free binding energy structure (∆G = 0) with conformation at RoG = 1.45387 nm and RMSD = 0.32822 nm. The local minima with RMSD = 0.3822 nm were captured at 240.75 ns, indicating structural changes around amino acid residues 200–211, which are not within the crucial binding pocket area. This change starts occurring at 90–100 ns, consistent with the RMSD profile. Meanwhile, the rest of the protein structure remains stable, including the binding pocket area, indicating that the introduction of the ligand did not significantly alter the protein structure. This stability is observed until the end of simulation. All trajectories, including secondary structure analyses and FEL interpretations, confirmed the stability of the C4-Bcl2 complex compared with the C4-MDM2 complex.

3.6. MM-PBSA calculation

The MM-PBSA method, known for its advanced capabilities in calculating free binding energy, was employed to validate the results obtained from molecular docking. Using this methodology, we computed the free binding energy (∆G_bind) of the anticancer C4-MDM2 and C4-Bcl2 complexes. Then, the total energies of these complexes were computed, and the results are presented in Table 4. The calculation was taken from the last 50 ns of the entire simulation. The binding energy of the C4-MDM2 complex was determined to be −21.26 kcal/mol. This finding falls within the range of free-binding energies reported by Li et al. [77] for the p53-MDM2 complex (−28.93 kcal/mol) and another peptide inhibitor (−28.22 kcal/mol), as calculated using a similar MM-PBSA method. In comparison to that of the C4-MDM2 complex, the free binding energy of the C4-Bcl-2 complex was calculated to be −15.70 kcal/mol. This value is close to the reported free binding energy of a small molecule inhibitor, i.e., α-muurolene, investigated for its inhibitory action on Bcl-2, which was reported to be −10.33 kcal/mol [78]. The difference indicates that the binding site may favor the formation of complexes with the ligands.

3.7. DFT and MEP analyses

The DFT calculations were performed in the gas phase using the B3LYP method with a 6-31G basis set. These calculations provided crucial parameters for analyzing the chemical reactivity of the synthesized compounds, such as energy gap, electronegativity (χ), chemical hardness (η), and electrophilicity index (ω). Frontier molecular orbitals encompass the E_HOMO and E_LUMO, which are pivotal quantum chemical parameters governing the reactivity of molecules and are utilized to compute various chemical reactivity descriptors. A high difference (∆E_gap) between the HOMO and LUMO energy values indicates low reactivity of the compound, resulting in higher stability. By computing the energy gap, C6 was determined to have the most substantial gap with a value of ∆E_gap = 4.0129 eV, indicating its superior stability relative to all other compounds examined. Following that, C1, C7, and C4 were ranked in decreasing order of stability. Interestingly, the trend indicates that the incorporation of a methoxy substituent in the phenyl group decreases the energy gap. The addition of another chloro group to thiophene, as observed in C8 and C9, further decreases the energy gap, resulting in the lowest energy gap in these two compounds. This finding indicates that the additional chloro group in thiophene destabilizes the system and renders the compound more reactive. The electronegativity (χ) tendencies indicate that C6 and C7 exhibit the highest electron-attracting capabilities, whereas C5 exhibits the lowest electronegativity. The additional chemical reactivity descriptors are summarized in Table S3.

To validate the evidence regarding the inhibitory properties of the compound, the MEP/ESP contour plot provides insights into the size, shape, and identification of reactive sites based on the electrostatic potential of the molecule. Figure S3 depicts the molecular electronic potential of all proposed compounds. Across all examined compounds, the red regions encircling the oxygen of the α, β unsaturated system, signify areas with the most substantial negative charge, indicative of high electron density available for electrophilic attack. Interestingly, a subtle increase in electron density is noted at sites where chlorine is substituted in the thiophene ring, indicating an enhanced capacity for interaction with electrophiles because of the addition of chlorine. However, oxygen from the carbonyl group remains the central focus of the interaction. The addition of one methoxy group results in a negative electron density around the oxygen. However, no significant increase in these red areas using dimethoxy is detected, as observed in C5 and C9. This finding confirms that the best compound design involves one substituent in the phenyl moiety (refer to Figure 3).

3.8. ADMET estimation

Drug-likeness is a critical attribute for compounds to become effective therapeutic candidates. Two approaches were employed to assess the drug-likeness of the compounds, namely, Lipinski’s rule of five and Veber’s criteria. Lipinski’s rule of five is based on five fundamental principles, such as (i) MW < 500 g/mol, (ii) lipophilicity or partition coefficient (Log P) < 5, (ii) TPSA < 140 Ų, (iv) number of HBD < 5, and (v) number of HBA < 10. Molecules meeting at least three rules were anticipated to exhibit favorable absorption properties. According to Veber’s criteria, compounds possessing ≤10 rotatable bonds and a polar surface area of ≤140 Ų are predicted to exhibit good oral bioavailability. This study used data from SwissADME predictions to generate information on several drug-likeness parameters.

Table S4 presents calculated molecular descriptors of all nine compounds. Remarkably, no violations were detected for either Lipinski’s rule of five or Veber’s criteria in any of the compounds. MW ranged from 248.73 to 343.22. However, the two comparison ligands, i.e., nutlin-3a and THIQ-PP, exhibited the maximal violation of Lipinski’s rule of five (>500). The number of HBD, number of HBA, TPSA, and Log P values all met the criteria. The Log P values ranged from 2.69 to 3.62. The observed changes in values indicate that the addition of a methoxy group reduces lipophilicity, making the compound more polar than the unsubstituted version (C1). The introduction of a chloro group (C6 and C7) increases lipophilicity. Based on this investigation, we can predict that the compounds may have better bioavailability than nutlin-3a and tetraisohydroquinoline phenyl pyrazole (THIQ-PP), which exhibit some violations and indicate poor bioavailability.

The evaluation of the ADMET parameters is a crucial step in drug development, ensuring the suitability of inhibitors for administration in biological systems. These assessments provide insights into how drugs are absorbed, distributed throughout the body, metabolized, and excreted and their potential for toxicity. To perform these evaluations, we utilized pkCSM, a comprehensive online tool known for its reliability in predicting ADMET properties. The results of these evaluations are summarized in Table S5, shedding light on various aspects of the pharmacokinetics and potential safety profiles of the compounds.

Initially, we evaluated GIA to assess the suitability of the compound for oral administration. The evaluation of the GIA thoroughly explains how a drug dissolves, traverses the gut/lumen interface, diffuses across membranes, and eventually gets eliminated from the bloodstream. Orally administered drugs primarily traverse gastrointestinal membranes via passive diffusion, a process influenced by the physicochemical properties of the molecules. Our compound exhibited a high GIA value (>90%), indicating that it will be absorbed into the bloodstream at a rate exceeding 90% from the intestinal tract.

The inadvertent penetration of the brain by drugs, especially those with CNS activity, can pose significant risks. Hence, the potential permeability of chemical compounds across the BBB needs to be evaluated. Our ligand exhibits high values for BBB and CNS permeability, ranging from 0.26 to 0.49, which are different from the two comparison ligands, with values of −1.59 and −0.92 for nutlin-3a and THIQ-PP, respectively. The addition of a methoxy group decreases these permeability values.

In drug metabolism, our compound is predicted to be a potential inhibitor of certain CYP isoenzymes, namely, CYP1A2, CYP2C19, and CYP3A4. Interestingly, some compounds with methoxy groups still exhibit inhibition activity toward CYP2D6, although none show such activity specifically for CYP2D6. The excretion parameters (total clearance) of our compound ranged from 0.05 to 0.42. Compounds C6 and C7 exhibit significantly lower total clearance than the other compounds, including the comparison ligands. As a result, compounds C6 and C7 are expected to necessitate shorter dosing intervals. In addition, we assessed the potential hepatotoxic effects to evaluate liver damage upon administration. Encouragingly, none of our compounds exhibited potency for this effect, indicating their overall safety profile.