Exploration of the mechanism of Ziyin Tongluo Formula on preventing and treating postmenopausal osteoporosis based on the network pharmacology and molecular docking

Background: Postmenopausal osteoporosis (PMOP) is a global chronic and metabolic bone disease, which poses huge challenges to individuals and society. Ziyin Tongluo Formula (ZYTLF) has been proved effective in the treatment of PMOP. However, the material basis and mechanism of ZYLTF against PMOP have not been thoroughly elucidated. Methods: Online databases were used to identify the active ingredients of ZYTLF and corresponding putative targets. Genes associated with PMOP were mined, and then mapped with the putative targets to obtain overlapping genes. Multiple networks were constructed and analyzed, from which the key genes were selected. The key genes were imported to the DAVID database to performs GO and KEGG pathway enrichment analysis. Finally, AutoDock Tools and other software were used for molecular docking of core compounds and key proteins. Results: Ninety-two active compounds of ZYTLF corresponded to 243 targets, with 129 target genes interacting with PMOP, and 50 key genes were selected. Network analysis showed the top 5 active ingredients including quercetin, kaempferol, luteolin, scutellarein, and formononetin., and the top 50 key genes such as VEGFA, MAPK8, AKT1, TNF, ESR1. Enrichment analysis uncovered two signicant types of KEGG pathways in PMOP, hormone-related signaling pathways (estrogen , prolactin, and thyroid hormone signaling pathway) and inammation-related pathways (TNF, PI3K-Akt, and MAPK signaling pathway). Moreover, molecular docking analysis veried that the main active compounds were tightly bound to the core proteins, further conrming the anti-PMOP effects. Conclusions: Based on network pharmacology and molecular docking technology, this study initially revealed the mechanisms of ZYTLF on PMOP, which involves multiple targets and multiple pathways.


Background
Postmenopausal osteoporosis (PMOP) refers to a global chronic metabolic bone disease with increased fragility and fracture susceptibility occurring in postmenopausal women, who are confronted with the drop of ovarian function and estrogen levels. It is characterized by greater bone resorption than bone formation, low bone mass and bone microstructure destruction. Osteoporotic fractures are one of the main causes of post-menopausal women's disability and death. Therefore, PMOP has now become a major issue of life quality worldwide, and its prevention research is one of the hot spots in modern medical research [1][2]. The current drug treatment on PMOP can be divided into bone resorption inhibitors (estrogen, bisphosphonates, calcitonin), bone formation promoters (parathyroid hormone, uoride), and bone mineralization (calcium, active Vitamin D). Estrogen replacement therapy has achieved good results in the clinic, but its cost is expensive, and the side effect such as endometrial hyperplasia, increased breast cancer risk [3], gastrointestinal reactions [4] are still one of the clinical problem.
Based on network pharmacology and molecular docking analysis, this study tries to screen the main active ingredients of ZYTLF in the prevention and treatment of PMOP, its targets and signaling pathways, and provides directions for elucidating the material basis and mechanism of action of ZYTLF in the prevention and treatment of PMOP. (Fig 1).

Identi cation of ZYTLF active ingredients
A total of 14 traditional Chinese herbs were imported into the Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP http://tcmspw.com/tcmsp.php) [8], BAT-MAN TCM (http://bionet.ncpsb.org/batman-tcm/) [9] and Traditional Chinese Medicine Information Database(TCMID, http://www.megabionet.org/tcmid/) [10] to obtain corresponding certain or potential compounds. The ADME (absorption, distribution, metabolism and excretion) properties of a compound is of importance in drug discovery and development [11]. Therefore, oral bioavailability (OB) and druglikeness (DL) , the two most important indicators of bioinformatics to evaluate the characteristics of ADME, were used as standards for screening active compounds. The TCMSP database utilizes Obioavail 1.1 to calculate the OB value [12], and applies the Tanimoto coe cient, see formula 1 in the supplementary les section.
to calculate the DL value. Compounds that meet the requirement OB ≥ 30% and DL ≥ 0.18 were regarded as biologically active ingredients. Because Radix Glycyrrhizae Preparata is less important in this formula, the standard was set as OB ≥ 60% and DL ≥ 0.36 to avoid greater bias. Based on the above conditions, the active ingredients of ZYTLF can be screened by TCMSP for subsequent analysis.

