Introduction

Lung cancer remains an aggressive malignancy with high mortality and incidence rate because of its aggressive characteristics1. Among the common subtypes of lung cancer, non-small cell lung cancer (NSCLC) is the primary subtype, accounting for 85% of the cases2. The standard therapy for NSCLC includes surgery, radiation, chemotherapy, targeted therapy, and immunotherapy3. In the past decade, despite advances in therapeutic approaches, the number of deaths has gradually increased due to NSCLC’s ability to metastasize, acquisition of chemotherapeutic resistance, and high recurrence rate; in addition, drug efficacy is limited by significant side effects4,5. Particularly, a five-year survival rate of only 7% was reported for patients with lung cancer at the metastatic stage1. Therefore, continuous drug discovery for lung cancer is necessary for improving the clinical outcome.

Plants are a versatile source of biologically active compounds. For example, the plants of the Aspidistra genus, discovered in Vietnam6, contain many active components, such as saponin, coumarin, and isoflavones7,8, with various demonstrated pharmacological activities, including antibacterial, antifungal, antitumor, and antiviral7,8,9,10. A previous study reported that Aspiletreins, mainly three steroidal saponins, exhibited cytotoxicity in breast, cervical, hepatocellular, gastric, and lung cancer cells8,11. These three steroidal saponins contain different sugar moiety side chains located at C3, suggesting distinct potencies (Table 1). The mechanistic investigation demonstrated that Aspiletrein A (AA) could suppress lung cancer metastasis by inhibiting protein kinase A (Akt) signaling11. However, the molecular targets of Aspiletreins B and C (AB and AC, respectively) have not been elucidated.

Table 1 Compound information.
Full size table

On the other hand, network pharmacology-based strategies, integrating complex systems of biology and the network analysis of multiple drug targets, have been established as a new paradigm for drug discovery12. This type of approach is a powerful tool for predicting the molecular targets of new chemical entities as well as disease pathways. In addition, molecular docking analysis reveals the potential intermolecular interactions between a target and compounds13. Furthermore, pharmacokinetics and pharmacodynamics were analyzed according to chemical properties and structures using the algorithm-based method14. These comprehensive approaches can help reduce the time and cost of preclinical experiments. In this study, we explored the molecular mechanism of Aspiletreins against NSCLC using network pharmacology and computation-based approaches.

Results

Identification of drug targets in non-small cell lung cancer

The experiment procedures are presented in Fig. 1. Targets in NSCLC were searched in the GeneCards Human database, DisGeNET, OMIM, and TTD using “non-small cell lung cancer” as a keyword. After removing duplication among these databases, a total of 6431 NSCLC targets were found based on the relevance scores (Fig. 2). The targets of each Aspiletrein were fetched from the Swiss Target Prediction. The predicted targets identified included 28 for AA, 32 for AB, and 30 for AC (Table 2). The 17 intersecting targets of Aspiletreins and NSCLC were in a Venn diagram (Fig. 2) and listed (Table S1). Furthermore, the compound-target network plot was constructed using Cytoscape 3.9.0 (Fig. 3). Among the targets, STAT3 and IL2 for AA and AB, GLRAs and TRPV1 for AB and AC, whereas TACR2 was a target only for AC.

Figure 1
figure 1

The framework of this study.

Full size image
Figure 2
figure 2

Venna diagram representing the overlapping of NSCLC targets (yellow) and compound targets (blue).

Full size image
Table 2 Lists of compound targets obtained from the Swiss Target Prediction.
Full size table
Figure 3
figure 3

Compound-target-NSCLC network constructed using Cytoscape v_3.9.0. The grey rectangles represent the compounds (AA, AB, and AC), and the white ovals represent the hub protein-target interaction. AA Aspiletrein A, AB Aspiletrein B, AC Aspiletrein C.

Full size image

Analysis of target protein–protein interaction (PPI) network

