







The mechanical properties of vertebrate bone are largely determined by a process which involves the complex interplay of three different cell types. This process is called bone remodeling and occurs asynchronously at multiple sites in the mature skeleton. The cells involved are bone resorbing osteoclasts, bone matrix producing osteoblasts, and mechanosensing osteocytes. These cells communicate with each other by means of autocrine and paracrine signaling factors and operate in complex entities, the so-called bone multicellular units (BMUs). To investigate the BMU dynamics in silico, we develop a novel mathematical model resulting in a system of nonlinear partial differential equations (PDEs) with time delays. The model describes the osteoblast and osteoclast populations together with the dynamics of the key messenger molecule RANKL and its decoy receptor OPG. Scaling theory is used to address parameter sensitivity and predict the emergence of pathological remodeling regimes. The model is studied numerically in one and two space dimensions using finite difference schemes in space and explicit delay equation solvers in time. The computational results are in agreement with in vivo observations and provide new insights into the role of the RANKL/OPG pathway in the spatial regulation of bone remodeling.
Mathematical Modeling of Spatio-Temporal Dynamics of a Single Bone Multicellular Unit
Abstract During bone remodeling, bone-resorbing osteoclasts and bone-forming osteoblasts are organized in bone multicellular units (BMUs), which travel at a rate of 20–40 μm/d for 6–12 mo, maintaining a cylindrical structure. However, the interplay of local BMU geometry with biochemical regulation is poorly understood. We developed a mathematical model of BMU describing changes in time and space of the concentrations of proresorptive cytokine RANKL and its inhibitor osteoprotegerin (OPG), in osteoclast and osteoblast numbers, and in bone mass. We assumed that osteocytes surrounding a microfracture produce RANKL, which attracted osteoclasts. OPG and RANKL were produced by osteoblasts and diffused through bone, RANKL was eliminated by binding to OPG and RANK. Osteoblasts were coupled to osteoclasts through paracrine factors. The evolution of the BMU arising from this model was studied using numerical simulations. Our model recapitulated the spatio-temporal dynamics observed in vivo in a cross-section of bone. In response to a RANKL field, osteoclasts moved as a well-confined cutting cone. The coupling of osteoclasts to osteoblasts allowed for sufficient recruitment of osteoblasts to the resorbed surfaces. The RANKL field was the highest at the microfracture in front of the BMU, whereas the OPG field peaked at the back of the BMU, resulting in the formation of a RANKL/OPG gradient, which strongly affected the rate of BMU progression and its size. Thus, the spatial organization of a BMU provides important constraints on the roles of RANKL and OPG as well as possibly other regulators in determining the outcome of remodeling in the BMU.

<i>In silico</i> models of bone remodeling from macro to nano—from organ to cell
Abstract Computational modeling is a tool through which researchers can achieve a greater understanding of the mechanisms governing biological systems. In the field of bone biology, a plethora of models exist which attempt to replicate and investigate bone's dynamic behavior at different scales. At organ level, models are continuum based and describe the variation of bone's apparent density as a function of both biological and external mechanical stimuli. At tissue level, models include bone microarchitecture and more descriptive parameters such as trabecular thickness, osteoclast resorption depth, and activation frequency. Finally, at cell level, models employ partial differential equations to describe complex cellular interactions in the temporal domain. Although informative, these models exist in isolation. Consequently, their interpretation is limited. In this review, we present an overview of the organ‐, tissue‐, and cell‐level models and assess their ability to reflect bone's metabolic processes reliably. Existing interscale synergies are then presented along with a computational framework which could be exploited to achieve a fully integrated, multiscale modeling approach. WIREs Syst Biol Med 2011 3 241–251 DOI: 10.1002/wsbm.115 This article is categorized under: Analytical and Computational Methods > Computational Methods Models of Systems Properties and Processes > Organ, Tissue, and Physiological Models

