Abstract
Ruoshi Yuan1,2,*, Shengfei Ma3,*, Xiaomei Zhu5, Jun Li4, Yuhong Liang4, Tao Liu4, Yanxia Zhu4, Bingbing Zhang4, Shuang Tan4, Huajie Guo4, Shuguang Guan3, Ping Ao1,2 and Guangqian Zhou4
1 Key Laboratory of Systems Biomedicine, Ministry of Education, Shanghai Center for Systems Biomedicine, Shanghai Jiao Tong University, Shanghai, China
2 School of Biomedical Engineering, Shanghai Jiao Tong University, Shanghai, China
3 Department of Physics, East China Normal University, Shanghai, China
4 The Center for Anti-Ageing and Regenerative Medicine, Shenzhen University, Shenzhen, China
5 GeneMath, N.E., Seattle, WA, USA
* These authors have contributed equally to this work
Correspondence to:
Ping Ao, email:
Guangqian Zhou, email:
Keywords: hormone, osteoporosis, regulatory network, strontium treatment, systems biology, Gerotarget
Received: August 30, 2015 Accepted: January 04, 2016 Published: January 15, 2016
Abstract
To develop and evaluate the long-term prophylactic treatment for chronic diseases such as osteoporosis requires a clear view of mechanism at the molecular and systems level. While molecular signaling pathway studies for osteoporosis are extensive, a unifying mechanism is missing. In this work, we provide experimental and systems-biology evidences that a tightly connected top-level regulatory network may exist, which governs the normal and osteoporotic phenotypes of osteoblast. Specifically, we constructed a hub-like interaction network from well-documented cross-talks among estrogens, glucocorticoids, retinoic acids, peroxisome proliferator-activated receptor, vitamin D receptor and calcium-signaling pathways. The network was verified with transmission electron microscopy and gene expression profiling for bone tissues of ovariectomized (OVX) rats before and after strontium gluconate (GluSr) treatment. Based on both the network structure and the experimental data, the dynamical modeling predicts calcium and glucocorticoids signaling pathways as targets for GluSr treatment. Modeling results further reveal that in the context of missing estrogen signaling, the GluSr treated state may be an outcome that is closest to the healthy state.
Introduction
Osteoporosis is characterized by a decrease in bone mass and density, leading to increased risk of fracture [1, 2]. Primary osteoporosis is aging related [3], while secondary osteoporosis is loss of bone caused by an agent or disease process such as corticosteroids, endocrine disorders, or an inflammatory process and malignant diseases [4, 5]. In addition to the apparent differential disease progress, osteoporosis is a multiple cell type disease including bone-forming osteoblasts and bone-resorbing osteoclasts that determine the bone remodeling balance [6-11]. Osteoporosis as a multi-caused disease at molecular level is also demonstrated by animal models [12]. The ovariectomized rat model is commonly used for research on postmenopausal osteoporosis [13]. Glucocorticoids feeding is another way to induce osteoporosis in animal models [14]. Mice receiving glucocorticoids for 7 days showed an early increase in bone resorption and exhibited decreased bone mineral density at week 4 together with changes in numbers of osteoblasts and osteoclasts; progenitors in the bone marrow; poor mineral depositional rate; and a dramatic reduction in cancellous bone mass [14]. Osteoporosis can also be induced by retinoic acids [15]. These animal models highlights critical molecular pathways for the disease, implicating a core structure for the underlying molecular mechanism, owing to the facts that hormones usually regulate a wide range of biological processes. The most common prophylactic methods are thus far limited to calcium and vitamin D supplements. Although a number of potential preventative and treatment targets have been investigated, the long safety is difficult to estimate. Strontium supplement has been used for a decade so far, without established relation to serious side effects [16].
The current method for osteoporosis study is often limited to assigning pro- or anti- osteoporosis role for each risk factor, which may lead to contradictions due to inability to account for multiple or even opposing effects simultaneously [17-19]. Estrogens were assigned anti-osteoporosis for their anti-inflammation function [20-23]. Glucocorticoids induced osteoporosis was suggested as a result of decreased bone formation, impaired osteoblast differentiation and increased osteocytic apoptosis [14, 24-26]. However, glucocorticoids are also anti-inflammation agents, which may reduce bone absorption. We found that estrogen loss in OVX rats lead to abnormal osteoblast phenotype, which was reversed by GluSr treatment, indicating complexity of the disease. We carried out osteoporosis study beyond individual molecular functions to the network level. It is a necessary step because risk factors and molecular pathways do not act independently. When combined with dynamical properties, a network has been shown to be able to evaluate risk factors and their contribution to pathogenesis [27]. In this particular case, the effect of a drug on the network level is an indication for long-term safety.
Our proposal has two key biological components. First, instead of focusing on quantities of osteocytes versus osteoclasts, we propose that loss of function of osteoblasts, which determine mineralization, is important for this disease. The recent experiment on the Senescence-Accelerated Mouse-Prone 6 (SAMP6) showed an impaired and likely depressed production of osteoclast in the animal [19], providing at least one example of osteoporosis not caused by overproduction of bone destructing osteoclasts. SAMP6 mouse stroma showed significant gene expression of osteoblast but poor mineralization of bone matrix [19]. We conducted our own experiments on OVX rats. The results also demonstrate clear abnormality of osteoblasts, as shown in Figure 2, exhibiting patterns somewhat similar to that of the SAMP6 mouse model [28]. These experiments suggest an important role of osteoblast in the disease.
Second, we propose to understand molecular mechanism of osteoblast function and its contribution to the disease through hierarchical network approach based on the facts that hormones and other signaling molecules do not work separately. Estrogens, glucocorticoids, retinoic acids all exert their regulatory effect through the overlapping molecular machinery of signal transduction, molecular transcription and translation. It is inevitable that their pathways cross-talk. For example, previous works suggested that glucocorticoids could attenuate estrogen response [29, 30] while estrogen can prevent glucocorticoid-induced apoptosis in osteoblast [25]. In accordance, we constructed a tightly linked core network including several pathways critical for bone formation using previous biological and biochemical studies reported in literature. In this network, hormones are assumed to be at the top of hierarchical structure of molecular regulations and their signaling pathways treated as interacting hubs at the network level.
Results
The network is presented in Figure 1 and details shown in Supplementary Materials. The summary of the nodes and a brief description of their interactions are the followings.
Figure 1: Schematic diagrams of the proposed molecular mechanism for osteoporosis due to osteoblast function loss. a. Molecular signaling pathways included in the core regularity network. b. Examples for cross-talks between pathways. Left: Ligated RARα was found to co-immunoprecipitate with ligated GR and to physically interact directly and proposed to be a mechanism for their mutual transcriptional repression. In addition, in vitro studies confirmed that RA abolished GR-mediated glucocorticoid-induced suppression of CRH expression, indicating a negative cross-talk between RARα and GR signaling pathways. Middle: At enzyme level, estrogen sulfotransferase, an enzyme important for the metabolic deactivation of estrogens, is a transcriptional target of GR. Right: Calcium regulates glucocorticoid receptor in mouse corticotrope tumor cells by reversible conversion of the receptor to a non-binding form. See Supplementary Materials for more details. d. The core molecular network. Transcriptional upregulation/activation by mechanisms such as phosphorylation is represented by activation with green line and arrow. Transcriptional downregulation/deactivation is represented by inhibition with red line and dot. f. Calculated attractors for the core network using Boolean method. c. and e. Calculated attractive basins for the attractors corresponding OVX rats, OVX rats treated with GluSr and control. The normal osteoblast phenotype associates with attractor (4), with ESR(1), GR(0), RAR(0), PPAR(1), VDR(1), GRP120(1), Ca2+ (1), NCOA(1), according to microarray data shown in Figure 3. The overactive glucocorticoid pathway in OVX is apparent from Figure 3. The OVX osteoporotic osteoblast phenotype associates with attractor (14), with ESR(0), GR(1), RAR(0), PPAR(1), VDR(0), GRP120(0), Ca2+ (0), NCOA(1). The GluSr treated phenotype associates with attractor (5), ESR(0), GR(0), RAR(1), PPAR(1), VDR(1), GRP130(1), Ca2+ (1), NCOA(1). g. Calculated transition stages after GluSr treatment, representing shortest path of switching from OVX attractor (14) to treated attractor (5). Biologically, Ca2+ activation and GR signaling pathway inhibition are targets for GluSr. h. Calculated transition stages after ovariectomization. Biologically, it presents loss of estrogen signaling and subsequent GR signaling pathway activation. Note that the colors in the barcode of a state used here denote different activated signaling pathways, not the difference of their activity levels.
Glucocorticoid receptor signaling pathway (GR)
Glucocorticoids are a class of steroid hormones that bind to the glucocorticoid receptors presented in almost every vertebrate animal cell. They inhibit proliferation and differentiation of osteoblast. Glucocorticoids also regulate ion channels.
Retinoic acids signaling pathway (RAR)
Retinoic acids act by binding to retinoic acid receptors (RARs), which are members of the nuclear receptor superfamily. The RAR family comprising 3 isotypes: RARα, RARβ, RARγ. RARs act in heterodimeric combinations with retinoid X receptors (RXRs). Retinoic acid is essential for physiological regulation of a wide range of biological processes including development, differentiation, proliferation, and apoptosis.
Estrogen receptor signaling pathway (ESR)
Estrogen is required for bone growth and development and for the maintenance of bone health. The cellular responses of osteoblasts and osteoclasts to estrogen are initiated via two high-affinity receptors ERα and ERβ.
The peroxisome proliferator-activated receptors signaling pathway (PPAR)
PPARs are members of nuclear receptor family that can form a heterodimeric complex with RXR and function as transcription factors regulating the expression of genes. PPARs are implicated in major metabolic and inflammatory processes, in the control of cell proliferation, differentiation and survival.
Vitamin D receptor signaling pathway (VDR)
Vitamin D is a prohormone with a key role in calcium and phosphate balance and bone structure. The primary molecular action of 1,25(OH)2D is to initiate or suppress gene transcription by binding to the vitamin D receptor (VDR), which belongs to the family of trans-acting transcriptional regulatory factors similar to the steroid and thyroid hormone receptors.
Ca2+ signaling pathway (Ca2+)
Calcium signaling is an important cellular process [31]. For example, Orai2 is a plasma membrane protein forming Calcium Release-Activated Channels (CRACs) which are specialized plasma membrane Ca2+ ion channels [32, 33]. Orai2 increased significantly in GluSr-treated rats, and reduced in OVX rats.
Nuclear receptor coactivators (NCOA)
Nuclear receptor coactivators directly bind nuclear receptors and stimulate the transcriptional activities in a hormone-dependent fashion. They introduce additional complexity into the network.
Figure 2: Ultrastructure of osteoblasts, microCT features, and the bone mineral density of the tibia of OVX rats and the OVX rats treated with strontium gluconate in comparison with healthy controls. The ultrastructural features of OVX rat bone tissues a. under transmission electron microscopy (TEM) show similar pattern as that of SAMP6 mouse model. In contrast to the abundant cisternae of the granular endoplasmic reticulum and well developed Golgi complex found in control rats b. the OVX rats have strikingly fewer organelles and rounder shape, without numerous long cytoplasmic projections. Myelin-like structures or amorphous-like materials were sometimes observed. It might be further surmised that the cell bodies in OVX rats are swollen and contain large amorphous vacuoles that push the decreased sub-cellular organelles to the nuclei. The osteocytes in control rats b. were encompassed by mineralized matrix in contrast to those from OVX rats by amorphous materials without collagenous fibrils a. c. Interestingly, osteocytes from OVX rats treated with strontium (GluSr), an element that has been used as for its anti-osteoporosis effects, appears close to control rats. d. Representative 3D images of trabecular bone microstructure in the proximal tibia of controls, OVX rats and OVX rats treated with vehicle or GluSr. e. Changes of body weight and bone mineral density with time of controls, OVX rats and OVX rats treated.
Omega 3 fatty acids signaling pathway (GPR120)
We separated omega 3 fatty acids signaling from fatty acids signaling through PPARs signaling for its potential role in reduction of oxidative stress.
Estrogen and glucocorticoid pathway interact at multiple cellular levels. At enzyme level, estrogen sulfotransferase, an enzyme important for the metabolic deactivation of estrogens, is a transcriptional target of GR. Overall, GR and ESR attenuate each other’s effects. ESR expression can alter the transcriptional regulation of PPARγ target gene expression. The interaction between retinoic acid signaling and estrogen signaling is complex and perhaps context dependent. RARα in the presence of its ligand can antagonize estrogen-ER function, and vice versa, because RARα and ESR can, in some cases, share common cis-regulatory elements. We adopted their anatomization tentatively. Estrogen upregulates the expression of VDR and increases the responsiveness to 1,25(OH)2D. The binding rate of the glucocorticoid receptor is higher in vitamin A-deficient rats than in controls and restored by retinoic acid supplementation. Calcium downregulates glucocorticoid receptor in mouse corticotrope tumor cells not due to a decrease in GR protein but reversible conversion of the receptor to a non-binding form. The details of molecular interactions obtained from literature are presented in Supplementary Materials.
Since the network was entirely constructed by using well-documented molecular pathway interactions, the microarray data in this study serve as independent test and validation of the network. We found that the overall profiling is consistent with the prediction of this network. Interestingly, gene expression profiling demonstrated that among abundant differentially expressed genes encoding enzymes, shown in Figure 3, many of which appear to be glucocorticoids regulated or related. St3gal3, a sialyltransferase, is an acute phase reactant that is transcriptional regulated by glucocorticoids. Sulfotransferase, Sult1c3, and N-acetylglucosaminidase, Hyal4, are glucocorticoids inducible genes. O-GlcNAc transferase, Aer61, may participate in glucocorticoids signaling by associating with ligand bound glucocorticoid receptor in a multi-protein repression complex, mediating glucocorticoid-induced apoptosis. Tktl2 is an enzyme in pentose phosphate pathway whose flux associates with interconversion of biologically inactive 11-keto derivatives to active glucocorticoids. Pcsk2 is a prohormone convertase that is regulated by both glucocorticoids and estrogens. Plcb2, a gene encodes phospholipase C that mediates activation by both estrogens and glucocorticoids. CYP3A family of enzymes is responsible for degradation of glucocorticoids. UDP-glucuronosyltransferase and Pck1 that encodes phosphoenolpyruvate carboxylase are inducible by glucocorticoids. The genes relevant to the network are summarized in Figure 3(a). References are in Supplementary Materials.
The modeling of the network is performed by Boolean method. Among 256 total combinations of the On (1), Off (0) states of the nodes, we identified 19 attractors with balanced interaction and locally stable. The attractors are shown in Figure 1. The normal osteoblast phenotype associates with attractor (4), with ESR(1), GR(0), RAR(0), PPAR(1), VDR(1), GRP120(1), Ca2+(1), NCOA(1), according to microarray data shown in Figure 3. The OVX osteoporotic osteoblast phenotype associates with attractor (14), with ESR(0), GR(1), RAR(0), PPAR(1), VDR(0), GRP120(0), Ca2+(0), NCOA(1). The GluSr treated phenotype associates with attractor (5), ESR(0), GR(0), RAR(1), PPAR(1), VDR(1), GRP130(1), Ca2+(1), NCOA(1). The overactive glucocorticoid pathway in OVX is apparent from Figure 3. The fatty acids related genes, especially those directly participate in fatty acids metabolism such as Baat and Hmgcs1, are high in OVX rats. Overall, PPAR pathway is less active in osteoblasts of normal rats. In normal osteoblast cells, genes encode enzymes involved in the synthesis of retinoic acid Adh4 and Aldh1a1 are highly expressed. Meanwhile Crabp2, which encodes a Cellular Retinoic Acid-Binding Protein (CBABP), responsible for retinoic acids sequestering is also high. Rarg, which encodes RARγ, is low in comparison to OVX rat. Thus retinoic acid signaling in normal osteoblast cells does not appear to be overactive but rather well buffered. This buffering is partially intact in osteoblast cells in OVX rats but lost in GluSr treated rats. We used buffering and Rarg as criteria to assign values to retinoic pathway.
Using these phenotypes as starting point, by performing Boolean analysis, we predicted that Ca2+ and glucocorticoid receptors are key targets for GluSr treatment, as shown in Figure 1. As glucocorticoid receptor is also a target of Ca2+, equivalently, we may interpret that Ca2+ and its downstream effectors as direct GluSr targets, although GluSr influence propagates to systems level through network interactions as indicated by microarray data. Since GluSr mechanism is regarded as unclear [34, 35], our prediction from the core network suggests a viable direction for future investigation. The attractors of the network already determine that if estrogen signaling is missing, the GluSr treated state is already the closest outcome to normal state (in Figure 1). The treated state is not totally normal because of the difference in RAR signaling.
Figure 3: a. Differentially expressed genes in microarray experiment for osteoblasts in controls (blue), OVX rats (red) and OVX rats treated with GluSr (green) for relevant genes regulating or targeted by the pathways in the core network. Crabp2, Pck1, S100g are targets for more than one signaling pathways. b. List of topmost differentially expressed genes with the largest fold changes and those in a..
Discussion
Proper mineralization of extracellular matrix has at least three requirements: adequate production of matrix vesicles [36], calcium binding proteins in the vesicles, and calcium signaling that regulates vesicle formation [37-39]. Abnormal calcium binding protein expression and calcium signaling are factors that may lead to poor mineralization. Vitamin D signaling is responsible for the transcription of calcium binding proteins. Calcium signaling pathway is regulated by the hormonal network. Inflammation may be an indirect result of the modified peroxisome proliferator-activated receptor signaling and omega-3 fatty acids signaling, possibly contributing to the increase of matured osteoclast [40]. For the OVX rats, these factors interplay, as shown by the gene expression profile in Figure 3, but are corrected after GluSr treatment. Since glucocorticoids increase the death of osteocytes [41], likely through GR signaling pathway, our proposal does not conflict with the current understanding on TGF-β, LPR5, and PPARγ signaling to increase the number of osteoblasts, or to decrease the number of osteoclasts [6, 42-44]. Instead, this core network approach provides additional targets for treatment to improve the quality of osteoblast, as well as enables the evaluation of long-term risk. In the GluSr case, the activation of RAR pathway after treatment may be a potential side effect.
In perspective, osteoporosis is a pathological condition of lost bone density, which can be induced by a range of conditions. The tightly cross-linked core molecular network structure uncovered in this report might serve as a unified molecular mechanism for different phenotypes of osteoblast. At this stage, the network is not meant to be extensive but to initiate study on hormones and hormone-like molecules at the systems and network level. With this network, we predict that GluSr treatment can give an optimized outcome when estrogen is deficient, although the outcome may still not be ideal owing to alternation in RAR signaling [45]. Other treatments for osteoporosis, such as stem cell therapy [46], which might be restricted by the uncertainty of the transplanted stem cell fate and related side effects [47], may be studied and improved by the network approach. It may further be employed to explain complex biological phenomena at systems level such as developmental processes [48], aging [49], and cancers [50].
Materials and Methods
Animal procedures
All experimental procedures were approved by the Animal Care and Use Committee of Medical School, Shenzhen University. A total of 36 female Sprague-Dawley rats (Guangdong Medical Laboratory Animal Center, China), aged between 7-8 weeks, were induced for postmenopausal osteoporosis by bilaterally ovariectomized under intraperitoneal Pentobarbital Sodium (Sigma, USA) anesthesia (40 mg/kg). An additional 6 rats were left intact. All the animals were housed in groups of six per cage with free access to food and water, kept in Animal Hall with an average temperature of 22°C and under a 12hr light/dark cycle, with lights on at 7:00 A.M. Ovariectomized (OVX) animals were randomly divided into 2 groups of 6 rats each. Totally there were 3 groups in this study: 1) the sham-operated group, 2) the OVX group, and 3) OVX treated with Strontium Gluconate (GluSr). Oral treatments were administered daily on the 4th week after ovariectomy and lasted for 20 weeks. Within 20 weeks, animals were weighed monthly. Bone mineral density (BMD) of each animal was monitored monthly by an in vivo micro-CT scanner (Sky-scan model 1076, Skyscan, Kontich, Belgium).
Tissue procedures
After 20 weeks post ovariectomization, animals were decapitated under Pentobarbital Sodium anesthesia (40 mg/kg). Blood samples were collected and centrifuged and the serum was stored at -20°C until biochemical analysis. The bone marrows from the femora were harvested for BMSC extraction. In addition, articular cartilage from femur and tibia was collected for chondrocyte extraction. Compact bone of long bone with removal of periosteum and endosteum was cut into small pieces and store in liquid nitrogen for RNA extraction using Trizol (Invitrogen, USA).
Electron microscopy procedure
Tissues processed for electron microscopy (EM) viewing were rinsed in 0.1 M sodium cacodylate buffer (pH 7.2), postfixed for 1 hour in 2% osmium tetroxide (OsO4) in 0.1 M sodium cacodylate buffer, dehydrated in a graded series of ethyl alcohols, impregnated with 1% uranyl acetate in 100% alcohol, and flat-embedded in Spurr resin (Electron Microscopy Sciences, Fort Washington, PA). For the flat-embedding, the sections were mounted on microslides pretreated with liquid releasing factor (Electron Microscopy Sciences). Pieces of embedded tissue were cut and glued to carrier blocks, and ultrathin sections were cut from these specimens with a Reichert ultramicrotome. The sections were mounted on mesh grids, stained with 0.4% lead citrate and 4.0% uranyl acetate using an LKB Ultrastainer, and finally viewed and images captured with a JEOL 2000EX electron microscope.
MicroCT measurement
The microarchitecture of the trabecular bone from the proximal tibiae was analyzed with microCT scanner (SkyScan 1172, Kontich, Belgium). The scanning was performed with the following parameters: 18 µm isotropic voxel size, 55 keV voltage, 109 µA current, 200 ms integration time and 4000 projections. A global threshold was chosen based on a subset of images by visual inspection and comparison of 2D and 3D binarized images with the original gray-scale images. The segmentation parameters for bone from background voxels were fixed at Sigma=1.2, Support=2, and a global threshold corresponding to 495 mmHA/cm3. The trabecular bone within tibiae was extracted from cortex with semi-automatically drawn contour at each two-dimensional (2-D) section. A total of 100 slices corresponding to a 1.5 mm region from the center of the tibiae was evaluated with the manufacturer-provided micro-CT software to determine trabecular microstructure parameters. The trabecular region was isolated from the cortical region in the 2D images by manual contouring analysis. Direct 3D measurement methods were used to calculate trabecular bone volume fraction, trabecular thickness, trabecular number and trabecular separation.
Gene expression profiling procedures
Compact bone of long bone (femora from 3 animals of each group) with removal of periosteum and endosteum was cut into small pieces and store in liquid nitrogen for RNA extraction using Trizol (Invitrogen, USA). Total RNA sample was quantified by the NanoDrop ND-1000 and RNA integrity was assessed by standard denaturing agarose gel electrophoresis. For microarray analysis, Agilent Array platform was employed. The sample preparation and microarray hybridization were performed based on the manufacturer’s standard protocols. Briefly, total RNA from each sample was amplified and transcribed into fluorescent cRNA with using the manufacturer’s Agilent’s Quick Amp Labeling protocol (version 5.7, Agilent Technologies). The labeled cRNAs were hybridized onto the Whole Genome Oligo Array (4x44K, Agilent Technologies). After having washed the slides, the arrays were scanned by the Agilent Scanner G2505C. Agilent Feature Extraction software (version 11.0.1.1) was used to analyze acquired array images. Quantile normalization and subsequent data processing were performed using the GeneSpring GX v11.5 software package (Agilent Technologies). After quantile normalization of the raw data, genes that 2 out of 3 samples have flags in Detected (All Targets Value) were chosen for further data analysis. Differentially expressed genes were identified through Fold Change filtering. Pathway analysis and GO Analysis were applied to determine the roles of these differentially expressed genes played in these biological pathways or GO terms. Finally, Hierarchical Clustering was performed to show the distinguishable gene expression pattern among samples.
Network modeling
We wish to understand how this network makes its own decision from dynamical analysis. Due to interactions among the nodes, which representing molecular pathways and modules, only in limited number of combinations the network can balance those interactions. We identified such combination by Boolean model. The interactions in the network are not real chemical reactions with known kinetic parameters since the nodes are consolidated pathways and modules. Boolean models are parameter free and can serve as a suitable starting point for modeling such biological systems. In particular, attractor analysis of Boolean models may provide insights into the long-term behaviors, i.e. observed phenotypes, of the network in response to stimulation [50]. We used the Boolean networks with threshold functions [51] to analyze this network. The attractors are listed in Table 3. Annotation for gene expression in different pathways is made mostly from literatures listed in Supplementary Materials and additional information provided by KEGG pathway [52] and www.genecards.org [53].
Boolean analysis
Boolean networks have been used for the modeling of large regulatory systems such as gene regulation. It has been test on real biological genetic networks available such as genetic network of flower development in Arabidopsis thaliana, the cell cycle networks in Saccharomyces cerevisiae and Saccharomyces pombe. Here we used Boolean networks with threshold functions to analyze this network. In the model, each node i has only two states, Si = 1 and Si = 0, representing the active state and the inactive state. The node states in the next time step are determined by the following rule:
Where aij = 1 for a green arrow from node j to node i, representing that the node j activates the node i and aij = -1 for a red arrow from node i to node i, representing that the node j inhibits the node i.We use this dynamic model to study the time evolution of the protein states. First, we study the attractors of the network dynamics by starting from each of the 28 = 256 initial states in the 8-node network. We find that all of the initial states eventually flow into one of the 19 stationary states (attractors) shown in Fig. 1(f). The calculated attractive basins are shown in Fig. 1(c), (e), (g) as well as Supplementary Figures 1 and 2.
The transition path in Fig. 1 was calculated based on biological knowledge as well as an optimized search on the network dynamics. Transition path from OVX osteoporotic osteoblast phenotype to GluSr treated state Fig. 1(g) was obtained by: starting from the OVX attractor with ESR(0), GR(1), RAR(0), PPAR(1), VDR(0), GRP120(0), Ca2+(0), NCOA(1), GluSr treatment will active Ca2+ pathway, and then the system falls into the attractive basin of an intermediate stable state (ESR(0), GR(0), RAR(1), PPAR(1), VDR(0), GRP120(1), Ca2+(1), NCOA(1)) and flows to it, by glucocorticoid receptor signaling pathway (GR) inhibition; finally by activating VDR pathway, the system reaches the GluSr treated one ESR(0), GR(0), RAR(1), PPAR(1), VDR(1), GRP120(1), Ca2+(1), NCOA(1). Osteoporosis genesis process in Fig. 1(h) was recognized as 3 steps: 1. From the normal state with ESR(1), GR(0), RAR(0), PPAR(1), VDR(1), GRP120(1), Ca2+(1), NCOA(1) to an estrogen loss situation, in which ESR is continuously lowly expressed, denoted by a grey block with slash inside; 2. By activating GR pathway, the system falls into the attractive basin of an intermediate state (ESR(0), GR(0), RAR(1), PPAR(0), VDR(0), GRP120(0), Ca2+(1), NCOA(1)) and flows to it; 3. Reactivating GR, the system enters the basin of the OVX state and goes to it.
Acknowledgments
This work was supported in part by the National 973 Project No.2010CB529200; by the Natural Science Foundation of China No.91029738 and No.81472126; and by the grants from the State Key Laboratory of Oncogenes and Related Genes (No.90-10-11) and Shenzhen Key Laboratory for Anti-aging and Regenerative Medicine (ZDSYS20140509100746679). S.M. and S.G. are supported by the NSFC grants 11135001, 81472126, and the Open Project Program of State Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, China (No. Y4KF151CJ1).
Authors’ contribution
X.Z., R.Y., and S.M. constructed molecular network, performed data and model analyses; G.Z. designed experiments; J.L., Y.L., T.L., Y.Z., B.Z., S.T., and H.G., conducted the experiments; S.G. commented and edited the manuscript; X.Z., R.Y., S.M., P.A., and G.Z. conceived the study and wrote the manuscript.
Conflicts of Interests
The authors declare no competing financial interests.
References
1. Rachner TD, Khosla S and Hofbauer LC. Osteoporosis: now and the future. Lancet. 2011; 377:1276-1287.
2. Wolf RL, Stone KL and Cauley JA. Update on the epidemiology of osteoporosis. Curr Rheumatol Rep. 2000; 2:74-86.
3. Lucas TS and Einhorn TA. Osteoporosis: The Role of the Orthopaedist. J Am Acad Orthop Surg. 1993; 1:48-56.
4. Lane JM and Nydick M. Osteoporosis: current modes of prevention and treatment. J Am Acad Orthop Surg. 1999; 7:19-31.
5. Inzerillo AM and Epstein S. Osteoporosis and diabetes mellitus. Rev Endocr Metab Disord. 2004; 5:261-268.
6. Manolagas SC. Birth and death of bone cells: basic regulatory mechanisms and implications for the pathogenesis and treatment of osteoporosis. Endocr Rev. 2000; 21:115-137.
7. Ducy P, Schinke T and Karsenty G. The osteoblast: a sophisticated fibroblast under central surveillance. Science. 2000; 289:1501-1504.
8. Zaidi M. Skeletal remodeling in health and disease. Nat Med. 2007; 13:791-801.
9. Canalis E, Giustina A and Bilezikian JP. Mechanisms of anabolic therapies for osteoporosis. N Engl J Med. 2007; 357:905-916.
10. Long F. Building strong bones: molecular regulation of the osteoblast lineage. Nat Rev Mol Cell Biol. 2012; 13:27-38.
11. Manolagas SC, O’Brien CA and Almeida M. The role of estrogen and androgen receptors in bone health and disease. Nat Rev Endocrinol. 2013; 9:699-712.
12. Lelovas PP, Xanthos TT, Thoma SE, Lyritis GP and Dontas IA. The laboratory rat as an animal model for osteoporosis research. Comp Med. 2008; 58:424.
13. Wronski T, Cintron M and Dann L. Temporal relationship between bone loss and increased bone turnover in ovariectomized rats. Calcif Tissue Int. 1988; 43:179-183.
14. Weinstein RS, Jilka RL, Parfitt AM and Manolagas SC. Inhibition of osteoblastogenesis and promotion of apoptosis of osteoblasts and osteocytes by glucocorticoids. J Clin Invest. 1998; 102:274.
15. Wu B, Xu B, Huang T and Wang J-R. A model of osteoporosis induced by retinoic acid in male Wistar rats. Acta Pharm Sin. 1995; 31:241-245.
16. Meunier PJ, Roux C, Ortolani S, Diaz-Curiel M, Compston J, Marquis P, Cormier C, Isaia G, Badurski J, Wark JD, Collette J and Reginster JY. Effects of long-term strontium ranelate treatment on vertebral fracture risk in postmenopausal women with osteoporosis. Osteoporos Int. 2009; 20:1663-1673.
17. Imai Y, Youn M-Y, Inoue K, Takada I, Kouzmenko A and Kato S. Nuclear receptors in bone physiology and diseases. Physiol Rev. 2013; 93:481-523.
18. Manolagas SC. Steroids and osteoporosis: the quest for mechanisms. J Clin Invest. 2013; 123:1919-1921.
19. O’Sullivan RP, Greenberger JS, Goff J, Cao S, Kingston KA, Zhou S, Dixon T, Houghton FD, Epperly MW, Wang H and Glowacki J. Dysregulated in vitro hematopoiesis, radiosensitivity, proliferation, and osteoblastogenesis with marrow from SAMP6 mice. Exp Hematol. 2012; 40:499-509.
20. Pfeilschifter J, Koditz R, Pfohl M and Schatz H. Changes in proinflammatory cytokine activity after menopause. Endocr Rev. 2002; 23:90-119.
21. Ralston SH, Russell RG and Gowen M. Estrogen inhibits release of tumor necrosis factor from peripheral blood mononuclear cells in postmenopausal women. J Bone Miner Res. 1990; 5:983-988.
22. Pacifici R, Brown C, Puscheck E, Friedrich E, Slatopolsky E, Maggio D, McCracken R and Avioli LV. Effect of surgical menopause and estrogen replacement on cytokine release from human blood mononuclear cells. Proc Natl Acad Sci U S A. 1991; 88:5134-5138.
23. Shanker G, Sorci-Thomas M and Adams MR. Estrogen modulates the expression of tumor necrosis factor alpha mRNA in phorbol ester-stimulated human monocytic THP-1 cells. Lymphokine Cytokine Res. 1994; 13:377-382.
24. Delany AM, Gabbitas BY and Canalis E. Cortisol downregulates osteoblast α1(I) procollagen mRNA by transcriptional and posttranscriptional mechanisms. J Cell Biochem. 1995; 57:488-494.
25. Gohel A, McCarthy M-B and Gronowicz G. Estrogen Prevents Glucocorticoid-Induced Apoptosis in Osteoblasts in Vivo and in Vitro. Endocrinology. 1999; 140:5339-5347.
26. Eijken M, Koedam M, van Driel M, Buurman C, Pols H and Van Leeuwen J. The essential role of glucocorticoids for proper human osteoblast differentiation and matrix mineralization. Mol Cell Endocrinol. 2006; 248:87-93.
27. Ao P, Galas D, Hood L and Zhu X. Cancer as robust intrinsic state of endogenous molecular-cellular network shaped by evolution. Med Hypotheses. 2008; 70:678-684.
28. Chen H, Shoumura S and Emura S. Ultrastructural changes in bones of the senescence-accelerated mouse (SAMP6): a murine model for senile osteoporosis. Histol Histopathol. 2004; 19:677-692.
29. Gong H, Jarzynka MJ, Cole TJ, Lee JH, Wada T, Zhang B, Gao J, Song W-C, DeFranco DB, Cheng S-Y and others. Glucocorticoids antagonize estrogens by glucocorticoid receptor—mediated activation of estrogen sulfotransferase. Cancer Res. 2008; 68:7386-7393.
30. Zhang Y, Leung DY, Nordeen SK and Goleva E. Estrogen inhibits glucocorticoid action via protein phosphatase 5 (PP5)-mediated glucocorticoid receptor dephosphorylation. J Biol Chem. 2009; 284:24542-24552.
31. Clapham DE. Calcium signaling. Cell. 2007; 131:1047-1058.
32. DeHaven WI, Smyth JT, Boyles RR and Putney JW, Jr. Calcium inhibition and calcium potentiation of Orai1, Orai2, and Orai3 calcium release-activated calcium channels. J Biol Chem. 2007; 282:17548-17556.
33. Mercer JC, Dehaven WI, Smyth JT, Wedel B, Boyles RR, Bird GS and Putney JW, Jr. Large store-operated calcium selective currents due to co-expression of Orai1 or Orai2 with the intracellular calcium sensor, Stim1. J Biol Chem. 2006; 281:24979-24990.
34. Kawai M, Modder UI, Khosla S and Rosen CJ. Emerging therapeutic opportunities for skeletal restoration. Nat Rev Drug Discov. 2011; 10:141-156.
35. Peng S, Liu XS, Huang S, Pan H, Zhen W, Zhou G, Luk KD, Guo XE and Lu WW. Intervention timing of strontium treatment on estrogen depletion-induced osteoporosis in rats: bone microstructure and mechanics. J Orthop Res. 2014; 32:477-484.
36. Anderson HC and Reynolds JJ. Pyrophosphate stimulation of calcium uptake into cultured embryonic bones Fine structure of matrix vesicles and their role in calcification. Dev Biol. 1973; 34:211-227.
37. Genge BR, Wu L and Wuthier R. Identification of phospholipid-dependent calcium-binding proteins as constituents of matrix vesicles. J Biol Chem. 1989; 264:10917-10921.
38. Zhou X, Cui Y, Luan J, Zhou X, Zhang G, Zhang X and Han J. Label-free quantification proteomics reveals novel calcium binding proteins in matrix vesicles isolated from mineralizing Saos-2 cells. Biosci Trends. 2013; 7:144-151.
39. Kapustin AN, Davies JD, Reynolds JL, McNair R, Jones GT, Sidibe A, Schurgers LJ, Skepper JN, Proudfoot D, Mayr M and others. Calcium regulates key components of vascular smooth muscle cell—derived matrix vesicles to enhance mineralization. Circ Res. 2011; 109:e1-e12.
40. Oh DY, Talukdar S, Bae EJ, Imamura T, Morinaga H, Fan W, Li P, Lu WJ, Watkins SM and Olefsky JM. GPR120 is an omega-3 fatty acid receptor mediating potent anti-inflammatory and insulin-sensitizing effects. Cell. 2010; 142:687-698.
41. O’Brien CA, Jia D, Plotkin LI, Bellido T, Powers CC, Stewart SA, Manolagas SC and Weinstein RS. Glucocorticoids act directly on osteoblasts and osteocytes to induce their apoptosis and reduce bone formation and strength. Endocrinology. 2004; 145:1835-1841.
42. Harada S-i and Rodan GA. Control of osteoblast function and regulation of bone mass. Nature. 2003; 423:349-355.
43. Crane JL and Cao X. Bone marrow mesenchymal stem cells and TGF-β signaling in bone remodeling. J Clin Invest. 2014; 124:466-472.
44. Marciano DP, Kuruvilla DS, Boregowda SV, Asteian A, Hughes TS, Garcia-Ordonez R, Corzo CA, Khan TM, Novick SJ, Park H, Kojetin DJ, Phinney DG, Bruning JB, Kamenecka TM and Griffin PR. Pharmacological repression of PPARγ promotes osteogenesis. Nat Commun. 2015; 6:7443.
45. Frankel SR, Eardley A, Lauwers G, Weiss M and Warrell RP. The retinoic acid syndrome in acute promyelocytic leukemia. Ann Intern Med. 1992; 117:292-296.
46. Romeo F, Costanzo F and Agostini M. Embryonic stem cells and inducible pluripotent stem cells: two faces of the same coin? Aging (Albany NY). 2012; 4:878-886.
47. Antebi B, Pelled G and Gazit D. Stem cell therapy for osteoporosis. Curr Osteoporos Rep. 2014; 12:41-47.
48. Kaech SM and Ahmed R. Memory CD8+ T cell differentiation: initial antigen encounter triggers a developmental program in naive cells. Nat Immunol. 2001; 2:415-422.
49. Scherfer C, Han VC, Wang Y, Anderson AE and Galko MJ. Autophagy drives epidermal deterioration in a Drosophila model of tissue aging. Aging (Albany NY). 2013; 5:276-287.
50. Zhu X, Yuan R, Hood L and Ao P. Endogenous Molecular-Cellular Hierarchical Modeling of Prostate Carcinogenesis Uncovers Robust Structure. Prog Biophys Mol Biol. 2015; 117:30-42.
51. Bornholdt S. Boolean network models of cellular regulation: prospects and limitations. J R Soc Interface. 2008; 5:S85-S94.
52. Kanehisa M and Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000; 28:27-30.
53. Rebhan M, Chalifa-Caspi V, Prilusky J and Lancet D. GeneCards: a novel functional genomics compendium with automated data mining and query reformulation support. Bioinformatics. 1998; 14:656-664.
