Network Pharmacology and Molecular Docking-Based Prediction of the Mechanism of Qianghuo Shengshi Decoction against Rheumatoid Arthritis

Qianghuo Shengshi decoction (QHSSD) is a classical Chinese medicine formula, which is used in clinical practice for the treatment of rheumatoid arthritis (RA) in China. However, the pharmacological mechanism of QHSSD on RA has remained unclear by now. We collected and screened active compounds and its potential targets by the pharmacology platform of Chinese herbal medicines. In addition, the therapeutic targets of RA were obtained and selected from databases. Network construction analyzed that 128 active compounds may act on 87 candidate targets and identified a total of 18 hub targets. GO annotation and KEGG enrichment investigated that the action mechanism underlying the treatment of RA by QHSSD might be involved in cell proliferation, angiogenesis, anti-inflammation, and antioxidation. Finally, molecular docking verification showed that TP53, VEGFA, TNF, EGFR, and NOS3 may be related to the RA treatment and molecular dynamics simulation showed the stability of protein-ligand interactions. In this work, QHSSD might exert therapeutic effect through a multicomponent, multitarget, and multipathway in RA from a holistic aspect, which provides basis for its mechanism of action and subsequent experiments.


Introduction
Rheumatoid arthritis (RA) is a chronic autoimmune disease with the following symptoms: tender, warm, swollen, and stiff joints in the morning and as RA progresses, resulting in bone erosion, joint deformity, and physical disabilities [1][2][3][4]. RA affects approximately 1% of the worldwide population and mostly occurs in middle-aged women [5]. The drug treatment of RA mainly focuses on nonsteroidal antiinflammatory drugs (NSAIDs), steroids, disease-modifying antirheumatic drugs (DMARDs), and glucocorticoids. However, these drugs can only delay the progression of RA, which may also increase the risk of infection. Long-term use of these drugs will cause drug dependence and a cascade of side effects, such as digestive tract irritation reaction, rashes, liver toxicity, and hair loss [2,6,7].
There is a rising trend of the application of traditional Chinese medicine in the treatment of RA, and the RA patients are more inclined to choose complementary medicine as a medical treatment [8][9][10]. Traditional Chinese medicine (TCM) has a long history of clinical use in China, used as complementary and alternative medicines in treatment [11]. Qianghuo Shengshi decoction (QHSSD) is a classic clinical prescription from internal and external causes (Nei Wai Shang Bian Huo Lun). Accordingly, QHSSD can dispel wind and dampness and sweat and alleviate pain, consisting of seven herbs: Notopterygii rhizoma et radix (Qianghuo, QH), Angelicae pubescentis radix (Duhuo, DH), Saposhnikoviae radix (Fangfeng FF), Ligustici rhizoma et radix (Gaoben, GB), Chuanxiong rhizoma (Chuangxiong, CX), Viticis fructus (Manjingzi, MJZ), and Glycyrrhizae radix et rhizoma (Gancao, GC). QHSSD has been used to treat influenza, allergic purpura, and inflammatory diseases, such as RA [12][13][14][15]. Unlike the single-target agent with a clear mechanism, however, the mechanism of QHSSD underlying the treatment of RA still remains a blur.
Network pharmacology is a novel approach to understand drug's pharmacological mechanism in the network perspective and to find out the interaction between the drug and the body on the basis of biological networks [16]. There is a growing number of researches that applied the network pharmacology to elucidate the mechanism between TCM/herbs and disease treatment, involving multicompound and multitarget ones [17][18][19][20][21][22][23][24]. We therefore tried to explore the mechanism of action of QHSSD on RA by network pharmacology and to further predict the interaction modes between QHSSD and the predicted targets by molecular docking and molecular dynamics. The flowchart of our study is presented in Figure 1.

