The Role of Osteocytes in Targeted Bone Remodeling: A Mathematical Model

Abstract

Until recently many studies of bone remodeling at the cellular level have focused on the behavior of mature osteoblasts and osteoclasts, and their respective precursor cells, with the role of osteocytes and bone lining cells left largely unexplored. This is particularly true with respect to the mathematical modeling of bone remodeling. However, there is increasing evidence that osteocytes play important roles in the cycle of targeted bone remodeling, in serving as a significant source of RANKL to support osteoclastogenesis, and in secreting the bone formation inhibitor sclerostin. Moreover, there is also increasing interest in sclerostin, an osteocyte-secreted bone formation inhibitor, and its role in regulating local response to changes in the bone microenvironment. Here we develop a cell population model of bone remodeling that includes the role of osteocytes, sclerostin, and allows for the possibility of RANKL expression by osteocyte cell populations. We have aimed to give a simple, yet still tractable, model that remains faithful to the underlying system based on the known literature. This model extends and complements many of the existing mathematical models for bone remodeling, but can be used to explore aspects of the process of bone remodeling that were previously beyond the scope of prior modeling work. Through numerical simulations we demonstrate that our model can be used to explore theoretically many of the qualitative features of the role of osteocytes in bone biology as presented in recent literature.

About the Authors:

Jason M. Graham

* E-mail: jason.grahamscranton.edu

Affiliation: Department of Mathematics, University of Scranton, Scranton, Pennsylvania, United States of America

Bruce P. Ayati

Affiliations Department of Mathematics, University of Iowa, Iowa City, Iowa, United States of America, Program in Applied Mathematical and Computational Sciences, University of Iowa, Iowa City, Iowa, United States of America, Department of Orthopaedics and Rehabilitation, University of Iowa, Iowa City, Iowa, United States of America

Sarah A. Holstein

Affiliation: Department of Internal Medicine, University of Iowa, Iowa City, Iowa, United States of America

James A. Martin

Affiliation: Department of Orthopaedics and Rehabilitation, University of Iowa, Iowa City, Iowa, United States of America

Introduction

The skeleton is a highly dynamic organ that undergoes constant remodeling throughout life. Remodeling is involved in systemic mineral homeostasis and is required for the routine repair of microfractures and for adaptation to mechanical stress. Even after the skeleton reaches its adult size, bone remains metabolically active. Bone remodeling involves the coupling of bone resorption and bone formation, and diseases are characterized by imbalances between formation and resorption. This includes such common disorders as osteopenia and osteoporosis [1], as well as malignant processes such as multiple myeloma and other cancers that metastasize to the bone [2], [3]. As these disorders contribute to morbidity and mortality in millions of people worldwide, there is great interest in improving our understanding the processes that regulate bone remodeling.

Bone remodeling is characterized by complex mechanical and biochemical signaling pathways. Remodeling typically occurs for one of three reasons: to maintain mineral homeostasis; to adapt to mechanical changes; and to repair damage [4], [5]. Remodeling related to mineral homeostasis need not occur at a specific site, and is therefore referred to as non-targeted remodeling [4], [5]. In contrast, remodeling due to mechanical changes or to damage is site-specific and is termed targeted remodeling [4], [5]. The key cellular players in targeted bone remodeling include the osteoblasts and osteoclasts. There is also increasing evidence to support the role of osteocytes as active participants in this remodeling process [6], [7].

Osteocytes are mature osteoblasts that become embedded within the bone matrix during the formation phase of bone remodeling. During this process osteocytes develop cytoplasmic processes which run through the canaliculi, forming a communication network that can convert mechanical signals into biochemical signals [6]–[8]. Osteocytes are the most abundant and long-living bone cells and are found throughout the skeleton, in contrast to the much shorter lived osteoblasts and osteoclasts which are only transiently found on the surface of the bone. Osteocytes actively secrete growth factors that stimulate bone formation, as well as sclerostin, which inhibits bone formation [7]–[9]. The complex relationship between osteocytes, osteoblasts, osteoclasts, and their precursors, is continuing to be elucidated.

A number of mathematical and computational models of various aspects of bone physiology have been proposed. Many of these models focus on either bone fracture healing, [5], [10]–[14], or bone remodeling, [1], [15]–[31]. These models take different approaches, and focus on various aspects of the remodeling process. The early modeling approach described in [21], [22] makes use of the biochemical systems analysis formalism of Savageau (e.g., [32]–[36]) developed to represent bone remodeling as the effect of changes in bone cell populations; a so-called “cell population” model. These models are comprised of a system of two ordinary differential equations (ODEs) representing the cell population dynamics of osteoclasts and osteoblasts. In addition, they are chemically implicit in that they do not track chemical quantities as functions of time or location, but instead incorporate the biochemical mechanisms in an abstract way. Furthermore, these models do not attempt to represent biomechanical mechanisms of remodeling. Instead, the authors relate remodeling cell populations to the change in bone mass, or bone volume, thus allowing a way to track the physical effects of resorption and formation. This approach was used in [37] to model osteomyelitis. Others, e.g., [25], [26], have modeled only the resorption phase of remodeling. These models do not track remodeling cell populations, but rather cellular activity represented by enzyme reaction equations. All of these models, as well as most others, incorporate the actions of the most prominent signaling pathways involved in bone remodeling, as well as various other autocrine and paracrine signaling pathways depending on the specific considerations.

