A novel bioinformatics strategy to uncover the active ingredients and molecular mechanisms of Bai Shao in the treatment of non-alcoholic fatty liver disease

Introduction: As a new discipline, network pharmacology has been widely used to disclose the material basis and mechanism of Traditional Chinese Medicine in recent years. However, numerous researches indicated that the material basis of TCMs identified based on network pharmacology was the mixtures of beneficial and harmful substances rather than the real material basis. In this work, taking the anti-NAFLD (non-alcoholic fatty liver disease) effect of Bai Shao (BS) as a case, we attempted to propose a novel bioinformatics strategy to uncover the material basis and mechanism of TCMs in a precise manner. Methods: In our previous studies, we have done a lot work to explore TCM-induced hepatoprotection. Here, by integrating our previous studies, we developed a novel computational pharmacology method to identify hepatoprotective ingredients from TCMs. Then the developed method was used to discover the material basis and mechanism of Bai Shao against Non-alcoholic fatty liver disease by combining with the techniques of molecular network, microarray data analysis, molecular docking, and molecular dynamics simulation. Finally, literature verification method was utilized to validate the findings. Results: A total of 12 ingredients were found to be associated with the anti-NAFLD effect of BS, including monoterpene glucosides, flavonoids, triterpenes, and phenolic acids. Further analysis found that IL1-β, IL6, and JUN would be the key targets. Interestingly, molecular docking and molecular dynamics simulation analysis showed that there indeed existed strong and stable binding affinity between the active ingredients and the key targets. In addition, a total of 23 NAFLD-related KEGG pathways were enriched. The major biological processes involved by these pathways including inflammation, apoptosis, lipid metabolism, and glucose metabolism. Of note, there was a great deal of evidence available in the literature to support the findings mentioned above, indicating that our method was reliable. Discussion: In summary, the contributions of this work can be summarized as two aspects as follows. Firstly, we systematically elucidated the material basis and mechanism of BS against NAFLD from multiple perspectives. These findings further enhanced the theoretical foundation of BS on NAFLD. Secondly, a novel computational pharmacology research strategy was proposed, which would assist network pharmacology to uncover the scientific connotation TCMs in a more precise manner.


