Abstract
Osteoporosis is a systemic skeletal disease that can easily lead to bone fractures. Berberine has been shown to be effective in treating osteoporosis. This study was conducted to identify the potential mechanism of berberine in treating this complaint. We screened potential targets of berberine and identified the osteoporosis-related differentially expressed genes (DEGs) in the microarray dataset GSE56815. Protein–protein interaction (PPI) network construction, hub targets identification, and pathway enrichment were carried out to find the potential targets. Molecular docking and molecular dynamics studies were performed to verify the combination of berberine with its treatment-related central targets. In addition, SwissADME preliminarily evaluated the physicochemical properties of berberine. Through data mining, 23 osteoporosis-related targets of berberine were selected. PPI and module analyses suggested that AKT1, MAPK1, ESR1, AR, TP53, and PTGS2 are the core targets of berberine. Docking and molecular dynamics studies showed that berberine could stably bind to core proteins to form a protein–ligand complex. The enrichment analysis showed that the estrogen signaling pathway and thyroid hormone signaling pathway play important roles in curing osteoporosis. To sum up, berberine primarily acts on AKT1, MAPK1, ESR1, AR, TP53, and PTGS2, mainly regulating the estrogen and thyroid hormone signaling pathways to treat osteoporosis in a multi-target, multi-pathway, and multi-system manner.
Keywords
Introduction
Osteoporosis is a systemic osteopathy characterized by micro-architectural deterioration of the bone tissue and low bone mass, with a consequent increase in susceptibility to bone fragility and fracture. 1 The main manifestations of osteoporosis are fracture (particularly hip and vertebral fractures), 2 hunchbacks, and diminished respiratory function, which typically occur in elderly and postmenopausal women. Osteoporosis is caused by factors such as micro-damage, genetic factors, estrogen deficiency, aging, and oxidative stress. 3 Currently, drug treatment is the principal treatment for osteoporosis, with the options including not only estrogen and a selective estrogen receptor modulator, 4 but also bisphosphonate, 5 the parathyroid hormone analog teriparatide, 6 and a human monoclonal antibody to the receptor activator of nuclear factor-κB ligand, denosumab. 7 Although considerable progress has been made in the development of drugs for treating osteoporosis, their side effects are not negligible. For example, estrogen increases the risk of cardiovascular events and breast cancer, 8 and bisphosphonates cause atypical femur fractures and osteonecrosis of the jaw. 9 Therefore, a safe and effective alternative for osteoporosis treatment is required (Figure 1).

Flowchart of this study.
Traditional Chinese medicine (TCM) has multiple pharmacological effects and mild adverse effects. Modern pharmacological studies have revealed several TCMs with good curative effects for osteoporosis via various mechanisms, including antioxidant, anti-inflammatory, and estrogenic-like activities. 10 Therefore, active compounds in TCM may be useful for treating osteoporosis. The phytochemical berberine (Figure 2) is derived from herbs such as Coptis chinensis and Radix Berberidis. 11 Recent studies showed that berberine is effective against many diseases, such as diabetes, hypercholesterolemia and inflammation, 12 as well as osteoporosis. 13