The models presented in [16], [18], [28], [29] are also cell population models that describe bone remodeling cell dynamics as controlled by autocrine and paracrine signaling. These models use a biochemically explicit approach, in that the model equations contain terms corresponding to interactions of remodeling cells with specific cytokines. Another distinction is that while the models represented in [21], [22], [25], [26] represent remodeling on the scale of an individual basic multicellular unit (BMU), the models in [17], [23], [27], [28] average over a volume of bone. By a BMU we mean a localized collective of bone cells working as an organized unit to carry out a single, complete cycle of the remodeling process, see e.g., [38]. We note that the models found in [1], [5], [6], [18], [44], [45] consider remodeling of individual BMUs, but in a spatially explicit context, which is beyond the scope of this paper. However, the model developed herein can be easily modified to obtain a spatially explicit model.

The complexity of the skeletal system and its interactions with the rest of the body limits the ability of a single model to capture all of the relevant biological, biochemical, and biophysical mechanisms of remodeling on all scales simultaneously. Moreover, parameterization of a given model is difficult, given the dependency on accuracy and availability of data. Here we introduce a novel mathematical model for bone remodeling that maintains those fundamental mechanisms captured in previous models, while incorporating biological aspects of bone remodeling that have not previously been considered. In particular, this model includes osteocytes, their biochemical processes, and their interactions with other bone remodeling cells, and thus represents a significant advance to the field.

Materials and Methods

Mathematical Modeling

In developing a new mathematical model for bone remodeling we have chosen to employ the biochemical systems analysis formalism used in e.g. [1], [15], [20]–[22], [30], but extended to incorporate further biological detail as summarized in Figure 1. We chose this approach in anticipation of the type of data to which we may eventually have access, and the types of questions which are of concern to us. For example, explicit chemical data are not obtained in current clinical practice. Our approach does represent many of the mechanisms contained in [17], [23], [28], [29], albeit implicitly. In particular, we incorporate the actions of an important signaling pathway based on the specific molecules: receptor activator of nuclear factor -B (RANK); and osteoprotegerin (OPG) [2], [3], [39]–[42]. These cytokines, plus the RANK ligand, form a pathway commonly known as the RANK/RANKL/OPG pathway. We also incorporate the actions of growth factors such as transforming growth factor (TGF-) [43], [44], and other cytokines on bone remodeling cells. It is well known that RANKL is a key cytokine in the differentiation process of osteoclast cells, while OPG, which is produced by differentiated osteoblastic cells, has been shown to function as an inhibitory factor for osteoclastogenesis [38].

[Figure omitted. See PDF.]

Figure 1. Interactions between bone cell populations.

https://doi.org/10.1371/journal.pone.0063884.g001

We have developed a cell population model for osteocyte-induced targeted bone remodeling. This model consists of the osteocyte, pre-osteoblast, osteoblast, and osteoclast cell populations; the interactions of these cells with one another, and through the power law formalism the autocrine and paracrine signaling among these cells. Figure 2 shows the equations that make up the model, these are described in detail below. In this model we explicitly include the population of pre-osteoblasts to emphasize the switch of cells in the osteoblast lineage from being osteoclastogenic, i.e., osteoclast generating, to being osteogenic, i.e., bone forming, see e.g., [45], [46]. Specifically, there is a clear distinction between the signaling behaviors of pre-osteoblasts and osteoblasts that is important for the considerations taken up in this work. On the other hand, for the purposes we have in mind there is no similar a priori reason to explicitly represent a population of pre-osteoclasts.

[Figure omitted. See PDF.]

Figure 2. System of ordinary differential equations constructed, using the biochemical systems analysis formalism [32]–[36], to model osteocyte-induced targeted bone remodeling.

https://doi.org/10.1371/journal.pone.0063884.g002

A general assumption is that there is a large pool of mesenchymal stem cells available to differentiate into pre-osteoblasts [28]. Similarly, we assume there is a large pool of osteoclast progenitor cells which are available to differentiate to fully committed mature osteoclasts [28]. Such cell differentiation is determined by autocrine and paracrine signaling discussed in more detail below, see also [8]. We assume that some percentage of pre-osteoblasts differentiate under the influence of autocrine and paracrine signaling while some percentage undergo apoptosis. We also assume that some percentage of mature osteoblasts will undergo apoptosis, and some percentage of osteoblasts will become embedded in the bone matrix as osteocytes. Figure 1 outlines the assumptions that we make to construct the mathematical model. The details of the assumptions and how they influence the development of a mathematical model is described in detail below.

In the following we denote by S(t), or simply , the osteocyte cell population at given time t. Sclerostin is produced by osteocytes and inhibits the Wnt/β-catenin pathway, [6], [8]. Wnt is known to promote osteoblastic proliferation and differentiation [7]. We incorporate the effects of sclerostin and the Wnt/β-catenin pathway into the mathematical model through a term of the form where and Ks is a parameter that describes the relation between osteocyte apoptosis and decrease in sclerostin inhibition. The idea is that, for a threshold level Ks of osteocytes, there is sufficient sclerostin production to inhibit local Wnt signaling. When osteocytes die, the sclerostin level decreases, [6], [8]. This releases osteoblast precursor cells from Wnt inhibition, thereby initiating a cycle of targeted bone remodeling. We note that the term is a dimensionless quantity that is used to modify the standard biochemical systems analysis formalism in order to accurately capture the action of sclerostin signaling as previously described. It is also important to note that in the biochemical systems analysis formalism the basic relations used to represent signaling are nonlinear. However, it may be the case that the exponents used take on numerical values equal or close to unity. This should not be misconstrued as an a priori assumption of linearity within the model.

