Unlocking the therapeutic potential of Origanum majorana through GC-MS and computational analysis to identify EGFR inhibitors for breast cancer.

Introduction
In the relentless battle against disease, cancer remains one of humanity’s most formidable adversaries, with breast cancer casting a particularly long shadow. It is a diagnosis that echoes through millions of lives each year, holding the grim distinction of being the world’s most prevalent cancer and a leading cause of mortality among women1. This global burden continues to grow, fueled by complex modern lifestyles and unequal access to care, turning the search for effective treatments into one of modern medicine’s most urgent quests2,3. The enemy, however, is cunning, often hijacking the body’s own cellular machinery to fuel its destructive campaign. A key conspirator in this insurgency is a protein known as the EGFR. When overexpressed, it acts like a rogue commander, relentlessly signalling for cell division and tumor growth, especially in aggressive forms of the disease like triple-negative breast cancer (TNBC), a subtype with a critical lack of targeted treatment options4,5. Blocking this commander has become a critical strategic goal6,7.
Yet, sometimes, the most promising solutions are not forged in a laboratory furnace but whispered by nature itself. For centuries, traditional healers have turned to the botanical world, and today, science is lending a rigorous ear. Among these natural apothecaries is Origanum majorana L. (sweet marjoram) a humble herb native to the sun-drenched lands of the Mediterranean and Middle East, long revered not just for its aroma but for its healing virtues8,9. Could this plant, thriving under the Saudi Arabian sun, hold within its essence a key to confronting a modern plague? Early investigations suggest it might; its essential oil has shown a remarkable ability to disrupt cancer cells10–12. The source of this power is believed to lie in its main chemical constituents: a fleet of volatile monoterpenes like sabinene, α-terpinene, and limonene molecules known for their potent biological activities, including the ability to induce cancer cell death13,14. Therefore, the first crucial step in unlocking the anticancer potential of Saudi Arabian O. majorana was a detailed chemical inventory. Gas Chromatography-Mass Spectrometry (GC-MS) analysis was employed to definitively characterize the complete phytochemical composition of the essential oil15. This analytical technique provides a precise qualitative and quantitative profile, identifying each constituent compound based on its unique mass-to-charge ratio and retention time. The resulting chromatogram served as the definitive source from which the major and biologically relevant compounds were selected for further investigation. This targeted approach moves beyond a generic study of “marjoram oil” to a precise analysis of its specific, dominant chemical entities. From the comprehensive list of compounds identified by GC-MS, a selection of the most abundant and pharmacologically promising monoterpenes and monoterpenoids was made.
This narrative leads us from the field to the computer terminal. To uncover this hidden potential, we turned to the digital realm of computational science. In modern drug discovery, computational screening methods such as molecular docking and molecular dynamics (MD) simulations are now integral, offering an efficient strategy to elucidate the structural fundamentals of ligand-target recognition16,17. Analogous to designing a key for a complementary lock, molecular docking predicts the binding interactions between small molecules and proteins and provides estimates of binding affinity18–22.
The rationale for focusing on EGFR in this study is based on its well-established and critical role in the molecular pathogenesis of breast cancer17,23–25. The receptor’s aberrant activation and overexpression are strongly linked to aggressive tumor progression, increased metastasis, and poor clinical prognosis, particularly in the TNBC subtype where effective therapeutic options are limited6,7. EGFR acts as a pivotal regulator of oncogenic signaling pathways, driving dysregulated cell proliferation, survival, and invasive potential4. Recently, computational work has applied structure-based design to discover EGFR inhibitors. For example, Mishra et al. (2025)23 performed 100 ns MD simulations to validate EGFR–ligand complex stability and Khan et al. (2025)24 reported strong binding affinities and stable interactions for novel EGFR ligands using docking and MD. Patil et al. (2025)25 similarly employed QSAR, docking, and MD to identify quinazoline-based EGFR inhibitors. Therefore, targeted inhibition of EGFR tyrosine kinase activity represents a rationally designed and promising strategy to hinder cancer progression. In this study, we employed computational screening to evaluate natural compounds derived from marjoram against the detailed structure of the EGFR tyrosine kinase domain (PDB ID: 2J6M). Molecular dynamics simulations26–28 complement docking analyses by assessing the stability of protein-ligand complexes over time under simulated physiological conditions, thereby probing the robustness of these interactions in a near-biological environment29,30.
This study characterized the phytochemical composition of wild Origanum majorana aerial parts harvested from Saudi Arabia’s Jazan region. The predominant compounds identified through this analysis were subsequently evaluated for their binding potential to EGFR via molecular docking to assess their affinity and mode of interaction within the kinase’s active site. The ligand erlotinib31, a known EGFR inhibitor was redocked into EGFR (2J6M), to confirm the reliability of our docking protocol. Based on their superior docking scores, the three highest-ranking ligands were chosen for further in-depth analysis using molecular dynamics simulations. This approach was implemented to elucidate the precise mechanisms underlying EGFR inhibition by these natural monoterpenes and to identify the most stable and potent candidate for future development as a breast cancer therapeutic. Overall, the results point to the promise of O. majorana metabolites as natural EGFR inhibitors and identify a lead scaffold that may be further explored for novel anti-breast cancer drug discovery.

