Systems pharmacology dissection of action mechanisms for herbs in osteoporosis treatment.

Objective: Osteoporosis has become the biggest cause of non-fatal health issue. Currently, the limitations of traditional anti-osteoporosis drugs such as long-term ill-effects and drug resistance, have raised concerns toward complementary and alternative therapies, particularly herbal medicines and their natural active compounds. Thus, this study aimed to provide an integrative analysis of active chemicals, drug targets and interacting pathways of the herbs for osteoporosis treatment.
Methods: Here, we introduced a systematic pharmacology model, combining the absorption, distribution, metabolism, and excretion (ADME) screening model, drug targeting and network pharmacology, to probe into the therapeutic mechanisms of herbs in osteoporosis.
Results: We obtained 86 natural compounds with favorable pharmacokinetic profiles and their 58 targets from seven osteoporosis-related herbs. Network analysis revealed that they probably synergistically work through multiple mechanisms, such as suppressing inflammatory response, maintaining bone metabolism or improving organism immunity, to benefit patients with osteoporosis. Furthermore, experimental results showed that all the five compounds (calycosin, asperosaponin VI, hederagenin, betulinic acid and luteolin) enhanced osteoblast proliferation and differentiation in vitro, which corroborated the validity of this system pharmacology approach. Notably, gentisin and aureusidin among the identified compounds were first predicted to be associated with osteoporosis.
Conclusion: Herbs and their natural compounds, being characterized as the classical combination therapies, might be engaged in multiple mechanisms to coordinately improve the osteoporosis symptoms. This work may contribute to offer novel strategies and clues for the therapy and drug discovery of osteoporosis and other complex diseases.

Introduction

Osteoporosis is a common chronic skeleton disease contributed by polygenic and multiple environmental factors, which has been aroused considerable attention all over the world with the growing number of osteoporosis patients, especially in postmenopausal women and the aged (Bilezikian et al., 2018). Given the laid serious health problem and social economic burden worldwide, proper therapeutic strategies seem particularly important for the treatment of osteoporosis. At present, monotherapy is the capital therapeutics and the common drugs for osteoporosis fall into two main categories: (1) bone formation-accelerating drugs, such as fluoride, low dose or intermittent administration of parathyroid hormone (PTH), growth hormone (GH) and isoflavonoids; (2) bone resorption-inhibiting drugs, including calcitonin (CT), bisphosphonates and estrogen (Khosla & Hofbauer, 2017). Despite good clinical curative effects of these medicines in alleviating osteoporosis, the long-term ill-effects, such as renal impairment, dyspepsia and nervous lesion, as well as drug resistance still exist. Therefore, comparing with monotherapy, combination and alternative therapies have been identified as more promising strategies for osteoporosis improvement and management.

Factually, combination and alternative therapies, especially the traditional Chinese medicine (TCM), have been widely used for numerous chronic diseases worldwide. TCM is a holistic medical system with thousands of years of clinical practice (Normile, 2003). Herbal medicines and natural products in TCM always display unique advantages in early intervention, combination therapies and personalized medicine over single-drug treatment for their multi-ingredients, multi-targets, multi-pathways and less toxicity characteristics (Zhao et al., 2015). Preclinical observations and clinical practices increasingly manifest that the nourish kidney herbs and natural products usually receive satisfactory curative effects for bone diseases, which is well consistent with the theory that “deficiency of kidney essence, reduction of marrow and flaccidity of bones” (Zhang et al., 2008, Cai et al., 2015, Murray, 2002). For example, Liu et al. pointed that Drynariae Rhizoma exerted anti-osteoporosis effects through intervening antioxidant-oxidation balance, tryptophan metabolism and phenylalanine metabolism (Liu et al., 2012). Epimedii Herba and its potential active ingredients, such as icariin, could promote the osteogenic action of BMP2 by activating the cAMP signaling pathway and are effective for the prevention and treatment of estrogen deficiency-induced bone loss (Chen et al., 2019, Nian et al., 2009). In addition, the active natural product anemonin, isolated from various Chinese natural herbs, has been revealed to attenuate RANKL-induced osteoclastogenesis and ameliorate LPS-induced inflammatory bone loss in mice through modulation of NFATc1 (Hou et al., 2019). Similarly, Sonchus oleraceus Linn and its main components were reported to protect against LPS-induced sepsis and inhibits inflammatory responses in RAW264.7 cells, thus possessing the potential to improve inflammatory bone loss (Chen, Cui et al., 2019). Therefore, we believe that reinforcing kidney herbal medicines and their active natural products might constitute a safe and important source of drug development for osteoporosis treatment. However, the multi-component, multi-target and synergistic interactions characteristics of herbs make it still a conundrum that how to dissect the multi-scale action mechanisms of herbs and develop novel natural products in osteoporosis treatment at a holistic level. Moreover, the conventional drug development methods usually failed to face the challenges to rapidly develop new drugs.

Fortunately, the advent of systems pharmacology has provided the opportunity and strategy for the investigation of the action mechanisms of herbs and the novel drug discovery. That is, the essence of systems pharmacology in Chinese medicine was to develop a mathematical and computational model for the analysis of complex herbal medicine system. Systems pharmacology was also applied in the discovery of bioactive molecules, the identification of new drug targets, the exploration of therapeutic mechanisms and the exploitation of novel drugs from a whole systematic level (Zhang et al., 2019). In recent years, systems pharmacology has been widely applied in exploration of the multi-scale mechanisms of Chinese traditional herbs in various complex diseases, such as cardiovascular diseases (Zhang et al., 2016), nerve system diseases (Zhang et al., 2019), psychiatric disorders (Wu, et al., 2019), rheumatoid arthritis (Li et al., 2015), depression (Huang et al., 2014). Herein, we adopted a novel systems pharmacology-based approach integrating ADME pharmacokinetics screening, drug targeting and network analysis to explore the therapeutic mechanisms of herbs and their natural products in osteoporosis treatment (Fig. 1). In brief, herbs most associated with osteoporosis were included based on the screening criterion: P < 0.01. Then, the potential active compounds of these herbs with favorable pharmacokinetic properties were screened out by the ADME system. Thirdly, the drug targets of these compounds associated with osteoporosis were predicted through systematic drug-target (SysDT) identification model and database mining. Next, the pivotal disease-relevant biological processes were obtained by the functional enrichment analysis. Meanwhile, the network analysis was implemented to interpret the multi-mechanisms of these herbs in the treatment of osteoporosis. Finally, some key active compounds were selected to verify the analysis results of systems pharmacology.

Fig. 1

Workflow of systems pharmacology approach in treating osteoporosis with herbal medicines. TCMSP: Traditional Chinese Medicine Systems Pharmacology Analysis Platform; ADME: Absorption, Distribution, Metabolism and Excretion; DL: Drug-Likeness; OB: Oral Bioavailability; Caco-2: Caco-2 permeability; C-T network: compound-target network; T-P network: target-pathway network.

Materials and methods

Identification of herbs associated with osteoporosis