Introduction
Non-alcoholic fatty liver disease (NAFLD) is defined as unusual hepatic fat accumulation in the absence of excessive alcohol consumption and any other conditions that may lead to hepatic steatosis (Chalasani et al., 2012).Patients suffered from NAFLD always presented liver dysfunction and are confronted with higher risk of liver-related morbidity (Mantovani et al., 2020).It was reported that NAFLD has been a leading cause of cirrhosis, hepatocellular carcinoma, and some other irreversible liver diseases (Powell et al., 2021).In addition, NAFLD also increased the risk of cardiovascular diseases (Yoshitaka et al., 2017), chronic kidney diseases (Mantovani et al., 2022), and malignancy (Rinella and Sanyal, 2016) significantly.In fact, NAFLD is a newly discovered disease.It has not been described until 1980s.However, in just a couple of decades, it has become a global health issue because of its high prevalence (Siegel et al., 2012;Perumpail et al., 2017).Epidemiological research showed that the morbidity of NAFLD in many countries was already close to 30% (Younossi et al., 2016;Estes et al., 2018).Although the liver disease specialists have made a tremendous amount of effort to seek effective therapeutic drugs for NAFLD during the past decades, there is no approved therapy until now.Therefore, scientific researches focused on exploring the candidate drugs of NAFLD deserve to be paid a great deal of attention.
Paeoniae Radix Alba, the dried roots of Paeonia lactiflora Pall., is a well-known herbal medicine and widely used in Asia for a long history.Owning to it presented very obvious effect on some complex diseases, it has been recorded in the Pharmacopoeias of many countries, including but not limit to China and Japan.In China, Paeoniae Radix Alba is known as Bai Shao (BS) and mainly used to regulate menstruation, antiperspirant, relieve pain and protect liver (Chinese Pharmacopeia Commission, 2020).In recent years, with the social and economic burdens caused by NAFLD increasing continually, the hepatoprotective effect of BS gradually gained more and more attention.It has been demonstrated that BS and its extracts have multiple beneficial effects on liver, including drug/ chemical-induced liver injury attenuation, anti-fatty liver, antihepatitis, anti-hepatic fibrosis, cholestasis alleviation, and hepatocellular carcinoma inhibition (Ma et al., 2020).Especially for the action of preventing and treating NAFLD, a great deal of evidence has been provided in previous publications.Sun et al. established fatty liver model by intervening HepG2 cells with oleic acid and explored the anti-fatty liver action of BS.As a result, they found that compared with the model group, the levels of TC, TG, ALP, and ALT of the pretreatment group and the treatment group were decreased significantly (Sun et al., 2022).The anti-fatty liver effect of BS in animal models has also been confirmed (Zhang et al., 2015;Ma et al., 2016).For example, Ma et al. found that BS treatment can obviously ameliorate the histopathological and biochemical changes in NAFLD rats induced by high-fat diet via inhibiting lipid ectopic deposition (Ma et al., 2017).In addition, some clinical scholars have attempted to treat NAFLD by combining polyene phosphatidylcholine and the total glycosides of BS.Excitedly, the total effective rates of the combination group increased by 47.06% in comparison with the polyene phosphatidylcholine treatment group.Besides, more significant improvements on clinical symptoms, liver function, and blood lipids were also observed in the combination group (Tian et al., 2014).In summary, the anti-NAFLD effect of BS has been well documented.However, which compounds are the active ingredients and how they produce these beneficial effects have not been fully elucidated.Therefore, more research focused on revealing the mechanisms of BS against NAFLD is needed.
Different from the western medicine of "one target, one drug," as complex mixtures of many bioactive ingredients, Traditional Chinese Medicines (TCMs) always produce their effects via multi-component and multi-target, making it to be a challenging work to unveil the material basis and mechanism of TCM systemically (Li, 2016).In 2007, network pharmacology, a new discipline, was proposed by Hopkins AL (Hopkins, 2008).Unlike the reductionism research strategies which treat both drugs and targets in isolation, network pharmacology attempts to illustrate the interactions between drugs and biological systems in a systematic manner (Kibble et al., 2015).The perspective of network pharmacology is in accord with the holistic theory of TCM, making it possible to deeply understand the scientific connotation of TCM.Furthermore, compared with the conventional methods, network pharmacology has the advantages of high-efficacy, effort-saving, and low cost.During the past decade, network pharmacology has attracted increasing attention and has become a popular tool in the field of TCM research (Li and Zhang, 2013).However, in the current state of the art, there still existed some problems in the field of network pharmacology.For example, in a research from Henan University of Traditional Chinese Medicine, the authors attempted to illustrate the material basis of Psoraleae Fructus-induced liver injury.Network pharmacology analysis indicated that a total of 25 compounds were associated with the hepatotoxicity of Psoraleae Fructus (Cao et al., 2022).However, we found that not all of these substances were harmful to the liver.Inversely, quite a few of them were liver-friendly and considered to be promising hepatoprotectors, such as daidzein, biochanin A, and delphinidin.Numerous data showed that these compounds have the effects of anti-NAFLD (Kim et al., 2011;Fan et al., 2021), anti-hepatic fibrosis (Domitrović and Jakovac, 2010;Breikaa et al., 2013b), anti-hepatocellular carcinoma (Xiao et al., 2020;Zakaria et al., 2021;Sun et al., 2023) and attenuating liver injury induced by multiple chemical agents or drugs (Domitrović and Jakovac, 2010;Breikaa et al., 2013a;Yu et al., 2020).In addition, both in-vivo and in silico toxicity studies indicated that daidzein is a safe natural substance without hepatotoxicity (Laddha et al., 2022).In another network pharmacology study from the University of Hong Kong, kaempferol and thymol exhibited the largest number of liver injury targets connections and considered to be play a crucial role in Xiao Chai Hu Tang-induced liver injury (Hong et al., 2017).However, both kaempferol (Ren et al., 2019;Xiao et al., 2022;Akiyama et al., 2023) and thymol (Jafari et al., 2018;Guo et al., 2023;Lahmi et al., 2023) are promising and safe agents for treating multiple liver diseases.They have never been implicated in the adverse events of drug-induced liver injury.In fact, the similar phenomena were broadly existed in network pharmacology researches (Wang et al., 2017;Zheng et al., 2020;Li et al., 2021).Cases mentioned above suggested that not all of the compounds connected with liver injury targets were hepatotoxic substances.Inversely, they could be liver-friendly substances.Likewise, we speculated that the compounds connected with liver disease targets were the mixtures of hepatotoxic and hepatoprotective substances.Therefore, in the network pharmacology researches focused on revealing the material basis and mechanism of TCMinduced hepatoprotection, establishing a method to identify the hepatoprotective ingredients may help to improve the precision of the results.
In our previous studies, a great deal of efforts have been paid to TCM-induced hepatoprotection (He et al., 2022).Firstly, a large scale of dataset, including 677 hepatoprotective phytoconstituents and 205 hepatoprotective TCMs, was established by conducting a Flowchart of the overall methodology.BS, Bai Shao; DL, drug likeness; OB, oral bioavailability; NAFLD, non-alcoholic fatty liver disease; PPI, proteinprotein interaction.CPM-TCMIHP indicated the computational pharmacology method for identifying the material basis of TCM-induced hepatoprotection which were detailed in Figure 15.No report × M1, M2, and M3 were the abbreviations of module 1, module 2, and module 3, respectively.The potential hepatoprotective ingredients identified by each module were marked with "√," and the unique hepatoprotective ingredients identified by each module were highlighted in red.In addition, those compounds which were predicted as non-hepatoprotective ingredients by module three were marked with "×."The substructures highlighted in red were the representative substructures for hepatoprotective activity.
Frontiers in Pharmacology frontiersin.orgcomprehensive literature retrieval.Then molecular network technique was adopted to construct a hepatoprotective "TCM-ingredient" network.In addition, by incorporating the use of multiple machine learning algorithms, an in silico model for predicting the liver protecting activity of phytoconstituents was developed for the first time.All of these studies mentioned above laid a foundation for identifying the material basis of TCM-induced hepatoprotection.In present study, we attempted to establish a computational pharmacology method to identify hepatoprotective ingredients from TCMs by integrating our previous studies.Then, the proposed method was utilized to explore the material basis and mechanism of BS against NAFLD by combining with network pharmacology, molecular docking, and molecular dynamics simulation techniques.We hope this work would be helpful to understand the scientific connotation of Bai Shao on NAFLD in a more precise manner.The detailed research project was displayed in Figure 1.

Candidate compounds of BS
A total of 97 unduplicated compounds were extracted from TCMSP, ETCM, and TCMID databases.After investigating the   The common targets between BS and NAFLD.Tissue-specific expression pattern analysis of the common targets between BS and NAFLD.If a gene expressed highest in the liver among the 84 organs, it was ranked 1. Inversely, if it expressed lowest in the liver, it was ranked 84.
pharmacokinetic parameters of DL and OB, twenty-three compounds were found to satisfy the inclusion criteria (Table 1, Compound 5-Compound 27).These compounds were defined as candidate compounds and further analyzed.