A review of mathematical modeling of bone remodeling from a systems biology perspective
Bone remodeling is an essential physiological process in the adult skeleton. Due to the complex nature of this process, many mathematical models of bone remodeling have been developed. Each of these models has unique features, but they have underlying patterns. In this review, the authors highlight the important aspects frequently found in mathematical models for bone remodeling and discuss how and why these aspects are included when considering the physiology of the bone basic multicellular unit, which is the term used for the collection of cells responsible for bone remodeling. The review also emphasizes the view of bone remodeling from a systems biology perspective. Understanding the systemic mechanisms involved in remodeling will help provide information on bone pathology associated with aging, endocrine disorders, cancers, and inflammatory conditions and enhance systems pharmacology. Furthermore, some features of the bone remodeling cycle and interactions with other organ systems that have not yet been modeled mathematically are discussed as promising future directions in the field.
A review of mathematical modeling of bone remodeling from a systems biology perspective
Bone remodeling is an essential, delicately balanced physiological process of coordinated activity of bone cells that remove and deposit new bone tissue in the adult skeleton. Due to the complex nature of this process, many mathematical models of bone remodeling have been developed. Each of these models has unique features, but they have underlying patterns. In this review, the authors highlight the important aspects frequently found in mathematical models for bone remodeling and discuss how and why these aspects are included when considering the physiology of the bone basic multicellular unit, which is the term used for the collection of cells responsible for bone remodeling. The review also emphasizes the view of bone remodeling from a systems biology perspective. Understanding the systemic mechanisms involved in remodeling will help provide information on bone pathology associated with aging, endocrine disorders, cancers, and inflammatory conditions and enhance systems pharmacology. Furthermore, some features of the bone remodeling cycle and interactions with other organ systems that have not yet been modeled mathematically are discussed as promising future directions in the field.

Mathematical modeling of postmenopausal osteoporosis and its treatment by the anti‐catabolic drug denosumab
SUMMARY Denosumab, a fully human monoclonal antibody, has been approved for the treatment of postmenopausal osteoporosis. The therapeutic effect of denosumab rests on its ability to inhibit osteoclast differentiation. Here, we present a computational approach on the basis of coupling a pharmacokinetics model of denosumab with a pharmacodynamics model for quantifying the effect of denosumab on bone remodeling. The pharmacodynamics model comprises an integrated systems biology‐continuum micromechanics approach, including a bone cell population model, considering the governing biochemical factors of bone remodeling (including the action of denosumab), and a multiscale micromechanics‐based bone mechanics model, for implementing the mechanobiology of bone remodeling in our model. Numerical studies of postmenopausal osteoporosis show that denosumab suppresses osteoclast differentiation, thus strongly curtailing bone resorption. Simulation results also suggest that denosumab may trigger a short‐term bone volume gain, which is, however, followed by constant or decreasing bone volume. This evolution is accompanied by a dramatic decrease of the bone turnover rate by more than one order of magnitude. The latter proposes dominant occurrence of secondary mineralization (which is not anymore impeded through cellular activity), leading to higher mineral concentration per bone volume. This explains the overall higher bone mineral density observed in denosumab‐related clinical studies. Copyright © 2013 John Wiley & Sons, Ltd.

