Molecular docking and dynamics analysis of selected phytocompounds against multi-targeted hepatocellular carcinoma

Abstract

Introduction: Hepatocellular carcinoma (HCC) is the fourth leading cause of cancer-related deaths worldwide, with a five-year survival rate of only 19%. Although Sorafenib is the primary systemic therapy, its limited efficacy and complex interactions with signaling pathways highlight the need for multi-target drugs.
Methods: This study evaluates the anti-cancer properties of selected phytochemicals against six key HCC target proteins, utilizing sorafenib tosylate as a positive control. Molecular docking was performed to evaluate binding affinities and interactions, and ADME/T predictions were generated to estimate drug-like properties. The topranking candidates were further evaluated using 100-ns molecular dynamics simulations to analyze conformational stability, protein-ligand interactions, and residue mobility.
Results: Silymarin (SA) emerged as the most effective compound, demonstrating greater predicted inhibitory activity than Sorafenib. SA showed high binding affinity for target proteins 6HH1 (−9.9 kcal/mol) and 1CM8 (−9.6 kcal/mol). Molecular dynamics simulations also revealed increased stability of the SA-protein complexes, particularly for the 1CM8-SA complex, which maintained high conformational stability. The root-mean-square deviation (RMSD) value was found to be around 2.1 Å, and the root-mean-square fluctuation (RMSF) values were
below 3 Å, indicating lower protein flexibility compared to both the native and sorafenib-bound complexes.
Conclusion: These computational findings provide a strong theoretical basis for Silymarin’s efficacy as a highly potent, multi-targeted therapeutic agent against HCC. The improved stability and binding properties of Silymarin compared with Sorafenib provide a strong rationale for advancing this compound into preclinical and clinical studies.

Introduction

Cancer is a significant global health concern and a leading cause of death. According to a WHO 2020 report, cancer is the primary or secondary cause of death in 112 countries, accounting for approximately 10 million deaths annually. Among these, liver cancer alone caused an estimated 8,30,000 deaths in 2020. Additionally, nearly 4,00,000 children are diagnosed with cancer each year. Hepatocellular carcinoma (HCC) is the most prevalent type of liver cancer, representing approximately 75% of cases, and is the fourth leading cause of cancer-related mortality, with a global survival rate of only 19%. Although the incidence of HCC has declined in some high-risk areas, it remains a significant health issue in many regions (1).

Primary liver cancer encompasses three main subtypes: HCC, ICC, and combined cHCC-CC, which account for approximately 75-85%, 10-15%, and 1-4.7% of cases, respectively. ICC is a malignant tumor that arises from the biliary epithelium and is the second most common primary liver cancer. Combined cHCC-CC is a rare form of liver cancer, characterized by features of both hepatocytes and bile duct epithelial cells, indicating bidirectional differentiation (2,3).

Plant-derived medicinal products have been pivotal in cancer treatment for centuries, with many modern drugs tracing their origins to natural sources. Ongoing research focuses on harnessing the potential of these plant-based compounds for cancer prevention, treatment, and management. Secondary metabolites, such as flavonoids, tannins, alkaloids, and terpenoids, exhibit significant anti-cancer properties. Approximately 60% of plant-derived products are effective lead molecules in drug discovery, often demonstrating synergistic effects. Some plant compounds inhibit tumor growth, induce cancer cell death, and enhance immune response, contributing to their therapeutic potential (4).

Paclitaxel (Taxol) is widely used to treat breast, ovarian, lung, and other cancers due to its ability to stabilize microtubules and prevent their disassembly, thereby halting cell division and inducing apoptosis in rapidly dividing cancer cells. Vincristine and vinblastine, vinca alkaloids, interfere with microtubule formation, disrupt mitosis, and trigger programmed cell death (5, 6). Camptothecin derivatives, such as irinotecan and topotecan, inhibit topoisomerase I, an enzyme essential for DNA replication, impeding cell proliferation in colon, ovarian, and lung cancers (7). Natural compounds such as curcumin, resveratrol, pterostilbene, and ginseng have also demonstrated anti-cancer effects comparable to those of synthetic agents. In this study, a multi-targeted approach involves concurrent inhibition of interconnected signaling pathways to prevent compensatory signaling that drives drug resistance. Six key target proteins, namely EGFR, RTKs, BRAF, MAPK, P53, and Caspase-9, were identified to have priority over other possible target proteins because they are the most critical and overlapping targets that control the survival, proliferation, and apoptosis of HCC cells. Although monotherapies have been shown to fail because HCC can evade inhibition by a single target protein, inhibiting this six-protein network is the key mechanism to overcome drug resistance. Based on a literature review and database availability, the identified phytocompounds are RA, SA, HMF, LU, and IC. These compounds will be screened to assess their anti-cancer properties, particularly their ability to target multiple signaling pathways involved in HCC progression. The chosen phytocompounds are rich in phenolic and flavonoid compounds, and several studies have examined their hepatoprotective properties and therapeutic potential.

