A Preliminary Inquiry Into the Potential Mechanism of Huang-Lian-Jie-Du Decoction in Treating Rheumatoid Arthritis via Network Pharmacology and Molecular Docking

Huang-Lian-Jie-Du decoction (HLJDD) has been widely applied to treat inflammation-associated diseases for thousands of years in China. However, the concrete molecular mechanism of HLJDD in the treatment of rheumatoid arthritis (RA) remains unclear. In this work, network pharmacology and molecular docking were applied to preliminarily analyze the potential active ingredients, drug targets, and related pathways of HLJDD on treating RA. A total of 102 active compounds with corresponding 189 targets were identified from HLJDD, and 41 common targets were further identified by intersecting with RA-related targets. Functional enrichment analysis was performed to screen the biological pathways associated with RA. Ten hub targets were further identified through constructing the protein–protein interaction (PPI) network of common targets, which were mainly enriched in the interleukin-17 (IL-17) signaling pathway, tumor necrosis factor (TNF) signaling pathway, and Toll-like receptor signaling pathway. Furthermore, a complex botanical drugs-ingredients-hub-targets-disease network was successfully constructed. The molecular docking results exhibited that these vital ingredients of HLJDD had a stable binding to the hub targets. Among these ingredients, quercetin (MOL000098) was the most common molecule with stable binding to all the targets, and PTGS2 was considered the most important target with multiple regulations by the most active ingredients. In vitro, we successfully validated the inhibitory role of quercetin in the cellular proliferation of human RA fibroblast-like synoviocyte cell line (MH7A cells). These findings indicated that the potential mechanisms of HLJDD for RA treatment might be attributed to inhibiting the immune-inflammatory response, reducing the release of chemokines, and alleviating the destruction of extracellular matrix (ECM) in the synovial compartment.


INTRODUCTION
Rheumatoid arthritis (RA) is a chronic inflammatory autoimmune disease characterized by synovitis, pannus formation, and cartilage erosion, eventually leading to progressive joint dysfunction, deformity, and increased mortality risk (Smolen et al., 2016). The prevalence rate of RA is approximately 1% in the world's population, with higher morbidity among elderly women (Abbasi et al., 2019). So far, the treatment of RA was mainly dependent on the appropriate combination of physical, pharmacological, and surgical approaches, of which nonsteroidal anti-inflammatory drugs (NSAIDs) and disease-modifying antirheumatic drugs (DMARDs) were widely used to lighten inflammation and improve joint function (Wilsdon and Hill, 2017;Conigliaro et al., 2019). However, these drugs may result in some adverse effects, including infection, gastrotract reaction, skin rashes, bone marrow suppression, and liver and kidney toxicity (Chen et al., 2019).
Recent advancements in pharmacy have attracted increasing attention on the use of traditional Chinese medicine (TCM) for the treatment of RA due to its significant therapeutic efficacy and mild side effects . As a well-known classic TCM formula for clearing heat and detoxification, Huang-Lian-Jie-Du decoction (HLJDD) is composed of Scutellaria baicalensis Georgi (Huangqin, HQ), Phellodendron amurense Rupr. (Huangbai, HB), Coptis chinensis Franch. (Huanglian, HL), and Gardenia jasminoides J.Ellis (Zhizi, ZZ) with a ratio of 2:2:3:3 (Supplementary Table S1). HLJDD has been widely applied to the treatment of inflammationassociated diseases including hepatitis (Wei et al., 2016), pneumonia , Alzheimer's disease (AD) (Gu et al., 2018), inflammatory bowel disease (IBD) (Yuan et al., 2019), and RA (Hu et al., 2013). Furthermore, numerous studies have shown that HLJDD could significantly reduce the levels of inflammatory cytokines such as TNF-α and inflammatory mediators such as prostaglandin E2 (PGE2) to realize the anti-inflammatory efficacy (Fang et al., 2004;Lu et al., 2011). In vitro, Sung et al. (2012) demonstrated that quercetin could significantly inhibit unstimulated and IL-1β induced proliferation of rheumatoid synovial fibroblasts (RASFs) and the messenger ribonucleic acid and protein expression of MMP-1, 3, COX-2. In animal models, Hu et al. also successfully demonstrated the treatments of HLJDT on collagen-induced arthritis in rats (Hu et al., 2013). During the period of COVID-19, HLJDD was also found to play a therapeutic role in COVID-19 through regulating multiple signaling pathways based on targeting genes such as IL6, IL10, MMP9, NOS2, VEGF, and TGFβ1 (Liu et al., 2021). Nevertheless, the concrete molecular mechanism of HLJDD in the treatment of RA was still unclear.
As a novel, promising, and cost-effective drug research approach based on bioinformatics and pharmacology, network pharmacology has been widely used in studies on the molecular mechanism of drugs with molecular docking. In previous studies, Yang et al. (2013) identified berberine, baicalin, and geniposide as three major ingredients of HLJDD and further demonstrated their anti-inflammatory effects through inhibiting NF-κB and MAPK signaling pathways in a dose-dependent manner. In addition, using network pharmacology and metabolomics analysis, Qu et al. found the antidepressant effects of HLJDD through candicine, berberine, tetrahydroberberine, and baicalein to regulate SLC6A4 and MAOA in tryptophan metabolism (Qu et al., 2021). In Li's study, the network pharmacology and molecular docking method was also well conducted to decode the potential mechanism of HLJDD in treating pneumonia and construct a complex botanical drugs-ingredients-targets-disease network . These studies provided a convincing possibility that network pharmacology and molecular docking analysis can be effectively applied to investigate potential therapeutic targets of HLJDD and facilitate a deep understanding of its underlying mechanism in RA.
In the present study, we investigated the potential targets and pathways of HLJDD in treating RA and constructed the botanical drugs-ingredients-targets-disease network through network pharmacology analysis and molecular docking validation. We were convinced that our results would help in illuminating HLJDD's possible mechanism of the treatment in RA, thus improving the curative effect and prognosis for RA. The workflow of the whole process of our study is shown in Figure 1.

