Network Pharmacology-Based Prediction and Verification of the Targets and Mechanism for Panax Notoginseng Saponins against Coronary Heart Disease

Coronary heart disease (CHD) is the worldwide leading cause for cardiovascular death. Panax notoginseng saponin (PNS), which is the main bioactive compound of panax notoginseng, has been generally accepted to exert a remarkable effect on CHD for a long time. However, to reveal the underlying treatment target and corresponding mechanism of PNS against CHD is still a substantial challenge. In this work, the targets and mechanism of PNS against CHD were successfully achieved by pharmacology-based prediction and experimental verification. 36 common targets were screened out through integrating the gene expression profile of CHD and the chemical-protein data of PNS. Then, two key nodes were further selected for verification by experiment after analyzing GO function, KEGG pathway, coexpression, and topology analysis. Results showed that PNS has protected the human umbilical vein endothelial cells from H2O2-induced oxidative stress by inhibiting early cell apoptosis via upregulating VEGFA mRNA expression. Therefore, our research has successfully pointed out one treatment target and apoptotic inhibition caused by PNS with method of integrating bioinformatics prediction and experimental verification, which has partially explained the pharmacological mechanism of PNS against CHD.


Introduction
Cardiovascular disease (CVD), as the leading reason in death worldwide, would increasingly result in over 23.6 million deaths by 2030 [1]. Coronary heart disease (CHD) is the most common cause of cardiovascular death, which accounts for about 44% [2]. 2018 AHA statistical data indicated that 16.5 million people over the age of 20 had CHD, and myocardial infarction (MI) occurred per 40 seconds approximately in America [2]. Moreover, silent MI or undiagnosed MI may account for almost 50% of the whole MIs [3]. This nonmanifestation MI is also an important risk factor for development of future heart failure [4]. Therefore, it could be speculated that the actual incidence of CHD is possibly higher than the reported data. Currently, the treatment strategy of CHD has achieved much progress involving new oral anticoagulant/antiplatelet drugs and drug eluting stent. However, recent treatment methods such as antiplatelet therapy, lipid regulation, and percutaneous coronary intervention (PCI) are still unsatisfactory for lots of patients worldwide due to difficulties in overcoming adverse effects. For example, as one cornerstone of CHD treatment, antiplatelet therapy encountered challenges of drug resistance [5,6] and inappropriate use [7,8] in clinical practice. Statins for lipid regulation sometimes induced the damage of liver and muscular cells [9]. Regarding PCI, there was increasing death rate observed from 2004 to 2014 [2]. As a result, it is noted that more effective and safer treatments for CHD are still urgently needed and continue to be challenging.
Panax notoginseng, as one of the traditional Chinese medicines (TCM) with effects of promoting blood circulation and stopping bleeding, has attracted enormous researches in the treatment of CVD. Panax notoginseng saponin (PNS) is the main bioactive compound in panax notoginseng, exerting cardiovascular protection. PNS mainly contains five compounds, namely, notoginsenoside R1, ginsenoside Rg1 (also called sanchinoside C1), ginsenoside Re, ginsenoside Rb1, and ginsenoside Rd (Figure 1), according to Pharmacopoeia of the People's Republic of China [10] and Hong Kong Chinese Materia Medica Standards [11]. Currently, PNS has been made into several preparations for CHD therapy. A meta-analysis revealed that panax notoginseng preparations are beneficial for unstable angina pectoris patients [12], possibly owing to the vital importance of PNS in antiplatelet aggregation [13], anti-inflammation [14], antioxidation [15], antiatherosclerotic [16], antiapoptosis [17], and lipid regulation [16]. However, essential compounds and specific targets of PNS are not clearly determined and the underlying mechanism against CHD remains poorly defined. Therefore, to reveal and understand the treatment targets and mechanism of PNS against CHD is a subject which is worthy of study.
Since network pharmacology is an effective approach for drug target prediction, it has been increasingly growing in popularity [18][19][20]. It could offer complementary support for new drug development and better understanding of complex mechanisms for drugs. Moreover, it can also provide systematic and holistic methods to deeply reveal the druggene-disease interactions which are especially suitable for understanding the multitarget, complementary, and synergic essence of TCM in molecular networks. Nevertheless, analysis of drug target based on network pharmacology is still in its infancy stage [21]. Predicted targets still need further confirmation from experiments. In this study, we predicted and verified the targets and mechanism for PNS against CHD through coupling the network pharmacology approach and experimental method. By integrating bioinformation from databases involving TCMSP, PharmMapper, STITCH, Gene Expression Omnibus (GEO) DataSets, WebGestalt, and GeneMANIA, key targets and corresponding mechanism of PNS against CHD were predicted. Then, cell experiment was applied to verify two key genes. As a result, one verified gene was identified as the treatment target of PNS against CHD, and the potential mechanism was regulation of apoptosis. The whole flow chart of our network prediction and experimental verification was shown in Figure 2. Topology analysis Figure 2: The whole framework of targets prediction and verification. ADME properties of PNS were first evaluated and then potential treatment targets of PNS against CHD were predicted through identifying common genes. Next, GO function, KEGG pathway, gene coexpression, and topology analysis were performed to predict key targets and their gene function and pathways. Finally, cell experiment was conducted to verify the selected targets and underlying mechanism.

