Multiple myeloma is a hematologic malignancy associated with the development of a destructive osteolytic bone disease.

Results

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.

Conclusions

The current model accurately reflects myeloma bone disease and illustrates how treatment approaches may be investigated using such computational approaches.

Reviewers

This article was reviewed by Ariosto Silva and Mark P. Little.

Background

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 [1]. 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) [2] 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 [7], 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 [15]. 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.

Results

We model the influence of tumor growth on bone remodeling, and in particular how the tumor influences autocrine and paracrine signaling in the osteoclast and osteoblast cell populations (see Fig. 1). Autocrine signaling represents the feedback from osteoclasts and osteoblasts to regulate their respective formation. Paracrine signaling represents the factors produced by osteoclasts that regulate osteoblast formation, and vice versa. We use the underlying model of bone remodeling in the absence of tumor presented in [22] and explored further in [23–25]. This model is a dynamical system with zero explicit space dimensions, but with a dependent variable that records bone mass as a function of time. If we interpret the bone mass equation as one for localized trabecular mass (spongy bone found within the bone marrow) underneath a point on the surface of the bone, we obtain a representation of one spatial dimension. We then present one-dimensional spatial models that suggest how we may incorporate additional spatial dimensions.

The organization of the paper is as follows: In Section 2 we first present the tumor-free zero-dimensional model in [22] 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.

Zero-dimensional Bone Model without Tumor

A model of normal bone remodeling at a single discrete site is developed in [22]. The model consists of a system of ordinary differential equations describing the bone cell populations in a BMU. These populations are the osteoclasts, which resorb bone, and osteoblasts, which form bone. The variables of the model are the density of osteoclasts C(t) and the density of osteoblasts B(t) at time t. The equations of the model in [22] are

(1)

(2)

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 system (1) has a unique nontrivial steady state

(3)

(4)

where

(5)

The system (1) has periodic solutions when

(6)

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).

An additional variable z(t) for the bone mass is obtained in [22] by assuming bone mass is determined by the extent to which normalized values of C(t) and B(t) exceed nontrivial steady state levels. We reinterpret the variable z(t) as representing localized trabecular mass beneath a point on the bone surface. Having found the solutions of system (1), we define the following equation for the bone mass z(t) at time t:

(7)

with the initial condition z(0) = z_{0}. In the case that C(t) and B(t) have periodic solutions, k_{1} and k_{2} are chosen so that the normal bone mass oscillates with the same periodicity as the osteoclast and osteoblast populations with an experimentally determined amplitude of the oscillations. The assumption is that bone mass decreases or increases cyclically according to the net effects of resorption (C(t) >) and formation (B(t) >). The constants k_{1} and k_{2} satisfy k_{1} = rR, and k_{2} = r, where

(8)

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) ≠ .

We reproduce here two examples of single-site normal bone remodeling given in [22]. The parameters are as in [22], with time unit in days. Fig. 2 simulates a targeted event corresponding to an external stimulus, taken as a perturbation of the nontrivial steady state (, ) by a momentary increase in the number of osteoclasts. The result is a single remodeling cycle (Ψ < 0, but not sufficiently close to 0 to yield damped oscillations). Fig. 3 corresponds to a series of internally initiated regular cycles over an extended period of time, with osteoclast and osteoblast populations exhibiting regular periodic cycles about their nontrivial steady state values (Ψ = 0). The distinction of the two behaviors is the value of the osteoclast autocrine parameter g_{11}, which is greater in Fig. 3. In [22] it is claimed that this parameter is the primary factor in the regulation of bone remodeling dynamics.

Zero-dimensional Bone Model with Tumor

We find that the regular cycles of the normal bone model above are perturbed by the presence of myeloma. The tumor parameters yield perturbations of the normal cycles that result in limit cycles, which are either damped oscillations that converge to the nontrivial steady state or undamped oscillations that diverge away from the nontrivial steady state. The variables for the model with tumor are C(t), B(t) as before, and T(t) = density of tumor cells at time t. Our equations are

(9)

(10)

(11)

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 [22] 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.

The system (7) has nontrivial steady states

(12)

and

(13)

(14)

(15)

Where

(16)

The system (7) has stable cycles when

(17)

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).

