Garuga pinnata and Mitragyna speciosa were selected based on their well-documented ethnomedicinal use among tribal populations of the Narmada Valley region (Alirajpur, Jhabua, and Dhar, Madhya Pradesh, India), where they are traditionally employed in the management of inflammation, pain, and fatigue-associated conditions15. Ethnopharmacological knowledge provides a biologically relevant framework for identifying candidate molecules with potential activity against complex, multi-pathway diseases such as cancer16. A systematic phytochemical profiling approach was performed using curated chemical databases, including PubChem, IMPPAT, and ChEMBL17. Compounds were prioritized based on reported pharmacological relevance, structural diversity, and the availability of validated chemical structures. A total of fifteen phytochemicals was selected, comprising seven compounds from G. pinnata and eight from M. speciosa. These compounds represented distinct chemical classes, including diterpenoids and indole alkaloids, which are known for their diverse pharmacological activities. A summary of the selected compounds, their classification, and selection rationale is presented in Table 1.

All ligand structures were retrieved in SDF format and subjected to geometry optimization using Open Babel (v2.4.1) to obtain energetically stable conformations prior to docking18. The optimized structures were converted into PDBQT format to ensure compatibility with docking algorithms. Drug-likeness evaluation was performed according to Lipinski’s Rule of Five, assessing molecular weight, lipophilicity (logP), hydrogen bond donors and acceptors, and rotatable bond count19. This screening ensured that the selected phytochemicals possessed favorable physicochemical properties consistent with oral bioavailability and pharmacokinetic suitability.

Four metastasis-associated protein targets were strategically selected based on their well-established roles in regulating epithelial–mesenchymal transition (EMT), tumor invasiveness, and metastatic progression: AXL kinase (PDB ID: 5U6B), TGF-β1 (PDB ID: 3KFD), TGF-β receptor II (PDB ID: 3DKC), and TGF-β receptor I kinase (PDB ID: 3ODU)20. The three-dimensional structures of these receptor proteins are illustrated in Figure 1. The simultaneous targeting of these receptors enabled the exploration of a multi-pathway inhibition strategy within the AXL–TGF-β signaling axis, which is known to synergistically drive metastatic phenotypes. Three-dimensional crystal structures were retrieved from the Protein Data Bank (PDB) and subjected to standardized preprocessing to ensure structural reliability for docking simulations. This included the removal of crystallographic water molecules, heteroatoms, and co-crystallized ligands to eliminate potential steric and electrostatic interference21. Polar hydrogen atoms were subsequently added, and structural optimization was performed using Discovery Studio Visualizer to stabilize protein geometry and prepare the receptors for interaction analysis. Binding pocket identification was performed using complementary computational tools, including CASTp and CB-Dock, allowing accurate detection of ligand-accessible cavities based on geometric and topological features22. The identified binding sites exhibited significant variability in pocket volume and residue composition, ranging from 630 ų (TGF-βRII) to 3503 ų (TGF-β1), highlighting differences in ligand accommodation capacity and physicochemical environments across targets. A comparative summary of receptor characteristics, including structural properties and binding pocket features, is presented in Table 2, providing a foundation for understanding differential ligand binding behavior across the selected targets.

Molecular docking was performed using CB-Dock (v2.0), which integrates cavity detection with Auto_Dock Vina scoring to predict ligand-binding conformations23. A blind docking approach was employed to identify both canonical and potential allosteric binding sites without predefined bias. The docking search space for each receptor was defined based on the top-ranked predicted cavities, and detailed grid parameters including center coordinates and box dimensions were provided in Table 3 to ensure reproducibility. Binding affinities were expressed as estimated Gibbs free energy (ΔG, kcal/mol), where more negative values indicated stronger predicted interactions. Post-docking interaction profiling, including hydrogen bonding, hydrophobic contacts, and π-mediated interactions, was performed using Discovery Studio. Ligand efficiency was calculated to normalize binding affinity relative to molecular size.