Results and discussion

Chemical constituents profiling of O. marijuana using GC-MS and GC-FID analysis
Hydro-distillation of the fresh aerial parts of O. majorana afforded a yellow essential oil (OMEO) in a yield of 0.2% (w/w). The yield obtained in the present study is comparable to those reported previously for O. majorana oils from Greece and Iran32. In contrast, it was 2–5 times higher than the yields documented for Tunisian accessions33. Reported yields of O. majorana essential oils vary considerably across different geographic origins, ranging from as low as 0.04% to as high as 12.8%. The maximum yield (12.8%) was obtained from plants cultivated in Cyprus, while the lowest yield was recorded for Tunisian material (Table 1). On average, however, most studies report yields in the range of 0.5–3%34. Such pronounced variability in oil yield among different populations of O. majorana can be attributed to multiple factors, including climatic and edaphic conditions, genetic differences among chemotypes, and the developmental stage of the plant at the time of harvest35. These findings highlight the importance of both environmental and genetic influences on the quantitative expression of essential oils in aromatic plants.

The as-obtained OMEO was analyzed using GC–MS and GC–FID in combination with linear retention indices (LRI), database searches, and comparison with literature data. This approach enabled the identification of 26 volatile constituents, representing 99.4% of the total composition (Table 2). The compounds are listed in order of elution on the HP-5MS nonpolar column. A representative gas chromatography (GC) chromatogram of the OMEO is shown in Fig. 1.

The oil was revealed to be a complex mixture dominated by oxygenated monoterpenes (57.8%), followed by monoterpene hydrocarbons (38.3%), while sesquiterpene hydrocarbons (3.3%) were detected in smaller amounts Fig. 2. Within the oxygenated monoterpene fraction, terpinen-4-ol (26.1%) was the principal constituent, accompanied by trans-sabinene hydrate (19.0%) and cis-sabinene hydrate (4.3%). Other noteworthy members of this group included α-terpineol (3.4%), cis-p-menth-2-en-1-ol (1.6%), and piperitone (1.2%).

The monoterpene hydrocarbon fraction was dominated by γ-terpinene (12.0%), sabinene (8.2%), and α-terpinene (6.8%), with additional contributions from limonene (3.6%), α-terpinolene (2.5%), and β-myrcene (2.4%). Sesquiterpene hydrocarbons were represented primarily by β-caryophyllene (1.8%) and bicyclogermacrene (1.3%). The Jazan province of Saudi Arabia, from which our sample was collected, is characterized by a hot, arid climate with high solar irradiance and limited water availability. Such environmental stresses are known to significantly influence the biosynthesis of secondary metabolites in aromatic plants. In particular, drought and heat stress can alter the activity of key enzymes in the terpenoid biosynthesis pathway (e.g., terpene synthases and cytochrome P450s), often favoring the accumulation of specific oxygenated monoterpenes. The high abundance of terpinen-4-ol and trans-sabinene hydrate, both oxygenated derivatives in our sample may represent a biochemical adaptation to oxidative stress, as these compounds possess established antioxidant properties.
Several of the identified constituents, including terpinen-4-ol, trans-sabinene hydrate, γ-terpinene, sabinene, and α-terpinene, limonene are known to impart the characteristic aroma of O. majorana oil and are frequently linked with diverse pharmacological activities, such as antimicrobial, antioxidant, and anti-inflammatory effects45–47.
Essential oils of O. majorana are known to exhibit chemo-typic variation, with previously described chemotypes including terpinen-4-ol/γ-terpinene, terpinen-4-ol/cis-sabinene hydrate, terpinen-4-ol/trans-sabinene hydrate, terpinen-4-ol/linalool, and carvacrol/thymol. Among these, the terpinen-4-ol/γ-terpinene and terpinen-4-ol/sabinene hydrate types (cis- and trans-) are most frequently reported (Fig. 3). Comparison of the present findings with earlier studies indicates that the investigated oil belongs to the terpinen-4-ol/trans-sabinene hydrate chemotype, similar to oils previously characterized from Albania and Egypt34.