Equation (1) in Figure 2 describes the dynamics of the osteocyte cell population. This equation simply states that osteocytes are mature osteoblasts that become embedded in extra cellular matrix at a given rate . We further note that there is no death term for osteocytes in equation (1). This is due to our assumption that, over the time scale of a single event of targeted remodeling considered here, the most significant influence on osteocyte apoptosis is the initial biomechanical action that begins remodeling [9], [47]. This appears in the mathematical model as an initial condition.

The pre-osteoblast cell population at a time t is denoted by P(t), or simply P. Pre-osteoblasts are differentiated mesenchymal stem cells. We assume that this differentiation is controlled by osteocytes through the sclerostin, Wnt/β-catenin pathways, and various growth factors. The effectiveness of sclerostin regulations of the differentiation of mesenchymal stem cells to become pre-osteolbasts is represented mathematically by Where g22 is a dimensionless parameter. Thus when osteocytes undergo apoptosis due to microdamage, local mesenchymal stem cells differentiate to pre-osteoblasts due to the resulting signaling. Moreover, the pre-osteoblasts are free to proliferate and differentiate to mature osteoblasts since they have been released from Wnt inhibition. Equation (2) in Figure 2 describes the dynamics of the pre-osteoblast cell population. This equation states that pre-osteoblasts are differentiated from a large pool of stem cells at a rate in response to signaling molecules produced by osteocytes, that pre-osteoblasts proliferate at a rate under the influence of autocrine signaling provided this is not inhibited by sclerostin. Furthermore, pre-osteoblasts differentiate to become mature osteoblasts at a rate . This is under the influence of autocrine and osteoclast regulated paracrine signaling. Finally, pre-osteoblasts undergo apoptosis at a rate .

The osteoblast cell population at a time is denoted by , or simply . Equation (3) in Figure 2 describes the dynamics of the mature osteoblast cell population. This equation states that osteoblasts are differentiated pre-osteoblasts, that osteoblasts undergo apoptosis, and that some osteoblasts are embedded in the extra cellular matrix during formation to become osteocytes. Notice that the term also appears in equation (1) of Figure 2, representing the embedding of osteoblasts that become osteocytes, and the term , also appears in equation (2) of Figure 2, which corresponds to the differentiation of pre-osteoblasts to become mature osteoblasts. Here the parameter describes the pre-osteoblast autocrine signaling. The parameter represents the effects of osteoclast derived paracrine signaling on pre-osteoblasts. This could represent, for example, the effects of TGF- on pre-osteoblasts as described for example in [28].

The osteoclast cell population at a time is denoted by , or simply . The equation (4) in Figure 2 describes the dynamics of the osteoclast cell population. This equation states that mature osteoclasts come from the differentiation of a large pool pre-osteoclasts at a rate . This differentiation is influenced essentially by the RANK/RANKL/OPG pathway. Thus the term , describes the effects of this pathway. The dimensionless parameter models the effectiveness of sclerostin regulation of osteoclastogensis. One novel feature of this model is that we have included osteocytes as a source of RANKL as discussed in recent literature [6]–[9], [48]. The parameter represents the effect of osteocyte derived RANKL signaling on osteoclastogenesis. We also retain pre-osteoblast derived RANKL signaling via the parameter . While we have explicitly included the pre-osteoblast cell population dynamics, we have neglected to explicitly include a dynamic cell population for pre-osteoclasts. This is based on a simplifying assumption that, over the relevant time span, there is a steady pool of circulating pre-osteocalsts that may be recruited for differentiation to active osteoclasts. This simplifying assumption is common the literature on mathematical modeling of bone remodeling, e.g., [15]–[18], [21]–[23].

We note that either of these, or , could take on the value . Thus advances in the understanding of the relative role of osteoblasts and osteocytes in RANKL production would result in parameter changes in our model, but not a change in the model structure itself.

The term represents the effect of OPG acting as a decoy receptor for RANKL. Typically the parameter takes on negative values, and since has as a steady state value, we add a sufficiently small number to avoid dividing by zero. This represents the factor of production when . We also have a term for the osteoclast cell death at a rate .

We denote by , or simply , the bone mass at a given time . We follow [22], [28] to develop an equation for the change of bone mass (or bone volume depending on the scaling) over time. The rate of change of bone mass has the form

As in e.g., [28], we assume that the amount of bone being resorbed is proportional to the osteoclast population, while the amount of bone being formed is proportional to the osteoblast population.

We can summarize the mathematical model with the following verbal description:

1. change in osteocytes = increase due to embedded osteoblasts

2. change in pre-osteoblasts = increase due to differentiation of stromal cells (released from sclerostin or exposure to growth factors) + proliferation of pre-osteoblasts (autocrine signaling of Wnt and growth factors) – differentiation to osteoblasts (growth factors) - apoptosis

3. change in osteoblasts = increase due to differentiation of pre-osteoblasts (growth factors) – apoptosis – embedding as osteocytes

4. change in osteoclasts = increase due to differentiation of pre-osteoclasts (due to RANKL and limited by OPG) – apoptosis

5. change in bone mass = increase due to activity of osteoblasts – activity of osteoclasts

The corresponding full mathematical model is the set of equations shown in Figure 2. The definitions of the variables and parameters appearing in the model equations of Figure 2 is summarized in table 1.

[Figure omitted. See PDF.]

Table 1. Definitions of Symbols Used in the Paper.

https://doi.org/10.1371/journal.pone.0063884.t001

Initial Conditions