Identification of the potential hepatoprotective ingredients in BS
As shown in Table 1, a total of 23 potential hepatoprotective ingredients (Compound 1-Compound 23) were identified by the computational pharmacology method proposed.Module 1, module 2, and module 3 separately identified 12, 6, and 10 potential hepatoprotective ingredients.The unique hepatoprotective ingredients provided by these three modules were 7, 1, and 10 in number, respectively.Such a phenomenon indicated that there existed complementary effects among module 1, module 2, and module 3.In other words, it was necessary to adopt a multi-module fusion strategy to clarify the material basis of TCM-induced hepatoprotection.In fact, we also attempted to collect the hepatoprotective ingredients in BS by retrieving concerned literature.As a result, only three hepatoprotective ingredients, including paeoniflorin (Kaur et al., 2018), albiflorin (Jeong et al., 2017), and palbinone (Li et al., 2023), were collected.Obliviously, compared with the method of retrieving literature, the computational pharmacology method proposed in the current study was able to identify the material basis of TCM-induced hepatoprotection more comprehensively.
Considering that the data included in module one and module 2 has been well documented in previous studies, the reliability of the results provided by these two modules was not discussed here.To investigate the reliability of the results provided by module 3, we have attempted to find some important evidence from the literature.
As a result, among those 10 potential hepatoprotective ingredients identified by module 3, a total of six ingredients (compound 14-compound 18 and compound 22) have been confirmed to be liver-friendly substances.A range of liver-protecting activities were involved, including against chemical and drug induced liver injury (Van Puyvelde et al., 1989;Oh et al., 2002;Chang et al., 2020;Wang F. et al., 2023), anti-proliferation (Park et al., 2020), anti-lipid peroxidation (Vinholes et al., 2014), and anti-hepatic fibrosis (Li et al., 2023).For those four non-hepatoprotective ingredients (compound 24-compound 27) identified by module 3, although there was no clear evidence that they are non-hepatoprotective ingredients, no one of them was reported to be potential liver protective drugs.In addition, in our prior study (He et al., 2022), we identified a series of representative substructures (RSs) for the hepatoprotective activity of phytoconstituents.Here, structure matching between the RSs and the potential hepatoprotective ingredients was performed.In consequence, among those 10 potential hepatoprotective ingredients identified by module 3, a total of 3 (compound 14, 17, and 21) ingredients were found to contain RSs.Conversely, for those four non-hepatoprotective ingredients identified by module 3, no one of them was found to contain RSs.In summary, both the literature retrieval and the RSs matching results demonstrated that the data provided by module three was highly reliable.The structural classification of the potential hepatoprotective ingredients in BS was also investigated.As illustrated in Figure 2, monoterpene glucosides, phenolic acids, flavonoids, and triterpenes were found to be the major structure types.

"Ingredient-target" network of BS on NAFLD
Firstly, Venn diagram analysis was performed to identify the shared targets between BS and NAFLD.As demonstrated in Figure 4, a total of 78 common targets were identified.Tissuespecific expression pattern analysis showed that there were Subcellular localization analysis of the common targets between BS and NAFLD.
36 targets that have the liver expression levels ranked in the top 10 among the 84 organs of human.Fifty-one targets have the liver expression levels ranked among the top 20 (Figure 5).In summary, most of the shared targets were highly expressed in the liver tissue, which may lay a foundation for BS alleviate NAFLD.To further explore the potential biological function of the common targets, subcellular localization analysis was performed.As shown in Figure 6, the common targets were enriched in a variety of cellular compartments, and the top five cellular compartments were nucleoplasm, cytosol, vesicles, plasma membrane, and golgi apparatus, respectively.In fact, all of these five cellular compartments have been demonstrated to be closely related to the development of NAFLD (Qiu et al., 2013;Wang Y. et al., 2023;Lu et al., 2023;Sherman et al., 2023).For example, Lipin proteins, including Lipin1, Lipin2, and Lipin3, play crucial roles in lipid metabolism.It has been reported that Lipin3 heterozygous knockout mice was more easily to suffered from NAFLD.Further mechanistic study suggested that such a situation was associated with the abnormal distribution of Lipin1 in cytosol and nucleoplasm (Wang F. et al., 2023).
Finally, the "ingredient-target" network of BS on NAFLD was constructed by inputting the common targets and the common targets related compounds into the network visualization tool Gephi (version 0.9.2).As shown in Figure 7, the network consisted of 90 nodes and 147 edges.A total of 12 potential hepatoprotective ingredients were involved, including four flavonoids (kaempferol, (+)-catechin, epigallocatechin, and astragalin), three triterpenes (oleanolic acid, β-sitosterol, and betulinic acid), three phenolic acids (gallic acid, benzoic acid, and paeonol), and two monoterpene glucosides (paeoniflorin and benzoylpaeoniflorin).The anti-NAFLD effect of BS may be largely attributed to these ingredients.In addition, a total of 78 NAFLD-related targets were included in the "ingredienttarget" network.To distinguish the known targets from the putative targets, we highlighted the edges by grey and yellow, respectively.It is not difficult to find that more than 90% of the "Ingredient-target" network of BS on NAFLD.The red nodes represented the potential hepatoprotective ingredients in BS, and the green nodes represented the NAFLD-related targets.In addition, the known targets and the putative targets were connected with the potential hepatoprotective ingredients by grey and yellow edges, respectively.edges were showed in grey, indicating that the "ingredient-target" network developed has a high reliability.

