Multi-Target Inhibition of Breast Cancer Receptors by Pomegranate Peel Phytochemicals: An Integrated Docking, MD Simulation, and DFT Study

INTRODUCTION

Breast cancer is the most frequently diagnosed malignancy among women worldwide and remains a leading cause of cancer-related mortality. Recent global statistics estimate approximately 2.3 million new cases annually, accounting for nearly one-quarter of all cancers in women, with over 680,000 deaths each year (Breast cancer statistics and resources, 2025). Despite significant advances in early diagnosis and targeted treatment, clinical management of breast cancer remains challenging due to pronounced molecular heterogeneity, hormone dependence, aggressive subtypes, and the frequent emergence of therapeutic resistance (Lei et al., 2026). These limitations underscore the need for novel treatment strategies capable of simultaneously modulating multiple oncogenic pathways. Breast cancer progression is driven by complex signaling networks involving receptor tyrosine kinases and nuclear hormone receptors. Epidermal Growth Factor Receptor (EGFR) and Human Epidermal growth factor Receptor-2 (HER2) are well-established oncogenic drivers that regulate tumor proliferation, survival, and metastasis through pathways such as PI3K/AKT and MAPK. EGFR overexpression is particularly associated with aggressive subtypes, including triple-negative breast cancer, whereas HER2 amplification occurs in approximately 15-20% of cases and is linked to poor prognosis. In parallel, Insulin-like Growth Factor-1 Receptor (IGF-1R) signaling promotes tumor growth and therapeutic resistance through extensive crosstalk with EGFR, HER2, and hormone receptor pathways, reinforcing its importance as a therapeutic target (Kavarthapu et al., 2021). Hormone receptors further contribute to breast cancer complexity. The Progesterone Receptor (PR) serves as a key prognostic marker in hormone receptor-positive tumors, while the Androgen Receptor (AR) is expressed across multiple subtypes, including luminal and triple-negative breast cancers, where its role is context-dependent. The functional interplay between growth factor receptors and hormone receptors limits the efficacy of single-target therapies and highlights the therapeutic potential of multi-target inhibition. Although synthetic inhibitors and monoclonal antibodies targeting these receptors have improved patient outcomes, their long-term effectiveness is often constrained by systemic toxicity, high cost, and acquired resistance. Consequently, increasing attention has been directed toward natural product-based therapeutics, which offer structural diversity, multi-target activity, and improved safety profiles. Punica granatum (pomegranate) is a rich source of bioactive phytochemicals, particularly within its peel, which contains high concentrations of polyphenols, flavonoids, and tannins. Pomegranate peel-derived compounds have demonstrated antiproliferative, pro-apoptotic, anti-angiogenic, and hormone-modulatory effects in breast cancer models (Dathan et al., 2025; Rauf et al., 2025). In this study, Plantaginin, Brevifolin, Brevifolin carboxylic acid, and Olivetolide were selected as representative pomegranate peel phytochemicals based on their structural features, biological relevance, and pharmacological suitability. Their polyphenolic and heterocyclic frameworks support stable non-covalent interactions with both receptor tyrosine kinases and nuclear hormone receptors, suggesting potential multi-target inhibitory activity. Accordingly, a comprehensive in silico strategy was employed to evaluate the therapeutic potential of these ligands against EGFR, HER2, IGF-1R, PR, and AR. Molecular docking was used to predict binding affinity and interaction patterns, molecular dynamics simulations assessed complex stability, Density Functional Theory (DFT) analysis elucidated electronic properties, and ADME profiling evaluated drug-likeness. Collectively, this integrated approach aims to provide molecular-level insights into the multi-target anticancer potential of pomegranate peel-derived phytochemicals and support their further development as promising therapeutic candidates.

MATERIALS AND METHODS

Molecular Docking

Ligand and Protein Preparation

Four bioactive phytochemical compounds-Plantaginin, Brevifolin, Brevifolin carboxylic acid, and Olivetolide were selected from Punica granatum (pomegranate) peel based on their reported pharmacological relevance and anticancer potential. The three-dimensional structures of all ligands were obtained and subjected to geometry optimization to achieve energetically stable conformations. Ligand preparation was carried out using AutoDock Tools (ADT), during which Gasteiger–Marsili partial charges were assigned, non-polar hydrogen atoms were merged, and rotatable bonds were defined to allow conformational flexibility. All prepared ligands were saved in PDBQT format for molecular docking analysis (Dathan et al., 2025; J Renukadevi et al., 2024; Roy & Preetha, 2024).

Five breast cancer–related molecular targets were selected, including epidermal growth factor receptor (EGFR; PDB ID: 3W2S), human epidermal growth factor receptor-2 (HER2; PDB ID: 3RCD), progesterone receptor (PR; PDB ID: 1A28), androgen receptor (AR; PDB ID: 2AM9), and insulin-like growth factor-1 receptor (IGF-1R; PDB ID: 2OJ9). The three-dimensional crystal structures of all target proteins were retrieved from the Protein Data Bank and prepared using AutoDock Tools by removing crystallographic water molecules, heteroatoms, and co-crystallized ligands. Polar hydrogen atoms were added, Kollman united atom charges were assigned, and all receptors were converted into PDBQT format for docking simulations.

Grid Box Definition and Docking Protocol

Docking grid boxes were defined individually for each receptor based strictly on the grid parameters used during docking simulations. For EGFR (3W2S), the grid box center was set at X = 6.30750036239624, Y = 1.012500286102295, and Z = 8.815500259399414 Å, with grid dimensions of 33.95100021362305 × 27.179000854492188 × 30.28499984741211 Å. For HER2 (3RCD), the grid box center was defined at X = 12.480, Y = 2.964, and Z = 28.015 Å, with grid dimensions of 50 × 50 × 48 Å. For PR (1A28), the grid box center was set at X = 22.656999588012695, Y = 9.395999908447266, and Z = 60.33250427246094 Å, with grid dimensions of 23.743999481201172 × 33.1619987487793 × 31.045001983642578 Å. For AR (2AM9), the grid box center was defined at X = 26.961, Y = 5.275, and Z = 4.748 Å, with grid dimensions of 50 × 50 × 50 Å. For IGF-1R (2OJ9), the grid box center was set at X = 6.269000053405762, Y = −6.880500316619873, and Z = 21.01849937438965 Å, with grid dimensions of 27.29400062561035 × 28.43899917602539 × 31.940998077392578 Å. These grid parameters were maintained exactly as used during docking to ensure full reproducibility of the study.Molecular docking simulations were performed using AutoDock software employing the Lamarckian Genetic Algorithm (LGA). Docking parameters were standardized for all ligand–receptor combinations, with the number of output binding modes set to ten and the energy range fixed at 4 kcal/mol. Ligands were treated as flexible, whereas receptors were kept rigid throughout the docking process. Multiple ligand conformations were generated for each docking run and ranked according to predicted binding free energy.

Post-Docking Analysis and Visualization

The docked ligand–receptor complexes were analyzed based on predicted binding energy values and interaction stability. Protein–ligand interaction profiling was carried out using the Protein–Ligand Interaction Profiler (PLIP) to identify and classify non-covalent interactions, including hydrogen bonds, hydrophobic contacts, π–π stacking, π–cation interactions, salt bridges, and water-mediated interactions at the residue level. Three-dimensional visualization and structural analysis of docked complexes were performed using PyMOL molecular visualization software to examine ligand orientation, binding pocket occupancy, and key interacting residues. The combined use of PLIP and PyMOL enabled detailed interpretation of binding mechanisms and facilitated comparative evaluation of the multi-target inhibitory potential of pomegranate peel–derived ligands.