A simulation of the bone model with tumor is given in Fig. 4 and Fig. 5 with Φ < 0. The osteoclast and osteoblast populations exhibit damped oscillations converging to the nontrivial steady state in Fig. 4, the tumor grows to maximum capacity, and the bone mass converges with oscillations to 0 (Fig. 5). We can see that as the tumor burden increases, there is an initial increase in both osteoclast and osteoblast number, reflecting the attempt of the system to maintain the normal coupling of bone resorption to bone formation even in the presence of tumor. Although the amplitude of osteoclast oscillations decreases over time, this reflects the decrease in bone mass. Importantly, the amplitude of osteoclast oscillations in the presence of tumor is higher than when tumor cells are not present (e.g., Fig. 3), reflecting the overall increase in osteoclasts associated with myeloma bone disease. In contrast, the amplitude of osteoblast oscillations is dramatically decreased when compared with the non-tumor simulation (Fig. 3), indicating the decrease in osteoblasts which is a characteristic feature of myeloma bone disease.

Another simulation of the bone model with tumor is given in Fig. 6 and Fig. 7 with the only change from Fig. 4 and Fig. 5 being that r_{11} = .02 instead of .005. In this case the osteoclast and osteoblast populations exhibit unstable oscillations (Φ > 0), whose amplitude grows with time, with a concomitant decrease in bone mass and increased growth of the tumor to its limiting size (Fig. 7).

As in [22] we can analyze the behavior of the model in terms of selected parameters. Here we investigate the dependence of the solution behavior on the parameters r_{11} and r_{22} and fix all the other parameters. The model has a unique nontrivial steady state as in (9a),(9b),(9c), which depends on r_{11} and r_{22}. The behavior of the solutions of the model as a function of r_{11} and r_{22} depends on the eigenvalues of the Jacobian at the nontrivial steady states (, ). If the maximum of the real parts of the eigenvalues is less than 0, then the solutions converge to the nontrivial steady state with damped oscillations, and if the maximum is greater than 0, then the solutions have unstable oscillations. These cases are illustrated in Fig. 8 as a function of r_{11} and r_{22}. The parameters r_{11} and r_{22} correspond to the alteration of the normal osteoclast-osteoblast regulation due to tumor burden. The relative values of r_{11} and r_{22} yield the two types of instability, that is, either decreasing amplitude oscillations or increasing amplitude oscillations, both departing from normal stable periodic oscillations.

Zero-dimensional Bone Model with Tumor and Drug Treatment

Many therapeutic approaches for the treatment of myeloma have the potential to affect directly both myeloma cells and cells of the bone marrow microenvironment, including osteoclasts and osteoblasts. Therefore, it is difficult to predict from in vivo studies the overall response to drug treatment in myeloma. We present a framework which provides us with the opportunity to model, for the first time, the effects of drug treatment on oscillatory bone remodeling. We have chosen to model the effects of proteasome inhibition in myeloma bone disease. Proteasome inhibitors are known to have direct anti-myeloma effects, and to have direct effects on osteoblasts to stimulate osteoblast differentiation and bone formation [26–30]. The equations are as before except that the treatment promotes osteoblast production and inhibits tumor growth:

(18)

(19)

(20)

The time-dependent treatment functions V_{1}(t) and V_{2}(t) in (12b) and (12c) are

(21)

(22)

(23)

(24)

where t_{start} is the starting time of treatment and v_{1} and v_{2} are the intensity parameters of treatment. A simulation of the effect of a proteasome inhibitor is given in Fig. 9 and Fig. 10, and corresponds to the untreated tumor simulations illustrated in Fig. 4 and Fig. 5. The tumor is extinguished, the osteoclast and osteoblast populations recover regular cycles, and the bone mass recovers to normal.

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.

One-dimensional Bone Model without Tumor

The normal bone model in [22] is a discrete site model for single event remodeling and internally regulated cycles of remodeling. We reinterpreted bone mass as implicitly providing one spatial dimension, such as trabecular mass beneath a point on the bone surface. To incorporate additional dimensions of spatial variability, we develop a diffusion model in a second spatial domain Ω. For convenience here we take a one-dimensional region Ω = [0, 1]. We assume that both osteoclasts and osteoblasts are diffusing in Ω. The variables of the model are C(t, x) = density of osteoclasts and B(t, x) = density of osteoblasts at time t with respect to x ∈ Ω. The equations are

(25)

(26)

with boundary conditions

(27)

(28)

and initial conditions C(0, x) = C_{0}(x) and B(0, x) = B_{0}(x) (see Fig. 11).

The nontrivial steady states and in the zero-dimensional case in Section are also steady states of the system (1) and viewed as constant functions and on Ω. Further, the change in bone mass z(t, x) as a function of x as well as t about a normalized value of 100 is given by the following equation: The value of 100 is arbitrary. Any other value, such as 1, could have also been used.