When evaluated as individual marker constituents, terpinen-4-ol has consistently emerged as the predominant component in O. majorana essential oils from a wide range of geographic regions34. Other compounds, including carvacrol, linalyl acetate, sabinene hydrates, γ-terpinene, 1,8-cineole, and pulegone, have also been reported as dominant constituents, though only in a limited number of studies (Fig. 4). Consistent with these observations, the present Saudi Arabian sample was likewise characterized by a high abundance of terpinen-4-ol. Collectively, these findings support the designation of terpinen-4-ol as the principal marker compound of O. majorana essential oil on a global scale.

The six compounds (α-Terpinene, Limonene, γ-Terpinene, Sabinene, Terpinen-4-ol, α-Terpineol) were selected based on their high relative abundance (> 3%) in the GC-MS profile, their documented pharmacological activities in prior literature and their structural diversity within the monoterpene class, allowing for meaningful structure-activity relationship analysis34.
This targeted approach ensures we focus on the most representative and bioactive constituents of O. majorana.

Analysis from molecular docking study
The docking analysis of six phytocompounds against the EGFR revealed differential binding affinities and interaction profiles (Table 3). Among the tested ligands, α-Terpinene exhibited the most favorable docking score (–5.98 kcal/mol), followed by Limonene (–5.32 kcal/mol) and Sabinene (–5.12 kcal/mol). In contrast, α-Terpineol showed the weakest binding affinity (–3.12 kcal/mol). These results suggest that structural variations within the terpene framework substantially influence binding efficiency towards EGFR. The docking score of the reference ligand erlotinib (–5.99 kcal/mol) is comparable to our efficient compounds. Erlotinib is a known EGFR inhibitor, and it proves the docking reliability.

Interaction profiling indicated that hydrophobic contacts, particularly van der Waals and alkyl interactions, were predominant in stabilizing the ligand–protein complexes. For instance, α-Terpinene established van der Waals interactions with LEU, ALA, and VAL, alongside alkyl contacts involving GLY, MET, THR, and LYS. Similarly, Limonene demonstrated extensive hydrophobic interactions with MET, THR, GLU, and multiple LEU residues, consistent with its relatively high docking score. These observations are in agreement with prior studies emphasizing the importance of hydrophobicity in small-molecule binding to EGFR’s active site48. Interestingly, hydrogen bonding was largely absent across most complexes, with the exception of α-Terpineol, which formed a single H-bond with ASP, albeit accompanied by an unfavorable acceptor–acceptor interaction with the same residue. This unfavorable electrostatic contact may contribute to the reduced binding score observed for α-Terpineol. Overall, the docking results highlight that α-Terpinene, Limonene, and Sabinene are the most promising candidates for further molecular dynamics studies, as they not only displayed superior binding energies but also engaged in multiple stabilizing hydrophobic interactions with key residues of EGFR (Fig. 5).

Notably, the docked ligands form critical contacts with the EGFR catalytic hinge. Interactions with critical EGFR catalytic residues (e.g., LEU, VAL, MET, LYS) (Fig. 5) are part of the ATP-binding site and are known to be essential for kinase activity. These hydrophobic pocket residues are hotspot interactions49. The persistent hydrophobic interactions with these residues, especially for limonene, provide a clear mechanistic basis for its inhibitory potential. α-terpinene, limonene, and sabinene were chosen for MD as the top-tier binders (docking scores < − 5.0 kcal/mol) to assess dynamic stability, a standard protocol in computational screening.

Analysis from molecular dynamics study