EGFR and RTKs are overexpressed in more than 60% of HCC patients. EGFR is an integral molecule that contributes to both tumor neovascularisation and metastasis. The clinical validation of these targets is demonstrated by Sorafenib and Lenvatinib, which inhibit RTKs, as RTKs are a broad family of cell-surface receptors. Multi-kinase inhibitor targets RTKs, such as VEGFR and PDGFR, to halt tumor angiogenesis, which in turn slows the progression of HCC.

The MAPK and BRAF pathways, as well as the entire Ras/Raf/MEK/ERK signaling pathway, are constitutively activated. Although less frequently mutated in HCC than in melanoma, BRAF mutations, when present, act as an essential mechanism for HCC to escape standard chemotherapy.

The two primary gatekeepers of cell death are P53 and Caspase-9. p53 is the most commonly mutated gene in HCC (30-40% of patients) and, as a result, allows cells to evade the apoptotic death pathway. Caspase-9 initiates the intrinsic apoptotic pathway. Plant-based compounds that reconstitute the activity of Caspase-9 can induce cancer cells to undergo programmed cell death.

Approach for Multi-Target Inhibition

A significant limitation of monotherapy for HCC is the rapid development of compensatory signaling when one target (such as the epidermal growth factor receptor) is inhibited, leading to upregulation of another signaling pathway (e.g., mitogen-activated protein kinase). The concurrent inhibition of multiple targets within the oncogene network through multi-target inhibition confers a synergistic effect. It provides a greater “genetic barrier” to resistance since the likelihood that a tumor cell would develop simultaneous bypass mutations across multiple tumorigenic signaling pathways is very low.

Toxicity and Safety Profiles

The concurrent suppression of multiple signal pathways, as a therapeutic measure for cancer treatment, raises the point of concern over the possibility of toxicity to normal, healthy hepatocytes. However, the secondary metabolites of plants, such as flavonoids, tannins, and phenolic compounds, exhibit selective toxicity as hormetic agents, causing significant metabolic stress to cancer cells and inducing antioxidants to protect normal cells. The reduced dosage of multi-targeted compounds can be effective at much lower concentrations than fully targeted synthetic inhibitors, as they act together, thereby reducing systemic toxicity and the many “off-target“ side effects associated with aggressive chemotherapeutics.

This study employs an in-silico approach to evaluate the binding affinities and pharmacological effects of various plant compounds for the previously discussed HCC target protein. By integrating molecular docking and dynamics alongside ADMET profiling, we aim to identify a potent drug that can disrupt the robust signaling pathway associated with HCC. This may contribute to the creation of more effective treatment options.

Materials and methods

The 3D structure of the proteins, phosphorylated map kinase P38-GAMMA (PDB ID: 1CM8), Dimeric Caspase-9 (PDB ID: 2AR9), EGFR (PDB ID: 4LQM), BRAF in Complex with RAF265 (PDB ID: 5CT7), p53 cancer mutant Y220C (PDB ID: 5G4N), and RTK (PDB ID: 6HH1), were obtained from the PDB online. The validity of the three-dimensional structures of the proteins was assessed using the Swiss-Model server online and docked using AutoDock Vina 4. The target proteins and their functions are listed in Supplementary Table 1. The grid box parameters for each protein were meticulously calibrated to encompass the active sites, ensuring that the grid dimensions and center coordinates accurately covered the critical residues (8,9) implicated in binding interactions (Supplementary Table 6). The docking results were visualized by using Maestro v14.5 from the Schrödinger suite (10).