(29)

with initial condition z(0, x) = z_{0}(x), and boundary condition

(30)

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).

We first simulate the normal bone model with the initial distribution of osteoclasts C(0, x) elevated above = 1.16 near one site in Ω, as graphed in Fig. 12. The initial distribution of osteoblasts B(0, x) is taken as constant at = 231.72. The density plots of the osteoclasts and osteoblasts, and the change in bone mass, are graphed in Fig. 13, where it is seen that the solutions sustain spatial and temporal cycles about these steady states.

One-dimensional Bone Model with Tumor

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 ∈ Ω.

The equations are

(31)

(32)

(33)

with boundary conditions

(34)

(35)

(36)

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.

We illustrate the tumor model with an additional space dimension by adding the tumor population to the simulation in Fig. 12 and Fig. 13. Tumor parameters are γ_{
T
}= .004, L_{
T
}= 100, r_{11} = .005, r_{21} = 0.0, r_{12} = 0.0, r_{22} = 0.2 (as in Fig. 3), the additional spatial parameters are σ_{1} = .000001, σ_{2} = .000001, σ_{4} = .000001, and the bone mass parameters are σ_{3} = .000001, k_{1}(x) ≡ 0.45, k_{2}(x) ≡ .0048, with z(0, x) ≡ 100.0. The zero-dimensional nontrivial steady state is = 1.16, = 231.72. In Fig. 14, 15, and the left side of Fig. 16 we simulate the model. The tumor is initially small and located on the right side of Ω = [0, 1]. Over time, as the tumor grows from the right side of Ω = [0, 1], the regular spatial and temporal cycles of the osteoclast and osteoblast populations are disrupted with these solutions ultimately approaching the zero-dimensional tumor model nontrivial steady state = 5.0, = 316.0, as in Fig. 4 (Fig. 14), the bone mass is depleted throughout Ω (Fig. 15), and the tumor grows to carrying capacity (Fig. 16, left side).

One-dimensional Bone Model with Tumor and Treatment

The one-dimensional model with tumor and treatment is obtained by adding the time-dependent treatment functions V_{1}(t) and V_{2}(t) as in (13). The equations are

(37)

(38)

(39)

with boundary conditions

(40)

(41)

(42)

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 illustrate the model in this section by adding treatment to the simulation in Fig. 14 and Fig. 15. All parameters and initial conditions are as in Fig. 11, 14, and 15. The treatment parameters are v_{1} = .0001 in (13b) and v_{2} = .006 in (13d) and t_{
start
}= 600. This form of the treatment corresponds to drugs such as proteasome inhibitors, which promote osteoblast production and inhibit tumor growth. The graphs of osteoclast and osteoblast populations are given in Fig. 17, where both populations recover to regular cycles after the start of treatment at t_{
start
}= 600. The bone mass is graphed in Fig. 18 and the tumor population is graphed on the right side of Fig. 16. The osteoclast and osteoblast populations recover normal cycling (Fig. 17), the bone mass recovers partially on the left side, but not on the right of Ω (Fig. 18), and the tumor is extinguished by the treatment (Fig. 16, right side).

Discussion

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 [15]. 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.

Conclusions

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.

Methods

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.

Authors Contributions

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.

Reviewers Comments

Comments from Ariosto Silva

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 [35] and pH in the activity of osteoblasts [36]. 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 [37] and Multiple Myeloma [38]. 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 [22] 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.

Comments from Mark P. Little

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. [22]).

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.

Declarations

Acknowledgements

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.

Authors’ Affiliations

(1)

Department of Mathematics, University of Iowa

(2)

Vanderbilt Center for Bone Biology, Department of Cancer Biology, Vanderbilt University Medical Center

(3)

Department of Mathematics, Vanderbilt University

(4)

Departments of Biomedical Engineering, Molecular Physiology and Biophysics, and Physics and Astronomy, Vanderbilt University

(5)

Vanderbilt Institute for Integrative Biosystems Research and Education

References

Parfit AM: Osteonal and hemi-osteonal remodeling: the spatial and temporal framework for signal traffic in adult human bone.J Cell Biochem 1994, 55:273–86.View Article

Martin MJ, Buckland-Wright JC: Sensitivity analysis of a novel mathematical model identifies factors determining bone resorption rates.Bone 2004, 35:918–928.PubMedView Article

Lekszycki T: Functional Adaptation of Bone as an Optimal Control Problem.J Theoretical Applied Mechanics 2005,43(3):555–574.

