Deciphering the Molecular Targets and Mechanisms of HGWD in the Treatment of Rheumatoid Arthritis via Network Pharmacology and Molecular Docking

Background Huangqi Guizhi Wuwu Decoction (HGWD) has been applied in the treatment of joint pain for more than 1000 years in China. Currently, most physicians use HGWD to treat rheumatoid arthritis (RA), and it has proved to have high efficacy. Therefore, it is necessary to explore the potential mechanism of action of HGWD in RA treatment based on network pharmacology and molecular docking methods. Methods The active compounds of HGWD were collected, and their targets were identified from the Traditional Chinese Medicine Systems Pharmacology Database (TCMSP) and DrugBank database, respectively. The RA-related targets were retrieved by analyzing the differentially expressed genes between RA patients and healthy individuals. Subsequently, the compound-target network of HGWD was constructed and visualized through Cytoscape 3.8.0 software. Protein-protein interaction (PPI) network was constructed to explore the potential mechanisms of HGWD on RA using the plugin BisoGenet of Cytoscape 3.8.0 software. Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) were performed in R software (Bioconductor, clusterProfiler). Afterward, molecular docking was used to analyze the binding force of the top 10 active compounds with target proteins of VCAM1, CTNNB1, and JUN. Results Cumulatively, 790 active compounds and 1006 targets of HGWD were identified. A total of 4570 differentially expressed genes of RA with a p value <0.05 and |log 2(fold change)| > 0.5 were collected. Moreover, 739 GO entries of HGWD on RA were identified, and 79 pathways were screened based on GO and KEGG analysis. The core target gene of HGWD in RA treatment was JUN. Other key target genes included FOS, CCND1, IL6, E2F2, and ICAM1. It was confirmed that the TNF signaling pathway and IL-17 signaling pathway are important pathways of HGWD in the treatment of RA. The molecular docking results revealed that the top 10 active compounds of HGWD had a strong binding to the target proteins of VCAM1, CTNNB1, and JUN. Conclusion HGWD has important active compounds such as quercetin, kaempferol, and beta-sitosterol, which exert its therapeutic effect on multiple targets and multiple pathways.


Introduction
Rheumatoid arthritis (RA), a common chronic systemic autoimmune disease, is characterized by synovial hyperproliferation and inflammatory/immune cell infiltration [1,2]. Typical clinical symptoms of RA include tender and swelling of the joints, accompanied by morning stiffness of the affected joints. In severe cases, large joints can be injured leading to joint deformity and loss of function. e current global incidence of RA is 0.24%, and it is expected to increase [3][4][5][6]. Besides, RA is ranked 74 th as a social burden and 42 nd as a disability [3][4][5][6]. Current international guidelines for the management of RA recommend the use of disease-modifying antirheumatic drugs (DMARDs), with methotrexate (MTX) as the first-line drug [7,8]. However, MTX is not an ideal therapeutic agent. It has damaging side effects on the neuronal, gastrointestinal, and immune systems [9]. erefore, it is imperative to explore safe and effective clinical treatment for RA. e traditional Chinese medicine (TCM) is based on the theory of syndrome differentiation and has long been established as an effective treatment of RA. Huangqi Guizhi Wuwu Decoction (HGWD), a classical prescription described in Jingui Yaolue, has been used over 1000 years in China to treat joint pain. HGWD was formulated by Zhang Zhongjing during the Han dynasty. It is a mixture of five Chinese medicines including Astragalus membranaceus (Huangqi), Cassia Twig (Guizhi), Radix Paeoniae Alba (Baishao), ginger (Shengjiang), and Chinese date (Dazao). HGWD is considered to have the ability to reinforce Qi and nourish blood, warm and smoothen meridian, and dredge arthritic pain traditionally [10]. In a study of rats with RA, HGWD administration reduced joint inflammation, synovial hyperplasia, and cartilage damage [11]. Previous studies have revealed that HGWD improves the clinical symptoms and signs, as well as laboratory indices, in RA patients [12]. Moreover, a systematic review and metaanalysis showed that when HGWD was combined with Western medicine therapy to treat RA, better efficacy, improved morning stiffness, reduced C-reactive protein, and rheumatoid factor content were achieved compared to Western medicine only [13]. However, the mechanisms of HGWD in RA treatment are not clear, a major limiting factor for its extensive application.
Network pharmacology is partly bioinformatics and was first proposed by Hopkins [14]. It has been successfully used to study complex TCM formulations. is is because it not only combines system network analysis and pharmacology but is also based on the connotation of holistic theory, multicomponents, multitargets, and multipathways of Chinese medicine [15,16]. Network pharmacology can elucidate the mechanisms of HGWD in the RA treatment at the molecular level via compoundcompound network, compound-target network, and target-disease network.
In the present study, network pharmacology was used to establish a compound-target-disease network for exploring the potential HGWD mechanism of action in RA treatment.
is study provides a reference for future pharmacological studies and clinical applications. e flow diagram of the network is shown in Figure 1.