The initial conditions for the mathematical model are given by , where is a constant that represents osteocyte apoptosis, which corresponds to a drop in sclerostin levels, , and . Observe that the system of model equations (1)-(4) has as a steady state, the values , , , and . The significance being that this steady state allows one to model the initiation of an event of targeted bone remodeling. In particular, our model does not assume the presence of a population of committed pre-osteoblasts or active osteoclasts in order to begin a remodeling cycle. This is in contrast with previously published mathematical models of bone remodeling which require either a population of active osteoclasts, or a population of pre-osteoblasts to initiate dynamics corresponding to physiological bone remodeling. However, it should be noted that physiological bone remodeling also involves the movement of bone lining cells away from the bone surface. This aspect of bone remodeling is complicated and not yet fully understood. The steady state value for bone volume, denoted depends significantly on the parameter values. For instance, if the parameter values are such that there is an imbalance of osteoblasts and osteoclasts, in favor of osteoblasts, then there will be over remodeling.

Numerical Simulations

The mathematical model is a system of coupled nonlinear ordinary differential equations (ODEs). To use the mathematical model to simulate bone remodeling we must solve this system along with the specified initial conditions. In order to do so we use the MATLAB software for scientific computing. In particular, we use the MATLAB ODE suite [49], and the program ode23s to integrate the problem (1)–(5) with the above initial conditions. For a further discussion of the numerical techniques underlying the codes in the MATLAB ODE suite see [50], [51].

Parameters

Table 2 lists the values of the parameters in the model equations (1)–(5) that we have taken as baseline values. These particular values are chosen so that a cycle of normal, targeted bone remodeling is completed, i.e., bone cell populations reach a steady state value, after approximately one hundred days.

[Figure omitted. See PDF.]

Table 2. Parameter Values.

https://doi.org/10.1371/journal.pone.0063884.t002

Results

Baseline Case

Figure 3 shows the simulation results for bone cell populations during a single cycle of normal targeted bone remodeling. Here, normal targeted bone remodeling is defined to be a complete remodeling cycle that takes place over a period of one hundred days, in which the amount of new bone formed is equal to the amount of old or damaged bone resorbed. One characteristic of normal bone remodeling that manifests itself quantitatively is that the steady state value of bone volume, here and below denoted by , is equal to one hundred percent. The simulations shown in Figure 3 result from solving the model equations (1)–(5) with the above initial conditions, and parameter values shown in table 2. Figure 3(a) shows the dynamics of the osteocyte cell population during an event of targeted bone remodeling. Initially there is a decrease from the steady state value in the osteocyte cell population. This corresponds to local osteocyte apoptosis, and results in a decrease in local sclerostin levels. This releases stromal cells from sclerostin inhibition, allows for Wnt signalling, and results in proliferation and differentiation of pre-osteoblast cells, Figure 3(b). There follows an increase in osteoblast cell number due to differentiation of pre-osteoblasts, and differentiation of osteoclast pre-cursors to mature osteoclasts, Figure 3(c) and Figure 3(d). Remodeling ceases once the local osteocyte cell population is replenished, and the osteocyte cell network is reestablished, returning sclerostin expression back to sufficient levels thereby inhibiting bone turnover.

[Figure omitted. See PDF.]

Figure 3. Dynamics of bone cells during a single event of targeted bone remodeling.

The dynamics of osteocyte (a), pre-osteoblast (b), osteoblast (c), and osteoclast (d) populations during an event of targeted bone remodeling.

https://doi.org/10.1371/journal.pone.0063884.g003

Figure 4 shows the dynamics of bone volume during a single event of targeted bone remodeling. An increase in the osteoclast cell numbers results in bone resorption and the decrease in bone volume, while increases in osteoblast cell number results in bone formation. In this case bone turnover is completely balanced in the sense that the amount of new bone formed equals the amount resorbed. Mathematically this corresponds to a steady state value of bone mass, .

[Figure omitted. See PDF.]

Figure 4. Dynamics of bone volume during a single event of targeted bone remodeling.

https://doi.org/10.1371/journal.pone.0063884.g004

Role of Osteocyte RANKL Production

At the scale of a BMU, bone cells are easily affected by changes in the microenvironment. Thus we expect that small changes in the model parameter values could significantly impact bone turnover during remodeling. In this section we study the influence of expression of RANKL by osteocytes on targeted bone remodeling. In particular, we are interested in the question of whether local RANKL expression by pre-osteoblasts, or by osteocytes has the greater impact. Following [28], we take as a control parameter the change in bone volume. This provides a measurable quantity by which we can assess the reasonableness of the mathematical model. If then there is over-remodeling. If then there is under-remodeling.

That RANKL/OPG expression by cells in the osteoblastic lineage is fundamental for the activation and control of bone remodeling is by now well established. Until recently, the paradigm was that pre-osteoblasts serve as the primary source for RANKL and active osteoblasts serve as the primary source of OPG [46]. More recent research highlights that osteocytes also express RANKL at significant levels, and have a strong influence on the activation and regulation of remodeling [6]–[9], [48], [52]. To explore this we vary the parameters and in equation (4) over the range from zero to two. These parameters correspond to the effectiveness of the expression of RANKL by osteocytes, and pre-osteoblasts respectively. By varying these parameter values simultaneously we gain insight into how RANKL expression by osteocytes and pre-osteoblasts impacts the regulation of a cycle of targeted bone remodeling. We note that the system is very sensitive to changes in the effectiveness of RANKL expression. In particular, small perturbations to the amount of RANKL can lead to pathological remodeling. First we compute as a simultaneous function of and , and observe that values of and that are simultaneously large, in this case , result in physiologically unrealistic over resorption (under remodeling). Next we compute as a function of each of and individually with the other fixed. We observe that it is possible to have values such that either , or and resulting in . Table 3 lists some specific values of , and that result in normal targeted bone remodeling. Thus, given cell population data, this model can be used, with parameter fitting, to determine whether osteocytes or pre-osteoblasts provide the dominant RANKL source during targeted bone remodeling. This allows for more thorough testing of hypotheses such as in [48], which suggests that RANKL produced by osteoblast precursors contribute little to adult bone remodeling.

