Exploring the potential therapeutic effect of Eucommia ulmoides–Dipsaci Radix herbal pair on osteoporosis based on network pharmacology and molecular docking technology

Eucommia ulmoides–Dipsaci Radix (EU–DR) is a commonly used herbal pair for the treatment of osteoporosis (OP) in China. The purpose of this study was to investigate the potential mechanism of EU–DR on OP through network pharmacology and molecular docking approaches. Combining data from multiple open-source databases and literature mining, the active compounds and potential targets of EU–DR were screened out. The OP related targets were identified from the interactive web tool GEO2R. The shared targets were obtained by intersecting the targets of EU–DR and OP. The protein–protein interaction (PPI) network was constructed via the STRING database and Cytoscape 3.7.2 software. GO and KEGG enrichment analysis were conducted using R 3.6.3 software with adjusted p-value < 0.05. Sybyl-x 2.1.1 and Autodock Vina 1.1.2 software were used to cross validate the affinity between active compounds and target proteins. Our results showed that a total of 50 active compounds were screened, corresponding to 895 EU–DR targets, 2202 OP targets and 144 shared targets. The flavonoids in EU–DR played an important role in anti-OP. The enrichment analysis of GO and KEGG suggested EU–DR exerted a therapeutic effect on OP mainly by regulating the osteoclast differentiation related signaling pathway. Meanwhile, molecular docking results showed that most active compounds in EU–DR had strong binding efficiency to the target proteins. In conclusion, this study elaborated the multi-component, multi-target, and multi-pathway interaction mechanism of the EU–DR herbal pair in the treatment of OP for the first time, which also provided a pharmacological basis for treating OP.


Introduction
Osteoporosis (OP) is a chronic metabolic bone disease characterized by bone loss, especially in elderly people and menopausal women. 1 It has become one of the most common bone diseases in the world and has caused a heavy economic burden on society. The causes of OP are very complex, including the interaction of endocrine, nutritional status, genetic, physical and immune factors. 2 The disorder of dynamic balance between bone formation of osteoblasts and bone resorption of osteoclasts is the main cause of its pathological changes. An increase in the number of osteoclasts may lead to excessive bone resorption, while impaired differentiation of osteoblasts may reduce the formation of new bone, disrupt the balance of bone reconstruction, and ultimately lead to bone loss. 3 Currently, the treatment strategies for OP mainly use drug therapy, including calcitonin, bisphosphonates and uorides. 4,5 These drugs can alleviate bone loss and improve the clinical symptoms to a certain extent, but their long-term clinical application is limited by low tolerance, severe side effects and high cost. 6 Meanwhile, these drugs cannot fundamentally improve bone metabolism and restore the dynamic balance of osteogenesis and osteoclastic activity. 7 Thus, it is very important to develop more effective and safer therapeutic strategies for OP.
Traditional Chinese Medicine (TCM) has been used to prevent and treat various diseases in a multi-target and multicomponent manner for thousands of years. 8 In the theoretical system of TCM, OP is dened as bone accidity or arthralgia syndrome caused by kidney deciency. 9 In the treatment of OP, TCM formulas can not only reduce bone loss, but also comprehensively regulate body functions and relieve back pain. Herbal pairs are common clinical TCM formulas, which generally consist of two kinds of herbs. Compared with TCM formulas containing more herbs, it is simpler and more bene-cial to elucidate the complex pharmacological mechanism of a herbal pair. 10 Eucommia ulmoides (EU, Chinese name: Duzhong) is the bark of Eucommia ulmoides Oliv., a typical Chinese herb for tonifying the kidney, while Dipsaci Radix (DR, Chinese name: Xuduan) is the dry root of Dipsacus asper Wall., another natural herb for strengthening the bone. 11,12 Previous studies have shown that EU and DR are currently the most widely used TCM for the treatment of OP. 13 EU contains various chemical components, including lignans, steroids, terpenoids and avonoids. EU extracts can inhibit bone loss and maintain metabolic balance. 14 A variety of chemical components have also been identied in DR, especially saponins, iridoid glycosides and sterols. The crude extracts of DR can increase bone density and change bone histomorphology. 15 Moreover, both of EU and DR are generally recognized as a classic herbal pair and play a synergistic effect in anti-OP. 13 However, the potential pharmacological mechanism of the EU-DR herbal pair and its interaction with OP related targets and pathways are still unclear, and further exploration is needed.
Modern pharmacological studies have demonstrated that TCM has the characteristics of complex mechanism and multipathway interaction in the treatment of diseases. 16 Based on the interaction network of multi-compound, multi-target and multi-pathway, network pharmacology can elucidate the mechanism of TCM through constructing a "drug-targetdisease" network. 17 This method conforms to the fundamental principles of TCM and has the characteristics of holistic and systematic. Molecular docking is a virtual screening technology that explores the interaction between receptors and drugs by simulating the behavior of small molecule ligands at the binding sites of receptor proteins. 18 In recent years, the application of molecular docking to explain the relevant mechanism has become a trend in new drug research and development. Therefore, in this study, network pharmacology and molecular docking were carried out to clarify the potential mechanism of the EU-DR herbal pair for the treatment of OP. The ow diagram of the research is shown in Fig. 1.