Chemical structure of berberine.
In ovariectomized-periodontitis rats, berberine protects the intestinal barrier through increasing butyrate, further preventing interleukin-17A-associated immune imbalance in alveolar bone and serum, consequently refraining estrogen deficiency-related alveolar bone loss. 14 Supplementation of 100 mg berberine sulfate trihydrate to the diet of healthy postmenopausal women can promote bone formation and reduce bone turnover by reducing the concentration of serum osteocalcin (OC) and increasing serum 25 (OH) D levels. 15 By treatment with berberine (100 mg/kg) for 12 weeks, the levels of bone turnover markers ALP, TRAP, and OC could be restored by activating the AMPK pathway, and the degree of oxidative stress in rats could be reduced, thus improving, to some extent, bone loss in diabetic rats induced by streptozotocin plus high-fat diet. 16 A recent study showed that berberine (30 or 60 mg/kg) treatment for 2 months promoted bone marrow mesenchymal stem cell osteogenesis and inhibited lipid differentiation through the cAMP/PKA/CREB signaling pathway, alleviating the deterioration of bone structure in 20-month-old male C57BL/6J mice (senile osteoporosis mouse model). 17 The above studies mainly focused on the treatment efficacy of berberine, but the potential targets and underlying molecular mechanism of berberine in the treatment of osteoporosis remain unclear. Because several special targets and signaling pathways were found to be involved in the activity of berberine against osteoporosis by cell, animal, and human experiments, we hypothesized that berberine might exert its effects by acting on multiple potential targets and signal pathways related to osteoporosis.
Network pharmacology analysis can provide a network target and multicomponent approach to clarify the multi-target and multi-pathway characteristics of berberine in treating osteoporosis. In this study, we used network pharmacology, microarray analysis, data identification, molecular docking, and molecular dynamics simulation to screen the core target proteins of berberine in osteoporosis and identified the signaling pathways related to its activity. We also identified differentially expressed genes (DEGs) to enrich the potential targets and compared the efficacy of a clinical drug (raloxifene) with berberine at the molecular level. Further, we used classical Newtonian mechanics to simulate the movement of berberine–core protein complexes at the micro level to verify the stability and flexibility of their binding. Additionally, computational chemistry tools predicted the pharmacokinetics of berberine so as to guide its future application.
Results
Potential Targets of Berberine
From the Traditional Chinese Medicine Systems Pharmacology (TCMSP) and STITCH databases, 17 and 10 prediction targets were identified. A total of 15 prediction targets were obtained from the PharmMapper database. After combining the results from the above 3 databases and removing duplicates, we obtained 39 prediction targets, including NOS2, MAPK1, F10, and PTGS2 (Table 1).
Potential Targets of Berberine.
Potential Targets of Osteoporosis and Intersection Targets
We screened 450 DEGs (176 upregulated, 274 downregulated) from the GEO database (Figure 3A), which were considered as osteoporosis-associated targets; 1314 targets of osteoporosis were identified using Online Mendelian Inheritance in Man (OMIM), DrugBank, GeneCards, and Therapeutic Target Database (Figure 3B). After combining these DEGs and removing the duplicates, 1738 potential targets were obtained. The predicted targets included PTPN4, CALB1, PDCD6, and MAPK1. Venn diagrams revealed 23 intersecting targets (Figure 4A) between berberine and osteoporosis (Table 2).

Potential targets of osteoporosis. (A) Differentially expressed genes of osteoporosis. Blue and red nodes represent downregulated and upregulated genes, respectively. (B) Total targets of osteoporosis from 4 databases and DEG 56815.

(A) Total targets of berberine and osteoporosis, and 23 intersection targets of berberine against osteoporosis were identified using the Venn diagram and PPI network. The right chart is the PPI network of 23 intersection targets of berberine against osteoporosis. Nodes represent the intersection targets of berberine and osteoporosis. The node color and size represent the strength of the interaction. A stronger intersection between targets has a color closer to orange. In contrast, a node closer to the blue color indicates a weaker intersection between targets. Moreover, the line width and transparency represent the strength of the interactions. (B) Most significant module in the PPI network. Nodes represent hub genes. (C) Degree of hub genes.
Intersection Targets of Berberine Against Osteoporosis.
PPI Network Results and Molecular Analysis of Berberine Activity Against Osteoporosis
The PPI network had 23 nodes and 86 edges (Figure 4A). The action mechanism of berberine in treating osteoporosis was predicted to be related to AKT1, ESR1, AR, MAPK1, PTGS2, and other targets. We obtained a significant protein cluster (11 nodes, 48 edges) through the MCODE function (Figure 4B), indicating that these targets are involved in the same processes. The results of biological process analysis (Figure 5) showed that these hub targets might be associated with mammary gland development and the response to oxidative stress. The Kyoto encyclopedia of genes and genomes (KEGG) signal pathway analysis (Table 3) indicated that the core genes play a therapeutic role through the AGE-RAGE signaling pathway and TNF signaling pathway.

Gene ontology enrichment analysis of the most significant module with ranking.
KEGG Signal Pathway of hub Genes.
Pathway Enrichment Analysis to Explore the Therapeutic Mechanisms
To clarify the mechanism of berberine against osteoporosis, Gene Ontology (GO) and KEGG pathway analyses of 23 intersection targets were performed using the Metascape database. GO enrichment analyses revealed the biological mechanism from the following 3 perspectives: biological processes, cellular components, and molecular functions (Figure 6). Based on the results, berberine was predicted to be involved in the estrogen metabolism process, steroid hormone metabolism process, and oxidative stress process through binding to protein kinase, protein domain-specific binding, and nitric-oxide synthase regulator activity.

