Network pharmacology of iridoid glycosides from Eucommia ulmoides Oliver against osteoporosis

Eucommia ulmoides Oliver is one of the commonly used traditional Chinese medicines for the treatment of osteoporosis, and iridoid glycosides are considered to be its active ingredients against osteoporosis. This study aims to clarify the chemical components and molecular mechanism of iridoid glycosides of Eucommia ulmoides Oliver in the treatment of osteoporosis by integrating network pharmacology and molecular simulations. The active iridoid glycosides and their potential targets were retrieved from text mining as well as Swiss Target Prediction, TargetNet database, and STITCH databases. At the same time, DisGeNET, GeneCards, and Therapeutic Target Database were used to search for the targets associated with osteoporosis. A protein–protein interaction network was built to analyze the interactions between targets. Then, DAVID bioinformatics resources and R 3.6.3 project were used to carry out Gene Ontology enrichment analysis and Kyoto Encyclopedia of Genes and Genomes pathway analysis. Moreover, interactions between active compounds and potential targets were investigated through molecular docking, molecular dynamic simulation, and binding free energy analysis. The results showed that a total of 12 iridoid glycosides were identified as the active iridoid glycosides of Eucommia ulmoides Oliver in the treatment of osteoporosis. Among them, aucubin, reptoside, geniposide and ajugoside were the core compounds. The enrichment analysis suggested iridoid glycosides of Eucommia ulmoides Oliver prevented osteoporosis mainly through PI3K-Akt signaling pathway, MAPK signaling pathway and Estrogen signaling pathway. Molecular docking results indicated that the 12 iridoid glycosides had good binding ability with 25 hub target proteins, which played a critical role in the treatment of osteoporosis. Molecular dynamic and molecular mechanics Poisson–Boltzmann surface area results revealed these compounds showed stable binding to the active sites of the target proteins during the simulations. In conclusion, our research demonstrated that iridoid glycosides of Eucommia ulmoides Oliver in the treatment of osteoporosis involved a multi-component, multi-target and multi-pathway mechanism, which provided new suggestions and theoretical support for treating osteoporosis.