Identication of active compounds in EU-DR
In order to obtain comprehensive information, the compounds of EU and DR were derived from four databases, including Traditional Chinese Medicine Systems Pharmacology (TCMSP) database (http://tcmspw.com/tcmsp.php), Traditional Chinese Medicine Information Database (TCM-ID) (http:// www.megabionet.org/tcmid/), Encyclopedia of Traditional Chinese Medicine (ETCM) database (http://www.tcmip.cn/ ETCM/index.php/Home/Index/index.html) and SymMap database (https://www.symmap.org/). The active compounds in EU-DR were screened based on the ADME related model of candidate compounds retrieved from TCMSP and Swiss ADME databases (http://www.swissadme.ch/). 19 The screening criteria in the TCMSP database were compounds fullling both oral bioavailability (OB) value $ 30% and drug-likeness (DL) value $ 0.18, and the screening criteria in the Swiss ADME database were as follows: gatrointestinal absorption: high, Lipinski: Yes, Veber: Yes, Ghose: Yes, Muegge: Yes and Egan: Yes. 9 In addition, it was also found that some compounds with signicant biological activity were removed because their values were lower than the screening criteria. Considering their potential druggability, we used large-scale text mining to supplement compounds with notable pharmacological activity. Finally, we obtained the PubChem ID, 2D structure and canonical SMILES les of the active compounds from PubChem database (https:// pubchem.ncbi.nlm.nih.gov/).

The potential targets of active compounds in EU-DR
To identify the potential targets related to the active compounds, we collected information from three databases, including SymMap database, Swiss Target Prediction database (http://www.swisstargetprediction.ch/) and TargetNet database (http://targetnet.scbdd.com/home/index/). All chemical structures were prepared and converted into canonical SMILES using PubChem. The SMILES les of the compound were put into Swiss Target Prediction database and TargetNet database, then the result with probability greater than 0.9 was output. 20 The targets related to the active compounds were collected through the molecular names in SymMap. Data were further analysed using UniProt database (http://www.uniprot.org/) to conrm that these target proteins had been veried and set as "Homo sapiens".

OP related targets and shared targets of EU-DR against OP
The interactive web tool GEO2R was used to retrieve the OP related targets. GEO2R is an analytical tool provided by Gene Expression Omnibus (GEO) Dataset (https:// www.ncbi.nlm.nih.gov/gds), which aims to compare different sample groups in the GEO Series to identify differentially expressed genes (DEGs) between normal samples and model samples. GSE35958 is a classical gene chip of OP model and is frequently cited. 21,22 In order to identify the genes that were signicantly differentially expressed in each sample, the Benjamini-Hochberg false discovery rate (FDR) method was employed to adjust P-values in GEO2R analysis. The cutoff criteria for screening DEGs was FDR < 0.05 & jlog FCj > 1. 23 Moreover, according to the results of DEGs analysis in GEO2R, a volcano diagram was drawn using the ggplot2 program package in R 3.6.3 soware. Subsequently, a venn analysis was performed through the VennDiagram package in R 3.6.3 soware to overlap the OP related targets and the potential targets of EU-DR.

Protein-protein interaction (PPI) and hub target analysis
PPI analysis is helpful to identify the hub target related to EU-DR on OP. The shared targets of EU-DR relevant to OP were retrieved to construct the PPI network via STRING database (https://string-db.org/) with "Homo sapiens" setting. Cytoscape 3.7.2 soware was used to visualize the interconnection network, and "Analysis Network" function was performed to analyze the network properties. In the PPI network, the degree of freedom (DOF) is an important index, since the points above the average DOF usually play an important role in the network. 24

GO and KEGG pathway enrichment analysis
In order to explore the core functions and biological pathways related to EU-DR in treatment of OP, we used the ClusterProler package in R 3.6.3 soware to perform Gene ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. GO functions include the biological process (BP), cellular component (CC), and molecular function (MF) terms, which are used to analyze the potential mechanism. KEGG was applied to identify the biological pathways. The enrichment results with FDR less than 0.05 were selected for further analysis.

Network construction
To better characterize the "multi-compound, multi-target and multi-pathway" therapeutic features of EU-DR, Cytoscape 3.7.2 soware was used to construct a primary regulatory network to intuitively display the "drug-target-disease" relationship between EU-DR and OP. In this bilateral network, the nodes represent active compounds, shared targets and pathway, while the edges represent the interactions between compounds, targets and pathway. The DOF of a node represents the number of targets associated with the node in the network.

Molecular docking simulation
The Surex-Dock module in Sybyl-X (version 2.1.1, Tripos Inc.) soware was used to evaluate the binding ability of active compounds with target proteins. The crystal structures of candidate target proteins were downloaded directly from RCSB Protein Data Bank (PDB, http://www.rcsb.org/) and modied by Sybyl-X soware, including water removal, hydrogen addition,    side chain repair, amino acid optimization and active pocket construction. The 3D chemical structures of the active compounds were obtained from the PubChem database. Subsequently, the Surex-Dock docking mode was used for molecular docking and the total score of binding affinity was obtained. In general, total score >4.0 indicates certain binding activity, total score >5.0 indicates good binding activity, while total score >7.0 indicates strong binding activity (the higher the

Collection of active compounds in EU-DR
According to the screening criteria of ADME, a total of 34 EU-DR active compounds were retrieved from the above databases. During the ADME screening process, we noticed that some compounds, such as Aucubin, Geniposide, Geniposidic acid and Asperosaponin VI had high pharmacological activities, but they did not meet all the criteria. 26-28 Therefore, we have also retained them for subsequent analysis. Aer deleting repetitions, a total of 50 active compounds were screened, including 34 active compounds in EU, 20 active compounds in DR, and 4 shared compounds in EU and DR. The details are shown in the Table 1.

Screening of potential targets in EU-DR
In this study, 833 targets corresponding to EU and 467 targets corresponding to DR were identied. However, most of the targets were repetitive, revealing that the active compounds of EU and DR may have the same biological effects and achieve synergistic effects when combined. Finally, through further integration of all targets, 895 targets corresponding to EU-DR were obtained.

Related targets of OP and shared targets of EU-DR against OP
The gene chip GSE35958 contains 4 normal bone marrow samples and 5 OP bone marrow samples to analyze the effect of OP on the transcriptome of human mesenchymal stem cells (hMSCs) from human bone marrow. Following the above screening criteria of FDR < 0.05 & jlog FCj > 1, a total of 2202 DEGs between normal samples and model samples were iden-tied as shown in the volcano diagram ( Fig. 2A). As shown in the Venn diagram (Fig. 2B), there were 144 shared targets between 2202 DEGs and 895 potential targets to further clarify the mechanism of EU-DR in OP treatment.

PPI network construction
A total of 144 shared target genes were imported into STRING database for the construction of PPI network. In this study, the 0.4 condence (medium) was set as the default threshold. In order to obtain a better understanding, the results were further analyzed by Cytoscape 3.7.2 soware, and a PPI network with 144 nodes and 1103 edges was obtained (Fig. 3A). The hub target possessed a higher DOF and was more likely to play a critical role in the network of EU-DR on OP, such as RACalpha serine/threonine-protein kinase (AKT1), cellular tumor antigen p53 (TP53), mitogen-activated protein kinase 3 (MAPK3), mitogen-activated protein kinase 1 (MAPK1), epidermal growth factor receptor (EGFR) and signal transducer and activator of transcription 3 (STAT3).The average DOF of PPI network was 17, and the top 42 hub targets above the average DOF were shown in Fig. 3B.

GO functional annotation and KEGG pathway enrichment analysis
In GO functional annotation analysis, the 144 shared targets were highly enriched in 2115 BP terms, 49 CC terms, and 125 MF terms, with FDR < 0.05. The top 10 enrichment results of terms and the proportions of each part were shown in Fig. 4. Among them, BP terms account for the largest proportion, accounting for 91.39%, followed by MF terms and CC terms, accounting for 6.18% and 2.42% respectively. Peptidyl-serine phosphorylation, peptidyl-serine modication, response to antibiotic, peptidyl-tyrosine phosphorylation and peptidyltyrosine modication were closely related to BP. In regard to CC, higher enrichment was found in membrane ra, membrane microdomain, membrane region, RNA polymerase II transcription factor complex and transcription factor complex. The main terms related to OP in MF included protein  tyrosine kinase activity, protein serine/threonine kinase activity, phosphatase binding, nuclear receptor activity and transcription factor activity, direct ligand regulated sequence-specic DNA binding. The analysis results showed that these targets were closely related to the processes of regulating steroid hormone activity and promoting steroid hormone receptor activity and binding. The focus of KEGG enrichment analysis of 144 shared targets was the anti-OP pathway of EU-DR, with FDR < 0.05. The results showed that 145 pathways were enriched, and the top 5 pathways with the largest number of enriched targets were shown in Fig. 5. Among them, MAPK signaling pathway (hsa04010), proteoglycans in cancer (hsa05205), PI3K-Akt signaling pathway (hsa04151), human T-cell leukemia virus 1 infection (hsa05166) and osteoclast differentiation (hsa04380) were highly associated with the occurrence and development of OP. In addition, other KEGG pathways were also reported to be relevant to OP, such as estrogen signaling pathway (hsa04915), VEGF signaling pathway (hsa04370) and Wnt signaling pathway (hsa04310).
Moreover, gene targets that act on multiple pathways at the same time were found in KEGG enrichment analysis, indicating the importance of specic targets in the entire OP bioinformatics network. For instance, AKT1, MAPK1 and MAPK3 appeared in the 5 pathways with the largest number of enriched targets. These also reected the "multi-target and multipathway" regulation mechanism of EU-DR.

Network construction analysis
As shown in Fig. 6, a regulatory network of "drug-targetdisease" was constructed, showing the interactions among 50 compounds, 144 targets, 30 KEGG pathways and 30 GO terms. This network was consisted of 256 nodes and 1848 interaction edges, which indicated that the active compounds of EU-DR act on the whole biological network system through multi-target and multi-pathway. According to the DOF, we found that avonoids played a greater role in the regulatory network. Among these compounds, Quercetin (degree ¼ 48), Genistein (degree ¼ 42) and Kaempferol (degree ¼ 29) were the most potential compounds. At the same time, the hub targets predicted by PPI also possessed higher DOF in this network, such as MAPK3 (degree ¼ 49), MAPK1 (degree ¼ 49), EGFR (degree ¼ 48), AKT1 (degree ¼ 47), RELA (degree ¼ 39), indicating that they were involved in more signaling pathways.

Molecular docking verication
In molecular docking, the ligand binds to one or more amino acid residues with hydrogen bonds in the active pocket, and participates in the process of conformational changes and energy complementation. 29 The binding sites and binding score values can intuitively reect the interaction and stability of the docking model. Here, we select the rst 15 targets with higher DOF (except for receptors lacking the original ligand in the PDB) to construct protein receptors based on hub targets of PPI. According to degree of the "drug-target-disease" regulatory network, the rst 6 compounds with higher degree values were selected from EU-DR herbal pair for molecular docking verication, including quercetin, genistein, kaempferol, ursolic acid, dehydrodiconiferyl alcohol and caffeic acid. The docking results showed that most of docking scores were greater than 4, indicating that the binding between ligand and receptor was relatively stable (Fig. 7A). Besides, the cross validation between Sybyl-X and Autodock Vina shown AUC ¼ 90.2709, indicating the docking results is reliable (Fig. 7B).
The binding modes of 6 targets with their corresponding original ligands and the above 6 active compounds of EU-DR were shown in Fig. 8. It was observed that there was a high resemblance between active compounds and original ligands of protein, and they nearly occupied the same active pocket. The binding mode of quercetin at the active site of GSK3B had been represented in its three-dimensional mode in Fig. 8A. Four amino acid residues formed hydrogen bond interactions with quercetin as follows: the phenolic hydroxyl group on the benzene ring of quercetin linked with Val135, Asp133, Asp200 and Lys85 residue respectively. The binding mode of genintein at the active site of AKT1 was shown in Fig. 8B. The molecule was located in the binding pocket, surrounded by Asn279, Leu156, Thr291 and Asp292 residue. The binding mode between kaempferol and ESR1 was shown in Fig. 8C. The residues of Leu346, Leu387 and Arg394 form hydrogen bond interactions with the hydroxyl groups of kaempferol. Fig. 8D showed the binding mode between ursolic acid and TP53, including a carboxyl group linked with Arg202 residue and the hydroxyl group linked with Pro223 and Val225 residue. Dehydrodiconiferyl alcohol had a strong binding ability with the receptor protein MAPK1 by forming four hydrogen bond interactions. One was the phenolic hydroxyl group connected to Ile31 residue from the benzene ring, and the other three were the hydroxyl group connected to Asp167, Tyr36 and Ala35 residues (Fig. 8E). Caffeic acid and the receptor protein EGFR showed good binding ability. Meanwhile, amino acid residues with caffeic acid form hydrogen bond interactions to maintain conformational stability, including the carboxyl group connected to Lys745, Asp855 residue and phenolic hydroxyl group linked with Gln791 and Met793 residue (Fig. 8F).
Taken together, these docking results provided evidences that the active compounds of EU-DR could stably bind to the

Discussion
Previous studies have shown that EU and DR are currently the most widely used TCM in treating OP. 13 However, due to the complex chemical compositions of EU-DR herbal pair, it is still difficult to elucidate its exact pharmacological mechanism. In this study, we systematically revealed mechanism of EU-DR in the treatment of OP for the rst time by using network pharmacology method.
Based on this method, we found that most of the compounds of EU-DR affected multiple targets, and many repetitive targets were identied from different compounds, revealing that EU-DR exerted a therapeutic effect on OP through the synergistic effects of various active components. The results indicated that avonoids might play a greater role in anti-OP as they affected more targets relevant to OP. Quercetin was the most important active compound, which was involved in the regulation of 48 target genes. Following quercetin, genistein and kaempferol regulated 42 and 29 target genes, respectively. Flavonoids have a variety of benecial biological activities due to their antioxidant, anti-inammatory and estrogenic effects, such as quercetin and kaempferol are substances that can both reduce the activity of osteoclastic and increase the activity of osteoblastic. 30 Genistein is a natural isoavonoid that can stimulate the growth of osteoblasts. 31 Moreover, we found that AKT1, MAPK1 and MAPK3 were the core target proteins of EU-DR for the treatment of OP. They not only exhibited abundant interactions with other target proteins, but also participated in more signaling pathways. Maintenance of bone homeostasis is a dynamic process regulated by osteoblast-mediated bone formation and osteoclast-mediated bone resorption. AKT1 is considered to be a unique inuence factor that mediates osteoblasts and regulates the differentiation of osteoblasts and osteoclasts. 9,32 During osteoclast maturation, MAPK1 and MAPK3 regulate osteoclast maturation by up-regulating the expression of osteoclast-related genes. 33 The main compounds in EU-DR were identied to bind to these core target proteins, and further conrmed by molecular docking.
Additionally, according to the results of the KEGG pathway analysis, we found that EU-DR mainly involved in osteoclast differentiation related signaling pathways in treating OP (Fig. 9). Osteoclast differentiation is stimulated by macrophage colonystimulating factor (M-CSF) and the receptor activator of nuclear factor-kB ligand (RANKL). 34 M-CSF bind to its receptor c-Fms to induce the proliferation and survival of osteoclast precursor cells through activation of the ERK and PI3K/Akt pathways. The binding of RANKL to RANK recruits TNF receptor-associated factor 6 (TRAF6) to activate MAPK, PI3K/Akt, and NF-kB pathway, and promote the differentiation of osteoclast precursors into osteoclasts. TRAF6 is recruited to RANK to form a signaling complex containing TGF-b-activated kinase (TAK) 1 and TAK-1-binding protein (TAB) 2 to activate all three MAPK pathways, including extracellular signal-regulated kinase (ERK), c-Jun N-terminal kinase (JNK), and p38, and induce the activation of their downstream targets such as c-Fos, AP-1 transcription factors. 35 Kaempferol inhibits RANKL-mediated osteoclast differentiation through inhibiting RANKL-mediated phosphorylation of ERK, JNK and p38, and expressions of c-Fos and nuclear factor of activated T cells 1 (NFATc1). 36,37 Several studies illustrated that quercetin suppressed the phosphorylation of Akt and affected nuclear factor kappa B (NF-kB), AP-1 and NFATc1 to inhibit osteoclast differentiation. 38,39 Activation of NF-kB signaling pathway by RANKL, TNF-a or IL-1 can induce osteoclast differentiation to extend osteoclast survival and increase bone resorption. 40 Genistein directly suppressed osteoclastic differentiation induced by TNF-a and suppressed the expression of two transcription factors c-Fos and NFATc1 caused by NF-kB up regulation to inhibit bone resorption. 41 In addition, robust NFATc1 induction can be provided through activation of calcium signaling. 42 Genistein induce apoptosis of osteoclasts through calcium signaling. 43 According to literature and network pharmacological analysis, we believed that avonoids, the active ingredient of EU-DR herbal pair, could affect the differentiation and survival of osteoclasts, reduce bone resorption, and achieve the purpose of treating OP.
Drug therapy is the usual method in treating OP. Regardless of the benets, most patients treated with these drugs have severe side effects. 44 Therefore, a mild treatment strategy is very important for the treatment of OP. In addition to treatment, the prevention of OP is also very important. Regarding preventive care, drugs with less side effects are needed. In this respect, TCM such as EU and DR have great advantages. With more research on EU-DR herbal pair, its medicinal substances and pharmacological mechanisms will become more and more clear.

Conclusions
In conclusion, this study elaborated the multi-component, multi-target, and multi-pathway interaction mechanism of EU-DR herbal pair in the treatment of OP through network pharmacology and molecular docking method for the rst time. The results indicated that EU-DR herbal pair can regulate the osteoclast differentiation related signaling pathways to affect the proliferation and survival of osteoclasts. The avonoids in EU-DR played a greater role in anti-OP. The research provided a theoretical basis for the effective ingredients and mechanism of EU-DR herbal pair for treating OP.

Author contributions
Conceptualization, S. F. and T. W.; investigation, L. F.; methodology, M. W. and X. Z.; supervision, X. A. and X. D.; visualization, S. F. and Y. L.; writing-original dra, S. F.; writing-review and editing, Y. C., J. H. and Y. L. All authors have read and agreed to the published version of the manuscript.

Conflicts of interest
There are no conicts of interest.