Gene ontology terms and signal enrichment analysis of candidate targets of berberine against osteoporosis.
KEGG pathway enrichment analysis revealed the major transmembrane signaling pathways involved in drug treatment of diseases. The analysis revealed 145 enriched KEGG pathways, which were related to thyroid hormone signaling pathway, estrogen signaling pathway, AGE-RAGE signaling pathway in diabetic complications, ovarian steroidogenesis, fluid shear stress, atherosclerosis, and bile secretion (Figure 6). Based on the number of targets involved in the pathways, a target-pathway Sankey plot (Figure 7) was constructed using the “ggalluvial” package in RStudio software. Most targets were involved in the estrogen signaling pathway, thyroid hormone signaling pathway, and AGE-RAGE signaling pathway in diabetic complications. Based on the KEGG enrichment analysis results, we predicted that the activity of berberine against osteoporosis, particularly postmenopausal osteoporosis and secondary osteoporosis caused by endocrine disease (hyperthyroidism), was mainly related to hormone secretion.

Main target-pathway Sankey plot of berberine for treating osteoporosis.
Molecular Docking Results of Hub Genes
Molecular docking was conducted to identify the interaction between berberine and its target proteins at the molecular level. Six intersection targets with the highest degree from the PPI network were selected for molecular docking analysis with berberine (Figure 8). The binding energies, binding sites, and modes of hub targets with berberine are shown in Table 4. The binding energy was less than 0 kcal/mol, indicating that the ligand and receptor can bind spontaneously. Moreover, lower binding energy indicated greater stability of binding of the active compound to the targets. The results showed that berberine could bind to these 6 core targets spontaneously.

Molecular docking between berberine and 6 intersection proteins of berberine and osteoporosis. (A) AKT1, (B) MAPK1, (C) ESR1, (D) AR, (E) TP53, (F) PTGS2, and (G) Virtual docking of Raloxifene and ESR1.
Molecular Docking Results.
The docking result of berberine and ESR1 was compared with the positive control compound raloxifene (Table 4). The binding energy of raloxifene to ESR1 was lower than that of berberine, and raloxifene targeted more binding sites, indicating that the therapeutic effect of berberine through binding with ESR1 was not as strong as that of raloxifene. Therefore, the clinical effect of berberine alone may not be significant, and it may be necessary to develop compounds with other TCMs to enhance the curative effects.
Molecular Dynamics Simulation Analysis of Hub Genes
We performed a molecular dynamics simulation to demonstrate further the stability of berberine binding to hub genes. As shown in Figure 9, the composite proteins of berberine and TP53, MAPK1, AR, and PTGS2 were simulated for as long as 20 ns, and berberine and AKT1 and ESR1 were simulated for as long as 50 ns. Based on the stability of the skeleton structure of the protein complex, the complex of berberine and AR tended to be stable after 1.5 ns with an average root mean square deviation (rmsd) of 0.18 nm. Berberine and PTGS2 and MAPK1 tended to be stable after 5 ns, with corresponding average rmsd values of 0.28 and 0.27 nm, respectively. Berberine and TP53 were particularly unstable initially and began to stabilize at 7 ns and 0.29 nm. Because of the large molecular structures of AKT1 and ESR1 proteins, the time required for stabilization was increased accordingly. AKT1 tended to be stable at 30 ns, with a stable rmsd of around 1.10 nm. ESR1 tended to be stable at 17 ns, with an average rmsd of around 1.7 nm. These results showed that the complex structures of berberine and core proteins had relatively stable conformations during molecular dynamics simulation.

Root mean square deviation (rmsd) plot of berberine-target complex during molecular dynamics simulation.
Absorption, Distribution, Metabolism, Excretion, Toxicology (ADMET) Analysis of Berberine
A brief analysis of the physicochemical properties and pharmacokinetics of berberine was performed using the SwissADME website to explore its possible beneficial effects. Figure 10 shows the bioavailability radar of berberine, which evaluates oral availability in terms of lipophilicity, size, polarity, solubility, flexibility, and saturation. The colored area in the figure represented the appropriate range of oral bioavailability. Based on the analysis, there were no out-of-range results, which indicates that berberine is expected to be orally bioavailable.

