Towards a New Spatial Representation of Bone Remodeling

Irregular bone remodeling is associated with a number of bone diseases such as osteoporosis and multiple myeloma. Computational and mathematical modeling can aid in therapy and treatment as well as understanding fundamental biology. Different approaches to modeling give insight into different aspects of a phenomena so it is useful to have an arsenal of various computational and mathematical models. Here we develop a mathematical representation of bone remodeling that can effectively describe many aspects of the complicated geometries and spatial behavior observed. There is a sharp interface between bone and marrow regions. Also the surface of bone moves in and out, i.e. in the normal direction, due to remodeling. Based on these observations we employ the use of a level-set function to represent the spatial behavior of remodeling. We elaborate on a temporal model for osteoclast and osteoblast population dynamics to determine the change in bone mass which influences how the interface between bone and marrow changes. We exhibit simulations based on our computational model that show the motion of the interface between bone and marrow as a consequence of bone remodeling. The simulations show that it is possible to capture spatial behavior of bone remodeling in complicated geometries as they occur \emph{in vitro} and \emph{in vivo}. By employing the level set approach it is possible to develop computational and mathematical representations of the spatial behavior of bone remodeling. By including in this formalism further details, such as more complex cytokine interactions and accurate parameter values, it is possible to obtain simulations of phenomena related to bone remodeling with spatial behavior much as \emph{in vitro} and \emph{in vivo}. This makes it possible to perform \emph{in silica} experiments more closely resembling experimental observations.


Introduction
Bones of the mature skeleton undergo constant remodeling through resorption and replacement of old matrix.The repair of micro-cracks caused by normal mechanical stresses entails stimulation of local remodeling by tissue-intrinsic factors [3,19,25].Bone remodeling also plays a role in systemic Ca ++ homoeostasis and is globally regulated by humoral factors (e.g.parathyroid hormone).Because the major mechanisms regulating bone metabolism are well-established and the links between abnormal metabolism and bone diseases are well-documented, the remodeling process is a potentially fertile ground for mathematical and computational analyses.Computer models that simulate the complexly regulated biology of bone remodeling can be used to rapidly and inexpensively screen drugs to estimate their disease-modifying activity.By minimizing the need for much more lengthy and expensive in vivo testing, such models can greatly accelerate the development of novel treatments for bone diseases.
Bone remodeling occurs simultaneously at various locations.However, remodeling from one site to another is not necessarily synchronous.Remodeling is localized to a particular remodeling site with teams of cells acting together as an individual structure commonly referred to as a basic multicellular unit (BMU) [3,19,24,25].Remodeling is driven locally by signals from osteocytes, which secrete chemotactic factors in response to mechanically induced microfractures and in response to bone strain [3,19,25].With the release of factors by osteocytes stromal stem cells are activated and begin to produce factors, such as macrophage colonystimulating factor (M-CSF), that promote osteoclastogenesis [3,19].Additional stromal cells are created and differentiate into pre-osteoblasts.These cells begin producing factors such as receptor activator of nuclear factor-B ligand (RANKL), which binds to receptors present in pre-osteoclasts and inhibits apoptosis [3,19].This is a case of osteoblast-derived paracrine signaling and is featured in many of the mathematical treatments of remodeling discussed below.In the presence of the factors M-CSF and RANKL pre-osteoclasts fuse to become mature osteoclasts.Osteoclasts bind to bone surfaces and resorb the matrix by secreting acids and matrix proteases.Mature osteoclasts in turn produce growth factors and other signaling molecules, an example of osteoclast-derived autocrine signaling.Resorption also releases growth factors such as IGF and TGF-β that act on cells of the osteoblast lineage, an example of osteoclast-derived paracrine signaling [3,19].
When resorption is nearly complete pre-osteoblasts mature into osteoblasts and stop secreting RANKL.Instead they produce OPG, a decoy receptor that blocks RANKL binding to osteoclasts.As a result, osteoclasts cease resorptive activity and undergo apoptosis.During formation some osteoblasts become trapped in the matrix as osteocytes, while some undergo apoptosis [3,19].This completes the description of a typical scenario for remodeling by a BMU in response to mechanical stress.This description is helpful in the development of the mathematical and computational models discussed below.
A first step towards formulating a theoretical representation of bone remodeling is to incorporate the mechanisms that regulate the process.Here this is the autocrine and paracrine biochemical signaling amongst cells.Komarova et al [11,12] develop a temporal model consisting of ordinary differential equations (ODEs) describing the local dynamics of the remodeling process.This model uses a power law approximation to capture the interaction of cell populations through autocrine and paracrine signaling.An alternative approach to that in [11,12] also capturing the cytokine signaling regulating bone remodeling is given by Pivonka et al in [23].The model presented in [23] incorporates signaling through Hill functions rather than power laws.Both approaches model the changes in osteoclast and osteoblast populations in response to RANKL/OPG signaling.Ayati et al [2] and Ryser et al [26,27] build on the model from [12] for the local dynamics to incorporate some of the spatial effects of bone remodeling.In Ayati et al [2], in the context of multiple myeloma, diffusion of cells is added to incorporate some spatial effects of bone remodeling, while in [26] Ryser et al consider chemotaxis of osteoclasts in the presence of remodeling promoting factors.
In this article we employ a geometric method developed by Osher and Sethian [17] to present another approach to capturing spatial effects of bone remodeling.We represent, via the so called level-set equation, the evolution of the interface between bone and marrow regions in response to remodeling.The interface is described geometrically as a closed curve in the plane.The interior of the curve models region of bone while the exterior to the curve represents the marrow region.A level-set function, which is determined by solving the level-set equation, provides a convenient way to represent a geometric object such as a curve.Through this geometric approach it is possible to consider effects of remodeling for complex geometries and spatial scales not easily represented with previous spatial models.This is particularly relevant for representing the remodeling of cancellous or trabecular bone tissue.By interpreting the change in total bone mass, at each point in space, during remodeling as the speed of the interface motion (normal, i.e. perpendicular, to the interface) 1 we obtain a novel representation of the spatial behavior of bone remodeling.We note that while in many other applications of this geometric approach, the motion of the interface is determined by some aspect of the local geometry such as curvature, this is not the case for the application we give here.The motion of the interface is a result of the behavior of the cell populations at each point in space, in other words, the changes in geometry are directly determined by biological phenomena alone.We find that our modeling and simulation efforts are representative of observations in animal models.Incorporating pathological behavior associated with abnormal bone remodeling in future work is straightforward and of clinical interest.

