Mechanisms of Spica Prunellae against thyroid-associated Ophthalmopathy based on network pharmacology and molecular docking

Background Thyroid-associated ophthalmopathy (TAO) is an autoimmune inflammatory disorder, which lacks effective treatment currently. Spica Prunellae (SP) is popularly used for its anti-inflammatory and immune-regulating properties, indicating SP may have potential therapeutic value in TAO. Therefore, the purpose of this study is to identify the efficiency and potential mechanism of SP in treating TAO. Methods A network pharmacology integrated molecular docking strategy was used to predict the underlying molecular mechanism of treating TAO. Firstly, the active compounds of SP were obtained from TCMSP database and literature research. Then we collected the putative targets of SP and TAO based on multi-sources databases to generate networks. Network topology analysis, GO and KEGG pathway enrichment analysis were performed to screen the key targets and mechanism. Furthermore, molecular docking simulation provided an assessment tool for verifying drug and target binding. Results Our results showed that 8 targets (PTGS2, MAPK3, AKT1, TNF, MAPK1, CASP3, IL6, MMP9) were recognized as key therapeutic targets with excellent binding affinity after network analysis and molecular docking-based virtual screening. The results of enrichment analysis suggested that the underlying mechanism was mainly focused on the biological processes and pathways associated with immune inflammation, proliferation, and apoptosis. Notably, the key pathway was considered as the PI3K-AKT signaling pathway. Conclusion In summary, the present study elucidates that SP may suppress inflammation and proliferation and promote apoptosis through the PI3K-AKT pathway, which makes SP a potential treatment against TAO. And this study offers new reference points for future experimental research and provides a scientific basis for more widespread clinical application.


Background
Thyroid associated ophthalmopathy (TAO), an inflammatory disorder affecting the orbit and ocular adnexa, is occurred in autoimmune thyroid disease, mainly Graves' disease [1]. A prospective study determined that the overall incidence rate of TAO was 10 per 10,000 persons in Europe [2]. And it is reported to occur in 20-30% of patients with Graves' disease and 2% of patients with Hashimoto's thyroiditis [3][4][5]. The pathogenesis of TAO is related to many factors, including genetic factors, environmental factors and autoimmunity, in which autoimmune disorder plays an essential role. An intricate interplay between the potential pathogenic autoantigens and the autoantibodies leads to the activation of autoimmune response and releasing related inflammatory factors, resulting in inflammation of the orbital tissues, abnormal fibroblast proliferation, and extraocular muscle thickening [5][6][7]. Clinically, it is characterized by proptosis, upper eyelid retraction, restrictive strabismus, exposure keratopathy and other presentation. As a result, its impact on vision and appearance in daily living is significantly associated with the quality of life and psychosocial problem [8][9][10][11]. However, current treatments are inadequate due to the incompletely understood pathophysiology of TAO. For mild disease, most patients experience spontaneous recovery within about 6 months [12]. Although helping with symptom management, 13.5% of patients had a progressively clinical deterioration [13]. For active severe phase, corticosteroids remain the mainstay treatment and orbital radiotherapy commonly serves as an adjuvant therapy [14][15][16]. However, these therapeutic approaches are often accompanied with limited efficacy and safety concerns, such as hormone resistance and complications [5]. Therefore, there is an urgent demand for novel effective therapy to improve the progression of TAO.
Spica Prunellae (SP), the fruit-spikes of the Prunella vulgaris L., is a traditional antipyretic Chinese herb with a wide distribution in Northeast Asia. According to the theory of traditional Chinese medicine (TCM), SP has long been thought to have benefits involving clearing the fire and eyesight, relieving edema and removing stasis. Moreover, SP has been extensively used in food additives and pharmaceuticals [17]. Previous pharmacologic studies indicated that the SP has multiple biological functions such as anti-inflammatory [18], anti-cancer [19,20] and immunomodulatory activities [21]. Nowadays, SP has a wide range of applications in thyroid disease, for instance, thyroid gland swell and subacute thyroiditis [22,23]. Indeed, some Chinese herbal formulations applied to TAO therapy employ SP as an important ingredient [24]. However, the underlying mechanism is yet to be completely elucidated.
Based on the development of bioinformatics, network biology and pharmacology analysis, network pharmacology becomes a novel strategy to understand the law of interactions between multicomponent and multitarget systematically and comprehensively [25,26]. In the latest years, network pharmacology has been successfully used to predict the mechanism of TCM in a variety of diseases treatment. Therefore, we adopted the network pharmacology approach to investigate the potential therapeutic targets and pathways of SP in TAO treating, and performed molecular docking studies to further predict the recognition and interaction modes between SP and its predicted targets.