Materials and Methods
2.1. Collecting and Screening of Active Components of QHSSD. We collected the chemical compounds of QHSSD from the traditional Chinese medicine systems pharmacology database and analysis platform (TCMSP) (https:// tcmspw.com/tcmsp.php), a common tool for network pharmacology research of TCM [25]. Oral bioavailability (OB) and druglikeness (DL) were set as the main parameters to screen the active components according to the absorption, distribution, metabolism, and excretion (ADME) process. We set OB ≥ 30% and DL ≥ 0:18 as the screening threshold values to select the candidate components through literature search [23,26,27].

Potential Targets of QHSSD.
Considering that the target screening criteria differ between different databases (including TCMSP, PubChem, and Swiss Target Prediction databases), therefore, only the TCMSP database was applied in this work. We predicted relevant targets (DrugBank, https://go.drugbank.com/) of active components in QHSSD on the TCMSP and then transformed the target name to a standard gene name on Uniport (https://www.uniprot.org/) database and removed the duplications.

Compound-Target Network Construction and Analysis.
Venny 2.1.0 (https://bioinfogp.cnb.csic.es/tools/venny/) [30] was used to find out the overlapping targets between compound targets and disease targets. To explore the relationship between compounds and disease more reasonably, Cytoscape 3.7.2 (https://cytoscape.org/) [31] with a visualized tool and all node degrees of networks calculated was used to construct a compound-potential target network.
2.5. Protein-Protein Interaction (PPI) Network Construction and Analysis. The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (https://string-db.org/) was used to predict the interaction between proteins [32]. According to the previous steps, drug-associated genes and disease targets were analyzed on STRING to create a protein-protein interaction network (PPI network) by    [36]. In this work, the GROMACS Program package [37] and CHARMM36 force field [38] were used to simulate the molecular dynamics of the complexes in the TIP3P water model. Ions (Na + or Cl − ) were added to neutralize charges. The energy of systems was minimized through a 4000-step steepest descent method and a 2000-step conjugate gradient method to eliminate high-energy collisions between molecules. Then, the systems were heated from 0 K to 300 K within 100 ps and simulated as an NPT ensemble with normal temperature (300 K) and normal pressure (101 kPa) for 100 ps; unconstrained MD simulation was carried out for 10 ns after the system reached dynamic equilibrium. Periodic boundary conditions and particle mesh Ewald (PME) methods were used to deal with long-range electrostatic interactions. The short-range Coulomb interaction was 0.9 nm, and the SHAKE algorithm was used to restrict the bond to the hydrogen atom. The temperature and pressure of the system were controlled by V-rescale and Parrinello-Rahman. The stability of systems was evaluated using the root-mean-square deviation (RMSD) of the aligned protein-ligand coordinate set calculated against the initial frame.

Results
3.1. Compound-Target Network and Analysis. We collected 153 compounds and 152 putative targets of QHSSD (Table S1 and Figure S1) on the screening criteria from the TCMSP. Among these compounds, we excluded 25 compounds without RA-related targets and 128 active compounds were obtained in this work. 1846 targets were selected from 29368 genes related to RA by the screening threshold of inference score ≥ 24:46 on the CTD database. Meanwhile, 306 genes were obtained from 4465 targets with a relevance score ≥ 19:94 on GeneCards. Targets from the above database merged and removed the overlapping parts, and at last, we received a total of 1992 RA-related targets. Based on the Venn diagram, we obtained 87 overlapping targets from 152 putative targets of QHSSD and 1992 RA-related targets. Consequently, 87 targets ( Figure 2 and Table 1) were considered as potential targets for the treatment. To further investigate the interaction between the compounds and targets, we constructed a compound-target network ( Figure 3) with 215 nodes and 1239 edges. The network results manifested that the average degree value of the compounds was 9.68; we received 66 active compounds. Five compounds were considered to have high connection with putative targets of RA: quercetin (degree = 57), kaempferol (degree = 31), luteolin (degree = 27), wogonin (degree = 21), and isorhamnetin (degree = 20) (Table 2).

PPI Network and Analysis.
Based on the previous results, 87 potential targets were imported in STRING with the organism set as "homo sapiens." Those data with a confidence level higher than 0.4 constructed a complete PPI network including 87 nodes and 1062 edges. We chose three important parameters, including "degree (DC)," "betweenness centrality (BC)," and "closeness centrality (CC)," as the thresholds to screen the key genes to establish major hub nodes by Cytoscape. First, we set the medium of DC, BC, and CC as thresholds, DC ≥ 21, BC ≥ 0:003, and CC ≥ 0:555, respectively, and obtained 37 nodes and 472 edges. Then, the 37 key nodes were further filtered by the second threshold of DC ≥ 27, BC ≥ 0:006, and CC ≥ 0:800 and we got 18 nodes and 153 edges at last ( Figure 4).     and atherosclerosis (hsa05418) and relaxin signaling pathway (hsa04926). The compounds of QHSSD may act on these signaling pathways and yield a therapeutic effect.

Molecular Docking Analysis.
To further elucidate the mechanism of action, 18 hub targets with their related compounds were docked and screened by the binding affinity on AutoDock Vina based on the compound-target network. We investigated a total of 39 pairs of ligand-protein interaction (Table 3), with binding affinity < −5 kcal/mol as the screening threshold. The lower the value of the docking affinity, the stronger the binding ability between the compound and the active site of the target. Accordingly, JUN and TNF get a low binding ability with their related compounds, as well as JUN-quercetin, JUN-kaempferol, and TNFquercetin. Therefore, 3 pairs of combinations were excluded, and at last, we retrieved 47 compound-target pairs. In addition, cognate docking was used to verify the rationality of molecular docking. In order to verify the rationality of the docking for the TP53 and quercetin system, we investigated the docking of IRAK4 (PDB ID: 2a9i) and quercetin; the binding affinity of IRAK4-quercetin was −6.6 kcal/mol. Small-molecule ligand potentially could fit into the interface port formed by the interaction of amino acid residues in protein. Overall, these complexes provide further evidence that these proteins could act as therapeutic target for the compound of QHSSD.

Molecular Dynamics Simulations.
To further research the stability of the molecular dynamics trajectory, we took a complex (TP53-quercetin) as an example to analyze the root mean square deviation (RMSD) of the TP53 backbone atoms in this complex, as shown in Figure 6. The RMSD of the TP53 backbone atom began to show an upward trend, whose trajectory was stable and reached an equilibrium state. The RMSD remained stable at about 0.2 nm after 4 ns, but in 7 ns, the trajectory fluctuated slightly. In addition, after 2 ns of the molecular dynamics simulation, the potential energy of the system also tends to be the minimum and stabilizes. The total interaction energy value was − 136:2 ± 19:0 kJ/mol. The 4 ns simulated trajectory was extracted by GROMACS as the representative frame of the interaction between quercetin and TP53 ( Figure 7); we could find that the amino acid ARG15 of TP53 and quercetin have hydrophobic interactions, resulting in quercetin and TP53 more closely. Moreover, quercetin and TP53 could form three hydrogen bonds (the hydrogen bond between the oxygen atom of quercetin and the hydrogen atom of amino acid residue GLY204 and the hydrogen bond between the hydrogen atom of quercetin and the nitrogen atom of amino acid residues GLU209 and GLN28) and caused the complex to bind more tightly.

Discussion
There has been increased researches on how traditional Chinese herbal medicine or formulae work on various diseases based on network pharmacology [21,22,24]. In our study, we first collected the putative targets of compounds in QHSSD from the database. And then, the compoundtarget network was constructed, which revealed that one herb could correspond with multiple components and multiple targets. The traditional Chinese medicine integrative database (TCMID) and a Bioinformatics Analysis Tool for Molecular Mechanism of Traditional Chinese Medicine (BATMAN-TCM) database were used to predict the targets of Chinese herbal medicine [27,28], whereas we found that the results of these databases were different because of the screening criteria. To ensure the consistency of data, we only used the TCMSP database to collect active compounds and related targets. In compound-target network analysis, the top five compounds (quercetin, kaempferol, luteolin, wogonin, and isorhamnetin) were considered to be crucial compounds in QHSSD. Surprisingly, these compounds are both flavonoids which have antioxidant and antiinflammatory effects [39][40][41][42]. For instance, researches indicated that quercetin mediate the anti-inflammatory mechanism of reducing acute laryngitis by restoring the Th17/Treg balance and activating heme oxygenase 1 and isorhamnetin increased the levels of TNF-α, IL-1β, IL-6,   IL-17A, and IL-17F as well as decreased levels of IL-35 and IL-10 in the CIA mice, resulting in an anti-inflammatory effect [43,44]. Luteolin prevented cartilage destruction in OA rats in vivo through inhibiting IL-1β-induced phosphorylation of NF-κB [45]. In conclusion, these active compounds are crucial material basis of the potential therapeutic effect of QHSSD on RA. Further, we selected 18 hub targets from the PPT network, as follows: IL-6, TP53, VEGFA, JUN, TNF, MAPK1, MAPK8, EGF, PTGS2, MAPK3, EGFR, ESR1, IL-1B, CAT, CCL2, MAPK14, MMP2, and NOS3. Combined with GO annotation and pathway enrichment analysis, we preliminarily considered that the treatment of RA by QHSSD may participate in biological processes of proliferation, inflammatory response, and angiogenesis. Previous studies indicated that IL-6 and TNF antagonists decrease the inflammatory response of RA patients and IL-1B stimulates inflammation and degradation of the bone and cartilage [46][47][48][49]. EGF is an important growth factor, regulating the proliferation of cells. MAPK enables the growth and proliferation of cells by binding with other protein kinases or factors [50]. In addition, the formation of new blood vessels supplied nutrients and oxygen to the augmented inflammatory cell mass, contributing to perpetuation of joint disease [51]. VEGFA, a proangiogenic molecule, promotes the angiogenic phenotype of RA and is upregulated in RA [52,53]. These literatures are consistent with our prediction. Therefore, these hub targets are considered as crucial proteins in RA therapy.
To further verify the hypothesis of action mechanism, we performed molecular docking on the hub target binding with crucial compounds. Surprisingly, we found that quercetin, luteolin, and wogonin could both fitly bind with TP53, a transcription factor, regulating cell cycle. Moreover, quercetin and luteolin also act on VEGFA. Previous studies revealed that TP53 mutations are linked to the VEGF   Figure 6: RMSD change of TP53 backbone atoms in MD simulation. 8 BioMed Research International pathway and showed that the transcriptional target of TP53 inhibited proliferation and angiogenesis via VEGFA [54][55][56]. We therefore hypothesized that quercetin and luteolin positively worked on TP53 to regulate cell cycle, and meanwhile act on VEGFA to suppress vascular formation. Additionally, luteolin could simultaneously act on EGFR and TNF, decreasing cell proliferation and the level of proinflammatory cytokines. This is in line with the researches before; the result of which manifested that luteolin inhibited IL-1β-induced expression of NO, PGE2, TNFα, and MMP-2 [45,50,57,58]. On the one hand, the synergistic action of quercetin and luteolin produces the effect of inhibiting cell proliferation and angiogenesis by multiple proteins and also has an anti-inflammatory effect. On the other hand, kaempferol acting on NOS3 may increase NO levels and prevent oxidative damage [59]. From a holistic aspect, QHSSD could inhibit the cell proliferation and angiogenesis and yield anti-inflammatory and antioxidant effects through regulating TP53, VEGFA, EGFR, TNF, and NOS3. MD is crucial to the discovery of new drug [60] and is a powerful method for studying protein-ligand interactions [61]. In contrast to experimental methods alone, MD simulations can address diseases caused by protein misfolding and virtual screening; it also identifies the stability of protein-ligand complexes and ligand binding kinetics [62]. Previous studies found that quercetin can play an important role in inflammation by acting on the TP53 target [54][55][56]. In our study, MD simulation was used to investigate the interactions of TP53-quercetin. This complex was simulated within 10 ns and reached an equilibrium state. We found that ligand has three hydrogen bonds with protein; a study [63] had shown that hydrogen bonding has an advantage over shapeless forces if it can form hydrogen bonds. In this study, MD simulation revealed that the structure of TP53quercetin was very stable when the simulation temperature was 300 K. It is important to explore the application of new drug by molecular docking and MD simulations [61]; however, it is still necessary to perform in vitro experiment to verify the prediction.

Conclusion
To conclude, we could learn that multicompounds act on multitargets in TCM therapy. We speculated that the therapeutic effect on the treatment of RA by QHSSD may be focused on antiangiogenesis, meanwhile, anti-inflammation, especially TP53 and VEGFA targets, plays a vital role during the treatment. This study was based on the combination of network pharmacology and molecular docking to predict and identify the mechanism of QHSSD in the treatment of RA, which provides the foundation for subsequent experimental investigation and offers ideas for the multidimensional and multilevel research of TCM formulae.

Data Availability
The data used to support the findings of this study are included within the supplementary information files.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this article.

Authors' Contributions
Xiaoli Bi designed this study and provided technical support toward the study. Zhihao Zeng, Guanlin Xiao, Ruipei Yang, Huajing Huang, and Huixian Zhong collected the data. Zhihao Zeng and Jiaoting Hu analyzed the data and wrote the manuscript. Jieyi Jiang, Sumei Li, and Yangxue Li revised the manuscript.

Acknowledgments
Thanks are due for Warren L. DeLano, the original PyMOL author of the open-source software. Our research was supported by the "Significant News Drugs Creation" Science and Technology Major Project of China (grant number 2018ZX09721005).