A mathematical model of bone remodeling dynamics for normal bone cell populations and myeloma bone disease
© Ayati et al. 2010
Received: 14 January 2010
Accepted: 20 April 2010
Published: 20 April 2010
Skip to main content
© Ayati et al. 2010
Received: 14 January 2010
Accepted: 20 April 2010
Published: 20 April 2010
Multiple myeloma is a hematologic malignancy associated with the development of a destructive osteolytic bone disease.
Mathematical models are developed for normal bone remodeling and for the dysregulated bone remodeling that occurs in myeloma bone disease. The models examine the critical signaling between osteoclasts (bone resorption) and osteoblasts (bone formation). The interactions of osteoclasts and osteoblasts are modeled as a system of differential equations for these cell populations, which exhibit stable oscillations in the normal case and unstable oscillations in the myeloma case. In the case of untreated myeloma, osteoclasts increase and osteoblasts decrease, with net bone loss as the tumor grows. The therapeutic effects of targeting both myeloma cells and cells of the bone marrow microenvironment on these dynamics are examined.
The current model accurately reflects myeloma bone disease and illustrates how treatment approaches may be investigated using such computational approaches.
This article was reviewed by Ariosto Silva and Mark P. Little.
Bone is continually renewed throughout the skeleton in a process known as remodeling. The bone remodeling process is spatially heterogeneous, with regular but asynchronous cycles at multiple sites that can occupy 5-25% of bone surface . In this way every part of the skeleton is remodeled periodically over time. Mathematical modeling of bone remodeling has focused on various aspects of this process. These approaches include models of how Michaelis-Menten-like feedback mechanics affect bone resorption (the process by which osteoclasts break down bone, resulting in bone loss)  and how biomechanical stress induces bone formation [3–6]. Other modeling efforts have examined the signaling pathways between osteoclasts and osteoblasts involved in bone remodeling , or accounted for the activity of both osteoclasts and osteoblasts in a microenvironment known as a basic multicellular unit (BMU) [8–12].
In this paper we develop mathematical models of myeloma bone disease. Multiple myeloma is a hematological malignancy associated with clonal expansion of malignant plasma cells within the bone marrow. One of the major clinical features of myeloma is the development of a progressive and destructive osteolytic bone disease, associated with severe bone pain, pathological fractures, osteoporosis, hypercalcemia and spinal cord compression. Interactions between myeloma cells and cells of the bone marrow microenvironment are critical for myeloma growth and survival and for the development of the osteolytic bone disease [13, 14]. The destructive nature of myeloma bone disease is increased by the vicious cycle that develops between myeloma cells and the bone marrow microenvironment. Although the precise molecular mechanisms responsible for the bone destruction in multiple myeloma are not completely understood, it is known that the bone destruction is primarily mediated by osteoclasts, and that this destruction is exacerbated by a reduction in osteoblastic bone formation. Patients with multiple myeloma have abnormal bone remodeling, where resorption and formation become uncoupled, with the end result being an increase in bone resorption and a decrease in bone formation. Myeloma cells are found in close association with sites of active bone resorption, and their ability to stimulate osteoclast formation and activity has been well characterized [14–17]. Histological studies have demonstrated that in the early stages of myeloma, bone formation is actually increased, and this is thought to reflect the attempt to compensate for the increase in osteoclastic resorption . However, as the disease progresses, bone formation is rapidly decreased [15, 18, 19]. This has been confirmed in studies which demonstrate that markers of bone formation are decreased in patients with multiple myeloma [20, 21]. Despite many significant advances in the understanding of the biology of multiple myeloma, it remains an incurable malignancy, and the destructive osteolytic bone disease is a major cause of morbidity in patients with multiple myeloma.
The organization of the paper is as follows: In Section 2 we first present the tumor-free zero-dimensional model in  and illustrate its simulation of normal bone modeling dynamics. In Section 3 we add the tumor cell population to the model of normal bone in Section 2, and in Section 4 we simulate treatment for the model in Section 3. In Section 5 we add an explicit spatial dimension to the normal bone model. This additional dimension allows for heterogeneity on the surface of the bone along one axis. In Sections 6 and 7 we add tumor and treatment to the spatial model. In Section 8 we conclude with a discussion of the results. Computations were conducted using the Mathematica function NDSolver. All models in this paper use dimensionless variables and parameters, including for the cell populations.
with initial conditions C(0) = C 0, B(0) = B 0. The power law nonlinearities in (1a) and (1b) are approximations for the interactions of the osteoclast and osteoblast populations in the proliferation terms of the equations. In (1a) autocrine signaling has a positive feedback on osteoclast production (g 11 > 0), and paracrine signaling has a negative feedback on osteoclast production (g 21 < 0). In (1b) autocrine signaling has a positive feedback on osteoblast production (g 22 > 0), and paracrine signaling has a positive feedback on osteoblast production (g 12 > 0).
The solutions exhibit limit cycles as Ψ passes through 0: Ψ < 0 yields damped oscillations converging to the nontrivial steady state , , and Ψ > 0 yields unstable oscillations diverging away from the nontrivial steady state (where Ψ is sufficiently close to 0).
is the period of the cycles of C(t) and B(t), and r is determined by the amplitude of the oscillations in the bone mass. The value R is well-defined as long as B(0) ≠ .
with the initial conditions C(0) = C 0, B(0) = B 0, and T (0) = T 0. The equation for bone mass z(t) is (5), as before. The tumor equation (7c) is of Gompertz form with growth constant γ T > 0 and maximum tumor size L T . We have taken γ T to be independent of bone loss. Future work will include models, simulations, and biological experiments to determine the dependence of γ T on bone loss. In (7) the tumor parameters r 11, r 12, r 21, r 22 are all nonnegative.
The tumor presence alters (1a) as follows: autocrine promotion of osteoclasts is increased ( , since g 11 > 0); paracrine inhibition of osteoclasts is reduced ( , since g 21 < 0); paracrine promotion of osteoblasts is reduced ( , since g 12 > 0); and autocrine promotion of osteoblasts is reduced ( , since g 22 > 0). The key difference between our model and that of  in (1) is the addition of the terms r ij T (t)/L T that couple the tumor density and maximum size to the power laws for the osteoclast/osteoblast interactions.
If Φ < 0, then the system exhibits damped oscillations converging to the nontrivial steady state, and if Φ > 0, then the system exhibits unstable oscillations converging away from the nontrivial steady state (where Φ is sufficiently close to 0).
Remark 1A similar analysis can be made of other model parameters, as well as other treatment functions.
The tumor is introduced at time 0, whereas treatment with a proteasome inhibitor is started at t start = 600. At this time, there is an increase in tumor mass and osteoclast number, and a decrease in osteoblast number, indicating the development of myeloma bone disease. Treatment with a proteasome inhibitor decreases tumor burden from the time of treatment, whereas there is a delay in the recovery of the bone mass. Analysis of the individual cell types (osteoclasts and osteoblasts) indicates a decrease in osteoclasts from time of treatment, but a delay in the expected increase in osteoblasts. This suggests that the reduction in osteoclast number is dependent on tumor burden, but is not sufficient to increase bone mass. The increase in bone mass appears to parallel the increase in osteoblast number. Taken together, this model suggests that proteasome inhibition has a dramatic effect to reduce tumor burden, and to prevent bone loss as well as increase bone formation to steady-state levels in multiple myeloma. These results are consistent with observations from both clinical studies and preclinical murine models, where tumor burden is decreased and bone volume or markers of bone formation are increased in response to bortezomib treatment [26, 30–32]. Furthermore, by enabling analysis of tumor burden, bone volume and cell number over time, the results from this modeling provide important insights into the biology behind the clinical response to bortezomib which cannot be obtained from current in vivo studies.
where σ3 is the diffusion coefficient for the bone mass, and k 1(x) and k 2(x) depend on C(0, x) and B(0, x). We remark that σ3 is typically small and represents stochasticity in the bone dynamics, not actual migration of bone stromal cells.
We give examples of the normal bone model for two cases of initial spatial inhomogeneity in the osteoclast population C(0, x). In both examples the nonspatial parameters are α 1 = 3.0, α 2 = 4.0, β 1 = 0.2, β 2 = .02, g 11 = 1.1, g 22 = 0.0, g 12 = 1.0, g 21 = -0.5 (as in Fig. 3, Part I), the spatial parameters are σ1 = .000001, σ2 = .000001, and the bone mass parameters are σ 3 = .000001, k 1(x) ≡ 0.45, k 2(x) ≡ .0048, with z(0, x) ≡ 100.0. For these parameters the constant functions (x) ≡ 1.16, (x) ≡ 231.72 are steady-state solutions (see Fig. 3, Part I).
For the spatial model with tumor we assume the tumor cells are also diffusing in Ω. The variables of the model are C(t, x) and B(t, x) as before, and T(t, x) = density of tumor cells at time t with respect to x ∈ Ω.
and initial conditions C(0, x) = C 0(x), B(0, x) = B 0(x), T (0, x) = T 0(x). The bone mass equations for z(t, x) are as in (14e) and (14f). The diffusion coefficient for the tumor is σ4, which allows for spatial growth.
and initial conditions C(0, x) = C 0(x), B(0, x) = B 0(x), T (0, x) = T 0(x). The bone mass equations for z(t, x) are again as in (14e) and (14f).
We have presented a dynamic model of spatially heterogeneous bone remodeling that incorporates the interaction of osteoclasts (bone resorption) and osteoblasts (bone formation) subject to myeloma bone disease and its treatment. Our model is a system of nonlinear partial differential equations for osteoclast-osteoblast interactions driven by autocrine-paracrine signaling, which allows for interpretation of the corresponding spatial changes in bone mass and tumor growth. In the case of normal bone the model views remodeling as stable regular oscillations in a spatially heterogeneous microenvironment. In the case of myeloma, these regular cycles are destabilized with an increase in osteoclasts typical of this disease, and with corresponding destruction of bone mass and progressive tumor growth. Furthermore, additional features of myeloma bone disease are reproduced by this model, including the initial increase in osteoblast number, which is thought to reflect the attempt to maintain normal coupling of bone resorption to bone formation . This is followed by a decrease in osteoblast oscillation below that of the non-tumor simulation, reflecting the decrease in osteoblasts which leads to reduced bone formation and subsequent bone loss in multiple myeloma. Of interest is the observation that osteoclast oscillations remain elevated above the non-tumor simulation at all times, despite an overall decrease reflecting the gradual bone loss. Again, this is highly representative of myeloma bone disease, which is associated with an increase in both osteoclast number and activity. We incorporated treatment into the model with a drug, such as a proteasome inhibitor, which promotes osteoblast production and inhibits tumor growth.
Proteasome inhibitors were initially identified by their dramatic effects to reduce myeloma tumor burden [28–30, 32]. Subsequent studies suggested that this class of drugs may have additional effects in myeloma, due to their direct effects to promote osteoblast differentiation and subsequent bone formation [26–30]. In addition, it has been suggested that proteasome inhibitors may have effects on osteoclastic bone resorption [31, 33, 34]. In the current simulation, treatment with a proteasome inhibitor which had direct effects on myeloma cells and osteoblast formation was found to significantly reduce tumor burden and prevent myeloma bone disease, in agreement with reported in vivo preclinical studies. Of interest, effects on osteoclast number were also observed, suggesting that proteasome inhibitors may have indirect effects on osteoclasts. Osteoclasts and tumor burden were reduced from time of treatment, whereas there was a delay in the increase in osteoblasts and bone mass. This suggests firstly that the reduction in osteoclast number is dependent upon the tumor burden, and secondly that the increase in bone mass is a result of the increase in osteoblasts. Such insights into potential cellular mechanisms cannot easily be gained from in vivo experiments where cellular effects can typically only be assessed at endpoint, thus highlighting the value of this mathematical model. Furthermore, although the simulated proteasome inhibitor treatment was not able to completely restore bone mass, this is most likely a reflection of the rates chosen for the model and it is expected that bone mass would recover over a longer time period. This may suggest that the effects of proteasome inhibition to replace bone already lost in myeloma may take longer than the effects to reduce tumor burden, and have implications for duration of treatment.
The model we have developed is one-dimensional in space, and is thus idealized with respect to the geometry of the bone microenvironment. In future work we will develop a more realistic higher dimensional description of the bone microenvironment, and the dynamics of bone myeloma and its treatment in this spatial microenvironment. The major problem in the identification of therapeutic approaches for the treatment of multiple myeloma is the complex relationship between the multiple types of cells distributed throughout the bone marrow microenvironment, which is almost impossible to recreate in vitro. Therefore, in order to determine the effect of a drug, the optimal dose, treatment regimen or combination, and relative timing of the delivery of drug combinations, in vivo studies are required. Although the 5T murine model of myeloma is very effective and reliable for such preclinical studies, it is impossible to evaluate all drugs or combinations in vivo. Those that are selected for in vivo evaluation are often chosen based upon in vitro assays, which do not reflect the interactions within the bone marrow microenvironment.
The frequencies of the osteoblast/osteoclast oscillations in the mouse are not yet known. Experimental determination of local trabecular density and osteoclast and osteoblast distributions in the 5T murine model requires histological examination of fixed sections of bone, precluding serial studies of the progression of multiple myeloma in an individual animal. The mathematical model should prove extremely useful in both minimizing the number of animals that have to be sacrificed to obtain statistically significant data regarding the time course of the disease, and in rationalizing inter-animal differences. Furthermore, the coupling between osteoblasts and osteoclasts, and bone remodeling cycle has yet to be recreated in vitro, either for normal bone or for diseased bone in the presence of myeloma cells. The current model accurately reflects myeloma bone disease and illustrates how treatment approaches may be investigated using such computational systems. Our ultimate goal is to quantify, with experimental support, the dynamics of bone remodeling in health and disease, and gain insight into the design and optimization of new therapeutic approaches for the treatment of myeloma bone disease.
The mathematical analyses in this paper use accepted methodology from the theories of partial differential equations and dynamical systems. The computations use standard algorithms as incorporated into the software Mathematica.
All authors contributed to the development and interpretation of the mathematical model and the design of the figures. BPA and GFW developed the mathematical parts of the manuscript and conducted the simulations. All authors contributed to writing the manuscript, which was approved by all authors.
I would suggest the authors to consider in their future additions to the model the effect of hypoxia in the bone formation equilibrium, mainly in the myelomatous bone marrow. Works have described the influence of hypoxia in the activity of osteoclasts  and pH in the activity of osteoblasts . Considering that the bone marrow vascularization is not evenly distributed as in other tissues, gradients of oxygen concentration (and supposedly also of pHe) could induce niches where bone degradation is more favorable. This effect could also be exacerbated by the proliferation of MM cells which are known to have increased glucose consumption, as seen in PET scans.
Response: We agree that hypoxia is an important aspect of the bone marrow microenvironment that will affect both bone cell activity and myeloma growth and survival, and as such, fully intend to incorporate hypoxia into future models.
In the conclusions, could you elaborate more on the reason of the increase in osteoblast activity in the early stages of MM?
Response: During normal bone remodeling, osteoclastic and osteoblastic activity is tightly coupled, with an increase in osteoclast activity followed by an increase in osteoblast activity. In the early stages of myeloma, myeloma cells within the bone marrow stimulate an increase in osteoclast activity. Osteoblast activity is also increased in an attempt to maintain the normal coupling of bone remodeling. As the disease progresses, osteoblast formation and activity are inhibited, leading to dramatic bone loss.
Does the delay after the reduction of tumor burden and recovery of bone mass depend on the parameters of the model? Is it found in vivo and in clinical treatment as well? Could the delay in recovery of bone mass in patients be mapped into parameters from the model? Could this possibly be used as a prognosis for relapse or patient recovery?
Response: The delay is dependent on model parameters. There is evidence to suggest that proteasome inhibition can increase bone formation in vivo, in addition to decrease tumor burden, however the time course of these events is unknown. This is an advantage of the mathematical model, that allows us to look in detail at all time points, whereas in vivo studies are limited to select time points which may not provide all information.
"Those that are selected for in vivo evaluation are often chosen based upon in vitro assays, which do not reflect the interactions within the bone marrow microenvironment". Some in vitro models have been used to assess the effect of fibronectin and contact with stroma cells in chemoresistance in leukemia  and Multiple Myeloma . It would be interesting indeed if data could be obtained from these in vitro models and ported into computational models like the one described in this work. Response: There is substantial evidence to demonstrate that interactions between bone marrow stromal cells and myeloma cells are important for promoting tumor growth and survival, both alone and in response to chemotherapeutic agents. Future studies will incorporate bone marrow stromal cells into the mathematical models.
"The mathematical model should prove extremely useful in both minimizing the number of animals that have to be sacrificed to obtain statistically significant data regarding the time course of the disease, and in rationalizing inter-animal differences". Depending on how you parameterize the model and which inputs you require from experiments, you may also created personalized models for each patient and even help on the prognosis of the disease or suggest more promising treatments for every patient.
Response: The long-term goal of this modeling is to create a mathematical model which closely resembles human myeloma, and with which we can explore multiple treatment protocols and select the one that has the optimal outcome. This would represent a major step forward in individualized cancer therapy - quantitative optimization of a combination of drugs with a possibly complicated delivery schedule, all based upon a validated mathematical model.
Even though this work is an extension of a previous one already published  it would be useful to describe the values used for the parameters, if they were estimated or obtained from literature and if the authors have plans on how to obtain these form experiments for the future versions of the model.
Response: The parameters for normal bone were taken from [1, 22]. The parameters for myeloma bone with and without treatment were estimated. In future work the parameters will be identified using experimental data.
This is a generally well-written paper, describing a schematic model of bone remodelling in the presence or absence of myeloma. It was not immediately obvious what the relevance of this article is to the Biology Direct special issue, but perhaps this doesn't matter too much. Arguably more important are the specific modelling assumptions, in particular the power law ODEs and PDEs relating osteoclasts and osteoblasts. While these appear plausible, a little more biological justification of these functional forms could be usefully provided here (e.g., as is given in ref. ).
Response: The theme of this Biology Direct special issue is Mathematics and the Evolution of Cancer. The manuscript outlines the development of a mathematical model of myeloma bone disease. Due to the interactions between myeloma cells and cells of the bone marrow microenvironment, the osteolytic bone disease associated with myeloma is inextricably linked with tumor progression, and therefore is relevant to this Biology Direct special issue. Power law nonlinearities have been used extensively in modelling biochemical kinetics in interacting systems. The parameters in the power law terms correspond to feedback and other regulatory mechanisms. In future work the power law parameters will be estimated using experimental data.
This work is supported by the National Science Foundation grant DMS-0914514 (BPA), the National Cancer Institute grant P01 CA-40035 (CME), the International Myeloma Foundation (CME), the Elsa U. Pardee Foundation (CME), the National Institutes of Health grant U54CA113007 (GFW, JPW), the Vanderbilt Institute for Integrative Biosystems Research and Education (JPW), and the Vanderbilt Integrative Cancer Biology Center (VICBC). The project was initiated during the 3rd Annual Workshop of the VICBC, July 15-19, 2007. We thank Allison Price for editorial assistance. As a one-time exception to the publishing policy of Biology Direct the articles in this series are being published with two reviewers.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.