Given the blind docking framework employed by CB-Dock, conventional validation through redocking of co-crystallized ligands was not uniformly feasible across all receptor systems. To address this limitation, a residue-based validation strategy was implemented, wherein predicted ligand-binding poses were evaluated against literature-reported active site residues of AXL kinase and TGF-β receptors26. The consistent involvement of key residues—such as Asp70, Glu74, Arg163, and Tyr185—supports the biological relevance and structural plausibility of the predicted binding modes. Furthermore, docking poses and interaction geometries were independently examined using Discovery Studio Visualizer to verify hydrogen bonding patterns, steric compatibility, and overall binding orientation within the receptor cavities. The observed interaction profiles, including hydrogen bonding and π-mediated contacts within known functional domains, were consistent with established ligand–receptor interaction principles. Although this approach provided indirect validation, it is consistent with accepted practices in blind docking workflows where experimentally resolved reference ligands are unavailable27. Nevertheless, it is acknowledged that docking predictions represent static approximations; therefore, future studies will incorporate re-docking protocols, consensus docking approaches, and molecular dynamics simulations to enhance predictive reliability and dynamic accuracy.

To complement docking-based affinity predictions, binding free energy calculations were performed using MM/GBSA and MM/PBSA approaches, which integrate molecular mechanics energy terms with implicit solvation models to evaluate ligand–receptor interaction energetics28. These methods were applied to comparatively assess binding favorability and to identify key residue-level energetic contributions involved in stabilizing the docked complexes.

The calculated free energy profiles generally supported the molecular docking results, in which top-ranked phytochemicals demonstrated favorable energetic stabilization within the receptor binding pockets. However, it is important to emphasize that these calculations were performed on energy-minimized static receptor–ligand complexes rather than on conformational ensembles generated through molecular dynamics simulations. Therefore, the obtained ΔG values should be interpreted as relative and comparative estimates of binding favorability rather than absolute thermodynamic parameters.

In addition, the present approach did not account for time-dependent conformational flexibility, explicit solvent dynamics, or entropic contributions, all of which can influence the accuracy of free energy estimation under physiological conditions29. Accordingly, the MM/GBSA and MM/PBSA analyses in this study were used to support docking-derived interaction trends and to prioritize structurally compatible phytochemicals for further investigation, rather than to infer definitive binding stability. Future studies incorporating trajectory-based molecular dynamics simulations and advanced free energy calculations will be necessary to provide deeper insight into kinetic stability and thermodynamic behavior of the identified ligand–receptor complexes.

Root Mean Square Deviation (RMSD) analysis was performed using Discovery Studio to evaluate the structural consistency and conformational deviation of receptor–ligand complexes following energy minimization. RMSD values below 2.5 Å were considered indicative of limited structural fluctuation and consistent ligand accommodation within the predicted binding pockets. The analysis was used to assess the reliability of docking-derived binding geometries and overall pose consistency across receptor systems. However, it is important to note that these RMSD measurements were derived from static, energy-minimized complexes and therefore reflect structural consistency rather than time-dependent dynamic stability under physiological conditions30.

Receptor–ligand complexes were visualized using Discovery Studio Visualizer. Two-dimensional and three-dimensional interaction maps were generated to illustrate hydrogen bonds, hydrophobic interactions, and electrostatic contacts within the binding pockets. Structural visualization facilitated the identification of key residues involved in ligand stabilization and provided insights into binding orientation and interaction networks31.

The computational framework employed in this study comprised a structured, multi-step pipeline integrating receptor selection, ligand preparation, blind molecular docking, interaction profiling, binding free energy estimation, and structural deviation analysis (Figure 2, 3). Ethnomedicinally relevant phytochemicals were initially identified and prepared through energy minimization and format conversion, while target receptor structures were retrieved and processed to ensure suitability for docking simulations.

Binding pockets were identified using cavity detection tools, followed by blind docking to explore ligand–receptor interactions across multiple potential binding sites. Post-docking analyses included interaction mapping to characterize hydrogen bonding, hydrophobic contacts, and π-mediated interactions. Binding free energy estimation using MM/GBSA and MM/PBSA methods was performed to evaluate relative interaction favorability. Structural analysis using RMSD was conducted to assess conformational deviation of the docked complexes. This workflow enabled the systematic evaluation of ligand–receptor interactions and facilitated the identification of phytochemicals with potential multi-target binding affinity toward components of the AXL–TGF-β signaling axis.