Models
The complex nature of the interactions of osteoclasts and osteoblasts and the chemical signaling that regulates them makes a detailed theoretical and experimental study difficult.Mathematical modeling may be used to tie existing knowledge about a system with new ideas.Simulations based on mathematical models can provide a useful description of observed or hypothesized results otherwise difficult with ordinary language.
In this section we present and discuss some mathematical models of bone remodeling with a focus on the ones that have influenced the development of our approach.The models presented here describe the change in osteoclast and osteoblast populations throughout the remodeling process.Throughout we denote by C the osteoclast population and B the osteoblast population.The model which forms the foundation for all that follows is the temporal model of Komarova et al [12] that describes the local dynamics of bone remodeling at a single remodeling site.The dynamics of the cell populations is determined by the interaction of cells through autocrine and paracrine signaling.In [12] the authors employ a power law approximation to describe the overall effects of the biochemical factors involved in signaling.This type of modeling, while simple to describe, is effective in capturing the nonlinear behavior of biochemical signaling.The exponents in the power law approximation describe the effectiveness of the autocrine and paracrine signaling on each cell type.From [12] the dynamics of the cell populations are given by Here α 1 and α 2 are constants that reflect the rate of cell production while while β 1 and β 2 are constants reflecting the rate of cell removal of osteoclasts and osteoblasts respectively.As mentioned above the exponents g ij describe the effectiveness of the different types of signaling on each cell type with g 11 describing the effectiveness of osteoclast derived autocrine signaling, g 21 the effectiveness of osteoblast derived paracrine signaling, g 12 the effectiveness of osteoclast derived paracrine signaling and g 22 the effectiveness of osteoblast derived autocrine signaling.The authors use this model to simulate the behavior of remodeling and find that the value of the parameter g 11 for the effectiveness of osteoclast derived autocrine signaling most strongly effects the mode of remodeling behavior.The authors also calculate the change in total bone mass due to remodeling via the equation Here z is the bone mass k 1 , k 2 are constants describing the activity of bone resorption and formation respectively.The constants C and B denote the stable steady state values of osteoclasts and osteoblasts respectively and here are given by In the level-set representation that follows we take the change in bone mass as giving the speed of motion of the interface in the normal direction.We note that in the work of Pivonka et al [23] there is an equation for the change in bone volume similar to (1c).Thus the equations from [23] can be used in the following in place of (1a) and (1b).However, for computational simplicity we use the power law approximation since the resulting interface dynamics should be the same, for the same z field, regardless which model is used for the cell populations is used to obtain such a z field.For example, in our case we can obtain the same z field from a system with and without explicit chemical concentrations.For convenience we reproduce (see figure 1) the change in bone mass computed in [12] as a percentage of initial bone mass using the system (1a),(1b), and (1c).Figure 1 can then be used for qualitatively comparison with the results in this paper.
In [2] the authors give an implicit spatial representation of bone remodeling by interpreting the bone mass equation as reflecting bone thickness.We use this interpretation with the level set method discussed below to give a first study of the dynamic behavior of the bone marrow interface due to remodeling.This assumes an unrealistic degree of homogeneity but is an instructive first step in developing the level set approach.The authors in [2] go on to incorporate explicitly spatial effects through the addition of linear diffusion terms to the system (1).This is all done in the context of multiple myeloma bone disease since the goal of their work was to show how perturbation of the mathematical models for bone remodeling can exhibit pathological behavior as is seen in association with metabolic bone diseases.The modeling approach developed below can also be used to study how pathological behavior can arise from perturbations representing irregular remodeling.
In [26] the authors develop a spatial representation of bone remodeling by modifying the system (1) to include motility of osteoclasts and explicit representation of the RANKL and OPG fields.The model is particularly successful in capturing the motion and steering of the BMU and simulations in one and two dimensions show the distinct features of the cutting cone of osteoclast cells.Their results are obtained by including the chemotaxis of osteoclasts towards increasing concentrations of RANKL and explicitly the dynamics of the RANKL/OPG pathway.The equations for RANKL and OPG concentrations include porous flow of the chemicals, activation of osteclastogenesis by RANKL and binding of RANKL/OPG.The two dimensional simulations based on the model elegantly show steering, as a chemotactic response, of a remodeling BMU along a microfracture .Also of note are the mathematical models of Martin and Buckland-Wright [13,14].The model in [13] simulates erosion depth due to resorption by modeling the changes in cellular activity due to cytokine signaling.This is done using the equations for enzyme kinetics of Michaelis and Menten.The model in [13] differs from the mathematical models in [2,11,12,23,26] in that it does not describe changes cell populations while the model in [14] does include changes in the osteoblast population.The models found in [13,14] provide another useful method for studying how changes in signaling are manifested as physical changes in the bone.The recent review paper of Pivonka and Komarova [22] gives a nice overview mathematical modeling in bone biology.
In the models discussed so far bone remodeling is regulated by biochemical signaling alone.Mechanical signaling also plays an important role in the regulation of cells involved in bone biology.There is a substantial body of work devoted to combining the biochemical and mechanical stimuli in mathematical models to study the effects that the two different types of stimuli have together [4][5][6][7][8]20].
For example in [8] the authors couple mechanical stimuli through finite element analysis with partial differential equation models describing spatio-temporal changes in cell populations.These changes are due to cell migration (chemotaxis, haptotaxis), proliferation and differentiation.They also include equations describing angiogenesis.Biochemical signaling is incorporated by including interactions of cell populations with growth factors.This all leads to a mathematical framework that is successful in simulating the combined effects of biochemical and mechanical regulatory mechanisms.An interesting result is the prediction of overload induced nonunion formation.The modeling approach just discussed is well suited to the study of wound and fracture healing and issues such as angiogenesis that arise under such situations.This is the case for example with larger fractures in or on cortical bone, rather than microfractures.Here we leave the role of angiogenesis to future considerations and focus on the spatial manifestations of other features of the remodeling process.
The geometric approach mentioned, was initially developed by Osher and Sethian in [17] to study geometry driven motion of an interface and is particularly well suited for problems in which there is a definite "inside" and "outside" such as in the modeling of soap bubbles.The method has since been further developed and extended, see for example [16,18,29,31,32], and has found many applications in various areas of applied mathematics and engineering [16,29,32].Since the surface of bone forms a distinct, sharp interface the levelset function seems well suited for representing physical changes to the surface of bone due to remodeling.
We present here some basic background on the level-set function approach specific to representing the spatial effects of bone remodeling in trabecular bone tissue.For a more in depth treatment and other applications see the basic texts by Osher and Fedkiw [16] and Sethian [29].As stated previously, we would like to describe the interface between bone and marrow regions as a closed2 curve in the plane.The basic idea is to represent an interface, geometrically a curve, surface, etc. implicitly as some level set, typically the zero level set, of a higher dimensional object described by a so called level-set function.For example consider the function φ(x, y) := 1 − x 2 − y 2 .The zero level set of φ(x, y) is the set of all points (x, y) in the plane satisfying φ(x, y) which describes a circle of radius 1 centered at the point (0, 0).Alternatively one could represent an interface locally by explicit functions, but this may require many different such function, or parametrically.However, the approach to representing geometry via a level-set function provides some advantages over these alternatives that make it well suited for applications.One advantage to representing an interface, in this case a curve, implicitly as above is that most of the geometric information one would like to have about the interface can be obtained from the level set function.For more details on this see the first chapter of the book by Osher and Fedkiw [16].What really makes the level-set function appealing for applications such as the changes to bone surface due to remodeling, is that, given an initial interface represented by a specified level set function, there is a simple way to evolve the interface, by evolving the level set function in time.More precisely, to represent the evolution of an interface we use a family of level-set functions φ(x, t) parameterized by time t so that at each time t the zero level set {x : φ(x, t) = 0} represents the interface at time t.Here x ∈ R n for some positive integer n, the case we are concerned with has n = 2. Then following Osher and Sethian in [17] we can derive a differential equation for the motion of the interface in the normal direction due to normal velocity with speed given by a function a(x, t).The equation is and is known as the level-set equation.Here φ t denotes the derivative with respect to time, ∇φ(x, t) is the gradient of φ in the spatial variables and |ȳ| is the length of a vector ȳ ∈ R n .In a sense the function a(x, t) gives the rule telling what drives the motion of the interface and solving the level set equation describes what the interface looks like at any given time.Notice that here the level set equation does not contain a curvature term so that local geometry is not used to move the interface.The ODEs (1a) and (1b) effectively describe the local population dynamics of the primary cells involved in remodeling.Alternatively, these ODEs can represent a larger, but well-mixed spatial domain.We have discussed previous approaches to incorporating spatial heterogeneity into mathematical descriptions of bone remodeling.Even with the success of these models it is valuable to move towards new and different spatial representations of bone remodeling.In that direction we employ the level-set function approach to develop such a representation.The interface between bone and marrow is distinct and sharp and we use this to motivate representing the system as one with a moving interface.We use the level-set equation discussed previously to describe the motion of this interface.
Since the power law approximation of [12] successfully captures the local dynamics of osteclast/osteoblast cell populations we elaborate (1a) (1b) by adding explicit space dependence and coupling with the level-set equation to represent spatial manifestations of the remodeling process.In other words, we (locally) couple the biological aspects, i.e. population dynamics of remodeling cells, of bone bone remodeling with the change in geometry of the interface between bone and marrow regions.This is done by assuming that the speed in the direction normal to the interface, a(x, t) in ( 4), is given by the change in bone mass.More precisely, the speed in the normal direction is given by a constant multiple of the change in bone mass, where ν is a constant.We note that here there is a (potentially) different population of each cell type at each location x in space.This makes the speed a(x, t) different at each point in space.In [26] chemotaxis plays the role of a steering mechanism, determining the spatial behavior of the osteoclast and osteoblast populations along an existing microfracture on the order of a few microns.In our model the spatial behavior, including steering of cells, is determined by the geometry via the level set method.This is appropriate for spatial scales on the order of a thousand microns or larger and is well suited for geometrical configurations common in trabecular bone.If necessary, slight modification of our approach allows for the inclusion of chemotaxis or other steering mechanisms.