[Figure omitted. See PDF.]

Table 3. Effectiveness of RANKL Expression.

https://doi.org/10.1371/journal.pone.0063884.t003

Signaling of Osteocytes and Pre-osteoblasts

In this section we explore the behavior of the system as it is affected by changes in the effectiveness of osteocyte paracrine, and pre-osteoblast autocrine signaling. In particular, we regard how these actions influence the differentiation of stromal cells into pre-osteoblasts, and the proliferation of pre-osteoblast cells. These aspects of BMU remodeling are strongly influenced by sclerostin signaling. While previous mathematical models have included pre-osteoblast proliferation, explicit inclusion of osteocyte signaling, and the effects sclerostin, has been outside the scope of those works and thus is another feature of novelty in this work.

We compute the steady state bone volume, , as a simultaneous function of the effectiveness of osteocyte paracrine signaling on stromal cell differentiation, , and pre-osteoblast autocrine signaling for pre-osteoblast proliferation, . The results are shown in Figure 5, Figure 6, and Figure 7. The reason to consider these parameters simultaneously is that these terms both influence the “production” of pre-osteoblasts in (2). First, we observe that the parameter is considerably more sensitive to variation than , this is seen in Figure 5 and in a comparison between Figure 6 and Figure 7. In the context of the model, this suggests that the influence of osteocyte signaling, and sclerostin specifically, on stromal cell differentiation into pre-osteoblasts is key to the model. In fact, we can obtain a steady state value for bone volume, , corresponding to normal bone remodeling for a relatively wide range of values, provided is at or near baseline value. The significance of the role of osteocyte signaling is particularly strong with regard to the initiation of a cycle of targeted remodeling, since as pointed out above, the steady state values for all cell populations, except osteocytes, is zero. That is, we do not have a presence of cells fully committed to either the osteoblast or osteoclast lineage when there is no active remodeling. This is a point of departure with previous mathematical models where positive steady state values of mature osteoblasts and osteoclasts occur.

[Figure omitted. See PDF.]

Figure 5. The steady state bone volume, , as a simultaneous function of the effectiveness of osteocyte paracrine signaling on stromal cell differentiation, , and pre-osteoblast autocrine signaling for pre-osteoblast proliferation, .

https://doi.org/10.1371/journal.pone.0063884.g005

[Figure omitted. See PDF.]

Figure 6. The steady state bone volume, , computed as a function of the effectiveness of osteocyte paracrine signaling on stromal cell differentiation, , with all other parameters held at baseline values.

https://doi.org/10.1371/journal.pone.0063884.g006

[Figure omitted. See PDF.]

Figure 7. The steady state bone volume, , computed as a function of the effectiveness of pre-osteoblast autocrine signaling for pre-osteoblast proliferation, , with all other parameters held at baseline values.

https://doi.org/10.1371/journal.pone.0063884.g007

Application to Study Anti-sclerostin Drugs

We can apply the model system to study the effects of using an anti-sclerostin drug. Many drugs for treating bone disease, such as biphosphonates, are designed to inhibit resorption, and to decrease the loss of bone mass. On the other hand, anti-sclerostin drugs purport to promote bone formation, thereby increasing bone mass [53], [54]. Sclerostin inhibits Wnt/-catenin, so anti-sclerostin would promote osteoblast differentiation. In addition, Wnt signaling increases OPG expression, which inhibits osteoclasts, so anti-sclerostin would be expected to increase OPG release too. Sclerostin has been shown to downregulate OPG and increase levels of RANKL. By modifying the parameters and which represent the effectiveness of sclerostin regulation of osteoblastogenesis and osteoclastogenesis respectively, albeit in an abstract manner, we can study the effects of an anti-sclerostin agent on bone remodeling in a situation where there is known to be an imbalance in resorption/formation (e.g., abnormal osteoclast activity).

To model the inclusion of an anti-sclersotin drug we proceed as follows. In general, the parameters and now become time dependent, that is and . The form of the time dependence is determined by the frequency of dosing, the half life of the drug, and potentially the method of drug delivery. We note that since and are exponents in a power law approximation they must be dimensionless quantities, this remains the case even when they are made to depend on time. We take the function to be a constant value over a specific time interval to , and the baseline value otherwise. We take the function to be a constant over another time interval to , and the baseline value otherwise. The constants an represent the perturbation of the baseline parameter values for and respectively. In order to get an increase in bone formation we must have be less than the baseline value, and be greater than the baseline value. The time intervals to , and to represent the period of activity of the given dose of the drug. The time interval to is the time over which the drug effects Wnt signaling. The time interval to is the time over which the drug effects OPG signaling. We note that these two time intervals may be the same.

Figure 8 shows the change in bone mass as a function of time in the case of simulated pathological remodeling. Here we simulate the results of a hypothetical bone disease that results in over resorption. This is modeled by increasing the value of from its baseline value as listed in table 2 to a value of . Figure 8 shows the bone mass as a function of time under conditions where there is over resorption and no treatment, that is all parameters other than are set to baseline values, and, in particular we do not yet modify and .

[Figure omitted. See PDF.]

Figure 8. We simulate the loss of bone mass associated with over resorption in conjunction with a bone degenerative disease, modeled by a slight increase in the value of with all other parameters set to the baseline values listed in table 2.

https://doi.org/10.1371/journal.pone.0063884.g008