Maldonado S, Borchers S, Findeisen R, Allgöwer F: Mathematical Modeling and Analysis of Force Induced Bone Growth.Proceedings of the 28th IEEE EMBS Annual International Conference, New York 2006.

Martínez G, Aznar JMG, Doblaré M, Cerrolaza M: External bone remodeling through boundary elements and damage mechanics.Mathematics and Computers in Simulation 2006, 73:183–199.View Article

Tezuka KI, Wada Y, Takahashi A, Kikuchi M: Computer-simulated bone architecture in a simple bone-remodeling model based on a reaction-diffusion system.J Bone Miner Metab 2005, 23:1–7.PubMedView Article

Lemaire V, Tobin FL, Greller LD, Cho CR, Suva LJ: Modeling the interactions between osteoblast and osteoclast activities in bone remodeling.J Theor Biol 2004, 229:293–309.PubMedView Article

García-Aznar JM, Rueberg T, Doblare M: A bone remodelling model coupling microdamage growth and repair by 3D BMU-activity.Biomech Model Mechanobiol 2005, 4:147–167.PubMedView Article

Moroz A, Wimpenny DI: Allosteric control model of bone remodelling containing periodical modes.Biophysical Chemistry 2007, 127:194–212.PubMedView Article

Pivonka P, Zimak J, Smith DW, Gardiner BS, Dunstan CR, Sims NA, Martin TJ, Mundy GR: Model structure and control of bone remodeling: A theoretical study.Bone 2008, 43:249–263.PubMedView Article

Restrepo JM, Choksi R, Jiang JMHY: Improving the damage accumulation in a biomechanical bone remodelling model.Computer Methods in Biomechanics and Biomedical Engineering 2009,12(3):341–352.PubMedView Article

Ryser MD, Nigam N, Komarova SV: Mathematical Modeling of Spatio-Temporal Dynamics of a Single Bone Multicellular Unit.J Bone Mineral Research 2009,24(5):860–870.View Article

Mundy GR: Myeloma bone disease.European J Cancer 1998, 34:246–251.View Article

Mundy GR, Raisz LG, Cooper RA, Schechter GP, Salmon S: Evidence for the secretion of an osteoclast stimulating factor in myeloma.New England J Medicine 1974, 291:1041–1046.View Article

Bataille R, Chappard D, Marcelli C, Dessauw P, Baldet P, Sany J, Alexandret C: Recruitment of New Osteoblasts and Osteoclasts Is the Earliest Critical Event in the Pathogenesis of Human Multiple Myeloma.J Clin Invest 1991, 88:62–66.PubMedView Article

Valentin OA, Charhon SA, Meunier PJ, Edouard CM, Arlot ME: Quantitative histology of myeloma-induced bone changes.British J Haematology 1982, 52:601–610.View Article

Taube T, Beneton MNC, McCloskey EV, Rogers S, Greaves M, Kanis JA: Abnormal bone remodelling in patients with myelomatosis and normal biochemical indices of bone resorption.European J Haematology 1992, 49:192–198.View Article

Evans CE, Galasko CS, Ward C: Does myeloma secrete an osteoblast inhibiting factor?J Bone Joint Surg Br 1989, 71:288–290.PubMed

Bataille R, Delmas PD, Chappard D, Sany J: Abnormal serum bone Gla protein levels in multiple myeloma. Crucial role of bone formation and prognostic implications.Cancer 1990, 66:167–172. PublisherFullTextPubMedView Article

Abildgaard N, Rungby J, Glerup H, Brixen K, Kassem M, Brincker H, Heickendorff L, Eriksen EF, Nielsen JL: Long-term oral pamidronate treatment inhibits osteoclastic bone resorption and bone turnover without affecting osteoblastic function in multiple myeloma.Eur J Haematol 1998, 61:128–134.PubMedView Article

Woitge HW, Horn E, Keck AV, Auler B, Seibel MJ, Pecherstorfer M: Biochemical markers of bone formation in patients with plasma cell dyscrasias and benign osteoporosis.Clin Chem 2001, 47:686–693.PubMed

Komarova SV, Smith RJ, Dixon SJ, Sims SM, Wahl LM: Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone remodeling.Bone 2003, 33:206–215.PubMedView Article

Akchurin T, Aissiou T, Kemeny N, Prosk E, Nigam N, Komarova SV: Complex Dynamics of Osteoclast Formation and Death in Long-Term Cultures.PLoS ONE 2008,3(5):e2104.PubMedView Article