Molecular Dynamics Simulation

All molecular dynamics (MD) simulations were performed using the GROMACS simulation engine to investigate the conformational stability, interaction dynamics, and collective motions of protein–ligand complexes in an explicit solvent environment (Mohan et al., 2025; Saini et al., 2024). Docked protein–ligand structures were parameterized using the AMBER99SB-ILDN all-atom force field, which incorporates refined backbone and side-chain torsional potentials to improve protein conformational accuracy. Ligand parameters, including atom types, bonded terms, and restrained electrostatic potential–derived partial charges, were generated using ACPYPE based on the Antechamber framework to ensure compatibility with AMBER force-field conventions. Each protein–ligand complex was centered in a cubic simulation box under periodic boundary conditions, maintaining a minimum solute–boundary distance of 1.0 nm. Systems were solvated using the SPC water model, and electrostatic neutrality was achieved by adding appropriate counter ions placed via electrostatic potential–guided solvent replacement. Energy minimization was performed using the steepest descent algorithm for up to 50,000 steps until the maximum residual force fell below 1000 kJ·mol⁻¹·nm⁻¹, eliminating steric clashes and unfavorable contacts. The minimized systems were equilibrated under constant number, volume, and temperature conditions using a leap-frog integrator with a 2 fs time step. Positional restraints were applied to protein and ligand heavy atoms during equilibration to allow solvent relaxation. Temperature was maintained at 300 K using a velocity-rescaling thermostat. Non-bonded interactions were treated using the Verlet cutoff scheme with a 1.0 nm cutoff, while long-range electrostatics were computed using the Particle Mesh Ewald method. All covalent bonds involving hydrogen atoms were constrained using the LINCS algorithm, enabling stable integration. Unrestrained production MD simulations were subsequently performed under periodic boundary conditions, and trajectories were saved at regular intervals for analysis. Trajectory analyses included backbone and ligand RMSD to assess global stability, residue-wise RMSF to evaluate local flexibility, radius of gyration to monitor global compactness, and hydrogen bond analysis using geometric distance and angle criteria to quantify interaction persistence. Principal Component Analysis (PCA) was conducted on equilibrated trajectories after removal of overall translational and rotational motions. The covariance matrix of backbone atomic fluctuations was diagonalized to extract dominant collective motions, and projections along the principal components were used to characterize ligand-induced modulation of conformational dynamics and conformational space sampling.

Density Functional Theory

Density Functional Theory (DFT) calculations were performed to evaluate the electronic structure and global reactivity descriptors of the ligand molecules.(Murugesan et al., 2025) All quantum chemical computations were carried out using the PySCF package within the Kohn–Sham DFT framework, employing a hybrid exchange–correlation functional in combination with the split-valence 6-31G(d,p) basis set, which incorporates polarization functions on both heavy atoms and hydrogen atoms. Initial three-dimensional ligand geometries were fully optimized without symmetry constraints to eliminate unfavorable conformations and ensure convergence to stationary points on the potential energy surface. Self-Consistent Field (SCF) calculations were converged using stringent energy and density matrix thresholds to ensure numerical stability and reproducibility. Following SCF convergence, frontier molecular orbital analysis was conducted by extracting the energies of the Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO). The HOMO and LUMO energies and the HOMO–LUMO energy gap were used as key descriptors of electronic stability, chemical reactivity, and charge-transfer potential. Within the Koopmans’ approximation, global conceptual DFT reactivity descriptors were derived, including ionization potential and electron affinity, estimated from the negative HOMO and LUMO energies, respectively. Electronegativity and chemical potential were calculated from the mean of the frontier orbital energies, while chemical hardness was defined as half of the HOMO–LUMO gap and chemical softness as its reciprocal. All orbital energies and derived reactivity indices were reported in electron volts after appropriate unit conversion. Collectively, these DFT-derived parameters provided a quantum-chemical foundation for assessing ligand stability and reactivity and complemented the molecular docking, molecular dynamics, and principal component analyses performed in this study.

ADME and Toxicity Analysis

In silico ADME and toxicity studies were performed to evaluate the pharmacokinetic properties and preliminary safety profile of compounds identified as promising candidates from molecular docking and molecular dynamics simulations. Based on their binding affinity, interaction stability, and conformational behavior observed during the simulation studies, Plantaginin and brevifolin carboxylic acid were selected for further pharmacokinetic and toxicity evaluation. ADME properties were predicted using the SwissADME web-based platform. Canonical SMILES representations of the selected compounds were used as input for the analysis. Physicochemical descriptors, including molecular weight, lipophilicity, topological polar surface area, number of hydrogen bond donors and acceptors, and number of rotatable bonds, were calculated to assess drug-likeness and oral bioavailability. Absorption-related properties, such as gastrointestinal absorption, were estimated using predictive models implemented in SwissADME. Distribution and metabolism were evaluated through lipophilicity-based descriptors and predicted interactions with cytochrome P450 enzymes. Excretion-related tendencies were inferred from molecular size and polarity parameters. Drug-likeness was assessed using established rule-based filters, including Lipinski’s rule of five and Veber’s criteria. In silico toxicity assessment was carried out using the ProTox-3 web server to predict potential toxicological effects of the selected compounds. ProTox-3 utilizes machine learning models trained on experimental toxicity data to estimate acute oral toxicity and classify compounds according to the Globally Harmonized System. Additional toxicity endpoints, including hepatotoxicity, carcinogenicity, mutagenicity, immunotoxicity, and cytotoxicity, were evaluated to provide a preliminary safety profile. The predicted ADME and toxicity results were analyzed in conjunction with molecular docking and molecular dynamics findings to support an integrated evaluation of pharmacokinetic feasibility and safety.

RESULTS

Molecular Docking Analysis