PPI network analysis
By mapping the 78 common targets identified in Section 2.5 into the STRING database, we acquired the PPI network of BS on NAFLD (Figures 8A).This network consisted of 76 nodes and 1,111 edges.The red, purple, and green nodes represented the targets that have the liver expression levels ranked in the top 10, top 20, and top 84 among the 84 organs of human, respectively.Of note, the tissue-specific expression data of CD274 was unavailable.Therefore, here, we highlighted it by blue.Modular analysis found that the PPI network can be divided into three functional modules (Figure 8B-D).As the major functional module of BS on NAFLD, functional module one consisted of 34 nodes and 494 edges.Functional module three is the minimal functional module which consisted of three nodes and three edges.In addition, functional module two including 11 nodes and 32 edges.

KEGG pathway enrichment analysis
In the previous section, we identified three functional modules of BS on NAFLD.Herein, by mapping the targets included in each functional module into the DAVID platform, we obtained 23 NAFLD-related KEGG pathways in total.All of these 23 pathways have been well documented to be associated with the development of NAFLD.
As shown in Figure 9A-C, a total of 12 KEGG pathways were enriched by functional module 1. Functional module two and functional module 3 separately enriched nine and 2 KEGG pathways.Then we sorted these pathways based on the adjusted p-value.In consequence, we found that the top five pathways were IL-17 signaling pathway, TNF signaling pathway, Non-alcoholic fatty liver disease, NF-kappa B signaling pathway, and Insulin resistance, respectively.Interestingly, the pathway of Nonalcoholic fatty liver disease was enriched significantly with the adjusted p-value ranked third, confirming the moderating effect of BS on NAFLD to a great extent.In addition, it's worth noting that all of the top four pathways were enriched by functional module 1, suggesting that it may be the major functional module of BS on NAFLD.
To elucidate the regulation effects of those three functional modules on NAFLD from the perspective of biological process.Firstly, we summarized the major biological processes associated with NAFLD by conducting a systematical literature retrieval.As showed in Figure 9F, we collected 12 NAFLD-related biological processes in total.Then, for each KEGG pathway enriched above, the related biological processes were annotated based on KEGG pathway database.In addition, we also retrieved the scientific literature database, by which we attempted to obtain the latest and most complete biological process data.Finally, Chord Plot was plotted to display the interactions between the KEGG pathways and the NAFLD-related biological processes (Figure 9E,F).As a result, we found that the major biological processes regulated by functional module one were inflammation, apoptosis, and lipid metabolism.For functional module 2, a total of four major NAFLD-related biological processes were observed, including glucose metabolism, lipid metabolism, mitochondrial dysfunction, and apoptosis.Functional module three was found to involve in the biological process of glutathione metabolism.It is not difficult to found that these three functional modules regulated NAFLD through different biological processes in a synergetic manner.
By mapping the BS-related targets into the Non-alcoholic fatty liver disease pathway, we attained a detailed mechanism map of BS against NAFLD.As illustrated in Figure 10, BS may relieve NAFLD through three major paths listed as follows.Firstly, it's well known that abnormal glucose and lipid metabolism was one of the major inducements of NAFLD (Chao et al., 2019).IRS-1/2, PPAR-α, and PPAR-γ have been proved to be essential regulatory factors of glucose and lipid metabolism (Chao et al., 2019).Therefore, ameliorating glucose and lipid metabolism via regulating IRS-1/2, PPAR-α, and PPAR-γ may be one of the major mechanisms of BS on NAFLD.Secondly, as one promising therapeutic strategy of NAFLD, decreasing hepatocyte apoptosis has gained increasing attention in recent years (Canbay et al., 2005).The significant positive regulatory effect of TNF-α on hepatocyte apoptosis has been well documented by many pharmacologists (Ezquerro et al., 2019).Thus, inhibiting the pro-apoptosis effect of TNF-α and decreasing hepatocyte injury may be another important mechanism of BS on NAFLD.Thirdly, inflammation is an important pathological change from steatosis to hepatitis (Patel and Mueller, 2024).This pathological change was significantly associated with the increase of expression levels of proinflammatory cytokines and inflammatory cytokines (Rao et al., 2018).Therefore, inhibiting the overexpression of IL-1, IL-6, IL-8 and TNF-α may be the third path of BS on NAFLD.In addition, it has been reported that TGFβ1 is a powerful fibrogenic cytokine (He et al., 2016).Depending on the data displayed in Figure 10, we speculated that some ingredients in BS may be TGFβ1 agonists which can alleviate liver fibrosis to a certain extent.In summary, we believed that the anti-NAFLD effect of BS should be attributed to multi-target, multi-pathway, and multi-biological process.
In addition, we investigated the reliability of the mechanisms of BS on NAFLD attained above via a literature-based method.Firstly, we summarized the molecular mechanisms of BS extracts induced hepatoprotection by conducting a systemic literature retrieval (Table 2).After that, consistency estimate between our results and the literature reports was performed.As a result, in the level of pathway, among those 23 pathways enriched by the DAVID platform, at least 11 pathways have been proved to be associated with BS-induced hepatoprotection, including P4, P6, P7, P11, P12, P13, P14, P15, P17, P18, and P22.Besides, at least five pathways (P4, P7, P13, P14, and P17) were reported to play essential roles in the process of BS alleviates NAFLD.In the level of biological process, among those 12 biological processes identified in Figure 9D,E, at least eight biological processes (F1, F3, F4, F5, F6, F7, F9, and F11) have been demonstrated to be significantly associated with the anti-NAFLD effect of BS.In the level of gene, among the pathway of NAFLD, a total of nine genes were found to be the potential targets of BS (Figure 10).According to the available literature, we found that eight out of these nine genes can be significantly regulated by BS, including IL-6, TNF-α, IRS1, PPAR-γ, IL-1β, JUN, TGF-β1, and PPAR-α (Table 2).To be specific, PPAR-α can be upregulated by BS, and the other seven genes can be downregulated by BS.All of above mentioned indicated that our results were highly in accord with the The mechanism map of BS relieve NAFLD.The genes regulated by BS were highlighted in pink explosive shape.FFAs, free fatty acids.
Frontiers in Pharmacology frontiersin.org13 He et al. 10.3389/fphar.2024.1406188data recorded in the literature, suggesting a higher reliability of our method.