Next, we simulate the treatment of the pathological remodeling, such as is show in Figure 8, via an anti-sclerostin agent. Figure 9 shows the results obtained by including treatment with an anti-sclerostin drug, implemented by modifying the parameters and as described above. Specifically, all the parameters are the same as those used to obtain the results show in 8, but now we modify and to be the functions and , previously defined, to model the effects of treatment with an anti-sclerostin antibody. We observe a dose dependent increase in bone mass as described in [53], [54].

[Figure omitted. See PDF.]

Figure 9. This figure shows simulation results of treating pathological bone remodeling, as simulated in figure 5, via the addition of an anti-sclerostin drug.

This results in a dose-dependent increase in bone mass. Treatment with an anti-sclerostin drug is modeled by modifying the appropriate signaling mechanisms, that is, by modifying the appropriate exponents and in the power law approximation.

https://doi.org/10.1371/journal.pone.0063884.g009

There are several notable features regarding the treatment of bone loss using an anti-sclerostin antibody within the context of this model. One is that the time of the dosing matters. The closer to initial stages of osteoclastogenesis that the dose is given, the better the response of increased bone formation. This is particularly the case with respect to OPG signaling related parameter . Another feature we observed is that the parameter , which corresponds primarily to Wnt signaling, is sensitive to the perturbations that model dosing with an anti-sclerostin drug. That is, it takes only a small decrease in the value over a short period of time for a marked increase in bone formation.

Discussion

We have developed a cell population mathematical model in the form of a system of nonlinear ordinary differential equations for the principal bone cells involved in local targeted bone remodeling. The primary goal for this work was to include features of major focus in some of the most recent experimental research on bone remodeling at the cellular level. In particular we have aimed to present a model that extends and compliments current mathematical models and allows for a theoretical testing of hypotheses put forth in recent experimental work that is difficult to test in vivo or in clinical settings.

The most prominent consideration in the development of this model, that differs significantly from previous work, has been the explicit inclusion of osteocyte cell populations, and the effects of osteocyte signaling on bone turnover during targeted bone remodeling. First, in contrast to previous mathematical models that required perturbation away from a positive steady state to begin a remodeling cycle, our model allows for a more natural (physiological) initiation of a cycle of targeted bone remodeling. Second, the steady state values for forming or resorbing bone cells is in better agreement with what is physiologically expected, since the presence of mature osteoclasts and osteoblasts is not expected in the absence of a remodeling BMU. Finally, the approach taken here allows for a natural representation of the action of sclerostin which inhibits Wnt signaling and pre-osteoblast proliferation, and may also promote osteoclast activity.

We have used, following [28], the steady state bone volume as a control to examine aspects of targeted bone remodeling such as osteocyte RANKL expression, initiation of targeted bone remodeling with osteocyte apoptosis, and the role of sclerostin that are of current interest to the larger bone remodeling community. We have shown that this model is appropriate for determining whether osteocyte or pre-osteoblast RANKL expression plays the more prominent role in osteoclastogenesis for adult bone remodeling. Furthermore, we have applied this model to study the effects of an anti-sclerostin drug, such as an anti-sclerostin monoclonal antibody described in [53], on bone turnover.

One important feature of bone remodeling that is beyond the scope of the model presented in this work is spatial bone remodeling. We point out that the main thrust of spatially explicit mathematical models of bone remodeling have been spatial extensions of purely temporal models, see e.g., [1], [15]. That is, the ordinary differential equations used to model local bone remodeling are extended to partial differential equations.

However, a different framework for spatial extensions based on level set methods (LSM) is also possible, see e.g., [20], [55]. The models used in [20], [55] did not include the role of osteocytes, or differentiate between precursor cells and active osteoblasts and osteoclasts. Nonetheless the works [20], [55] provide a proof of concept for a new spatial representation of bone remodeling, and demonstrate the sort of pitfalls that might occur when spatially extending the model to scales of interest for bone marrow biopsy cores. Quite importantly, we did not see a change in the local behavior based on embedding into the LSM spatial representation; in fact local behavior translated into spatial behavior well.

Another feature, somewhat related to spatial behavior and not considered in this work are the mechanical mechanisms involved in the initiation and regulation of bone remodeling. This is primarily due to the fact that this work considers local remodeling at the scale of an individual BMU but without treating the problem of BMU steering.

Decades of experimental studies on bone remodeling reveal a dynamic process that is sensitive to multiple local and humoral factors. In silico modeling based on this accumulated knowledge offers a means to assess systematically the effects of perturbing individual parts of the system at a rate that cannot be matched by experiments conducted in vivo. As such, model-based predictions are particularly useful for the development of new treatments for bone disorders. This is somewhat illustrated by results from our model, which underscore the decisive influence of osteocytes and the potential for sclerostin-targeted treatments to reduce bone loss in pathologic conditions ranging from osteoporosis to osteolytic cancers.

It has been noted that dysfunctional or apoptotic osteocytes contribute to post-menopausal, steroid-, and immobilization-induced osteoporosis (reviewed by [56]). High circulating levels of sclerostin have been associated with osteoporosis, Paget's disease, and metastatic bone disease [57], [58]. Intriguingly, sclerostin has also been found to be expressed by malignant plasma cells in multiple myeloma [59]. Agents targeting sclerostin are under development, representing a novel therapeutic strategy compared with such currently available options as the bisphosphonates which directly disrupt osteoclast function. Disorders of bone remodeling affect millions of individuals and are responsible for considerable morbidity and mortality, resulting billions of dollars in healthcare costs each year. Given the complexity of the cellular and biochemical processes which mediate bone remodeling, the availability of a mathematical model which can capture all of these elements is an exceedingly useful tool for understanding bone remodeling in the setting of disease. Our model, which, to our knowledge, is the first to incorporate osteocytes and sclerostin, can be used to answer questions regarding the relative impact of cell type-specific activity. Not only can the effects of existing therapeutic agents be modeled, but this model could be used to explore the effects of modulating the activity of new targets thus providing rationale for new drug development. Furthermore, this model could be extended to multiple myeloma by the incorporation of plasma cells and the resulting interactions between the bone and stromal cells.