Effect of ligands of protein conformation
MD simulations were performed to investigate the conformational stability of EGFR in both its apo form and when bound to three phytochemicals, namely α-terpinene, limonene, and sabinene. Structural integrity was evaluated by monitoring root mean square deviation (RMSD), while conformational flexibility was analyzed via residue-based root mean square fluctuation (RMSF). Molecular compactness was quantified using the radius of gyration (Rg) (Fig. 6)50. The RMSD analysis demonstrated that ligand binding imparted greater stability to the EGFR structure relative to its unbound state. Specifically, the α-terpinene complex exhibited moderate structural fluctuations during the initial 30 ns of simulation, subsequently achieving stability with RMSD values ranging from 0.25 to 0.40 nm. Limonene binding resulted in a more stable trajectory, with backbone deviations consistently remaining in the range of 0.20–0.35 nm, suggesting higher conformational stability. Similarly, the sabinene complex maintained lower deviations than the apo structure throughout the 100 ns simulation. Overall, these observations suggest that ligand binding stabilizes EGFR mostly with limonene than sabinene and these two exerting stronger effects than α-terpinene.
The RMSF results further confirmed these observations, demonstrating decreased residue-level flexibility in the ligand-bound forms of EGFR relative to the unbound receptor. While α-terpinene binding decreased fluctuations across most residues, some loop regions remained moderately flexible. Limonene exhibited the greatest reduction in RMSF, suggesting enhanced rigidity in both functional and non-functional regions of EGFR. RMSF peaks occur in the activation loop region (within the C-lobe) and residues in the N-lobe show much lower fluctuations. Sabinene also reduced fluctuations, although its stabilizing influence was slightly weaker than that of limonene. The reduction in flexibility indicates restricted conformational mobility of EGFR upon phytochemical binding, which is often associated with improved structural stability.
The Rg plots demonstrated the effect of ligand binding on protein compactness. Apo-EGFR maintained an average Rg between 2.05 and 2.10 nm, while complexes with limonene and sabinene displayed slightly lower values (~ 1.95–2.05 nm), suggesting enhanced compactness and structural integrity. In contrast, α-terpinene binding did not induce notable changes, maintaining values comparable to apo-EGFR. These results indicate that while α-terpinene stabilizes EGFR to some extent, limonene and sabinene promote a more compact and stable protein conformation.
Taken together, the analyses of RMSD, RMSF, and Rg confirm that all three phytochemicals contribute to EGFR stability, but limonene and sabinene have stronger effects compared to α-terpinene. The enhanced stability and molecular compaction appear to result from optimal hydrophobic contacts within the binding pocket, which restrict conformational dynamics and preserve structural cohesion and previous studies also supports this type of observations17.

Analysis from distance plot with hydrophobic interaction
The distance plot (Fig. 7) illustrates the temporal evolution of the ligand–receptor center-of-mass distance for EGFR in complex with three different ligands: α-terpinene (black), limonene (red), and sabinene (green) over a 100 ns molecular dynamics trajectory. The distance metric provides insight into the stability of ligand binding and potential dissociation events throughout the simulation.
The EGFR–limonene complex (red trace) demonstrates notable structural stability, sustaining a consistent distance of approximately 1.0 nm throughout the simulation duration. This molecular behavior indicates that limonene remains firmly anchored within the active site without significant dissociation, indicative of a high-affinity binding mechanism. In contrast, the EGFR–α-terpinene complex (black trace) exhibits pronounced variations in intermolecular distance, most notably during the 10–30 ns interval, where the separation surpasses 5.0 nm. Deviations of this magnitude are indicative of transient dissociation or low binding affinity. Although the ligand intermittently returns closer to the protein, the lack of stability indicates transient interactions rather than sustained binding. The EGFR–sabinene complex (green trace) demonstrates intermediate behavior. The ligand initially remains close to the binding pocket but undergoes repeated fluctuations after 20 ns, reaching distances above 4.0 nm before returning to near-bound states. This suggests dynamic association–dissociation events, reflecting moderate binding stability relative to limonene and weaker retention compared to a strong binder. Dissociation events observed for α-terpinene and sabinene are indicative of weaker binding affinity and transient interactions. The pocket accommodates all ligands, but only limonene maintains stable binding due to optimal hydrophobic complementarity.
Overall, among the three ligands, limonene shows the strongest and most stable binding to EGFR, while α-terpinene appears least stable, frequently dissociating from the receptor. Sabinene demonstrates partial stability but undergoes significant fluctuations, suggesting less consistent binding.

Hydrophobic contact analysis
Hydrophobic contacts (Fig. 7) are crucial for stabilizing ligands within protein binding pockets, a feature especially relevant for small nonpolar molecules like terpenes51. The hydrophobic contact profiles further support the distance-based findings.
EGFR–Limonene: The complex exhibited consistently high numbers of hydrophobic contacts (~ 15–20) throughout the 100 ns trajectory. This stable interaction profile indicates that van der Waals and nonpolar contacts contribute significantly to the retention of limonene in the binding pocket.

EGFR–α-terpinene
The number of hydrophobic contacts fluctuated strongly, with periods of reduced interactions coinciding with the dissociation events observed in the distance plot. This instability highlights a weaker hydrophobic anchoring of α-terpinene within EGFR.