The Role of Osteocytes in Targeted Bone Remodeling: A Mathematical Model
Until recently many studies of bone remodeling at the cellular level have focused on the behavior of mature osteoblasts and osteoclasts, and their respective precursor cells, with the role of osteocytes and bone lining cells left largely unexplored. This is particularly true with respect to the mathematical modeling of bone remodeling. However, there is increasing evidence that osteocytes play important roles in the cycle of targeted bone remodeling, in serving as a significant source of RANKL to support osteoclastogenesis, and in secreting the bone formation inhibitor sclerostin. Moreover, there is also increasing interest in sclerostin, an osteocyte-secreted bone formation inhibitor, and its role in regulating local response to changes in the bone microenvironment. Here we develop a cell population model of bone remodeling that includes the role of osteocytes, sclerostin, and allows for the possibility of RANKL expression by osteocyte cell populations. We have aimed to give a simple, yet still tractable, model that remains faithful to the underlying system based on the known literature. This model extends and complements many of the existing mathematical models for bone remodeling, but can be used to explore aspects of the process of bone remodeling that were previously beyond the scope of prior modeling work. Through numerical simulations we demonstrate that our model can be used to explore theoretically many of the qualitative features of the role of osteocytes in bone biology as presented in recent literature.
Spatio-temporal simulations of bone remodelling using a bone cell population model based on cell availability
Here we developed a spatio-temporal bone remodeling model to simulate the action of Basic Multicelluar Units (BMUs). This model is based on two major extensions of a temporal-only bone cell population model (BCPM). First, the differentiation into mature resorbing osteoclasts and mature forming osteoblasts from their respective precursor cells was modelled as an intermittent process based on precursor cells availability. Second, the interaction between neighbouring BMUs was considered based on a ``metabolic cost'' argument which warrants that no new BMU will be activated in the neighbourhood of an existing BMU. With the proposed model we have simulated the phases of the remodelling process obtaining average periods similar to those found in the literature: resorption ($\sim22$ days) - reversal ($\sim$8 days) - formation ($\sim$ 65 days) - quiescence (560 - 600 days) and an average BMU activation frequency of $\sim$1.6 BMUs/year/mm$^3$. We further show here that the resorption and formation phases of the BMU become coordinated only by the presence of TGF-$\upbeta$ (transforming growth factor $\upbeta$, i.e. a major coupling factor stored in the bone matrix. TGF-$\upbeta$ is released through resorption so upregulating osteoclast apoptosis and accumulation of osteoblast precursors, i.e. facilitating the transition from the resorption to the formation phase at a given remodelling site. Finally, we demonstrate that this model can explain targeted bone remodelling as the BMUs are steered towards damaged bone areas in order to commence bone matrix repair.

Mechanobiological osteocyte feedback drives mechanostat regulation of bone in a multiscale computational model
Significant progress has been made to identify the cells and signaling molecules involved in the mechanobiological regulation of bone remodeling. It is now well accepted that osteocytes act as mechanosensory cells in bone expressing several signaling molecules such as nitric oxide (NO) and sclerostin (Scl) which are able to control bone remodeling responses. In this paper, we present a comprehensive multiscale computational model of bone remodeling which incorporates biochemical osteocyte feedback. The mechanostat theory is quantitatively incorporated into the model using mechanical feedback to control expression levels of NO and Scl. The catabolic signaling pathway RANK–RANKL–OPG is co-regulated via (continuous) PTH and NO, while the anabolic Wnt signaling pathway is described via competitive binding reactions between Wnt, Scl and the Wnt receptors LRP5/6. Using this novel model of bone remodeling, we investigate the effects of changes in the mechanical loading and hormonal environment on bone balance. Our numerical simulations show that we can calibrate the mechanostat anabolic and catabolic regulatory mechanisms so that they are mutually exclusive. This is consistent with previous models that use a Wolff-type law to regulate bone resorption and formation separately. Furthermore, mechanical feedback provides an effective mechanism to obtain physiological bone loss responses due to mechanical disuse and/or osteoporosis.

Spatial extension of a bone remodeling dynamics model and its finite element analysis
Abstract There are many works dealing with the dynamics of bone remodeling, proposing increasingly complex and complete models. In the recent years, the efforts started to focus on developing models that not only reproduce the temporal evolution, but also include the spatial aspects of this phenomenon. In this work, we propose the spatial extension of an existing model that includes the dynamics of osteocytes. The spatial dependence is modeled in terms of a linear diffusion, as proposed in previous works dealing with related problems. The resulting model is then written in its variational form, and fully discretized using the well‐known finite element method and a combination of the implicit and explicit Euler schemes. The numerical algorithm is then analyzed, proving some a priori error estimates and its linear convergence. Finally, we extend the examples already published for the temporal model to one and two dimensions, showing the dynamics of the solution in the spatial domain.

Dynamic Biochemical and Cellular Models of Bone Physiology: Integrating Remodeling Processes, Tumor Growth, and Therapy
Bone is an activeCoelho, Rui M. Neto, Joana P. connective tissue composed of different types of cells. The dynamic behavior of bone remodeling processes is typically represented throughValério, Duarte differential equations, which represent the physiological phenomena occurring in this organ. These models takeVinga, Susana into account the tight biochemical regulation between osteoclasts and osteoblasts and have also been enriched with variables and parameters related to bone pathologies and treatment. This chapter reviews some of the more recent models describing bone physiology, focusing on those that include the main cellular processes, along the biochemical control, and also the pharmacokinetic/pharmacodynamic (PK/PD) of the most common treatments for diseases such as cancer. These models are then compared in terms of the simulations obtained and, finally, some highlights on integrating them with the biomechanical component of the system which will be given. These models are expected to provide a valuable insight into this complex system and to support the development of clinical decision systems for bone pathologies.

Bone remodeling: A tissue-level process emerging from cell-level molecular algorithms
The human skeleton undergoes constant remodeling throughout the lifetime. Processes occurring on microscopic and molecular scales degrade bone and replace it with new, fully functional tissue. Multiple bone remodeling events occur simultaneously, continuously and independently throughout the body, so that the entire skeleton is completely renewed about every ten years.Bone remodeling is performed by groups of cells called Bone Multicellular Units (BMU). BMUs consist of different cell types, some specialized in the resorption of old bone, others encharged with producing new bone to replace the former. These processes are tightly regulated so that the amount of new bone produced is in perfect equilibrium with that of old bone removed, thus maintaining bone microscopic structure.To date, many regulatory molecules involved in bone remodeling have been identified, but the precise mechanism of BMU operation remains to be fully elucidated. Given the complexity of the signaling pathways already known, one may question whether such complexity is an inherent requirement of the process or whether some subset of the multiple constituents could fulfill the essential role, leaving functional redundancy to serve an alternative safety role. We propose in this work a minimal model of BMU function that involves a limited number of signals able to account for fully functional BMU operation. Our main assumptions were i) at any given time, any cell within a BMU can select only one among a limited choice of decisions, i.e. divide, die, migrate or differentiate, ii) this decision is irreversibly determined by depletion of an appropriate internal inhibitor and iii) the dynamics of any such inhibitor are coupled to that of specific external mediators, such as hormones, cytokines, growth factors. It was thus shown that efficient BMU operation manifests as an emergent process, which results from the individual and collective decisions taken by cells within the BMU unit in the absence of any external planning.
A review of recent developments in mathematical modeling of bone remodeling
In this article, we summarize the developments in the mathematical modeling of the mechanics of bone and related biological phenomena. We will devote special attention to the results of the last 10–15 years, although we will cover some relevant classical work to better frame the more recent researches. We will propose a division of the literature based on the main aim of the model (mechanical/biomathematical) and the type of biological phenomena considered (stimulus, growth, cell population dynamics). Finally, we will suggest some possible directions for future investigations.