Results and Discussion
In what follows, we denote by Γ t = {x : φ(x, t) = 0} the interface at time t, where φ is the level-set function.
We solve the following system of equations to simulate the motion of the interface between bone and marrow due to bone remodeling, where a(x, t) is given by ( 5), and C(x, t) and B(x, t) are the number of osteoclasts and osteoblasts respectively, recruited from the marrow to a point x on interface at time t.Thus we have a system of equations where C(x, t), B(x, t) depend on the interface, i.e. level set function φ(x, t), since at each time, the population at the interface Γ t is updated according to (6d) and (6e).The rapid motion of the cells from deeper within the marrow to a remodeling site is captured by the local recruitment terms (6d) and (6e).When remodeling is initiated the cells are recruited quickly and directly to the remodeling site.As a result, diffusion of these cells is not relevant for the purposes we are considering.The domain for this problem is a square.To mimic a large domain, for computational simplicity we assume periodic boundary conditions for this domain.The initial interface is described by the zero level set of an initial level-set function which we denote by φ 0 (x, y).
For this discussion we take the initial interface to be a circle but note that any other shape could be taken by prescribing an appropriate level-set function.It is even possible to use, as a level set function, images taken from experiments, see for example [30].For the numerical simulations presented below we use parameter values as in [2,12,27], these are C = 1.06 B = 212.13 All of our simulations have been carried out using MATLAB.The implementation of our model required solving a coupled system of differential equations at each point in the domain.To carry this out we first solved the equations for the cell populations using implicit time stepping and a Krylov iterative method [9,10] to solve the nonlinear system.This provided the normal speed for the motion of the interface.Evolving the interface required solving the level-set equation.For this we used the level set toolbox for MATLAB [15] developed by Ian Mitchell.This toolbox contains many useful tools for implementing the level-set numerical method in MATLAB.
We begin remodeling at three distinct sites by perturbing the osteoclast population away from steady state, figure 2(a).This initiates bone resorption.The effects of remodeling on the interface after 5 days of bone resorption is shown in figure 4(b).Note that the interface moves inward (toward the interior region on bone) due resorption since there is an increase in osteoclast population away from steady state.The osteoclast population quickly returns to steady state, figure 2(f).In accordance with equation 1c when the osteoclast population returns to steady state the motion should cease to be inward.This can be seen in figures 4(e) and 4(f).In the meantime the osteoblast population increases from steady state, figures 3(b), 3(c), and new bone is being formed.The effect of remodeling on the interface during the early stages of formation of new bone by osteoblasts is shown if figures 4(f) and 4(g) where the indentions become filled in.After reaching approximately maximum population density, figures 3(e), 3(f), the osteoblasts return to steady state by which time bone has been fully remodeled, figures 3(i), 4(i).
Comparing the results in the three figures 2,3 and 4 at corresponding times shows how the changes in the cell populations are manifested as physical changes in the bone surfaces at three different remodeling sites.The interface moves in the normal direction as a result of remodeling in accord with (1).While here we have initiated remodeling at each of the three sites at the same instant this is not a requirement of this modeling approach.Remodeling at different sites can easily be taken to be out of phase.
We have described in figure 5 quantitatively the change in the bone by plotting the area (in square microns) interior to bone, this can be compared with figure 2, which shows the bone mass as a function of time.We remark that bone mass resorption and formation occurs somewhat faster in figure 5 as we add more spatial dimensions due to the geometry, (i.e.proportional to radius squared versus proportional to local thickness).We note however, the qualitative agreement with figure 2 in [12].
As discussed previously, in the work [26], a mathematical model including explicit spatial terms is considered.To illustrate the flexibility of the level-set approach we discuss how this can be done here as well.The first step is the inclusion of signaling pathways not through exponents in a power law approximation but explicitly.For example the RANKL/OPG pathway can be represented explicitly by modifying (1a) to  include the effect of RANKL on osteoclast arriving at Then the RANKL/OPG field is represented by equations for the concentrations of each chemical species.The population dynamics of osteoblast cells at three separate remodeling site.The osteoblast population starts at steady state, then increases in response to an increase in osteoclast population, and finally returns to steady in accord with (1a) and (1b).Bone is formed in the presence of an osteoclast population above steady state levels.These are The evolution of the interface between bone and marrow regions due to remodeling at three distinct remodeling cites.The speed in the normal direction is given by the change in the bone mass (5). Figure 4(a) shows the initial interface before any remodeling has begun.The interior of the circle represents the bone region while the exterior is marrow.The interface between bone and marrow moves in and out, i.e. in the normal direction, due to changes in the populations of osteoclasts and osteoblasts.
where u denotes the concentration of RANKL and v the concentration of OPG.The equations for the osteoblast population and change in bone mass remain as in (1b) and ( 5) respectively.Note that while in [26] steering is controlled by the chemotaxis of osteoclasts toward RANKL concentrations, here it is controlled by the geometry through the level-set function as above.Our computational results including explicit spatial terms for the chemical species showed an evolution of the interface similar to figure 4. The amount of change due to remodeling.The units are days for the horizontal axis and mm 2 for the vertical axis.We arrived at this by computing the interior area at each time in the evolution of the interface.This can be compared with figure 1.The differences in geometry (i.e.velocity times radius squared here, and local bone mass in 1) results in a slightly faster, but qualitatively similar, cycle than in figure 1.
Based on the local dynamics of bone remodeling as developed in [11,12,26,27] we have developed a new approach to modeling the spatial effects of bone remodeling.Through the use of the level-set equation we link local dynamics with geometric effects and spatial behavior such as steering of a BMU in a way that is appropriate on many different spatial scales.The approach taken is useful but not limited to describing the effects of bone remodeling in complex geometries such as is seen in the remodeling of trabecular bone.The extension of this model to various other situations is straightforward and software for implementing the level-set equation approach numerically is widely available (see e.g.[15] and references therein).
The development of mathematical and computational models of biological phenomena is valuable in design and analysis of experiments and aids in building intuition and conceptual understanding of biological processes.These models can also be useful to the biomedical community.Analysis and simulation based on mathematical and computational models can suggest how and where pathological behavior occurs.Examples of this in bone remodeling can be seen in [2,11,12] and the modeling approach in this paper can be used in the same way by including the appropriate modifications.Modeling can also aid in considering treatments and therapies for osteopathies (see e.g.[1]).We feel the type of modeling presented in this work can be especially helpful in ruling out ineffective treatments and therapies.
For future work more sophisticated numerical schemes for solving the level-set equation and computing solutions of the model equations will be employed.In particular since remodeling occurs primarily near a small region of the interface, a local level-set method such as in [21,28] can be used to speed up computations.It may also be appropriate in the future to include further details of the cytokine signaling, something that is likely to involve coupling additional differential equations with the existing model.This could call for a need of more efficient and accurate numerical algorithms to carry out simulations, particularly important in moving to three dimensional problems.While the level-set equation scales well in higher dimensions, numerical solutions of other differential equations become much more expensive from a computational standpoint.We feel the level-set approach taken here will provide a supplement or alternative to other models of the fundamental process of bone remodeling.Further development in this direction will allow for applications to a variety of biomedical issued related to bone remodeling at a variety of scales.