Structural and chemical details for a library of five phytochemicals were retrieved from the Ethnobotanical and Dr. Duke’s Phytochemical databases online. The 3D conformer of the ligands-SA (PubChem Id: 5213), HMF (PubChem Id: 150893), LU (PubChem Id: 5280445), RA (PubChem ID: 5281792), IC (PubChem ID: 5318528), and Sorafenib (PubChem ID: 406563) the 3D conformer of the ligands were obtained from the NCBI PubChem database in SDF format Online (Supplementary table 2) and were subsequently categorized based on their toxicity class. Further, the phytocompounds were energy-minimized and optimized to a low-energy state, and their flexibility was ensured using Chimera 1.18 and Avogadro, and they were converted to the Sybyl Mol2 format Online.

Molecular Dynamics (MD) simulation

MD simulation is a computational tool for modeling protein-ligand interactions, examining atomic and molecular vibrations, and tracking conformational changes over a specified period (11). This study provides insights into the stability of protein-ligand complexes within a solvent system. MD simulations for the two highest-ranking proteins, namely MAP kinase P38 gamma (PDB ID: 1CM8) and RTK (PDB ID: 6HH1), were conducted using the Schrödinger suite of Maestro and Desmond v6.5, developed by D.E. Shaw Research. The MD simulation used an orthorhombic simulation box measuring 10 Å × 10 Å × 10 Å, determined using the buffer-size approach, and the box volume was subsequently minimized. The MD simulations were performed in three separate phases. The first phase was the system preparation, energy minimization, and a 100 ns simulation, which provided sufficient time to observe the complex’s equilibration and stability. The second phase, the solvated system with SPC water, was neutralized with 0.15 M NaCl and pre-relaxed. The simulation was then carried out in the NPT ensemble at 300 K and 1.01325 bar. In the third phase, the final complex with a loaded trajectory was evaluated using a 10-ps block length for averaging. The data was generated in a time frame of a 100 ns simulation, and graphs are plotted using total energy (kcal/mol), potential energy (kcal/mol), temperature (K), pressure (bar), and volume (A³) parameters. The stability of the complex’s structure was determined by calculating the complex RMSD, defined as the average change in the displacement of certain atoms in the defined frame relative to the reference frame, using the formula (12).

Where N is the number of atoms in the atom selection.

tref is the reference time, and r’ is the position of the selected atoms in frame X

Results

Evaluation of Pharmacokinetic and Pharmacological Properties

ADME/T analysis was performed to predict pharmacokinetic and pharmacodynamic properties (Fig. 1). ADMET properties included molecular weight, dipole moment, number of H-bond donors and acceptors, log P values, violation of the rule of five, and oral absorption (13) were listed in (supplementary table 4). Toxicity analysis of shortlisted molecules was performed using the PROTOX-III web tool Online (14) The parameters used for evaluating the compounds were cytotoxicity, carcinogenicity, hepatotoxicity, immunotoxicity, mutagenicity, predicted LD₅₀, and toxicity (oral toxicity) were listed in (supplementary Table 5) out of all the ligands RA and Sorafenib fell in the toxicity class 5 (3000 < LD₅₀ ≤ 5000 mg/kg) indicating their possibility of being toxic if swallowed and all other ligands fell in the toxicity class 4.

FIGURE 1 -. Bioavailability Radar plot of each shortlisted ligand.

Molecular docking

The present study conducted site-specific docking analysis to evaluate the binding affinities of six key targeted proteins associated with HCC against a selection of phytocompounds, including Sorafenib. The findings are summarized in supplementary Figures (1-6), demonstrating the docking scores and differential interactions between all phytocompounds and the target proteins. The binding affinity scores of each ligand and protein are presented in supplementary Table 7. SA showed the strongest binding affinity among the compounds tested, especially in targets 6HH1 and 1CM8. It demonstrated a significant ability for hydrogen bonding, which contributed to its overall stability within the binding sites of these proteins (supplementary Table 8). IC also exhibited moderate binding affinity for all the targeted proteins, characterized by substantial hydrogen bonding and hydrophobic interactions. In contrast, LU exhibited moderate binding affinity, accompanied by significant hydrogen bonding and other interactions. Conversely, HMF and RA exhibited the lowest binding affinities, while maintaining moderate hydrophobic and hydrogen bonding interactions. Overall, Sorafenib revealed its superior binding affinity across all targets, particularly 6HH1, 1CM8 and 5CT7 as detailed in (supplementary Tables 9-14). Overall, the phytocompound SA and Sorafenib showed stronger interactions with the target proteins 6HH1 and 1CM8, prompting further analysis by molecular dynamics to examine interaction stability, docking poses, and binding states.