Author Contributions

Contributed reagents/materials/analysis tools: JMG BPA SAH JAM. Wrote the paper: JMG BPA SAH JAM.

Citation: Graham JM, Ayati BP, Holstein SA, Martin JA (2013) The Role of Osteocytes in Targeted Bone Remodeling: A Mathematical Model. PLoS ONE 8(5): e63884. https://doi.org/10.1371/journal.pone.0063884

References

1. Ryser M, Komarova S, Nigam N (2010) The cellular dynamics of bone remodeling: a mathematical model. SIAM J Appl Math 70: 1899–1921.

2. Hofbauer L, Kuehne C, Viereck V (2004) The opg/rankl/rank system in metabolic bone diseases. J Musculoskelet Neuronal Interact 4: 268–275.

3. Hofbauer L, Schoppet M (2004) Clinical implications of the osteoprotegerin/rankl/rank system for bone and vascular diseases. JAMA 292: 490–495.

4. Burr D (2002) Targeted and nontargeted remodeling. Bone 30: 2–4.

5. Parfitt A (2002) Targeted and nontargeted bone remodeling: relationship to basic multicellular unit origination and progression. Bone 30: 5–7.

6. Bonewald L (2011) The amazing osteocyte. J Bone Miner Res 26: 229–238.

7. Neve A, Corrado A, Cantatore F (2012) Osteocytes: central conductors of bone biology in normal and pathological conditions. Acta Physiol 204: 317–330.

8. Kular J, Tickner J, Chim S, Xu J (2012) An overview of the regulation of bone remodelling at the cellular level. Clin Biochem In Press.

9. Heino T, Kurata K, Higaki H, Väänänen H (2009) Evidence for the role of osteocytes in the initiation of targeted remodeling. Technol Health Care 17: 49–56.

10. Geris L, Gerisch A, Maes C, Carmeliet G, Weiner R, et al. (2006) Mathematical modeling of fracture healing in mice: comparison between experimental data and numerical simulation results. Med Biol Eng Comput 44: 280–289.

11. Geris L, Gerisch A, Schugart R (2010) Mathematical modeling in wound healing, bone regeneration and tissue engineering. Acta Biotheor 58: 1–13.

12. Geris L, Gerisch A, Sloten JV, Weiner R, Oosterwyck HV (2008) Angiogenesis in bone fracture healing: a bioregulatory model. J Theor Bio 25: 137–158.

13. Geris L, Reed AA, Vander Sloten J, Simpson AH, Van Oosterwyck, et al. (2010) Occurrence and treatment of bone atrophic non-unions investigated by an integrative approach. PLoS Comput Bio 6: 189–193.

14. Geris L, Vander Sloten J, Van Oosterwyck H (2010) Connecting biology and mechanics in fracture healing: an integrated mathematical modeling framework for the study of nonunions. Biomech Model Mechanobiol 9: 1–12.

15. Ayati B, Edwards C, Webb G, Wikswo J (2010) A mathematical model of bone remodeling dynamics for normal bone cell populations and myeloma bone disease. Biology Direct 5.

16. Buenzli P, Jeon J, Pivonka P, Smith D, Cummings P (2012) Investigation of bone resorption within a cortical basic multicellular unit using a lattice-based computational model. Bone 50: 378–389.

17. Buenzli P, Pivonka P, Gardiner B, Smith D (2012) Modelling the anabolic response of bone using a cell population model. J Theor Biol 307C: 42–52.

18. Buenzli P, Pivonka P, Smith D (2011) Spatio-temporal structure of cell distribution in cortical bone multicellular units: a mathematical model. Bone 48: 918–926.

19. Defranoux N, Stokes C, Young D, Kahn A (2005) In silico modeling and simulation of bone biology: a proposal. J Bone Miner Res 20: 1079–1084.

20. Graham J, Ayati B, Ramakrishnan P, Martin J (2012) Towards a new spatial representation of bone remodeling. Math Biosci Eng 9: 281–295.

21. Komarova S (2005) Mathematical model of paracrine interactions between osteoclasts and osteoblasts predicts anabolic action of parathyroid hormone on bone. J Endocrinol 146: 3589–3595.

22. Komarova S, Smith R, Dixon S, Sims S, Wahl L (2003) Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone remodeling. Bone 33: 225–234.

23. Lemaire V, Tobin F, Greller L, Cho C, Suva L (2004) Modeling the interactions between osteoblast and osteoclast activities in bone remodeling. J Theor Biol 229: 293–309.

24. Liotier P, Rossi J, Wendling-Mansuy S, Chabrand P (2012) Trabecular bone remodelling under pathological conditions based on biochemical and mechanical processes involved in bmu activity. Comput Methods Biomech Biomed Engin In Press.

25. Martin M, Buckland-Wright J (2004) Sensitivity analysis of a novel mathematical model identifies factors determining bone resorption rates. Bone 35: 918–928.

26. Martin M, Buckland-Wright J (2005) A novel mathematical model identifies potential factors regulating bone apposition. Calcif Tissue Int 77: 250–260.

27. Pivonka P, Komarova S (2010) Mathematical modeling in bone biology: From intracellular signaling to tissue mechanics. Bone 47: 181–189.