Structural characterization of the selected metastasis-associated receptors revealed distinct variations in binding pocket architecture, residue composition, and physicochemical properties that influenced ligand accommodation and interaction behavior (Table 4). Pocket volumes ranged from 630 ų (TGF-βRII) to 3503 ų (TGF-β1), indicating substantial variability in cavity size and flexibility. The AXL kinase domain (5U6B) exhibited a well-defined ATP-binding cleft enriched with acidic residues (Asp70, Glu74, Asp179) and aromatic residues (Phe36, Tyr187), forming a favorable electrostatic and hydrophobic environment for ligand binding. In contrast, TGF-β1 presented the largest and most aromatic-rich cavity, characterized by residues such as Trp30, Phe54, and Tyr69, supporting π–π interactions and preferential binding of aromatic ligands. TGF-βRII displayed a comparatively compact, amphipathic pocket dominated by Asp119, Arg163, and Tyr185, suggesting selective ligand accommodation. Meanwhile, TGF-βRI kinase exhibited a balanced donor–acceptor distribution (Arg80, Asp84, Tyr401), supporting stable ligand anchoring within its catalytic domain.

Overall, the presence of conserved acidic, basic, and aromatic residues across all receptors established a structurally compatible framework for multi-target ligand binding within the AXL–TGF-β signaling axis.

Molecular docking analysis demonstrated that phytochemicals derived from both Garuga pinnata and Mitragyna speciosa exhibited favorable binding affinities across all selected receptors, with clear differences observed based on ligand class and receptor topology (Table 5). Among Garuga pinnata compounds, 21-hydroxydammar-24-en-3-one showed consistently strong binding across multiple targets, with the highest affinity observed against TGF-βRII (ΔG = −16.8 kcal/mol) and AXL kinase (ΔG = −16.3 kcal/mol). Other diterpenoids, including Garuganin-V1_Et and Garuganin-V1_Me, also demonstrated favorable binding (ΔG ≈ −14.5 to −13.5 kcal/mol), indicating broad receptor compatibility. For Mitragyna speciosa, alkaloid derivatives such as 7-acetoxymitragynine and mitraphylline-26 exhibited strong binding interactions, particularly against AXL kinase (ΔG = −12.7 kcal/mol) and TGF-βRII (ΔG = −14.6 kcal/mol). These compounds showed consistent affinity across kinase and receptor targets, suggesting potential multi-target engagement. RMSD values for most ligand–receptor complexes remained within a low deviation range, indicating stable docking conformations and minimal structural perturbation within the predicted binding pockets.

Post-docking visualization of the receptor–ligand complexes (Figures 4 and 5) revealed that the selected phytochemicals adopted well-defined orientations within the binding pockets of AXL kinase and TGF-β receptors, demonstrating consistent spatial accommodation and alignment with key functional regions. In AXL kinase, ligands were predominantly positioned within the ATP-binding cleft, exhibiting a stable orientation along the catalytic groove, where their hydrophobic cores aligned with aromatic residues while polar functional groups extended toward solvent-accessible regions, facilitating hydrogen bonding. Similarly, in TGF-β receptors, ligands occupied central cavity regions with orientations that maximized surface complementarity, allowing simultaneous engagement of polar and hydrophobic residues within the pocket. The visual representations further indicated that high-affinity phytochemicals maintained compact and well-fitted conformations, with minimal steric clashes and favorable pocket occupancy. Diterpenoid scaffolds tended to occupy deeper regions of larger cavities, whereas indole alkaloids aligned along aromatic-rich surfaces, supporting π-mediated interactions. Overall, the observed ligand orientations and binding site conformations demanded (or demonstrated) a high degree of geometric compatibility between phytochemicals and receptor pockets, reinforcing their ability to effectively engage multiple targets within the AXL–TGF-β signaling axis.