The 17 intersected targets, including VEGFA, FGF1, FGF2, HPSE, CDK1, HSP90AA1, ADRA2B, DRD2, CYP2D6, LGALS3, RORC, IL2, ADRA1A, STAT3, TRPV1, SLC6A2, TRPV1, and GLRA1, were imported to the STRING database, and the potential relationships among them were investigated. The PPI network with a confidence score of 0.9 was constructed (Fig. 4A). The STRING database analysis revealed that the average node degree, defined as the average number of interactions of a protein in a network, was 1.06 and the local clustering coefficient, defined as the wellness of the connected nodes in a network, was 0.312. The interactions between the targets comprised 17 nodes and 9 edges, with each edge representing the association between nodes. Cytoscape analysis showed that there was one main cluster associated in the PPI network (Fig. 4B). According to the scores of degree, closeness centrality, betweenness centrality, and clustering centrality, signal transducer and activator of transcription 3 (STAT3), heat shock protein HSP 90-alpha (HSP90AA1), vascular endothelial growth factor A (VEGFA), fibroblast growth factor-2 (FGF2), and interleukin-2 (IL-2) were identified as the top 5 intersecting targets of Aspiletreins and NSCLC interaction (Fig. 4C and Table S2).

Figure 4
figure 4

Protein–protein interaction (PPI) analysis. (A) Protein–protein interaction network of Aspiletreins and NSCLC targets obtained from STRING v_11.5 database. (B) The PPI network was constructed using the plug-in targets from the STRING database and imported into Cytoscape. (C) The top 5 targets in the PPI network as ranked using the cytoHubba plug in network analyzer. The higher degree value is represented by colors ranging from red to yellow.

Full size image

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis

GO enrichment analysis was conducted to analyze the impact of the targets in NSCLC. The classification of GO was set for three criteria, including biological processes (Fig. 5A), GO molecular functions (Fig. 5B), and GO subcellular localizations (Fig. 5C). First, the data with the candidate target and KEGG pathway analysis were interpreted using the R software with a ggplot2 plug-in, uncovering 245 biological processes, 9 molecular functions, and 3 subcellular localizations. Then, the top 10 biological processes of Aspiletreins were identified by sorting them by the degree of significance, including response to ethanol, regulation of protein modification process, regulation of phosphate metabolic process, positive regulation of signal transduction, positive regulation of protein phosphorylation, positive regulation of multicellular organismal process, positive regulation of kinase activity, positive regulation of intracellular signal transduction, positive regulation of catalytic activity, and cellular response to chemical stimulus (Fig. 5A). Molecular functions were mainly enriched in growth factor activity, chemoattractant activity, growth factor receptor binding, protein binding, receptor-ligand activity, molecular function regulator, sulfur compound binding, catecholamine binding, and alpha-adrenergic receptor activity (Fig. 5B). Finally, subcellular localization includes extracellular space, platelet-derived growth factor complex, and VEGF-A complex (Fig. 5C). Furthermore, the KEGG pathway suggested that the predicted targets were components of the pathways involved in oncogenesis, particularly in the PI3K-Akt, Rap1, Ras, and MAPK signaling pathways, and resistance to EGFR tyrosine kinase inhibitors (Fig. 5D, Fig. S1).

Figure 5
figure 5

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The biological processes (A), molecular functions (B), subcellular localization (C), and KEGG terms (D) distributed in the ordinate and the degree of enrichment were analyzed. The size of the dots represents the gene count. Adjusted P-value indicates the importance of enrichment in which the blue-to-red represents the high-to-low value of the enrichment.

Full size image

Compounds-target interaction analysis by molecular docking and molecular dynamic simulation

Among the 17 main targets of Aspiletreins in NSCLC, the top 5 potential targets, including STAT3, VEGFA, HSP90AA1, FGF2, and IL-2, were investigated for possible interaction with Aspiletreins. The Aspiletreins were molecularly docked using the PyRx Virtual Screening Tool, and their interactions with the highest affinity in each target were presented (Fig. 6, Fig. S2). The major binding interactions included hydrogen bonding and Van der Waals interactions. The binding energy score during docking indicates the affinity of a component for the target protein15. Here, all of the binding energy scores analyzed were less than 0, suggesting high-affinity interactions among the targets and the Aspiletreins (AA, AB, and AC) (Table. 3).

