In order to collect information on SAFI main components, “salvianolic acid for injection” was used as the key word to obtain Chinese-language literature through CNKI (https://www.cnki.net/) and English-language literature through PubMed (https://pubmed.ncbi.nlm.nih.gov/). Then, the components contained in SAFI were collated manually and standardized through PubChem database (https://pubchem.ncbi.nlm.nih.gov/).

In order to obtain SAFI target information, the term “salvianolic acid for injection” was used to obtain related literature through CNKI and PubMed databases.

The potential targets of five important SAFI components-salvianolic acid B, rosmarinic acid, lithospermic acid, salvianolic acid D, and salvianolic acid Y, were collected using HERB database (Fang et al., 2020).

IS-associated genes were collected from the DisGeNET database (Piero et al., 2019) by searching the term “Ischemic stroke”. Differentially expressed genes in IS were obtained from GSE22255 dataset (downloaded from Gene Expression Omnibus (GEO), https://www.ncbi.nlm.nih.gov/geo/).

Anti-IS drugs collection were referred to “Chinese guidelines for diagnosis and treatment of acute ischemic stroke 2018” (Neurology and Society, 2018), and their targets were collected from the Drugbank database (Wishart et al., 2017).

To evaluate the efficacy of a drug on disease, the correlation of the drug targets and the genes associated with the disease from the perspective of network propagation was used (Yang et al., 2020). Briefly, the drug targets and disease genes were used as seed genes to run Random Walk with Restart (RWR) algorithm (Köhler et al., 2008) in the background network STRING, which was performed in R package dnet (version 1.1.7) with 0.75 of the restart probability. In this way, the influence score vector of the two sets of seed nodes on all nodes in the background network was obtained. Pearson correlation coefficient of the two score vectors (Cor) was then calculated, and Z-score was used to evaluate the significance of the correlation using the formula z − score = C o r − E ( C o r ) δ ( C o r ) where E(Cor) and δ(Cor) are the mean and standard deviation of the Pearson correlation coefficients between the influence score vector of drug targets and those of 1000 groups of random contrast disease genes, each of which contained the same number of randomly selected proteins as the disease seed nodes.

The similarity of two drugs was evaluated based mainly on chemical structure, comprehensive targets or cellular function fingerprints. Drug clusters with similar features were obtained using hierarchical clustering.

The chemical structure similarity was measured based on the Morgan fingerprint using RDKit in python. Molecular fingerprints encode molecular structure in a series of binary digits (bits) that represent the presence or absence of particular substructures in the molecule. Comparing fingerprints enable us to determine the similarity between two molecules. Once SMILES strings are converted to scalar fingerprints, the Tanimoto coefficient was used as the similarity score to measure the absolute similarity between two molecules(Hakime et al., 2016).

The target similarity between two drugs was measured based on comprehensive targets of compounds. Two drugs acting on same targets could be considered to have the same effect. For multi-target drugs, drugs whose targets are very close in the PPI network show similar effects (Barabási et al., 2011). Here, the network proximity index proposed by Barabasi (Joerg Menche et al., 2015) was used to explore the similarity between two drugs. Briefly, the network proximity of drug target module A and drug-target module B is defined using the separation measure as follows: S A B = d A B − d A A + d B B 2 which compares the mean shortest distance within the interactome of each targets module, ⟨d _AA ⟩ and ⟨d _BB ⟩, to the mean shortest distance ⟨d _AB ⟩ between targets module A and B. The smaller the S _AB , the closer the topological distance between the two drugs, that is, the more similar the functions of the two drugs. BNet was used as the background human PPI network.

The cellular function fingerprint similarity of two drugs was measured based on a “compound-target-cellular function” heterogeneous network using the PathSim method (Yizhou Sun et al., 2011). As described in the literature (Fg et al., 2020), in the “compound-target-cellular function” heterogeneous network, the metapath “compound-target-cellular function target-compound” of two drugs was considered to describe the linkage between two drugs. In this instance, cellular function fingerprints of compounds were described using Gene Ontology (GO) biological processes term of targets of compounds. Under the metapath framework, PathSim was developed to find peer objects in the network and to measure the similarity of peer objects based on metapaths. The “compound-target-cellular function” heterogeneous network consists of drug-comprehensive target interactions and target-cellular function relations from the Gene Ontology database (The Gene Ontology, 2017).

Given a set of N compounds to be clustered, an N × N similarity matrix was generated. Finally, the similarity matrix was used to perform hierarchical cluster, which was executed by the R package hClust.

The SRWR algorithm, which was proposed in the research of social networks for measuring trust and distrust in signed social networks (Jung et al., 2016), was used here to measure how the activation or inhibition of a seed node corresponding to a drug target, leading to the activation or inhibition of other nodes in the human signaling network.

As described in the literature (Zhao et al., 2019), suppose a signed random walker starts at one of the seed nodes s and walks on a signed network. The sign of the walker is either positive or negative, meaning that it exerts activation or inhibition to a node, respectively. At each step, the walker either moves to a randomly chosen neighbor u of the current node v or it jumps back to its starting seed node s and restarts its random walking. When the random walker goes through a negative edge, it changes its sign from positive to negative, or vice versa. Otherwise, it keeps its sign. Once the walker jumps back to its starting node, it regains its original sign. The SRWR can simulate the process that the activation or inhibition of a seed node propagates to other nodes in a signaling network.

The hypergeometric cumulative distribution was applied to quantitatively measure whether a gene set associated with IS (validation gene set) was more enriched with the genes significantly inhibited by SAFI (identified by SRWR) than would be expected by chance. We used IS-associated disease genes collected from DisGeNET as distinct validation gene sets. The detailed methods to define the p-value have been previously described (Zhao et al., 2019).

For the purpose of exploring the biological function of genes correlated with IS, GO, Kyoto Encyclopedia of Gene and Genomes (KEGG) and tissue enrichment analysis were performed using DAVID. Significance of each term was assessed with a p-value, and a term with p-value <0.05 was considered being significant.

The dataset of key common genes in IS and SAFI was obtained from the ChEMBL database (Anna et al., 2012), and processed in the Knime workflow. IC₅₀ values in Schrödinger’s Ligprep module were converted to predicted IC₅₀ values (PIC₅₀). The QSAR model based on traditional methods was generated for key common genes through Schrodinger’s AutoQSAR module.

AutoQSAR is an automated machine learning application for building, validating, and generating quantitative structure-activity relationship (QSAR) models. The process of descriptor generation, feature selection, and creation of a large number of QSAR models has been automated into a single workflow in AutoQSAR (Dixon et al., 2016). These models are built using various machine learning methods (Bayes, KPLS, MLR, PCR, PLS, and RP), and each model is scored using a novel method, considering R² and Q² respectively (de Oliveira and Katekawa, 2018). AutoQSAR randomly selects a combination of training set and test set, and allocates data to the learning set according to the same pattern. Then, it calculates physicochemical and topological descriptors and 2D fingerprints (linear, dendritic, radial or MOLPRINT2D) through the Canvas module. These methods are used in combination with machine learning methods and can be used to build various predictive models.

To evaluate the influence on PTGS1 and PTGS2 by SAFI and obtain IC₅₀ values, Pharmaron (Beijing) was commissioned to perform in vitro PTGS1 and PTGS2 enzymatic assays. The detailed information about the assays, such as the reagents, instruments, assay procedure, data analysis, and calculation is presented in Supplementary Methods.

In order to validate the influence of SAFI on PTGS2 activity, RAW264.7 macrophages were plated in 96-well plate and pretreated with different concentrations of SAFI for 1 h. Then cells were stimulated by lipopolysaccharide (LPS, 0.5 μg/ml) for 18 h to induce PTGS2 expression. Similarly, BV-2 microglia were treated by LPS(0.1 μg/ml) and different concentrations of SAFI for 24 h to induce PTGS2 expression. The prostaglandin E2 content in supernatants, generated by cyclooxygenase 2 (COX2, encoded by PTGS2) in the conversion of arachidonic acid, was assayed with enzyme linked immunosorbent assay according to the protocol of manufacturer(JINGMEI II, Lot 202107).

Values are represented as mean ± SD. Statistical significance was determined by One-way ANOVA in GraphPad Prism 8. p-value <0.05 was considered statistically significant.

In order to identify the important targets of SAFI in the treatment of IS, we first extracted its components through literature mining. It was found that salvianolic acid B, rosmarinic acid, lithospermic acid, salvianolic acid D, and salvianolic acid Y were major chemical components of SAFI, and four constituents present in the rat blood including salvianolic acid B, rosmarinic acid, lithospermic acid and salvianolic acid D were determined as effective Quality Markers of SAFI, identified by LC-MS/MS (Li et al., 2018) (Figure 2A; Table 1, Supplementary Figure S1). Subsequently, the targets of SAFI were collected from both literature mining and database searches. As a result, a total of 79 targets were collected from both CNKI and NCBI databases after searching with the keyword “SAFI” (Supplementary Table S1). Potential targets of SAFI’ main components were obtained from HERB database, and after manual removal of duplicates, 227 genes were obtained (Supplementary Table S2). The relationship between components and targets is shown in Figure 2B. Importantly, salvianolic acid B and rosmarinic acid possessed the most targets, whereas no targets were obtained from salvianolic acid Y after literature mining or database searches.

We also collected the IS-associated disease genes from DisGeNET database, and obtained a total of 1159 disease genes (Supplementary Table S3). Figure 2C indicate the overlaps of SAFI and its main components and IS-associated genes. Intersection analysis highlighted 38 common genes (Supplementary Table S4) between IS, SAFI’s components and SAFI, representing 3.28% of IS genes, 48.1% of SAFI putative targets and 16.74% of main components datasets. The PTGS2 gene was targeted by three main components of SAFI, salvianolic acid B, rosmarinic acid, and salvianolic acid D, suggesting an important role in the treatment of IS. Functional annotation of the 38 common genes was performed by enrichment analysis. As showed in Figure 2D, top GO terms were enriched in the BP category like negative regulation of apoptotic process and inflammatory response. Annotated pathways indicated that inflammatory associated pathways such as TNF signaling pathway, HIF signaling pathway, and NF-kappa B signaling pathway were significantly enriched in KEGG pathways (Figure 2E). In fact, apoptosis and inflammatory response are key characteristics of IS progression (Datta et al., 2020; Maida et al., 2020). In summary, we identified the putative target genes of SAFI correlated with IS and their association with inflammation.

In order to investigate whether SAFI has a similar mechanism of action with current IS marketed drugs, we compared them using an unsupervised clustering evaluation method based on similarity. The “Chinese guidelines for diagnosis and treatment of acute ischemic stroke 2018” (neurology and society, 2018) recommends 16 drugs to treat IS, representing five different mechanisms of action, including antiplatelet (tirofiban, clopidogrel, ticagrelor, aspirin), anticoagulant (heparin, argatroban), antilipemic (lovastatin, pitavastatin, fluvastatin, rosuvastatin, pravastatin, simvastatin), antihypertensive (labetalol, nicardipine) and neoroprotectant (citicoline, edaravone). Next, the main compounds in SAFI were compared with the recommended anti-IS drugs based on similarity of chemical structure, functions and targets. The results showed that the main compounds of SAFI were not clustered with the five types of anti-IS drugs based on chemical structure (Figure 3A), indicating that the main compounds of SAFI do not share a similar structure with the recommended anti-IS drugs. However, hierarchical clustering of compounds showed that rosmarinic acid and salvianolic acid B were clustered with aspirin both based on functions and targets (Figure 3B,C), and salvianolic acid B was also clustered with simvastatin based on functions (Figure 3B). This indicates that SAFI may share similar effects with aspirin and simvastatin.

In order to evaluate and quantify the correlations between SAFI and IS-associated genes collected from DisGeNET, a parameter named Z-score was applied. Z-score, which provides a quantitative measure for the significance of the correlation between two genes, were set as greater than 3. Z-score indicated that the targets of both SAFI and its main components have a close link with IS-associated genes, especially salvianolic acid B (Z-score = 20.994) and rosmarinic acid (Z-score = 19.023), with a total of more than 200 IS-associated genes (Table 2).

Although it is known which targets obtained from literature mining are inhibited or activated by SAFI, the directionality of action of targets obtained from databases is unknown. The inhibition or activation of a target by a drug propagates to the other nodes in the human signaling network, and this can be measured by SRWR algorithm (Zhao et al., 2019), helping in the assessment on the potential impact of the drug. To construct the human signaling network, 7118 genes were used, and 4317 were used as the giant strongly connected component (GSC) in the background of the human signaling network. IS-associated genes, IS-GEO-UP genes and IS-GEO-down genes were used as the validation sets (Supplementary Tables S3, S5, S6). The method has been previously described (Zhao et al., 2019). We defined IS-associated genes significantly inhibited by SAFI as those with negative activation scores and with absolute values that ranked within the top 10% of all the nodes obtained with SRWR algorithm. In this way, the top 362 inhibited genes were identified. Table 3 shows the number of genes in common between the top inhibited genes by each drug and the three validation gene datasets. In most cases, the p-values were much smaller than 0.05, indicating the statistically significant enrichment of inhibited genes in the validation sets. The p-values in IS-GEO-UP dataset (the upregulated genes in GSE22255 dataset, Supplementary Table S5) were similar to those in IS disease genes (from DisGeNET), suggesting that SAFI may work as an inhibitor of IS-associated disease genes. It is worth mentioning that PTGS1 is in the top 10% inhibited by SAFI via SRWR algorithm, so it is also considered one of the most important targets of SAFI.

To investigate whether the top inhibited genes by SAFI were expressed in IS target tissues, tissue enrichment analysis was conducted using DAVID. As shown in Table 4, a large number of the inhibited genes are expressed in whole brain and blood in the category of GNF_U133A_QUARTILE and UP_TISSUE, respectively (p-value < 0.05). These results suggest that SAFI could inhibit target tissue proteins to exert its anti-IS effects.

For the 362 top inhibited genes predicted by SRWR algorithm, DAVID Functional Annotation Clustering tool was used to conduct the functional annotation from GO and KEGG analysis. The results showed that the top inhibited genes are related with processes such as inflammatory response, platelet activation and cell adhesion, which are all anti-IS associated process. KEGG enrichment analysis indicated that the top inhibited genes are involved in PI3K-Akt, TNF, NF-kappa B, and HIF signaling pathways, among others, which are associated with inflammatory response (Du et al., 2020; Minghua et al., 2021); Supplementary Figure S2). These results indicate that SAFI significantly inhibits a large number of genes associated with IS and inflammation, further validating its effect on IS.

Based on the intersection analysis and SRWR results, PTGS1 and PTGS2 show big potential as important genes regulated by SAFI. We predicted the potential binding abilities between SAFI/main components and PTGS1/2 using the QSAR model. The PTGS1/2 dataset was downloaded from the ChEMBL database and processed in the Knime workflow. We prepared molecules with IC₅₀ values in Schrödinger’s Ligprep module, and converted IC₅₀ values to PIC₅₀ values. The QSAR model based on traditional methods was generated for PTGS1/2 through Schrodinger’s AutoQSAR module. R² (correlation coefficient) and Q² (cross validation coefficient) represent the availability of the model (the better model the closer to 1). By default, the dataset was divided into a 75% training set (PTGS1: 1217 molecules, PTGS2: 2634 molecules) and a 25% test set (PTGS1: 406 molecules, PTGS2: 879 molecules). We used two models with Q² > 0.5 in top10 to predict the activity of the compound on the target. The model score of PTGS1/2 with predictive ability is shown in Table 5, and the scatter plot in Figure 4 shows the performance of the model in predicting the experimental binding affinity of the learning set. We also validated the model with positive molecules, and the results are shown in Table 6. The prediction results showed that our model can identify molecules that are active to the corresponding target, and the difference in activity is no more than an order of magnitude. The activity values of SAFI predicted by the QSAR model are shown in Table 7. The results indicated that salvianolic acid B, lithopermic acid, salvianolic acid Y, and salvianolic acid D may possess inhibitory effect on PTGS1 and PTGS2 like two anti-IS drugs, aspirin and NS-398.

Arachidonic acid metabolism plays an important role in acute ischemic syndromes affecting the coronary or cerebrovascular territory, whereas, cyclooxygenase, including PTGS1 and PTGS2, is the key enzyme of the arachidonic acid metabolism (Santovito et al., 2009; Cipollone et al., 2010). Our computational prediction results show that SAFI has an inhibitory effect on these two enzymes. To validate the effect of SAFI on PTGS1 and PTGS2, we determined how their enzymatic activity was modulated. The enzymatic assay showed that SAFI had a strong inhibition on PTGS1 and PTGS2, with IC₅₀ of 0.04 μg/ml (Figure 5A) and 0.03 μg/ml (Figure 5B) respectively. The inhibition curves of positive controls SC-560 and celecoxib on PTGS1 and PTGS2 can be seen in Supplementary Figure S3. In order to validate the anti-inflammation effect of SAFI, we determined the content of prostaglandin E2 induced by LPS in RAW264.7 macrophages and BV-2 microglia, respectively. As shown in Figure 5C, SAFI, at a concentration of 250 μg/ml, had a significant inhibition on the production of prostaglandin E2 both in the RAW264.7 macrophages and BV-2 microglia. These indicate that SAFI inhibits the production of prostaglandin E2 in a cell-based assay, and these offer a mode of action for its anti-inflammation effect.

Intravenous administration of tissue plasminogen activator and endovascular treatment are currently used to recanalize the blood flow in acute IS. However, reperfusion leads to a highly harmful reactive oxygen species production, generating oxidative stress, which is responsible for most of the ischemia-reperfusion injury and brain tissue damage (Sofía et al., 2020). The positive effect of SAFI on ischemia reperfusion injury has been widely reported in the literature (Tang et al., 2015; Wang et al., 2017; Huang et al., 2019), although its targets and mechanism of action remain to be discovered.

In this study, we have systemically analyzed the mechanism of action of SAFI on IS using the network pharmacology combined computational prediction and experimental validation. We identified 38 genes shared by datasets of SAFI targets, SAFI main components targets and IS-associated genes. GO and KEGG enrichment analysis of this common shared gene set indicates their association with inflammation. Two genes associated with inflammation (PTGS1 and PTGS2) were preliminarily validated in vitro. SAFI inhibited PTGS1 and PTGS2 activity in a dose-dependent manner and inhibited the production of prostaglandin E2 induced by LPS in RAW264.7 macrophages and BV-2 microglia. Cytokines, like tumor necrosis factor (TNF), modulate tissue injury in experimental stroke and affect infarct evolution, making them good targets for potential future stroke therapy (Lam Be Rtsen et al., 2012). Nuclear factor kappa B (NF-κB) is a key regulator of the inflammatory process, regulating expression of proinflammatory and proapoptotic genes in ischemia-reperfusion injury of the brain and representing a therapeutic target of IS (Howell and Bidwell, 2020). HIF-1α and Notch-1 engage pro-inflammatory and apoptotic signaling pathways, thus promoting neuronal cell death in IS (Cheng et al., 2014). SAFI and its main component, Salvianolic acid B, have been reported to inhibit IL-1β, IL-6 NF-κB, TNF, and HIF-1 (Ming-chao et al., 2016; Wang et al., 2016; Wang et al., 2017; Zhuang et al., 2017; Fan et al., 2018; Zhao et al., 2020). Interestingly, as one of the main components, no potential targets of salvianolic acid Y were obtained from HERB database and literature mining. The only report related to salvianolic acid Y was about its antioxidant effect (Jun et al., 2015). In view of the similar structure with salvianolic acid B, the potential anti-inflammatory activity of salvianolic acid Y deserves to be mined in the future. QASR model prediction analysis indicates that SAFI has the potential to bind to PTGS1 and PTGS2, with binding as good as the positive drugs, aspirin and NS-398. Drug similarity analysis also indicates that SAFI and aspirin may share a common set of targets. This highlights once more that the effect of SAFI on IS may be due to its anti-inflammatory effect by inhibiting the activity of PTGS1 and/or PTGS2. Aspirin has been frequently prescribed to prevent cardiovascular disease due to its analgesic, anti-pyretic, anti-inflammatory and antithrombotic qualities. However, its side effects, such as internal bleeding and gastrointestinal damage, should also be noted when long-term medication is used (Sostres and Lanas, 2011; Fiala and Pasic, 2020). SAFI has been widely used to treat IS clinically in China, mostly for mild to moderate cases, with no gastrointestinal side effects. Although SAFI has been reported to inhibit the expression level of PTGS2, but not PTGS1, the relationship between SAFI, IS and PTGS1/2 was not confirmed (Wang et al., 2017). Therefore the interaction between SAFI and PTGS1/PTGS2 in IS are still not clear in the current literature. Overall, our network pharmacology approach indicates that the beneficial effect of SAFI on IS is likely due to its anti-inflammatory properties via direct inhibition of PTGS1 and PTGS2. It should be noted that PTGS1 and PTGS2 are involved in inflammation not only via the arachidonic acid (AA) metabolic pathway, but also through the maintenance of cerebral blood flow, synaptic plasticity, and cerebrovascular regulation (Niwa et al., 2000; Chen and Bazan, 2005; Jayaraj et al., 2020). Due to the advantages of multiple components, multiple targets and multiple pathways, other potential effects of SAFI in the treatment of IS should be further explored.

In this study, we have used network pharmacology combined with computational prediction to elucidate the molecular mechanisms of SAFI acting on IS, and demonstrated the feasibility of the approach in investigating its mode of action. A similar approach can be used for a variety of herbal medicine used in TCM for a variety of diseases, thus contributing to its modernization.

The emergence of TCM databases such as HERB, TCMSP and others makes relatively easy to obtain the chemical components and corresponding targets of TCM, saving a lot of time and cost of experiments (Ru et al., 2014; Fang et al., 2020). However, the chemical components in the prescription and their pharmacological characteristics, such as their absorption, have not been fully evaluated in many of them, which may lead to some bias when using enrichment analyses.

Propagation-like algorithms, such as SRWR, have been used in the identification of proteins or genes influenced by drugs or diseases, and can simulate the spread of active compound-induced activation or inhibition to a group of targets on a particular signaling network(Zhao et al., 2019). Using this algorithm, we have identified the top 10% genes inhibited by SAFI, and enriched them by pathways and biological process analysis. The results illustrate the mechanism by which SAFI likely acts when used to treat IS: inhibition of a group of proteins associated with inflammation. However, the human signaling network may not include some unknown protein interactions with orientation and/or patterns, and this may lead to some loss of information in subsequent identification analysis of inhibited and activated genes. Given that, more studies using high-throughput data such as RNA sequencing, single-cell sequencing and other computational models to elucidate the mechanisms of SAFI treating IS are warranted, which may provide us more rich and accurate information on IS(Olsen and Baryawno, 2018; Peng et al., 2021; Zheng et al., 2022).

QSAR has been applied for decades in the development of relationships between physicochemical properties of chemical substances and their biological activities in order to obtain a reliable statistical model for prediction of the activities of new chemical entities (Verma et al., 2010). Since the study of QSAR is based on a series of assumptions using known compounds, this methodology does not consider the possible interactions of complex components in SAFI, or other TCMs, and there may be deviations between predicted activities and actual conditions.

In summary, in this study we deduce the mode of action of SAFI on IS from the levels of the components, targets, pathways, network, and activity prediction. SAFI putative targets are significantly enriched in several pathways associated with inflammation, which is critical in IS. PTGS1 and PTGS2 were found in a dataset of genes shared between SAFI and its components’ targets and IS-associated genes, and their binding activities were further predicted by QSAR. The effect of SAFI on the enzymatic activities of PTGS1 and PTGS2 confirmed the data deduced from network pharmacology combined with computational prediction. Due to the limitations of network pharmacology and the complicated pathogenesis progress of IS, more extensive and in-depth research is still needed to corroborate these data to elucidate other action mechanisms of SAFI in the treatment of IS.