Binding free energy calculations using MM/GBSA and MM/PBSA were performed to refine docking-based rankings and to quantify residue-level contributions to ligand stabilization (Table 6). Across targets, the calculated ΔG values were consistently negative, indicating energetically favorable ligand–receptor associations and supporting the trends observed in molecular docking. Among Mitragyna speciosa alkaloids, 7-acetoxymitragynine exhibited the most favorable binding against AXL kinase (ΔG ≈ −45.2 kcal/mol), characterized by a combination of hydrogen bonding and electrostatic interactions within the catalytic domain. Mitraphylline-26 also demonstrated strong stabilization, primarily driven by cation–π and electrostatic contributions within receptor cavities. For Garuga pinnata, diterpenoid scaffolds showed consistent energetic stabilization across receptors, with key contributions from hydrogen bonding and hydrophobic interactions involving residues such as Asp161, Asp179, and Phe36 within the kinase-associated regions. Importantly, these energy estimates were derived from energy-minimized static complexes and therefore represented relative binding favorability rather than absolute thermodynamic parameters. Within this context, the MM/GBSA and MM/PBSA results served to corroborate docking-derived trends and supported the prioritization of phytochemicals exhibiting stable and energetically favorable interaction profiles.

Structural consistency of the top-ranked phytochemical–receptor complexes was evaluated using RMSD analysis derived from docking-based energy-minimized structures. As summarized in Table 7 and illustrated in Figure 6, most ligand–receptor complexes exhibited low RMSD values, generally below 3.0 Å, indicating limited conformational deviation and consistent ligand accommodation within the predicted receptor binding pockets. These findings supported the structural reliability of the docking-derived binding geometries and suggested favorable spatial orientation of the phytochemicals within metastasis-associated receptor cavities.

Among the evaluated compounds, 21-hydroxydammar-24-en-3-one from Garuga pinnata demonstrated the lowest RMSD value against TGF-βRII (3DKC) (0.586 Å), reflecting highly consistent pocket accommodation and favorable structural stabilization within the compact receptor cavity (Table 7; Figure 5, 6). Similarly, the indole alkaloids 7-acetoxymitragynine and 7-hydroxymitragynine from Mitragyna speciosa exhibited low RMSD values against AXL kinase (5U6B), supporting favorable ligand orientation and structurally coherent interaction geometry within the catalytic region. Other phytochemicals, including Garuganin-V1_Et, Garuganin-V1_Me, and Mitraphylline-26, also maintained relatively low conformational deviation, suggesting consistent docking conformations across multiple receptor environments. Comparative RMSD analysis further indicated scaffold-dependent accommodation patterns, in which diterpenoid compounds preferentially aligned within deeper hydrophobic cavities, while indole alkaloids demonstrated effective adaptation within aromatic-rich binding regions, as represented in Figure 5, 6.

Importantly, the structural consistency observed through RMSD evaluation correlated well with the molecular docking, interaction mapping, and MM/GBSA–MM/PBSA analyses, collectively supporting the structural compatibility of these phytochemicals with key components of the AXL–TGF-β signaling axis. Within the limitations of a static computational framework, these findings provided additional structural validation for the predicted ligand–receptor interactions and supported the prioritization of selected ethnopharmacological scaffolds for future experimental investigation.

Collectively, the integrated computational analyses demonstrated that phytochemicals derived from Garuga pinnata and Mitragyna speciosa exhibited consistent binding trends, favorable interaction profiles, and structurally coherent docking conformations across multiple metastasis-associated targets. The convergence of docking affinities, residue-level interaction mapping, and MM/GBSA–MM/PBSA energy estimates indicated that these compounds were capable of engaging conserved functional regions within both AXL kinase and TGF-β receptor systems.

Notably, the complementary binding behavior observed between diterpenoid and indole alkaloid scaffolds suggested a scaffold-dependent multi-target interaction pattern, enabling effective accommodation within diverse receptor pocket environments. Within the constraints of a static computational framework, these findings supported the prioritization of selected phytochemicals as candidate multi-target ligands with potential relevance to the modulation of metastatic signaling pathways.