Molecular docking studies were performed to evaluate the binding affinity, binding orientation, and intermolecular interaction profiles of four selected ligands against multiple cancer-associated receptor targets, including IGF-1R, Androgen Receptor, HER2, EGFR, and the Progesterone Receptor. Comparative evaluation of docking scores, binding poses, and interaction stability revealed that, among the four screened ligands, brevifolin carboxylic acid and plantaginin emerged as the most favorable candidates, exhibiting stable binding conformations and well-defined interaction networks within the active sites of different receptors (as shown in Table 1). Consequently, the docking interaction profiles discussed below correspond exclusively to these two top-performing ligands in their respective high-affinity receptor complexes. Docking of brevifolin carboxylic acid within the IGF-1R binding pocket as shown in Figure 1 B demonstrated a stable and well-oriented binding pose, primarily stabilized by hydrophobic interactions involving Leu975, Val983, and Lys1003, with interaction distances ranging from 3.61 to 3.82 Å. These residues contributed to anchoring the ligand within the hydrophobic core of the binding cavity. Additional stabilization was achieved through directional hydrogen bonding interactions with Gln977, Gly978, and Thr1127, exhibiting hydrogen–acceptor distances between 2.24 and 3.63 Å. The presence of both backbone- and side-chain-mediated hydrogen bonds indicates favorable polar complementarity and supports strong ligand accommodation within the IGF-1R active site. Within the androgen receptor ligand-binding domain, plantaginin adopted a stable binding (as shown in Figure 1 A) orientation characterized by extensive hydrophobic contacts with Pro682, Val715, Leu744, Ala748, and Lys808, with interaction distances spanning 3.36–3.98 Å. These hydrophobic interactions played a central role in ligand stabilization within the predominantly nonpolar environment of the AR binding pocket. Polar contributions were provided by hydrogen bonding interactions involving Val685, Gln711, and Pro801. In addition, aromatic π–π stacking interactions with Trp718, characterized by a centroid distance of approximately 5.28 Å, and π–cation interactions involving Arg752, with distances ranging from 4.72 to 4.84 Å, contributed significantly to binding stabilization, highlighting the role of aromatic and electrostatic interactions in AR recognition. Docking of plantaginin to HER2 as shown in Figure 2 A revealed strong ligand accommodation within the receptor binding site, supported by a combination of hydrophobic, polar, and electrostatic interactions. Hydrophobic contacts with Leu726, Ala751, Leu800, and Leu852 were observed at distances between 3.34 and 3.58 Å, contributing to ligand anchoring within the hydrophobic region of the pocket. A dense hydrogen-bonding network involving Gly727, Gly729, Ala730, Lys753, Asp845, Arg849, and Asn850 was detected, with hydrogen–acceptor distances as short as 2.15 Å, indicating strong and directional interactions. Notably, a salt bridge interaction between the ligand carboxylate group and Lys753, with an interaction distance of approximately 3.99 Å, provided additional electrostatic stabilization, reinforcing ligand binding within the HER2 active site.

Click to enlarge

Click to enlarge

For EGFR, plantaginin exhibited a favorable binding mode characterized by hydrophobic interactions with Val726, Lys745, and Thr790, with distances ranging from 3.49 to 3.94 Å as shown in Figure 2B. These interactions were complemented by hydrogen bonds involving Gly721, Ala722, and Phe723, exhibiting hydrogen–acceptor distances between 2.09 and 3.79 Å. Furthermore, the formation of a salt bridge between the ligand and Lys745, with an interaction distance of approximately 5.04 Å, contributed to electrostatic anchoring and enhanced stability of the ligand within the EGFR catalytic pocket. Docking analysis of the progesterone receptor as shown in Figure 3 demonstrated that brevifolin carboxylic acid adopted a stable binding conformation within the ligand-binding domain, supported by hydrophobic interactions with Leu718 and Leu763, with interaction distances of 3.67 and 3.94 Å, respectively. Polar stabilization was achieved through hydrogen bonding interactions with Asn719, Gly722, Gln725, Met759, and Arg766, with hydrogen–acceptor distances ranging from 1.84 to 3.63 Å. Additionally, aromatic π–π stacking interactions with Phe778, characterized by a centroid distance of approximately 5.00 Å, further contributed to ligand stabilization within the PR binding pocket. Collectively, the docking results demonstrate that out of four screened ligands, brevifolin carboxylic acid and plantaginin exhibited receptor-specific favorable binding profiles, with brevifolin carboxylic acid showing stable interactions with IGF-1R and the progesterone receptor, and plantaginin displaying consistent and robust binding across the androgen receptor, HER2, and EGFR. The presence of strong hydrophobic anchoring residues, directional hydrogen bonds, aromatic stacking interactions, and electrostatic contacts provides a robust structural basis for ligand selection and justified the progression of these compounds to subsequent molecular dynamics simulations and pharmacokinetic evaluation.

Click to enlarge

Molecular Dynamics Results: Comparative Protein Backbone RMSD Analysis

The conformational stability of all protein–ligand complexes was assessed by monitoring the backbone Root Mean Square Deviation (RMSD) after least-squares fitting to the initial backbone coordinates over the course of the molecular dynamics simulations (Figure 4). Comparative analysis of RMSD profiles revealed ligand- and receptor-dependent differences in backbone stability and dynamic behavior. The AR–plantaginin complex exhibited a rapid equilibration phase, with RMSD values increasing to approximately 0.10–0.12 nm within the first 5 ns. Following equilibration, the RMSD stabilized predominantly within the range of ~0.13–0.15 nm for the majority of the simulation, with limited fluctuations. Toward the final segment of the trajectory, transient increases reaching ~0.18–0.21 nm were observed; however, no sustained drift was detected, indicating preservation of the overall androgen receptor fold. For the EGFR–plantaginin complex, the backbone RMSD converged early to ~0.15–0.17 nm and remained largely confined within ~0.14–0.18 nm throughout the simulation. Occasional short-lived excursions approaching ~0.20–0.22 nm were observed during the later stages but did not persist, suggesting stable global conformation accompanied by moderate intrinsic flexibility of the EGFR backbone. In the HER2–plantaginin complex, a higher RMSD magnitude was observed relative to AR and EGFR. Following rapid equilibration to ~0.20–0.22 nm, the RMSD fluctuated primarily within ~0.21–0.25 nm across the production phase. Transient deviations approaching ~0.27–0.29 nm occurred intermittently but were followed by rapid re-stabilization, indicating increased flexibility without loss of structural integrity of the HER2 receptor. The IGF-1R–brevifolin carboxylic acid complex displayed the largest RMSD values among the studied systems. After initial equilibration, the RMSD fluctuated within ~0.30–0.38 nm during the early production phase and increased further to ~0.38–0.46 nm between approximately 45 and 85 ns, with occasional peaks approaching ~0.50 nm. Notably, the RMSD decreased during the final phase of the simulation, stabilizing around ~0.28–0.32 nm, indicating convergence toward an alternative conformational state rather than progressive destabilization. In contrast, the PR–brevifolin carboxylic acid complex exhibited the lowest and most stable RMSD profile. Following rapid equilibration to ~0.10–0.12 nm within the first 5 ns, the RMSD remained predominantly within ~0.11–0.14 nm for the majority of the simulation. A gradual increase to ~0.14–0.17 nm was observed during the later stages, without abrupt transitions or directional drift, reflecting minimal ligand-induced perturbation of the progesterone receptor backbone. Overall, the comparative RMSD analysis indicates that all protein–ligand complexes achieved structural equilibration and remained dynamically stable throughout the simulations. The PR–brevifolin carboxylic acid and AR–plantaginin complexes exhibited the highest backbone stability, followed by EGFR–plantaginin and HER2–plantaginin, while IGF-1R–brevifolin carboxylic acid showed enhanced conformational flexibility. Importantly, none of the systems displayed sustained RMSD drift indicative of unfolding or structural collapse. These results support the structural robustness of the selected complexes and justify further analysis of residue-level flexibility, compactness, intermolecular interactions, and essential dynamics.

Click to enlarge

Comparative Analysis of Ligand Positional Stability Relative to Protein Frameworks

The dynamic stability and positional behavior of ligands relative to their respective protein frameworks were evaluated by analyzing ligand Root Mean Square Deviation (RMSD) following least-squares fitting of the protein coordinates throughout the molecular dynamics simulations (Figure 5). Ligand RMSD serves as a robust indicator of ligand mobility, conformational adaptability, and persistence of association with the protein environment under dynamic conditions.

Click to enlarge