EGFR–sabinene
The contact profile displayed irregular patterns, with several intervals where hydrophobic interactions dropped to near zero. This correlates with the increased distances observed in the trajectory, reinforcing that sabinene undergoes frequent unbinding and rebinding events, resulting in less stable complex formation.
The three top-ranked ligands (limonene, α-terpinene, and sabinene) were found to bind within the ATP-binding cleft of EGFR and interact with key catalytic residues essential for kinase activity. All ligands established stable hydrophobic contacts with hinge and pocket-forming residues, notably MET, LEU, VAL, and the catalytic LYS. Limonene showed the most persistent interactions with these residues during molecular dynamics simulations, leading to sustained occupation of the ATP-binding site.
Taken together, the distance and hydrophobic contact analyses indicate that limonene exhibits superior stability and affinity toward EGFR compared to α-terpinene and sabinene. The persistent hydrophobic interactions and minimal fluctuation in distance suggest that limonene forms a robust binding mode17,24,25, whereas α-terpinene and sabinene interact less consistently. The ability of limonene to maintain these interactions could therefore explain its enhanced binding stability.

Analysis of ligand position diagrams
The positional snapshots (Fig. 8) of EGFR in complex with α-terpinene, limonene, and sabinene at different time intervals provide a visual confirmation of the binding stability inferred from the distance and hydrophobic interaction analyses. For the EGFR–α-terpinene complex, the ligand shows clear displacement from the binding pocket as early as 25 ns, with subsequent frames (30 ns and 100 ns) confirming its unstable binding and frequent dissociation events. This observation aligns with the distance plot, where α-terpinene exhibited large fluctuations exceeding 5.0 nm, indicating weak retention. The EGFR–limonene complex maintains a consistent ligand orientation within the binding site across all snapshots (0–100 ns). The ligand remains stably positioned without significant drift, corroborating the minimal fluctuations observed in the distance trajectory and the persistent hydrophobic contacts. This strongly suggests that limonene establishes robust and stable interactions with EGFR throughout the simulation. In contrast, the EGFR–sabinene complex displays intermittent displacement, with the ligand leaving the active site around 33 ns and 80 ns before partially re-establishing contact at 100 ns. This behavior supports the earlier findings of dynamic association–dissociation events, reflected by unstable hydrophobic contact patterns and fluctuating distances.
Overall, the positional diagrams visually reinforce the quantitative analyses: limonene demonstrates the highest binding stability with EGFR, while α-terpinene and sabinene exhibit weaker and transient interactions.

Binding energy calculation
The binding affinity was quantified as the binding energy between Limonen and EGFR using MMGBSA method [Table 4]. ‘EGFR + α-Terpinene’ and ‘EGFR + Sabinene’ systems were not considered as they deviated from their initial binding pose. Binding energy results for the EGFR–limonene complex provide deeper insight into the driving forces behind ligand binding. The van der Waals contribution (ΔE_vdw = − 21.28 ± 2.82 kJ/mol) is the dominant stabilizing factor, highlighting the importance of hydrophobic packing between limonene and the EGFR active site. Electrostatic interactions contribute minimally (ΔE_elec = − 1.92 ± 1.55 kJ/mol), consistent with the largely nonpolar nature of limonene. The overall binding free energy (ΔE_binding = − 13.64 ± 2.94 kJ/mol) indicates that limonene forms a stable complex with EGFR, primarily driven by hydrophobic interactions. This finding is in strong agreement with the distance, hydrophobic contact, and positional analyses, all of which demonstrated stable retention of limonene within the EGFR binding pocket.

DFT-based analysis of limonene
To further interpret the binding characteristics of limonene with EGFR, DFT calculations were carried out. Optimized molecular structure confirmed the characteristic cyclic monoterpene framework with conjugated double bonds, a structural feature that contributes to its rigid and hydrophobic nature (Fig. 9). This property explains why limonene displayed strong van der Waals interactions and stable retention within the EGFR binding site in the molecular dynamics (MD) simulations.

Mulliken charge distribution and electrostatic potential
The Mulliken charge distribution (Fig. 5) highlights electron-rich (negative) regions around the double bonds and electron-deficient (positive) regions near hydrogen atoms. This distribution aligns with the electrostatic potential map, which shows localized areas of partial negative potential that could facilitate weak dipole-driven interactions with polar residues, although the dominant interaction mode remains hydrophobic. These electronic features explain why limonene maintains a strong and consistent orientation in the EGFR binding pocket, as seen in the trajectory snapshots.

Electronic properties
Frontier molecular orbital (FMO) analysis indicated that the HOMO is primarily distributed over the conjugated double bonds, whereas the LUMO is delocalized across the same π-framework. The calculated HOMO–LUMO energy gap of 6.87 eV suggests high electronic stability, with relatively low reactivity toward charge transfer processes. This wide gap supports limonene’s preference for nonpolar stabilization through van der Waals and hydrophobic interactions rather than strong polar or electrostatic contributions.