Screening of active compounds of SP
To collect the pharmacologically active compounds of SP, we searched Traditional Chinese Medicine System Pharmacology Database [27] (TCMSP, http://lsp.nwu. edu.cn/tcmsp.php), a systematic pharmacology platform of Chinese herbal medicine. Based on their ADME (absorption, distribution, metabolism and excretion) features, Oral bioavailability (OB) and Drug-likeness (DL) were set as the key parameters [28]. OB refers to the percentage of which the pharmaceutical agents is used and unchanged throughout systemic circulation [29]. And DL is a concept of optimize pharmacokinetics and properties [30]. Then, based on the threshold values of OB 30% and DL 0.18 [31], eleven active compounds were selected. In addition, five researches [32][33][34][35][36] using high-performance liquid chromatography (HPLC) to identify active compounds of SP were screened through literature search. After integrating with the results obtained from TCMSP, total twenty active compounds of SP were used to the subsequent network pharmacology study.
Analysis of overlapping targets of SP in TAO Target mapping To obtain the candidate targets responsible for TAO therapy, the putative targets of SP mapped to the TAO related targets to obtain overlapping targets. For overlapping targets, numbers were visualized with a Venn diagram generated in R version 3.2.0, a freely available, open-source statistical programming language and environment for statistical computing [45], with VennDiagram package.

Protein-protein interaction (PPI) analysis
These related targets of SP in treating TAO were then put in STRING tools [46] (https://string-db.org/), an online platform that critically predict functional protein association networks, to generate protein interaction network. Network data were used for the topological properties analysis to identify important targets in anti-TAO system, and Table S4 for more information.

Gene ontology (GO) and KEGG pathway enrichment analysis
To further reveal the potential mechanism of SP in TAO treating, we imported overlapping targets into the Functional Annotation tool of Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 [47] (https://david.ncifcrf.gov/) to perform GO functional analysis and KEGG pathways analysis. The enrichment P values of functional annotations were corrected by both Bonferroni (P < 0.05) and Benjamini (P < 0.05) correction. The plots of GO enrichment and KEGG pathway enrichment were then carried out with R version 3.2.0.

Network construction
The Cytoscape3.6.1 [48] (http://cytoscape.org/), a software environment for integrated models of biomolecular interaction networks, was utilized to construct the following networks:(1) Compound-target network; (2) PPI network of TAO targets; (3) PPI network of compound-TAO targets; (4) Compound-target-pathway network. The Cytoscape plug-in Network Analyzer was applied to analyze the topological properties of the PPI network. Three indices were used to assess the topological properties of every node. (1) "Degree" refers to the number of nodes that directly interact with that node in a network [49], (2) "Betweenness Centrality": the proportion of the number of nodes that have passed through the shortest path in a network [50], with important influence due to their control over information passing between other nosdes. Thus, the betweenness reflects its global importance, while the degree measures the local importance of a node; (3) "Closeness Centrality" reflects the closeness of a node to other nodes [51]. The level of these three paraments reflects the topological importance of the node in a network, which the higher the values are, the more important the nodes are. The targets with degree > twofold median were selected as the hub genes and used for further molecular docking.

Molecular docking simulation Target protein preparation
The crystal structures of hub genes were downloaded from RCSB Protein Data Bank [52] (http://www.pdb.org/ ). The downloaded complexes were embellished by PyMol2.3.0 (the PyMol Molecular Graphics System) to remove original ligand, water molecules and phosphates. Moreover, the AutoDock Tools1.5.6 [53] (http:// mgltools.scripps.edu/documentation/links/autodock) was used to prepare receptors, including adding hydrogen and seting docking paraments. "Grid box" was set to maximum to perform the blind docking.

Ligand preparation
Before docking, the 2D structures of compounds were drawn using the ChemBioDraw13.0 software. Then the ChemBio3D13.0 software was used to transfer the 2D structure to 3D chemical structure and make energy minimizing for the further docking.

Molecular docking
All ligand and receptor files were saved as pdbqt format. Then we used Autodock Vina [54], a freely available open-source packages, to evaluate and verify the binding affinity of compound-target relationship, and the prediction results from network pharmacology. The binding models were visualized by PyMol2.3.0 software and Discovery Studio3.5 software.

Compound-target network
A flowchart of the network pharmacology study is summarized in Fig. 1. Based on TCMSP, 11 active compounds were selected by the screen criteria of OB and DL. Considering the large number of researches performed HPLC to identify the major active components of SP, we conducted the literature search in Pubmed and China National Knowledge Infrastructure (CNKI). And nine compounds were supplemented accounting for the relatively large proportion and proven active role in SP, including Betulinic acid, Caffeic acid, Corosolic acid, Euscaphic acid, Maslinic acid, Rosmarinic acid, Rutin, Oleanolic acid and Ursolic acid (Table 1). Then we identified 610 putative targets from TCMSP, Pubchem, STITCH and SwissTargetPrediction, which the detailed information is described in Table S1. To obtain further understanding of the relationship between the 20 compounds and their corresponding targets on a system level, we conducted the compound-target network (Fig. 2a). In a network topology analysis, we identified top six compounds adopting the average degree (103.94) as the threshold value: Quercetin (degree = 298), Ursolic acid (degree = 155), Luteolin (degree = 149), Kaempferol (degree = 146), Morin (degree = 135) and Delphinidin (degree = 113).

PPI network of compound-TAO targets
Based on the above results, the 610 putative targets of SP mapped to the 813 TAO related targets to obtain overlapping targets. Consequently, 127 targets were identified as the candidate targets responsible for TAO therapy (Fig. 3a). A PPI network was then constructed to evaluate the role of the targets in complex disease and discover the interactive effects. After analyzing the topological feature of the abovementioned PPI network, 127 genes were sorted in descending order by degree in Table S4 with an average degree of 36.74. As shown in Fig. 3b, nodes arranged in a concentric circle according to the degree and the innermost center circle is composed of 18 hub nodes, including AKT1, ALB, CASP3 (caspase-3), CXCL8 (interleukin-8), EGF (epidermal growth factor), IL10 (interleukin-10), IL1B (interleukin-1 beta), IL6, INS, JUN (transcription factor AP-1), MAPK1 (mitogen-activated protein kinase 1), MAPK3 (mitogen-activated protein kinase 3), MMP9 (matrix metallopeptidase 9), PTGS2 (prostaglandin-endoperoxide synthase 2), STAT3 (signal transducer and activator of transcription 3), TP53 and VEGFA (vascular endothelial growth factor A). These hub genes are of significance in the network and used for the following molecular docking study.

GO and KEGG pathway enrichment analysis
To further elucidate the mechanism of drug treatment systematically, the enrichment analysis was performed with DAVID on 127 genes. The top 10 GO items and KEGG pathways are selected based on counts of hit genes and P values (Fig. 4). For biological processes, it can be found that the targets were mainly  For pathway analysis, most involved in immune and inflammatory related pathways such as PI3K-AKT signaling pathway (hsa04151), NOD-like receptor signaling pathway (hsa04621), Toll-like receptor signaling pathway (hsa04620), NF-kappa B signaling pathway (hsa04064) and TNF signaling pathway (hsa04668) ( Table 2). These results demonstrated that the main action mechanism underlying TAO treatment.

Molecular docking analysis
In the present studies, the possible interaction activity between 18 hub genes and their corresponding  Fig. 5A (a)). It showed that six hydrogen bond formation between ligand and residues invloved in LEU 188, ALA 189, GLN 227, TYR 248 and MET 247. Three hydrogen bond interactions were considered as strong interaction with the distance of 3.2, 3.3, 3.5 Å respectively. Moreover, HIS 226 residue contributed to form two pi-pi stacking interactions with two benzene rings of ligand. The other essential residues (GLY 186, LEU 187, VAL 223, ARG 249, LEU 222, LEU 243, ALA 242, TYR 245 and PRO 246) interacted with Luteolin through electrostatic forces, van der Waals forces, etc. (Fig. 5A(b)). Consequently, Luteolin stably binded to MMP9 protein bonds through mutiple interaction forces. Overall, we found that hydrogen bond and electrostatic forces were the main forms of interactions of the top 10 docking complexes, and the diverse of the interaction forms determined the ability of binding affinity.

Compound-target -pathway network
As shown in Fig. 6A, a compound-hub gene-pathway network was established by connecting potential pathways, compounds and hub genes. To find the major nodes, we took the average of respective degree as the threshold values to identify the significant compounds and targets. Three compounds and 8 targets were selected as core nodes in anti-TAO system: Quercetin (degree = 14), Ursolic acid (degree = 11) and Rutin (degree = 8), PTGS2 (degree = 17), MAPK3 (degree = 14), AKT1 (degree = 13), TNF (degree = 13), MAPK1 (degree = 10), CASP3 (degree = 9), IL6 (degree = 9) and MMP9 (degree = 9). Subsequently, for a more detailed explanation of the associations described above, we extracted a subnetwork to exhibit the core nodes and key mapping pathways (Fig. 6b).

Discussion
In recent years, there has been growing researches devoted to drug discovery and combined treatment with TCM against complex disease such as TAO, in which the network pharmacology approach plays an essential role. In the present study, we constructed the following network to reveal the potential targets and pathways of SP in TAO treatment: compound-target network, PPI network of TAO targets, PPI network of compound-TAO targets and compound-target -pathway network. After integrating and consolidating information from diverse sources of available databases, 20 active compounds of SP acted on 127 different targets associated with TAO. According to the compound-target network, six compounds and putative targets are highly connected and can be defined as vital compounds in SP. Previous studies have indicated the anti-inflammatory, proapoptosis and immunomodulation effects of Quercetin [55,56], Ursolic acid [57][58][59], Luteolin [60][61][62], Kaempferol [63][64][65], Morin [66,67] and Delphinidin [68,69] in various diseases, suggesting SP may exhibit beneficial effects in TAO which is characterized by immunoinflammatory response. Actually, the PPI network of TAO  targets did reflect the underlying pathogenesis of TAO, which are consistent with prior researches. Commonly reported mechanism [6] is the promotion and mutual interaction of these including that cellular immunity and humoral immunity mainly caused by massive activation of T lymphocytes, inflammatory response and proliferation of orbital fibroblasts. This manifests as an uncontrolled secretion of inflammatory cytokines and growth factors by activated orbital fibroblasts, leading to orbital tissue remodeling and enlargement. Further, to better elucidate the mechanism of SP in treating TAO, the 127 targets were screened for 18 hub genes in the PPI network of compound-TAO targets. Coincidentally, there are 6 targets considered as key pathogenic gene in the PPI network of TAO targets, which are involved in the 18 hub genes, namely TNF, ALB, IL6, INS, AKT1 and TP53, reconfirming SP might possess good effects on TAO.
According to the compound-hub gene-pathway network, Quercetin, Ursolic acid and Rutin interacted with the large number of targets, indicating the important roles in the anti-TAO system. Quercetin is a flavonoid phytoestrogen exhibiting antioxidant and anti-inflammatory properties and reducing proliferation in orbital fibroblasts [70][71][72]. In addition to the above effects, Ursolic acid [57][58][59] and Rutin [73,74] are reported to promote apoptosis and regulate immune in other cell systems and animal models. Eight important therapeutic targets with the degree > 7.83 are PTGS2, MAPK3, AKT1, TNF, MAPK1, CASP3, IL6 and MMP9. Meanwhile, we performed molecular docking simulation between 18 hub genes and 20 active compounds to complement key targets. The results suggested that all of the target-compound pairs possess good docking affinity, especially, MMP9, PTGS2, CASP3 and ALB. ALB is the most abundant protein in human plasma, to which drugs usually bind, playing a crucial importance to understand the pharmacodynamics and pharmacokinetics [75]. Beyond that, the drop in ALB level reflects the oxidative stress damage [76], albeit not confirmed in TAO pathogenesis. We therefore speculate that the good affinity of ALB and Betulinic acid may refer to the general drug response rather than the treatment mechanism of SP. Consequently, SP plays a beneficial role in TAO, at least in part, through modulating PTGS2, MAPK3, AKT1, TNF, MAPK1, CASP3, IL6 and MMP9. GO functional annotation and pathway enrichment analysis further supported the potential therapeutic mechanism, the results of which was predominated by immune inflammation, proliferation and apoptosis.
In terms of immune inflammation, proinflammatory cytokine PTGS2, IL6 and TNF were validated the participation in TAO onset. PTGS2 is involved in prostaglandin biosynthesis, which have been ascribed key role in inflammation [77]. Prior studies have noted that PTGS2 diminished with a decrease in clinical activity score in TAO [78,79], is currently believed to be pivotal to inflammatory process in TAO patients. IL6 has been implicated in the pathogenesis of autoimmune disease [80]. It has been reported that AKT1/NF-κB signaling pathway contributed to the production of IL6 in the retrobulbar space in active phase of TAO [81]. Likewise, research showed that the elevation levels of TNF in serum is mediated by AKT1/NF-κB signaling pathway in inflammatory response of TAO [82]. At present, there has been reported some inhibitors of TNF obtained promising results in patients with TAO, irrespective of rare adverse effects [83]. Then SP with TNF as a key therapeutic target may have potential to be an effective agent for TAO. Hence, it is possible that SP could suppress the inflammatory cytokines secretion of T lymphocytes and orbital fibroblasts through AKT1/NF-κB signaling pathway.  Regarding whether SP could modulate immune like corticosteroids exerts [83], for instance, decrease antigen presentation to T lymphocytes in TAO, it requires further laboratory investigation. For proliferation and apoptosis, PI3K-AKT signaling pathway seems to be significant in mediating cell growth and death in TAO. Recently studies demonstrated that orbital fibrocytes express thyrotropin receptor (TSHR) [84] and activation of TSHR leads to increased inflammatory gene expression and proliferation through PI3K-AKT pathway [85][86][87]. Besides this, the increased cell proliferation of orbital preadipocytes is also contributed in TAO progression [88]. CASP3 activation is one of the final steps in the execution of apoptosis, which is activated in orbital preadipocytes by some compounds in TAO [88]. And SP exhibits pro-apoptotic actions via activation of CASP3 in various cancer [89]. Although there is a lack of relevant studies of MAPK1, MAPK3 in TAO, they act in a signaling cascade that regulates various cellular processes such as proliferation and apoptosis [90]. MMP9 as a kind of matrix metalloproteinases has been proved to participate in inflammation of multiple tumors through TNF signaling pathway [91]. However, it still requires further researches to elucidate the role of MMP9 in TAO pathogenesis. Therefore, it can be summarized from the results that whether in immune inflammation or proliferation and apoptosis, the PI3K-AKT signaling pathway plays a key role in TAO, and inhibition of this process could be an effective therapeutic target for SP against TAO (Fig. 7).

Conclusion
In conclusion, up till now, TAO treatments have been an intractable problem, but TCM may be a potential alternative treatment to a certain extent. In the present study, we predicted 8 key targets from complex networks through the integration of network pharmacology and molecular docking, and performed comprehensive explanation of the therapeutic mechanism of SP in TAO, which is likely to focus on three aspects: immunoinflammatory suppression and regulation of proliferation and apoptosis. In addition, we found that PI3K-AKT signaling pathway may occupy core status in anti-TAO system. Our research would be enlightening for developing therapeutic methods targeted TAO. However, further extensive experiments will be necessary to unravel the above possibility, such as in vivo and in vitro efficacy experiments of SP in TAO therapy.
Additional file 1: Table S1 The putative targets of SP. Fig. 7 The PI3K-AKT signaling pathway plays a central role in anti-TAO system of SP. The red nodes represent key targets, the yellow nodes represent overlapping targets of SP and TAO targets, and the green nodes represent the other targets in PI3K-AKT signaling pathway. Grey ellipses represent important pathways associated with TAO and mediated by PI3K-AKT pathway as well