Introduction

Breast cancer (BC) is the predominant form of cancer in women worldwide. In 2020, it exceeded lung cancer to become the primary cause of cancer occurrence worldwide, with around 2.3 million new cases, accounting for 11.7% of all cancer occurrences1. In solid tumors such as BC, the usual oxygen level is less than 2%, creating an area with low oxygen known as the hypoxic region2. Hypoxia enhances tumor plasticity and heterogeneity, leading to a more aggressive and metastatic phenotype that causes resistance to therapy and worsens the prognosis for patients. The hypoxia-inducible factor (HIF) family, which is essential for responding to low oxygen levels, comprises transcription factors vital for regulating cellular responses to hypoxic stress3,4. Currently, three forms of HIF have been discovered. HIF is a complex made up of two different types of subunits: an oxygen-sensitive α subunit (HIF-1α, HIF-2α, and HIF-3α) and an oxygen-insensitive β subunit (HIF-1β)5. Specifically, HIF-1α, also called the “master regulator,” is responsible for initiating genes linked to angiogenesis, cellular proliferation and viability, invasion and proliferation of cancer, glucose metabolism, evasion of the immune system, and resistance to multiple cancer treatments6,7. Multiple studies have consistently shown that hypoxia-inducible factor 1-alpha (HIF-1α) is excessively produced in breast cancer. During normoxia, Prolyl hydroxylase (PHD) post-translationally alters HIF-1α in the presence of O2, allowing it to interact with the VHL complex comprising elongin-B, elongin-C, CUL2, RBX1, and a ubiquitin-conjugating enzyme (E2)8. This complex, in conjunction with an E1 ubiquitin-activating enzyme, causes the ubiquitylation of HIF-1α and subsequent proteasomal degradation. During hypoxia, PHD cannot change HIF-1α without oxygen, and the protein remains stable and enters the nucleus and interacts with cofactors such as aryl hydrocarbon receptor nuclear translocator (ARNT), CBP/p300, and DNA polymerase II complex to bind to hypoxia-responsive elements (HREs) and stimulates transcription of target genes such as pro-angiogenic factor Vascular Endothelial Growth Factor (VEGF) a primary catalyst of angiogenesis and glucose transporters namely GLUT1 that promotes the transition of tumor cells from oxidative to glycolytic metabolism in low-oxygen environments9,10. As cancer cells undergo a metabolic shift to aerobic glycolysis, they replace pyruvate with lactate, which is then released into the tumor microenvironment (TME). Aerobic glycolysis promotes angiogenesis by generating lactate, which lowers the pH of the surrounding area and increases the expression of VEGF11. In addition, the final products of glycolysis, lactate and pyruvate, affect the expression of VEGF by increasing the levels of HIF-1α12. To summarize, this interplay between HIF-1α,VEGF and GLUT1 has a pivotal function in controlling invasion and metastasis and the dissemination of cancer cells from the original tumor location to distant organs, resulting in the development of additional tumors. Thus, targeting the HIF-1α and its key effectors offers a promising avenue, potentially disrupting the tumor’s hypoxia-driven survival mechanism. However, developing effective inhibitors that selectively target HIF-1α-related pathways remains a challenge.

Phytochemicals and their derivatives have been used as treatments since ancient times. Phytochemicals and their metabolites are present in various parts of the plant, such as the root, leaf, flower, stem, and bark, possessing diverse pharmacological properties13. Moringa oleifera (Mo), commonly referred to as the drumstick tree and also known as the “wonder tree” or “Miracle Tree”, is native to the foothills of the Himalayan ranges in the Indian subcontinent. Mo is extensively utilized in traditional and ayurvedic medicine, with its diverse components believed to possess the ability to prevent 300 diseases14,15. This is attributed to the presence of numerous bioactive compounds such as phenolic acids, vitamins, glucosinolates, carotenoids, polyphenols, alkaloids, flavonoids, saponins, isothiocyanates, tannins, oxalates, and phytates16. Bioactive compounds of Mo, such as Apigenin, could modulate HIF-1α levels via regulation of AKT and NF-kB pathways, and inhibit PI3K/Akt signaling pathway, resulting in simultaneous inhibition of HIF-1α via blocking HIF-1α–Hsp90 interaction both under normoxia and hypoxia17,18. Another isoflavone of Mo namely Daidzein, can downregulate NF-κB in neoplastic cells and inhibit Src/STAT3/HIF-1α pathway activation, inhibit angiogenesis, with other pleiotropic effects on tumor cells19. Another study reported that Isorhamnetin could inhibit the hypoxia response elements, resulting in downregulation of the HIF-1α expression and transcriptionally regulated genes like GLUT1, pyruvate dehydrogenase kinase-1, lactate dehydrogenase-A, and carbonic anhydrase-IX20. In one in silico study, quercetin was found to decrease the binding energy between HIF-1α with FIH and PHD2, thus resulting in an increase in hydroxylation of HIF-1α by FIH and PHD2, which is responsible for proteasomal degradation of HIF-1 and/or transactivation inhibition21. Sulforaphane inhibits neoangiogenesis, cancer cell migration, and adhesion by reducing STAT3, HIF-1, and VEGF expression and decreases the HIF-1α synthesis via reverse regulation of the JNK and ERK signaling pathways22. It was also found to suppress a key regulator of glycolysis pathways, PFKFB4 (6-phosphofructo-2-kinase/fructose 2,6-biphosphatase4), as well as expression of HIF-1, which tightly influences PFKFB4 expression23.