Screening Active Ingredients and Potential Targets in HLJDD
In this study, the Traditional Chinese Medicine Systems Pharmacology (TCMSP) (Ru et al., 2014) database was utilized to screen active ingredients of HLJDD by using "Huangbai," "Huanglian," "Huangqin," and "Zhizi" as keywords with the filtration criteria of bioavailability (OB) ≥30% and druglikeness (DL) ≥0.18, respectively (Song et al., 2018). The corresponding molecular structure, structural information, and PubChem CID of those active ingredients were further obtained from the "PubChem" online tool (Wang et al., 2017). Subsequently, using the above structural information of ingredients, we utilized the "BATMEN-TCM" online tools to identify potential targets in HLJDD with a score ≥20 and adjusted p ≤ 0.05 (Liu et al., 2016).

Identification of Target Genes in Rheumatoid Arthritis
score ≥0.4. Then, the PPI network was constructed and visualized by the Cytoscape software (Shannon et al., 2003), and hub targets were detected according to levers of degree value (the number of interactions for each node) in the PPI network. In the further cluster analysis, the plug-in MCODE (Deng et al., 2019) was used to identify significant modules using the filter conditions: k core FIGURE 1 | The summary and description of the study workflow in the potential mechanism of HLJDD in treating RA. The active ingredients and corresponding targets of HLJDD were obtained through TCMSP and BATMEN-TCM database. The RA-related targets were downloaded from four different databases, and the common targets were identified by intersecting ingredient targets and RA-related targets. GO and KEGG enrichment analysis was conducted, and PPI network with cytoHubba plug-in was used to select hub targets in the common targets. Finally, the complex botanical drugs-ingredients-hub-targets-disease network was constructed and validated by molecular docking and normalized expression of RA-related target genes from GEO datasets.

GO and KEGG Functional Enrichment Analysis
To investigate the biological function of potential targets in RA, we performed gene ontology (GO) analysis, including biological processes (BP), molecular functions (MF), and cellular components (CC), and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis. These analyses were conducted through the online tool "g:Profiler" (Raudvere et al., 2019) with Bonferroni correction and adjusted p ≤ 0:05. Moreover, the ClueGo (Bindea et al., 2009) and CluePedia (Bindea et al., 2013) plug-in of Cytoscape software were further applied to visualize the network of significant pathways and their corresponding target genes.