A mechano-chemo-biological model for bone remodeling with a new mechano-chemo-transduction approach
Bone remodeling is a fundamental biological process that develops in bone tissue along its whole lifetime. It refers to a continuous bone transformation with new bone formation and old bone resorption that changes the internal microstructure and composition of the tissue. The main objectives of bone remodeling are: repair of the internal microcracks; adaptation of the macroscopic stiffness and strength to the actual changing mechanical demands; and control of the calcium homeostasis. Understanding this process and predicting its evolution is critical to reduce the effects of long-term disuse as happens during periods of reduced mobility. It is also important in the design of bone implants to avoid long-term stress shielding. Many mathematical models have been proposed from the earliest purely phenomenological to the latest that include biological knowledge. However, there still exists a lack of connection between the mechanical driving force and the biochemical and cell processes it triggers. Here, and following previous works that model independently the mechanobiological and biochemical processes in bone remodeling, we present a more complete model, useful for both cortical and trabecular bone, that uses a new mechanotransduction approach based on the effect of strains onto the bonding–unbonding rate of RANK/RANKL/OPG receptor–ligand reactions. We compare the results of this model with previous ones, showing a good agreement in similar conditions. We also apply it to realistic situations such as a femoral bone after implantation of a hip prosthesis, getting similar results to the clinical ones in the final bone density distribution. Finally, we extend this approach to the anisotropic case, getting not only the mean density, but also the directional homogenization of the microstructure. This biochemical approach permits, not only to predict the bone evolution under changes in the mechanical loads, but also, to consider anabolic and catabolic drugs to control bone density, such as those used in osteoporosis.