Key targets of BS on NAFLD
Here, we ranked the genes in the PPI network based on CytoHubba.A total of three algorithms were adopted, including Degree, MCC, and EPC.For each algorithm, the top 10 genes were extracted, and the common genes among these three algorithms were defined as hub genes.As demonstrated in Figure 11A-D, a total of seven genes were identified as hub genes, including hypoxiainducible factor 1-alpha (HIF1A), interleukin-1 beta (IL1-β), interleukin-6 (IL6), proto-oncogene c-Jun (JUN), peroxisome proliferator-activated receptor gamma (PPAR-γ), toll-like receptor 4 (TLR4), and tumor necrosis factor (TNF).
To further investigate the relationship between the hub genes and NAFLD, microarray data analysis was performed (Figure 11E,F).Consequently, a total of 637 differentially expressed genes were identified in the group of nonalcoholic steatohepatitis, including 271 downregulated and 366 upregulated genes.In the simple steatosis group, we observed 197 downregulated ID, of the pathway and biological process corresponded to the ID, in Figure 9. DILI, drug-induced liver injury; TGP, total glucosides of paeony; WEP, water extract of paeony; HCC, hepatocellular carcinoma.
and 407 upregulated genes, respectively.These differentially expressed genes may play critical roles in the occurrence and development of NAFLD.Further analysis found that only three hub genes named IL1-β, IL6, and JUN were included in the datasets of differentially expressed gene.These three hub genes were defined as the key genes of BS on NAFLD and further analyzed.In fact, it has been demonstrated that inflammation and apoptosis were important pathogenesis of NAFLD which can be activated by the overexpression of IL1-β, IL6, and JUN (Li et al., 2013;Zhang J. et al., 2022;Deng et al., 2022).Interestingly, as shown in Table 2, the inhibitory effect of BS on these three genes has been well documented in previous studies, suggesting the importance of these three genes in the process of BS against NAFLD.
To answer the question that whether there exists significant interaction between the key genes identified above and the potential hepatoprotective ingredients in BS, molecular docking analysis was carried out.Here, the index of binding energy was calculated to evaluate the binding capacity between the ligand and the receptor.The specific docking parameters were provided in Table 3.Generally, a lower value of binding energy indicates a stronger interaction between the ligand and the receptor.It has been demonstrated that there may exist significant binding capacity between the ligand and the receptor when the binding energy between them is less than −1.2 kcal/mol (Ma et al., 2022).As shown in Figure 12, the binding energies between the potential hepatoprotective ingredients and the key genes ranged from −9.7 kcal/mol to −4.3 kcal/mol.Obviously, the binding energies of all studied docking complexes were within the acceptable range.Therefore, we can claim that the hepatoprotective ingredients have great potential to bind with the key genes and affect their biological effects.Among all of the potential hepatoprotective ingredients, monoterpene glucosides, triterpenes, flavonoids, and coumarins were found to have lower values of binding energy with the key genes, indicating that there may exist stronger binding affinity between these ingredients and the key genes.In addition, compared with the genes of IL1-β and  IL6, JUN exhibited stronger binding affinity with the hepatoprotective ingredients.So, we speculated that JUN may play a non-negligible role in the process of BS alleviate NAFLD.
The stability of the protein-ligand complexes formed by the hepatoprotective ingredients and the key genes were investigated based on the method of MD simulation.Here, for each key gene, the docking complexes with the lowest binding energy was selected to perform MD simulation (Figures 13A,C,E).Considering that paeoniflorin is the major active ingredient of BS.Therefore, the docking complexes containing paeoniflorin were also subjected to MD simulation (Figures 13B,D,F).As a result, although there existed some small fluctuations, the RMSD profiles of all studied complexes were stable in nature (Figure 14A-C)).Specifically for the complexes of JUN-benzoylpaeoniflorin and JUN-paeoniflorin, the RMSD always kept at the value of around 2 Å over the entire course of MD simulation.For the RMSF profiles (Figure 14D-F), except for the end regions and several loop regions, the RMSF values for most of the amino acid residue of the studied genes were within acceptable range with fluctuation range less than 3 Å.Furthermore, we also calculated binding free energy for each selected docking complex by utilizing the method of MMPBSA.As shown in Table 4, the binding free energy of the docking complexes were found to lie between −31.66 and −17.15 kcal/ mol.In summary, all of the results of the MD simulations mentioned above demonstrated that the studied docking complexes has higher stability.