Evaluation of ADME-Related Properties for PNS.
PNS mainly contains notoginsenoside R1, ginsenoside Rg1, ginsenoside Re, ginsenoside Rb1, and ginsenoside Rd. The last compound is often exhibited as the metabolite, ginsenoside Rd qt. In order to assess druggable characters of the above five compounds, ADME (abbreviation of absorption, distribution, metabolism, and excretion) properties are evaluated by using TCMSP server (http://lsp.nwu.edu.cn/tcmsp.php).
Particularly, important parameters of oral bioavailability (OB), drug-likeness (DL), molecular weight (MW), bloodbrain barrier (BBB), intestinal epithelial permeability (Caco-2), fractional negative surface area (FASA-), the polar surface area (TPSA), and RBN [22] are used. Among all the descriptive parameters, OB and DL are the most two important indicators [23,24] to assess and identify candidate compounds or drugs. OB is closely related to drug absorption and delivery of an orally administered dose. Drugs with high OB value mean they are more likely to reach the potential targets and can be therapeutic agents. DL can be used to estimate the druggability of a drug qualitatively, which helps to optimize pharmacokinetic and pharmaceutical properties. Commonly, OB⩾30% and DL⩾0.18 are suggested as drug screening criteria in the traditional Chinese herbs [22]. In the TCMSP server, OB value is calculated by an in-house model OBioavail1.1 [22]. Besides, the Tanimoto coefficient [25] is used to assess DL indexes of each compound in PNS. The formula is shown as follows [26]: Evidence-Based Complementary and Alternative Medicine

Targets Prediction of PNS against CHD.
Pharmmapper (http://lilab-ecust.cn/pharmmapper/) and STITCH (http:// stitch.embl.de/) were used to identify targets related to PNS in the study. Pharmmapper is a reverse docking server which identifies potential protein targets for small molecule compounds via a pharmacophore mapping approach. By calculating the fit scores and z'-scores (normalized fit scores) of a compound with different targets respectively, the correlation degree between the query compound and its targets are predicted. Generally, large positive z'-score indicates high significance [27], while STITCH provides a collection of protein-chemical interactions that integrates many sources of experimental and manually curated lines of evidence with text-mining information and interaction predictions [28,29]. It calculates a score for each pair of proteinchemical interaction. Usually, interactions with score ⩾0.4 are considered as medium confidence and those with the score⩾0.7 are high confidence [30].
In addition, gene expression series of CHD were downloaded from GEO DataSets [31] (http://www.ncbi.nlm .nih.gov/geo/). Four researches with GSE20686, GSE71226, GSE18612, and GSE120774 accession numbers were offered in the platform of GPL4133, GPL570, and GPL6244, respectively. The expression profiles were all received by array. After unifying different gene names as Gene Symbol by UniProtKB (http://www.uniprot.org), the differential expression genes (DEGs) of CHD were obtained. Furthermore, Venn diagram (http://bioinformatics.psb.ugent .be/webtools/Venn/) was drawn to identify the intersection nodes of PNS and CHD. Finally, 36 common targets were predicted.

GO Function, KEGG Pathway, and Key Targets Analysis.
GO function and KEGG pathway of 36 common targets were analyzed with WebGestalt (http://www.webgestalt.org/). Then, coexpression network of common targets was established by GeneMANIA (http://pages.genemania.org/) which was usually applied to generate hypotheses about gene function, analyze gene lists, and prioritize genes for functional assays [32]. In addition, topology analysis was conducted by CentiScape and MCODE to calculate the degree centrality and extract subnetwork, in order to identify key nodes in the coexpression network. Ultimately, two key targets which participate in apoptotic signaling pathway were selected to perform further cell experiment.

Experimental Materials and Cell
Culture. Human umbilical vein endothelial cells (HUVECs) were purchased from Bei Na Chuanglian Biotechnology Research Institute, Beijing, China. PNS was purchased from National Institutes for Food and Drug Control, Beijing, China. PNS was dissolved with Phosphate Buffered Saline (1X) (Hyclone, Lot AC12557265) before experiment. Dulbecco's modified Eagle's medium (DMEM, Hyclone, Lot AD19822332) supplemented with 10% fetal bovine serum (Gibco, Lot 1739463) and 1% penicillinstreptomycin (Hyclone, Lot J160016) was used to culture HUVECs. The cultures were maintained in a humidified atmosphere containing 5% CO 2 at 37 ∘ C. Then cells were inoculated onto 6-well plates to perform experiments and each group had 3 samples. Treatment groups were pretreated with 30 g/mL PNS for 21 h and then added into 100 mol⋅L −1 H 2 O 2 for 3 h. In contrast, model groups were only treated with 100 mol⋅L −1 H 2 O 2 in the last 3 h. Cells in control groups were cultured in medium added with equal volume of DMEM.

Cell Apoptosis Analysis.
HUVECs were seeded in 6well plates (1×10 6 cells/well) and incubated with 30 g/mL PNS and/or 100 mol⋅L −1 H 2 O 2 . After the culture for 24 h, cells were harvested and washed with cold PBS twice. Then cells were resuspended in 1×Binding Buffer and stained with FITC Annexin V and Propidium Iodide (PI) supplied in the FITC Annexin V Apoptosis Detection kit I (BD Pharmingen, USA). Next, the cells were incubated for 15 min at RT (25 ∘ C) in the dark. Finally, 1×Binding Buffer was added to cells and those cells were analyzed by a flow cytometer (BD FACSCalibur, USA) within 1 h. All experiments were performed in triplicate.

Statistical Analysis.
The quantitative data were presented as mean±standard deviation (±SD). One-way analysis of variance (ANOVA) was applied to compare multiple means. When the variance was equal, the least-significant difference (LSD) method was used. Dunnett's T3 test was conducted when the variance was not equal. Significance was defined as P<0.05. All statistical analysis was performed with SPSS 17.0 and data were presented as two decimals. In addition, all statistical column charts were drawn with Origin 8.0.  DL value above 0.18, which indicated that they presented excellent druggability such as better solubility and chemical stability [26]. Similarly, their OB indices which were estimated to be 10.04% and 12.23%, respectively, demonstrated that they can be efficiently absorbed to systemic circulation, compared with other three compounds. As for MW, TPSA, and RBN parameters, ginsenoside Rd qt showed better properties than others as well (see Table 1). Therefore, it could be seen that ginsenoside Rd qt was more druggable and may be the main bioactive compound in oral management.

Common Targets of PNS and CHD. CHD-related genes
were searched in GEO database and the series (GSE20686, GSE71226, GSE18612, and GSE120774) were selected. The former two contained DEGs of 399 samples totally from peripheral blood, while the latter two involved DEGs of 49 samples totally from epicardial and subcutaneous adipose tissue which were closely related to coronary atherosclerosis   [ [33][34][35]. CHD is diagnosed when the luminal stenosis⩾50% in at least one major vessel. Therefore, patients who met the above diagnostic criteria were assigned to the CHD group, and the rest belonged to non-CHD group, when analyzing GSE20686 and GSE71226 data. Regarding GSE18612 and GSE120774, patients who underwent cardiac surgery for coronary artery disease were assigned to the CHD group; and those undergoing cardiac surgery for other reasons were considered as non-CHD group. By comparing CHD and non-CHD groups with GEO2R, there were 1775 DEGs finally identified in the above four series with selection criteria of FC⩾1.5 and P<0.05. Moreover, further comparison between 1775 CHD-related genes and 679 PNS targets identified 36 common genes that may become the predicted treatment targets of PNS against CHD (Figure 4).

GO Function and KEGG Pathway Analysis of 36 Common
Targets. To investigate the 36 common genes of PNS and CHD, GO function and KEGG pathway were analyzed with WebGestalt. GO enrichment was analyzed from the perspective of biological process, cellular component, and molecular function. It was noted that the common genes were mainly responsible for protein binding, biological regulation, response to stimulus, and metabolic process ( Figure 5). Importantly, two apoptotic signaling pathways have been enriched in the top 10 GO functions (as shown in Table 2). Additionally, KEGG pathways of the common genes were analyzed and shown in Figure 6. Among all the pathways, hsa04926, hsa05219, hsa04915, hsa05418, hsa05205, hsa04370, hsa04211, hsa04657, hsa01522, and hsa04933 enriched the most genes and ranked the top 10 levels (Table 3). Particularly, relaxin signaling pathway, fluid shear stress and atherosclerosis, VEGF signaling pathway, and IL-17 signaling pathway were reportedly related to CHD [36][37][38].   pathways, and colocalization, those 34 nodes established a coexpression network with other 20 nodes in GeneMANIA (Figure 7). In addition, topology analysis was performed to identify key nodes in the coexpression network, involving calculation of degree centrality and k-core decomposition by CentiScape and MCODE, respectively, in Cytoscape 3.6.1. Degree centrality refers to the degree to which a node in a network is associated with all other nodes. Node with higher degree usually means more interaction edges and more significance in the network. The k-core of a network is defined as the largest subnetwork where every node has at least k links. Commonly, the extracted subnetwork is considered to be more stable and often used for hub nodes screening [39]. In this study, we set the parameters of degree>15 and k-core=2 to help identify key nodes in   Figure 7: The coexpression network of common targets for PNS against CHD. Circle nodes mean genes, and edges with different connecting colour represent different correlations. Genes in black circles with white stripes were submitted as query terms, while pure black circles were genes associated with query terms. Genes with larger circle indicate more correlations with other genes in the network.

Coexpression and Topology Analysis of 36 Common
the coexpression network. Results indicated that common genes involving BCL2A1, MMP9, MAPK14, S100A12, and VEGFA have degree above 15 (Figure 8(a)). Furthermore, three subnetworks were extracted from the coexpression network by k-core decomposition (Figure 8(b)). Interestingly, it is shown in our analysis that VEGFA also constituted one part of a subnetwork, showing more possibility as a key node in the whole network. More importantly, VEGFA participates in regulation of apoptosis and BCL2A1 is a wellknown antiapoptotic gene [40,41]. As a result, we selected the two key nodes of VEGFA and BCL2A1 as potential targets for next experimental verification. Additionally, cell apoptosis was also observed to reveal VEGFA and BCL2A1-involved pathway in the treatment.

PNS Upregulated VEGFA mRNA Expression in H 2 O 2 -Treated HUVECs.
To verify the predicted targets of VEGFA and BCL2A1 in PNS against CHD, 100 mol⋅L −1 H 2 O 2induced HUVECs injury was performed as the oxidative stress model in cardiovascular diseases [42,43]. In our experiment, the mRNA expressions of VEGFA and BCL2A1 were detected by real-time PCR in control, model, and treatment groups, respectively. According to the analyzing method of real-time PCR data [44], we performed the calculation of 2 −△△CT values and analyzed corresponding FC and expression level. The detail values were described in Table 4. Furthermore, statistical differences of relative expression of VEGFA and BCL2A1 in each group were also analyzed, respectively. Results showed that H 2 O 2 significantly reduced VEGFA expression compared with the control group (P=0.041); moreover, PNS promoted VEGFA markedly compared to model group (P=0.007) (Figure 9(a)). With respect to BCL2A1, the mRNA expression of BCL2A1 had the trend to be downregulated by H 2 O 2 and upregulated by PNS, but current results showed no statistical significance (P>0.05) (Figure 9(b)).

PNS Inhibited Early Apoptosis in H 2 O 2 -Treated HUVECs.
Previous reports have pointed out that VEGFA can participate in suppression of cell apoptosis [45,46]. Thus, we further carried out apoptosis analysis to find out whether VEGFA overexpression induced by PNS could inhibit cell apoptosis in HUVECs. Our results showed that HUVECs in the model group were obviously reduced and disorderly arranged, while treatment group cells grew well compared with the model group ( Figure 10). Apoptosis analysis further indicated that H 2 O 2 induced cell apoptosis, but PNS pretreatment could protect cell from apoptosis ( Figure 11). Particularly, the early apoptosis rate in model group was significantly higher than that in control group (P=0.004, Figure 12(a)). In contrast, PNS pretreatment markedly relieved H 2 O 2 -induced damages and protected HUVECs from early apoptosis when compared with model group (P=0.009, Figure 12(a)). However, no significant difference existed in the late apoptosis, although declined rates were also observed in treatment group (Figure 12(b)). These results suggested that PNS inhibited HUVECs apoptosis mainly ascribing to the suppression of early phase.

Discussion
In this study, we revealed that Ginsenoside Rd qt had better pharmacokinetics properties than the other four compounds in PNS. Thus, Ginsenoside Rd qt may be the main bioactive      Figure 12: The early (a) and late (b) apoptosis rate of HUVECs in different groups. Results were presented as mean ± SD (n = 3). Significance was defined as p < 0.05 ( * p < 0.05; * * p < 0.01).
compound that exerted protective effect in the treatment of CHD. Notably, all compounds had relatively low OB value compared with the recommended screening criteria by TCMSP server [22], which indicated that PNS was not suitable for oral administration directly. More preparations involving injection [12,47,48] and spray [49] were therefore widely applied in practice or studied for new drug development. Although PNS has shown positive druggability and efficacy in clinical application, current researches about its specific targets and mechanism against CHD are still insufficient. Targets prediction based on network pharmacology have been reported [50,51], but those predicted results are not convincing enough. Experimental verification is further needed to confirm the authenticity and reliability. As a result, we combined network pharmacology and cell experiment to predict and verify target and underlying mechanism of PNS in the treatment of CHD in our study. The results showed that VEGFA and BCL2A1 had great potential as target genes. However, only VEGFA was proved to be able to act as the target and participate in inhibition of apoptosis by cell experiment. Previous studies reported consistent results, indicating that VEGFA was an important protective factor in apoptosis of cardiomyocytes [45] and microvascular endothelial cells [46]. Moreover, VEGFA also took part in lncRNA NEAT1/miR-377/VEGFA axis to gain the protective effect in cell apoptosis [52]. It shows that VEGFA overexpression induced by PNS which inhibits cell apoptosis is reasonable. Notably, the inhibitory effect focused on the early apoptosis, rather than late phase. Therefore, our study revealed that PNS could upregulate the mRNA expression of VEGFA to participate in the suppression of apoptosis in early phase.
VEGFA was identified as a treatment target of PNS against CHD by integrating bioinformatics prediction and experimental verification. But beyond that, some other nodes in the coexpression network may also play important roles in the treatment of CHD, especially for MMP9 and MAPK14. A meta-analysis indicated that MMP9-1562C/T polymorphism is markedly associated with CHD [53]. Thus, the circulating level of MMP9 was upregulated in coronary artery disease, which may act as a risk factor for evaluating the development and severity of disease [54]. A mechanistic study further revealed that MMP-9 activation could promote cell migration to stimulate angiogenesis and preserve cardiac function [55]. As for MAPK14, also called p38MAPK , it played negative role in expression of vascular smooth muscle cell (VSMC) contractile genes and differentiation of VSMC [56]. Besides, it promoted platelet activation and served as a potential indicator of highly thrombotic lesions and no-reflow after myocardial infarction [57]. As a result, experimental verification for more potential targets needs to be performed in the future. In addition, because of the fact that the expression level of mRNA and protein may be inconsistent after treatment [58], the detection of the protein expression from VEGFA should be carried out by Western blot assay in subsequent research. Furthermore, there have been some reports [46,52] about the modulation of VEGFA by miRNA and lncRNA. Therefore, to reveal the regulatory mechanism of PNS against CHD based on lncRNA-miRNA-VEGFA axis in the apoptosis pathway is also a subject which is worthy of exploration.

Conclusion
In conclusion, our study predicted and verified that PNS can promote mRNA expression of VEGFA for the suppression of early apoptosis in vascular endothelial cells to exert curative effect for CHD. The results helped partially explain the pharmacological mechanism of PNS against CHD. Moreover, our research provided an integrated method for discovery of drug targets and corresponding mechanism.

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

Conflicts of Interest
The authors declare no conflicts of interest.