Molecular dynamics (MD)

RMSD and RMSF plots

RMSD values over time are pivotal for assessing protein compactness during ligand binding and serve as indicators of the overall stability of the docked complexes (15). Moreover, it reflects the overall stability of the docked complexes (16). The minimal fluctuation observed in the RMSD values indicated enhanced stability of both the protein and the protein-ligand complex. The RMSD plots for MAP kinase and its complex over a 100-ns time scale are presented in Figures 2A and B. During the initial 0-15-ns time scale, the RMSD of 1CM8 fluctuated within a range of 1.7 to 2.8 Å during the 10-ns period. After 25 ns, it became stable at 3.2-4.7 Å until the end of the simulation. Additionally, an increase in RMSD was observed between 60 and 70 ns, fluctuating from 3.5 to 4.1 Å before restabilizing. In contrast, the RMSD of the 1CM8-SA complex (Fig. 2B) initially fluctuated within a range of 1.5-2.0 Å during the period of 20 ns. After 20 ns, it remained highly stable at ~2.1 Å, followed by minor fluctuations between 3.8 and 4.2 Å until 100 ns. In Sorafenib (Fig. 4A), the protein and ligand complex fluctuated within a range of 3.5 Å during the 10-30 ns period, after which it became stable.

The RMSD values of 6HH1 and 6HH1-SA complex are shown in (Fig. 3A and B), respectively. The RMSD of 6HH1 in (Fig. 3A) fluctuated from 1.3 to 2.4 Å during the first 10 ns. Then it stabilized at ~2.3 Å between 20 and 60 ns, increased to 3.3-3.7 Å until 70 ns, and finally stabilized at ~3.1 Å until the end of the simulation. Conversely, the RMSD of the 6HH1-SA complex (Fig. 3B) decreased initially, stabilizing in the 2.1-2.7 Å range up to 40 ns. After 40 ns, the complex fluctuated between 3.2 and 3.5 Å until 65 ns, then stabilized at 3.8-4.1 Å until 100 ns. In the sorafenib complex (Fig. 4B), structural fluctuations occurred within a 3.5 Å range during the first 10-30 ns, after which it remained stable.

The RMSF is a numerical measure of individual-residue flexibility; lower coordinate fluctuations indicate greater stability. The RMSF of the 1CM8 and 1CM8-SA complexes are illustrated in (Figs 2C and D). The active site residues MET109 and PHE111 exhibited significant hydrogen bonding throughout the analysis, whereas PRO156 and LEU159 demonstrated notable hydrogen bonding and hydrophobic interactions. As depicted in Figure 2D, residues 30-65 and 100-160 engaged in multiple interactions with the SA ligand (indicated by green lines representing bond formation). In contrast, Sorafenib showed greater residual fluctuation (Fig. 4C) due to fewer hydrogen bonds, with only ASP115 maintaining consistent binding interactions throughout the run. The RMSF data show that the 1CM8-SA complex was more stable than 1CM8 alone. The green lines indicate the amino acid residues that participate in the bond formation within a molecular fragment. As depicted in (Fig. 2D), residues ranked from 30 to 65 and 100 to 160 engaged in multiple interactions with the ligand.

FIGURE 2 -. RMSD of 1CM8 (A) and 1CM8-SA (B), RMSF of 1CM8 (C), 1CM8-SA, (D) at 100 ns.

FIGURE 3 -. RMSD of 6HH1 (A) and 6HH1-SA (B), RMSF of 6HH1, (C) 6HH1-SA, (D) at 100 ns.

FIGURE 4 -. RMSD of 1CM8—Sorafenib (A) and 6HH1-Sorafenib (B), RMSF of 1CM8—Sorafenib (C) and 6HH1—Sorafenib (D) at 100 ns.

FIGURE 5 -. Heat map and protein—ligand contact histogram of 1CM8-SA (A), and 6HH1-SA (B) at 100 ns.