Discussion
In this work, an advanced network pharmacology approach was proposed to uncover the scientific connotation of BS against NAFLD.Compared with the method of conventional network pharmacology, the research strategy adopted in the current study has the advantages of more precision.Generally, network pharmacology identifies the material basis of drugs based on the connections between the components and the disease targets, and the components exhibited higher number of disease targets connections were defined as the potential material basis.However, numerous network pharmacology research cases showed that the components connected with the disease targets were the mixtures of beneficial and harmful substances rather than the real material basis.Herein, taking TCM-induced hepatoprotection as a case, we established a novel computational pharmacology-based method to identify the material basis of TCMs, by which we attempted to improve the precision and actual application value of network pharmacology.The proposed computational pharmacology method consisted of three modules.These three modules were developed based on molecular network, database retrieval, and structure activity relationship techniques, respectively.They possess different characteristics and complement each other.Compared with the method of literature retrieval, our computational pharmacology-based method could identify the material basis of TCM-induced hepatoprotection more comprehensively and more efficiently.Based on the computational pharmacology method proposed, a total of 23 ingredients in BS were identified as liver-friendly substances, among which 19 ingredients have been demonstrated to benefit to the liver by animal or cell experiments.For the other four ingredients, although direct evidence focused on their liverprotecting effect was not found, no one of them was implicated by drug-induced liver injury.In addition, RSs analysis found that there existed abundant hepatoprotection-related RSs among the potential hepatoprotective ingredients identified.Data mentioned above indicated that our computational pharmacology method was highly reliable.We believed that our method would provide a strong technical support to disclose the material basis of TCM-induced hepatoprotection.Then based on the hepatoprotective ingredients Molecular docking score of the protein-ligand complexes.Compound 20: 11alpha, 12alpha-epoxy-3beta-23-dihydroxy-30-norolean-20-en-28,12beta-olide.
discovered above, we investigated the material basis and mechanism of BS against NAFLD comprehensively.
Consistency analysis between the NAFLD-related targets and BS-related targets identified 78 common targets.Both the tissue-specific expression pattern analysis and subcellular localization analysis indicated that these common targets were significantly correlated with NAFLD.Molecular network analysis found that at least 12 ingredients in BS could regulate these common targets.Docking modes of the protein-ligand complexes which were selected to perform MD simulation.Figures (A, C, and E) separately represented the docking complexes with the lowest binding energy for IL6, IL1β, and JUN.The docking complexes containing paeoniflorin were display in figures (B, D, and F), respectively.
Flavonoids, triterpenes, monoterpene glucosides, and phenolic acids were found to be the major structure types.Our prior studies suggested that phenolic acids, flavonoids, and triterpenes were important natural sources of liver protectants (He et al., 2022).In addition, monoterpene glucosides were reported to be the major active components of BS (Jiang et al., 2020).Therefore, we speculated that these 12 ingredients may be the major material basis of BS against NAFLD.Interestingly, it seemed that molecular docking analysis also verified such a speculation.Compared with the other structure types, flavonoids, triterpenes, monoterpene glucosides exhibited stronger binding affinity with the key targets of BS on NAFLD.
Obviously, molecular docking analysis further indicating the importance of these ingredients in the process of BS against NAFLD.In fact, oxypeucedanin, eugenitin, α-cedrene and αcedrol also showed stronger binding affinity with the key targets.We speculated that the lack of target data may be an important reason for these ingredients escaped from the identification of molecular network analysis.After all, for these four ingredients, only two targets were collected.
Focused on the PPI network of BS on NAFLD, modular analysis identified three functional modules in total.By mapping the targets included in each functional module into the DAVID platform, a total of 23 NAFLD-related KEGG pathways were enriched, among Frontiers in Pharmacology frontiersin.org18 He et al. 10.3389/fphar.2024.1406188which not less than 11 pathways have been reported to be associated with BS-induced hepatoprotection.Besides, at least five pathways have been demonstrated to play essential roles in the process of BS alleviates NAFLD.Further analysis found that a total of 12 NAFLDrelated biological processes were involved by the KEGG pathways enriched.It's worth noting that not less than eight out of these biological processes have been proved to be involved in the process of BS against NAFLD.Except for identifying the key pathways and the major biological processes, we also identified the key targets based on hub gene analysis and microarray data analysis.As a result, a total of three genes, including IL1-β, IL6, and JUN, were identified.Through retrieving previously published literature, we found that these three key genes indeed can be significantly downregulated by BS.In summary, our findings were highly consistent with reports in the literature, indicating that the bioinformatics strategy adopted in the current study were reasonable and reliable.However, we must acknowledge that there still existed some limitations in our research.Although we have validated our results through molecular docking analysis, molecular dynamics simulation analysis and literature analysis, but more in-depth experimental verification is required to further confirm our findings.In addition, the incomplete and missing of BSrelated targets may also affect our results.After all, some potential hepatoprotective ingredients were found to lack of the target data when we collected BS-related targets.With the constant perfection of relevant data in the future, we believed that our bioinformatics strategy would achieve a more satisfactory performance.

Collection of candidate compounds of BS
The chemical components of BS were extracted from three typical TCM databases, including TCMSP (Ru et al., 2014), TCMID (Xue et al., 2013), and, ETCM (Xu et al., 2019).The duplicates and the compounds without structures were removed.For the remaining components, the pharmacokinetics parameters, including drug likeness (DL) and oral bioavailability (OB), were investigated.Those components satisfy at least one of the following two conditions were defined as candidate compounds: (1) OB ≥ 30%, DL ≥ 0.18 (data from TCMSP database), (2) both QED (quantitative estimate of DL) and OB are good or moderate classes (data from, ETCM database).
4.2 A computational pharmacology method to identify the material basis of TCMinduced hepatoprotection In this section, a computational pharmacology method was proposed to screen hepatoprotective ingredients from TCMs, by which we attempted to provide valuable clues for elucidating the material basis of TCM-induced hepatoprotection.The computational pharmacology method consisted of three modules detailed as follows.
Module 1: Module one is a hepatoprotective "TCM-ingredient" network which consisted of 638 nodes and 2,262 edges.A total of 433 hepatoprotective ingredients and 205 hepatoprotective TCMs were involved.This network intuitively displayed which compounds were associated with the liver protection of the TCMs existed in the network.
Module 2: Module two is a large scale dataset of TCM-induced hepatoprotection, including 677 hepatoprotective phytoconstituents.For each phytoconstituent, we provided the English name and the simplified molecular input line entry system (SMILES) information, which enables researchers to retrieve this dataset by the method of name or structure matching.
Module 3: Module three is an in silico model which aimed at predicting the hepatoprotective activity of phytoconstituents derived from TCMs.This in silico model was established based on 709 phytoconstituents by integrating seven types of machine learning algorithms.Both the 5-fold crossvalidation and the external validation produced accuracy greater than 85%, indicating that this model is reasonably successful.
Details on these three modules mentioned above can be found in our prior published work (He et al., 2022).It's worth noting that the hepatoprotective ingredients identified by the computational pharmacology method can be classified into two categories, including the known hepatoprotective ingredients and the putative hepatoprotective ingredients.The former was provided by module one and module 2, and the latter was provided by module 3.They constitute together the material basis of TCM-induced hepatoprotection.The detailed workflow of the computational pharmacology method was provided in Figure 15.