In the AR–plantaginin complex, ligand RMSD increased rapidly during the equilibration phase and subsequently stabilized within a narrow range of approximately 0.16–0.22 nm for the majority of the simulation. Only minor, transient deviations were observed, with RMSD values occasionally approaching ~0.25–0.30 nm. The persistence of a well-defined RMSD plateau and absence of sustained directional drift indicate strong positional stability of plantaginin relative to the androgen receptor, with limited conformational flexibility. The EGFR–plantaginin complex exhibited a distinctly higher ligand RMSD profile. Following equilibration, the ligand RMSD fluctuated predominantly within ~0.35–0.50 nm, with intermittent excursions reaching ~0.55–0.65 nm during intermediate simulation stages. These fluctuations reflect increased ligand mobility and adaptive reorientation in response to the intrinsically flexible EGFR framework. Importantly, the RMSD remained bounded and did not display progressive increase, supporting stable ligand association despite enhanced conformational freedom. For the HER2–plantaginin complex, ligand RMSD values demonstrated moderate-to-high variability, primarily ranging between ~0.28 and 0.40 nm. Short-lived increases approaching ~0.45–0.50 nm were observed but were followed by periods of re-stabilization. Notably, a late-stage reduction in RMSD suggested partial convergence toward a more restricted ligand–protein configuration. This behavior indicates reversible ligand flexibility consistent with the dynamic nature of HER2, rather than ligand displacement. The IGF-1R–brevifolin carboxylic acid complex displayed the highest ligand RMSD among all studied systems. After rapid equilibration to ~0.45–0.55 nm, the ligand RMSD fluctuated extensively within ~0.50–0.65 nm, with intermittent excursions approaching ~0.70–0.80 nm during the mid-to-late stages of the simulation. Toward the final phase, a clear reduction in RMSD to ~0.45–0.50 nm was observed, indicating convergence toward a dynamically stable ligand–protein configuration following substantial ligand reorientation. These elevated RMSD values are indicative of pronounced conformational adaptability rather than ligand dissociation. In contrast, the PR–brevifolin carboxylic acid complex exhibited comparatively low ligand RMSD values and the most restrained fluctuation profile. After equilibration, RMSD values remained largely confined within ~0.22–0.30 nm for the majority of the simulation, with only moderate late-stage increases approaching ~0.35–0.38 nm. The limited amplitude of RMSD fluctuations and absence of abrupt transitions indicate high positional stability of the ligand relative to the progesterone receptor. Collectively, comparative ligand RMSD analysis demonstrates that all ligands remained dynamically associated with their respective protein targets throughout the simulations. PR–brevifolin carboxylic acid and AR–plantaginin complexes exhibited the highest ligand positional stability, whereas EGFR– and HER2–plantaginin complexes showed moderate ligand flexibility, and IGF-1R–brevifolin carboxylic acid displayed the greatest conformational adaptability. Importantly, none of the systems exhibited sustained RMSD drift indicative of ligand dissociation or loss of association. The observed differences reflect receptor-specific dynamic environments and ligand-dependent conformational responses, providing a consistent and robust dynamic basis for subsequent interaction, flexibility, and essential dynamics analyses.

RMSF Analysis

Root Mean Square Fluctuation (RMSF) analysis was performed to evaluate ligand-induced modulation of protein backbone flexibility and to understand how local conformational dynamics relate to the biological function of each receptor system. RMSF values were calculated for backbone atoms following least-squares removal of global translational and rotational motions, thereby isolating intrinsic protein flexibility and ligand-dependent dynamic effects (Figure 6). In the AR–plantaginin complex, the protein backbone exhibited predominantly low RMSF values across the structured core, with fluctuations largely confined below ~0.08–0.10 nm. Such restrained flexibility is indicative of a conformationally stable receptor architecture, which is critical for maintaining the integrity of the ligand-binding domain and associated regulatory regions involved in androgen receptor signaling. Elevated RMSF values were primarily restricted to terminal regions, with fluctuations exceeding ~0.25 nm at the N-terminus and ~0.50 nm at the C-terminus, consistent with intrinsically disordered or regulatory segments that are not directly involved in ligand recognition. The preservation of low flexibility within the core suggests that plantaginin binding stabilizes the androgen receptor in a biologically relevant conformation conducive to sustained ligand engagement. The EGFR–plantaginin complex displayed comparatively higher RMSF values in selected regions, reflecting the inherently dynamic nature of EGFR, which undergoes conformational transitions during activation and signal transduction. While the structured core retained RMSF values predominantly below ~0.10 nm, several flexible regions exhibited increased fluctuations (~0.20–0.25 nm). Such localized flexibility is biologically relevant, as conformational adaptability in EGFR is required for allosteric communication and regulatory control of kinase activity. Terminal regions displayed pronounced fluctuations exceeding ~0.50 nm, consistent with regions involved in regulatory interactions rather than ligand stabilization. These results indicate that plantaginin binding preserves core stability while allowing biologically necessary flexibility. In the HER2–plantaginin complex, RMSF analysis revealed a generally stable backbone with most regions fluctuating below ~0.10 nm, suggesting that ligand binding does not disrupt the structural framework of HER2. Moderate RMSF peaks (~0.20–0.30 nm) were localized to loop regions, which are often implicated in receptor dimerization and signaling regulation. The controlled flexibility observed in these regions supports a model in which plantaginin binding maintains HER2 structural integrity while permitting conformational plasticity essential for biological function.

Click to enlarge

The IGF-1R–brevifolin carboxylic acid complex exhibited a largely rigid backbone across the structured domains, with RMSF values remaining below ~0.10–0.12 nm for most regions. Localized increases in flexibility (~0.20–0.30 nm) were observed in surface-exposed segments that are frequently associated with ligand-induced conformational adaptation and downstream signaling interactions. Elevated fluctuations at the N-terminal region (>~0.50 nm) reflect intrinsic flexibility rather than ligand-induced destabilization. Overall, the RMSF profile suggests that brevifolin carboxylic acid binding preserves the conformational stability required for IGF-1R functional regulation while allowing adaptive local motions. The PR–brevifolin carboxylic acid complex demonstrated the lowest overall backbone flexibility among all systems. RMSF values across the majority of the protein remained below ~0.08 nm, indicating a highly rigid and conformationally constrained receptor. Such reduced flexibility is biologically significant, as progesterone receptor activity is closely linked to the stabilization of specific ligand-bound conformations that regulate transcriptional activity. Moderate RMSF peaks (~0.12–0.20 nm) were limited to flexible segments, while terminal regions displayed expected higher mobility (~0.35 nm). The overall RMSF profile supports a ligand-stabilized receptor state with minimal dynamic perturbation. Collectively, RMSF analysis indicates that ligand binding does not induce global structural destabilization in any of the studied complexes. Instead, ligand-induced effects are manifested as localized modulation of flexibility in regions associated with regulatory function, while the structural cores of the receptors remain dynamically stable. AR–plantaginin and PR–brevifolin carboxylic acid complexes exhibit enhanced backbone rigidity, consistent with stabilization of biologically relevant receptor conformations, whereas EGFR– and HER2–plantaginin complexes retain moderate localized flexibility, reflecting the dynamic requirements of receptor signaling mechanisms. These findings, together with RMSD analyses, provide a coherent dynamic framework supporting stable ligand engagement and biologically meaningful protein–ligand interactions.

Radius of Gyration (Rg) Analysis