Despite the extensive documentation of Moringa oleifera’s (Mo) therapeutic potential in traditional medicine, most studies to date have focused on crude extracts rather than isolating and characterizing specific bioactive constituents responsible for its anticancer effects. Mo-derived phytochemicals have demonstrated pleiotropic modulation of hypoxia-driven pathways, particularly the HIF-1α/VEGF/GLUT1 axis, which has emerged as a promising target in solid tumors such as breast cancer. However, current literature predominantly explores these compounds in the context of hepatocellular and pancreatic carcinomas, with a notable paucity of studies elucidating their role in suppressing hypoxia-mediated oncogenic signaling in BC. Moreover, many potentially active constituents of Mo remain undiscovered or uncharacterized at the molecular level. This study addresses this critical knowledge gap by employing a systems-level in silico approach to simulate and analyze the interactions between characterized Mo phytocompounds and key hypoxia-regulated targets (HIF-1α, VEGF, and GLUT1). Through molecular docking, pharmacokinetic profiling, long-timescale molecular dynamics simulations, and density functional theory, we aim to identify and prioritize promising candidates for future in vitro and in vivo validation. The overarching goal is to establish a mechanistically sound and computationally validated framework for integrating plant-based HIF-1α inhibitors into modern BC therapeutic strategies.

Methodology

Preparation of extracts and isolation

The seeds of Moringa oleifera were collected from herbarium at AIIMS Bhopal campus, Madhya Pradesh, and it was authenticated by Senior Medical Officer Dr. Danish Javed (Ayurveda), Department of AYUSH (Ayurveda, Yoga and Naturopathy, Unani, Siddha, and Homeopath), AIIMS Bhopal and deposited with the accession number AIIMS/MO/23/001 for reference and traceability. Seeds were subjected to washing using Tween 20, air-dried for 2–3 days, and finely ground into powder using a 100 mm mesh sieve. The extraction was done by cold sequential Maceration Extraction- 500 g of crushed powder was extracted in solvents of varied polarity using petroleum ether, n-butanol, n-hexane, chloroform, ethyl acetate, acetone, ethanol (HPLC Grade), and water for 3–4 days, with continuous shaking at 200 rpm on a magnetic stirrer at 37 °C. Liquid was strained off. Solid residue (marc) pressed (recovered as much as an occluded solution). Strained and expressed liquids were mixed. Liquids were filtered using Whatman filter No. 4, concentrated using a rotary vacuum evaporator, and stored at 4 °C until further use.

Isolation of benzyl isothiocyanate

A Glass Column was used, and packing was done using silica gel of mesh size 60–120 μm. The column was pre-equilibrated with a non-polar solvent, namely n-hexane, to remove impurities. 10–20 g of the alcoholic extracts of seeds were bound to silica and loaded on top of the column. The column was then eluted using a gradient solvent system with varying concentrations of n-hexane + Ethyl acetate, Ethyl acetate + Acetone, and Acetone + Methanol (95:5, 90:10, 80:20, 50:50, 20:80, 10:90 v/v ratios). The eluted fractions were collected in test tubes. Fractions were then spotted on the TLC plate using a capillary tube. The fractions with the same Rf values were mixed and collected. The isolates were again concentrated using a rotary vacuum evaporator and stored at 4 °C.

Characterization

The isolated fractions were analyzed using GC-MS (Gas Chromatography-Mass Spectrometry), equipped with a split/split-less injector and an auto-sampler, connected to an apolar 5-MS capillary column (30 m × 0.25 mm i.d. with a 0.25-µm film thickness) and a mass detector. Mass spectra were acquired in scan mode (70 eV) over a range of 50 to 550 m/z. Each fraction was diluted by adding 20 µl to 2 ml of methanol, and 1 µl of this diluted sample was injected for GC-MS analysis. Mass spectra interpretation was conducted using the National Institute of Standards and Technology (NIST, USA) Mass spectral library, which is a validated reference database with 62,000 known compound spectra. The spectra of the isolates were compared and matched with the known compounds stored in the NIST library based on fragmentation patterns and molecular ion peaks.

Synthesis of aurantiamide acetate from aurantiamide [(-)-2]

figure a

In a flame-dried round-bottom flask, aurantiamide (200 mg, 0.497 mmol, 1 equivalent) was taken, and 5 mL of pyridine was added. To this reaction mixture, acetic anhydride (52 µL, 0.546 mmol, 1.1 equivalent) was added dropwise via a syringe at room temperature. After complete addition, the resulting mixture was stirred for 2 h under the same conditions. After complete consumption of starting materials (as monitored by TLC analysis), the reaction mixture was quenched with 5 mL of 4 N Hydrochloric acid (HCl) and extracted with EtOAc (5 mL × 3). The combined organic layers were washed with brine (5 mL), dried over Na2SO4, and concentrated in vacuo under reduced pressure to get a crude product. The crude materials were purified by flash column chromatography using 40% EtOAc in n-hexane. A pure product as a white solid was obtained (188 mg, 85% yield).

Characterization of aurantiamide acetate

1HNMR and 13CNMR spectra of the compound were recorded using an Advanced III Bruker (500 MHz). The compound was dissolved in deuterated chloroform CDCl3 (Chloroform-d) and poured into a clean NMR tube using a pipette to avoid air bubbles. The measuring conditions for 1HNMR were Frequency – 500 MHz, a pulse angle of 30 degrees, and a pulse repetition time of 1–3 s. 13CNMR spectra of the compound were recorded at 125 MHz, a pulse angle of 30–45 degrees, and a pulse repetition time of 4–7 s. Later, HRMS (High-Resolution Mass Spectrometry) was recorded using micrOTOF-Q- Q II (MaXis Impact Bruker). The compound was dissolved in methanol at a concentration of 5 µg/mL. The solution was filtered using a 0.22 μm syringe filter to remove any particulate matter. The spectra were recorded using the Electrospray Ionization (ESI) technique. The observed m/z values were compared with theoretical values to confirm the identity of the compound.

Target proteins data set