Vibrational (IR) spectrum
The simulated IR spectrum displayed absorption peaks between 3043 and 3227 cm⁻¹, characteristic of C–H stretching vibrations. A sharp band at 1725 cm⁻¹ corresponded to C = C stretching, while additional signals around 1506 and 922 cm⁻¹ were assigned to bending and out-of-plane motions of the double bonds. The observed vibrational frequencies align with previously published infrared spectra for limonene and analogous monoterpenes52,53 thereby corroborating the precision of the computational models employed.

Integration with MD finding
The combination of DFT and MD analyses provides a coherent explanation for limonene’s superior binding behavior. Its stable π-electron framework, large HOMO–LUMO gap, and predominantly hydrophobic character all contribute to strong van der Waals stabilization inside the EGFR active site. These results help explain why limonene, compared to α-terpinene and sabinene, maintained the most stable binding mode across the trajectory simulations.

Global reactivity descriptor (GRD) analysis of limonene
To complement the frontier orbital and electrostatic potential results, global reactivity descriptors (GRDs) were evaluated from Fukui function analysis. These descriptors provide critical information regarding the compound’s inherent reactivity and molecular stability, which correlate with its binding behavior toward biological targets54,55. Fukui function analysis of limonene’s global reactivity descriptors offers deeper insight into its chemical behavior and stability. The compound exhibited an electronic potential of − 0.0877 au and a global hardness of 0.3162 au, suggesting low charge transfer susceptibility. A global softness value of 3.1629 au suggests a moderate capacity for electron redistribution, while an electrophilicity index of 0.3309 eV implies a limited tendency for electrophilic behavior. These values, together with the large HOMO–LUMO gap (6.87 eV), confirm that limonene is electronically stable and chemically less reactive. Thus, the GRD analysis reinforces the conclusion that limonene’s affinity for EGFR arises mainly from hydrophobic stabilization rather than polar interactions54,55.

Integration with binding behaviour
The GRD results confirm that limonene is chemically stable, electronically inert, and unlikely to engage in strong polar interactions with EGFR residues. This characteristic explains the low electrostatic binding energy (ΔE_elec = − 1.92 kJ/mol) calculated in MMGBSA and reinforces the view that its retention in the EGFR active site is primarily due to hydrophobic forces. Collectively, the Fukui-based descriptors and frontier orbital analysis provide a coherent explanation for limonene’s molecular stability and preference for nonpolar binding, making it the most promising of the tested ligands.

ADME and drug-likeness analysis of limonene
The pharmacokinetic and drug-likeness properties of limonene were evaluated (Table 5) using the SWISS ADME web server to assess its potential as a lead-like compound.

Physicochemical properties
Limonene (C₁₀H₁₆, MW = 136.23 g/mol) is a small non-aromatic monoterpene with a high fraction of sp³ carbons (0.60) and minimal conformational flexibility (1 rotatable bond). It lacks hydrogen bond donors and acceptors, which explains its zero topological polar surface area (TPSA = 0.0 Ų). A low TPSA generally favors passive membrane diffusion but limits hydrogen bonding with target proteins56.

Lipophilicity and solubility
Consensus log P was estimated at 3.37, placing limonene in the moderately lipophilic range. Such lipophilicity is consistent with the hydrophobic stabilization observed in MD simulations and MMGBSA binding energy decomposition. Predicted water solubility varied between soluble and moderately soluble classes depending on the model, with ESOL predicting 4.33 × 10⁻² mg/ml (classified as soluble) and Ali predicting lower solubility (moderately soluble).

Pharmacokinetics
Limonene is predicted to have low gastrointestinal absorption but is capable of crossing the blood–brain barrier (BBB), reflecting its high lipophilicity and small molecular size. It is not a substrate for P-glycoprotein (P-gp), reducing the risk of efflux-related bioavailability issues. Among cytochrome P450 enzymes, only CYP2C9 inhibition was predicted, while no inhibitory activity was indicated for CYP1A2, CYP2C19, CYP2D6, or CYP3A4.

Drug-likeness
Limonene satisfies Lipinski’s rule of five with no violations, indicating acceptable oral drug-likeness. It also passes the Veber and Egan filters, but fails the Ghose and Muegge criteria due to its small molecular weight (< 200 Da) and low heteroatom count. Despite these limitations, limonene achieved a bioavailability score of 0.55, suggesting moderate oral bioavailability.