To obtain the herbs for the treatment of osteoporosis, a wide-scale text mining was conducted on PubMed and CNKI, using “osteoporosis” and “herbs name” as search terms. Besides, we screened some herbs associated with osteoporosis based on the herb-disease interactions presented in the Traditional Chinese Medicine System Pharmacology Database and Analysis Platform (TCMSP, http://tcmspnw.com): a database of systems pharmacology for drug discovery from herbal medicines (Ru et al., 2014). After removing the duplicates, a list of anti-osteoporosis herbs was constructed preliminarily. Owing to these herbs with different research extents, a statistical index, i.e., P value, the ratio of the number of osteoporosis-herb-related articles/the number of herb-related articles (as displayed in Eq. (1)), was calculated to further assess this bias and further appraise the co-occurrence probability of each herbs and osteoporosis (Zhang et al., 2014).where, N represents the total number of articles published in database PubMed and CNKI, K is the number of articles related to osteoporosis, n shows the number of papers about each single herb and k is the number of papers about the effects of corresponding herbs on osteoporosis. P value indicates the correlations between the herbs and osteoporosis (significant when P < 0.01).

Database construction of herbal ingredients

All ingredients of these herbs associated with osteoporosis were extracted from TCMSP. The chemical structures, drug properties and the SDF format files of these compounds were obtained from online database PharmGKB (https://www.pharmgkb.org/annotated Drugs). Given that the deglycosylation of glycosides by colonic bacteria in humans, the corresponding aglycones of these ingredients were also added into the database for following studies.

Bioactive compounds screening by ADME system

Generally, ADME system as an in silico integrative model is used to predict the pharmacokinetic behaviors and potential drug-drug interactions, which are critical procedures in drug discovery and development. In the present work, three in silico prescreening models, i.e., PreOB (predict oral bioavailability), PreCaco-2 (predict Caco-2 permeability) and PreDL (predict drug-likeness) were employed to screen the bioactive compounds of these included herbs from the TCMSP. Specifically, PreOB, a robust mathematical model, first integrated the main line of defense of limiting the oral bioavailability (OB) of drugs: P-glycoprotein (P-gp) and cytochrome P450s into construction of QSAR modeling for human OB based on 805 structurally diverse drug and drug-like molecules (Xu et al., 2012). This model was verified by the linear (multiple linear regression: MLR, and partial least squares regression: PLS) and nonlinear (support-vector machine regression: SVR) methods with five-fold cross-validation and independent external tests. Compared with the previous prediction models of OB, this PreOB model possesses an optimal predictive ability (R2 = 0.80, SEE = 0.31 for the training set, Q2 = 0.72, SEP = 0.22 for the independent test set) and has been widely and successfully applied for material-based analysis of various Chinese medicines (Wu et al., 2020, Yuan et al., 2020). Caco-2 permeability (Li, Li, Wang, Zhang & Yang, 2007) was introduced to evaluate the absorption rates of the ingredients across the intestinal epithelial barrier. Another ADME model PreDL (Willett, Barnard & Downs, 1998) was developed to discriminate between drug-like and nondrug-like chemicals based on the molecular descriptors and Tanimoto coefficient using the Tanimoto coefficient (as displayed in Eq. (2))where the A shows the molecular properties of ingredients in herbs, and B represents the average drug-likeness index of all compounds in DrugBank database (http://www.drugbank.ca/).

Here, the ingredients matching the criteria: OB ≥ 30%, Caco-2 ≥  − 0.4 and DL ≥ 0.18, were considered as potential bioactive compounds for further analysis. Chemical structures of these compounds were obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov/compound) and the online databases Chemical Book (http://www.chemicalbook.com).

Drug target identification

To obtain the targets of these active natural compounds, we carried out a novel computational model termed SysDT based on Random Forest (RF) and Support Vector Machine (SVM) methods, which integrates large scale information of chemistry, genomics and pharmacology (Yu et al., 2012). The obtained targets were subsequently input into Uniprot (http://www.uniprot.org) database to normalize their name and organisms. All the initially obtained drug targets and their corresponding active compounds were listed as a one-to-one mapping (compound-target list). The gene list of osteoporosis-related therapeutic targets was obtained from the GeneCards and PharmGKB databases which are widely used to search the disease-target associations. After intersecting and matching the gene list of osteoporosis-related therapeutic targets with the compound-target list, the overlapped targets were remained as the potential drug targets, meanwhile, the compounds without targets or osteoporosis-associated targets were removed.

Network construction and functional enrichment analysis

To characterize the multi-component therapies of the anti-osteoporosis herbs, three networks including compound-target (C-T), compound-target-function (C-T-F) and target-pathway (T-P) networks were generated and visualized by an open source of bioinformatics package Cytoscape v3.6.0 (Shannon et al, 2003). The key topological parameter degree is the number of edges associated to the node, which represents the importance of the node in a network.

To further explore the biological activities of the predicted targets in osteoporosis, all the obtained targets were mapped onto DAVID (http://david.abcc.ncifcrf.gov) database for Gene Ontology (GO) enrichment analysis and KEGG pathway analysis. The terms with P < 0.05 were selected to further analysis.

Molecular docking

Molecular docking is one of the common approaches to illuminate the binding modes between the small molecules and their targets. In this section, three pivotal targets (CD40 ligand, interleukin 6 and androgen receptor) for osteoporosis drugs in the C-T network were selected to perform molecular docking simulations by program MOE (Molecular Operating Environment, version 2015.10). The X-ray crystal structures of the three targets were extracted from RCSB Protein Data Bank.

Experimental validation

Materials and reagents

(1) Reagents preparation

Mouse osteoblast cell line MC3T3-E1 cell line was generously provided by Dr. Hong Zhou (University of Sydney, Sydney, Australia). Alpha-Modified Eagle’s Medium (α-MEM) were purchased from ThermoFisher (Gibco, Carlsbad, USA). The fetal bovine serum (FBS) was purchased from Biological Industries (Kibbutz BeitHaemek, Israel). Penicillin, streptomycin and trypsin were purchased from Amresco (Washington, USA). L-glutamine, β-glycerophosphate, L-ascorbic acid, Alizarin red S and Dimethyl sulfoxide (DMSO) were purchased from Sigma-Aldrich (St. Louis, USA). The microplate reader was purchased from Molecular Devices (California, USA). BCIP/NBT Alkaline Phosphatase Color Development Kit was purchased from Beyotime (Shanghai, China). Total RNA Kit was purchased from Omega Bio-tek (USA). PrimeScript RT reagen Kit (Perfect Real Time) and SYBR Premix Ex Taq II (Tli RnaseH Plus) were purchased from Takara (Takara, Dalian, China).

(2) Test sample preparation

To validate the accuracy and efficiency of the systematic pharmacological screening models, some active ingredients were selected for cell experiments. Following the rules of randomness and availability, we sifted five active compounds to quantify the effects of compounds on mouse osteoblast cell line MC3T3-E1 cell. The five selected active compounds were commercially available and represented different pharmacological properties, i.e., relatively optimal pharmacological properties (calycosin and betulinic acid), general pharmacological properties (hederagenin and luteolin) and poor pharmacological properties (asperosaponin VI). These five compounds were purchased from Chroma Biotechnology Co., Itd (Chengdu, China) with purities (≥98%). All the test samples were dissolved in DMSO to make a stock solution of 400 μmol/L, and the final concentrations of DMSO presented in the culture media (<0.1%) had no effect on cell viability.

Cell culture

The MC3T3-E1 cell line was cultured in α-MEM solution supplemented with 10% fetal bovine serum (FBS), 1% L-glutamine and 1% penicillin/streptomycin as well as 1% β-glycerophosphate and L-ascorbic acid. Cell cultures were maintained at a humidified, 37 °C, 5% CO2 incubator (Thermo Fisher Scientific, Waltham, MA).

Cell viability measured by MTT assay

MC3T3-E1 cells were seeded into 96-well plates at a density of 5000 cells/cm2. After incubation for 24 h, the medium was changed into fresh medium containing active compounds with different concentrations from 0.1 µmol/L to 100 µmol/L accordingly, DMSO served as controls. The cells were then incubated at 37 °C for 24 h. For MTT assay, the medium was replaced with 20 µL of MTT (5 mg/mL) and left to incubate for 4 h in dark at 37 °C. After discarding the cultural supernatant, 150 µL/well DMSO was added to dissolve the formazan. Finally, the OD values were read at the wavelength of 490 nm on a microplate reader.

Alkaline phosphatase (ALP) and Alizarin red S (ARS) staining

To examine the effects of these test ingredients on osteogenic differentiation, MC3T3-E1 cells were washed and seeded onto 24-well plates and incubated in complete α-MEM for 24 h for cell adherence and growth. Then, the medium was replaced with the osteoblast inducing conditional media (complete α-MEM supplemented with 1% β-glycerophosphate and 1% ascorbic acid). At the same time, the selected ingredients were added into the osteogenic medium at three concentration gradients (three repeat wells per concentration). The cultures were maintained at 37 °C with 5% CO2, and the medium was replaced every two days. ALP staining was monitored using a BCIP/NBT Alkaline Phosphatase Color Development Kit. Cells were fixed by immersion 4% paraformaldehyde (PFA) solution for 10 min and rinsed in PBS for 5 min 3 times. The samples were then placed in an alkaline phosphatase staining solution (BCIP/NBT solution) for 30 min. The whole procedure was protected from light. After discarding the solution, cells were rinsed in deionized water 2 min 2 times and scanned with a CanoScan 9000F Mark II scanner (Canon, Tokyo, Japan). Alizarin red staining of mineralized osteoblast nodules was carried out. Briefly, cells were fixed and washed as above and then stained in 0.5% Alizarin red S (pH 4.0) for 30 min. After washing, the plates were scanned with a CanoScan 9000F Mark II scanner.

RNA extraction, reverse transcription and quantitative real-time polymerase chain reaction (qRT-PCR)

Total RNA was extracted from MC3T3-E1 cells using TRIzol reagent after 4 d, 7 d and 10 d of the selected compounds treatment. Total RNA (1 μg) was used as a template for double-stranded cDNA synthesis. The SYBR® Premix Ex Taq™ II was applied for qRT-PCR. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as endogenous controls for normalization. Data were analyzed using the comparative Ct method (2-△△Ct) and expressed as fold changes compared to the corresponding control (GAPDH). Primers (sequences listed in Table S3) were synthesized by Tsingke (Xi’an, China).

Statistical analysis

All experiments were independently repeated at least three times with each done in triplicate. Statistical analyses of the data were performed using the GraphPad Prism 6 software (GraphPad Software, La Jolla, CA). All data were reported as the mean ± standard deviation, and P < 0.05 were considered statistically significant for all comparisons.

Results and discussion

Herbal medicines for treatment of osteoporosis

As displayed in Table 1, a total of seven herbs were identified as the most well studied anti-osteoporosis medicines (P < 0.01). All botanical plant names of the included herbs have been matched to the latest name revision in “The Plant List” (www.theplantlist.org). Among them, Drynariae Rhizoma possesses the highest ratio and favorable P value (9.47%, P < 0.01), suggesting its crucial roles in osteoporosis treatment. Drynariae Rhizoma has been widely used as an effective anti-osteoporosis medicine by directly promoting osteoblastic bone formation and inhibiting osteoclastic bone resorption (Liu et al., 2012, Jeong et al., 2005). Moreover, Epimedii Herba, with a well ratio and P value (4.50%, P < 0.01), has a long history of thousands of years in the treatment of osteoporosis in China. Epimedii Herba and its constituents were found to exhibit dual actions in maintaining bone remodeling and skeletal integrity through promoting bone formation and suppressing bone resorption (Wang et al., 2016). The third popular herb is Dipsaci Radix (with P < 0.01), followed by Eucommiae Cortex (with P < 0.01), Dogwood (with P < 0.01) and so forth. Dipsaci Radix and its active ingredients could induce the differentiation of bone marrow mesenchymal stem cells (BMMSCs) into osteoblasts both in vivo and in vitro experiments (Niu et al., 2011). Eucommiae Cortex and Dogwood have been broadly used in TCM prescriptions for tonifying kidney and strengthening bones (Park et al., 2017, Huang et al., 2018). In the following sections, we aimed to investigate the underlying mechanisms of these herbs in the treatment of osteoporosis.

Table 1

Detail information of herbs and their correlations with osteoporosis.

TCMs	Latin scientific names	Genus	Number of articles
			Total (n)	Related to OP (k; ratio)	P value
Drynariae Rhizoma (GSB)	Davallia mariesii T. Moore ex Baker.	Davallia	1847	175 (9.47%)	< 0.01**
Epimedii Herba (YYH)	Epimedium brevicornu Maxim.	Epimedii	8082	364 (4.50%)	< 0.01**
Dipsaci Radix (XD)	Dipsacus asperoides C. Y. Cheng et T. M. Ai	Dipsacus	2363	61 (2.58%)	< 0.01**
Eucommiae Cortex (DZ)	Eucommia ulmoides Oliv.	Eucommia	8229	132 (1.63%)	< 0.01**
Dogwood (SZY)	Cornus officinalis Sieb. et Zucc.	Cornus	4849	43 (0.89%)	< 0.01**
Cistanches Herba (RCR)	Cistanche deserticola Y.C. Ma	Cistanche	3071	27 (0.88%)	< 0.01**
Astragali Radix (HQ)	Astragalus propinquus Schischkin	Astragalus	46,412	215 (0.46%)	< 0.01**

Note: For all the herbs, the abbreviated Chinese names were listed in the brackets after the English name. GSB: Gusuibu; YYH: Yinyanghuo; XD: Xuduan; DZ: Duzhong; SZY: Shanzhuyu; RCR: Roucongrong; HQ: Huangqi. All the names were consistent to that in the Plant List (www.theplantlist.org). n: Total number of researches associated with the corresponding herbs; k: Total number of papers about the effects of corresponding herbs on osteoporosis; ratio: k/n; OP: osteoporosis; ** P < 0.01.

Active compound identification

Bioactive compounds of herbal medicines are considered as the significant contributors of their pharmaceutical effects. In this study, a total of 767 compounds were initially captured from TCMSP and incorporated into the ingredient database. Given that only a few herbal ingredients with favorable ADME properties could exert pharmacological effects, it is requisite to filter the potential bioactive compounds with satisfactory pharmacokinetic properties. As a result, 82 bioactive compounds (10.8%) of the seven herbs fulfilled the following three filter criteria of the reliable in silico ADME model: OB ≥ 30%, Caco2 ≥  − 0.4 and DL ≥ 0.18 simultaneously. In addition, to ensure the integrity of the data, some excluded compounds with relatively poor pharmacokinetic properties were also available for further analysis since they are the well-known main active ingredients in some herbs. For instance, though both chlorogenic acid (OB = 11.93%, Caco-2 =  − 1.03, DL = 0.33) and rutin (OB = 3.20%, Caca-2 =  − 1.93, DL = 0.68) possess poor ADME properties, they are the main ingredients of Eucommiae Cortex and have been confirmed to exhibit potent anti-osteoporosis effects (Wang et al., 2017). Likewise, asperosaponin VI (OB = 1.67%, Caco-2 =  − 3.02, DL = 0.07) and cornin (OB = 12.69%, Caco-2 =  − 1.42, DL = 0.44), as the principal active components of Dipsaci Radix and Dogwood, respectively, have been reported to improve bone loss effectively, thus are regarded as the candidate anti-osteoporosis compounds (Huang et al., 2018, Ke et al., 2016). Therefore, these four compounds were added into the bioactive compound database. As a result, a total of 86 ingredients were considered as the potential active compounds of these herbs (Table S1). Deserved to be mentioned, the predicted ADME properties based on the in silico integrative model were well consistent the reported experimental results. For instance, the bioavailability of orally caffeine in rat was approximately 100%, which was consistent with that predicted by the preOB model (90%) (Samojlik et al., 2016). Besides, OB value of quercetin-3-O-gentiobioside (3.5%) approached to that reported in the literature (3%) (Makino et al., 2009). Yeleswaram et al. reported that the OB of melatonin in rat was 53.5%, which was consistent with the predicted value (53.0%) based on this novel preOB model (Yeleswaram, McLaughlin, Knipe & Schabdach, 1997). Notably, when compared the predicted results from our PreOB model with that from the prediction tools of oral bioavailability (OB) developed by the other groups, we concluded that the trend of the predicted values is basically the same as others (Yoshida and Topliss, 2000, Lipinski et al., 2001). Thus, we highly recommended this in silico ADME screening model as complementary tools in “screening prior to synthesis” procedures for drug discovery. Most of the obtained active compounds have been demonstrated to be involved in osteoporosis treatment, which illustrated the reliability of this system pharmacology-based approach. Some newly identified potential anti-osteoporosis natural compounds, such as gentisin and aureusidin, needed further clinical investigations to confirm their effects and provided a source of phytomedicines as new therapeutics for osteoporosis and bone-related chronic diseases.

Furthermore, the average number of bioactive compounds per herb was 12.3, displaying the multi-component characteristics of herbs. Specifically, these potential active compounds were classified into several categories by their structural elucidation (Fig. 2A) and the ingredients which cannot be assigned into any of these categories were grouped into “others” (Table S1). Notably, flavonoids represent 33% of the active compounds, which is consistent with the fact that flavonoids are biologically major and chemically diverse groups of secondary metabolites and have been implicated having beneficial dietary effects on human health (Wu et al., 2019). Other types including eaters (9.3%), triterpenoids (4.7%), sterols (4.7%) and alkaloids (4.7%) also play important roles in the treatment of osteoporosis. Equally important, compounds of different structures have different biological activities and properties. Take quercetin (MOL57) as an example, it is one of the most commonly polyphenols in human clinical studies, whose absorption and corresponding glycosidesis are associated with the cleavage and release of its aglycones. Quercetin has been proved to exert multiple pharmacological actions including anti-inflammation, oxidation resistance and immune regulation (Yang et al., 2019). While, only the aglycone form of quercetin is available to manufacturers for the supplementation as food products (Teng & Chen, 2018). Here, as an illustration, three representative herbs, i.e., Drynariae Rhizoma, Epimedii Herba and Dipsaci Radix were specified in detail to interpret these filtering principles.

Fig. 2

Classification of active compounds and their targets. A. Classification of compounds; B. Distribution of drug targets according to their biochemical criteria; C. Classification of targets in enzyme; D. Average degree of five main kinds of protein targets, i.e., nuclear receptor, transcription factor, enzymes, cytokine and glycoproteins.

Drynariae Rhizoma

Drynariae Rhizoma is the root of perennial fern Polypodiaceae, which contains various types of ingredients, mainly including flavonoids, triterpenes, phenolic acids, etc. (Song et al., 2017). Due to its favorable biological activities of healing-promotion, anti-osteoporosis and anti-inflammatory, it was generally used to reverse bone loss triggered by inflammation or other pathologic conditions (Lu et al., 2011, Lu et al., 2011). Based on ADME system, 16 bioactive compounds with satisfactory OB, DL and Caco-2 values were screened out from the 71 ingredients of Drynariae Rhizoma. Thereinto, three representative flavonoids, i.e., naringenin (OB = 59.29%, Caco-2 = 0.28, DL = 0.21), eriodictyol (OB = 71.79%, Caco-2 = 0.17, DL = 0.24) and luteolin (OB = 36.16%, Caco-2 = 0.19, DL = 0.25) possess favorable biological activities in the prevention and treatment of osteoporosis. For example, luteolin, a major flavonoid of Drynariae Rhizoma, has the potential to restrain bone loss effectively through inhibiting osteoclast differentiation in postmenopausal osteoporosis (Kim et al., 2011). β-sitosterol (OB = 36.91%, Caco-2 = 1.32, DL = 0.75) is one of the typical steroids of Drynariae Rhizoma, which promotes bone formation by increasing the ratio of OPG (osteoprotegerin) /ODF (osteoclast differentiation factor) in osteoblasts (Zeng et al., 2012).

Epimedii Herba

Epimedii Herba is a traditional tonifying kidney crude herb with significant pharmacological effects. In Chinese Pharmacopoeia, it was recorded that this herb was employed clinically for osteoporosis treatment and bone defect repairment in the TCM prescriptions (Zhao et al., 2016). Also, modern studies have demonstrated that active ingredients of Epimedii Herba exhibit multiple pharmacological activities, which mainly focus on bone system, immune regulation, reproductive system and so on (Zhai et al., 2013). Totally, 130 ingredients were identified from Epimedii Herba, 20 of which were identified as potential active compounds with favorable biochemical properties (OB ≥ 30%, Caco-2 ≥  − 0.4 and DL ≥ 0.18). Among these active compounds, icariin (OB = 41.58%, Caco-2 =  − 1.82, DL = 0.61) is the most abundant and active compound, which possesses favorable biological properties and pharmacological actions in the treatment of osteoporosis. For example, Yao et al. found that icaritin possessed the potential for enhancing the osteogenic differentiation and bone formation (Yao et al., 2012). Another study showed that icariin could promote osteoblasts and bone marrow stromal cells (BMSCs) differentiation via inhibiting the secretion of interleukin-8 (IL-8) and tumor necrosis factor-α (TNF-α), thus increasing bone mass for treating osteoporosis (Wang et al., 2018).

Dipsaci Radix

Dipsaci Radix is derived from the dry roots of Dipsacus asperoides C. Y. Cheng et T. M. Ai, which belongs to Dipsacaceae species (Park et al., 2019). Dipsaci Radix, as one of the high-grade herbs recorded in Sheng Nong’s Classic of Materia Medica, has been commonly applied in the treatment of various orthopedic diseases for several centuries in Asia, such as joint pain, low back pain, bruises, osteoporosis and fracture healing (Mukwaya, Xu, Wong & Zhang, 2014). In addition, Dipsaci Radix could raise the macrophagocyte phagocytosis to reinforce immune function in mice, which is of great importance to bone remolding (Niu et al., 2013). Overall, Dipsaci Radixis is one of anti-osteoporosis herbs with pharmacological activities, including tonifying liver and kidney, strengthening bones and muscles, regulating immune response and resisting aging. Six active compounds with good pharmacological properties (OB ≥ 30%, Caco-2 ≥  − 0.4 and DL ≥ 0.18) were picked out from 31 ingredients of Dipsaci Radix. Therein, asperosaponin VI, the most representative and abundant natural product of Dipsaci Radix, exerts pro-osteogenic, pro-angiogenic and anti-inflammatory response effects in bone fracture treatment (Peng et al., 2010). Besides, asperosaponin VI also acts as an induction of osteoblast maturation and differentiation via increasing bone morphogenetic proteins-2 (BMP-2) synthesis, thus increasing the bone formation (Huang et al., 2018). Moreover, our previous research has found that two active compounds ursonic acid and β-sitosterol of Dipsaci Radix inhibited the osteoclast differentiation in vitro (, Zhang et al., 2019). However, poor pharmacological properties of asperosaponin VI (OB = 1.67%, Caco-2 =  − 3.02, DL = 0.07) may be the most challenge issue and major reason for its ambiguous therapeutic effects as anti-osteoporosis agents in clinical trials. Other study revealed that Dipsaci Radix saponins inhibited osteoclastogenesis through decreasing the ratio of RANKL relative to its decoy receptor OPG, which controls the differentiation of osteoclast precursors (Kong et al., 2012). Indeed, RANKL and OPG are important molecules explicitly linking the bone and immune systems owing to that activated T cells also can secrete RANKL and OPG (Boyle, Simonet & Lacey, 2003).

Moreover, the other four herbs as well as their ingredients also have been widely reported to serve essential roles for osteoporosis treatment. For example, Li et al. found that Eucommiae Cortex was generally contained in the TCM prescriptions of kidney tonifying and bone nourishing (Li et al., 2016). Modern pharmacological studies also showed that Eucommiae Cortex could significantly reduce bone loss via accelerating osteoblastic bone formation and suppressing osteoclastic bone resorption (Ha et al., 2003). Besides, chlorogenic acid, the major constituent of Eucommiae Cortex, promotes the osteogenic differentiation of BMSCs via the Shp2/PI3K/Akt/cyclin D1 pathway in OVX-mice (Zhou et al., 2016). Astragali Radix has been reported to play important roles in various inflammatory diseases, including osteoporosis, periodontal disease and arthritis (Li et al., 2016, Rahman et al., 2018). Astragali Radix extracts could improve the level of sex hormones thus accelerating the differentiation of bone marrow mesenchymal stem cells towards osteoblasts and preventing bone loss in ovariectomized rats (Manolagas, 2010). Additionally, Li et al. uncovered that the oxidative stress in serum and bone tissue of ovariectomized rats were reduced with Astragali Radix, indicating that Astragali Radix may improve osteoporosis by suppressing oxidative stress (Li et al., 2012). Together, these results strongly demonstrated that the seven herbal medicines and their active compounds are effective for osteoporosis treatment through different mechanisms.

Intriguingly, some of the 86 active ingredients could be further metabolized into other biochemical components by gut microbiota. MOL49 (Luteolin) is transformed to baicalein 6-methylether, which processes an anti-inflammatory action and is used to reduce inflammatory bone loss (Amat, De Angelis, Sgamellotti & Fantacci, 2008). MOL57 (Quercetin) is converted into three biochemical components with favorable oxidation resistance and anti-inflammatory action by cleavage of its glycoside, including quercetin-glucuronide, sulfate, and methylated, which could be used to alleviate the bone loss induced by oxidative stress and inflammation (Cao et al., 2015, Biasutto and Zoratti, 2014). Metabolites of MOL54 (Naringenin) are mainly its O-glucoside, 5-O-glucoside, 5-O-glucoside, 7-O-glucoside, 5, 7-O-diglucoside, 7-O-diglucoside, 5-O-diglucoside, 7-O-gentiobioside and 7-O-rutinoside, which participate in the treatment of osteoporosis through their multiple bioactivities, such as immune regulation and antioxidation (Sordon et al., 2016). These results showed that all the metabolites of these compounds are effective in osteoporosis treatment.

Drug target identification and network construction

Drug target identification

Given the multi-component and multi-target characteristics of herbs, the drug target exploration conduces to elucidate the action mechanisms of anti-osteoporosis drugs (Zhou et al., 2016). In this section, a total of 291 targets were predicted for these active compounds initially based on the SysDT model. After being mapped into GeneCards and PharmGKB databases, 58 osteoporosis-associated therapeutic targets were finally remained (Table 2). The average number of targets per herb was 8.3, demonstrating their multi-target characteristics. Among these targets, IL-1β is an important pro-inflammatory cytokine of IL-1 family, which not only stimulates bone resorption but also inhibits bone formation in postmenopausal osteoporosis (Ruscitti et al., 2015). The compound rutin (MOL58) could suppress the expression of IL-1β to promote osteoblast proliferation and differentiation, thus preventing osteoporosis in OVX rats (Abdel-Naim et al., 2018). Moreover, osteopontin (OPN) is one of the abundant non-collagenous proteins in bone matrix secreted by osteoblasts and osteoclasts, which has been linked to many physiological and pathological events, including bone remolding and inflammation (Wein et al, 2019). For example, women with OPN overexpression show less resistance to postmenopausal osteoporosis, indicating that increased OPN is a risk factor for osteoporosis in menopausal women (Chang, Chiang, Yeh, Lee & Cheng, 2010). In addition, the absence of OPN rescued the imbalance of bone resorption and formation induced by the reduced mechanical stress (Ishijima et al., 2001). These results indicated that despite the important roles of OPN in bone remolding, the therapeutic effect of deficiency or neutralization of OPN in bone loss diseases is also remarkable.

Table 2

Detail information of osteoporosis-related targets of seven herbs.

Proteins	Gene symbols	Categories	Uniprot ID
Acetylcholinesterase	ACHE	Enzymes	P22303
Albumin	ALB	Albumin	P02768
Androgen receptor	AR	Nuclear receptor	P10275
Antileukoproteinase	SLPI	Enzymes	P03973
Calcium-activated potassium channel subunit alpha 1	KCNMα1	Ion channels	Q12791
Calmodulin	CALM	Ca-binding protein	P0DP23
Catalase	CAT	Enzymes	P04040
Cathepsin D	CTSD	Enzymes	P07339
CD40 ligand	CD40L	Cytokine	P29965
Collagen alpha-1(I) chain	COL1A1	Collagen	P02452
C-Reactive protein	CRP	Cytokine	P02741
Cyclin-D1	CCND1	Cytokine	P24385
Cyclin-dependent kinase inhibitor 1	CDKN1	Protein kinase inhibitor	P38936
Cytochrome P450 family 1 subfamily A member 2	CYP1A2	Enzymes	P05177
Cytochrome P450 family 3 subfamily A member 4	CYP3A4	Enzymes	P08684
D (2) dopamine receptor	DRD2	G-protein coupled receptor	P14416
Dipeptidyl peptidase IV	DPP4	Enzymes	P27487
E2F Transcription factor 1	E2F1	Transcription factor	Q01094
Epidermal growth factor receptor	EGFR	Cytokine	P01133
Estrogen receptor	ESR	Nuclear receptor	P03372
Estrogen receptor beta	ESRβ	Nuclear receptor	Q92731
Glutathione S-transferase Mu 1	GSTM1	Enzymes	P09488
Heat shock protein family B member 1	HSPβ1	Chaperone	P04792
Heme oxygenase 1	HMOX1	Enzymes	P09601
Insulin	INS	Hormone	P01308
Insulin receptor	INSR	Enzymes	P06213
Insulin-like growth factor II	IGF2	Enzymes	P01344
Insulin-like growth factor-binding protein 3	IGFBP3	Enzymes	P17936
Intercellular adhesion molecule 1	ICAM1	Glycoprotein	P05362
Interferon gamma	IFNγ	Cytokine	P01579
Interferon regulatory factor 1	IRF1	Transcription factor	P10914
Interleukin 1 alpha	IL-1α	Cytokine	P01583
Interleukin 1 beta	IL-1β	Cytokine	P01584
Interleukin-10	IL-10	Cytokine	P22301
Interleukin-2	IL-2	Cytokine	P60568
Interleukin-4	IL-4	Cytokine	P05112
Interleukin-6	IL-6	Cytokine	P05231
Low-density lipoprotein receptor	LDLR	Glycoproteins	P01130
Maltase-glucoamylase	MGAM	Enzymes	O43451
Mitogen-activated protein kinase 1	MAPK1	Enzymes	P28482
Mitogen-activated protein kinase 14	MAPK14	Enzymes	Q16539
Mitogen-activated protein kinase 3	MAPK3	Enzymes	P27361
Mitogen-activated protein kinase 8	MAPK8	Enzymes	P45983
Myeloperoxidase	MPO	Enzymes	P05164
Nitric oxide synthase	NOS	Enzymes	P29474
Nuclear receptor subfamily 1 group I member 2	NR1I2	Nuclear receptor	O75469
Nuclear receptor subfamily 1 group I member 3	NR1I3	Nuclear receptor	Q14994
Osteopontin	OPN	Cytokine	P10451
Oxidized low density lipoprotein receptor 1	OLR1	Glycoproteins	P78380
Peroxisome proliferator activated receptor gamma	PPARγ	Nuclear receptor	P37231
Progesterone receptor	PGR	Nuclear receptor	P06401
Runt related transcription factor 2	RUNX2	Transcription factor	Q13950
Selectin E	SELE	Cytokine	P16581
Activator protein 1	AP-1	Transcription factor	P05412
Transforming growth factor beta-1	TGFβ1	Cytokine	P01137
Tumor necrosis factor	TNF	Cytokine	P01375
Tyrosinase	TYR	Enzymes	P14679
Vascular endothelial growth factor A	VEGFA	Cytokine	P49767

To better clarify the functions of these targets in the pathogenesis of osteoporosis, we analyzed the functional distribution of these targets according to their biochemical classification (Fig. 2B). It was observed that these targets mainly consisted of enzymes, cytokines, nuclear receptors (NRs) and transcription factors (TFs). Of note, enzymes accounted for the largest proportion (35%), including eight transferases (seven of which were kinases), six oxidoreductases, five hydrolases (three of which were proteases) and one synthase (Fig. 2C). In fact, there has been growing interest in exploiting enzymes as the therapeutic targets in recent years. For example, protein methyltransferase (PMTs) is revealed as a novel, chemically tractable drug target class (Copeland, Solomon & Richon, 2009). Moreover, catalase (CAT) is a typical antioxidant in almost all living organisms, which has been shown to remove the excess deleterious reactive oxygen species (ROS) produced in bone metabolism (Zhou et al., 2016). The excess deleterious ROS is one of the important triggering factors for pathological bone loss (Chen et al., 2008). On the one hand, elevated ROS production increases RANKL expression and stimulates RANK-RANKL signaling in osteoblasts and osteoblast precursors, which further promoting osteoclast differentiation and bone resorption (Mercer et al., 2014). On the other hand, ROS also blocks Wnt-β-catenin signaling in human MSCs (HMSCs), thus restraining the osteogenic differentiation and resulting in excessive bone loss (Chen et al., 2010).

Cytokines including polypeptides and proteins were the second major category of drug targets (28%), which mostly expressed on various immune cells and bone cells (Azizieh, Raghupathy, Shehab, Al-Jarallah & Gupta 2017). An impressive number of studies have shown that certain cytokines produced by immune cells or bone cells, such as TNF-α, IL-1, IL-6, IFNγ, TGF-β and IL-10, are involved in the pathogenesis of osteoporosis (Abdel Meguid et al., 2013; Zupan et al., 2013; Takayanagi, 2009; Kapoor et al., 2011). For instance, the increased expression of TNF-α could upregulate RANKL and M−CSF generation, thus stimulating osteoclastic bone resorption and aggravating estrogen deficiency-induced bone loss (Zha et al., 2018). Additionally, the vital roles of interleukins in the bone destruction and bone loss have been demonstrated both in clinical trials and in animal models of osteoporosis (Kapoor, Martel-Pelletier, Lajeunesse, Pelletier & Fahmi, 2011; Meguid, Hamad, Swilam & Barakat, 2011). Take pleiotropic cytokine IL-17 as an example, it is secreted by Th17 cells and potently potentiates RANKL-induced osteoclastogenesis, thus intensifying bone loss (Li et al., 2015). Another pivotal cytokine is the T cell costimulatory factor CD40L, which has been reported to intensify the dysregulation of bone resorption and formation under both estrogen deficiency and continuous PTH treatment through stimulating the additional production of RANKL and inhibiting the secretion of OPG (Gao et al., 2008, Li et al., 2011). Overall, these observations strongly support the fact that the cytokines are key therapeutic targets for osteoporosis and offer proper proof for the intricate nexus between osteology and immunology (Arron & Choi, 2000).

NRs (12%) serve as potential promising drug targets for osteoporosis due to their regulatory roles in bone development and remodeling (Zuo & Wan, 2017). Drugs targeting nuclear receptors are gradually popular in clinical treatment of osteoporosis through modulating bone formation and resorption rates. For example, drugs like estrogen receptor (ESR) agonist remain important in the prevention of postmenopausal osteoporosis clinically (Riggs & Hartmann, 2003). Another important kind of drug target is TF (or sequence-specific DNA-binding factor), occupying 7% of all the targets, which controls the transcription rates of genes from DNA to mRNA and regulates the gene expressions by binding to a specific DNA sequence (Lambert et al., 2018). The average degree of these five types of drug targets was shown in Fig. 2D. Herein, NRs possesses significantly more connected components than the other types, suggesting that they can bind to various specific ligands and may exhibit multiple physiological effects.

Compound-target network analysis

After discarding the compounds without targets or with no osteoporosis-related targets, the remained 54 candidate compounds (blue hexagon) and their 58 targets (orange circle) were visualized in a directed bipartite compound-target (C-T) network (Fig. 3), consisting of 112 nodes and 238 interactions (edges). In this network, the node size was proportional to its degree, thus, the node with multiple join points was considered as the key node in the network. Specifically, each compound is connected to an average of 4.67 targets, suggesting their broad pharmacological properties. Note that the top three compounds were quercetin (MOL57, degree = 39), luteolin (MOL49, degree = 19) and kaempferol (MOL46, degree = 16), indicating that they might be the crucial compounds in osteoporosis treatment. Coincidently, we also found that several active ingredients could act together to one common target protein, which might exhibit additive effects for improving the osteoporosis outcome. Specifically, there were 24 (41.4%) targets attached to ≥ 2 ligands, of which CALM (degree = 28), AR (degree = 26), DPP4 (degree = 26) and PGR (degree = 21) linked to>20 compounds. These results suggested that the active compounds of these herbs may possess synergic combination effects in osteoporosis treatment. For example, luteolin, with specific anti-inflammatory properties, has the potential to decrease osteoclastogenesis and osteoclastic bone resorption to prevent bone loss by suppressing the expression of multiple cytokines, such as TNF-α, IL-6, IL-8, IFNγ and so on (Kim et al., 2011, Seelinger et al., 2008). Moreover, IL-6 induces osteoclast activity and increases bone resorption, resulting in excessive bone loss. Multiple researches reported that natural products such as quercetin (MOL57), rutin (MOL58) and calycosin (MOL21) could reduce the production of the osteoclastogenic cytokine IL-6 to treat osteoporosis effectively (Ivanova, Vasileva, Ivanova, Peikova & Obreshkova, 2015). In summary, as displayed in the C-T net, natural active compounds of botanical medicines and their targets may engage in the complicated interactions to exhibit the synergistic therapeutic actions for osteoporosis therapy and prophylaxis.

Fig. 3

Compound-target (C-T) network was constructed by linking the candidate compounds (blue hexagons) with their potential targets (orange circles). The node size was proportional to its degree.

Compound-target-function network analysis

Given the complicated pathophysiology of osteoporosis, the use of combination for osteoporosis therapy and prophylaxis is especially rational. Emerging evidence indicated that therapies that restore bone homeostasis, reduce inflammatory reactions or affect the immune responses exerted promising therapeutic potentials for osteoporosis (, Zhang et al., 2019, Wang et al., 2017). Coincidentally, the compound-target-function (C-T-F) network of anti-osteoporosis herbs (Fig. 4), consisting of 51 compounds, 57 candidate targets and their major functional annotations, showed that the potential active compounds and their associated targets were involved in three main fundamental processes: inflammation, metabolism and immunity. Specifically, 21 targets (accounting for 34.5% of the total targets) were mainly associated with inflammation (cytokines), following 18 targets and 14 targets were involved in immune and metabolism regulation, respectively. Interestingly, eight targets (IL10, IL-1, IL2, IL4, IL6, JUN, OPN and CRP) engaged in the regulation of both immune and inflammation processes, which might display a combination therapy for effective long-term treatment for osteoporosis patients. Here, some typical inflammation cytokines (factors) relevant to osteoporosis were screened out and listed in Table 3. For example, CD40L (CD40 ligand), an important T-cell costimulatory molecule, also known as CD154, exerts its function by binding to CD40, which is expressed on antigen-presenting cells (APCs), stroma cells (SCs) and osteoblasts (OBs) (Quezada, Jarvinen, Lind & Noelle, 2004). Several researches have reported the pivotal roles of CD40L in bone metabolism. For example, Li et al. pointed that CD40L combined with its costimulatory receptor CD40 could intensify the bone loss induced by ovariectomy (Li et al., 2011). Also, Gao and his colleagues uncovered that CD40/CD40L costimulatory system increased osteoclastogenic activity of SCs through stimulating the additional production of RANKL and inhibiting the secretion of OPG on SCs under continuous infusion of PTH (Gao et al., 2008). In addition, TNF-α, an important mononuclear-macrophage-derived cytotoxin, contributes to the pathological damage of some inflammatory diseases, such as osteoporosis (Bystrom et al., 2018). The increased expression of TNF-α could aggravate the bone loss through promote the osteoclastic bone resorption, which has been observed in humans and experimental models of osteoporosis (Zuo and Wan, 2017, Sang et al., 2017).

Fig. 4

Compound-target-function (C-T-F) network was constructed by ingredients (circles) and their corresponding protein targets (blue hexagons).

Table 3

Function of inflammation factors (cytokines) in bone-immune system.

Factors	Sources	Effects on immune system	Functions in bone Metabolism	References
IL-6	Dendritic cells (DCs), Macrophage	Pro-inflammation, Th17 induction	Activation of osteoclastogenesis	Song, Gao & Qian, 2014
IL-4	Th2	Humoral Immunity	Inhibits osteoclastogenesis	Mangashetti, Khapli & Wani, 2005
TGF-β	Multiple cell lines	Blocks activation of lymphocytes and monocytes derived phagocytosis	Indirect osteoclast activation. Inhibits osteoblast differentiation	Adamopoulos & Bowman, 2008
IL-1	Macrophage and DCs	Pro-inflammation	Directly activates RANK signaling to promote osteoclastogenesis	Adamopoulos et al., 2010
IL-10	Treg	Anti-inflammation	Suppress bone resorption	Wing, Yamaguchi & Sakaguchi, 2011
IFN-γ	Th1, NK cells	Cellular immunity	Inhibits osteoclastogenesis	Kotake et al., 2005
TNF-α	Th17, Macrophage DCs	Pro-inflammation	Indirect osteoclastic activation throughRANKL	Boyce & Xing, 2008
CD40L	antigen-presenting cell (APC), stroma cell; T cell; OB	Pro-inflammation	Indirect osteoclastic activation throughRANKL	Gao et al., 2008
IL-2	T cell	Pro-inflammation	Activation of osteoclastogenesis	Sun, Niu & Qi, 2006
JUN (transcription factor AP-1)	Th17, B cell and dendritic cell	Pro-inflammation	Activation of osteoclastogenesis	Wagner & Eferl, 2005
RUNX2	endothelial cells	Anti-inflammation	stimulate the differentiation of osteoblasts	Kawane et al., 2018
OPN	dendritic cells (DCs), monocytes, osteocyte	Anti-inflammation	stimulate the differentiation of osteoblasts	Filip, Radzki & Bieńko, 2018
CALM	chondrocytes and articular cartilage	Anti-inflammation	Inhibits osteoclastogenesis	Wu, Ahn, McKenna, Yeo & McDonald, 2005
CRP	lymphocytes	Pro-inflammation	CRP concentration was inversely associated with BMD	Pablo, Cooper & Buckley, 2012
PPARγ	T cell	Pro-inflammation	Inhibits osteoclastogenesis	Kawaguchi et al., 2005

In fact, osteoporosis has been recognized as one of the most common immune-related bone diseases under inflammatory conditions. On the one hand, osteoblasts and osteoblast-mediated bone formation could be affected by soluble factors such as cytokines in the immune system, including CD40L, TNF-α, IL-1, IL-6 and IL-4 (Takayanagi, 2007, Gilbert et al., 2000). On the other hand, osteoclasts and osteoclast-mediated bone resorption could be triggered by multiple factors such as RANKL, M−CSF, CD40L and IL-17 (Chen et al., 2008, Banuelos and Lu, 2016, Tyagi et al., 2012). Therefore, cytokines or immune factors as well as inflammatory responses generated by the aberrant activation of immune system may influence the bone remolding process via altering the delicate balance of osteoblastic bone formation and osteoclastic bone resorption (Dar, Azam, Anupam, Mondal & Srivastava, 2018). Intriguingly, Tyagi et al. found that both immune cells and bone-resorbing osteoclasts are derived from hematopoietic stem cells (HSCs) (Tyagi et al, 2012), which further provided reliable proof for the reciprocal regulation between bone and immune systems. Overall, these results indicated that anti-osteoporosis herbs and their active ingredients could synergistically improve and alleviate symptoms of osteoporosis by multiple mechanisms, such as suppressing bone resorption, promoting bone formation, preventing inflammatory reaction as well as modifying immune response and other osteoporosis risk factors.

Target-pathway network analysis

For better elaborating the holistic mechanisms of anti-osteoporosis medicinal herbs, we extracted the canonical pathways that are highly relevant to osteoporosis from KEGG database (http://www.genome.jp/kegg/), resulting in 37 pathways including Wnt signaling pathway, TNF signaling pathway, osteoclast differentiation. Subsequently, all the drug targets of these herbs were mapped onto these 37 pathways, generating a bipartite target-pathway (T-P) network graph as displayed in Fig. 5. The T-P network is constructed with 77 nodes (40 targets and 37 pathways) and 240 edges after discarding 17 targets nonparticipation of any pathways. We observed that more than one-third targets (14/40) participated in multiple pathways (degree ≥ 5), indicating that these targets may play a key role in the pathogenesis of osteoporosis through various biological processes. Note that the top three pathways were mitogen activated protein kinase 1 (MAPK1, degree = 28), tumor necrosis factor (TNF, degree = 18) and transcription factor AP-1 (JUN, degree = 17). For instance, MAPK1, also known as ERK2, participated in almost all physiological and pathological processes of the organism. In addition, TNF is mainly involved in pathways associated with immune system which is another crucial mechanism of osteoporosis. As mentioned above, overexpression of TNF-α could stimulate the production of RANKL and M−CSF, thus increasing osteoclastic bone resorption and aggravating estrogen deficiency-induced bone loss (Zha et al., 2018). JUN, a dimeric transcription factor complex composed of Fos, Jun and activating transcription factor (ATF) families of proteins, has attracted increasing attention as its crucial roles in pathological bone loss. Harada S et al. found that JUN can be activated by TGF-β, PTH and vitamin D, which in turn negatively regulates the differentiation and proliferation of osteoblasts, thus exacerbating the bone loss (Harada & Rodan, 2003). Moreover, JUN also acts as an osteoclastogenic transcription factor to modulate the osteoclastic bone resorption by its downstream target genes like NFATc1 (Matsuo et al., 2004).

Fig. 5

Target-pathway (T-P) network, where ellipse and round rectangle nodes represent osteoporosis-related targets and pathways, respectively. Color codes were given in the legend. Node size was proportional to its degree.

In addition, the KEGG pathway analysis showed that the drug targets were primarily implicated in processes associated with skeletal system, immune system, endocrine system and signal transduction, which might be the key mechanisms that drugs engender their anti-osteoporosis effects (Zheng and Spector, 2012). The main pathways included PI3K-Akt signaling pathway (degree = 14, P = 1.10E-14), Toll-like receptor signaling pathway (degree = 10, P = 1.93E-14), TNF signaling pathway (degree = 10, P = 2.84E-14), Th17 cell differentiation (degree = 9, P = 1.47E-12) and Osteoclast differentiation (degree = 7, P = 1.22E-08) (Fig. 5 and Supplementary Table S2). For example, Gu et al. has found that the activation of PI3K/Akt signaling pathway promotes the differentiation of MC3T3-E1 cell into osteoblasts, playing important roles in bone stability and bone reconstruction (Gu et al., 2013). Another study also showed that PI3K/Akt signaling pathway could accelerate the HMSCs to differentiate into osteoblasts (Baker et al., 2015). With respect to Th17 cell differentiation, it serves as a bridge between immune system and skeletal system. Th17 cells are recognized to be pathogenic subset of CD4+ T cells in osteoporosis owing to their potent osteoclastogenic activities (Wing, Yamaguchi & Sakaguchi, 2011). Both activated Th17 cells and Th17 cells-derived IL-17 participate in the regulation of osteoclastogenic activity and stimulate the abundant production of additional RANKL under inflammatory conditions, resulting in intensifying bone loss (Li et al., 2015). All these results exactly confirmed to the prior interesting observations that immune response maybe the pathogenesis mechanism of osteoporosis. Furthermore, osteoclasts play central roles in the physiological bone homeostasis and pathological bone loss. Osteoclast differentiation signaling pathway is closely associated with bone diseases such as rheumatoid arthritis and osteoporosis (Teitelbaum, 2000). Zhu et al. pointed that suppressed RANKL-induced osteoclast differentiation could reduce OVX-induced bone loss in mice (Zhu et al., 2013). Therefore, we speculated that herbal medicines probably mediate these pathways to exhibit the anti-osteoporosis properties, and thereby might provide a combining system for osteoporosis prevention and treatment.

Osteoporosis pathways analysis

Considering the pathogenesis and pathological state of osteoporosis, in this section, an integrated ‘Osteoporosis Pathway net’ that comprises of six signaling pathways such as osteoclast signaling pathway, T cell receptor signaling pathway, TNF signaling pathway, MAPK signaling pathway and PI3K-Akt signaling pathway were assembled. As shown in Fig. 6, 31 proteins (pink rectangles) located from upstream to downstream on the osteoporosis pathway can be linked with active ingredients in our work, indicating that these herbs and their compounds may antagonize osteoporosis by the regulation of these target proteins on the pathways. Besides, there exist multiple biological cross-talks among these pathways, where several mutual targets link these pathways together to achieve anti-osteoporosis effects. The most typical representative is the joint effects of osteoclasts differentiation, T cell receptor, MAPK and TNF signaling pathways. These pathways are bonded together to regulate MAPK14 (also called p38α MAPK), which is detected in activated immune cell macrophages and activated by some cytokines and growth factors involved in bone development and remodeling, such as BMP2, TGF-β and TNFα (Yamashita et al., 2008, Cuadrado and Nebreda, 2010). MAPK14 expressed on BMMSCs facilitates bone health by promoting osteogenic differentiation and bone formation (Cong et al., 2016). Moreover, PI3K was a key intersection for TNF, MAPK, estrogen and PI3K-AKT signaling pathways, while JUN was the cross target of osteoclast differentiation, T cell receptor, MAPK and estrogen signaling pathways.

Fig. 6

Distribution of target proteins on representative osteoporosis-related pathways (purple squares) and their therapeutic modules (blue oval). Arrows indicate activation, T-arrows indicated inhibition. Pink squares represent the predicted target proteins of these herbs; Yellow squares represent the proteins on these pathways.

Notably, target proteins in this integrated ‘Osteoporosis Pathway net’ take control of several therapeutic modules, such as immune response, cell proliferation, cell cycle progression, cell survival and so forth (Fig. 6). Here, we mainly focused on three representative modules to dissect the synergistic effects of these herbs in osteoporosis treatment.

Immune response module: As shown in Fig. 6, T cell receptor signaling pathway, TNF signaling pathway, PI3K/Akt pathway and osteoclast differentiation were involved in the modulation of immune responses, which is consistent with the results of C-T-F network and T-P network. In fact, some studies have pointed that immune cells and immune cell-derived cytokines broadly participant in the development of osteoporosis, therefore, the immune response maybe the mechanism of osteoporosis (Clowes, Riggs & Khosla, 2005). The immune cytokines involved in these pathways, such as TNFα, IL-2, IL-6, IL-10 and CD40L, play crucial roles in the pathogenesis of osteoporosis under inflammatory conditions.

Cell proliferation module: T cell receptor (TCR) signaling pathway and MAPK signaling pathway were found to be closely associated with the cell proliferation module. For example, MAPK signaling pathway participates in the proliferation of vascular endothelial cells, which further engages in the regulation of angiogenesis and osteoblasts differentiation, thus mediating bone metabolism (Zhang and Liu, 2002, Duttenhoefer et al., 2015). Moreover, Chang et al. have revealed that T cells activated by TCR signaling pathway is sufficient to stimulate T-cell proliferation, which feeds forward on enhancing the pathway of Th17 differentiation, thus leading to inflammatory bone loss (Chang et al., 2003).

Cell cycle progression module: Cell cycle is an ordered set of events that eventually leads to cell growth and division. As shown in Fig. 6, MAPK signaling pathway and PI3K/Akt signaling pathway were related to the cell cycle progression module. Some evidences have revealed that the PI3K/Akt pathway is associated with regulation of cell cycle progression (Sherr & Roberts, 1999). The cyclin-D1 (CCND1) and cyclin-dependent kinase inhibitor 1 (CDKN1/p21) located in PI3K/Akt pathway form complexes to regulate cell cycle progression through various cell cycle stages. Many investigations have reported the important roles of p21 (CDKN1) in adjusting cyclin-dependent kinase (CDK) activity. Under physiological conditions, p21 may potentially induce the activities of CDK4 and CDK6 from early G1 until the middle of S phase by binding to CDK4/CCND1 and CDK6/CCND1 complexes (Taylor et al., 2002).

Taken together, all these results indicated that the active compounds and drug targets of herbal medicines probably conjunctively act through multiple mechanisms, such as suppressing inflammatory response, maintaining the balance of bone metabolism as well as improving the organism immunity, to synergistically benefit patients with osteoporosis, which may provide a novel therapy strategy for the treatment of osteoporosis.

Molecular docking is of fundamental importance in modern structure-based drug design (Taylor, Jewsbury & Essex, 2002). To clarify the reliability and accuracy of the binding modes between the compounds and their targets in herbs, we selected four compounds (luteolin, quercetin, rutin and kaempferol) and their corresponding targets (CD40L, IL-6 and AR) for docking simulation. Presently, 5 ns MD simulations for all the docked complexes were carried out to observe the kinetic conformational changes between the compounds and the targets in aqueous solution. The snapshots of 3D binding conformations of each complex with the key amino acids of their target proteins at the last 1 ns of the MD simulations were shown in Fig. 7. Obviously, CD40L-luteolin and CD40L-quercetin (Fig. 7A and B) were stabilized by the interactions between the ligand and Asp26, Glu23 and Lys27. Fig. 7C, demonstrated that kaempferol is directed towards the binding pocket in the entrance cavity of AR, establishing interactions with residues His789, Trp796 and Gln792. Besides, kaemferol, rutin and quercetin were located within the binding cavity of IL-6 (Fig. 7D, E and F), all which form at least two hydrogen bonds and an additional water-mediated hydrogen-bonding to stabilize their binding sites. Therefore, the results obtained from the molecular docking showed that hydrogen bonds and water-mediated hydrogen-bonding play central roles in the protein–ligand recognition and stability, which may contribute to assessing the activity of anti-osteoporosis drugs.

Fig. 7

3D Binding conformations of compounds and their target proteins from MD simulation. (A) CD40L-luteolin, (B) CD40L-quercetin, (C) AR-kaempferol, (D) IL-6-kaemferol, (E) IL-6-quercetin, (F) IL-6-rutin. The molecules were present as ball and stick models.

Experimental verification

To verify the accuracy of the prediction results, five commercially available active compounds with different pharmacological properties (calycosin < MOL21>, asperosaponin VI < MOL16>, hederagenin < MOL40>, betulinic acid < MOL51 > and luteolin < MOL49 > ) were selected to investigate their effects on osteoblast precursor cells MC3T3-E1 cell line.

MTT assay

The effects of these five compounds on MC3T3-E1 cell proliferation were examined by MTT assay (Fig. 8A). The results showed that calycosin and asperosaponin VI did not cause cytotoxic responses on MC3T3-E1 cells even at high concentrations (up to 100 μmol/L) (Fig. 8A (a − b)). In addition, we observed that there was no obvious toxicity to MC3T3-E1 cell when hederagenin, betulinic acid and luteolin were used at doses<50 μmol/L (Fig. 8A (c − e)). Moreover, it was apparent that all five compounds could dramatically promote MC3T3-E1 cells proliferation at appropriate concentrations ranges. Based on these observations, the five active compounds were used at the dose range of 0.1 μmol/L to 10 μmol/L for the following evaluations.

Fig. 8

Effects of the active compounds on osteoblast proliferation and differentiation (mean ± SD, n > 3). A. The cell viability of MC3T3-E1 cells after treated by the active compounds for 24 h. B. ALP staining of MC3T3-E1 cells for 48 h. C. Alizarin red mineralized nodules staining of MC3T3-E1 cells for 14 d. *P < 0.05, **P < 0.01, ***P < 0.001 vs Blank group, ns: no statistical differences).

ALP staining and Alizarin red S (ARS) staining

ALP is a well-characterized early marker of osteogenic differentiation. Mineralization is a major performance index of osteogenic differentiation and bone regeneration. To detect the roles of the five active compounds on osteogenic differentiation of MC3T3-E1 cell, the dose-dependent effects and optimal concentrations need to be determined by ALP staining and Alizarin red mineralized nodules staining. The MC3T3-E1 cells showed general increase of ALP activity after treated with the active compounds (0.1, 1.0 and 10.0 μmol/L) and the number of ALP+ cells was significantly more than that in the blank controls. Dose-dependent study revealed that the optimal concentration of all the five compounds for increasing ALP activity was at 1.0 µmol/L (Fig. 8B (a − b)). In addition, alizarin Red S (ARS) staining showed that at day 14, there was more calcium nodules formation in the presence of active compounds compared with the controls. Largely consistent with the results of ALP staining, all these active compounds except luteolin increased calcium deposition with the maxima at the concentration of 1 µmol/L (Fig. 8C (a − b)). From the results mentioned above, we observed that all the five key compounds effectively enhanced the osteoblast proliferation and differentiation of MC3T3-E1 cells under the concentration gradients of 0.1 µmol/L, 1.0 µmol/L and 10 µmol/L, of which 1.0 µmol/L was the optimal concentration.

Marker gene expressions of osteoblast differentiation by qRT-PCR

To further ascertain the effect of these five active compounds on osteoblast differentiation, the expression levels of marker genes including ALP, RUNX2 and COL Iα1 were examined by qRT-PCR (Fig. 9A − E). In Fig. 9A, the expression of ALP was significant promoted at day 4 (P < 0.05) and 7 (P < 0.01) with calycosin treatment. RUNX2 was increased by about 34% (P < 0.05), 86% (P < 0.01) and 85% (P < 0.01) after treated with calycosin for 4, 7 and 10 d, respectively. Besides, COL Iα1 was significant up-regulated at day 7 (P < 0.05) and 10 (P < 0.01) after being treated with calycosin. The same alterations of these marker genes were also observed after asperosaponin VI treatment (Fig. 9B). In addition, hederagenin significantly up-regulated the expressions of ALP and RUNX2 at day 4, 7 and 10 (P < 0.01). Both at day 7 and 10 of hederagenin treatment, the expression of COL Iα1 was increased (Fig. 9C). As for betulinic acid (Fig. 9D), the expression of ALP (P < 0.01), RUNX2 (P < 0.01) were significant up-regulated at day 4 and 7. The COL Iα1 was up-regulated significantly only after 7 days of hederagenin treatment (P < 0.01). Luteolin could significantly increase the expression of ALP, RUNX2 and COL Iα1 at day 4 and 7 (P < 0.01), while only RUNX2 and COL Iα1 were up-regulated after treated with luteolin for 10 d (P < 0.01) (Fig. 9E).

Fig. 9

Effect of active compounds on osteoblast differentiation. The mRNA expression level of osteoblast differentiation marker genes: ALP, RUNX2 and COL Iα1 in MC3T3-E1 cells treated by the active compounds for 4, 7 and 10 d detected by qRT-PCR. The concentrations of the active compounds were at 1.0 μmol/L (mean ± SD, n > 3, *P < 0.05, **P < 0.01 vs Blank group, ns: no statistical differences). A. Calycosin; B. Asperosaponin; C. Hederagenin; D. Betulinic acid; E. Luteolin.

Coincidentally, all the five selected active compounds are proven to be dramatically efficient in promoting osteoblast proliferation and differentiation, which is quite consistent with the previous reported pharmacological action of them on curing osteoporosis. For instance, Kong and colleagues revealed that calycosin exhibited a positive effect regarding improving the osteogenic function of osteoblasts (Kong, Wang, Niu, Wu & Pan, 2018). In addition, both asperosaponin VI and betulinic acid were demonstrated to synergically enhance bone formation through stimulating BMP/Runx2/β-Catenin signals or Smad/p38 pathways, which are essential for the management of osteoporosis (Lo et al., 2010, Niu et al., 2011; Choi et al., 2016). Similarly, it was reported that luteolin promoted osteoblastic differentiation by regulating the ERK/Lrp-5/GSK-3β pathway to attenuate glucocorticoid-induced osteoporosis (Jing et al., 2018). Notably, an authorized patent by Chinese National Patent Office showed that hederagenin can significantly inhibit the decrease of bone mineral density in ovariectomized rats and improve the internal characteristics of femur (Li and Qi, 2018).

Results of the experiments showed that all the five selected compounds within appropriate concentrations could promote the osteoblast proliferation and differentiation, which were in good agreement with our theoretical predictions, indicating the reasonability of the system pharmacology-based approach.

Conclusion

Osteoporosis has long been a pervasive public health concern and caused significant economic losses worldwide, the preferable prevention and management of osteoporosis were particularly important. In China, herbs and their natural active compounds have historically been accepted as an important source of therapeutic agents for osteoporosis. In this study, we mainly presented an integrated system pharmacology-based approach to dissect and understand the multi-scale mechanisms that underlie the spread effects of Chinese herbs from molecular-level interactions to organismal-level phenotypes in osteoporosis therapy.

In summary, the main findings of this work are as follows:

After wide-scale text mining and the P value evaluation, seven herbal medicines like Drynariae Rhizoma, Epimedii Herba, Dipsaci Radix, etc. have been shown exhibiting significant correlations with osteoporosis.

Based on the ADME system and SysDT method, 86 compounds with the favorable pharmacokinetic profiles and their 58 targets from the seven herbs were identified to be implicated with osteoporosis probably. Notably, most the potential active compounds have been demonstrated to be involved in the mechanism of osteoporosis, which further illustrated the reliability of this system pharmacology-based approach. Some other obtained novel potential anti-osteoporosis natural compounds, such as gentisin and aureusidin, needed more experiments to prove their effects and provided a source of phytomedicines as new therapeutics for osteoporosis and bone-related chronic diseases.

The T-P network and the integrated “osteoporosis pathway net” indicated that the active compounds and drug targets of herbal medicines probably conjunctively work by multiple mechanisms, such as suppressing inflammatory response, maintaining the balance of bone metabolism as well as improving the organism immunity, to synergistically benefit patients with osteoporosis.

The experiment results showed that calycosin, asperosaponin VI, hederagenin, betulinic acid and luteolin could promote the osteoblast proliferation and differentiation at the proper concentrations, which strongly supported the potential bioactive compounds we identified based on systems pharmacology analysis.

Taken together, the natural active compounds identified based on the present novel system pharmacology approach may provide a source of phytomedicine for osteoporosis and bone-related chronic diseases, which would be of great value. Importantly, this system pharmacology approach contributes to understand the intricate associations among biological organisms, drugs and chronic diseases from a network perspective, as well as provides a novel approach to promote drug discovery. Despite this system pharmacology model has been widely applied and exhibited great influence in the development of novel drugs, it is still in its infant stage and showed a few flaws such as the insufficient accuracy and the limited herbs. Besides, more tests are needed to confirm the effects of the novel potential anti-osteoporosis compounds identified in our study, such as gentisin and aureusidin. Therefore, in the follow-up researches, more experimental data should be collected for the continuous model improvement and optimization.

Author contribution

Ying Huai wrote the first draft of the manuscript and drew the figures; Wen-juan Zhang and Wei Wang helped to do the data analysis; Ai-rong Qian, Yu Li and Wen-juan Zhang helped to design the experiment; Shan-feng Jiang, Kai Dang and Zhi-ping Miao guided the experiment; all the other authors revised the manuscript. All the authors proof read and approved the final manuscript.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.