RMSF data for the 6HH1 and 6HH1-SA complex are presented in (Figs 3C and D). In the 6HH1 structure (Fig. 3C), this complex shows structural fluctuations, particularly LEU595, VAL603, and VAL644, with variations reaching up to 5.8 Å. This complexity shows structural fluctuations in key interacting residues VAL668, TYR672, ASP810, and ILE789, which were notably restricted to a much narrower range of 1.2 to 1.8 Å and subsequently stabilized at 1.2-1.8 Å. Furthermore, the residues spanning from 502-630 and 670-800 regions exhibited robust interactions that contributed to the formation of a stable secondary structure. The sorafenib complex showed weaker interactions due to the lack of residue interactions (GLU640 and ASP810) and fewer H-bond interactions (Fig. 4D).

Protein-ligand contact histogram

The interaction between 6HH1 and SA was predicted to involve a combination of hydrogen bonds, hydrophobic interactions, water-mediated interactions and ionic interactions. Generally, hydrogen bonds and water bridges indicate that the amino acid residues ALA621, LYS623, ILE789, and HIS790 showed hydrophobic interactions. Concurrently, VAL668 and CYS673 formed robust hydrogen bonds, while TYR672, CYS788, and ASP810 participated in water-bridge interactions. The detailed protein-ligand contacts for 6HH1-SA are demonstrated in (Fig. 5B). In contrast, the sorafenib-6HH1 complex (Fig. 6B) showed only moderate H-bond interactions with a few residues, specifically GLU640 and ASP810, throughout the simulation.

Furthermore, the protein-ligand contact of 1CM8-SA (Fig. 5A) exhibited stable and strong hydrogen bonding with MET109, PHE111, MET112, and ASP171. Moderate hydrophobic interactions were observed with VAL33, ALA40, VAL41, ALA54, and ALA160, while water bridges were formed with THR114, ASP115, PRO156, and LEU159. Compared to the 6HH1-SA complex, 1CM8-SA reflected a highly stable interaction network characterized by dense hydrogen bonding and hydrophobic contacts. In contrast, the 1CM8-sorafenib complex showed weaker hydrogen bonding, particularly relying on a single residue, ASP115, which interacted with and stabilized the protein until 100 ns (Fig. 6A).

The heat map of the 1CM8-SA complex showed a high number of contacts with MET109 and PHE111, indicated by darker bands (Fig. 5A); however, these specific contacts frequently broke and reformulated between the 10-30 ns and 40-50 ns intervals. Similarly, the heat map for the 6HH1-SA complex (Fig. 5B) showed distinct dark bands for VAL668 and CYS673, indicating strong, consistent binding affinity. This was accompanied by more dynamic contact behavior from ILE789, ALA621, HIS790, and ASP810, which were observed to lose and gain contacts over shorter timescales (20-45 ns and 80-90 ns).

Discussion

In silico methodologies are crucial in identifying potential binding sites and facilitating the discovery of novel molecules that can interact with known sites. Virtual screening and active site-predicted docking are frequently used to develop new pharmacological agents. However, a limited body of research has focused on the in silico molecular docking of phytocompounds in the context of HCC (17,18). In this current research, site-specific docking was used to investigate receptor-ligand interactions. Five unique phytocompounds (SA, HMF, RA, LU, IC), in addition to sorafenib tosylate as a positive control, were tested against six key HCC targets: EGFR, Caspase-9, BRAF, RTK, P53, and MAPK. The docking results indicate that SA exhibits more favorable binding stability than other phytocompounds and Sorafenib, as evidenced by stronger hydrogen-bonding interactions with each targeted protein. Although Sorafenib exhibited higher overall affinity, stabilization was compromised owing to reduced hydrogen bonding and hydrophobic interactions. Overall, the study suggests that SA could be a potent multi-target therapeutic compound for HCC. Meanwhile, IC and RA also exhibited moderate binding affinity with all listed targets.

SA is a diverse combination of flavonoids extracted from the seeds of the Silybum marianum plant, commonly referred to as milk thistle. The main components of SA are silybin and its isomers, silicristin and silidianin. Studies suggest that SA possesses notable antioxidant capabilities and helps stabilize membranes. Additionally, it protects various tissues and is induced by chemicals, indicating its potential as a hepatoprotective agent (19).