The radius of gyration (Rg) was monitored throughout the molecular dynamics trajectories to assess global compactness, structural integrity, and folding stability of each protein–ligand complex (Figure 7). Rg profiles for all five systems were well converged, with narrow fluctuation ranges (typically ≤0.05–0.08 nm) and no long-term drift, indicating that ligand binding did not induce large-scale unfolding or structural destabilization. The AR–plantaginin complex exhibited a highly stable Rg profile, equilibrating around ~1.83–1.86 nm with minimal fluctuations, reflecting preservation of the helix-rich ligand-binding domain and reinforcement of a conformationally restrained receptor state. The EGFR–plantaginin system showed slightly higher absolute Rg values (~2.00–2.04 nm), consistent with its larger multi-domain architecture, yet maintained stable compactness with only low-amplitude oscillations attributable to physiological kinase domain breathing. Similarly, the HER2–plantaginin complex displayed rapid convergence and sustained Rg stability (~1.99–2.02 nm), suggesting ligand-induced conformational locking of this intrinsically rigid oncogenic receptor. The IGF-1R–brevifolin carboxylic acid complex demonstrated stable Rg values (~2.00–2.07 nm) with occasional short-lived decreases indicative of ligand-induced compaction rather than destabilization, consistent with controlled inter-domain ordering. The PR–brevifolin carboxylic acid complex showed one of the most tightly regulated Rg profiles (~1.83–1.86 nm), with minimal fluctuation, indicating enhanced structural rigidity and reduced conformational entropy within the ligand-binding domain. Collectively, the Rg analyses confirm that both plantaginin and brevifolin carboxylic acid preserve or enhance global protein compactness across all targets, with no evidence of unfolding or excessive flexibility. From a biophysical perspective, these findings demonstrate maintenance of structural integrity upon ligand binding, while biologically they support ligand-mediated stabilization of receptor conformations less permissive to activation. Overall, the Rg profiles provide robust quantitative evidence for the formation of structurally stable and functionally relevant protein–ligand complexes, reinforcing the therapeutic potential of the studied phytochemicals.

Click to enlarge

Solvent Accessible Surface Area (SASA) Analysis

The Solvent Accessible Surface Area (SASA) was monitored throughout the molecular dynamics trajectories to assess ligand-induced modulation of protein surface exposure, hydration behavior, and global compactness (Figure 8). SASA is a sensitive indicator of protein–solvent interactions and conformational stability, where stable, bounded profiles reflect preserved folding and controlled surface dynamics. Across all five protein–ligand complexes, SASA trajectories were well equilibrated with system-specific mean values and limited fluctuations, indicating that ligand binding did not induce excessive solvent exposure, unfolding, or destabilization. The AR–plantaginin complex exhibited a tightly regulated SASA profile (mean ~128 nm²), with moderate fluctuations corresponding to physiological surface breathing and no sustained expansion, indicating preservation of ligand-binding domain integrity. Similarly, the EGFR–plantaginin system showed higher absolute SASA values (mean ~160 nm²), consistent with the larger multi-domain architecture of EGFR, yet remained stable without progressive solvent exposure, supporting ligand-induced restraint of kinase-domain flexibility. The HER2–plantaginin complex displayed a well-defined and stable SASA distribution (mean ~151 nm²), reflecting maintenance of HER2’s intrinsically rigid architecture and suggesting suppression of conformational freedom associated with aberrant activation. The IGF-1R–brevifolin carboxylic acid complex showed bounded SASA fluctuations (mean ~155 nm²), with minor reductions at later simulation stages indicative of ligand-induced compaction rather than destabilization. Likewise, the PR–brevifolin carboxylic acid complex exhibited one of the most consistent SASA profiles (mean ~131 nm²), reflecting enhanced structural rigidity and preserved surface topology within the progesterone receptor ligand-binding domain. Collectively, the SASA analyses demonstrate that plantaginin and brevifolin carboxylic acid maintain native-like hydration shells while preventing pathological solvent accessibility across all systems. The absence of SASA trends associated with unfolding or destabilization confirms structurally coherent and dynamically stable protein–ligand complexes. From a functional perspective, controlled solvent exposure supports ligand-mediated stabilization of receptor conformations less permissive to activation, reinforcing the biological relevance and inhibitory potential of the studied phytochemicals.

Click to enlarge

Hydrogen Bond (H-Bond) Interaction Analysis

Hydrogen bond (H-bond) analysis was performed over the entire molecular dynamics trajectories as shown in Figure 9 to characterize the temporal persistence, interaction strength, and stabilizing contribution of intermolecular polar contacts formed between the ligands (plantaginin and brevifolin carboxylic acid) and their respective protein targets. Hydrogen bonds constitute a key energetic determinant of protein–ligand recognition, governing ligand orientation, binding specificity, and conformational restraint, particularly in nuclear hormone receptors and receptor tyrosine kinases where precise electrostatic complementarity regulates functional states. The number of intermolecular hydrogen bonds was monitored on a frame-by-frame basis, allowing quantitative evaluation of average H-bond occupancy, fluctuation amplitude, and contact persistence throughout the simulations.

Click to enlarge

In the AR–plantaginin complex, the hydrogen bond population was predominantly maintained between 1 and 2 hydrogen bonds, with occasional short-lived transitions to 3 hydrogen bonds. At least one hydrogen bond was present for the majority of simulation frames, indicating continuous polar anchoring of plantaginin within the androgen receptor ligand-binding pocket. The restricted fluctuation range reflects a geometrically optimized interaction network in which a limited number of well-oriented hydrogen bonds are sufficient to stabilize ligand positioning while maintaining the structural integrity of the helix-rich binding domain, thereby constraining conformational rearrangements associated with receptor activation. The EGFR–plantaginin complex exhibited a comparatively more dynamic hydrogen-bonding profile. The number of hydrogen bonds fluctuated primarily between 1 and 3, with recurrent transient events reaching 4 hydrogen bonds. Despite this variability, hydrogen bonds were continuously re-established following short disruptions, indicating adaptive polar interactions within the intrinsically flexible kinase domain. This behavior reflects accommodation of ligand binding within EGFR domain breathing motions while maintaining sustained electrostatic engagement that stabilizes conformational states less favorable for catalytic activation and autophosphorylation. In the HER2–plantaginin system, hydrogen bond analysis revealed a broader interaction capacity, with hydrogen bond counts spanning 1 to 3 bonds for most of the trajectory and clear transient formation of up to 4 hydrogen bonds. These higher-order hydrogen-bonding events, although short-lived, indicate the ability of plantaginin to transiently engage multiple polar residues within the HER2 binding environment. The overall interaction pattern remains controlled and well-distributed over time, consistent with HER2’s relatively rigid structural framework, while the presence of higher hydrogen-bond counts suggests enhanced polar stabilization without inducing excessive conformational flexibility. The IGF-1R–brevifolin carboxylic acid complex displayed moderate and tightly regulated hydrogen bond variability, with hydrogen bond numbers predominantly oscillating between 1 and 3 hydrogen bonds. These fluctuations reflect the multi-domain architecture and inherent flexibility of IGF-1R, while the consistent presence of hydrogen bonds across the trajectory indicates effective polar engagement mediated by the hydroxyl and carboxyl functionalities of brevifolin carboxylic acid. The interaction profile supports adaptive binding that accommodates local conformational adjustments without over-stabilizing the system. In the PR–brevifolin carboxylic acid complex, the hydrogen bond profile was the most restrained among all systems, with bond counts predominantly limited to a single hydrogen bond, and infrequent, short-duration transitions to two hydrogen bonds. Despite the lower absolute number, these interactions were persistent and repeatedly formed at consistent regions of the trajectory, indicating stable positional anchoring of the ligand. This pattern suggests that binding stabilization in PR is dominated by shape complementarity and hydrophobic packing, with hydrogen bonds acting as precise orientational constraints that limit receptor flexibility. Overall, hydrogen bond analysis across all five complexes demonstrates sustained and system-specific polar interaction networks throughout the molecular dynamics simulations. Plantaginin exhibits greater hydrogen-bonding adaptability in receptor tyrosine kinases, particularly EGFR and HER2, while brevifolin carboxylic acid maintains controlled and persistent hydrogen bonding in IGF-1R and PR. The observed differences in hydrogen-bond occupancy and fluctuation amplitudes directly reflect receptor architecture and flexibility, providing molecular-level insight into stable ligand engagement and modulation of receptor conformational dynamics.