Komarova SV: Mathematical Model of Paracrine Interactions between Osteoclasts and Osteoblasts Predicts Anabolic Action of Parathyroid Hormone on Bone.Endocrinology 2005,146(8):3589–3595.PubMedView Article

Komarova SV: Bone Remodeling in Health and Disease: Lessons From Mathematical Modeling.Ann NY Acad Sci 2006, 1068:557–559.PubMedView Article

Edwards CM, Lwin ST, Fowler JA, Oyajobi BO, Zhuang J, Bates AL, Mundy GR: Myeloma cells exhibit an increase in proteasome activity and an enhanced response to proteasome inhibition in the bone marrow microenvironment in vivo.Am J Hematol 2009,84(5):268–272.PubMedView Article

Garrett IR, Chen D, Gutierrez G, Zhao M, Escobedo A, Rossini G, Harris SE, Gallwitz W, Kim KB, Hu S, Crews CM, Mundy GR: Selective inhibitors of the osteoblast proteasome stimulate bone formation in vivo and in vitro.J Clin Invest 2003,111(11):1771–1782.PubMedView Article

LeBlanc R, Catley LP, Hideshima T, Lentzsch S, Mitsiades CS, Mitsiades N, Neuberg D, Goloubeva O, Pien CS, Adams J, Gupta D, Richardson PG, Munshi NC, Anderson KC: Proteasome Inhibitor PS-341 Inhibits Human Myeloma Cell Growth in Vivo and Prolongs Survival in a Murine Model.Cancer Research 2002, 62:4996–5000.PubMed

Oyajobi BO, Garrett IR, Gupta A, Flores A, Esparza J, Munoz S, Zhao M, Mundy GR: Stimulation of new bone formation by the proteasome inhibitor, bortezomib: implications for myeloma bone disease.British J Haematology 2007,139(3):434–438.View Article

Pennisi A, Li X, Ling W, Khan S, Zangari M, Yaccoby S: The proteasome inhibitor, bortezomib suppresses primary myeloma and stimulates bone formation in myelomatous and nonmyelomatous bones in vivo.Am J Hematol 2008, 84:6–14.View Article

Uy GL, Trivedi R, Peles S, Fisher NM, Zhang QJ, Tomasson MH, Dipersio JF, Vij R: Bortezomib inhibits osteoclast activity in patients with multiple myeloma.Clin Lymphoma Myeloma 2007, 7:587–589.PubMedView Article

Hideshima T, Richardson P, Chauhan D, Palombella VJ, Elliott PJ, Adams J, Anderson KC: The Proteasome Inhibitor PS-341 Inhibits Growth, Induces Apoptosis, and Overcomes Drug Resistance in Human Multiple Myeloma Cells.Cancer Research 2001, 61:3071–3076.PubMed

Terpos E, Heath DJ, Rahemtulla A, Zervas K, Chantry A, Anagnostopoulos A, Pouli A, Katodritou E, Verrou E, Vervessou EC, Dimopoulos , Meletios-Athanassios , Croucher PI: Bortezomib reduces serum dickkopf-1 and receptor activator of nuclear factor-[kappa]B ligand concentrations and normalises indices of bone remodelling in patients with relapsed multiple myeloma.British J Haematology 2006,135(5):688–692.View Article

Zavrski I, Krebbel H, Wildemann B, Heider U, Kaiser M, Possinger K, Sezer O: Proteasome inhibitors abrogate osteoclast differentiation and osteoclast function.Biochem Biophys Res Commun 2005, 333:200–205.PubMedView Article

Bozec A, Bakiri L, Hoebertz A, Eferl R, Schilling A, Komnenovic V, Scheuch H, Priemel M, Stewart C, Amling M, Wagner E: Osteoclast size is controlled by Fra-2 through LIF/LIF-receptor signalling and hypoxia.Nature 2008, 454:221–225.PubMedView Article

Kohn D, Sarmadi M, Helman J, Krebsbach P: Effects of pH on human bone marrow stromal cells in vitro: implications for tissue engineering of bone.J Biomed Mater Res 2002, 60:292–299.PubMedView Article

Fortney J, Zhao W, Wenger S, Gibson L: Bone marrow stromal cells regulate caspase 3 activity in leukemic cells during chemotherapy.Leuk Res 2001, 25:901–907.PubMedView Article

Kirshner J, Thulien K, Martin L, Debes Marun C, Reiman T, Belch A, Pilarski L: A unique three-dimensional model for evaluating the impact of therapy on multiple myeloma.Blood 2008, 112:2935–2945.PubMedView Article

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.