Active compounds and potential targets of GHQWG
In total, 135 chemical compounds from ten herbs in GHQWG were obtained from TCMSP databases and the literature, including 6, 11, 23, 9, 17, 18, 8, 23, 5 and 15 compounds from GMG, GHX, JYH, CZ, CH, FF, JJS, LQ, NH and QH, respectively. After removing redundancy, 90 chemical compounds were identified (Supplementary File Table 1). 312 probable targets of GHQWG active chemicals and their related gene symbols were identified by scanning the TCMSP target module (Supplementary File Table 2). The active compounds and common targets were loaded into Cytoscape 3.7.1 to create a herb-compound-target network with 411 nodes and 2134 edges (Fig. 2). The network sheds light on the multiple complicated impacts of GHQWG active chemicals on flu. We determine the key components based on the degree ranking of the active ingredients, and identify the top eight active ingredients with degree ranking as key components. A network analysis showed that quercetin (degree: 585), Kaempferol (degree: 175), Luteolin (degree: 163), β-Sitosterol (degree: 137), wogonin (degree: 127), Stigmasterol (degree: 91), Caffeic Acid (degree: 82), and Isorhamnetin (degree: 33) were key nodes.Figure 2Herb–compound–target network. Blue arrows represent the herbs in GHQWG, yellow circles node are QH compounds, green circles node are CH compounds, orange circles node are JYH compounds, deep purple circles node are GHX compounds, light purple circles node are intersection components, pinkish purple circles node are JJS compounds, dark green circles node are GMG compounds, navy blue circles node are NH compounds, bright yellow circles node are FF compounds, red–orange circles node are CZ compounds, powdered circles node are LQ compounds. The edges represent the interaction between compounds and targets, and the node size is proportional to the degree of interaction.Prediction of drug and disease related targetsWe collect target genes related to flu from five open-source databases, namely, GeneCards (1656), OMIM (5), DisGeNET (314), TTD (17), and DrugBank (232). After checking duplications, 1996 target genes were obtained (Supplementary File Table 3). Using the Venny online mapping tool, 1996 influenza-related targets and 331 GHQWG projected targets were loaded. Following mapping, 134 intersection targets of GHQWG and influenza were obtained (Fig. 3).Figure 3Disease–drug target Venn diagram.Construction and analysis of GHQWG-flu-related PPI network134 intersection targets were then imported into the STRING platform to build a PPI network, yielding 109 nodes and 397 edges. The size and color of the node denoted the size of the value. The intersection targets were screened using the double median of “Degree,” that is, “Degree 24.” As a result, the PPI network data from the STRING11.5 database was imported into the Cytoscape 3.9.1 software for visualization (Fig. 4). The PPI analysis revealed that the therapeutic targets of GHQWG had many networks and synergistic interactions. Topological analysis showed RELA proto-oncogene, NF-kBr (RELA), tumor protein p53 (TP53), mitogen-activated protein kinase 3 (MAPK3), tumor necrosis factor (TNF), AKT serine/threonine kinase 1 (AKT1), and mitogen-activated protein kinase 1 (MAPK1) occupy the core position in the PPI network and are considered hub genes of GHQWG against flu (Supplementary File Table 5).Figure 4Core target PPI network. As shown in the figure, the darker color of the circle is proportional to its importance in this network.Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway analysisTo study the biological activities and pathways of GHQWG against flu, GO and KEGG enrichment analysis on common targets were done. On the Metascape platform, the GO function enrichment study of the 134 core targets yielded 888 GO items, comprising 674 “Biological Processes (BP)”, 80 “Cellular Components (CC)”, and 134 “Molecular Functions (MF)”(Supplementary File Table 6). Based on the P value for visual analysis, the first 15 “Biological Processes” items, 10 “Cellular Components” items, and 10 “Molecular Functions” items were chosen (Fig. 5). Representative BP terms included the response to xenobiotic stimulus, positive regulation of gene expression, cellular response to cadmium ion, inflammatory response, positive regulation of the apoptotic process, positive regulation of transcription from RNA polymerase II promoter, lipopolysaccharide-mediated signaling pathway, positive regulation of cell proliferation, etc. Representative CC terms included extracellular space, macromolecular complex, membrane raft, caveola, an integral component of the plasma membrane, extracellular region, postsynaptic membrane, plasma membrane, cytosol, presynaptic membrane, etc. Representative MF terms included enzyme binding, identical protein binding, protein binding, protein homodimerization activity, heme binding, cytokine activity, G-protein coupled acetylcholine receptor activity, protein kinase binding, transcription cofactor binding, etc. In addition, 174 signal pathways were enriched (Supplementary File Table 7) by KEGG pathway analysis of the core targets using the DAVID platform (Fig. 6) representative pathways included the AGE-RAGE signaling pathway in diabetic complications, Lipid and atherosclerosis, Fluid shear stress and atherosclerosis, TNF signaling pathway, Pathways in cancer, IL- 17 signaling pathway, Chagas disease, C- type lectin receptor signaling pathway, Hepatitis B, Kaposi sarcoma-associated herpesvirus infection, Prostate cancer, Chemical carcinogenesis—reactive oxygen species, Human cytomegalovirus infection, etc. In conclusion, Go and KEGG analysis highlighted that GHQWG’s anti-inflammatory, antiviral, and anticancer properties are important targets/pathways in flu treatment.Figure 5The GO function analyzes the histogram. BP is marked in teal, CC in sienna, and MF in steel blue. The bar graph is obtained by the Bioinformatics Platform.Figure 6KEGG enrichment bubble diagram of the treatment of flu by GHQWG.Molecular dockingIn general, it is thought that a docking score of 0 kcal/mol or less suggests that the component and the target can interact spontaneously, a value of − 4.0 kcal/mol or less shows good docking affinity, a value of − 7.0 kcal/mol or less shows strong docking affinity, and the smaller the value, the higher the binding activity49,50,51. The docking targets with the 8 active components with the highest degree of Quercetin, Kaempferol, Luteolin, β-Sitosterol, Wogonin, Stigmasterol, Caffeic Acid, Isorhamnetin, and the top 6 targets of degree are RELA, TP53, MAPK3, TNF, AKT1, and MAPK1 were used for molecular docking analysis. The docking results have shown the main components and hub genes have good binding activity (affinity < − 6.00 kcal/mol) (Fig. 7) and the binding energy of quercetin and hub genes was the lowest, indicating that quercetin had the strongest binding ability. Detailed information on molecular docking is shown in Table 1.Figure 7Heat map of binding energies.Table 1 The molecular docking results of the hub genes with the main components of GHQWG.PyMoL-1.7.2.1 and Discovery Studio 2020 were used to depict the compound-target interactions with the greatest free binding energy scores, as well as their mechanisms of binding (Fig. 8A–L). Docking results had free binding energies ranging from − 3.505 to − 9.853 kcal/mol, indicating stable binding. The free binding energy of TNF with quercetin was − 9.853 kcal/mol. Binding affinities were attributed to hydrogen bonding with TYR-199 residues van der Waals forces with LEU-94, LEU-120, TYR-59, LEU-57, GLN-61, GLY-121, LEU-120, SER-60, TYR-119 and GLR-61 residues, as well as hydrophobic interactions with LEU-57 residues of TNF (Fig. 8I,J).Figure 8Molecular docking diagram of chemical composition to target: TNF to Quercetin (2D and 3D) (A,B); AKT1 to Kaempferol (2D and 3D) (C,D); MAPK1 to Kaempferol (2D and 3D) (E,F); MAPK3 to Kaempferol (2D and 3D) (G,H); RELA to Luteolin (2D and 3D) (I,J); TP53 to Caffeic-Acid (2D and 3D) (K,L).Molecular dynamics simulationMolecular dynamics simulations can be used to understand the stability of protein–ligand complexes. Based on the results of molecular docking, we selected the six compounds with the highest free binding energy scores as the targets of molecular dynamics simulations, namely AKT1-Kaempferol, MAPK1-Kaempferol, MAPK3-Kaempferol, RELA-Luteolin, TNF-Quercetin, and TP53-Caffeic Acid were subjected to simulation analyses of molecular dynamics for 100 ns to assess their motion, trajectory, structural features, binding potential, and conformational changes.The root mean square deviation (RMSD), which measures the degree of atom position departure from the initial position, is a useful indicator of the conformational stability of proteins and ligands. An improved conformational stability is shown by a decreased deviation. The complexes’ variations in RMSD values were examined. The RMSD of the AKT1–Kaempferol complex varied early on and stabilized after 0.4 nm, as seen in Fig. 9A. Although the MAPK1–Kaempferol complex’s RMSD trajectory varied between 86,000 and 90,000 ps and was largely smooth the remainder of the period (Fig. 9B), it did so. The trajectory became steady for the MAPK3–Kaempferol complex at 1 nm (Fig. 9C). The RELA–Luteolin complex’s RMSD fluctuated early on before stabilizing at 48,000 ps (Fig. 9D). RMSD of the TNF-Quercetin combination has been stabilized overall (Fig. 9E). In the beginning, the RMSD of the TP53–Caffeic Acid complex varied until stabilizing at 0.4 nm (Fig. 9F).Figure 9Root-mean-square deviation (RMSD) plots during of molecular dynamics simulation. (A) The RMSD of AKT1-Kaempferol. (B) The RMSD of MAPK1-Kaempferol. (C) The RMSD of MAPK3-Kaempferol. (D) The RMSD of RELA-Luteolin. (E) The RMSD of TNF-Quercetin. (F) The RMSD of TP53-Caffeic Acid.The active pockets of small molecules and proteins were found to be in a stable condition based on the RMSD values for the ligand and pocket. This shows that the protein’s structure does not change considerably following the interaction with the small molecule ligand, and the combination is reasonably stable.Together, these findings resoundingly validate the docking results. Temperature and pressure have no discernible impact on the structural conformation.The stability of the target protein at the residue levelThe vibrations of each residue following chemical binding were examined as root mean square fluctuations (RMSF) in order to investigate the local fluctuations of macromolecular proteins at the residue level. RMSF can be used in molecular dynamics simulations to depict the flexibility of proteins. The medication typically has the function of stabilizing the protein and exerting the effect of enzymatic activity by binding to the protein and causing a reduction in its flexibility. Protein amino acid flexibility and motion intensity are described by RMSF throughout the simulation. The medication attaches to the protein, stabilizing it and activating enzymes, as seen in the picture, which generally makes the protein less flexible. All compounds show the same tendency when compared to the respective reference ligand, i.e., they cause some fluctuation in the same protein region (Fig. 10A–F). However, all compounds generally have high rigidity and structural stability.Figure 10MD simulates RMSF trace values of protein–ligand complexes. (A) The RMSF of AKT1-Kaempferol. (B) The RMSF of MAPK1-Kaempferol. (C) The RMSF of MAPK3-Kaempferol. (D) The RMSF of RELA-Luteolin. (E) The RMSF of TNF-Quercetin. (F) The RMSF of TP53-Caffeic Acid.Analysis of the radius of gyrationThe radius of gyration (Rg) represents the compactness of the embodiment and can indicate the degree of system binding. Rg, as indicated in the figure, can be used to reflect the protein structure’s compactness. Throughout the simulation, the structural dynamics of AKT1–Kaempferol, MAPK1–Kaempferol, and TNF–Quercetin complexes remained relatively stable (Fig. 11A,B,E). MAPK3-Kaempferol, RELA-Luteolin, and TP53-Caffeic Acid fluctuated, and the Rg value was biased (Fig. 11C,D,F), indicating that kaempferol was strongly bound to most proteins. Meanwhile, the Rg of AKT1–Kaempferol complex, MAPK1–Kaempferol complex, MAPK3–Kaempferol complex and TNF–Quercetin complex were stabilised at 2.16–2.25 nm; that of RELA–Luteolin complex was stabilised at 2.80 nm and that of TP53-Caffeic Acid was stabilised at 1.67 nm. The Rg of the RELA–Luteolin complex eventually stabilised at 2.80 nm, and that of the TP53–Caffeic Acid complex eventually stabilised at 1.67 nm.Figure 11MD simulates Rg trace values of protein–ligand complexes. (A) The Rg of AKT1-Kaempferol. (B) The Rg of MAPK1-Kaempferol. (C) The Rg of MAPK3-Kaempferol. (D) The Rg of RELA-Luteolin. (E) The Rg of TNF-Quercetin. (F) The Rg of TP53-Caffeic Acid.Hydrogen bond analysisBecause hydrogen bonding is one of the most powerful non-covalent binding interactions, understanding the binding affinity between ligands and proteins is critical. The results indicated that the hydrogen bond numbers for the AKT1–Kaempferol, MAPK1–Kaempferol, MAPK3–Kaempferol, RELA–Luteolin, TNF–Quercetin, and TP53–Caffeic Acid complexes were 0–6, 0–10, 0–8, 0–8, 0–8, and 0–6, respectively (Fig. 12A–F). The number of hydrogen bonds produced by all protein–ligand complexes remains constant throughout the simulation. Furthermore, the presence of persistent amino acid residues at the active site contributes to the complex’s overall structural stability.Figure 12MD simulates the number of hydrogen bonds (HBond) of protein–ligand complexes. (A) The HBond of AKT1-Kaempferol. (B) The HBond of MAPK1-Kaempferol. (C) The HBond of MAPK3-Kaempferol.(D)The HBond of RELA-Luteolin.(E)The HBond of TNF-Quercetin. (F) The HBond of TP53-Caffeic Acid.Analysis of solvent accessible surface areaThe interface enclosed by solvent is used to compute the solvent accessible surface area (SASA)52. Because this solvent behaves differently under different conditions, it can be used to examine protein conformational dynamics in a solvent context. The contact area between the six complexes and water is comparable, and the small molecule has little effect on the protein and water effect, the results are shown in (Fig. 13A–F). The SASA values of AKT1–Kaempferol complex, TP53–Caffeic Acid complex, MAPK3–Kaempferol complex and AKT1–Kaempferol complex remained generally unchanged during the simulation, RELA–Luteolin complex decreased from the initial 230 nm2 to 195 nm2 and TNF–Quercetin complex decreased from the initial 215 nm2 to 185 nm2 during the simulation, indicating that the protein–protein interactions have little effect on the characterisation and stability of the surface of protein molecules.Figure 13MD simulates SASA of protein–ligand complexes. (A) The SASA of AKT1-Kaempferol. (B) The SASA of MAPK1-Kaempferol. (C) The SASA of MAPK3-Kaempferol. (D) The SASA of RELA-Luteolin. (E) The SASA of TNF-Quercetin. (F) The SASA of TP53-Caffeic Acid.
Mechanisms underlying the therapeutic effects of Gang Huo Qing wen granules in the treatment of influenza based on network pharmacology, molecular docking and molecular dynamics
