The increasing incidence of osteoporosis worldwide requires anabolic treatments that are safe, effective, and, critically, inexpensive given the prevailing overburdened health care systems. While vigorous skeletal loading is anabolic and holds promise, deficits in mechanotransduction accrued with age markedly diminish the efficacy of readily complied, exercise-based strategies to combat osteoporosis in the elderly. Our approach to explore and counteract these age-related deficits was guided by cellular signaling patterns across hierarchical scales and by the insight that cell responses initiated during transient, rare events hold potential to exert high-fidelity control over temporally and spatially distant tissue adaptation. Here, we present an agent-based model of real-time Ca2+/NFAT signaling amongst bone cells that fully described periosteal bone formation induced by a wide variety of loading stimuli in young and aged animals. The model predicted age-related pathway alterations underlying the diminished bone formation at senescence, and hence identified critical deficits that were promising targets for therapy. Based upon model predictions, we implemented an in vivo intervention and show for the first time that supplementing mechanical stimuli with low-dose Cyclosporin A can completely rescue loading induced bone formation in the senescent skeleton. These pre-clinical data provide the rationale to consider this approved pharmaceutical alongside mild physical exercise as an inexpensive, yet potent therapy to augment bone mass in the elderly. Our analyses suggested that real-time cellular signaling strongly influences downstream bone adaptation to mechanical stimuli, and quantification of these otherwise inaccessible, transient events in silico yielded a novel intervention with clinical potential.
Mechanical stimuli are anabolic for bone and hold promise to counteract skeletal fragility associated with bone loss pathologies [1], [2]. Exercise based strategies are especially attractive given the critical need for inexpensive options to treat osteoporosis worldwide [3]. However, readily complied and tolerated exercise has proved ineffective in enhancing bone mass in the elderly population most in need of such interventions [4], [5]. Studies that have examined the apparent ineffectiveness of mild exercise in the elderly using cell culture systems and animal models have sometimes led to conflicting outcomes. For instance, specific aspects of cell signaling pathways may/may not be altered with age [6], [7]. Furthermore, reports suggest that exercise training can elicit enhanced bone tissue responses in animals at advanced age [8], [9]. However, in studies where in vivo deformations and strains induced by mechanical stimuli are equivalently calibrated, observations suggest that aging markedly blunts the osteogenic response to mechanical stimuli [10], [11], and renders adaptation into a low-level binary off-on state [12]. To begin to explore underlying potential for age-related deficits in mechanotransduction function, we first focused upon observations that a single bout of mechanical stimuli (∼100 s) is sufficient to influence bone matrix secretion up to a week later [13]. Furthermore, brief stimuli (∼15 s) repeated every 24 hrs robustly enhances bone formation and bone mass [14]. These observations suggest that cell signaling activated during brief stimuli are focal events that guide unique downstream bone adaptation. A variety of second messengers (e. g. , Ca2+, NO, PGE2, cAMP, ATP) are acutely activated by mechanical stimuli. Of these, the Ca2+ ion second messenger system may be unique in its specificity given that all important aspects of mechanical stimuli that influence bone formation in vivo (e. g. , magnitude, strain rate, frequency, rest intervals) have been observed to provoke highly specific real-time Ca2+ oscillations in bone cells in vitro [7], [15]–[17]. Furthermore, while blockade of Ca2+ signaling disrupts mechano-responsive gene expression [18], [19], gap junctional communication integrates Ca2+ signaling within the bone cell network [20] and may be required to influence ‘effector’ osteoblast cell differentiation [21], [22]. Downstream of signaling through the Ca2+ ion system, activation of a variety of transcription factors (e. g. , NF-κB, JNK, NFAT) provides putative links between cell responses on the order of seconds to cell function over successive days of the week. Of these, the nuclear factor of activated T-cell family of transcription factors (NFAT c1–c4) may be unique. Specifically, NFAT activation dynamics (within minutes) is remarkably specific to Ca2+ amplitudes and frequencies known to influence distinct downstream cell functions, including proliferation, differentiation and apoptosis [23], [24]. Taken together, and given recent evidence of NFAT' s critical involvement in bone mechanotransduction [25], [26], characterization of Ca2+/NFAT signaling and age-related alterations in bone cells in vivo may prove useful. However, it is not currently possible to experimentally quantify real-time activation of Ca2+/NFAT signaling in situ within bone. Given this inaccessibility, we previously developed a model for Ca2+ signaling induced in simple networks of osteocytic cells subject to mechanical stimuli using a numerical/computational technique called agent-based modeling (ABM) [27]. ABMs have been used to explore bottom-up emergent phenomena in biological systems ranging from cell-cell interactions [28] to ecosystem dynamics [29]. In our prior study, we used the characteristic of the ABM technique to examine how local properties (e. g. , cell response thresholds) and dynamic interactions between local properties over time (e. g. , cell-cell communication) influenced the emergence of global properties (e. g. , collective Ca2+ signaling in cell ensembles). Based in part upon our original framework [27], here we present an ABM of Ca2+/NFAT signaling in situ within bone' s cellular syncytium and its relation to osteoblastic relative mineral apposition rates induced by mechanical loading (r. MAR; where r. MAR = MAR in loaded−MAR in contralateral bones; Fig 1). We hypothesized that the Ca2+/NFAT ABM would accurately simulate bone adaptation, predict aspects of the modeled pathway that were altered by aging and identify aspects amenable to intervention. To test our hypothesis, we examined whether the model could be calibrated to accurately simulate relative periosteal bone formation rates (rp. BFR) induced in vivo by a variety of mechanical stimuli in young and senescent female C57BL/6 mice [12], [30]. We subsequently used the model to predict age-related deficits within the pathway. Similar to the literature, the model predicted age-related deficits in the ability of Ca2+ oscillations to dephosphorylate NFAT [31] and in NFAT DNA binding capacity [32], as factors underlying the muted adaptation observed in senescent animals [12]. We then surveyed the literature to examine pharmaceutical agents that could modulate these deficits and identified Cyclosporin A (CsA) as one possibility [33]–[35]. Lastly, to validate model predictions, we examined whether the use of low-dose CsA as a supplement could restore the bone formation response to loading when implemented in vivo in senescent animals.
The ABM of the Ca2+/NFAT pathway incorporated sufficient parametric complexity to simulate rp. BFR induced by a variety of loading protocols in young adult [30] and senescent animals [12]. We identified an optimal ABM parameter vector by maximizing likelihoods via an optimization procedure, called simulated annealing that explored the model parameter space (Table 1). At the young and aged maximum likelihood estimate (MLE, Table 2), simulation of Ca2+ oscillations in individual cells (Fig 2,3) was qualitatively similar to that observed in in vitro systems [15], [17]. Further, NFAT signaling downstream of Ca2+ oscillations conferred a high-fidelity memory, and the resulting relative mineral apposition rates in individual cells retained information regarding the distinct osteogenic potentials of the loading protocols (Fig 2,3). Importantly, the consequent tissue level rp. BFR simulated for animal specific strains was not significantly different from in vivo data in young adult and in senescent animals (p = 0. 67,0. 57 respectively; Fig 2,3). To test for redundant model parameters, we next examined whether the ABM would remain ‘true’ upon further simplification via parameter ‘knock-outs’ (i. e. , constraining parameter values = 0, one at a time). ‘Knock-out’ of parameters resulted in simulations that were significantly different from both the experimental rp. BFR data (p≤0. 02; Table S1) [12], [30] and from rp. BFR simulated by the unconstrained ABM (p≤0. 01; Fig 4). Of note, removing parameters ER0 or r. MARx always resulted in rp. BFRi = 0 (regardless of the loading protocol ‘i’ and whether the remaining parameters were optimized). Post-hoc analysis of ABM simulations suggested that ‘knock-out’ of the remaining parameters de-couples Ca2+ oscillation from induced strains (Tεx = 0), inhibits secondary Ca2+ transients (R0 = 0), decreases Ca2+ induced dephosphorylation of NFAT (α = 0) and renders mineral apposition rates into a binary ‘on/off’ state in osteoblasts (NFATnx = 0; Fig 4). As such, attempts to simplify the model via functional ‘knock-outs’ significantly degraded model ability to simulate the in vivo data (Fig 4). The ability of the unconstrained ABM to accurately simulate the in vivo data from both young [30] and senescent animals [12] uniquely permitted an examination of aging related alterations in the model' s parameters. The null-hypothesis that aging does not alter the Ca2+/NFAT pathway ABM was rejected (6-df, p<0. 0001). Given this result and in the context of literature reports, analysis of null hypotheses that aging does not alter the parameters Tεx, ER0, r. MARx were accepted both when considered individually (p = 0. 31,0. 58 and 0. 75, respectively; 1-df) and synchronously (p = 0. 85; 3-df). Subsequently, when it was assumed that aging does not alter Tεx, ER0, and r. MARx, the null hypothesis that aging does not degrade components (1-df) was rejected for the parameter α (p = 0. 004) and NFATnx (p = 0. 002) but accepted for the parameter R0 (p = 0. 27). Lastly, we examined the influence of simulated interventions that restored (the identified) age-related deficits in components of the Ca2+/NFAT pathway. We found that synchronously restoring parameters with significant age-related deficits (i. e. , α, NFATnx) to their young optima values significantly increased rp. BFR induced by a variety of loading stimuli (p<0. 001; mean: +95%, range: +49% to +120%). We next performed in vivo experiments to validate model insights and to thereby evaluate the promise of this strategy (i. e. , interventions that could restore age-related deficits in α, and NFATnx). Specifically, senescent female C57BL/6 mice (22 Mo) were subjected to mechanical loading 3 d/wk for 3-wks with/without low-dose CsA supplements (0. 3 or 3. 0 mg/kg s. c.). Additionally, to examine the extent to which bone adaptation can be rescued at senescence, young female C57BL/6 mice underwent an identical loading protocol (without CsA supplementation). We found that in contralateral bones (not subject to exogenous loading), periosteal mineralizing surface (p. MS/BS, p = 0. 29), mineral apposition rate (p. MAR; p = 0. 82) and bone formation rate (p. BFR/BS; p = 0. 56) were not significantly different between aged animals with or without CsA supplements or compared to young animals (Fig 5). In these experiments, loading induced periosteal strains were not significantly different between senescent and young mice (p = 0. 27). While p. MS/BS was not significantly different in the loaded limbs across groups, we found that loading induced p. MAR (p<0. 01) and p. BFR/BS (p = 0. 04) were significantly lower in vehicle treated senescent mice compared with young mice (Fig 5). In contrast, loading supplemented with CsA at both 0. 3 and 3. 0 mg/kg significantly enhanced both p. MAR (p<0. 01 for both dosages) and p. BFR/BS (p<0. 01 for both dosages) compared to that in vehicle treated senescent mice and to levels not different from that in young animals (Fig 5, p>0. 55). Lastly, we examined whether the model could be used to predict rp. BFR induced in senescent animals supplemented with low-dose CsA. Specifically, to simulate the hypothesized mode of interaction between loading and low-dose CsA (Fig 6 a), model parameters significantly degraded by age (i. e. , α, NFATnx) were synchronously restored to their young optima values (Table 2). Simulated restoration of age-related deficits in α and NFATnx (that are downstream of Ca2+ signaling) did not differentially influence Ca2+ oscillations in cells around the surface in response to the loading protocol (1700 με, 50 c/d; not shown here for brevity). Restoring α and NFATnx resulted in increases in NFAT dephosphorylation in surface precursors (Fig 6 b, c) and consequent r. MAR (Fig 6 d) in differentiated osteoblasts around the bone surface compared with the case with aged model parameters (Table 2). Validating model predictions, we found that simulation of rp. BFR was not significantly different from in vivo data (Fig 5 f) obtained from senescent mice subject to loading without or with CsA, and from young mice subject to loading (p = 0. 87; Fig 6 f).
We present the development of an ABM that describes activation of Ca2+/NFAT signaling within and between cells in bone. The model incorporated sufficient complexity (i. e. , parameters and analytics) to accurately simulate relative periosteal bone formation rates induced by a variety of mechanical stimuli in young adult and senescent animals. Subsequent in silico ‘knock-out’ of ABM parameters indicated that attempts to further simplify the model significantly degraded the ability to simulate the in vivo data. Model predictions of age-related alterations (or lack thereof) in 5 of 6 ABM parameters were similar in direction to that reported individually in the literature. Model simulations suggested that restoring deficits in NFAT dephosphorylation/translocation and DNA binding capacity could substantially enhance bone' s responsiveness to loading at senescence. Finally, follow-up in vivo experiments validated the model predictions and confirmed that low-dose CsA, when used as an adjuvant to skeletal loading, can completely rescue loading induced bone formation in senescent mice. The broad assumption underlying our study is that activation of intracellular Ca2+ oscillations and subsequent activation of the NFAT pathway is a critical mechanism underlying bone mechanotransduction. As described in the introduction, the literature clearly supports this hypothesis [23]–[26]. Our model represents a first attempt at quantifying activation of the Ca2+/NFAT pathway in situ and its downstream consequences for bone tissue. However, given that a number of additional second messenger systems are also activated by brief mechanical stimuli (e. g. , cAMP, NO), and that transcription factors in addition to NFAT also modulate and influence cell and tissue adaptation downstream [1], [36], [37], we expect to explore these interacting pathways in future model refinements. Another limitation was associated with the mathematical formulations used to simulate loading induced activation of the Ca2+/NFAT pathway. The formulations were specified to efficiently (computationally) describe experimental observations and were similar to previous mathematical descriptions [17], [21]–[23], [27], [37]–[45]. While the current model parameterization undoubtedly represents a simplified description, the ABM technique would readily enable finer mathematical specifications of more fundamental mechanisms and processes via which Ca2+/NFAT signaling emerges in individual cells. One direct limitation of the ABM formulation is that model time increments are restricted to 1 s and greater (and hence, the current model cannot simulate adaptation induced by aspects of loading such as strain rates and frequencies) [2], [46]. While formulating the biological observations in a more refined manner (e. g. , differential equation based) would address this limitation, such an approach would also be expected to substantially increase the computational effort required. Also, the parameter optima for this ABM (Table 2) are relevant only to this current model. As such, any further refinement of the model would be likely to result in new optima. In the context of these observations, the accuracy of our model predictions (age-related pathway deficits, prediction of the effects of targeted interventions) does lend initial validity for the implemented formulations. We next consider some of other simplifications and limitations of the model. For instance, while the current model accurately simulates adaptation at the periosteal surface, it does not accurately predict adaptation at the endocortical surface. This is not entirely surprising given that the cells responsible for endocortical adaptation may be unique from periosteal precursors, and their distinct milieu can induce differential responses to mechanical stimuli [47]. Further, cells present at the endocortical surface may receive stimulation secondary to loading in a very different manner (e. g. , marrow pressurization [48]) compared to periosteal cells - a difference that may instead contribute to the differential response at the bone surfaces. As such, planned model expansions to simulate adaptation at the endocortical surface will have to address some of these possibilities. Another limitation is that the current model is a 2-D description of a biological process that undoubtedly unfolds in 3-D. While there is no feature within the mathematical model per se that precludes extension to 3-D, a representation of the cell network in the tibia diaphysis via imaging is first needed to incorporate 3-D network topology and anatomy within the model, and will be a focus of future model refinements. A third limitation is that the model does not explicitly incorporate parameters that distinguish between the different loading protocols in previous young (3 days/wk for 3-wks) vs senescent animal experiments (5 days/wk for 2-wks) [12], [30] from our group that were used in model calibration. The reasons for this omission were: 1) reports, including our preliminary data, which suggest that loading bone every 24 hrs vs 48 hrs does not differentially influence bone formation [49], 2) the time between fluorochrome labeling were the same in the two studies (9 days) [12], [30], 3) the current model can accurately simulate the previous data [12], [30] without an additional parameter, and 4) bone formation rates induced in our current in vivo experiments and model predictions in aged mice (1700 με, 50 c/d, 3 days/wk, 3-wks; Fig 6) were per expectations given our previous data in senescent mice (that underwent 50 cycles/d of loading for 5 days/wk over 2-wk protocols at both 1200 με and 2400 με). However, the lack of explicit model consideration of this loading bout scheduling aspect of protocols is a limitation that will need to be addressed if broader optimization questions are to be examined (e. g. , when and how often is it best to load bone?). Despite these and other limitations, analysis indicated that the model of the Ca2+/NFAT pathway was sufficient to simulate bone formation induced in both young adult and senescent animals in response to a wide variety of loading stimuli [12], [30]. In the statistical realm, refinement of our idealized model for Ca2+/NFAT signaling, parametric inclusion of other second messengers (e. g. , NO, ATP) and interacting pathways (e. g. , MAPK, Wnt signaling) [33], [50], or alternate mathematical forms will not result in significantly more accurate simulations of the data at hand. As such, we believe that this result establishes the current model as a critical first step in our exploration of bone mechanotransduction in silico. Furthermore, given the modularity of the ABM technique, our approach readily permits future ABM expansions that address limitations inherent to the current model. We expect that such iterative expansions could lead to more comprehensive models that provide increasingly nuanced representations of mechanotransduction while conferring the ability to predict and optimize bone adaptation induced by a substantially wider variety of mechanical stimuli. On the other hand, attempts to simplify the current model of the Ca2+/NFAT pathway via parameter ‘knock-outs’ significantly degraded the accuracy of simulations. Similar to literature reports [40], knock out of the parameter specifying the maximal capacity of ER Ca2+ stores (ER0 = 0) abolished intracellular Ca2+ oscillations and eliminated all further downstream outcomes. Knock-out of the parameter specifying the maximal rate of osteoblastic mineral apposition (r. MARx = 0) represented the trivial case where rp. BFR was always zero. As the ABM was unable to adapt to these knock-outs (via optimization of the remaining parameters), the parameters ER0 and r. MARx were therefore critical components of our model of the Ca2+/NFAT pathway. While knock-outs of the remaining parameters (Tεx, R0, α and NFATnx) were not synthetically lethal for bone formation (i. e. , rp. BFR ≠ 0), they nevertheless caused significant loading protocol specific declines in the accuracy of ABM simulations of the in vivo data (Fig 4). This was despite the flexibility of the model and the ability of the remaining intact parameters to adapt (via optimization and location of alternate MLEs under the knockout constraint). Taken together, these results further suggest that while the 6-parameter baseline model was sufficient and accurate, further simplification of the model (via parameter ‘knock-outs’) significantly degraded model ability to simulate the data at hand [12], [30]. The highly specific adaptation that resulted from simulated knock-out of ‘sub-critical’ components of the pathway had an alternate utility as they further elucidated the functioning of the pathway (as modeled). Specifically, knock-out of parameter Tεx decoupled strain amplitudes from the Ca2+ dynamics and caused highly unusual Ca2+ oscillations in cells (Fig 4 a). While these unusual Ca2+ oscillations occurred directly as a consequence of compensatory adaptation of remaining ABM parameters, it is unclear whether they are representative of what might occur in vivo (or in vitro) should such a knock out prove feasible. Ultimately, knockout of parameter Tεx resulted in significant differences in simulated rp. BFR in loading protocols where increasing strains induce differential adaptive responses. Knock-out of the parameter R0 prevented re-filling of stores, eliminated secondary intracellular Ca2+ transients, and despite model compensation, significantly degraded rp. BFR simulated when bones were subject to the extended rest-inserted loading protocol (that was most impacted by the lack of secondary Ca2+ transients). Knock-out of the parameter α decoupled the ‘memory’ conferred to NFAT dephosphorylation events from upstream Ca2+ oscillation histories, and impacted the simulated rp. BFR responses to cyclic loading stimuli (the secondary Ca2+ transients induced by rest-inserted loading shielded this stimulus from the effects of knocking out parameter α). However, the effect of knocking out parameter NFATnx was more ubiquitous and significantly altered simulated rp. BFR for both cyclic and rest-inserted loading protocols. Of note, experimental exploration of some of these knock-outs is possible in vivo and could form the basis for additional, future model validation. For example, knock-out of model parameter R0 could be partially achieved via functional knock-out of proteins regulating SERCA pumps [51]. Additionally, knock-out of the parameter α could be achieved, for instance, by modulating expression of negative feedback elements of the Ca2+/NFAT pathway such as RCAN1 [52]. Each of these functional knockout experimental models could provide extremely valuable tools to further validate the ABM, and more importantly, to investigate bone mechanotransduction function in finer detail. The ability to accurately simulate data in both young adult and senescent animals (without resorting to additional qualifiers) enabled us to examine age-related alterations in our model for the Ca2+/NFAT pathway. Of note, parameter alterations were not specified a priori. Similar to literature, analysis suggested that aging does not significantly alter three components of the modeled pathway: 1) the induced strain magnitude that maximizes Ca2+ amplitudes (Tεx) [7], 2) the Ca2+ sequestration capacity of the ER stores (ER0) [53], and 3) the maximal rates of mineral apposition by osteoblasts (r. MARx) [54]. Subsequent analysis predicted that aging significantly decreased Ca2+ induced dephosphorylation/translocation of NFAT (α) [31], and decreased NFAT DNA binding capacity (NFATnx) [32]. However, contrary to literature [55], predicted alterations in Ca2+ store refilling rates (R0) did not attain statistical significance. It is possible that addressing model limitations (e. g. , addressing our lack of consideration of cell viability at senescence [56]) could resolve this one apparent inconsistency. However, it also remains to be determined whether the reported deficit in R0 occurs synchronously with correctly predicted age-related deficits (i. e. , in α, NFATnx) [31], [32] and lack of deficits in other parameters (i. e. , in Tεx, ER0, r. MARx) [7], [53], [54], and in the specific context of bone response to mechanical stimuli. Taken together, we believe that the validity of the model is augmented by its ability to correctly predict the direction of age-related alterations in 5 of 6 components of the modeled Ca2+/NFAT pathway. By describing pathway alterations and their relation to bone adaptation within a consolidated framework, the ABM provided a unique tool to query the impact of targeted interventions. Interestingly, a hypothetical simulated intervention that fully and synchronously restored parameters significantly degraded by aging (α, NFATnx) to their young ‘states’ was found to nearly double bone formation response to loading at senescence. To explore and exploit this promising prediction, we considered the use of low-dose CsA supplements to mechanical stimuli. CsA is a currently approved immunosuppressant normally regarded as an inhibitor of NFAT signaling. There is substantial conflict regarding the influence of CsA upon bone adaptation, with reported effects dependent upon a number of cofactors, including experimental system and dosage [57], [58]. A recent study clarifies these inconsistencies [35], and suggests that CsA has bi-phasic effects in the trabecular compartment of bone in young adult mice (anabolic at low dosages and inhibitory at high dosages). However, to our knowledge, no previous study has examined whether and how CsA and exogenous mechanical stimuli might interact and modulate bone adaptation. Here, we hypothesized that supplementation with low-dose CsA would mitigate predicted deficits in loading (and Ca2+) induced dephosphorylation and nuclear translocation of NFAT (α) in part via CsA suppression of negative regulators such as p38 [34] and/or GSK-3β [33] (Fig 6 a). We additionally speculated that CsA supplementation could also effectively counteract predicted deficits in NFAT-DNA binding (NFATnx) by enhancing cooperative binding interactions between loading induced translocation of NFAT and CsA enhanced activation of transcription factor families such as AP-1 [35] (Fig 6 a). We found that supplementation with CsA (at both 0. 3 and 3. 0 mg/kg) significantly enhanced periosteal mineral apposition rates (by over 50%) and bone formation rates (by over 80%) compared to that induced by loading alone, and to levels not different from that observed in young animals. While the mechanisms underlying this interaction are not known at this time, we expect that a combination of RT-PCR, immunohistochemistry (for NFAT proteins and targets of the pathway such as RCAN1), and primary culture in vitro experiments might elucidate how CsA and loading interact in the aged skeleton. In this context, while the agreement between model simulations and experimental data in senescent mice supplemented with CsA validates model predictions (Fig 6), it also provides the rationale for initially investigating CsA-loading interactions as hypothesized (Fig 6a). Importantly, the low-dose CsA supplements tested here were observed to completely rescue loading induced bone formation in senescent animals, a deficit that, to our knowledge, has not previously been overcome. These pre-clinical data suggest that CsA in combination with readily complied, mild physical exercise could represent an extremely low-cost anabolic option to augment bone mass. Specifically, we estimate that low-dose CsA could be utilized as an adjunct with mild physical activity for less than $1 US/month [59]. If optimized, and efficacy is borne out in clinical trials, this intervention would enable an inexpensive, anabolic treatment option that could be readily implemented in elderly populations at risk in developed and developing countries alike. In summary, we present an agent-based model that describes, for the first time, how bone tissue formation in response to brief loading might be controlled by real-time signaling interactions within the bone cell syncytium. Despite the relative simplicity of our approach compared to complexity of bone mechanotransduction, the model was sufficiently sophisticated to enable successful calibration and hence, was able to quantitatively simulate bone formation rates induced by a variety of mechanical stimuli in both young adult and senescent animals. Post calibration, the model was found to correctly predict the directions of age-related alterations in 5 of 6 modeled components of the Ca2+/NFAT pathway and identified deficits that were amenable to therapeutic manipulation. Validating model predictions, supplementing mild mechanical stimuli with low-dose CsA was observed to completely restore the bone formation response to loading at senescence. Our model-enabled discovery provides the direct rationale to consider this currently approved pharmaceutical alongside mild physical exercise as an inexpensive, yet potent therapy to augment bone mass in the elderly. These results also further emphasize that intracellular events initiated during transient/infrequent stimuli (and deficits within) may have profound and lasting influences upon cells, tissues, and organs. Given the inaccessibility of such real-time signaling pathways in biological systems and the success of an agent-based modeling approach in this context, we speculate that an analogous approach may prove successful in understanding and predicting the behavior of other primarily homeodynamic systems facing sudden perturbations (e. g. , adaptation following musculoskeletal trauma).
The ABM was developed to examine signaling interactions between cells at the murine tibia mid-shaft, the site where experimental data had previously been measured (Table S1) [12], [30]. To determine cell network topology at the tibia mid-shaft (Fig 1), we obtained thin (5 µm) decalcified sections from the tibia mid-shaft of young adult (4 Mo, n = 6) and senescent female C57BL/6 mice (22 Mo, n = 5). Sections were Toluidine blue stained per standard protocols [60], imaged under white light excitation (200 X), combined into composites and oriented anatomically. For each composite image, the cortical bone was sub-divided and osteocytic cell lacunae present within each sub-sector (∼100 µm2,384 sub-sectors) were counted for each bone and averaged across bones (young, senescent separately). This provided an estimate of the ‘average’ count of osteocytic cells present within each of 384 anatomically oriented sub-sectors. To define bone anatomy, we also imaged un-decalcified sections separately obtained from young adult (90 µm, 4 Mo, n = 5) and senescent female C57BL/6 mice (22 Mo, n = 5). The cross-section chosen for ABM implementation was representative of the average anatomy (i. e. , in cortical area, thickness, moments of inertia). Upon this cross-section, each of 384 sub-sectors was ‘seeded’ with the average number of cells previously determined within the sub-sector. Cell locations within each sub-sector were then specified via a non-overlapping, random walk process. With regard to the effector cells (collectively termed ‘precursors’ in the model), per literature [61], we assumed that a single layer of precursor cells were present around the cortex and were positioned at the endocortical and periosteal surfaces (along 96 equal angle sectors). While immediately adjacent precursors were assumed to be functionally coupled, osteocyte - osteocyte coupling and osteocyte-precursor couplings were specified, per literature [62], if cell-cell separation was within a canalicular length of 50 µm (Fig 1; where ‘functional coupling’ simulates an ability to pass Ca2+ ions between coupled cells). Similar to literature [63], in the idealized network for young animals, the mean (± S. D.) number of functional connections per periosteal precursor was 4. 0±1. 6 and that for osteocytes was 7. 9±2. 3, and the average canalicular length was 34±12 µm. In the network for senescent animals, the mean (± S. D.) number of functional connections was 4. 2±1. 5 per periosteal precursor and 10. 7±2. 7 for the osteocyte (p<0. 001 vs young), and the average canalicular length was 34±12 µm. Of note, osteocytic lacunar density was significantly increased in these aged vs young female C57BL/6 mice (940±80 vs 626±90/mm2; p = 0. 001). However, this data is contrary to declines in lacunar density typically observed with age in healthy human bone [64], but similar to increases observed in osteoporotic vs control bones [65]. Given that, to our knowledge, there is no comparable report on age-related alterations in lacunar density in mice, the increase in lacunar density observed here may be specific to the female C57BL/6 mice at this advanced age (22 Mo) and remains to be further investigated. Given the cell network topology (Fig 1), we next determined localized tissue strains at the location of each cell induced by loading of the tibia. Using a combined strain gauging and finite element modeling approach, we have previously calibrated the force, moment and tissue level strain environment that is induced at the tibia mid-shaft by loading [12], [30]. As such, for each independent loading protocol (Table S1) [12], [30], we utilized beam theory along with known force and moment boundary conditions at the tibia mid-shaft to determine tissue level strains at the location of each mechanosensory osteocytic cell over the time duration of the loading protocol. The ABM is a biophysical model that relies on experimentally identified ‘building blocks’ to algorithmically describe how activation of the Ca2+/NFAT pathway by loading of bone influences mineral apposition rates induced in individual osteoblasts (Fig 1, Table S2, Protocol S1). In the model, intracellular Ca2+ oscillations induced in individual cells (osteocytes and precursors) at time ‘t’ during application of loading, [Ca2+]i (t), emerged via an interplay between Ca2+ influx due to loading induced strains at the location of the osteocytic cell (εi (t) ), cell-cell gap junctional exchange of Ca2+ ions and efflux from ER Ca2+ stores as follows: (1a) where, (1b) (1c) and, (1d) where, Tεx is the parameterized threshold strain magnitude required to maximize the amplitude of Ca2+ oscillations in mechanosensory osteocytes (normalized to 100% values), ni is the number of cells (osteocytes and/or precursor cells) functionally coupled to the ith cell, j is an element of the cell set ni, TCn represents threshold ‘incoming’ Ca2+ ion levels required to initiate oscillations in the recipient cells (specified via preliminary simulations to be 2. 5%), and ER0 and R0 are the parameterized maximal Ca2+ store capacity and maximal store recovery rates, respectively. The specific mathematical formulations for osteocyte sigmoidal dose-response relation to strain (eqn 1b) were based upon observations of Ca2+ oscillations induced in bone cells by mechanical stimuli [17], [38]. Further, cell-cell communication (eqn 1c) is a simple description (via averaging) of Ca2+ ion propagation between cells based upon reports [38], [39] and is similar to our previous descriptions [27]. The formulation for ER Ca2+ store depletion (eqn 1d) is similar to previously [27]. Lastly, the formulation for ER Ca2+ store recovery upon cellular quiescence (eqn 1d), models the action of the SERCA pump working against a Ca2+ ion concentration gradient and is similar to literature reports [40]–[42]. Given observations that Ca2+ oscillations in osteocytic networks hold potential to regulate surface osteoblastic cells [21], [22], we next consider the downstream consequences of Ca2+ oscillations upon NFAT transcription factor dynamics. Specifically, downstream of increased intracellular Ca2+ oscillations, cytoplasmic NFAT dephosphorylation and nuclear translocation in surface precursor cells was modeled to be reflective of a ‘memory’ of Ca2+ oscillation characteristics such as amplitude, frequency and durations, and included an implicit consideration of negative feedback elements [23], [37] as follows: (2a) where, the parameter ‘α’ controls the amount of NFAT dephosphorylated based upon prior Ca2+ oscillations. The quantity of NFAT translocated into the ith cell' s nucleus () over the duration of an entire loading bout (tbout) was modeled simply as: (2b) Finally, given observations that the activation of the Ca2+ signaling and the NFAT pathway code for highly specific cell proliferation, differentiation and apoptotic events [23], [37], [43], we assume that nuclear NFAT induced by a given loading bout controls the differentiation of precursor cells and ultimate, the relative mineral apposition rates by differentiated osteoblasts () as follows: (3) where, NFATnx is the maximal nuclear NFAT DNA binding capacity, and r. MARx, the maximal relative mineral apposition rate by osteoblasts. The mathematical form (eqn 3) reflects considerations that transcriptional regulation is likely to be a threshold-driven process [44], [45]. Simulations with the ABM involved solving equations (1–3), given osteocyte level strains derived over each second of a loading bout in the context of the ABM' s 6 independent parameters (i. e. , specific values assigned to Tεx, ER0, R0, α, NFATnx, r. MARx prior to execution of the simulations). To initialize the ABM (at t = 0), strains induced by loading at osteocyte cell locations were 0 (i. e. , εi (t = 0) = 0), and the stores of all cells were assumed to be filled to capacity (i. e. , ERi (t = 0) = ER0). Based upon preliminary studies, we simulated intracellular Ca2+ oscillations for an extra 200 seconds beyond cessation of a given loading bout to permit the unfolding of cell-cell cross talk, and NFAT dephosphorylation for an extra 100 s beyond cessation of Ca2+ oscillations to more fully account for prior Ca2+ oscillation histories. We implemented standard dynamic histomorphometry methods to determine bone formation indices simulated by the model at the tissue level [66]. We first made the simplifying assumption that repeated loading of bone over days of the week does not adaptively alter cells (i. e. , make them more/less sensitive to repeated bouts of loading). Therefore, r. MAR induced in individual osteoblasts over a loading protocol (e. g. , loading for 9 bouts provided over 3 days/wk for 3 weeks) was simply determined as the average of induced by each bout of loading. Second, surface referent relative mineralizing surface (r. MS/BS) at the tissue level was determined simply as the proportion of surface osteoblasts with non-zero. Tissue level r. MAR was then determined as the average of non-zero at each bone surface. Finally, tissue level relative bone formation rates (surface referent) was determined as r. BFR/BS = r. MS/BS×r. MAR at each of the surfaces. The ABM was implemented in C++ on an Apple Mac Pro (2×quad core processors, 3. 0 GHz, 8 Gb RAM), and simulations were performed at 1-Hz (i. e. , smallest time unit was 1 s). The ABM simulates activation of the Ca2+/NFAT pathway in cells and ultimately its relation to relative bone formation rates at the tissue level. However, given that the ABM is an idealized biophysical model, the parameter values must be estimated from observed rp. BFR induced by a variety of mechanical stimuli [12], [30]. We modeled the observed rp. BFRi from the ith loading protocol as normally distributed with population mean μi and variance σi2. The ABM with values specified for the 6 parameters θ = (Tεx, ER0, R0, α, NFATnx, r. MARx) at the initiation of the simulations then predicts mean values μi (θ) = rp. BFRi for each loading protocol ‘i’. Under the normal error model, the log-likelihood function for the ABM is then [67]: (4) where ‘yij’ is the rp. BFR measured in the jth animal for the ith loading protocol, μ (θ) is the vector of means given by the ABM and σ2 is the vector of population variances for rp. BFR under the different loading protocols, ‘Nprotocol’ is the number of independent loading protocols, ‘nai’ is the number of animals that underwent the ith loading protocol, ‘na’ is the total number of experimental observations. The maximum likelihood estimate (MLE), for the vector of means, μ (θ) and variances σ2, is simply the set of model parameters that maximizes L (μ (θ), σ2). For a given vector of mean parameters μ (θ) of σ2, maximization with respect to the variances σ2 is achieved by setting: (5) However, maximization with respect to the ABM parameters θ is more challenging because the likelihood function contains local maxima. To address this, we applied the simulated annealing procedure (SA) [68], with random restarts to find the MLE for the ABM parameters. Briefly, the SA algorithm is based on annealing in metallurgy and is a probabilistic heuristic applied to a stochastic process (random walk) that seeks to determine a material configuration where energy is minimized [68]. In the context of our ABM, the SA procedure explores the 6-D model parameter space seeking to identify parameter values that would maximize the log-likelihood function L (μ (θ), σ2). Using this procedure, separate ABM models were calibrated to the rp. BFR data obtained from young animals subjected to 10 protocols and senescent animals subjected to 7 protocols (Table S1) [12], [30]. We denote the resulting maximum likelihood estimates of the parameter vector for each ABM by and respectively (Table 2). To analyze model simulations, we addressed the following questions (described in more detail in supporting information, Text S1): (a) Are the calibrated ABM models for young and senescent animals compatible with the observed data [12], [30]? (b) What are 95% confidence intervals for the parameters θyng and θaged? c) Is there evidence against the null hypothesis that each of the ABM parameters takes value zero (i. e. , simulating parameter ‘knock-outs’)? (d) Is there evidence in the observed data against the null hypothesis that the population values of the ABM parameters are the same for young and senescent animals? (e) If null hypothesis in (d) is rejected, is there evidence against the null hypotheses that specific model parameters or combination of parameters are unchanged by age? (f) Do the calibrated ABMs make predictions about the results of hypothetical interventions that would restore specific ABM parameters in senescent animals to their young values? (g) Do ABM simulations of loading induced rp. BFR, upon complete restoration of parameters significantly degraded by aging, predict data from senescent animals subject to loading supplemented with CsA? We examined these questions using a series of likelihood ratio tests (within MATLAB) and factorial ANOVAs (within SPSS) where appropriate; please see supporting information for additional details (Text S1).
Post-menopausal and age-related osteoporosis afflicts large segments of the population and can markedly increase skeletal fragility. Bone fractures that occur as a consequence substantially increase health care expenditures and raise levels of morbidity. While strategies that prevent further loss of bone exist, options that compensate for bone loss accrued over age are less numerous. Physical exercise holds promise in this realm. However, in part due to deficits in how cells within bone respond to skeletal loading, readily complied exercise has proved ineffective in enhancing bone mass in the elderly. In this study, we examined whether the ability of physical exercise to increase bone mass can be restored at advanced age. To this end, we developed a computational model describing how a specific aspect (or pathway) activated in bone cells by skeletal loading may be altered with age. Our model proved successful in describing age-related pathway alterations and identified specific deficits that were amenable to therapeutic manipulation. We subsequently discovered that when an extremely inexpensive, currently approved pharmaceutical is used as a supplement, bone response to skeletal loading was completely rescued in aged animals. We believe that this result provides the rationale to consider this approach as a means to increasing bone mass in the elderly.
physiology/cell signaling
computational biology/transcriptional regulation
biophysics/theory and simulation
computational biology/signaling networks
Rescuing Loading Induced Bone Formation at Senescence