Figure 6
figure 6

Molecular docking between Aspiletreins and the top 5 targets that have the highest affinity. (A) The interaction between AA and STAT3 (PDB: 6NJS). (B) The interaction between AA and VEGFA (PDB: 4KZN). (C) The interaction between AB and HSP90AA1 (PDB: 4BQG). (D) The interaction between AB and FGF2 (PDB: 2FGF). (E) The interaction between AA and IL2 (PDB: 1M48). Hydrogen bonds were displayed in green, and Van der Waals interactions were displayed in light green. Alkyls were displayed in pink, and π–sigma was displayed in purple. AA Aspiletrein A, AB Aspiletrein B.

Full size image
Table 3 Binding energy among the compounds and the top 5 potential targets.
Full size table

Since our finding demonstrated that STAT3 exhibited the most potential target of compounds in lung cancer cells, molecular dynamic simulation between STAT3 and AA or AB was then performed. A RMSD plot generated by R-Studio was used to evaluate the residual deviations in the complexes. The data shows that both AA and AB have comparable consistent interaction with STAT3 over the time after their binding (Fig. 7, Videos S1 and S2). However, the ligand movement in the binding pocket side of STAT3 showed that AB was more stable due to its longer time in the aspect of equilibrium, which average RMSD of AB was 3.07 ± 0.49 Å, whereas that of AA was 6.87 ± 2.4 Å. Apart from RMSD plot, other factors such as the amount of hydrogen, hydrophobic, and Van der Waals interaction between ligand and protein are important to determine the efficacy of ligand–protein interaction16. In this case, AB has more significant hydrogen interaction (Table S3) compared to AA, in which hydrogen bonds play a major role in the stabilization of ligand–protein complexes17. These results suggest the stable behaviors of the complexes formed between AA or AB with STAT3.

Figure 7
figure 7

Molecular dynamic interaction between AA and AB with STAT3 during 0–10 ns. The RMSD for ligand movement plot for interaction of AA (red) or AB (green) with STAT3 complex.

Full size image

In vitro cytotoxicity and target validation

Cytotoxicity of the tested compounds was performed in NSCLC-H460 cells by MTT assay. The results demonstrated that IC50 of AA, AB, and AC in H460 cells were 13.25 ± 0.77, 6.82 ± 0.87, and 15.76 ± 0.17 μM respectively (Fig. S3A), indicating that AB has the most potent compound. Next, an in vitro experiment for validating the molecular target of the compounds was examined. Since STAT3 plays an important role in cancer survival and apoptosis, and our finding demonstrated that STAT3 was the most relevant core target (Fig. 4, Table S2), the effect of AB on STAT3 activity was then investigated. The data demonstrated that phosphorylated STAT3 (pSTAT, an active form) was significantly decreased in response to AB, whereas its total form was unchanged (Fig. S3B), confirming the molecular mechanism identified. This suggests that network pharmacology in combination with molecular docking approaches is a beneficial platform for the identification of new drugs.

Discussion

Despite the advances in cancer treatment, such as chemotherapy, radiotherapy, targeted therapy, and immunotherapy, lung cancer continues to be one of the main causes of cancer-related deaths worldwide and a major global health problem due to its low five-year overall survival rate and poor prognosis1. Natural compounds are a promising source of drug candidates for lung cancer treatment18; therefore, they are investigated for their potential in improving clinical overcome. Recently, our group discovered Aspiletreins A, B, and C, promising natural compounds from Aspidistra letreae. In addition, these steroidal saponins reportedly exhibited cytotoxicity in various cancer cell lines8,11. In this study, we have identified the molecular mechanism between the Aspiletreins and their predicted targets using the network pharmacology approach in combination with in silico molecular docking.