Collection of Active HGWD Compound Information.
Information on the HGWD compounds was retrieved from the Traditional Chinese Medicine Systems Pharmacology Database (TCMSP, http://tcmspw.com/tcmsp.php) [17]. TCMSP is a unique system pharmacology database of Chinese herbal medicines with data on absorption, distribution, metabolism, and excretion (ADME)-related parameters of herbal ingredients as well as the relationships between diseases, targets, ingredients, and drugs [18]. e active compounds of HGWD were primarily screened based on oral bioavailability (OB) and drug-like (DL) properties, the two important indicators for bioinformatics evaluation of ADME characteristics [19]. e OB is a major pharmacokinetic parameter of oral drugs and is the proportion of oral drug dose in the systemic circulation [20]. DL properties are the physical and chemical properties that qualitatively evaluate whether a compound is similar to existing drugs [21]. Subsequently, the compounds with OB ≥ 30% and DL ≥ 0.18 were chosen as the candidate compounds for further analysis.

Identification of HGWD Potential
Targets. Target identification is an important aspect of drug exploration [22]. e target candidate compounds were imported into the DrugBank database (https://www.drugbank.ca/) to identify the corresponding potential targets of HGWD [23].

Known RA-Related Targets.
e differentially expressed genes in RA were retrieved from the GEO database (https:// www.ncbi.nlm.nih.gov; series: GSE21959; samples: there were 36 samples, 18 for healthy individuals and 18 for those with rheumatoid arthritis). e genes with a p value <0.05 and |log 2(fold change)| > 0.5 were regarded as differentially expressed genes and RA-related targets.

Network Construction.
e compound-target network of HGWD was constructed and visualized via Cytoscape 3.8.0 software [24]. Protein-protein interaction (PPI) networks were constructed to explore the potential mechanisms of HGWD on RA using the plugin BisoGenet of Cytoscape 3.8.0 software [25]. Afterward, the topological importance of each node by calculating degree centrality (DC), betweenness centrality (BC), closeness centrality (CC), eigenvector centrality (EC), local average connectivity-based method (LAC), and network centrality (NC) was evaluated with a Cytoscape plugin CytoNCA. eir definitions and computational formulas have been reported and are the topological importance representative of each node [26].

Gene Ontology and Pathway Enrichment Analysis.
Gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) are important methods that describe the features of candidate targets. e two were performed in R software (Bioconductor, clusterProfiler) with the standard p value cutoff of 0.05 and the q value of 0.05 [27,28]. e GO analysis was applied for target protein analysis, and the top 20 functional categories in biological process (BP), cellular component (CC), and molecular function (MF) were chosen. Based on the targets of HGWD in RA treatment, KEGG enrichment analysis was carried out. e top 20 KEGG pathways were selected for plotting a histogram. Meanwhile, the network diagram of the gene pathway was drawn. Evidence-Based Complementary and Alternative Medicine 2.6. Molecular Docking of the Main Active Constituents of HGWD and Core Proteins. e 3D protein structure of the three core proteins corresponding to the core targets, VCAM1, CTNNB1, and JUN, was downloaded from the UniProt database. Subsequently, the structure of the active ingredient in HGWD (the top 10 places in the number of targets) was downloaded from the PubChem database and saved in the PDB format. Using PyMOL software, the three proteins were virtually dehydrated and hydrogenated, the original ligands were extracted in each protein, and then were stored separately. AutoDockTools 1.5.6 was utilized to convert compounds, ligands, and proteins into the "pdbqt" format and to define if the location of each protein or its ligands was the active pocket of the protein. Finally, Vina 1.5.6 was run to assess molecular docking. At a binding energy value <0, the molecular proteins were considered to spontaneously bind and interact with each other. Accordingly, the lower the energy is, the more stable the molecular conformation is.

Target Screening of HGWD and RA.
A total of 790 compounds of the five herb medicines in HGWD were retrieved from the TCMSP database.
is included 87 compounds in Huangqi, 220 in Guizhi, 85 in Baishao, 265 in Shengjiang, and 133 in Dazao. Among them, 74 compounds passed OB ≥ 30% and DL ≥ 0.18 filtering. Specifically, the numbers of candidate compounds in Huangqi, Guizhi, Baishao, Shengjiang, and Dazao were 20, 7, 13, 5, and 29, respectively. e candidate compounds in HGWD used for further analysis are shown in Table 1. e DrugBank database retrieval predicted a total of 1006 potential targets. e potential targets linked to Huangqi, Guizhi, Baishao, Shengjiang, and Dazao were 405, 57, 104, 60, and 380, respectively. Moreover, a total of 4570 differentially expressed genes in RA were collected from the GEO database. As shown in Figure 2, a volcano plot was drawn to show the distribution of the differentially expressed genes. e genes are represented by the red and green dots in the plot. We compared the target genes regulated by the active compounds in HGWD, and different genes in RA were compared, obtaining 49 common target genes. ese 49 target genes were found to be regulated by 28 active compounds.

Compound-Target Network Analysis.
e compoundtarget network of HGWD was established with the collected compounds and their targets as shown in Figure 3. e network contains 77 nodes (28 compounds in HGWD and 49 compound targets) and 130 edges elucidating the compound-target interactions. Quercetin, kaempferol, and betasitosterol acted on 33, 14, and 8 targets, respectively. e OB of quercetin, kaempferol, and beta-sitosterol was 43.43, 41.88, and 36.91%, respectively. erefore, they are potential key active compounds due to their relative positioning in the network. Besides, it has been revealed that Guizhi, Baishao, and Dazao share a common component (+)-catechin (ID: MOL000492); Shengjiang and Dazao have the same ingredient stigmasterol (ID: MOL000449); Guizhi, Baishao, Shengjiang, and Dazao share a common ingredient beta-sitosterol (ID: MOL000358); Huangqi and Baishao share the same constituent kaempferol (ID: MOL000422); and Huangqi and Dazao share the same component quercetin (ID: MOL000098).

e Candidate Targets for HGWD against RA.
To elucidate the mechanism by which HGWD ameliorates RA, the potential target network was merged with the RA-related Evidence-Based Complementary and Alternative Medicine  (Figure 4(a)). e previous research showed that the median degree of all nodes was more than 74 degrees which were identified as significant targets [29]. Accordingly, a network of the significant targets for HGWD against RA with 426 nodes and 17,112 edges with a DC ≥ 74 was constructed ( Figure 4(b)). During the second screening, since the number of genes was limited, only BC average value was used. e candidate targets were further identified, and 111 targets with BC ≥ 348.07 (BC average value) were chosen (Figure 4(c)). Eventually, 111 HGWD target genes against RA were identified. VCAM1 (degree: 452), CTNNB1 (degree: 390), and JUN (degree: 274) might be the most important among 111 target genes for HGWD against RA.

GO and KEGG Pathway Enrichment Analysis.
To further confirm the biological responses from RA treatment with HGWD, GO analysis of the 49 RA-related potential therapeutic target genes was performed based on BP, CC, and MF. e analysis results revealed a total of 739 entries. In BP enrichment analysis, 666 entries were obtained including response to antibiotic, response to alcohol, and response to steroid hormone. In CC enrichment analysis, 34 entries involved in membrane raft, membrane microdomain, membrane region, etc. were obtained. e MF enrichment analysis revealed 39 entries, including kinase regulator activity, protein kinase regulator activity, and serine hydrolase activity. e top 20 terms are shown in Figure 5.
To further show the biological processes of these targets, the KEGG pathway analysis was performed.
e analysis results showed that these targets were enriched in 79 pathways with an adjusted p value <0.05. e top 20 KEGG pathways' enrichment analysis is shown in Figure 6. Among these pathways, the top three were fluid shear stress and atherosclerosis, TNF signaling pathway, and IL-17 signaling pathway based on the number of the pathway target genes.

e Gene-Pathway Network.
e gene-pathway network analysis was constructed according to the significantly enriched pathways and genes which regulate these pathways as shown in Figure 7. e diagram shows the network relationship between the top 20 pathways and their regulated target genes. According to the network analysis, JUN had the            Evidence-Based Complementary and Alternative Medicine 7 CTNNB1 (PDB: 1jdh), and beta-carotene and stigmasterol had the strongest binding force with JUN (PDB: 2g01). e corresponding molecular docking diagram of stepholidine, stigmasterol, and VCAM1, beta-carotene, quercetin and CTNNB1, and beta-carotene, stigmasterol, and JUN is illustrated in Figure 8.

Discussion
TCM, characterized by multicompound and multitarget medicines, cures diseases via multiple targets, multiple pathways, and multiple links. Due to the complex chemical ingredients of TCM, conventional research methods such as clinical and pharmacological research are not capable of fully elucidating the mechanism of action of TCM. Fortunately, network pharmacology provides a solution to this challenge since it is suitable for multicompound and multitarget research. In the present study, the mechanism of action of HGWD on RA was explored via network pharmacology methods. is provided a clear direction for further basic and clinical research. Modern pharmacology has shown that HGWD has a specific therapeutic effect on RA, where it elicits obvious anti-inflammatory and analgesic effects [30]. Specifically, Huangqi has an anti-inflammatory effect and improves the immunity [31]. e volatile oil of Guizhi has good antiinflammatory, immunological, and chondrocyte proliferation effects [32]. e total glycosides, the main components of Baishao, have anti-inflammatory, analgesic, and autoimmunological effects [33]. Shengjiang resists oxidation, scavenges free radicals, and relieves pain [34]. e polysaccharide of Dazao can significantly inhibit proinflammatory cytokines, such as IL-6 and TNF, with antiinflammatory activity [35]. It is suggested that this prescription has anti-inflammatory, antioxidation, analgesic, and autoimmune reaction pharmacological effects, thus providing specific pharmacological basis for its clinical function in the treatment of RA.
In this study, the compound-target network of HGWD was constructed using 28 compounds and their 49 targets. e network diagram demonstrated that most of the HGWD compounds affected multiple targets. For instance, quercetin, kaempferol, and beta-sitosterol acted on 33, 14, and 8 targets, respectively. erefore, they were the likely key active compounds for HGWD. Quercetin is a common active component of Huangqi and Dazao. Kaempferol is a common active component of Huangqi and Baishao. Betasitosterol is a common active component of Guizhi, Baishao,  Shengjiang, and Dazao. ese drugs have been mainly used together in Chinese medical history. It indicates that the compatibility between these drugs has a synergistic effect and increases their efficacy.
Quercetin is a flavonol that has unique therapeutic biological properties including anticarcinogenic, anti-inflammatory, antioxidant, antiviral, and psychostimulant activities [36,37]. Kaempferol, a representative flavonoid, has been known to exert a range of pharmacological actions including the mediation of antioxidant, antimicrobial, and anti-inflammatory effects [38,39]. Kaempferol is a potent immunosuppressant that reduces the harmful immune responses including autoimmunity and chronic inflammation [40]. Beta-sitosterol, a main constituent of plants and vegetables, has versatile activities that impact cell activities including anti-inflammatory effect [41,42]. In this study, quercetin, kaempferol, and beta-sitosterol regulated most RA-related targets and exhibited anti-inflammatory effects. Besides, they have a high oral bioavailability. erefore, they were identified as the representative compounds of HGWD. e PPI networks of HGWD targets and RA-related targets were constructed and merged to identify the candidate targets of HGWD against RA. To accurately obtain the targets, two parameters including DC and BC were set to construct a new network. rough PPI analysis, VCAM1, CTNNB1, and JUN were established as the important targets in RA treatment. VCAM-1 is a glycoprotein expressed in vascular endothelial cells, whose serum is positively correlated with RA [43]. Studies have found that [44] TNF-α upregulates the expression of VCAM-1 in the endothelial cell membrane at the small joint synovium. VCAM-1 binds to α4β1 on the surface of free immune cells in the blood vessels, inducing the migration of immune cells to inflammatory joints thus expanding the inflammatory response cascade and aggravating small joint injury. Proliferation and invasion of fibroblast synovial cells are important mechanisms that result in the thickened synovial lining, increased secretion of inflammatory cytokines, and cartilage and bone injury [45]. Studies have shown that inhibition of CTNNB1 transcription reduces the proliferation of fibroblast synovial cells and the levels of proinflammatory cytokines such as interleukin-6 (IL-6), IL-10, and TNF-α [46]. Elsewhere, it has been shown that CD44 can activate the transcription factor AP-1 (the whole is a protein) expressed by the gene JUN, thus promoting the activation of T cells and aggravating RA synovitis [47]. e targets of HGWD against RA were enriched in BP, CC, and MF through GO analysis. Results suggested that HGWD mainly regulated response to antibiotics, response to steroid hormones, and response to radiation. Furthermore, it was shown to affect some cellular components and molecular functions including membrane raft, membrane microdomain, membrane region, and kinase regulator activity. In the present study, 79 KEGG pathways including TNF signaling pathway and IL-17 signaling pathway were significantly enriched. IL-17 and TNF-α are classic proinflammatory cytokines. IL-17 is mainly secreted by helper T-cells type 17 and has strong proinflammatory effects. Analysis of synovial fluid in patients with RA found that IL-17 content was significantly increased and positively correlated with Disease Activity Score-28 (DAS28) [48]. is was a confirmation that IL-17 is involved in the occurrence of RA. Fischer et al. [49] found that, after IL-17 treatment of fibroblast synovial cells, the expression of proinflammatory cytokines such as IL-6, IL-8, and granulocyte-macrophage colony-stimulating factor (GM-CSF) increases, aggravating the inflammatory damage of fibroblast synovial cells. Besides, IL-17 can also induce joint bone and cartilage damages [50]. It has been revealed that IL-17 upregulates the levels of nuclear factor-kappa B (NF-κB) receptor activator ligand RANKL. Consequently, this increases the content of matrix metalloproteinases, possibly promoting cartilage degradation, inhibiting chondroprogenitor cell formation, and stimulating osteoclast bone damaging [51][52][53].

Stigmasterol-JUN
TNF-α plays an important role in the development of RA. On the one hand, it can bind to the TNFR1 receptor on fibroblast synovial cells, promoting the release of inflammatory cytokines such as IL-1, IL-6, and IL-8 and aggravating the damage of the articular cartilage and bone [54]. On the other hand, TNF-α can stimulate immune cells in blood to enter the joint through VCAM-1, aggravating the joint injury [44]. Besides, the other pathways in the first 20 pathways have not been reported to be related to RA, providing new ideas and clues for future research. e gene-pathway network was constructed to explore the main target genes of HGWD against RA. e results revealed that JUN, FOS, CCND1, IL6, E2F2, and ICAM1 are important target genes. A study found that the ERK-JNK-P38 signaling pathways in autoimmune disease are activated, leading to high levels of downstream JUN and Fos protein.
Subsequently, the two combine to form dimer transcription factor-activating protein 1 (AP-1), which is involved in the occurrence and development of RA. It regulates the transformation of initial T cells into effector T cells, hence regulating the immune response and inflammatory process [55][56][57][58]. Inflammatory response runs throughout RA. e proinflammatory cytokines such as TNF-α, IL-6, and IL-17 play a crucial role in synovial fibroblast inflammation. Studies have demonstrated that TNF-α upregulates the E2F2 expression by activating the nuclear transcription factor NF-κB pathway.
is enhances synovial cell proliferation and synovial tissue thickening leading to joint injury [59]. e serum level of intercellular adhesion molecule-1 (ICAM-1) is high in RA patients. It has been found that ICAM-1 induces inflammatory cell aggregation and promotes synovial cell inflammation once it combines with ligand lymphocyte function antigen-1 (LFA-1) [60].
rough the molecular docking analysis of the effects of the chief active components in HGWD and the core target proteins of RA, the molecular mechanism of the treatment of RA by the active components of TCM can be predicted. is is of great reference significance for subsequent research and development of targeted drugs. e results showed that the main compounds in HGWD have a strong binding force to the core proteins VCAM1, CTNNB1, and JUN. Among them, stigmasterol and stepholidine closely link to VCAM1 through hydrogen bonds and hydrophobic forces. Beta-carotene and quercetin closely link to CTNNB1 through hydrogen bonds and hydrophobic interactions. Stigmasterol closely links to JUN through hydrogen bonds and hydrophobic interactions. Beta-carotene closely links to JUN only through hydrophobic interactions. Based on this, it can be concluded that stigmasterol, stepholidine, beta-carotene, and quercetin are the key compounds in the treatment of RA.

Conclusion
Our study systematically elucidated the mechanisms of action and molecular targets of HGWD against RA via the network pharmacology approach. Quercetin, kaempferol, and betasitosterol regulated most of the targets related to RA. e HGWD can regulate gene function through their related pathways including TNF signaling pathway and IL-17 signaling pathway. e key target genes in the gene-pathway network of HGWD against RA were JUN, FOS, CCND1, IL6, E2F2, and ICAM1. Furthermore, according to molecular docking analysis, important compounds such as stepholidine, stigmasterol, beta-carotene, quercetin, and the core protein CTNNB1, VCAM1, and JUN all have good binding ability. e network pharmacology is a promising suitable approach for the study of TCM formulations.

Data Availability
e data supporting the findings of this study are available from the corresponding author upon request.

Disclosure
Wei Liu and Yihua Fan are the co-first authors.

Conflicts of Interest
e authors declare that they have no conflicts of interest.

Authors' Contributions
Wei Liu and Yihua Fan contributed equally to this work.