FIGURE 6 -. Heat map and protein-ligand contact histogram of 1CM8- Sorafenib (A), and 6HH1—Sorafenib (B) at 100 ns.

The varied effects of SA on different cellular pathways underscore its potential as a treatment for HCC. This compound has been shown to induce apoptosis, promote angiogenesis in cancer cells, and inhibit tumor development. Research studies conducted before clinical trials have indicated that SA protects hepatocytes from toxins and oxidative stress. Additionally, clinical studies have revealed its potential effectiveness in treating conditions like hepatitis, cirrhosis, and non-alcoholic fatty liver disease (20,21).

SA has been identified as a potent inhibitor of p38 MAP kinase, demonstrating anti-cancer, anti-diabetic, and antitumor effects. Evidence suggests that it has the potential to serve as an anti-cancer therapeutic agent. Silibinin inhibits p38 MAPK phosphorylation during the activation of macrophages and subsequent nitric oxide production in RAW 264.7 cells. The findings indicate that Silibinin effectively inhibits macrophage activation by obstructing the p38 MAPK signaling pathway (22).

Silybin has been shown to have cardioprotective properties (23). The study’s intricate network of multi-target interactions and numerous pathways through which silybin exerts its effects in managing MI. The findings suggest that silybin mitigates oxidative stress by increasing antioxidant levels via the PI3K-Akt/HIF-1 signaling pathway.

Time-Dependent Variation and Multi-target Activity

The time-dependent cytotoxicity of LU, RA, and SA varies across cell lines, with differing IC₅₀ values. LU exhibited IC₅₀ values ranging from 15.1 to 103.1 µM across various cell lines, including SMMC-7721, HepG2, and Bel7402, and its mechanism of action involves inducing apoptosis, autophagy, or metabolic disturbance by interacting with target proteins such as Caspase-8, LC3B, p53, and DPD (24,25) The cytotoxicity of RA shows a significant effect against HepG2 cells, with IC₅₀ values of 14 µM and 33 μg/mL, via mechanisms that inhibit glycolysis and the metastasis pathway, induce apoptosis via Fyn or PI3K/AKT signaling pathways, and affect GLUT-1, matrix metalloproteinases (MMPs), and various caspase signaling pathways (26,27).

In contrast, SA is identified as a potent anti-proliferative agent against HepG2 and Hep3B cells under hypoxic conditions, with IC₅₀ values ranging from 45.71 to 75.13 µmol/L (28,29) Its mechanism of action involves inhibiting the HIF-1α/VEGF pathway, which is crucial for angiogenesis, and inducing cell cycle arrest by modulating key proteins such as Cytochrome c, p53, and cyclin D1. The various phytocompound studies show the target metabolic reprogramming, autophagic survival mechanisms, and hypoxia-inducible pathways, suggesting these phytocompounds could be promising multi-target compounds for the development of novel therapeutics for HCC.

In addition, results from in silico studies indicate strong binding of these compounds to key HCC targets, including EGFR, MAPK, and RTKs; however, these in silico findings require experimental confirmation. In summary, the current findings of time-dependent inhibition of HCC signaling pathways provide evidence for the efficacy of SA, RA, and LU as potential multi-target anti-cancer agents. In contrast, experimental IC₅₀ values for HMF and IC in HCC cell lines are not yet available.

Numerous studies have shown that SA, known for its hepatoprotective properties, has demonstrated considerable therapeutic efficacy in managing HCC. Nonetheless, the precise molecular mechanisms underlying SA’s therapeutic effects remain only partially elucidated. Therefore, examining the interplay between the active constituents of SA and key molecular targets is essential to clarify its beneficial effects and enhance its therapeutic potential. In this pivotal study, molecular docking and MD simulations were used to analyze how SA affects HCC (20,30).