Minimum Distance Analysis

The minimum distance between protein and ligand atoms was monitored throughout the molecular dynamics trajectories to quantitatively assess binding persistence, interfacial compactness, and stability of protein–ligand interactions (Figure 10). This metric provides a sensitive indicator of binding-site retention and the continuity of short-range noncovalent interactions, including hydrogen bonding, van der Waals contacts, and electrostatic interactions. Sustained low minimum distances are indicative of robust ligand anchoring, whereas prolonged increases may reflect binding-site rearrangement or partial dissociation. Across all five systems, minimum distance trajectories remained well converged and consistently within sub-nanometer ranges (~0.15–0.25 nm), demonstrating persistent ligand engagement and stable protein–ligand interfaces throughout the simulation period. The AR–plantaginin complex exhibited a tightly confined distance profile (mean ~0.19 nm) with no systematic drift, indicating deep and stable ligand embedding within the androgen receptor ligand-binding cavity. Similarly, the EGFR–plantaginin system showed slightly broader but well-controlled fluctuations (mean ~0.18–0.19 nm), reflecting intrinsic kinase domain flexibility while maintaining continuous ligand association. The HER2–plantaginin complex displayed a uniformly stable minimum distance distribution (mean ~0.19 nm), consistent with sustained interfacial locking and restricted conformational freedom.The IGF-1R–brevifolin carboxylic acid complex exhibited modestly higher variability (mean ~0.20 nm), consistent with the flexible architecture of IGF-1R; however, transient expansions rapidly returned to low-distance values, confirming stable ligand tethering within the binding region. The PR–brevifolin carboxylic acid complex showed one of the most persistent and narrowly distributed distance profiles (mean ~0.21 nm), indicating strong and continuous ligand–receptor contacts and restricted binding-site flexibility. Collectively, all complexes maintained stable sub-nanometer contact distances with average values below ~0.21 nm and no evidence of ligand dissociation. These results confirm strong intermolecular coupling, effective binding-site occupancy, and adaptive ligand accommodation to protein motions without loss of interaction strength. The sustained proximity observed across systems supports conformational stabilization of receptor architectures and suppression of activation-associated structural transitions, reinforcing the dynamic stability and mechanistic relevance of the studied protein–ligand complexes.

Click to enlarge

Number of Protein–Ligand Contacts (< 0.6 nm) Analysis

The number of protein–ligand contacts within a 0.6 nm cutoff was monitored throughout the molecular dynamics trajectories to assess the persistence, density, and temporal stability of non-covalent interactions governing complex integrity (Figure 11). This metric captures cumulative short-range interactions, including hydrophobic, van der Waals, and polar contacts, and serves as a robust indicator of binding stability beyond individual hydrogen bonds. Across all five systems, contact profiles were well equilibrated and temporally consistent, indicating sustained ligand engagement without evidence of detachment or destabilization.

Click to enlarge

The AR–plantaginin complex exhibited a high and stable contact population (mean ~2,050 contacts), reflecting dense hydrophobic packing and persistent anchoring of plantaginin within the androgen receptor ligand-binding domain. Similarly, the EGFR–plantaginin system showed the highest absolute contact numbers (mean ~2,550 contacts), consistent with the larger kinase surface and extensive ligand–protein interface. Despite transient fluctuations associated with intrinsic kinase flexibility, contact density remained robust, indicating broad and stable engagement of the EGFR binding region. The HER2–plantaginin complex displayed a well-defined and stable contact profile (mean ~2,000 contacts), supporting persistent ligand binding and maintenance of a conformationally restrained receptor state. The IGF-1R–brevifolin carboxylic acid complex exhibited moderate but highly stable contact numbers (mean ~1,000 contacts), consistent with the flexible architecture of IGF-1R while maintaining secure ligand accommodation. Likewise, the PR–brevifolin carboxylic acid complex showed a consistently high and narrowly distributed contact population (mean ~1,520 contacts), indicating tight ligand packing and reduced conformational freedom within the progesterone receptor ligand-binding domain. Collectively, the contact number analyses confirm that both plantaginin and brevifolin carboxylic acid form dense, persistent, and well-maintained interaction networks with their respective targets. The absence of progressive contact loss or large-scale disengagement across all systems corroborates RMSD, radius of gyration, SASA, and minimum-distance analyses, reinforcing the formation of structurally coherent, dynamically stable, and energetically favorable protein–ligand complexes capable of effective receptor modulation.

Principal component analysis (PCA)

Principal component analysis (PCA) was employed to systematically characterize and compare the dominant collective motions governing the conformational dynamics of all five protein–ligand complexes AR plantaginin, EGFR plantaginin, HER2 plantaginin, IGF-1R brevifolin carboxylic acid, and PR brevifolin carboxylic acid over the full molecular dynamics trajectories as shown in Figures 12 and 13. PCA decomposes the correlated atomic fluctuations into orthogonal eigenvectors, where the first two principal components (PC1 and PC2) capture the majority of large-amplitude, low-frequency motions that are most relevant to functional conformational changes in receptor systems. Across all complexes, the projections along PC1 and PC2 display continuous, smooth fluctuations without abrupt discontinuities, indicating well-equilibrated trajectories and the absence of large-scale unfolding or destabilizing transitions. The temporal evolution of PC1 in each system is dominated by gradual drifts followed by stable oscillatory behavior, reflecting global collective motions such as breathing, hinge-like domain rearrangements, and ligand-binding pocket adjustments. PC2 consistently captures secondary motions with comparatively lower amplitudes, corresponding to localized flexibility and fine structural adaptations that complement the dominant PC1 dynamics. For the AR plantaginin complex, the PC1 PC2 projections are confined within a relatively narrow range, and the two-dimensional phase-space representation reveals compact clustering with limited dispersion. This indicates that plantaginin binding restricts androgen receptor dynamics to a stable conformational ensemble, suppressing large-amplitude transitions while permitting essential low-frequency motions required for structural accommodation within the ligand-binding domain. The EGFR plantaginin complex exhibits broader sampling along PC1 and PC2 compared to AR, consistent with the intrinsically higher flexibility of receptor tyrosine kinase domains. The 2D projection reveals multiple partially overlapping conformational basins connected through continuous trajectories, indicating dynamic interconversion among closely related substates. This behavior reflects adaptive collective motions of EGFR that accommodate ligand binding while maintaining overall structural integrity and avoiding transitions into highly active conformations. In the HER2 plantaginin system, PCA reveals a comparatively restrained dynamical profile. The PC1 and PC2 projections show reduced amplitude fluctuations and tighter clustering in the 2D conformational space, highlighting the conformational rigidity of HER2 relative to EGFR. The limited exploration of phase space suggests that plantaginin stabilizes HER2 in a restricted ensemble, dampening collective motions associated with receptor dimerization and activation. The IGF 1R brevifolin carboxylic acid complex displays moderate conformational heterogeneity along PC1, with PC2 capturing localized, reversible fluctuations. The corresponding 2D projection shows an asymmetric but continuous conformational distribution, indicating directional collective motions rather than random diffusion across phase space. This pattern reflects the multi-domain architecture of IGF-1R, where ligand binding allows coordinated domain rearrangements while preserving a stable global fold and sustained receptor engagement. For the PR brevifolin carboxylic acid complex, PCA reveals the most compact conformational landscape among all systems. The PC1 and PC2 time series are characterized by low-amplitude, tightly bounded oscillations, and the 2D projection shows a densely populated single basin with minimal dispersion. This indicates strong conformational confinement of the progesterone receptor upon ligand binding, with collective motions restricted to subtle correlated fluctuations rather than large-scale transitions.