Integration with binding studies
The ADME profile complements the molecular docking, MD, and DFT results by highlighting limonene’s hydrophobic nature, small molecular size, and passive diffusion capability. While its poor aqueous solubility and low GI absorption may limit oral administration, its favorable lipophilicity, BBB permeability, and metabolic selectivity indicates possible suitability for central nervous system oriented therapeutic use. Moreover, the lack of P-gp substrate recognition and limited CYP inhibition profile support its potential as a safe scaffold for drug design. BBB permeability could be relevant for targeting brain metastases but warrants careful preclinical toxicity assessment.

Discussion
In conclusion, this study successfully integrates phytochemical analysis with advanced computational biology to evaluate the anti-breast cancer potential of Origanum majorana essential oil constituents. GC-MS profiling provided a comprehensive chemical inventory, highlighting the abundance of bioactive monoterpenes. The in silico screening pinpointed α-terpinene, limonene, and sabinene as high-affinity ligands for the EGFR tyrosine kinase domain. Molecular dynamics simulations were pivotal in differentiating the stability of these complexes, with limonene emerging as the superior candidate. Its robust binding, driven predominantly by van der Waals forces, was characterized by minimal protein backbone deviation, sustained hydrophobic contacts, and a favorable binding free energy. DFT calculations further rationalized this strong interaction by elucidating limonene’s stable, hydrophobic electronic structure. Complementing these findings, ADME profiling suggested that limonene possesses promising drug-like properties, including high lipophilicity and BBB permeability. Therefore, this study provides a strong mechanistic foundation for limonene as a natural, bioavailable EGFR inhibitor. While limonene shows promising binding and drug-like properties, its simple structure makes it a fragment-like lead suitable for further optimization, such as chemical modification to enhance potency, selectivity, and pharmacokinetic properties.
Future work should focus on the experimental validation of these computational insights through in vitro assays to confirm the efficacy and therapeutic potential of limonene with its optimized derivatives.

Materials and methods

Collection of plant materials
Aerial parts of O. majorana cultivated in the Jazan Province of Saudi Arabia were collected in May 2023. The plant material was authenticated by a plant taxonomist, Dr. Rajakrishnan Rajagopal at the Herbarium Division, King Saud University, Riyadh. A voucher specimen (KSU-24648) was deposited in the university’s herbarium for the record. Collection of plant materials were carried out according to institutional, national, and international guidelines. All necessary permissions were obtained from the relevant authorities for the collection of plant materials. The fresh samples were subsequently processed for essential oil isolation.

Isolation of essential oil
The freshly collected aerial parts of O. majorana were finely chopped and subjected to hydrodistillation for 3 h using a conventional Clevenger-type apparatus, following previously described procedures57. The distillation afforded a light-yellow oil with a yield of 0.2% (w/w, fresh weight basis). The obtained oil was dried over anhydrous Na₂SO₄ and stored at 4 °C until further analysis.

GC–MS and GC–FID analysis
The essential oil was diluted in diethyl ether and analyzed by GC–MS and GC–FID using an HP-5MS nonpolar stationary phase column, as reported58. Detailed analytical conditions are provided in the Supporting Information (Scheme S1). The identified constituents and their relative abundances are listed in Table 2 in the order of elution on the HP-5MS column.

Determination of linear retention indices (LRIs)
The LRIs of the volatile components were calculated according to established procedures58, with full details described in the Supporting Information (Scheme S2).

Characterization of volatile constituents
The volatile constituents of O. majorana essential oil were characterized on the HP-5MS column following standard protocols58. A comprehensive description of the methodology is provided in the Supporting Information (Scheme S3).