28. Pivonka P, Zimak J, Smith D, Gardiner B, Dunstan C, et al. (2008) Model structure and control of bone remodeling: A theoretical study. Bone 43: 249–263.

29. Pivonka P, Zimak J, Smith D, Gardiner B, Dunstan C, et al. (2010) Theoretical investigation of the role of the rank-rankl-opg system in bone remodeling. J Theor Biol 262: 306–316.

30. Ryser M, Nigam N, Komarova S (2009) Mathematical modeling of spatio-temporal dyanmics of a single bone multicellular unit. J Bone Miner Res 24: 860–870.

31. Zumsande M, Stiefs D, Siegmund S, Gross T (2011) General analysis of mathematical models for bone remodeling. Bone 48: 910–917.

32. Savageau M (1969) Biochemical systems analysis. i. some mathematical properties of the rate law for the component enzymatic reactions. J Theor Biol 25: 365–369.

33. Savageau M (1969) Biochemical systems analysis. ii. the steady-state solutions for an n-pool system using a power-law approximation. J Theor Biol 25: 370–379.

34. Savageau M (1970) Biochemical systems analysis. iii. dynamic solutions using a power-law approximation. J Theor Biol 26: 215–226.

35. Savageau M (1976) Biochemical systems analysis. Reading, MA: Addison-Wesley.

36. Voit E (2000) Computational analysis of biochemical systems. Cambridge, UK: Cambridge University Press.

37. Lio P, Paoletti N, Atwell M, Merelli E, Viceconti M (2012) Modelling osteomyelitis. BCM Bioinformatics 13 Suppl 14.

38. Bilezikian J, Raisz L, Rodan G (2002) Principles of bone biology. Boston: Academic Press, second edition.

39. Bell N (2003) Rankl and the regulation of skeletal remodeling. J Clin Invest 111: 1120–1122.

40. Hofbauer L, Heufelder A (2001) Role of receptor activator of nuclear factor-b ligand and osteoprotegerin in bone cell biology. J Mol Med 79: 243–253.

41. Martin T (2004) Paracrine regulation of osteoclast formation and activity: milestones in discovery. J Musculoskelet Neuronal Interact 4: 243–253.

42. Simonet W, Lacey D, Dunstan C, Kelley M, Chang M, et al. (1997) Osteoprotegerin: a novel secreted protein involved in the regulation of bone density. Cell 89: 309–319.

43. Cohen M (1997) Transforming growth factor beta's and fibroblast growth factors and their receptors: role in sutural biology and craniosynostosis. J Bone Miner Res 12: 322–331.

44. Janssens K, ten Dijke P, Janssens S, Hul WV (2005) Transforming growth factor-beta1 to the bone. Endocr Rev 26: 743–744.

45. Boyle W, Simonset W, Lacey D (2003) Osteoclast differentiation and activation. Nature 423: 337.

46. Gori F, LC Hofbauer, Dunstan C, Spelsberg T, et al. (2000) The expression of osteoprotegerin and rank ligand and the support of osteoclast formation by stromal-osteoblast lineage cells is developmentally regulated. Endocrinology 141: 4768–76.

47. Cardoso L, Herman B, Verborgt O, Laudier D, Majeska R, et al. (2009) Osteocyte apoptosis controls activation of intracortical resoprtion in response to bone fatique. J Bone Min Res 24: 597–605.

48. Xiong J, Onal M, Jilka R, Weinstein R, Manolagas S, et al. (2002) Matrix-embedded cells control osteoclast formation. Nature Med 17: 1235–1242.

49. Shampine L, Reichelt M (1997) The MATLAB ODE suite. SIAM J Sci Comput 18: 1–22.

50. Shampine LF, Gladwell I, Thompson S (2003) Solving ODEs with MATLAB. Cambridge: Cambridge University Press, viii+263 pp. doi:10.1017/CBO9780511615542. Available: http://dx.doi.org/10.1017/CBO9780511615542.

51. Shampine L (1994) Numerical solution of ordinary differential equations. New York: Chapman & Hall, x+484 pp.

52. Giuliani N, Ferretti M, Bolzoni M, Storti P, Lazzaretti M, et al. (2012) Increased osteocyte death in multiple myeloma patients: role in myeloma-induced osteoclast formation. Leukemia 26: 1391–1401.

53. Lewiecki E (2011) Sclerostin: a novel target for intervention in the treatment of osteoporosis. Discov Med 12: 263–273.

54. Lim V, Clarke B (2012) New therapeutic targets for osteoporosis: Beyond denosumab. Maturitas 73: 269–272.

55. Ayati B, Graham J, Holstein S (2012) Using a level set to model multiple myeloma induced bone loss. Proceedings of the 6th European Congress on Computational Methods in Applied Sciences and Engineering.

56. Manolagas S, Parfitt A (2012) For whom the bell tolls: Distress signals from long-lived osteocytes and the pathogenesis of metabolic bone diseases. Bone.

57. Ardawi M, Rouzi A, Al-Sibiani S, Al-Senani N, Qari M, et al. (2012) High serum sclerostin predicts the occurrence of osteoporotic fractures in postmenopausal women: The center of excellence for osteoporosis research study. J Bone Miner Res 27: 2592–2602.

58. Yavropoulou M, va Lierop A, Handy N (2012) Serum sclerostin levels in paget's disease and prostate cancer with bone metastases with a wide range of bone turnover. Bone 51: 153–157.

59. Brunetti G, Oranger A, Mori G, Specchia G, Rinaldi E, et al. (2011) Sclerostin is overexpressed by plasma cells from multiple myeloma patients. Ann NY Acad Sci 1237: 19–23.