Click to enlarge

Click to enlarge

Taken together, the PCA results demonstrate that both plantaginin and brevifolin carboxylic acid modulate receptor dynamics in a receptor-specific manner, balancing stability and flexibility according to intrinsic protein architecture. Hormone receptors (AR and PR) exhibit highly confined conformational sampling upon ligand binding, whereas receptor tyrosine kinases (EGFR, HER2, and IGF-1R) retain controlled dynamic adaptability. The absence of widely separated clusters or abrupt transitions in any system confirms that all complexes remain dynamically stable throughout the simulations, with ligand binding biasing collective motions toward well-defined, energetically favorable conformational ensembles.

Density Functional Theory (DFT)

Density Functional Theory (DFT) calculations were carried out to elucidate the electronic structure, frontier molecular orbital characteristics, and global reactivity descriptors of the investigated ligands, plantaginin (Mol1) and brevifolin carboxylic acid (Mol2). The analysis focused on HOMO–LUMO energy levels, energy gaps, dipole moments, and derived conceptual DFT parameters, including ionization potential, electron affinity, chemical hardness, softness, and electronegativity, which collectively describe molecular stability, charge distribution, and reactivity (as shown in Table 2).

Plantaginin (Mol1) exhibited a highest occupied molecular orbital (HOMO) energy of −5.81 eV and a lowest unoccupied molecular orbital (LUMO) energy of −1.35 eV, resulting in a HOMO–LUMO energy gap of 4.46 eV as shown in Figure 14. This relatively wide energy gap indicates a high degree of electronic stability and a reduced propensity for spontaneous electronic excitation. The HOMO (as shown in Figure 15) density is predominantly localized over the aromatic rings and oxygen-rich functional groups, identifying these regions as primary electron-donating centers, whereas the LUMO (as shown in Figure 16) is mainly distributed across the conjugated π-system and carbonyl functionalities, suggesting preferred sites for electron acceptance. Using Koopmans’ approximation, the Ionization Potential (IP) and Electron Affinity (EA) of plantaginin were estimated to be 5.81 eV and 1.35 eV, respectively. The corresponding chemical hardness (η) was calculated as 2.23 eV, while the global softness (S) was 0.224 eV⁻¹, reflecting a chemically hard and structurally robust molecule with controlled reactivity. The electronegativity (χ) value of 3.58 eV indicates balanced electron-donating and electron-withdrawing tendencies. The calculated dipole moment of 3.24 D suggests moderate polarity and directional charge separation within the molecular framework. In contrast, brevifolin carboxylic acid (Mol2) displayed a HOMO energy of −6.41 eV and a LUMO energy of −2.22 eV, yielding a comparatively smaller HOMO–LUMO gap of 4.19 eV (as shown in Figure 17). The reduced energy gap indicates enhanced electronic polarizability and increased chemical reactivity relative to plantaginin. The HOMO (as shown in Figure 18) electron density is largely concentrated over the phenolic hydroxyl groups and conjugated aromatic system, while the LUMO (as shown in Figure 19) density is strongly localized around the carbonyl and carboxylic acid moieties, highlighting efficient intramolecular charge-transfer pathways. The ionization potential and electron affinity of brevifolin carboxylic acid were calculated as 6.41 eV and 2.22 eV, respectively. The chemical hardness was determined to be 2.10 eV, lower than that of plantaginin, while the global softness increased to 0.238 eV⁻¹, confirming a softer and more reactive electronic character. The electronegativity value of 4.32 eV reflects a stronger electron-withdrawing nature, consistent with the presence of multiple oxygenated functional groups. Notably, brevifolin carboxylic acid exhibited a significantly higher dipole moment of 8.46 D, indicating pronounced molecular polarity and substantial charge separation. Overall, the DFT results clearly differentiate the two ligands in terms of electronic behavior. Plantaginin demonstrates greater electronic stability, higher hardness, and moderate polarity, consistent with a rigid and well-stabilized electronic structure. In contrast, brevifolin carboxylic acid shows increased softness, stronger polarity, and enhanced charge-transfer capability, arising from its dense network of hydroxyl and carboxyl functionalities. These intrinsic quantum-chemical characteristics provide a robust electronic foundation for understanding their distinct interaction profiles and stability trends observed in subsequent computational analyses.

Click to enlarge

Click to enlarge

Click to enlarge

Click to enlarge

Click to enlarge

Click to enlarge

ADME and Toxicological Profiling (SwissADME and ProTox-III)