The RMSD and RMSF values for the protein and complex indicate that 1CM8 was more stable throughout the simulation, from 20 ns, with an RMSD range of 2.1-3.1 Å. Also, interactions among residues stabilize the protein’s fold throughout the simulation more than in 6HH1. 1CM8 shows greater stabilization relative to its unbound state than 6HH1. The binding of SA to 6HH1 successfully brought the complex to a local minimum energy state, ensuring structural integrity. These in silico molecular docking and dynamics analyses strongly suggest that the phytocompound SA effectively inhibits both 1CM8 and 6HH1, showing a promising candidate for modulating the target-mediated progression of HCC. Most importantly, extensive hydrogen-bond networks bind ligands to specific amino acid residues in target areas. By elucidating these molecular interactions, this study aims to inform the design of targeted treatments and to highlight SA’s potential to enhance current HCC therapies.

The 100-ns timeframe utilized for the MD simulations was used to validate the sustained stability of the docked complexes. While shorter simulations (< 50 ns) often capture only the initial relaxation of the protein-ligand system, the 100-ns simulation allowed thorough observation of dynamic equilibrium. The extended timeframes of RMSD and RMSF values confirmed the stabilization of the 1CM8-SA and 6HH1-SA complexes and indicated a sustained thermodynamic state rather than a transient local energy minimum. This timeframe duration was sufficient to identify the breaking and reforming of hydrogen bonds, reinforcing the reliability of the observed multi-point anchoring systems.

The dynamic stability of SA complexes, in contrast to Sorafenib, is quantitatively influenced by the density of their non-covalent interaction networks. SA establishes a highly redundant, multi-point anchoring system within the binding pockets of both targeted proteins. In the case of structure 1CM8, SA maintains a robust network comprising at least 13 distinct residue interactions, including four strong hydrogen bonds. In comparison, Sorafenib relies heavily on a single hydrogen bond interaction with ASP115. In the 6HH1 structure, the binding of SA to nine critical residues varies in binding types, whereas Sorafenib binds to only two residues. This extensive binding mode greatly reduces atomic fluctuations, thereby limiting the RMSF of critical residues like VAL668 and TYR672 in 6HH1 to the range of 1.2-1.8 Å. This results in greater structural compactness in SA, enabling it to act as an anchor to prevent conformational drifts that would have occurred if only single-point binding had occurred. The overall results, summarised from various in silico analyses, suggested that SA strongly inhibits MAP kinases and RTKs.

Conclusion

This study employed a combined In-silico method to explore the therapeutic efficacy of selected phyto-compounds and the standard drug (sorafenib tosylate) in treating HCC. We identified five selected phytocompounds from an ethnobotanical database and six common HCC targets. Among these, SA showed a potent therapeutic effect in molecular docking studies involving MAP kinase and RTK in combination with SA and Sorafenib, elucidating the binding modes and residue interactions between the target protein and ligands. SA exhibited a higher binding affinity for MAPK and RTK proteins than Sorafenib, as indicated by a lower free energy of binding. This binding interaction may inhibit phosphorylation events and disrupt essential signaling pathways critical for tumor growth, invasion, and angiogenesis. Therefore, SA can impede cancer progression by blocking cell proliferation, migration, and survival signals that are crucial for HCC development, thereby demonstrating its potential as a potent anti-cancer compound.

In addition, MD revealed that the RMSD and RMSF of 1CM8 and 6HH1 exhibited reduced deviation and fluctuation in SA compared to Sorafenib, demonstrating stability through hydrogen bonding and protein-ligand interactions involving interacting residues throughout the final run. These results provide a theoretical foundation for the possible therapeutic effects of SA in HCC. However, given the computational nature of this research, further experimental and clinical studies are imperative to validate these findings and assess the practical application of SA as a multi-target treatment for HCC.

Acknowledgments

All authors thank Sathyabama Institute of Science and Technology for their constant encouragement and support.

Other information

This article includes supplementary materials

Corresponding author:

Swarnalatha Yanamandala

email: swarnalatha.y@rajalakshmi.edu.in

Disclosures

Conflicts of interest: The authors declare no conflict of interest.

Financial support: This work is supported by the DST-SEED, Ministry of India, through a grant to Swarnalatha Yanamandala (DST sanction letter No. DST/SEED/TSP/STI/2021/466).

Data availability statement: The data are available from the corresponding author upon reasonable request.

Authors’ contributions: NC: Contributed to conceptualization and methodology and wrote the main manuscript. SY: oversaw supervision, investigation, and editing. SC, SS, SK, PX, KV, KAR & SBY prepared the figures and revised the manuscript. All authors reviewed the manuscript.