Prediction of ZYTLF targets
The next step after screening active ingredients was to predict the targets that might stimulate biological effects. The TCMSP analysis platform is loaded with prediction function of large-scale molecular network targets based on chemistry, genomics, pharmacology and system analysis technology, which is currently a commonly used target prediction platform [12]. After inputting the active ingredient into the TCMSP analysis platform, the corresponding target could be collected. Uniprot (https: //www.Uniprot.org) [13] is a database containing accurate annotations of proteins and other substances, where the genetic information of the target can be obtained. The target information were inputted into the Uniprot database to acquire the standard gene name.
Construction of a network model of "herb-active ingredient-overlapping gene" The overlapping target genes and their corresponding active ingredients and herbs were inputted into Cytoscape 3.7.1 software [17], to build a "herb-active ingredient-overlapping gene" network model. CytoNCA plug-in [18] was used to conduct network topology analysis. According to the Degree Centrality (DC) and Between Centrality (BC), the key nodes in the network were screened. The higher the node's degree value was, the more important it was in the network.

Construction of PPI network of overlapping genes
The STRING database (http://string-db.org/, ver.11.0) [19] was utilized to constructed protein-protein interaction (PPI) network of overlapping proteins, setting the condition of data analysis mode "Multiple proteins", the type "Homo sapiens" (human), and the minimum mutual threshold "high con dence (0.700)". The remaining parameters remained unchanged. The information of PPI network was obtained and then imported into Cytoscape 3.7.1 software. The MCC algorithm in the CytoHubba plug-in was applied to perform network analysis to screen out core genes.

Enrichment analysis of core genes
The Database Visualization and Integrated Discovery system (DAVID, https://david.ncifcrf.gov/) is an online biological information database, integrating biological data and analysis tools to provide annotation information of comprehensive biological function. Quantitative screening technology has its advantages on the biological function analysis of genes [20][21]. The species was limited to "Homo sapiens", and enrichment analysis of Gene Ontology (GO, http: //www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG, http: //www .genome.jp/kegg/) was conducted. R software was used to draw advanced bubble diagrams and gene-pathway mechanism diagrams was obtained from the KEGG database.

Molecular docking simulation
The computer simulation docking technology was utilized to verify the reliability of this study through evaluating the binding e ciency of the overlapping proteins and main active ingredients in ZYTLF. To obtain the SDF structure of top 10 core compounds in the network from the PubChem database (https://pubchem. ncbi.nlm.nih.gov/). The 2D structure was processed and transformed into PDB format through PyMOL, and they were saved in PDBQT format as docking ligands. Simultaneously, to collect all the crystal structures of the key proteins from the RCSB protein data bank (PDB, http://www.pdb.org/) and select the ones with relatively higher resolution and unique ligands. AutoDock Tools software was used to remove the water molecules, isolate proteins and save it as receptors. The receptors and ligands were processed with PyMOL software and Auto Dock software and then docked through Vina.

Results
Active ingredients and corresponding targets of ZYTLF Through TSMSP platform, BAT-MAN TCM and TCMID database, a total of 1,738 compounds were obtained, and 92 compounds that met the standard were regarded as active ingredients. The detailed information of active ingredients of ZYTLF are displayed in Additional le 1: Table S1. There were 242 corresponding target genes for these active ingredients.

PMOP-related genes and overlapping genes
From Genecard, Disgenet and OMIM databases, a total of 1113 PMOP related genes after removing duplication were obtained. Compared target genes of ZYTLF with PMOP-related genes, a total of 129 overlapping target genes were acquired. (Fig. 2).The detailed information of ZYTLF-genes, PMOP-genes, and ZYTLF-PMOP overlapping target genes are shown in Additional le 2: Table S2.

Network analysis
A network of "herb-active ingredient-overlapping gene" was constructed. (Fig. 3). There were 231 nodes and 1,072 edges in the network, and the detailed information of the network are shown in Additional le 3: Table S3. The results of network topology analysis showed that the active ingredients such as quercetin, kaempferol, luteolin, scutellarein, and formononetin had higher degrees, which played an important role in the network. The network revealed that one active compound could correspond with multiple genes, and one gene could also correspond with multiple active compounds. The complex interaction between active components and intersecting genes manifested the multi-component and multi-target features of ZYTLF.
The intersecting genes were inputted into the STRING database to built a PPI network after setting the relevant parameters. The network displayed the relevant information of the genes and the interaction relationship. (Fig. 4). Then Cytoscape 3.7.1 software was used to visually process and analyze the PPI network information. Based on the MCC algorithm of the CytoHubba plug-in, the top 50 key genes were selected. (Fig. 5).
GO and KEGG enrichment analysis DAVID database was used to perform GO enrichment and KEGG pathway enrichment analysis on 50 core genes and extract signi cant enrichment results (FDR < 0.05). A total of 62 GO-BP terms, 6 GO-CC terms, 11 GO-MF terms and 59 terms on the KEEG pathway were obtained. The enrichment results showed that the biological processes involved mainly included positive transcriptional regulation signals of RNA polymerase II promoter, negative regulation of apoptosis, positive regulation of transcription using DNA as a template, aging, positive regulation of gene expression, etc.. Main cell components included nucleus, cytoplasm, extracellular space, etc. Regulated molecular functions mainly included enzyme binding, transcription factor binding, and protein binding. (Fig. 6 a-c).
A total of 15 signal pathways related to PMOP were screened from the signi cant signal pathways. (Fig.  6 d). The signal pathways mainly involved in ammation, including TNF, HIF-1, PI3K-Akt, Toll-like receptors, MAPK and other signaling pathways, and hormonal pathways such as estrogen, prolactin, and thyroid hormone signaling pathways, indicating ZYTLF could play a synergistic role through many pathways. From the KEGG database, the mechanism of action of key genes on estrogen and TNF signaling pathways was obtained. The red marks in the gure represent potential targets for possible intervention by ZYTLF. (Fig. 7 and Fig. 8).The results and detailed information on the GO terms and pathways are shown in Additional le 4: Table S4.
The core compounds in ZYTLF are docked with key proteins Molecular docking simulation was applied to verify if the top 10 compounds have an pivotal role in regulation of the 50 key proteins. In molecular docking, if the ligand can bind to one or more amino acid residues in the active site (also called active pocket) of the receptor and participates in the process of conformation change, energy complementation, and so on, then the small molecule ligand can combine with the receptor to form a stable structure. This study indicated that the top 10 compounds had strong a nity with 23 core proteins,including ESR1 (PDB id: 4pxm), TNF (PDB id: 3m2w), AKT1 (PDB id: 6npz), PTGS2 (PDB id: 5f19), which have an important role according to the result of network pharmacology.
The quercetin can bind closely with the active pocket of ESR1 and TNF proteins. Quercetin developed four hydrogen bonds with the amino acid residues LEU-72, LYS-93, EU-141, VAL-78, making quercetin and TNF form a stable complex (Fig. 9). Meanwhile, Quercetin developed three hydrogen bonds with the amino acid residues GLY-521, GLU-353, ARG-394, making quercetin and ESR1 form a stable complex (Fig. 10). The binding energies of the various compounds are shown in Additional le 5:Table S5.

Discussion
At the molecular level, multiple factors participate in the development of PMOP, interacting with each other and regulating bone metabolism in a complex regulatory network. Estrogen de ciency is the pivotal factor in the development of PMOP, which causes to the disorder of bone metabolism. Bone resorption is greater than bone formation, and eventually osteoporosis is caused [22]. Compared with hormone replacement therapy, traditional Chinese medicine, as a welcome alternative therapy, has greater advantages possessing many merits of simplicity and low side effect. Modern research have found that traditional Chinese herb contains various bioactive ingredients, which can act on multiple targets through multiple pathways, validating the holistic view of traditional Chinese medicine and the theory of treatment based on syndrome differentiation [23]. Compared with the "single component-single targetsingle pathway" research model of traditional pharmacology, network pharmacology combines the concepts of bioinformatics and multi-directional pharmacology to analyze the relationship between biological systems, medicines, and diseases from a network perspective. It has opened a new research model of traditional Chinese medicine from empirical medicine to evidence-based medicine [24].

Key active ingredient of ZYTLF for anti-PMOP
According to the network analysis results of the "herb-active ingredient-overlapping gene", the key active ingredients of ZYTLF are quercetin, kaempferol, luteolin, rhubarbin, formononetin and other avonoid compounds. At the same time, they are also phytoestrogens, which has two-way regulating effects. It exerts the role of estrogen on target organ when estrogen de ciency occurs in body, and plays the role of anti-estrogen when the level of estrogen increases. It can alleviate a series of clinical symptoms caused by the postmenopausal estrogen de ciency and avoid the occurrence of tumors in the reproductive system caused by estrogen replacement therapy [25]. Quercetin with the highest degree value in the network is a typical avonoid compound. It has multiple pharmacological effects, including scavenging free radical, anti-cancer, anti-infection, and protecting cardiovascular [26][27][28][29]. Quercetin has also been proven to be an effective component against osteoporosis, with the dual effects of inhibiting osteoclastogenesis and osteoblast differentiation. Zhixing Li found that quercetin can relieve osteoporosis symptoms in ovariectomized rats, possibly by up-regulating ALP gene expression and inhibiting JNK, ERK, and p38 MAPK signaling pathways [30]. Quercetin can also directly or indirectly down-regulate the expression of RANKL, inhibit osteoclast differentiation, reduce bone resorption, and stimulate osteoblast activity through estrogen signaling pathway and ERK signaling pathway [31]. Kaempferol has also been shown to have bone protective effects on ovariectomized rats [32], possibly through estrogen receptor, MAPK, NF-κB and other signaling pathways [33]. Luteolin can prevent bone loss after osteoporosis by inhibiting osteoclast differentiation.

Key Genes of ZYTLF for Preventing PMOP
Combined with the analysis of PPI network and key gene network, TNF, JUN, MAPK8, AKT1, IL-6, MMP9, PTGS2, MAPK1, CASP3 were of great importance and all related to the in ammatory response process. Actually, increasing studies have found that the reduction of estrogen levels in postmenopausal women can stimulate the immune system to produce a large number of osteoclastogenic factors, which in turn activates related signaling pathways, further aggravating bone loss [33]. Zha Li and other scholars have found that postmenopausal women with osteoporosis have signi cantly increased levels of TNF-α. Invitro experiments shows that TNF-α and RANKL synergistically enhance bone resorption of osteoclasts through NF-κB and PI3K/Akt signaling pathways [35]. TNF-α and IL6 play an important role in the immune response and bone metabolism, mainly affecting the differentiation and proliferation of osteoclasts by regulating complex mechanisms, and they are important pathogenic factors for immune-mediated bone diseases [36]. Studies have shown that increased levels of in ammatory factor TNF-α in the serum of postmenopausal women may be one of the important causes of osteoporosis. TNF-α can activate the RANK/RANKL signaling pathway and induce the formation of osteoclasts [37]. We speculates that elevated TNF-α level in the serum of postmenopausal women with osteoporosis may be related to estrogen de ciency.

The signal pathways ZYTLF in preventing PMOP
The enrichment analysis of KEGG pathways for the key genes of ZYTLF in preventing PMOP showed that 15 signaling pathways were related to the development and progression of PMOP, including hormone pathways and in ammation-related pathways. Hormones in postmenopausal women, including estrogen [38], parathyroid hormone [39], and prolactin, if not normally secreted could affect bone metabolism [40]. Estrogen bound with the estrogen receptor in osteoblasts and osteoclasts to act on the OPG/RANK/RANKL signaling pathway, which further promoted OPG secretion, down-regulated the expression of RANKL, and inhibited the formation of osteoclasts [41]. After the combination of estrogen and estrogen receptor, it could also regulate the expression of various target genes through the estrogen signaling pathway, thereby activating downstream PI3K/Akt, MAPK, WNT and other signaling pathways to promote the differentiation and proliferation of osteoblasts [42][43]. The TNF signaling pathway had the largest number of enriched genes, which might be the key signaling pathway for ZYTLF to prevent PMOP. The TNF signaling pathway was mainly opened by TNF-α and interacted with multiple signaling pathways to synergistically inhibit osteoclast differentiation and bone resorption function [44].

Veri cation of molecular docking
In this study, network pharmacology was applied to elaborate on the potential anti-PMOP effect of ZYTLF and this effect was visualized by molecular docking for veri cation. In the result, ESR, TNF and other 21 proteins can combined with the top 10 active compounds, which have a signi cant role in regulation of the 23 proteins associated closely with PMOP. The result proved that ZYTLF could exert the treatment of PMOP by activating those receptors in the estrogen signaling pathway and TNF signaling pathway. It is consistent with the result of a vitro study showing that quercetin protected against TNF-α induced impairments in bone marrow mesenchymal stem cells [45]. In another study, quercetin can exert dual effect on estrogen, displaying anti-estrogenic effect at low dose (10mg/kg/day), and exaggerated estrogenic activities at high dose (100mg/kg/day) [46]. However, there are several limitations in this study. First, the compatibility and dosage of the formula are not taken into consideration. Second, the screened active components are different from the actual components in blood. Third, even if the results of network pharmacology and molecular docking simulation were combined, we still could not completely understand the accurate therapeutic mechanism of ZYTLF. Therefore, further experiments need to be conducted combining with Western Blot, Rt-PCR and multi-omics technology, etc. However, the network pharmacology methodology can save our time and energy by predicting the active ingredients, targets and signal pathways, and verifying it with molecular docking simulation.

Conclusions
Quercetin, kaempferol, luteolin and other core compounds of ZYTLF may be enriched in estrogen  signaling pathway and TNF signaling pathway through ESR1, TNF, IL6, MAPK8, STAT3, MMP9, PTGS2,  JUN, inhibiting osteoclast formation, promoting    There are 242 active compounds of ZYTLF and 1113 genes associated with PMOP. 129 overlapping genes were selected as potential targets for further study.   Key genes in overlapping target genes The color of the node changed from light yellow to dark red, indicating that the higher the MCC value was, the more signi cant the role it played in the network.  GO analysis and KEGG pathway enrichment analysis of Key gene In the bubble diagrams above (a-d), the ordinate represents the names of BP, CC, MF terms and pathways, respectively, and the abscissa represents the degree of enrichment. The smaller the FDR is, the higher the importance of enrichment, and the redder the color on the diagram. Speci c mechanism of key gene enrichment in estrogen signaling pathway Figure 8 Speci c mechanism of key gene enrichment in TNF signaling pathway