Identification of Hub Targets and Critical Network of HLJDD on RA
The common targets between HLJDD and RA were identified through taking the intersection, and the PPI network of these targets was also constructed and visualized in the Cytoscape software. Subsequently, we used plug-in cytoHubba (Chin et al., 2014) to explore important nodes in the above PPI network by 12 topological algorithms including MCC, MNC, DMNC, BottleNeck, EcCentricity, EPC, Radiality, Betweenness, Closeness, Clustering Coefficient, Degree, and Stress. The top 10 nodes identified by each topological algorithm were chosen to identify the shared genes more than six ways as the most pivotal hub genes in the network (Luo et al., 2020). Finally, a complex botanical drugsingredients-hub-targets-disease network was further constructed and visualized in the Cytoscape.

Visualization and Validation of Molecular Docking
To validate the binding of hub targets and active ingredients in HLJDD, the 3D molecular conformations of ingredients were retrieved from the PubChem Compound database. The crystal structures of target proteins were obtained from the RCSB Protein Data Bank (PDB) database (Burley et al., 2021). AutoDock 4.2.6 tool was used to remove the redundant structure, ligands, and water molecules and add polar hydrogen atoms and partial charges into protein crystal structures before running molecular docking (Li H. et al., 2018). The visualization of molecular docking was exhibited by PyMol software (Seeliger and de Groot, 2010). The transcriptomic dataset of peripheral blood mononuclear cells (PBMCs) (GSE17755, including 112 RA and 45 controls) was downloaded from Gene Expression.
Omnibus (GEO) database to validate the normalized expression of RA-related target genes. Values were presented as the mean ± standard deviation (SD), and the comparison between groups was performed using the Wilcoxon test by the "ggpubr" R package. p < 0.05 was considered statistically significant.

Cell Culture and Treatment
Human RA fibroblast-like synoviocyte cell line (MH7A cells) was purchased from the Riken cell bank (Ibaraki, Japan), and cells were maintained in DMEM supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin solution (HyClone, Shanghai). Subsequently, the MH7A cells were plated into a sixwell plate (5 × 10 3 cells/well) and treated with different concentrations of quercetin (0, 25, 50, 75, and 100 μM; Sigma, St Louis, MO, United States) in 37°C incubator with 5% CO 2 atmosphere for 24 h. Meanwhile, the culture medium with 10% DMEM and 0.05% DMSO but without MH7A cells was set as control groups and also plated in the same condition for 24 h.

Cell Proliferation Assay
Equal numbers of the above MH7A cells were plated in a 96-well plate (5 × 10 3 cells/well), and Cell Counting Kit-8 (CCK-8) (Dojindo, Tokyo, Japan) (10 μL/well) was applied to measure the cell proliferation at 24 h. Finally, the imaging of cell proliferation was recorded, and the absorbance value (OD value; A 450nm ) of MH7A cells was further detected using a Microplate Reader (Thermo Fisher, MA). Based on the OD value of different quercetin concentrations, we successfully drew the drug-concentration curve to identify the halfmaximal inhibitory concentration (IC50) value of quercetin in RA.