Network pharmacology-based assessments provide in-depth analyses on how various network factors interact with potential drug candidates19,20,21. This approach has emerged as a powerful tool for drug target identification based on a multidisciplinary concept integrating biology systems and polypharmacological models22. Aside from discovering new drugs and their targets, this tool enables the exploration of prospective target areas and the repurposing of existing drugs for diverse diseases23.

In this study, the potential targets of the compounds that overlapped with NSCLC targets were identified, including STAT3, VEGFA, HSP90AA1, FGF2, and IL2. These targets are highly relevant in the progression and metastasis in lung cancer24,25,26,27,28,29. STAT3, one of the main members of the STAT family, is responsible for the transcription of several genes related to oncogenesis, such as those encoding the proteins of the CDK family and pro-survival Bcl-2 family, VEGF, and epithelial-to-mesenchymal transition-related transcription factors SNAIL and SLUG25,30. Even though STAT3 expression was not significantly altered in tumor (Fig. S4A), the STAT3 mutation was strongly correlated to the overall survival rate of lung cancer patients (Fig. S4B), and hyperphosphorylation of STAT3 was extensively found in lung adenocarcinoma (Fig. S4G), suggesting that suppression of its phosphorylation provides promising therapeutic approach31,32. Meanwhile, VEGFA, a growth factor, functions particularly through VEGFR-1 and VEGFR-2 on endothelial cells. It has various modes of action, including forming new blood vessels and vascular networks, increasing vascular permeability, stimulating endothelial cell proliferation and migration, and preventing endothelial cell apoptosis26,33. On the other hand, HSP90AA1 is a gene that regulates the expression of heat shock protein 90α (HSP90α) in response to intracellular stress27. In addition to its intracellular roles, HSP90α acts as an inflammation-stimulated, secreted extracellular factor that promotes malignant cell motility and metastasis by activating NF-kB and STAT3 transcription programs27. In the case of growth factor FGF2, it possesses broad mitogenic functions under normal conditions, especially in embryonic development and tissue repairment. In contrast, in cancer, its gene becomes overexpressed, inducing uncontrolled proliferation and the metastasis of several malignant tumors29,34. In addition, IL2, a multifunctional glycoprotein of the interleukin family, is a growth factor involved in the immune system, promoting the growth of natural killer cells, B-cells, and T-cells35,36. IL2 also aids the immune system by improving the ability of specific white blood cells to identify and destroy cancer cells37. In lung cancer tissues, low expression of IL2 was associated with poor prognosis (Fig. S4C). Even though some of them, their expressions were not directly relevant to the overall survival rate of NSCLC patients (Fig. S4), they participate in cancer signaling as downstream molecules or ligands of growth factor receptors whose upstream regulators or receptors were aberrant expressed or overactivated in cancers38,39,40,41. Furthermore, targeting these core targets was widely reported as potential therapeutic strategy for cancer42,43,44. Taken together with previous study, the mechanisms of these molecular targets may inform on the mechanisms of Aspiletreins against cancer8,11.

These promising targets were investigated for their interactions with Aspiletreins using molecular docking analysis. Consistent with our hypothesis, the Aspiletreins bound to the targets with high affinity and likely inhibited the tumor-promoting functions of STAT3, VEGFA, HSP90AA1, and FGF2 or enhanced tumor-suppressing activity of IL2, or both. In addition, the cytotoxicity of these compounds was shown, which AB posed the most potent toxic to lung cancer cells, and STAT3 was validated as a potential target of this compound. Although advances in the field of biological systems and network pharmacology can aid rational drug design and development as well as target identification45,46, the information of new compounds including Aspiletreins are not completely available in several databases, the interactions between drugs and their targets need to be further elucidated, and pharmacokinetic profile and the safety and efficacy of the drugs need to be verified. This work, at least, full filled the scientific information of this compound suggesting their potential for anti-cancer drug research and development for clinical application.

Conclusion

