Abstract
Recurrent aphthous stomatitis (RAS) are complex inflammatory diseases caused by multi-factors, which severely impact patient quality of life. However, there is still no effective treatment method for RAS without side effects. Traditionally, Cortex Phellodendri known as “Huang Bai” was used to treat RAS for antibacterial and anti-inflammatory properties in China. Network pharmacology methods and bioinformatics analysis were utilized to search and fish incorporating target. Network analysis and silico validation were used to discover the pharmacological mechanisms of “Huang Bai” for the treatment of RAS. A total of 25 active ingredients in HB, 200 drug targets, and 578 differentially expressed genes (DEGs) between Recurrent aphthous stomatitis and normal samples were obtained. The Gene Ontology enrichment analysis revealed that the immune response was the most significantly enriched term within the DEGs. The KEGG pathway analysis identified 60 significant pathways, most of which involved in the inhibition of inflammation and regulation of immunological response. The functions are dependent on a multi-pathway, particularly the TNF signaling pathway and the HIF-1 signaling pathway. We identified six hub genes in the PPI network, most of which were validated as highly expressed in oral ulcers by DiseaseMeth databases. In addition, molecular docking displayed that the primary molecule combined well with the key targets. “Huang Bai” contains potential anti-RAS active compounds. This study reflects the multi-component multi-target multi-pathway action characteristics of “Huang Bai.” Our study provides potential biomarkers or treatment targets for further research.
Keywords
Introduction
Recurrent aphthous stomatitis (RAS), also called recurrent oral ulcer (ROU) or recurrent aphthous ulcer (RAU), is a chronic inflammatory oral ulcer disease, which is the most principal clinical disease in the oral cavity. 1 The typical symptom is a single or several painful erosions that recur on the lips, the tip of the tongue, edge of the tongue, cheek, and other areas. The disease can affect people of all ages, races, and regions. 2 RAS usually first appears in childhood or adolescence, frequently recurring and persistent lifetime morbidity, with incidence rates varying from 5 to 60% in different population. 3 In children, the prevalence of the disease is significantly influenced by one or both parents with the presence of RAS. 4 Children with parents with positive RAS have a 90% chance of developing RAS, while children with parents with negative RAS have a 20% chance. 5 A study of genetic testing for single-gene disease indicates that children with RAS have Mutations in TNFAIP3, DNASE1L3, PTPN22, and STAT1. 6 Therefore, researchers consider that individuals who are genetically predisposed to RAS are susceptible to RAS. In severe cases, RAS can lead to difficulty in eating and speech, causing pain to the patient. Scholars have explored multiple strategies to prevent RAS, including Laser therapy, probiotics and vitamins, mouthwash treatments, and immunotherapy. However, these therapies have adverse side effects and high recurrence rates. Effective treatment options with few side effects are urgently needed to treat RAS.
“Huang Bai”(HB) is one of the most frequently used traditional Chinese medicine, also known as Cortex Phellodendri. It is the dried bark of Phellodendron chinense Schneid. or Phellodendron amurense Rupr. HB is bitter, cold in nature, and functions on cleansing damp-heat and drying up the dampness, relieving fire and detoxification, remove swelling and decay. 7 HB has a combination of immunomodulatory, anti-inflammatory, antibacterial, anti-hypertensive, antioxidant, and antipyretic therapeutic effects.8–10 Modern research has shown that HB contains numerous chemical components, including alkaloids and terpenoids. Both of these compounds are effective infections and neurological diseases. 11 Some studies have proven the effects of HB against RAS. Yang et al. reported that HB promotes the healing of phenol-induced oral ulcer in rats by inhibiting the expression of TNF-α. 12 Shao et al. found berberine, as well as other compounds, could increase the expression of Sirt1 and Foxp3 and suppress the expression and acetylation of NF-κB in oral mucosa cells. 13 A randomized controlled trail found that HB has better effects than bingpeng powder in the ulcer healing in patients with RAS. 14 Although several studies have reported on the anti-RAS activity of HB, the mechanisms of HB contributing to therapeutic effects remain unclear.13,15
Network pharmacology is a useful bioinformatics tool to discover candidate targets, functions, and mechanisms of bioactive components to treat disease. In the research, we aimed to explore the potential therapeutic targets and mechanisms of HB by utilizing a network pharmacology technology.
Materials and methods
The data were obtained from public databases, so the ethics committee approval was unnecessary.
Active ingredients selection
In the TCMSP(https://tcmspw.com/tcmsp.php), using “Huang Bai” as the keyword, we determined the oral bioavailability (OB) ≥ 30% and drug likeness (DL) ≥ 0.18 as the selection criteria. “Oral bioavailability” (OB) refers to the percentage of unmodified drugs that enter the circulatory system after oral administration.16,17 OB is a useful indicator for objective evaluation of the internal quality of drugs. The higher OB of an ingredient, the higher is the likelihood that it can be used clinically. The term “drug likeness” (DL) refers to the similarity between ingredients and known drugs. Most ingredients in Chinese formulations have poor pharmacologic properties, so they can’t bind to specific protein targets on cells efficaciously. Several researchers18,19 have recommended that molecules with OB ≥ 30% or DL ≥ 0.18 are considered to have better pharmacologic effects and are selected as candidate ingredients for further analyses. Protein sequences of these drug targets were normalized to official gene symbols using the UniProt database (https://www.uniprot.org/).
Acquisition of differentially expressed genes
GEO database (http://www.ncbi.nlm.nih.gov/geo), which is a publically available database of gene/microarray profiles, was used to screen relative targets of RAS. The search strategy (“recurrent aphthous stomatitis” [MeSH Terms] OR “recurrent aphthous stomatitis” [All Fields]) AND (“Homo sapiens” [Organism] was adopted to search differentially expressed genes (DEGs). We downloaded the microarray dataset GSE37265 containing 14 normal samples and 14 RAS samples from the GEO. The datasets were quantile-standardized by using R script. Each gene expression was calculated based on the false discovery rate (FDR; p < .05) using the Benjamini–Hochberg method and t-test. DEGs with the threshold criterion of |log FC| >1 and p < .05 in samples were screened using the limma V3.42.0 package of the R software (version 4.1.2). The WGCNA package (version: 1.70–3) was employed to construct a gene co-expression network using a variant set of genes. The analysis was performed based on the package instructions. The minimum number of genes in each module was set to 30.
Construction of a network and analyses
We obtained the intersection of HB and RAS targets using the Venn diagram of Venny 2.1.0 (https://bioinfogp.cnb.csic.es/tools/venny/). The shared target genes were then imported into the STRING database 20 (https://string-db.org/) with the species defined as human as well as the confidence score >0.4. Cytoscape software 3.8.2 was used to visualize the PPI network for common targets. In the component-key target network, nodes represent components and targets, and edges represent the interactions between nodes. The targets without interaction were removed from the network. Afterward, the topological parameters of each node in the net were analyzed by using Network Analyzer tool Cytohubba pluginin in Cytoscape. Cytohubba is a plug-in of Cytoscape, was employed to study essential nodes in the network with 11 methods (MCC exhibits a satisfied comparative performance) The topological parameters including degree, betweenness centrality, and closeness centrality were used to evaluate the significance of the nodes in networks. The nodes with “Degree” values 2-fold greater than all the network nodes’ median value, “Closeness centrality” and “Betweenness centrality” values above all the network nodes’ median value were identified as the key targets.
GO function and KEGG pathway enrichment analysis
For Biological Processes (BP), Cellular Components (CC), Molecular Functions (MF), and KEGG Pathways, the R software analysis package’s clusterProfiler, Digest, and Goplot programs, as well as Cytoscape’s CluePledia plug-in and Cluego plug-in, were utilized. The threshold p < .05 was set. The enrichment score (ES) represents the degree to which the set S is overexpressed at the top or bottom of the sorted list L. The higher the enrichment score, the more relevant the pathway is for disease progression. The more genes contained in a pathway, the larger the bubble represented by that pathway. The top 10 GO results with the lowest p adjust values were displayed in the form of bar chart and cnet plot. The top 15 pathways were mapped to the oral ulcer disease related pathways in the CTD database (http://ctdbase.org/) to get common pathways.
Molecular docking of core targets to HB components
The molecular docking performance was operated by site-features directed docking for the study of the interactions between receptor and ligand. 21 The protein structure of hub genes was obtained from the PDB database (http://www.rcsb.org), The molecular structure of key components ligand files in mol2 formats were obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov/). CB-Dock (http://cao.labshare.cn/cb-dock/) was used for molecular docking. The CB-Dock server, developed by Dr Liu’s research team, is a user-friendly blind docking network server. Autodock Vina is used for docking, and it employs a unique curvature-based cavity detection technique. 22 This method had a success rate of more than 70%, outperforming even the most advanced blind docking tools. The atom positional root-mean-square deviation (RMSD) was developed as a metric for structural comparison by crystallographers. It has become a standard tool for comparing pairs of structures, characterizing ensemble. 23 RMSD ≤2 was considered the threshold for conformation of the ligand after docking to match the conformation of the original ligand.
Results
Composite ingredients of HB
The details of active ingredients of HB.
Identification of DEGs of RAS and HB-RAS intersection targets
The gene expression profiles GSE3726537.
24
were obtained from the GEO database. The platform for GSE37265 is GPL570, which includes 14 samples of RAS patients, and 14 samples from normal tissue. Data table header descriptions and series matrix files of GSE37265 were downloaded. The expression dataset is corrected as shown in Figure 1 with good consistency. Relative log expression (RLE) reflects the consistency of parallel experiments. Filtration based on principal component analysis (PCA) includes all samples. The first two principal components distinguished RAS and normal samples well, accounting for 23.36% (first component) and 18.65% (second component) of the observed difference (Figure 2(b)). The gene expression profiles of these 28 samples were used in subsequent analyses. We selected genes with variances greater than all quartiles of variance for further analysis (20,486 genes). The Pearson’s correlation coefficient was conducted to cluster the samples in GSE37265. After removing outliers, we drew a sample clustering tree. We set the soft threshold to 10 (R2 = 0.98) to construct a scale-free network. WGCNA analysis identified seven modules based on average hierarchical clustering and dynamic tree clipping (Figure 2(c) and (d)). The blue module was highly related to age of patients. A total of 578 DEGs were extracted after standardization of the datasets based on the defined criteria. The DEGs are shown in the volcano plots and the heatmaps whose clustering was performed with Euclidean distance. GSE37265 includes 535 up-regulated genes and 43 down-regulated genes (Figure 2, Supplementary Table 1). Analysis of dataset. (a) Box plot of RLE. The result indicates good quality control in normal and RAS groups. (b) Principal component analysis. Blue represents normal samples, red represents RAS samples. The points of the RAS group and the normal group are far apart in space, indicating that the principal components are different. (c) Heatmap of the correlation between the module eigengenes and clinical traits of RAS. (d) Clustering dendrogram of DEGs depending on the dissimilarity measure (1-TOM). Each branch represents a gene, and each color represents a co-expression module. The color band shows the results obtained from the automatic single-block analysis. The x-axis represents soft-thresholding power, the y-axis represents mean connectivity (degree). TOM, topological overlap matrix. Volcano plot and heatmap of DEGs. (a) Volcano plot of DEGs, Data points in left represent down-regulated, and right represent up-regulated genes. The differences are set as |log FC|>1 and p < .05. (b) Heatmap of top 10 up-regulated DEGs and top 10 down-regulated DEGs (red: high expression; blue: low expression).

We obtained 200 potential targets of active ingredients of HB. The co-expressed genes were integrated using the Venn Diagram online tool, including 27 genes in HB and RAS (Figure 3 and Table 2). Venn diagram of common differentially expressed genes of RAS in HB. Twenty-seven DEGs were screened out. Differentially expressed genes in HB and RAS.
Construction of HB-RAS interaction networks
PPI network of common targets for HB and RAS
We obtain the network of intersection targets by the STRING database with 27 targets, saved them as TSV format files, used Cytoscape software to draw network diagrams for the core targets. The hub PPI network comprises 27 nodes and 215 edges with the maximum degree of freedom being 24, the minimum degree of freedom is 1. Eight nodes with higher value (“Degree” > 19, “Betweenness centrality” > 0.01,595 and “Closeness centrality” > 0.73,958) calculated by cytoHubba were screened as the key targets. The top ten degree-ranked targets in the net are PTGS2, IL6, TNF, IL1B, CCL2, ICAM1, CXCL8, IFNG, VCAM1, and CXCL10, indicating their essential roles for the network. We adjusted the color and size of each node according to the size of the “Combined” value of each node, as in Figure 4. PPI network diagram of common targets for HB and RAS. The color scale is from orange to red, the darker the color, the higher the degree of freedom. The size of circles represents nodes degree value.
Network of drug-active components-key targets
The component-key target network including 73 nodes and 324 edges was constructed by using Cytoscape. We found that many targets are hit by different numbers of compounds and other targets, implying the multi-targets characteristics of HB as shown in Figure 5. A network analysis showed that quercetin (HB19, degree = 142), and rutaecarpine (HB7, degree = 102) had the highest number of connections to different targets. Network diagram of drug-active components-key targets. Diamonds, triangles, and circles represent “Huang Bai,” 25 compounds, and 27 intersection targets, respectively. Different colors represent the degree, the color scale is from yellow to red, the darker the color, the higher the degree of freedom. Node size is proportional to the degree of interaction.
Pathway enrichment analysis for key targets
The results of GO analysis show that 1332 of 2015 biological processes, 19 of 53 cell components, 58 of 108 molecular functions enriched for these targets, and 29 target-related pathways were recognized as p < .05 (Figure 6, Supplementary Table 2). Regarding BP, DEGs were significantly enriched in response to lipopolysaccharide, response to molecule of bacterial origin, cellular response to lipopolysaccharide, and cellular response to molecule of bacterial origin. This means lipopolysaccharide and bacterial production play key roles in RAS. For the CC, the DEGs were enriched in membrane raft, membrane microdomain, membrane region, and so on. MF analysis showed that the DEGs were significantly enriched in cytokine activity, cytokine receptor binding, receptor ligand activity, signaling receptor activator activity, and chemokine receptor binding. KEGG pathway-enrichment indicated that DEGs were mainly enriched in the TNF signaling pathway, AGE-RAGE signaling pathway in diabetic complications, Lipid and atherosclerosis, IL-17 signaling pathway, HIF-1 signaling pathway, Coronavirus disease—COVID-19, etc. The top 15 KEGG pathway enrichments are shown in Figure 7. Among these pathways, the HIF-1 signaling pathway and the TNF signaling pathway were associated with oral ulcer disease according to the CTD database. It is well known that a typical symptom of recurrent aphthous stomatitis is the inflammation of the oral cavity. Therefore, the mechanisms of HB’s anti-inflammatory function were investigated. Among the top 15 pathways, several pathways are associated with inflammation, and the TNF signaling pathway is the most distinct one. Actually, the immunologic disturbances play a crucial role in the etiopathogenesis of RAS. HIF-1 signaling pathway is an important regulator of the immune system. Increased serum NO level in HIF-1 pathway is thought to take part in the pathogenesis of RAS.
25
Hence, it is meaningful to identify the impact of HB on immune function. Gene Ontology (GO) enrichment analysis of the key targets. ( Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis of key genes. (a) Bubble chart shows enrichment of key genes in the top 15 signaling pathways. The x-axis label represents the gene ratio and the y-axis label represents pathway. Pathway enrichment results at p < .01. (b) Cnet plot shows the top 10 enriched KEGG pathways of DEGs. The yellow points represent the GO categories, the color of the line delivered by a point indicates the category of the point in the legend, the size of a point indicates the number of the genes it includes.

Compound-target docking
Vina score and RMSD of active chemical and key targets.

Ligand-receptor interactions. (a) TNF-rutaecarpine, (b) TNF-quercetin, (c) MMP3-rutaecarpine, (d) MMP3-quercetin.
Discussion
Recurrent aphthous stomatitis is a complex inflammatory disease caused by multiple factors, which severely influences the life of patients. The etiology of RAS is unknown, but some local trauma, systemic, immune, genetic, allergic, nutritional, and microbiological factors are considered to be causative factors. 27 However, there are no effective drugs with low side effects against RAS.
HB has several biological activities against many Gram-positive and some Gram-negative bacteria.28,29 Besides antibacterial activity, HB also has well-established anti-inflammatory agent. Previous studies have shown the effectiveness of HB against RAS.30–32
The chemical composition of HB mainly includes alkaloids, flavonoids, sterols, lactones, phenolic acids, terpenes, phenylpropanoids, and other chemical components using liquid chromatography-tandem mass spectrometry. Alkaloids are the main active ingredients in HB. About 42 alkaloids have been reported in the current studies. In the study, quercetin and rutaecarpine were among the important active compounds of HB. In Chinese herbs, quercetin is a flavonoid with therapeutic properties, including antioxidant and anti-inflammatory activities. 33 Quercetin has a long history of consumption as part of normal human diets. A prospective active control trial confirmed the effectiveness of quercetin in reducing ulcer size and pain in patients with RAS. 34 In two randomized controlled trials, quercetin was shown to be a highly effective and safe adjuvant therapy.35,36 Furthermore, 90% of patients responded that they appreciated the ease of application when using the topical quercetin, and they did not object to its consistency or taste. 35 It may be promising to formulate a quercetin-based drug treatment for RAS in the future. Rutaecarpine, a major indolequinazoline alkaloid, is the dried fruit of Evodia rutaecarpa (Juss.) Benth 37 . Several studies have revealed that rutaecarpine has extensive pharmacological effects, including anti-inflammatory, anti-oxidation and vasodilation.38–40 Evodiamine inhibits hypoxia induced inflammatory responses (COX-2 and iNOS expression) via suppressing the activation of hypoxia inducible factor-alpha (HIF-α) in macrophages. 41 Many studies have shown that rutaecarpine can reduce pro-inflammatory mediators such as TNF-α, IL-1β, IL-6, and NF-κB to improve the symptoms of inflammation and oxidative stress.42–44 Taken together, we believe that these compounds play important roles in the healing process of RAS.
In the present study, integrated bioinformatics methods assisted in analyzing how the expression of critical genes changed to reveal potential RAS pathways based on GEO datasets “GSE37265.” We identified 578 DEGs, including 535 up-regulated DEGs and 43 down-regulated DEG. Further, using the network pharmacology approach, we identified 27 intersection genes with HB treatment against RAS. PPI network analysis demonstrated that one active ingredient matches multiple targets, and one target matches multiple ingredients; reflecting the complex mechanism of multi-ingredient and multi-target against diseases in Chinese medicine. Gene pathway network analysis revealed that TNF, IL6, IL1B, PTGS2, CCL2, and IL10 are the core targets of the HB for RAS treatment. Previous studies suggested that TNF-α plays a vital role in the pathological process of RAS. Oliveira 45 et al. found in their experiment that reducing TNF-α expression optimized the healing of traumatic oral ulcers in diabetic rats. TNF-α, as a pro-inflammatory factor, is involved in the development of RAS and can be used as an indicator of the severity of the disease. PTGS2 encoding cyclooxygenase-2(COX-2) was associated with the development of RAS. Experiments in rats showed a decrease in PTGS2 activity and a decrease in COX-2 protein expression during the healing of oral ulcers. 46
GO enrichment analysis were mainly on biological processes, including response to lipopolysaccharide, response to molecule of bacterial origin, cellular response to lipopolysaccharide, and cellular response to molecule of bacterial origin. Cellular response to molecule of bacterial origin and reactive oxygen species metabolic process are important BPs involved in the development of RAS. Wound healing and aging are key contributors that cause RAS. Reactive oxygen species metabolic process may trigger the stress response by regulating gene expression to shift redox homeostasis to a more oxidative state. This process facilitates the onset of RAS. Positive regulation of protein phosphorylation is a central signaling pathway that regulates cell survival, playing significant roles in the inhibition of RAS47,48.
The mechanisms of HB for treating RAS are fulfilled at least two function modules, including the inhibition of inflammation and the regulation of immunological response. The implementation of the functions relies on a smooth run of the complex multi-pathways network, particularly TNF signaling pathway and HIF-1 signaling pathway. CCL2, CXCL10, IL6, MMP3, TNF, IRF1, SELE, VCAM1, IL1B, PTGS2, ICAM1, and FOS are involved in TNF-1 signaling pathway. Previous studies have been reported that excessive production of TNF-α, IL-6, and IL-1β increased the risk of RAS development.49–51 In addition, matrix metalloproteinases play a negative role in mucosal pathology of the oral cavity.52,53 In summary, the anti-inflammatory effect is an important factor for HB contributing to the treatment of oral ulcers. Three targets are involved in HIF-1 signaling pathway which is an important regulator of the immune system, including IL6, NOS2, and IFNG. NOS2 is the isoforms of NOS, which produce high-level sustained NO synthesis and strongly affect adaptive immune responses. 54 NO regulates the responses of many immune and inflammatory cell types like macrophages which are associated with some of the most important immune pathologies. 55 Increased serum NO level is thought to take part in the pathogenesis of RAS (25). Thus, NOS2 plays an important role in the HIF-1 signaling pathway for the regulation of anti-inflammation and immune protection functions of HB. Immunologic disturbances play a crucial role in the etiopathogenesis of oral ulcers. Hence, it is worthwhile to explore the impact of HB on immune function. According to the pathway enriched results, many target genes are enriched in the pathways related to immune inflammatory diseases. For example, nine targets are enriched in Toll-like receptor signaling pathway (p-value =1.45 × 10−10). Six targets are enriched in NF-kappa B signaling pathway (p-value =4.94× 10−6) and four targets in T cell receptor signaling pathway (p-value =1.16 × 10−6). Therefore, it can be presumed that one pharmacological function of HB is related to the regulation of immune response for RAS.
Molecular docking displayed components have good binding to core target genes. All selected core proteins had good affinity for the ligand (<−5 kcal/mol) and all effective docking fractions were > −7 kcal/mol, indicating that the compound has a strong binding affinity for the docked protein. 56 In silico validation of key targets for bioactive components offers an alternative path for the exploration of ligand-target interactions and action mechanisms. We could conclude that the different components display diverse interactions with these targets in the predicted pathways. The results demonstrated a synergistic effect mode of the HB in its oral ulcer treatment effects. PTGS2, which has the highest docking score, is a biomarker of iron death and inhibits the expression of inflammatory factors and apoptosis.57,58 Future research directions can start from iron death mechanism pathway.
The limitation of this study is the lack of clinical or animal experiments. Further studies will combine molecular biology and pathophysiology to validate the predicted potential key targets and pathways, and explore the mechanism of action of the active ingredients.
Conclusion
In this study, we used network pharmacology and molecular docking methods to investigate the mechanism of HB against RAS. HB exhibited the therapeutic effects on RAS probably by inhibiting inflammation and regulating immunological response. TNF signaling pathway associated with inflammation and HIF-1 signaling pathway associated with immune system may play crucial roles in the protection of HB against recurrent aphthous stomatitis. This study is expected to be useful for the research and development of the drug and further elucidates its mechanism.
Supplemental Material
Supplemental Material - Network pharmacology and bioinformatics analyses identify the intersection genes and mechanism of Huang Bai for recurrent aphthous stomatitis
Supplemental Material for Network pharmacology and bioinformatics analyses identify the intersection genes and mechanism of Huang Bai for recurrent aphthous stomatitis by Lulu Tang, ling Huang and yingtao Lai in International Journal of Immunopathology and Pharmacology
Supplemental Material
Supplemental Material - Network pharmacology and bioinformatics analyses identify the intersection genes and mechanism of Huang Bai for recurrent aphthous stomatitis
Supplemental Material for Network pharmacology and bioinformatics analyses identify the intersection genes and mechanism of Huang Bai for recurrent aphthous stomatitis by Lulu Tang, ling Huang and yingtao Lai in International Journal of Immunopathology and Pharmacology
Footnotes
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.
Ethical approval
The data were obtained from public databases, so the ethics committee approval was unnecessary.
Availability of data and materials
Supplemental Material
Supplementary material for this article is available on the online.
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