Comprehensive SwissADME and ProTox-III analyses revealed that both plantaginin and brevifolin carboxylic acid possess favorable safety profiles but exhibit polarity-driven pharmacokinetic constraints typical of polyphenolic natural products. Plantaginin, a high–molecular weight polyphenolic glycoside (MW = 448.38 g/mol), displayed an exceptionally elevated TPSA (~190.3 Ų) and a markedly low consensus logP (−0.19), confirming a strongly hydrophilic character. Consistent solubility predictions across ESOL, Ali, and SILICOS-IT models classified plantaginin as soluble to moderately soluble ; however, gastrointestinal absorption was predicted to be low due to excessive polarity, extensive hydrogen-bonding capacity (7 donors, 11 acceptors), and limited passive membrane permeability. Distribution analysis further indicated an absence of blood–brain barrier permeability and no predicted P-glycoprotein substrate liability, suggesting restricted CNS exposure and minimal efflux-related pharmacokinetic loss. Brevifolin carboxylic acid exhibited a lower molecular weight (292.20 g/mol) and reduced TPSA (~145 Ų), resulting in a more compact yet still hydrophilic architecture, with lipophilicity predictions converging on a near-neutral consensus logP (~0.04). SwissADME models consistently classified it as soluble to very soluble, supporting efficient dissolution and molecular availability. Despite these advantages, gastrointestinal absorption was likewise predicted to be low, reflecting persistent polarity and hydrogen-bonding capacity (4 donors, 8 acceptors). Similar to plantaginin, brevifolin carboxylic acid was predicted to be non-BBB permeant and not a P-glycoprotein substrate, indicating limited CNS exposure and negligible transporter-mediated efflux liability.Neither compound was predicted to inhibit major cytochrome P450 isoforms (CYP1A2, CYP2C9, CYP2C19, CYP2D6, or CYP3A4), indicating a low risk of metabolic interference and drug–drug interactions. Drug-likeness assessment showed that plantaginin violates Lipinski, Veber, and Muegge criteria primarily due to excessive TPSA and hydrogen-bonding features, resulting in a low predicted bioavailability score (0.17). In contrast, brevifolin carboxylic acid satisfied Lipinski, Ghose, and Muegge rules, with only Veber and Egan violations linked to TPSA, and exhibited a substantially higher predicted bioavailability score (0.56), indicating improved systemic exposure potential. Medicinal chemistry filters identified one PAINS and one Brenk alert for plantaginin and a single Brenk alert for brevifolin carboxylic acid, all associated with catechol-like substructures commonly observed in polyphenolic scaffolds and generally linked to redox behavior rather than intrinsic toxicity.ProTox-III predictions indicated low acute and chronic toxicity for both compounds. Plantaginin showed low predicted risks for hepatotoxicity, cardiotoxicity, mutagenicity, carcinogenicity, immunotoxicity, and cytotoxicity, while brevifolin carboxylic acid exhibited similarly low risk across hepatotoxic, nephrotoxic, cardiotoxic, mutagenic, carcinogenic, immunotoxic, and cytotoxic endpoints. Pathway-level analysis further indicated minimal interaction with nuclear hormone receptors, stress-response signaling pathways, and central nervous system ion channels, reducing concerns related to endocrine disruption, neurotoxicity, or off-target receptor engagement. Overall, both compounds combine high aqueous solubility, negligible CYP inhibition liability, and low predicted systemic toxicity, while oral absorption remains limited by elevated polarity. Plantaginin exhibits more pronounced ADME constraints due to molecular complexity, whereas brevifolin carboxylic acid demonstrates a superior balance between physicochemical properties, metabolic compatibility, bioavailability, and safety, positioning it as the more attractive lead scaffold for further optimization and pharmacological development.

DISCUSSION

The present study provides a mechanistically integrated in silico evaluation of pomegranate peel–derived phytochemicals as multi-target modulators of key breast cancer–associated receptors, addressing molecular heterogeneity, signaling redundancy, and adaptive resistance that limit the durability of conventional targeted therapies. By concurrently targeting receptor tyrosine kinases (EGFR, HER2, and IGF-1R) and hormone-responsive nuclear receptors (PR and AR), this work adopts a systems-level strategy that reflects the extensive pathway crosstalk characteristic of breast cancer. The integrated application of molecular docking, Molecular Dynamics (MD) simulations, Density Functional Theory (DFT), and ADME–toxicological profiling establishes a unified framework linking binding recognition, dynamic stability, electronic behavior, and translational relevance. Docking analyses revealed receptor-selective polypharmacology rather than nonspecific promiscuity. Plantaginin exhibited preferential binding toward EGFR, HER2, and AR, whereas brevifolin carboxylic acid showed stronger affinity for IGF-1R and PR. These complementary binding preferences arise from differences in ligand architecture, including hydrogen-bonding patterns, aromatic surface distribution, and electrostatic polarity, enabling effective accommodation within both kinase active sites and nuclear receptor ligand-binding domains. Such selective multi-target engagement is advantageous in breast cancer, as it allows concurrent modulation of proliferative signaling, hormone-driven transcription, and compensatory survival pathways associated with therapeutic resistance.MD simulations supported the stability and biological plausibility of the docked complexes. All ligand–receptor systems rapidly reached equilibrium and maintained structural integrity throughout the simulation period, indicating that ligand binding does not disrupt receptor architecture. Hormone receptor complexes displayed enhanced conformational stabilization consistent with restrained transcriptional activity, while receptor tyrosine kinases retained controlled flexibility characteristic of inactive or signaling-incompetent states. Ligand mobility observed within kinase binding pockets reflected adaptive intra-pocket rearrangement rather than dissociation, supporting sustained receptor engagement under dynamic conditions. Global and residue-level structural analyses further indicated that ligand binding reinforces receptor compactness without inducing destabilizing perturbations. Ligand-induced flexibility was confined primarily to surface and regulatory regions, while structured cores remained stable. These effects are particularly relevant for EGFR and HER2, where excessive conformational freedom underlies constitutive activation, suggesting that stabilization of compact receptor conformations contributes to functional inhibition. Persistent intermolecular interactions underpinned complex stability across all systems. Hydrogen bonding and non-covalent contact analyses demonstrated receptor-dependent interaction strategies, with plantaginin forming adaptable interaction networks in kinase domains and brevifolin carboxylic acid establishing fewer but more persistent polar interactions in hormone receptors. Such long-lived interaction profiles are consistent with effective multi-pathway modulation. DFT calculations provided an electronic rationale for the observed binding and dynamic trends. Plantaginin exhibited higher electronic stability, characterized by a larger HOMO–LUMO gap and greater chemical hardness, supporting sustained and stable receptor engagement. In contrast, brevifolin carboxylic acid displayed increased electronic softness and a higher dipole moment, indicative of enhanced electrostatic responsiveness and charge-transfer capability. These complementary electronic features rationalize their differential receptor selectivity and adaptive binding behavior. Pharmacokinetic and toxicological predictions further supported the translational potential of the identified leads. Both compounds demonstrated high aqueous solubility, minimal cytochrome P450 inhibition, and low predicted risks of major toxicological liabilities. Although elevated polarity may limit passive gastrointestinal absorption-particularly for plantaginin-this limitation is characteristic of polyphenolic natural products and can be addressed through rational optimization or formulation-based delivery strategies. Notably, brevifolin carboxylic acid exhibited a more favorable overall drug-likeness profile, identifying it as a promising candidate for further development. In summary, this study provides compelling molecular-level evidence that plantaginin and brevifolin carboxylic acid function as structurally stable, dynamically compatible, and electronically complementary multi-target modulators of key breast cancer–associated receptors. Their capacity to simultaneously engage growth factor and hormone signaling pathways, combined with favorable safety predictions, supports their potential to overcome limitations associated with single-target therapies. While experimental validation is required, the present findings establish a strong foundation for the rational optimization and future biological evaluation of pomegranate peel–derived phytochemicals as multi-pathway therapeutic candidates for breast cancer.

CONCLUSION

This study employed an integrated insilico approach to investigate pomegranate peel derived phytochemicals as potential modulators of breast cancer associated molecular targets. Molecular docking, molecular dynamics simulations, density functional theory calculations, and ADME–toxicological predictions were combined to assess ligand–receptor interactions, conformational stability, electronic properties, and preliminary pharmacokinetic suitability. Among the evaluated compounds, plantaginin and brevifolin carboxylic acid demonstrated favorable binding behavior across receptor tyrosine kinases (EGFR, HER2, and IGF-1R) and hormone-responsive nuclear receptors (PR and AR). Molecular dynamics analyses indicated stable protein–ligand complexes with preserved structural integrity and controlled conformational flexibility. Density functional theory calculations provided insight into the electronic characteristics of the ligands, while insilico ADME and toxicity analyses suggested acceptable safety profiles, with limitations related primarily to predicted oral absorption. Overall, the findings provide a molecular-level basis for the further investigation of pomegranate peel–derived phytochemicals as multi-target agents in breast cancer. The results support subsequent experimental validation and optimization to better define their therapeutic potential.