Physicochemical properties of berberine. The pink area represents the optimal range for lipophilicity, size, polarity, solubility, saturation, and flexibility. The 6 red dots represent the predicted values of berberine in these 6 aspects.
Discussion
Osteoporosis is a bone disorder characterized by decreased bone mineral density, low bone mass, and increased bone fragility. 18 Osteoporosis has become a major health problem for the aging population worldwide. 19 TCM has a long history of use of over 2000 years and is often considered an alternative pharmacotherapy for many intractable diseases because of its mild side effects. 20 Berberine, a key component of C. chinensis, has been shown to be effective in preventing bone loss. For example, a decade ago, Xu et al 21 discovered through animal mechanics experiments that berberine inhibits bone resorption and improves bone formation to prevent glucocorticoid-induced osteoporosis. Subsequent studies have shown that berberine may induce osteogenic differentiation of bone marrow mesenchymal stem cells through a typical Wnt signal pathway. 22 Moreover, berberine inhibited osteoclast differentiation and osteoclast marker gene expression induced by RANKL, proving that berberine can reduce bone loss in vivo.23,24 These in vitro and in vivo experiments demonstrated the therapeutic effect of berberine on bone remodeling exerted via the basic pathways, but the mechanism of action against osteoporosis has not been fully discussed. Computer-based methods such as data mining and bioinformatics can reveal potentially relevant mechanisms; thus, we used these computer-based methods to explore comprehensively the mechanism of action of berberine in treating osteoporosis.
Through network pharmacology, we selected the top 6 core genes as hub targets (ESR1, AKT1, MAPK1, AR, TP53, and PTGS2) predicted to play a crucial role in the effects of berberine against osteoporosis. Molecular dynamics simulation showed that berberine displayed suitable binding activity with the 6 key anti-osteoporosis targets. Thus, berberine may exert its effect on osteoporosis through these 6 targets.
Previous studies showed that berberine stimulates the differentiation of osteoblasts through activating the p38 MAPK pathway, which upregulates COX-2 (also known as PTGS2) expression by stimulating the activity of Runx2. 25 This finding is consistent with our results. AKT1 is a signaling intermediate that regulates osteoblasts and osteoclast differentiation. 26 Mukherjee suggested that AKT1 deficiency can promote osteoblast differentiation and inhibit osteoclastogenesis to increase further expression of the osteogenic transcription factor Runx2 and lead to net bone formation. 26 Moreover, androgens act on their target cells by binding to the AR, 27 and AR deficiency in male mice leads to increased bone resorption and decreased bone volume. 28
KEGG pathway analysis suggested that the thyroid hormone signaling pathway, estrogen signaling pathway, and AGE-RAGE signaling pathway in diabetic complications are the vital pathways associated with osteoporosis. The thyroid hormone signaling pathway plays an important role in bone formation. Studies have shown that thyroid hormone acts directly or indirectly on osteoblasts, osteocytes, chondrocytes, or other bone marrow cells 29 and regulates the signal pathways of many growth factors such as Wnt, insulin-like growth factor 1, PTH-related peptide, and fibroblast growth factor to regulate bone growth. 30 In the body, when thyroid hormone is lacking, the bone remodeling cycle is longer, and bone mass is increased. In contrast, in patients with hyperthyroidism, the bone remodeling cycle is shortened, and bone loss increases. 31 In the present study, the thyroid hormone signaling pathway was found to be enriched, which may be a direct or indirect effect of estrogen. Under the influence of estradiol, vital thyroid hormone receptors such as THα-1, THβ-1, and TSHR are highly expressed in different endometrial compartments. 32 Estrogen affects osteoporosis mainly via activated estrogen receptor alpha (Erα, encoded by ESR1). In the estrogen signal pathway, estrogen binds to ERα, then ERα translocates to the cell nucleus, binds DNA at the estrogen response elements, and activates the expression of estrogen response element-dependent genes, thus exerting an anti-osteoporosis effect. In osteoclasts, activated ERα positively regulates the expression of the FasL gene to cause apoptosis. 33 In osteoblasts, activated ERα enhances osteoblast activity by enhancing the paracrine signals of osteoblasts. 34 In addition, estrogen can inhibit the production of deubiquitinase Usp10, accelerating the degradation of senescence factor p53 (encoded by TP53), thus alleviating the senescence of osteoblasts. 35 RANKL induces the binding of c-Src and nuclear factor-κB receptor activator (RANK) to promote directly the bone resorption of osteoclasts 36 ; however, estrogen can form an ERα/SHP2/c-Src complex in osteoclasts, which attenuates c-Src activation upon RANKL stimulation to reduce its cytoskeletal recombination signal and osteoclast proliferation. 37 It is worth mentioning that estrogen deficiency can induce oxidative stress. 38 In the bone environment, the production of oxygen free radicals will lead to the up-regulation of RANKL and down-regulation of osteoprotegerin (OPG). 39 At the same time, oxidative stress can activate FoxO through JNK; FoxO enters the cell nucleus and activates β-catenin, resulting in bone formation deficiency. 40 In addition, estrogen can regulate the proliferation and apoptosis of osteoblasts and osteoclasts by interacting with the Wnt, β-catenin, and BMP-2 signaling pathways. 41 The KEGG analysis results indicated that the related signaling pathways have intricate interactions with other signaling pathways to a certain extent, showing that treatment of osteoporosis should not be limited to key pathways; instead, a multisystem intervention should be adopted by considering various factors such as disease causes, development status, and body environment.
This study had some limitations. Although we performed network pharmacology, microarray analysis, molecular docking, and molecular dynamics simulation, animal experiments and clinical trials are needed to verify the findings. Network pharmacology emphasizes the analysis of the molecular correlation between berberine and osteoporosis from the view of the system level and biological network as a whole, but it can only provide theoretical support; thus, in the future, we will explore the relationship between berberine and osteoporosis through in vivo and in vitro experiments, as well as clinical trials.
Although berberine has curative effects on osteoporosis of different causes (such as diabetic osteoporosis, senile osteoporosis, and postmenopausal osteoporosis), there are several differences in the potential targets and pathways of berberine regarding different osteoporosis types. When a more comprehensive discussion is carried out on a certain type of osteoporosis, corresponding conditions need to be added to the selection of disease targets, such as postmenopausal state, old age, or diabetes. In addition, we need to dig out the related DEGs. In the future, we can also use computer-related methods to study osteoporosis with specific causes and provide corresponding theoretical support for future biological experiments.
In conclusion, berberine can activate ESR1, AKT1, MAPK1, AR, PTGS2, and TP53, thereby regulating the thyroid hormone signaling pathway and estrogen signaling pathway for treating osteoporosis. The mechanism of action of berberine against osteoporosis is mainly associated with hormonopoiesis.
Conclusion
Through database target screening, DEG identification, PPI, GO, and KEGG analyses, molecular docking, and molecular dynamics simulation, we revealed a series of potential targets and signaling pathways which have not been widely examined in previous studies on the action of berberine against osteoporosis. We identified the key pathways related to the effects of berberine in osteoporosis: the estrogen signaling pathway and thyroid hormone signaling pathway. Our results provide a theoretical basis and reference for further studies of the anti-osteoporosis effects of berberine.
Materials and Methods
Screening of Potential Berberine Targets
Potential targets were obtained from the following 3 databases: TCMSP, 42 STITCH, 43 and PharmMapper.44–46 The details of the websites are shown in Table 5. Next, target names from the 3 databases were searched in the Universal Protein Database 47 for standardization, and the corresponding gene IDs were acquired.
Databases and Analyzing Platforms.
Screening of Potential Osteoporosis Targets
The following 4 databases were used to identify candidate targets of osteoporosis: OMIM, DrugBank, 48 GeneCards, 49 and Therapeutic Target Database. 50 The GSE56815 dataset obtained using Affymetrix HG-133A arrays was acquired from the NCBI-GEO database to analyze 40 Caucasian females with high and 40 with low hip bone mineral density characters of monocytes. Next, we applied the “limma” package in RStudio software to identify DEGs with an adjusted P-value of < .05 and |log2 fold-change| of > 0.1 as potential osteoporosis targets. The prediction targets were combined, and duplicates were removed.
Screening of Intersection Targets
The potential targets of berberine and osteoporosis were inputted into the Venn website to acquire the intersection targets.
Construction of PPI Network
The intersection targets of berberine in osteoporosis were applied to the STRING database. 51 We set the organism to “Homo sapiens,” selected hidden free points, and set the minimum required interaction score to 0.4. The data were imported into Cytoscape 3.8.0 software to construct a network of intersecting targets. The MCODE function in the Cytoscape 3.8.0 software was used for module analysis of the network, setting the K-Core value as > 4 and aiming to identify the key protein cluster and hub targets with high degrees. The key protein cluster was further inputted into the Metascape database, setting a P-value < .01 as the criterion for the statistical significance to perform the biological process and KEGG signal pathway analyses. Further, we sorted the intersection targets in the protein cluster by degree value, and those with a degree ≥ 12 were considered hub targets.
GO and KEGG Pathway Enrichment Analyses
The Metascape database 52 was used for GO and KEGG pathway enrichment analyses of the intersecting targets. The gene symbols were inputted into the Metascape database, with “Homo sapiens” selected as the organism, and a P-value < .01 was considered the cutoff criterion for statistical significance.
Molecular Docking
The 3D chemical structure of berberine was obtained from the ZINC database. 53 The crystal structures of the target proteins were acquired from the RCSB Protein Data Bank. Co-crystalized ligands and water from the targets were separated using PyMOL-2.1 software. Autodock-1.5.6 software was used to design the active pockets, perform molecular docking of berberine with the hub targets, and visualize the binding energies. To understand further the potential effect of berberine against osteoporosis, we used raloxifene as a positive control to compare its binding energy with the hub targets. According to the minimized binding energy, binding sites, and binding mode, the results obtained from Pymol-2.2 and Ligplot-2.2.4 were used for visualization and further analysis.
Molecular Dynamics Simulation
After docking, to verify the binding stability, we conducted molecular dynamics simulation for complexes of drugs and proteins using Gromacs software and the force field parameters with the ATB 54 website. The structures of core genes were obtained from the RCSB Protein Data Bank, and the required file types were constructed in GROMACS. The force field parameter of the small-molecule drug was obtained from the ATB website. To ensure that the system was electrically neutral, we added the number of sodium ions or chloride ions to the corresponding composite systems. The periodic boundary condition was set, and equilibrium phase treatment was carried out with energy minimization, NVT equilibrium, and NPT equilibrium. The energy minimization was generated using the steepest descent algorithm, and the maximum–minimization step was 50 000 steps. The NVT balance was set to 100 ps, with 2 fs at each step. The NPT balance was performed for a total of 100 ps. The systems entered the production phase, running for 20 to 50 ns to determine the stability of the combination.
SwissADME Analysis of Berberine
ADME properties are a key measure of a compound's ability to become a drug. To explore the possible beneficial effects of berberine, its physicochemical properties and pharmacokinetics were analyzed. SwissADME 55 is a commonly used web tool for performing ADME pharmacological analysis, allowing access to fast but robust predictive models of physicochemical properties, pharmacokinetics, and drug-like properties. The structural formula of berberine was input and analyzed in SwissADME.
Footnotes
Acknowledgments
We thank Yuchuan Traditional Chinese Medicine (Official bilibili id: 412608068) of Beijing University of Chinese Medicine for providing suggestions for this article.
Author Contributions
Jin-Ling Yu and Zhi-Hui Liu designed this study. Jin-Ling YU, Bo-Wei Wang, and Hui-Li Zhang conducted the experiments. Jin-Ling Yu and Liu-Qing Yang drafted the main manuscript. Jin-Ling Yu, Jing-Jing Yao and Han-Dan Huang revised the manuscript. Jin-Ling YU, Lu-TAO and Ying-GAO drew the picture in the manuscript. All the authors have read and approved the final manuscript.
Ethical Approval
Ethical approval is not applicable, because this article does not include any research on humans or animals.
Declaration of Conflicting Interests
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Key Scientific and Technological Research and Development Projects, Medical and Health Project (grant nos. 20180201056YY, 20190304032YY).
Informed Consent
Not applicable, because this article does not contain any studies with human or animal subjects.
Trial Registration
Not applicable, because this article does not contain any clinical trials.