Collection of targets for the potential hepatoprotective ingredients in BS
For each ingredient, the targets consisted of the known targets and the putative targets.We obtained the known targets by searching DrugBank (Wishart et al., 2008) and CTD (Wishart et al., 2008) databases using the ingredient name as input.Generally, data from DrugBank and CTD databases was highly reliable, because it always has been well documented by experiments.The putative targets were extracted from two frequently used database platforms, STITCH (Szklarczyk et al., 2016) and Swiss Target Prediction (Gfeller et al., 2014).The reliability of the putative targets is often evaluated by the index of confidence (range from 0 to 1).The closer the confidence score is to 1, the more reliable the putative target is.However, if the threshold of the confidence is set too high, it is difficult to obtain sufficient data.Therefore, here, to make a better compromise between the quality and quantity of the data, we set the threshold values of data from STITCH and Swiss Target Prediction as 0.70 and 0.85, respectively.After removing the duplicates and those non-human targets, the remaining targets were further analyzed.

Collection of NAFLD-related targets
The NAFLD-related targets were collected from three existing resources, including DisGeNET (Piñero et al., 2017), CTD (Davis et al., 2009), and NCBI-gene (https://www.ncbi.nlm.nih.gov/gene/).We searched these databases by the keyword of "Non-alcoholic Fatty Liver Disease," and only those human-related and non-repeated targets were taken into consideration.In addition, to improve the reliability of the data from CTD, only those targets curated by the experts were extracted.

Construction of the "ingredienttarget" network
Firstly, consistency analysis between the NAFLD-related targets and BS-related targets was conducted based on the online analysis tool Venn Diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/), by which we aimed at identifying the common targets between NAFLD and BS.Then network visualization tool Gephi (version 0.9.2) was utilized to construct the "ingredient-target" network by taking the common targets and the common targetsrelated compounds as input.

Tissue-specific expression pattern analysis
The tissue-specific pattern of mRNA expression can indicate important clues about gene function.BioGPS (Wu et al., 2016) is a centralized gene-annotation portal.One of its most important function is that enables researchers to access the mRNA expression data of genes in 84 organs.Herein, to investigate the distribution of the common targets attained in Section 4.5 within the liver, we ranked the target expression patterns based on their expression levels in the liver.To be specific, for certain target, if it expressed highest in the liver among the 84 organs, it was ranked 1. Inversely, if it expressed lowest in the liver, it was ranked 84.

Subcellular localization analysis
As we all known that protein activities are tightly linked to the cellular compartment and microenvironment where they are present.Elucidating the subcellular localization of proteins contributes to deeply understand their biological functions.Here, The Human Protein Atlas (Pontén et al., 2008), a large-scale biological database aiming at mapping the entire human proteome, was utilized to perform the subcellular localization analysis.Workflow of the computational pharmacology method for identifying the material basis of TCM-induced hepatoprotection.HISG, hepatoprotective ingredient subgroup; HIG, hepatoprotective ingredient group.

Construction of the protein-protein interaction (PPI) network
Protein-protein association network is of great value to understand the biological phenomena.As a free and publicly available source, STRING (Szklarczyk et al., 2023) database providing a comprehensive and objective global network where both the direct (physical) and indirect (functional) interactions between proteins were included.In the current study, the protein-protein interaction (PPI) information was extracted from this global network by using the common targets between NAFLD and BS as input.The PPI information with confidence score <0.4 was removed, and the remaining data was input into the network visualization tool Gephi (version 0.9.2) to obtain the PPI network map.

Identification of the hub genes
The importance of nodes in a network was often inferred based on their network topological features.As a powerful Cytoscape plugin, CytoHubba (Chin et al., 2014) provides at least 11 algorithms to rank nodes in a network, among which the Degree algorithm is most commonly used.Similar to the Degree algorithm, the Maximal Clique Centrality (MCC) algorithm was also proposed based on the local network topological features of nodes.Chin et al. have attempted to identify the essential proteins from the yeast PPI network by implementing multiple algorithms.As a consequence, the MCC algorithm was found to exhibit a better performance on the precision in comparison with other algorithms (Chin et al., 2014).To integrate the advantages of different algorithms, both Degree and MCC algorithms were adopted.Besides, Edge Percolated Component (EPC) algorithm was also implemented.Different from the Degree and MCC algorithms, EPC algorithm was a global network topological features-based algorithm.For each algorithm, the top 10 genes were collected, and the common genes of those three algorithms mentioned above were defined as hub genes.