Statistical Analysis
All statistical analyses were performed in R software (version 4.0.1, https://www.r-project.org/). Continuous variables were presented as mean ± standard deviation (SD), and Wilcoxon signed-rank test was used to compare continuous variables. The two-tailed P value <0.05 was considered statistically significant.

Identification of Potential Targets-Active Ingredients Network
Based on the criteria of OB ≥ 30% and DL ≥ 0.18, a total of 102 active ingredients were predicted in HLJDD through the TCMSP database, including 37 from HB, 14 from HL, 36 from HQ, and 15 from ZZ (Supplementary Figure 1A). Moreover, we also screened a total of 1008 corresponding targets of HLJDD using the TCMSP and BATMEN-TCM, of which 271 were from HB, 216 from HL, 226 from HQ, and 295 from ZZ (Supplementary Table S3). Combined with the 102 active ingredients, 189 targets of these ingredients were further identified, and the complex targets-active ingredients network is constructed in Figure 2A. Notably, by setting the degree value ≥ 45, the top 10 active ingredients were identified in the network, including quercetin (MOL000098), beta-sitosterol (MOL000358), stigmasterol (MOL000449), kaempferol (MOL000422), wogonin (MOL000173), baicalein (MOL002714), palmatine (MOL000785), berberine (MOL001454), coptisine Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 740266 (MOL001458), and 5-hydroxy-7-methoxy-2-(3,4,5trimethoxyphenyl) chromone (MOL003095). The detailed characteristics of these active ingredients and targets are shown in Supplementary Table S2 and Supplementary Table  S3. To further identify the changes of gene expression in the progress of RA, four databases, namely, OMIM, GenBank, PharmGKB, and TTD, were screened to obtain RA-related target genes, and a total of 440 targets of RA were obtained (Supplementary Figure 1B, Supplementary Table S4). Subsequently, comparing the target genes of HLJDD and RA, The active ingredients-targets network of HLJDD. (C) The PPI network of all RA-related targets. The cyan nodes represent the 10 hub targets obtained from the subsequent cytoHubba analysis. (D) The PPI network of 39 common targets between RA-related and HLJDD; The PPI networks of 10 vital targets were identified using Cytohubba plug-in with approval of more than 6 methods.
Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 740266 5 we found that RA shared 41 common targets with that of active ingredients (Figure 2A).

PPI Network and Cluster Analysis of Disease Targets in RA
To identify central target genes of RA in the physical interaction network and provide potential clues for the possible pathogenic mechanisms of RA, we successfully constructed an interconnected PPI network of 440 RA-related targets with 394 nodes and 9391 edges ( Figure 2C) based on the STRING database. Using the MCODE plug-in in Cytoscape with MCODE score greater than 5, we further identified five clusters with close interconnection (Supplementary Figure 1C, Table 1). The blue nodes represented the significant targets with high degree values in the network, which might play essential roles in the progression of RA.

Functional Enrichment Analysis of Disease Targets in RA
In order to further interpret biological processes associated with the RA-related targets, GO and KEGG enrichment were performed based on the 440 target genes. It turned out that these target genes were enriched in 1237 BPs, 81 MFs, 67 CCs, and 76 KEGG pathways (Supplementary Table S5). The results showed that the biological processes of RA were mainly related to the immune response process, such as "immune system process," "response to organic substance," "cellular response to cytokine stimulus," and "inflammatory response." The molecular functions of RA might be related to "signaling receptor binding," "molecular transducer activity," and "signaling receptor activity." The cellular components were principally associated with "cell periphery," "extracellular region," and "plasma membrane" ( Figure 3A). The KEGG pathway analysis also revealed these targets were significantly enriched in immune-activated associated diseases and pathways, including "rheumatoid arthritis," "inflammatory bowel disease," "cytokine-cytokine receptor interaction," "TNF signaling pathway," "NF-kappa B signaling pathway," and "Th17 cell differentiation" (Figure 3B). In addition, the concrete network of the KEGG pathways indicated that massive inflammatory cytokines participated in the progression of RA ( Figure 3C).

Identification of Hub Genes and Critical Network of HLJDD on RA
In order to identify the hub genes of HLJDD in treating RA, we also imported the 41 common targets into the STRING database and subsequently constructed their PPI network ( Figure 2D). The cytoHubba plug-in was further applied to calculate the significance of targets with 12 methods, and we found that there were 10 hub genes shared with more than 6 topological analyses (MMP9, IL1β, IFNG, IL10, ICAM1, CCL2, PTGS2, IL6, CXCL8, and JUN) by selecting top 10 genes of each method ( Table 2). The PPI network of the hub genes included 10 nodes and 45 edges, with an average node degree of 9 and P value of 3.22e-14 ( Figure 2D).

Validation of Hub Genes by Molecular Docking and mRNA Expression
To further evaluate and validate the effective binding of targets and active ingredients, ten hub genes and nine pivotal ingredients were used to perform molecular docking as receptors and ligands, respectively. The stability and strength of the binding were based on the binding energy. The lower the binding energy value was, the more stable the binding was. Setting the binding energy value < −5.0 kcal/mol as the threshold, we successfully demonstrated the ligands binding well to the target receptor with the lowest energy ( Table 3). Among these ingredients, quercetin (MOL000098) was the most significant molecule with stable binding to all the targets, and PTGS2 was considered the most important target with multiple regulations by the most active ingredients ( Figure 5). Normalized expression of 10 hub targets is also validated in GEO datasets, and it revealed that IL6, IL1β, JUN, CXCL8, MMP9, PTGS2, ICAM1, and CCL2 were significantly upregulated. In contrast, IFNG and IL10 were significantly downregulated in the PBMC of RA patients compared to controls (Figures 6A-J).

Quercetin Inhibited Cell Proliferation in MH7A Cells
Initially, we estimated the proliferation of MH7A cells after the treatment of series quercetin concentrations (0, 25, 50, 75, and 100 μM) and found the inhibitory effectiveness was generally increased with the enhancement of quercetin ( Figure 6K). However, the growth did not increase speedily when the concentration was more than 50 μM ( Figure 6L), and the IC50 value of quercetin was identified as 38.02 μM according to the drug-concentration curve ( Figure 6M).

DISCUSSION
As one of the most common systematic autoimmune disorders, RA has a high proportion of incidence and joint deformity with high disease activity. HLJDD is a classical formula from TCM and has been widely used to treat inflammation-related diseases for thousands of years, and recent preclinical studies also demonstrated its effective treatment in RA (Hu et al., 2013). However, its concrete mechanism of the clinical effects remains unclear. Therefore, we investigated the potential mechanisms of the therapeutic effect of HLJDD in RA employing network pharmacology and molecular docking analysis.
Moreover, 10 hub targets (MMP9, IL1β, IFNG, IL10, ICAM1, CCL2, PTGS2, IL6, CXCL8, and JUN) of HLJDD in treating RA were identified through constructing the PPI network of the 41 common targets from HLJDD and RA. These targets were significantly enriched in the biological process of the cellular inflammatory response to different substances such as organic substance and chemical and cytokine stimulus and associated with IL-17 signaling pathway (CCL2, CXCL10, CXCL8, IFN-γ, IL1β, IL6, JUN, MMP9, PTGS2), TNF signaling pathway (CCL2, CXCL10, ICAM1, IL1β, IL6, JUN, MMP9, PTGS2), rheumatoid arthritis pathway (CCL2, CXCL8, ICAM1, IFNG, IL1β, IL6, JUN), and Toll-like receptor signaling pathway (CXCL10,CXCL8,IL1β,IL6,JUN). Accumulating studies showed that IL-17 and IL-17-producing T helper (Th17) cells played a critical role during the development and progression of RA ) and Lee's study exhibited that IL-17 upregulated the expression of TLR3 by STAT3 pathways in fibroblast-like synoviocytes of RA (Lee et al., 2014). Tumor necrosis factor (TNF) signaling pathway has been reported to mediate leukocyte recruitment and inflammatory responses, and various anti-TNF inhibitors also acquired prominent potent effects on RA . In addition, Toll-like receptor (TLRs) and its corresponding downstream signaling pathways, such as Wnt, NF-κB, and MAPK signaling pathways, have been demonstrated in the biological process of synovial inflammation and bone remodeling of RA (Andreakos et al., 2005). Interestingly, low IL-10 and IFN-γ levels were found in RA, and this contradictory finding might be explained by published studies. Seung et al.  found IFN-γ regulated inflammatory cell death by targeting necroptosis, and IFN-γ deficiency increased the number of Th17 cells and upregulated the expression levels of IL-17 and TNF-α using experimental mice of collagen-induced arthritis (CIA). In addition, an authoritative study declared that IL-10 was also considered a highly promising treatment target for RA attributing to its unique capacity to inhibit cellular immunity and deactivate macrophages through downregulating the production of multiple proinflammatory cytokines, including IL1 and TNF-α (St Clair, 1999). Overall, the results suggested that complicated immuneassociated signaling pathways were activated in the progression of RA, and inhibiting associated targets would alleviate the disease activity and prognosis for RA.
As important inflammatory regulators, IL6 and IL1β belong to the family of cytokines and could stimulate the production of neutrophils and megakaryocytes, guide immune cell differentiation, and upregulate the Th17/Treg balance, which was pathologically associated with the development of chronic inflammatory diseases (Tanaka et al., 2018;Gomes da Silva et al., 2020). In addition, the clinical practice had authenticated that the therapeutic monoclonal antibody against IL6, such as tozumab, could acquire significant curative effects in the treatment of RA (Humby et al., 2021). The activity of macrophages in the synovial compartment was regulated by substantial inflammatory factors, especially chemokines families, which included four groups: the C, CC, CXC, and CX3C chemokine receptor families (Zlotnik and Yoshie, 2000). Interestingly, three hub targets belong to the CC (CCL2) and CXC (CXCL8, CXCL10) subfamilies, of which CXCL8 was produced constitutively by macrophages in the synovial tissues and demonstrated to induce synovial inflammation in the animal model (Endo et al., 1994), and the main function of CCL2 was responsible for the recruitment of macrophages (Hachicha et al., 1995). Interferon-gamma (IFN-γ) was predominantly produced by T lymphocytes and natural killer (NK) cells and played an essential role in autoimmune disorders through activating multiple immune-related pathways, especially the JAK-STAT pathway (Simon et al., 2021). In Masaru's study, a high level of IFN-γ was found associated with the severity and poor prognosis of RA, and the expression level of IFN-γ was significantly decreased after the treatment of JAK inhibitors (Kato, 2020). Matrix metalloproteinase-9 (MMP9) played a critical role in degrading components of extracellular matrix (ECM) and leukocyte migration in rheumatoid arthritis  patients (Alwan and Ghali, 2021). Prostaglandin-endoperoxide synthase 2 (PTGS2), also called COX-2, has been recognized as a significant proinflammatory target of RA, and COX-2 inhibitor (NSAIDs) was the preferred anti-inflammatory drug to control pain and stiffness (Crofford, 2013). Notably, Nicole et al. found that the AP-1 Transcription Factor (JUN) could directly activate proinflammatory factor COX-2 in macrophages to promote arthritis in rat models (Hannemann et al., 2017). Summarily, the potential mechanisms of HLJDD in the treatment of RA might be attributed to the following aspects: inhibiting the immune inflammatory response, reducing the release of chemokines, and alleviating the destruction of ECM in the synovial compartment.
However, there still are some limitations in our study. On the one hand, although network pharmacology provides a comprehensive approach based on the combination of pharmacology, biochemics, bioinformatics, and network biology, the data of active ingredients and targets were collectively obtained from multi-public databases, and still selection bias exists. On the other hand, in this study, we just validated quercetin's therapeutic role in inhibiting the proliferation of RA cells in vitro, but the ingredients-targets network and concrete molecular mechanism of HLJDD in the treatment of RA remain to be verified by in-depth in vivo and in vitro studies. Finally, the curative effect of quercetin might be non-specific due to artificial synthesis in this experiment and need more validation using extractive compounds in the future.

CONCLUSION
In this study, we successfully screened and obtained the pivotal ingredients of HLJDD in treating RA, and 10 hub targets were identified to construct the botanical drugs-ingredients-targetsdisease network. Functional enrichment analysis suggested HLJDD could regulate related targets through associated pathways, including TNF signaling pathway, Toll-like receptor signaling pathway, and IL-17 signaling pathway. Furthermore, based on the molecular docking analysis, important compounds such as quercetin, beta-sitosterol, wogonin, kaempferol, coptisine, and stigmasterol exhibited good binding ability with these core proteins. In conclusion, our study systematically expounded the mechanisms and molecular targets of HLJDD in the treatment of RA through the network pharmacology and molecular docking approach. The therapeutic action of HLJDD in RA was mainly by inhibiting the immune inflammatory response, reducing the release of chemokines, and alleviating the destruction of ECM in the synovial compartment.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
CL and JP contributed to the data analysis and drafting of the manuscript. CX contributed to data acquisition, figures presentation, and revision of the manuscript. ZJ and XC contributed to the design of the study. All authors contributed to the article and approved the submitted version.

FUNDING
This study was supported by the General Scientific Projects of the Zhejiang Education Department (Y202147905).