Using network pharmacology and molecular docking-molecular dynamic approaches, this study systematically analyzes the pharmacological mechanism of Aspiletreins in the treatment of NSCLC. Aspiletreins suppress NSCLC by targeting various proteins, including STAT3, VEGFA, HSP90AA1, FGF2, and IL-2. The network pharmacology approach has significant advantages for uncovering the mechanism of Aspiletreins, and narrowing down the number of targets, thus minimizing the time and cost of preclinical investigations. Further in vitro and in vivo experiments will be conducted to verify the mode of action of Aspiletreins against NSCLC.

Methods

Compound database

The compounds information of Aspiletrein A (AA), Aspiletrein B (AB), and Aspiletrein C (AC) were identified using the J-GLOBAL Database (https://jglobal.jst.go.jp/en), while the log P information (Table 1) was obtained from the pkCSM tools (http://biosig.unimelb.edu.au/pkcsm/), an online pharmacokinetic prediction model.

Compounds-target identification

The canonical SMILE of the compounds was uploaded into the Swiss Target Prediction database (http://www.swisstargetprediction.ch/)47. The compound-target network was visualized using Cytoscape (3.9.0)48.

Non-small cell lung cancer target identification

The NSCLC-related targets were retrieved by searching for “non-small cell lung cancer” in the GeneCards database (https://www.genecards.org/), DisGeNET (https://www.disgenet.org/), OMIM (https://www.omim.org/), and TTD (http://db.idrblab.net/ttd/). The targets of NSCLC and compounds (AA, AB, and AC) were integrated using VENNY (https://bioinfogp.cnb.csic.es/tools/venny_old/index.html) and the intersected targets were presented in a Venn diagram.

Construction of protein-PPI network

A PPI network was constructed by uploading the genes to the STRING v_11.0 database (https://string-db.org/)49. The settings for building the PPI network were established in accordance with the “Homo sapiens” model, and the confidence of the interaction between the targets was set at 0.950. The network nodes represented proteins, and the edges reflected the protein–protein interactions. For the core target protein identification, the Cytoscape data was sorted by using these parameters, including degree, closeness centrality, betweenness centrality, and clustering coefficient.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis

Gene Ontology (GO)51 and Kyoto Encyclopedia of Genes (KEGG)52 were analyzed by plotting the data analysis obtained from STRING 11.5 database as a scatterplot to display the relationship between two continuous variables of gene count and adjusted P-value. The R software with ggplot253 was applied to construct a bubble scatterplot of the GO (molecular function, biological function, and subcellular localization) and KEGG biological pathway. The target count and significance value (P < 0.05) were mapped to size and color, respectively.

Molecular docking

The possible interaction between compounds AA, AB, and AC and the targets, including STAT3 (PDB ID: 6NJS), VEGFA (PDB ID: 4KZN), HSP90AA1 (PDB ID: 4BQG), FGF2 (PDB ID: 2FGF), and IL-2 (PDB ID: 1M48), were modeled. Each protein was retrieved from RCSB Protein Data Bank (https://www.rcsb.org/pages/policies). The interaction between each compound with its target was predicted using PyRx Virtual Screening Tools (Version 0.8)54. All the interactions among the target proteins and compounds were constructed using PyMOL 2.5 in the PDBQT format file55. Lastly, the 2-dimensional visualization of the interaction between compounds and target proteins was investigated using Discovery Studio Visualizer 2021.

Molecular dynamic simulation

SwissParam web service56 was used to create topology and parameter in CHARMM36 force field (version Jul 21)57 for AA and AB ligands. The topology and parameter of STAT3 complexes were patched and generated using pfsgen (version 2.0). The complexes were solvated by TIP3 water58 in 15 Å for each direction and ionized using VMD’s solvate and autoIonize plugins (version 1.9.2)59. The simulations were performed using NAMD (version 2.14)60. The simulation system was minimized for 1000 steps the simulated for 50,000 steps, 2 fs for each step using particle mesh Ewald, Langevin dynamic simulation under 300 K and 1 atm. The rmsd of molecular dynamic simulation were calculated using VMD’s NAMD Energy and RMSD Trajectory Tool plugins (version 1.9.2)59.