Molecular docking study
The three-dimensional structures of the functionally relevant breast cancer-associated proteins were obtained from the Protein Data Bank (http://www.rcsb.org). The details of the receptor protein is presented in Table 6.

For proteins containing unresolved regions, missing residues were computationally modeled using the SWISS-MODEL homology modeling server17. The EGFR crystal structures (PDB IDs 1M17 and 1XKK) contained unresolved segments. These missing residues were reconstructed using homology modeling (SwissModel)24 These completed receptor structures were then prepared for molecular docking by stripping heteroatoms and incorporating polar hydrogens alongside Gasteiger partial charges, utilizing the AutoDock Tools software package50,60. The predominant compounds identified in the GC-MS analysis of O. majorana aerial parts were selected as ligands [Table 2] for molecular docking. These ligand structures were first subjected to geometry optimization and then docked into the binding site using AutoDock Vina61,62, which calculated docking scores to quantify binding affinity. The resulting binding poses were visualized and analyzed with Discovery Studio62 and PyMOL63 to characterize specific interaction types and binding modes. The co-crystallized ligand (erlotinib)31 was redocked into EGFR (2J6M), yielding an RMSD < 2.0 Å, confirming the reliability of our docking protocol.

DFT study
Geometry optimization of the chosen ligands [Table 7] was carried out using the B3LYP hybrid functional64 implemented in the Gaussian 16 software package65. All atoms were treated using the localized 6–31G(d, p) basis set. The optimized ligand geometries were subsequently employed in molecular docking and dynamics simulations. For the most active compound, limonene, several electronic descriptors such as Mulliken charge distribution, electrostatic potential (ESP) surface, and HOMO–LUMO energy gap were determined at the same theoretical level. The vibrational (IR) spectrum of Limonene was simulated with the same method66. Global reactivity descriptors were further obtained using Fukui function analysis to assess the compound’s chemical reactivity67.

Molecular dynamics simulations
Following molecular docking studies, molecular dynamics (MD) simulations were performed for the most efficient 3 ligands as per the docking scores with the same protein EGFR. All simulated systems, including the apo proteins and protein-ligand complexes, were prepared using the CHARMM-GUI web server68. Each system was solvated in a cubic box of TIP3P water69, with a minimum solute-to-box distance of 10 Å to mitigate finite-size effects. Biological ionic strength was maintained at 0.15 M through the addition of potassium and chloride ions. The CHARMM36 all-atom force field70 was applied to the proteins, and the CGenFF force field71 was used for the EGFR ligand. Prior to production dynamics, each system underwent energy minimization for 1000 steps using the Adopted Basis Newton-Raphson algorithm. Equilibration stage involved a 1 ns run under the NVT ensemble, followed by 5 ns under the NPT ensemble. Afterward, production MD simulations were performed for 100 ns using GROMACS 202230 in constant pressure (1 atm) and temperature (300 K). The Particle Mesh Ewald (PME) approach72 was employed to consider for long-range electrostatic interactions, while a cutoff distance of 12 Å was set for short-range van der Waals and Coulomb forces. An integration timestep of 1.0 fs was maintained, and trajectory coordinates were written every 20 ps for analysis. Structural properties were evaluated using native GROMACS tools, and the MM-GBSA method was applied to calculate the binding energy of each complex73.

Binding energy calculation
The binding affinities for the three protein–ligand complexes were determined using the gmx_MMGBSA too73 within the Amber software suite. The binding free energy (ΔE_binding) was computed based on the principles of free energy perturbation and is represented by the following equation:
The energy components were calculated as ensemble averages during the last 50 ns of the simulation trajectories. For EProtein–Ligand, the complex trajectory was generated by stripping water molecules and counter ions from the simulated system. Independent trajectories of the protein and ligand were then extracted from this desolvated complex trajectory to evaluate EProtein and ELigand. The molecular mechanics contributions namely Einternal, Evdw, and Eelec were determined using the same force field applied in the MD simulations70. The polar solvation free energy (Epolar solvation) was calculated using the Generalized Born (GB) implicit solvent model available in the gmx_MMGBSA module74, with dielectric constants assigned as 4 for the solute and 80 for the solvent75. The nonpolar solvation energy (Enonpolar_solvation) was calculated using a linear relation to the solvent-accessible surface area (SASA)76. Time-dependent changes in these energy components were analyzed and visualized using the XMGRACE software package77, while Visualization and structural analysis of the protein–ligand complexes was performed using the Visual Molecular Dynamics (VMD) software78.

Absorption, distribution, metabolism, and excretion (ADME) prediction methodology
The pharmacological profile and drug-like characteristics of the most potent compound, limonene, were assessed by forecasting its Absorption, Distribution, Metabolism, and Excretion (ADME) parameters via the SwissADME online platform79. For this computational assessment, the canonical SMILES notation of limonene was input into the platform, generating an extensive suite of calculated descriptors. The evaluation incorporated key physicochemical properties, including molecular weight, rotatable bond count, and hydrogen bond donor/acceptor capacity. Additional analyses encompassed lipophilicity (Log P o/w), solubility (Log S), and pharmacokinetic behavior such as gastrointestinal absorption, blood-brain barrier permeability, and CYP450 enzyme inhibition potential. Adherence to recognized drug-likeness criteria, notably Lipinski’s Rule of Five, was also evaluated79. These in silico predictions provide crucial preliminary insight into the compound’s potential bioavailability and pharmacokinetic profile before experimental studies.

Supplementary Information
Below is the link to the electronic supplementary material.