X-ray crystallographic structure of human HIF-1α, VEGF, and GLUT1 were retrieved from the protein data bank (https://www.rcsb.org/) PDB ID: 1LQB (resolution: 2.00 Å; R-value free: 0.277)24, 4ASD (resolution: 2.03 Å; R-value free: 0.230)25, 5EQI (resolution: 3.00 Å; R-value free: 0.293)26 in .pdb format. Protein preparation was done in AutoDock Tools v4.2.6; co-crystallized water molecules, heteroatoms, and nonstandard residues were deleted. Polar Hydrogens were added to provide protonation states at physiological pH. Partial charges were assigned by adding Kollman charges. Furthermore, incomplete side chains were also repaired27.

Ligand design

Canonical smiles (Supporting Information 2) and corresponding 3D chemical structures of the Aurantiamide acetate and Benzyl isothiocyanate were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/), in SDF format. The SDF files of the structure were changed to PDB format using the open Babel tool, rendering them suitable for docking studies. (https://www.cheminfo.org/Chemistry/Cheminformatics/FormatConverter/index.html). Subsequently, Gasteiger charges were added, and torsion requirements for proper binding were set using Autodock tools v4.2.628.

Molecular docking studies

The software AutoDock v4.2.6 was utilized to convert both receptor and ligand PDB structures into pdbqt format, which comprises atomic charges, atom-type definitions for protein, and topological information such as rotatable bonds for ligands28. Later, a grid was centered to adequately cover the binding site of the structures. Subsequently, receptor and ligand molecules were subjected to molecular docking by using the AutoDock Vina program to generate 9 conformations of the ligand in complex with the receptor29. The final ranking of these complexes was based on binding energy. The grid dimension for molecular docking was set with grid centers 44.812, 54.179, 66.219 Å; -21.981, -1.141, -3.791 Å; -41.587, 8.363, 10.095 Å (x, y, z) for HIF-1α, VEGF, and GLUT1, respectively. After the molecular docking was completed, the docked poses were analyzed by Discovery Studio 09 software30.

Molecular dynamics simulations

Molecular dynamics simulations were conducted using the Desmond simulation package LLC (Schrodinger). The OPLS-3e (optimized potential for liquid simulations- Schrodinger) force field’s protein preparation wizard was employed for the preparation, optimization, and minimization of the docked locations over time, assessing the stability of the MD simulations. The initial system setup encompassed the protein-ligand complexes within a TIP3P water model configured in an orthorhombic shape, constructed using a system builder with the dimensions of 20 Å X 20 Å X 20 Å. To achieve neutralization, an appropriate number of metal ions were added, maintaining a salt concentration of 0.15 M. Each complex’s-built system underwent a standard equilibration protocol. Simulations were initiated at 300 K and 1 bar pressure, pH 7.4, for a duration of 500 ns. The NPT ensemble was used for simulation. Each simulation was performed using periodic boundary conditions (PBC) in all three directions. The long-range electrostatic interactions were managed with the Particle Mesh Ewald (PME) summation. The Berendsen thermostat was used for maintaining the temperature. relaxation constant of 1.0 ps. The Parrinello-Rahman barostat was used for regulating the pressure at 1 bar with a 5 ps coupling constant and compressibility of 4.5 * 10− 5 bar. Trajectories were recorded over 500 ns to analyze the dynamic nature of the interactions and component stability. The results, including root mean square fluctuation (RMSF), root mean square deviation (RMSD), and hydrogen bonds, solvent accessible surface area (SASA), and radius of gyration (rGyr), were analyzed using the Desmond simulation interaction analysis tool31.

DFT studies were performed with Gaussian 09 software. All calculations were performed using the combination of the hybrid density functional B3LYP32 with the 6-311 + G basis set. The approach is furthermore referred to as B3LYP/6-311 + G. The charge analysis was performed using the natural bond orbital (NBO) method as implemented in the Gaussian 09 program33, using the B3LYP/6-311 + G approach with the implicit water effects. Frontier molecular orbitals (FMOs) and molecular electrostatic potential (MEP) were computed at the B3LYP/6-311 + G level with the implicit water effects as well. We considered the calculated structures of the AA and BITC, structural parameters for the lowest-lying isomer, its NBO charges, FMOs, and MEP. Furthermore, we used the values of the energies of the lowest-lying isomer, HOMO and LUMO, to compute its global reactivity parameters (GRP). Equations 1, 2, 3 were used to calculate the values of the band gap (B), ionization potential (IP), and electron affinity (EA):

$$\:B=\:{E}_{HOMO}-{E}_{LUMO}$$
(1)
$$\:IP=\:{-E}_{HOMO}$$
(2)
$$\:EA=\:{-E}_{LUMO}$$
(3)

For electronegativity χ, and global hardness η, Softness (S) values, we used Eqs. 4, 5 & 6:

$$\:\chi\:=\:-\mu\:$$
(4)
$$\:\eta\:=\:\:\frac{IP-EA}{2}$$
(5)
$$\:S=\:\:\frac{1}{\eta\:}$$
(6)

Also, the global electrophilicity index ω value was calculated by Eq. 7:

$$\:\omega\:=\:\frac{{\mu\:}^{2}}{2\eta\:}$$
(7)

And µ, the chemical potential of the system, was calculated by Eq. 8:

$$\:\mu\:=\:-\frac{IP+EA}{2}$$
(8)

ADME prediction by computational analysis

The pharmacokinetic properties, which include absorption, distribution, metabolism, excretion, and toxicity (ADMET) were computed using the SwissADME (http://www.swissadme.ch/) and ADMETSaR (http://lmmd.ecust.edu.cn/admetsar2/) web server. The main characteristics determinants of PK properties were analyzed namely the polar surface area PSA_2D (determines fractional absorption), atom-based LogP (AlogP98) for lipophilicity levels, membrane permeability [indicated by colon cancer cell line (Caco-2)], skin permeability, intestinal absorption, and P-glycoprotein substrate or inhibitor, distribution through blood-brain barrier (BBB), CYP models for substrate or inhibition (CYP3A4, CYP2D6, CYP1A2, CYP2C9, CYP2C19, CYP2D6, and CYP2C19), the total drug clearance via renal OCT2 inhibitor. Drug toxicity was predicted based on AMES toxicity, hERG inhibition, hepatotoxicity, and skin sensitization. These parameters were checked for compliance with their standard ranges.

Drug-likeness

The drug-likeness of the compounds was determined using SwissADME (http://www.swissadme.ch/). The canonical SMILES of the selected compounds were acquired from the PubChem database. The compounds underwent Lipinski, Ghose, Veber, Egan, and Muegge screening.

Bioavailability radar plot

The bioavailability of the compounds was analyzed employing an online ADMETlab 3.0 (https://admetlab3.scbdd.com/) web server. The canonical SMILES of the chosen compounds were obtained from the PubChem Database.

Results

Isolation and characterization of the benzyl isothiocyanate from Mo

We were able to isolate pure fractions from the ethyl acetate-ethanolic fraction of the seeds using column chromatography. The GC-MS analysis was performed to identify the compounds present in the fractions. The mass spectrum obtained for the Fraction is shown in Fig. 1. The molecular ion peak was observed at m/z 149. In addition to the molecular ion peak, the mass spectrum displayed several prominent fragment ions, m/z 91, m/z 105, and m/z 65. The calculated mass was 149 g/mol. To confirm the identity of the compound, the spectrum was compared against the NIST Mass Spectral Library using the NIST MS Search Program. Based on the mass spectrum analysis, calculated molecular weight, and comparison with known spectral data, the Fraction was identified as isothiocyanatomethylbenzene (Benzyl isothiocyanate). This identification is supported by the molecular ion peak and the characteristic fragmentation pattern observed in the GC-MS analysis (Table 1).

Table 1 Assigned identification of compounds from GC-MS analysis.
Full size table
Fig. 1
figure 1

GC-MS analysis of the extract showing the retention time (min) of the compounds in the X axis and percentage (%) of peak area in the Y axis. Benzyl Isothiocyanate (BITC).

Full size image

Characterization of aurantiamide acetate

Aurantiamide acetate was characterized by NMR and HRMS: (S)-2-((S)-2-benzamido-3-phenylpropanamido)-3-phenylpropyl acetate [(-)-2], using 1H NMR and 13C NMR spectra and High-Resolution Mass spectrometry (HRMS) (Fig. 2).

figure c

1H NMR (500 MHz, CDCl3) δ 7.74–7.66 (m, 2 H), 7.51 (t, J = 7.4 Hz, 1 H), 7.42 (t, J = 7.6 Hz, 2 H), 7.28–7.21 (m, 5 H), 7.17–7.08 (m, 3 H), 7.08–7.01 (m, 2 H), 6.80 (d, J = 7.9 Hz, 1 H), 6.08 (d, J = 10.6 Hz, 1 H), 4.78 (q, J = 7.5, 7.0 Hz, 1 H), 4.39–4.29 (m, 1 H), 3.92 (dd, J = 11.4, 5.0 Hz, 1 H), 3.82 (dd, J = 11.1, 4.3 Hz, 1 H), 3.21 (dd, J = 13.6, 6.0 Hz, 1 H), 3.06 (dd, J = 13.7, 8.3 Hz, 1 H), 2.80–2.68 (m, 2 H), 2.01 (s, 3 H).

13C NMR (125 MHz, CDCl3) δ 170.8, 170.3, 167.1, 136.7, 136.6, 133.7, 131.9, 129.3, 129.1, 128.7, 128.6, 128.6, 127.1, 127.1, 126.7, 64.6, 55.0, 49.5, 38.5, 37.5, 20.8.

Fig. 2
figure 2

(A) 1H NMR (500 MHz, CDCl3) (B) 13C NMR (125 MHz, CDCl3) of Aurantiamide acetate.

Full size image

The high-resolution mass spectrometry graph is shown in Fig. 3. The topmost graph shows the Base Peak Chromatogram (BPC), which has a prominent peak indicating the elution of Aurantiamide Acetate at around 0.8 min. The main ion of interest is at m/z 467.1961, supporting the identification of Aurantiamide Acetate.

Fig. 3
figure 3

High-resolution mass spectrometry analysis showing the retention time (min) of the compounds in the X axis and percentage (%) of peak area in the Y axis.

Full size image

Molecular docking studies of selected bioactive compounds AA and BITC showed good binding with TNBC targets HIF-1α, VEGF, and GLUT1

Table 2 shows the binding energy and hydrogen bond interactions of the two bioactive compounds of Mo namely Aurantiamide acetate (AA) and Benzyl isothiocyanate (BITC) simulated by AutoDock Vina ranging from − 6.8, -4.2 kcal/mol for HIF-1α, -7.9, -4.7 kcal/mol for VEGF, and − 7.3, -5.2 kcal/mol for GLUT1 respectively.

Table 2 Binding energies of the ligands with protein:
Full size table

In addition to conventional hydrogen bonding other stabilizing interactions such as van der Waal forces, pi bonds, and alkyl bonds were also seen between the ligands and amino acid residues present in the active site of the proteins (Figs. 4, 5 and 6). Based on the binding energies and stabilizing interactions, the compounds (ligands) were found to have strong affinity for the target proteins.

Fig. 4
figure 4

3D and 2D Docked pose and interacting residues of (A) AA and HIF-1α (B) BITC and HIF-1α.

Full size image
Fig. 5
figure 5

3D and 2D Docked pose and interacting residues of (A) AA and VEGF (B) BITC and VEGF.

Full size image
Fig. 6
figure 6

3D and 2D docked pose and interacting residues of (A) AA and GLUT1 (B) BITC and GLUT1.

Full size image

The absorption, distribution, metabolism, excretion, and toxicity profiles of the bioactive compounds showed good pharmacokinetic properties in silico

The ADMET properties of AA and BITC are presented in Table 3. The PSA_2D of AA and BITC were less than 100, indicating good oral absorption. Both compounds were predicted as having ideal lipophilicity (AlogP98 ≤ 5). When the Papp coefficient is > 8 × 10− 6, the predicted value is > 0.90; thus, the compound is said to have high Caco-2 permeability and is easy to absorb. BITC and AA were predicted to have high Caco-2 permeability, suggesting high intestinal absorption. If log Kp > − 2.5, the compound is considered to have relatively low skin permeability. All compounds were found to have low skin permeability. AA was a substrate as well as an inhibitor for P-glycoprotein, whereas BITC was a non-substrate and non-inhibitor. For blood-brain barrier membrane permeability, only BITC was found to be permeable. Neither compound was a substrate for CYP2D6; AA was a substrate for CYP3A4 and was a CYP2C19 and CYP2C9 inhibitor. The prediction results show that none of the compounds were OCT2 substrates. BITC may be hepatotoxic, AA was a hERG inhibitor, and may cause skin sensitization. Overall, the compounds were found to have good pharmacokinetic properties.

Table 3 Predicted pharmacokinetics properties of the ligands:
Full size table

In silico drug-likeness filters identified the bioactive compounds to be drug-like

The drug-likeness of the compounds, was assessed using filter variants such as the Lipinski rule of (≤ 5 hydrogen bond donors, ≤ 10 hydrogen bond acceptors, A molecular weight (MW) ≤ 500 Da, An octanol-water partition coefficient (log P) ≤ 5), Ghose filter (Atomic mass between 160 and 480 Da, Log P between − 0.4 and 5.6, Molar refractivity from 40 to 130, Total number of atoms from 20 to 70), Veber rule (Rotatable bonds ≤ 10, Polar surface area (PSA) ≤ 140 Å2), Egan rule, and Muegge rule. Aurantiamide acetate had one violation, whereas Benzyl Isothiocyanate had two violations, largely suggesting a strong drug-likeness Table 4.

Table 4 Drug-Likeness of the compounds.
Full size table

The bioactive compounds demonstrated good oral bioavailability in silico

The bioavailability of the drug was predicted in silico through bioavailability radar plots. Each axis represents a distinct attribute related to the compound’s bioavailability, as listed in Table 5 and illustrated in Fig. 7.

Fig. 7
figure 7

Bioavailability Radars Plots of Ligands (A) AA and (B) BITC. The green area represents the lower limit for each attribute, suggesting the minimum value, and the blue area represents the upper limit of each property a compound should have to be considered bioavailable. The yellow line and points represent the actual compound properties.

Full size image
Table 5 Bioavailability attributes and their acceptable reference ranges (RR) for oral drug-likeness.
Full size table

Aurantiamide acetate (AA) had two violations, namely 21 rigid bonds (≤ 10) and 13 rotatable bonds (≤ 10). Benzyl isothiocyanate (BITC) had one violation. Overall, both compounds showed minimal violations and were suggestive of good bioavailability.

Molecular dynamic simulations revealed ideal stability of the protein-ligand complex

The conformational stability and dynamic interaction of the protein HIF-1α, VEGF, and GLUT1, and the ligand AA and BITC over a 500-nanosecond (ns) molecular dynamics simulation are presented in Figs. 8, 9, 10, 11, 12 and 13 (A-F) in the form of Root mean square deviation (RMSD), RMSF (Root Mean Square Fluctuation), SASA (Solvent-Accessible Surface Area), and rGyr (Radius of Gyration) of protein-ligand.

The RMSD trajectory for the HIF-1α–AA complex stabilized after 100 ns, maintaining consistent protein RMSD values between 4.3 and 4.8 Å, indicating minimal backbone fluctuations. Whereas Ligand RMSD showed a monotonic increase up to 20 Å by the end of the run. In contrast, the BITC-HIF-1α complex had initial ligand mobility (< 100 ns), with BITC RMSD stabilizing around 3.5–5 Å. The protein RMSD, however, stayed comparatively constant at 2.5–3.5 Å with small oscillations. Per-residue flexibility and areas displaying local conformational alterations were assessed using RMSF. HIF-1α-AA showed prominent peaks at important binding site residues (avg. RMSF: 3–10 Å). In line with anticipated flexible regions, terminal loop sections displayed a slightly increased RMSF. The binding site residues for HIF-1α-BITC showed very low fluctuations, remaining flat for all 10 atoms. SASA analysis helped determine the extent of protein surface exposure to solvent and revealed insights into protein folding dynamics. Average SASA for HIF-1α-AA was ~ 140–160 Å2, with ~ 100–200 Å2 repeated peaks and for HIF-1α-BITC, it had transient spikes reaching ~ 300 Å2. rGyr analysis was used to assess the compactness of protein structure over the 500 ns simulation. The rGyr for HIF-1α-AA remained stable with mild fluctuations ~ 3.9–4.7 Å, whereas for HIF-1α-BITC it was stable and compact at ~ 2.5–3 Å.

Fig. 8
figure 8

HIF-1α-AA (A) RMSD plot. (B) RMSF profile. (C) Normalized stacked bar charts of binding site residues interacting with ligand. (D) Schematic representations of ligand interactions with protein. (E) Radius of gyration (rGyr) plot. (F) SASA profile.

Full size image
Fig. 9
figure 9

HIF-1α-BITC (A) RMSD plot. (B) RMSF profile. (C) Normalized stacked bar charts of binding site residues interacting with ligand. (D) Schematic representations of ligand interactions with protein. (E) Radius of gyration (rGyr) plot. (F) SASA profile.

Full size image

The RMSD of VEGF bound to AA remained consistent between 2 and 3.5 Å for the protein and 1.5–3.2 Å for the ligand. The VEGF-BITC complex also demonstrated low variations, with protein RMSD attaining 3.5 Å, while ligand RMSD, although it rises up to 5 Å but remains consistent. VEGF-AA exhibited low fluctuations in RMSF values. Active site residues exhibited stability within the range of approximately 1.0 to 1.7 Å. VEGF-BITC showed higher RMSF values around 3.5–4.5 Å, suggesting greater flexibility in binding site residues. For VEGF-AA, the rGyr stabilized after 50 ns with minor fluctuations approximately ~ 4.2–5.0 Å, for VEGF-BITC, rGyr was compact and stable (~ 2.4–3.0 Å). The SASA for VEGF-AA remained steady over simulation, ~ 80–160 Å2, indicating a properly folded structure. However, for VEGF-BITC, SASA drops sharply and stabilizes at ~ 20–60 Å2.

Fig. 10
figure 10

VEGF-AA (A) RMSD plot. (B) RMSF profile. (C) Normalized stacked bar charts of binding site residues interacting with ligand. (D) Schematic representations of ligand interactions with protein. (E) Radius of gyration. (rGyr) plot. (F) SASA profile.

Full size image
Fig. 11
figure 11

VEGF-BITC. (A) RMSD plot. (B) RMSF profile. (C) Normalized stacked bar charts of binding site residues interacting with ligand. (D) Schematic representations of ligand interactions with protein. (E) Radius of gyration (rGyr) plot. (F) SASA profile.

Full size image

The GLUT1–AA complex had fairly stable interactions, with protein RMSD typically ranging from 2.5 to 4.0 Å and ligand RMSD fluctuating between 4.5 and 7 Å. GLUT1 complexed with BITC exhibited minor fluctuations in ligand RMSD, mostly between 3.0 and 5 Å, indicating modest structural flexibility. Nonetheless, the protein RMSD remained steady at approximately 2.0–4.0 Å. The RMSF values for GLUT1-AA concerning transmembrane helices were notably low, approximately 1.2–3 Å, whereas a few regions showed higher fluctuations. The terminal sections and loop portions of GLUT1-BITC exhibited minimal fluctuations, while the core residues demonstrated stability, but the ligand itself showed structural rigidity. GLUT1-AA SASA was around ~ 50–200 Å2 with oscillations indicating an augmented spatial dispersion of structural mass distribution. GLUT1-BITC SASA remained steady at ~ 50–300 Å2, with occasional spikes. GLUT1-AA exhibited consistent compactness with a radius of gyration ranging from ~ 4.5 to ~ 4.9 Å. GLUT1-BITC exhibited a stable rGyr of 2.4–2.9 Å.

Fig. 12
figure 12

GLUT1-AA. (A) RMSD plot. (B) RMSF profile. (C) Normalized stacked bar charts of binding site residues interacting with ligand. (D) Schematic representations of ligand interactions with protein. (E) Radius of gyration (rGyr) plot. (F) SASA profile.

Full size image
Fig. 13
figure 13

GLUT1-BITC. (A) RMSD plot. (B) RMSF profile. (C) Normalized stacked bar charts of binding site residues interacting with ligand. (D) Schematic representations of ligand interactions with protein. (E) Radius of gyration (rGyr) plot. (F) SASA profile.

Full size image

DFT studies

The optimization of the compound BITC revealed that the isothiocyanate group (-N = C = S) does not remain planar to the plane of the benzene ring upon geometry optimization (Fig. 14). It rather stays in the other plane, where the connecting CH2 group acts as a pivot for the molecule.

Fig. 14
figure 14

Optimized structure of BITC in B3LYP/6-311 + G level of theory.

Full size image

The FMO analysis of BITC indicates the HOMO electron density of BITC mainly resides in the isothiocyanate group (-N = C = S) group and the connecting -CH2 group. But the LUMO electron density of the molecule is spread throughout the molecule (Fig. 15). The band gap of the molecule is also found to be ∆E = 5.61466 eV. The electrostatic potential analysis of BITC indicates that the negative potential is localized over the isothiocyanate group. This indicates the isothiocyanate possesses the affinity to accept as an electron donor, and it possesses the affinity to form hydrogen bonds in the protein structure.

Fig. 15
figure 15

The frontier molecular orbitals HOMO and LUMO and electrostatic potential map at the ground state at B3LYP/6-311 + G.

Full size image

The NBO analysis of the isothiocyanate bond indicates that the -N = C part of the -N = C = S bond consists of three parts. The σ bond is attributed to the overlap of the sp-sp-hybridized orbital of N and the sp-sp-hybridized orbital of C. The π bonds are attributed to the overlap of the remaining p orbitals of N and C atoms (Fig. 16). Whereas the C = S part of the -N = C = S bond consists of two parts. The σ bond is attributed to the overlap of the sp hybridized orbital of C and the sp3 hybridized orbital of S, and the π bond is attributed to the overlap of the remaining sp3 orbital of C and the p orbital of the S atom. The remaining orbitals of the S atom can readily donate the electron density and act as an acceptor in the formation of electrostatic bonds and hydrogen bonds.

Fig. 16
figure 16

The NBO analysis of the BITC isothiocyanate bond at B3LYP/6-311 + G level of theory.

Full size image

The optimization of the compound AA indicated that the three benzene groups present in the molecule upon geometry optimization reside in three different planes, and they are connected through the aliphatic chains (Fig. 17).

Fig. 17
figure 17

Optimized structure of AA in B3LYP/6-311 + G level of theory.

Full size image

The FMO analysis of AA indicates that the HOMO electron density of AA mainly resides in the aliphatic chain of the molecule (Fig. 18). But the LUMO electron density of the molecule is localized over the benzene ring, indicating an electron transfer from aliphatic carbons to aromatic carbons. The band gap of the molecule is also found to be ∆E = 5.142413 eV. The electrostatic potential analysis of AA indicates that the negative potential is localized over the carbonyl oxygen atoms. This indicates that the carbonyl oxygen atoms possess the affinity to accept electrons as donors, and they possess the affinity to form hydrogen bonds within the protein structure.

Fig. 18
figure 18

The frontier molecular orbitals, HOMO and LUMO, at the ground state at B3LYP/6-311 + G.

Full size image

The NBO analysis of the carbonyl oxygen bond indicates that the σ bonds are attributed to the overlap of the sp3 hybridized orbital of C and the sp hybridized orbital of O (Fig. 19). The π bonds are attributed to the overlap of the p orbitals of C and O atoms. Due to the electron density localization and presence of free p orbital over O atoms, they can readily donate the electron density and act as an acceptor in the formation of electrostatic bonds and hydrogen bonds.

Fig. 19
figure 19

The NBO analysis of the AA carbonyl bond at B3LYP/6-311 + G level of theory.

Full size image

The analysis of Global Reactivity Parameters (GRP) values presented in Table 6 indicates that AA possesses a lower band gap (0.54 eV) than BITC (0.89 eV), implying that AA is more chemically reactive and may exhibit greater efficiency in charge transfer processes within the binding pocket. The elevated electrophilicity index (ω = 59.22 eV) for AA indicates an enhanced ability to take electrons, signifying pronounced reactivity towards nucleophilic residues in the active site. Increased softness (S) and decreased hardness (η) for AA suggest enhanced compatibility with the electronic milieu of the target proteins, hence promoting complex formation. BITC, exhibiting lower electrophilicity (0.54 eV), demonstrates superior chemical stability, evidenced by a larger band gap (0.89 eV) and hardness (0.44 eV), suggesting it may enhance systemic half-life while perhaps exhibiting reduced reactivity at the target location.

Table 6 Calculated GRPs for the compounds (eV):
Full size table

Discussion

Molecular docking studies are a valuable in silico technique that can accurately predict prospective drug candidates for therapeutic development against specific target proteins41. Binding energy provides insights into the strength with which a compound (ligand) binds to the pocket of a target protein (receptor), and a compound with more negative binding energy is desirable42. The molecular docking studies, conducted using AutoDock Vina, evaluated the binding interactions of Aurantiamide acetate (AA) and Benzyl isothiocyanate (BITC) with three target proteins: HIF-1α, VEGF, and GLUT1. The binding energies and key residues involved in hydrogen bonding were as follows: HIF-1α: AA (-6.8 kcal/mol, GLU160, ARG200) and BITC (-4.2 kcal/mol). VEGF: AA (-7.9 kcal/mol, ASP1046) and BITC (-4.7 kcal/mol, ASP 1056). GLUT1: AA (-7.3 kcal/mol) and BITC (-5.2 kcal/mol). Besides conventional hydrogen bonding, other stabilizing interactions, including van der Waals forces, pi bonds, and alkyl bonds, were also observed between the ligands and the amino acid residues in the protein’s active site. AutoDock Vina’s scoring algorithm predicted a stronger binding affinity of these bioactive compounds to HIF-1α, VEGF, and GLUT1 binding pockets, indicating their potential as promising leads targeting breast cancer.

To further confirm the stability of AA and BITC compounds in the binding pocket of HIF-1α, VEGF, and GLUT1, we extracted Molecular dynamics simulations conformations at different intervals and visualized the interactions between protein and ligands. The molecular dynamics simulations provide an assessment of the binding stability of the protein-ligand complex on the atomic level43. The superposition of the ligand conformation at the beginning and end of the simulation gives us information about the stability of the complex44. Target proteins (HIF-1α, VEGF, GLUT1) exhibited significant structural and dynamic changes when complexed with AA and BITC in 500 ns molecular dynamics simulations. Both AA and BITC-bound complexes showed great conformational stability. Interestingly, the HIF-1α–AA complex reached equilibrium quickly and maintained structural integrity. While protein backbones remained below 3.4 Å, BITC complexes, especially HIF-1α–BITC, showed higher ligand RMSD values but still stable binding pose after equilibration45. RMSF analysis revealed that BITC had reduced residue mobility at binding interfaces (1.0-1.8 Å), indicating a rigid protein-ligand combination. In contrast, VEGF complexes showed higher fluctuations in loop and periphery domains, indicating reduced local stabilization46. SASA and Rg analyses revealed that AA-bound complexes (HIF-1α, VEGF, GLUT1) consistently maintained lower fluctuations, reflecting compact and stable protein conformations. In contrast, BITC-bound systems displayed greater SASA variability and mild rGyr perturbations, indicating higher solvent exposure and conformational flexibility. Overall, AA conferred greater structural stability across all targets, whereas BITC binding was associated with comparatively weaker stabilization and dynamic behavior47,48. These data suggest AA interacts with hypoxia-associated protein targets in a more stable, solvent-free, and conformationally restricted manner. AA outperforms BITC in structural performance, supported by docking scores and DFT-derived electronic descriptors, indicating its potential as a strong HIF-1α/VEGF/GLUT1 modulator in physiological circumstances.

Evaluating absorption, distribution, metabolism, excretion, and toxicology (ADMET) is crucial in drug discovery. Laboratory assessments of pharmacokinetic properties are costly and time-consuming, making in silico evaluations a more efficient and feasible alternative. This approach helps prevent drug failures during clinical trials49. The in silico ADMET profiling of AA and BITC revealed promising pharmacokinetic properties. Both compounds exhibited good oral absorption (PSA < 100 Å2) and ideal lipophilicity (AlogP98 ≤ 5). BITC showed high Caco-2 permeability (a model to study intestinal absorption), and AA also indicated good intestinal absorption50. Both compounds had low skin permeability (log Kp > -2.5). AA was both a substrate and inhibitor for P-glycoprotein, potentially affecting drug distribution51. Neither compound was a substrate for CYP2D6, suggesting a minimal impact on the metabolism of other drugs52. They were not inhibitors of renal OCT2, indicating a low risk of renal toxicity53. Neither compound caused genetic mutations (AMES test)54 or skin sensitization. However, AA was found to be a hERG inhibitor, suggesting potential cardiac toxicity55, and BITC indicated potential hepatotoxicity. Consequently, the compounds exhibit satisfactory ADMET properties, with few exceptions, indicating their potential as therapeutic agents. The drug-likeness of AA and BITC was assessed using various filters (Lipinski, Ghose, Veber, Egan, and Muegge rules), with AA showing only one violation and BITC showing 2 violations, indicating their potential as orally bioavailable drugs. The bioavailability radar plots showed that AA exceeded only two limits for nRig and nRot, whereas BITC met all criteria, with only one violation which further confirmed that these compounds fit within the desired pharmacokinetic parameters.

The electronic characteristics of AA and BITC were viewed quantum mechanically by DFT analysis. The Frontier Molecular Orbitals (FMO) showed that BITC’s HOMO is mostly on the isothiocyanate (-N = C = S) group, indicating its function in hydrogen bonding and nucleophilic interactions. Delocalizing the LUMO in the aromatic system allows for π-π interactions and dispersion forces with hydrophobic residues56. The HOMO was localized on the aliphatic chains, and the LUMO on the aromatic rings in AA, indicating intramolecular charge transfer from side chains to the core scaffold. Electrostatic potential mapping (ESP) confirmed nucleophilic areas on AA’s carbonyl oxygen atoms and BITC’s isothiocyanate group, which are ideal for the protein’s active site’s electrostatic potential and hydrogen bonding57. Natural Bond Orbital (NBO) research has shown substantial delocalization in π systems and lone pair donation in sulfur (BITC) and oxygen (AA), thereby promoting binding interactions. Critical analysis of Aurantiamide Acetate and Benzyl Isothiocyanate GRPs was also performed. GRPs indicate molecular stability, reactivity, and electrophilicity. Aurantiamide Acetate (AA) had significant chemical reactivity and polarizability due to its low bandgap (0.54 eV), high electrophilicity index (ω = 59.22), and high softness (S = 3.64). This indicates high electron acceptance and adaptation to the protein binding pocket’s electrical conditions. AA stabilizes negative charge distribution with polar residues due to its high chemical potential and electron affinity. Benzyl Isothiocyanate (BITC) has a higher bandgap (0.89 eV) and chemical hardness (η = 0.44), representing decreased reactivity but increased kinetic stability. Lower electrophilicity (ω = 0.54) indicates decreased charge transfer interactions, supporting its greater RMSD in MD simulations58. Due to its stability and cellular bioactivity, BITC may remain biologically relevant. DFT and GRP data confirm docking and MDS studies by showing AA’s greater electronic flexibility and reactivity, while BITC’s more stable but less dynamic interaction profile may contribute to biological consequences.

Overall, this study used an integrated strategy to identify and characterize natural inhibitors of the hypoxia-driven HIF-1α/VEGF/GLUT1 axis in breast cancer, focusing on bioactive compounds from Moringa oleifera. Several studies have shown the anticancer effects of Moringa crude extracts, but this is the first study to the best of our knowledge which systematically evaluated the compounds Aurantiamide Acetate (AA) and Benzyl Isothiocyanate (BITC) using a multi-parametric in silico pipeline, including molecular docking, ADMET profiling, long-timescale (500 ns) molecular dynamics simulations, and DFT studies, explicitly focusing on three hypoxia-associated proteins (HIF-1α, VEGF, and GLUT1), targeting multiple hallmarks of tumor aggressiveness, such as angiogenesis, metabolic adaptation, and invasion.

But despite the robustness of the computational framework, the study has a few limitations, such as in silico predictions have not been tested in vitro or in vivo. Lack of benchmarking with known HIF inhibitors hampers contextualization of compound performance. The tumor microenvironment, protein allosterism, and metabolic interactions in vivo are too complex for simulations. However, these limitations do not invalidate the findings but rather highlight important translational development directions and experimental validation.

Conclusion

This study utilized an extensive in silico framework to assess the therapeutic potential of two bioactive chemicals derived from Moringa oleifera—Aurantiamide Acetate (AA) and Benzyl Isothiocyanate (BITC), focusing on hypoxia-regulated molecular pathways critical to breast cancer progression. Molecular docking indicated that both AA and BITC displayed significant binding affinities for HIF-1α, VEGF, and GLUT1, with AA exhibiting notably robust and stable interactions. ADMET profiling validated their drug-likeness, bioavailability, and advantageous safety parameters. Extended molecular dynamics simulations (500 ns) confirmed the conformational stability and compactness of the protein-ligand complexes, while DFT-derived global reactivity parameters further validated the compound’s exceptional electronic adaptability and electrophilicity, underscoring its potential as a hypoxia-targeted therapeutic candidate.

Collectively, our findings designate AA and BITC as a potential natural inhibitor of the HIF-1α/VEGF/GLUT1 axis, providing the twin benefits of structural integrity and pharmacological safety. To translate these in silico discoveries into therapeutic significance, in vitro assessment of cytotoxicity, apoptotic induction, and hypoxia gene regulation in breast cancer cell lines is recommended, along with in vivo xenograft studies to evaluate tumor regression, pharmacokinetics, and safety profile in hypoxic tumor microenvironments. Investigation of formulation techniques, such as nanoparticles or liposomes, to improve the bioavailability of compounds should also be explored. Augmentation of the compound library derived from Moringa oleifera and structurally analogous flora to discover novel inhibitors of the HIF pathway is encouraged. This study ultimately seeks to connect computational predictions with experimental validation, thereby facilitating the advancement of phytochemical-based therapies for aggressive subtypes of breast cancer.