Figure 1 :
Figure1: Change in bone mass according to the model from[12].This figure was obtained by solving equations (1a),(1b), and (1c).Here the horizontal axis is time in days and the vertical is percentage of initial (100%) bone mass.

Figure 2 :
Figure 2:The population dynamics of osteoclast cells at three separate remodeling sites.Remodeling is initiated at three distinct sites by perturbing the osteoclast population away from steady state.The osteoclast population quickly returns to steady state.Bone is resorbed in the presence of an osteoclast population above steady state levels.

Figure 3 :
Figure 3:The population dynamics of osteoblast cells at three separate remodeling site.The osteoblast population starts at steady state, then increases in response to an increase in osteoclast population, and finally returns to steady in accord with (1a) and (1b).Bone is formed in the presence of an osteoclast population above steady state levels.

Figure 4 :
Figure 4:The evolution of the interface between bone and marrow regions due to remodeling at three distinct remodeling cites.The speed in the normal direction is given by the change in the bone mass(5).Figure4(a)shows the initial interface before any remodeling has begun.The interior of the circle represents the bone region while the exterior is marrow.The interface between bone and marrow moves in and out, i.e. in the normal direction, due to changes in the populations of osteoclasts and osteoblasts.

Figure 5 :
Figure 5:The amount of change due to remodeling.The units are days for the horizontal axis and mm 2 for the vertical axis.We arrived at this by computing the interior area at each time in the evolution of the interface.This can be compared with figure1.The differences in geometry (i.e.velocity times radius squared here, and local bone mass in 1) results in a slightly faster, but qualitatively similar, cycle than in figure1.