Microarray data analysis
Firstly, hepatic gene expression data (GSE89632) in patients with NAFLD and in healthy donors were downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).Chip GSE89632 was from platform GPL14951 and contained 24 healthy controls (HC), 20 simple steatosis samples (SS), and 19 nonalcoholic steatohepatitis (NASH) samples.Then GEO2R analysis (https://www.ncbi.nlm.nih.gov/geo/geo2r/) was performed to identify genes that are differentially expressed between the NAFLD patients and the healthy controls.The detailed screening criteria were set as p < 0.05 and | log(FC)| > 1.Finally, we compared the differentially expressed genes with the hub genes obtained in Section 4.9, and the common genes between these two datasets were further analyzed.

Molecular docking analysis
To investigate the interactions between the active compounds and the key targets, molecular docking analysis was performed based on AutoDock Vina (version 1.1.2).Here, the 2D structures of the active compounds were obtained from PubChem database, and the X-ray crystallography-based structures of the key targets were downloaded from RCSB Protein Data Bank.When we implemented the AutoDock program, the size of the Grid Box was set to 40 × 40 × 40, and the exhaustiveness value was set as 100.Finally, the docking complexes with the lowest binding energy were selected for further analysis.

Molecular dynamics (MD) simulations
MD simulation is an efficient method to evaluate the stability of the protein-ligand complexes.In present study, MD simulation was carried out based on AMBER 18 software package with the force fields of FF14SB and GAFF.For each studied docked complex, a 100 ns MD simulation was performed under the condition of constant temperature.Finally, several important indices, including Root Mean Square Deviation (RMSD), Mean Square Fluctuation (RMSF) and binding free energy were calculated, by which we attempted to evaluate the stability of the docked complexes comprehensively.Of note, here, the binding free energy was calculated based on the last 50ns trajectories by utilizing the MMPBSA.py module.

Modular analysis of the PPI network
Molecular network is a biological system which consisted of a large number of nodes and edges.The characteristic of high complexity making it to be a challenging work to elucidate its scientific connotation.Network module theory holds the view that biological system network has modularity, and the nodes in the same module always perform some biological functions cooperatively (Lorenz et al., 2011).Therefore, modular analysis is considered to be an effective method to analyze the complex biological network.During the past decade, dozens of modular analysis methods have been proposed.To solve the problem with optimization of the network module division, Gu et al. proposed a network structure entropy-based method to evaluate the effect of 11 commonly used module division methods.As a result, Molecular Complex Detection (MCODE) algorithm was found to be superior to other algorithms (Gu et al., 2018).Therefore, in present study, the MCODE algorithm was implemented to identify the functional modules of BS on NAFLD.

Kyoto Encyclopedia of genes and genomes (KEGG) pathway enrichment analysis
To identify the pathways involved in BS on NAFLD, the platform of Database for Annotation, Visualization and Integrated Discovery (DAVID) (Dennis et al., 2003) was utilized to perform Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, and only those pathways of which adjusted p-values <0.05 were taken into consideration.In addition, for each pathway attained, we annotated its biological functions systematically by summarizing the related information recorded in literature and databases.
All of the databases utilized in the current study were listed in Table 5 in detail, which would assist the reader to understand our research methods more easily and quickly.

Conclusion
In present study, to clarify the material basis and mechanism of BS against NAFLD in a precise manner, an original computational pharmacology method for identifying the hepatoprotective ingredient group of TCMs was proposed.Then by incorporating with the techniques of network pharmacology, molecular docking, and molecular dynamics simulation, the proposed computational pharmacology method was utilized to reveal the scientific connotation of BS on NAFLD from multiple perspectives, including active ingredients, key targets, key pathways and the major biological processes involved.As a result, a total of 12 ingredients, mainly including monoterpene glucosides, flavonoids, triterpenes, and phenolic acids, were found to be associated with the anti-NAFLD effect of BS.Hub gene analysis and microarray data analysis indicated that IL1-β, IL6, and JUN were the key targets of BS on NAFLD.The findings mentioned above were then further validated via molecular docking analysis, molecular dynamics simulation analysis, and literature analysis.In addition, the key KEGG pathways and the major biological processes of BS on NAFLD were also identified.It's worth noting that the NAFLD pathway was significantly enriched.Further analysis found that there was a great deal of evidence available in the literature to support the regulatory effect of BS on NAFLD pathway.In addition, inflammation, apoptosis, lipid metabolism, and glucose metabolism were found to play critical roles in the process of BS alleviate NAFLD.In summary, a novel and effective bioinformatics strategy was proposed to uncover the material basis and mechanism of TCM in this work, based on which the anti-NAFLD effect of BS was systematically investigated from multiple perspectives.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers.Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

FIGURE 2
FIGURE 2Structural classification of the potential hepatoprotective ingredients in BS.

FIGURE 8
FIGURE 8 PPI network (A) and functional modules (B-D) of BS on NAFLD.

FIGURE 9 KEGG
FIGURE 9KEGG pathway enrichment analysis of the major functional modules of BS on NAFLD.Figures (A-C) displayed the KEGG pathways enriched by functional modules 1, 2, and 3, respectively.In Figure (F), F1-F12 were the major biological processes which were involved in the occurrence and development of NAFLD. Figure (D) displayed the NAFLD-related biological processes enriched by functional module 1.The NAFLD-related biological processes enriched by functional module two and three were showed in Figure (E).

FIGURE 11
FIGURE 11 Identification of the key targets of BS on NAFLD.Figure (A) Venn diagram analysis.Figures (B-D) showed the top 10 genes in the PPI network identified by Degree, Stress, and EPC algorithms, respectively.The node size was proportional to the gene score.Higher score indicated a higher ranking.Figures(E, F) displayed the differentially expressed genes related to NAFLD.HC, healthy controls; SS, simple steatosis; NASH, non-alcoholic steatohepatitis.

FIGURE 14 MD
FIGURE 14 MD simulation.RMSD (A-C) and RMSF (D-F) profiles of the selected docking complexes.

TABLE 1
Potential hepatoprotective ingredients in BS.

TABLE 2
Molecular mechanisms of BS extracts-induced hepatoprotection.

TABLE 3
Molecular docking parameters.

TABLE 4
Binding free energies of the selected docking complexes.All of the energies were provided in the format of average ± standard deviation.

TABLE 5 A
summary of the databases and analysis tools utilized in the current study.