Osteoporosis (OP) is a chronic epidemic characterized by low bone mass and deterioration of bone microarchitecture, which leads to increased bone fragility and fracture risk, and caused a heavy economic burden to society 1 . The etiology of OP is very complex, including the interaction of endocrine, nutritional status, genetic, physiological and immune factors 2,3 . The imbalance between bone formation of osteoblasts and bone resorption of osteoclasts is the underlying cause of OP 4 . The treatment of OP depends on drug therapy, including bisphosphonate, selective estrogen receptor modulator, mixed steroid receptor agonist, monoclonal antibody against RANKL, parathyroid hormone analogue and so on 5 . These drugs can alleviate bone loss and improve clinical symptoms to a certain extent, but their long-term clinical application is limited by low tolerance, severe side effects and high cost 6 . Therefore, it is of great significance to develop more safe, effective and economical drugs for the treatment of OP.
Traditional Chinese medicine (TCM) has a long history in China. It is more and more popular with the advantages of good curative effect, few side effects and affordable price 7 . In the theoretical system of TCM, OP is recognized as bone atrophy or arthralgia syndrome caused by kidney deciency 8 . Eucommia ulmoides Oliver (EU) is one of the most important nourishing medicinal materials in TCM. It has been found that EU can Molecular docking. Molecular docking was used to confirm the interactions between the compounds and the hub targets of IGEUs in treating OP and to verify the accuracy of the network pharmacology prediction. The 3D structures of the target proteins were downloaded from the RCSB database (https:// www. rcsb. org/), and the MOL2 structures of the compounds were downloaded from the TCMSP database (https:// tcmspw. com/ tcmsp. ph) 52,53 . SYBYL-X 2.1.1 was used to carry out the molecular docking of core compounds (ligands) with hub targets (receptors) and determine their binding activity 54 . The Total score is a complex score obtained by docking the receptor and ligand with corresponding parameters using Surflex-Dock procedure. It is generally believed that when the conformation of the ligand and receptor complex is stable, the higher the Total score, the higher the affinity of the receptor and ligand 55 . Total Score > 4.0 indicates certain binding activity, Total Score > 5.0 indicates good binding activity, while Total Score > 7.0 indicates strong binding activity 56 . PyMOL 2.4 (https:// pymol. org) software was used to visualize the docking results 57 . Molecular dynamic simulation and binding free energy analysis. The selected compounds were prepared with ATB webserver (http:// atb. uq. edu. au/) prior to MD simulation to generate an initial topology for the ligands 58 . MD simulation of the complex was performed using GROMACS 2019.6 59 . GROMOS96 54a7 force field was applied to the system, and dodecahedron water box consisting of a TIP3P water model was used to solvate the system. Na + and Cl − ions were also solvated in the box to neutralize the system charge. Next, energy optimization process was performed using the steepest descent method. For equilibration simulation, 100 ps NVT equilibration was performed by using the velocity rescaling thermostat coupling method to keep the temperature constant at 300 K. Then, 100 ps NPT equilibration was performed. To treat the longrange coulombic interactions, the PME method was used. Production run of MD simulation was performed till 50 ns for each protein-ligand complex. Consequently, root mean square deviation (RMSD), root-mean-square fluctuation (RMSF), the radius of gyration (Rg), solvent-accessible surface area (SASA) and hydrogen bonds were calculated according to the trajectory for further analysis. The binding free energies (ΔG bind ) including electrostatic interactions (ΔE elec ), Vander Waals interactions (ΔE vdW ), non-polar solvation energy (ΔG SASA ) and polar solvation energy (ΔG polar) were calculated using the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) method implemented in GROMACS compatible tool "g_mmpbsa" 60,61 . In addition, the energy of each residue was decomposed, and the energy decomposition could be analyzed to determine the contribution of key residues to binding.

Results
Collection and screening of active IGEUs and construction of networks. After text mining and removal of non-target compounds reported in the literatures, 12 iridoid glycosides were obtained as active IGEUs (Table 1). In total, 161 targets of IGEUs were identified. The compound-target network is composed of 175 nodes and 554 edges ( Fig. 2A). By searching the disease-related database, a total of 4124 targets of OP were integrated. Finally, 97 common targets were identified as therapeutic targets for the anti-OP activity of IGEUs (Fig. 2B). Moreover, a compound-target-disease network with 111 nodes and 414 edges was constructed using IGEUs, osteoporosis and common targets (Fig. 2C).
PPI network of the anti-OP targets of IGEUs. The 97 common targets were imported into the STRING11.0 database to construct the PPI network. After removing 4 disconnected nodes, there were 93 nodes left, but the original PPI network was usually rough (Fig. 3A). Therefore, a second PPI network was constructed by Cytoscape 3.8.0 in order to obtain a better visualization and understanding. The results showed that the  (Fig. 4). The results of GO functional enrichment analysis showed that IGEUs in the treatment of OP were mainly regulated by response to oxygen-containing compound, regulation of cell proliferation and apoptotic, response to an organic substance, regulation of protein phosphorylation in BP, and mitochondrion, nucleus, cytoplasm and nucleoplasm in CC, and protein binding, enzyme regulator activity, protein phosphatase binding and protein phosphatase binding in FM. The above analysis suggested that the active compounds of IGEUs may exert anti-OP effects by participating in various biological regulatory processes. The KEGG enrichment analysis revealed 105 pathway items, and the top 10 pathways were shown in Fig. 5A. The detailed information of targets and pathways were listed in Table 4. The main pathway included pathways in    [62][63][64] . Therefore, the above three signaling pathways and related targets were considered as the candidate pathways for further validation. The compound-hub target-pathway network, including 46 nodes and 173 edges (Fig. 5B). The main active compounds of IGEUs were distributed in different pathways and played a coordinating role in the treatment of OP. The core compounds included aucubin (edgecount = 9), reptoside (edgecount = 9), geniposide (edgecount = 6) and ajugoside (edgecount = 6).

Molecular docking.
We screened 12 IGEUs and 25 hub targets for molecular docking verification, which played a more critical role in the treatment of OP. The results showed that the Total score between most target proteins and compounds were above 4, indicating that these ligands and receptors could bind stably (Fig. 6). Moreover, we selected 4 hub targets as representative examples to show their docking modes, namely AKT1(PDB: 3OS5), ESR1(PDB ID: 6PFM), MAPK1(PDB ID: 5K4I) and MAPK3(PDB ID: 2ZOQ). The binding mode of reptoside with AKT1 showed that reptoside formed 5 hydrogen bonds with Asn296, Tyr337, Glu356, Ser224 and Ser225, formed hydrophobic interaction with the Gly227, Asn265 and Tyr353 (Fig. 7A). The docking mode of aucubin and ESR1 had the highest Total score. We could observe the formation of 6 hydrogen bonds at www.nature.com/scientificreports/ different sites of ESR1, namely aucubin with Gly521, His524, Leu346, Arg394 and Leu387. We could also observe the hydrophobic interaction between aucubin with Trp383, Leu349 and Met343 (Fig. 7B). The binding between geniposide and MAPK1 included 4 hydrogen bonds with Lys151, Arg67, Tyr36 and Gly37, and hydrophobic interaction linked with Glu33, Ser153 and Asp167 (Fig. 7C). The ajugoside at the active site of MAPK3 formed 5 hydrogen bonds interactions with Gly49, Gly50, Asp184, Asn171 and Ser170. At the same time, ajugoside formed hydrophobic interactions with Ala69, Leu124, Met125 and Asp128 (Fig. 7D). www.nature.com/scientificreports/ Molecular dynamic simulation and binding free energy analysis. We selected the above 4 representative molecular docking complexes for further MD study. The RMSD analysis of MD trajectories showed that all systems reached equilibrium at 50 ns with minimal fluctuations (Fig. 8A). The plot of AKT1-reptoside complex showed a stable equilibrium approximately at 30 ns, then a stable RMSD could be seen approximately at ~ 0.27 nm. The ESR1-aucubin complex was stable from the beginning and had a consistent RMSD around ~ 0.2 nm. The MAPK1-geniposide complex stabilized at ~ 0.22 nm during 50 ns of MD simulation. Compared with other complexes, the MAPK3-ajugoside took a relatively slower time to reach a stable conformation. In summary, the stability of complexes was high, and the trajectories were suitable for further analysis.
In order to understand the structural stability of the protein-ligand complexes, we determined the compactness of the protein structure by computing the Rg (Fig. 8B). The Rg plots showed that the structural dynamics of AKT1-reptoside, ESR1-aucubin, MAPK1-geniposide and MAPK3-ajugoside complexes were quite stable throughout the simulation time, with mean Rg values of ~ 1.83 nm, ~ 1.64 nm, ~ 1.50 nm and ~ 2.19 nm, respectively. AKT1-reptoside complex had the highest Rg towards the whole simulation, and Rg for MAPK1-geniposide and MAPK3-ajugoside complexes were almost similar.
Moreover, the SASA of the complexes was analyzed to assess the complex volume change through the simulation trajectories (Fig. 8C). The complexes of AKT1-reptoside and ESR1-aucubin had a higher SASA profile in the  www.nature.com/scientificreports/ initial phase, followed by a lower descriptor, and then maintained a stable SASA until the end of the simulation. The SASA of MAPK1-genpiside and MAPK3-reptoside complexes were similar, with less fluctuation during the simulation, indicating that the stability of these two complexes were less affected by the solvent. Hydrogen bond interaction is one of the main parameters reflecting the stability of the ligand at the active site of the protein. Thus, we further investigated the changes in the numbers of hydrogen bonds (Fig. 8D). We found there was a maximum occupancy of 6 hydrogen bonds between reptoside and AKT1. Among them, 4 hydrogen bonds remained consistent until ~ 50 ns. There were up to 7 hydrogen bonds between aucubin and ESR1, which could be seen at ~ 20 ns, and 3 hydrogen bonds remained stable in the last 10 ns. MAPK1 formed up to 6 hydrogen bonds with geniposide, which were observed consistent until 50 ns. MAPK3 showed the possibility of forming up to 7 hydrogen bonds with ajugoside. The above analysis results revealed that the hydrogen bond interactions between these amino acid residues and compound were dynamic.
To better understand the molecular interaction and stability related to the complexes, the binding free energy was analyzed in detail. Results showed that all the binding free energy was less than zero, indicating the reaction can proceed spontaneously ( Table 5). The detailed decomposition of the energy components revealed that the van der Waals energy and electrostatic interaction energies were the major contributors to the binding free energy of the complexes. The nonpolar solvation energy played a supplement role in binding. To analyze the contribution of residues to protein ligand interaction, the free energy decomposition per residue was employed (Fig. 9). Residues with energy > 5.0 kJ/mol or < − 5.0 kJ/mol were considered to be the critical residues for ligand binding to protein 65 . The calculation results showed that Trp352 and Tyr335 in AKT1 had strong interactions with reptoside (Fig. 9A). Aucubin had the lowest interaction energy with Leu346. In addition, the binding energy of aucubin with Leu387, His524 and Leu525 were also low (Fig. 9B). The binding of geniposide to MAPK1 was mainly supported by the amino acids' residues Leu156, Val39, Ile31, Lys54 and Asp111 (Fig. 9C). Analysis of MAPK3-ajugoside complex showed that Leu346, Leu387, His524, Leu525 and Glu353 energetically favor the binding of ajugoside (Fig. 9D). Overall, the identification of critical residues in these proteins facilitated the discovery of new selective inhibitors against OP-related targets.

Discussion
OP is becoming a major health problem with increasing age and aging bones, placing a heavy economic burden on society and families 66 . The etiology and pathogenesis of OP are still unclear. Therefore, single-target drugs cannot fundamentally prevent the development of OP 67 . TCM has the characteristics of high safety and few side effects, and it has unique advantages for complex diseases 68 . EU is the top grade of TCM, which has the effect of strengthening muscles and bones, nourishing liver and kidney 69 . Numerous studies have reported the beneficial effects of EU on skeletal and renal diseases 69,70 . But research on IGEUs is very limited. Therefore, this study was the first to explore the mechanism of IGEUs in the treatment of OP through network pharmacology and molecular simulations.
In this study, we identified effective compounds, target proteins and important pathways for IGEUs in the treatment of OP. Network analysis showed that aucubin, geniposide, reptoside and ajugoside were the core compounds of IGEUs against OP. Aucubin had strong anti-OP activity, which can not only increased the differentiation of osteoblasts, promote the increase of cortical bone thickness and bone density, but also prevented www.nature.com/scientificreports/ the apoptosis of osteoclasts 71,72 . Moreover, we have also reported that the total glycosides in Eucommia ulmoides seeds contain high content of aucubin, which could enhance bone mineral density and bone strength, suggesting that it may be a potential alternative drug for the treatment of OP 73 . Geniposide significantly promoted the formation of calcified nodules and induced osteogenic differentiation 74 . Ajugoside could resist oxidative damage in some tissues by increasing the activity of SOD 75 . Reptoside exerted good anti-inflammatory activity through inhibiting COX-1 and COX-2 enzymes 76 . Through the construction of PPI network, we identified 25 hub targets of IGEUs in the treatment of OP. These hub targets showed rich interactions with other target proteins and were also involved in 105 biological pathways. In order to find key biological pathways, we analyzed the pathways with more annotation targets and lower p-value. Importantly, PI3K-Akt signaling pathway, MAPK signaling pathway and Estrogen signaling pathway were related to the OP system, which were also consistent with the results reported before [62][63][64] . PI3K-Akt signaling pathway is an important signaling pathway for increasing osteoblast   77 . AKT1 is a crucial signaling molecule in this signaling pathway and plays a significant role in osteogenesis 78,79 . Deficiency of AKT1 results in decreased bone mineral density throughout the body and in the femur, which is reflected in decreased femoral fracture resistance 80,81 . Many studies have shown that osteoclast generation and function can be inhibited by inhibiting the MAPK signaling pathway 82,83 . MAPKs, including p38, JNK and ERK, play an intermediary role in the regulation of bone formation. The p38, ERK1 (MAPK3) and ERK2 (MAPK1) promoted osteoblast differentiation via phosphorylation of Runx2 84 . Estrogen signaling pathway is essential for the development of OP. Estrogen is an important regulatory hormone in the human body, and its physiological role is mainly to regulate the transcription and translation of target genes by acting on the estrogen receptors of tissue cells 85 . As the major estrogen receptor subtype in bone tissue, ESR1 also plays an important role in regulating bone metabolism. Estrogen deficiency is one of the main causes of postmenopausal OP, which indicates the role of ESR1 in human bone homeostasis 86 . Studies have shown that the ratio of RANKL/OPG increased after estrogen cessation, which further leading to the increase of bone resorption 87 . Then, molecular docking technology was used to verify the binding of hub targets with core compounds. Results showed that most hub targets had a certain binding activity, and AKT1, ESR1, MAPK1 and MAPK3 had a strong binding activity, indicating they could tightly integrate with corresponding compounds. We further study the interactions between proteins and ligands using MD simulation, which can present the fluctuation and movement of residues at any specific time in motion 88 . The RMSD showed that all complexes reached equilibrium till 50 ns, and the RMSF analysis proved the flexibility of the active amino acid residues, which facilitated the ligand binding. Similarly, the Rg and SASA analysis showed that the complexes were less affected by the solvent. Furthermore, the hydrogen bond assessment results were consistent with the results of molecular docking. The binding free energy analysis revealed that all core IGEUs showed stable binding at the binding pocket of the target proteins during the simulation time.
Through network pharmacology, molecular docking and MD simulation analysis in this study, we believed that the active compounds of IGEUs could affect the differentiation and survival of osteoclasts and osteoblast, regulate the balance of bone resorption and osteogenesis, and achieve the purpose of treating OP.

Conclusion
The present study combined network pharmacology, molecular docking and molecular dynamic methods for the first time to reveal the pharmacological mechanism of IGEUs in the treatment of OP. Our data demonstrated that IGEUs could treat OP through a multi-component, multi-target and multi-pathway mechanism. In conclusion, this valuable finding may provide theoretical support for the clinical application of IGEUs and molecular design of therapeutic targets for OP.