Simulated Interventions to Ameliorate Age-Related Bone Loss Indicate the Importance of Timing
Bone remodeling is the continuous process of bone resorption by osteoclasts and bone formation by osteoblasts, in order to maintain homeostasis. The activity of osteoclasts and osteoblasts is regulated by a network of signaling pathways, including Wnt, parathyroid hormone (PTH), RANKL/OPG and TGF-β, in response to stimuli such as mechanical loading. During aging there is a gradual loss of bone mass due to dysregulation of signaling pathways. This may be due to a decline in physical activity with age and/or changes in hormones and other signaling molecules. In particular, hormones such as PTH have a circadian rhythm which may be disrupted in aging. Due to the complexity of the molecular and cellular networks involved in bone remodeling, several mathematical models have been proposed to aid understanding of the processes involved. However, to date there are no models which explicitly consider the effects of mechanical loading, the circadian rhythm of PTH and the dynamics of signaling molecules on bone remodeling. Therefore, we have constructed a network model of the system using a modular approach which will allow further modifications as required in future research. The model was used to simulate the effects of mechanical loading and also the effects of different interventions such as continuous or intermittent administration of PTH. Our model predicts that the absence of regular mechanical loading and/or an impaired PTH circadian rhythm leads to a gradual decrease in bone mass over time which can be restored by simulated interventions and that the effectiveness of some interventions may depend on their timing.

A new biological bone remodeling in silico model combined with advanced discretization methods
Abstract Bone remodeling remains a highly researched topic investigated by many strands of science. The main purpose of this work is formulating a new computational framework for biological simulation, extending the version of the bone remodeling model previously proposed by Komarova. Thus, considering only the biological aspect of the remodeling process, the action of osteoclasts and osteoblasts is taken into account as well as its impact on bone mass. It is conducted a spatiotemporal analysis of a remodeling cycle obtaining a dynamic behavior of bone cells very similar to the biological process already described in the literature. The numerical example used is based on bone images obtained with scanning electron microscopy. During simulation, it is possible to observe the variation of bone's architecture through isomaps. These maps are obtained through the combination of biological bone remodeling models with three distinct numerical techniques—finite element method (FEM), radial point interpolation method (RPIM), and natural neighbor radial point interpolation method (NNRPIM). A study combining these numerical techniques allows to compare their performance. Ultimately, this work supports the inclusion of meshless methods due to their smoother results and its easiness to be combined with medical images from CT scans and MRI.

Bone refilling in cortical basic multicellular units: insights into tetracycline double labelling from a computational model
Bone remodelling is carried out by ‘bone multicellular units’ ($$\text{ BMU }$$s) in which active osteoclasts and active osteoblasts are spatially and temporally coupled. The refilling of new bone by osteoblasts towards the back of the $$\text{ BMU }$$occurs at a rate that depends both on the number of osteoblasts and on their secretory activity. In cortical bone, a linear phenomenological relationship between matrix apposition rate and $$\text{ BMU }$$cavity radius is found experimentally. How this relationship emerges from the combination of complex, nonlinear regulations of osteoblast number and secretory activity is unknown. Here, we extend our previous mathematical model of cell development within a single cortical $$\text{ BMU }$$to investigate how osteoblast number and osteoblast secretory activity vary along the $$\text{ BMU }$$’s closing cone. The mathematical model is based on biochemical coupling between osteoclasts and osteoblasts of various maturity and includes the differentiation of osteoblasts into osteocytes and bone lining cells, as well as the influence of $$\text{ BMU }$$cavity shrinkage on osteoblast development and activity. Matrix apposition rates predicted by the model are compared with data from tetracycline double labelling experiments. We find that the linear phenomenological relationship observed in these experiments between matrix apposition rate and $$\text{ BMU }$$cavity radius holds for most of the refilling phase simulated by our model, but not near the start and end of refilling. This suggests that at a particular bone site undergoing remodelling, bone formation starts and ends rapidly, supporting the hypothesis that osteoblasts behave synchronously. Our model also suggests that part of the observed cross-sectional variability in tetracycline data may be due to different bone sites being refilled by $$\text{ BMU }$$s at different stages of their lifetime. The different stages of a $$\text{ BMU }$$’s lifetime (such as initiation stage, progression stage, and termination stage) depend on whether the cell populations within the $$\text{ BMU }$$are still developing or have reached a quasi-steady state whilst travelling through bone. We find that due to their longer lifespan, active osteoblasts reach a quasi-steady distribution more slowly than active osteoclasts. We suggest that this fact may locally enlarge the Haversian canal diameter (due to a local lack of osteoblasts compared to osteoclasts) near the $$\text{ BMU }$$’s point of origin.
