A Systems Approach to Understanding Bone Cell Interactions in Health and Disease

1. Introduction

Bone is an important organ performing three essential physiological functions: mechanical support, mineral homeostasis (such as calcium and phosphate) and support of haematopoiesis. In fact, bone diseases in the elderly are associated with high morbidity and increased mortality. Osteoporosis and related skeletal complications are amongst the most important diseases impacting both the quality of life of our aging population and contributing costs to our health care system.

These different physiological functions of bone involve complex regulatory mechanisms at different spatial and temporal scales. For example, calcium and phosphate homeostasis is controlled on the whole organism scale and involves several organs such as parathyroid glands, gut, kidney and bone (Figure 1). On the other hand, bone structure and its mechanical properties are controlled on the tissue scale by the processes of bone modelling and bone remodelling. Bone modelling enables bones to adapt their shape during growth and in response to the prevailing mechanical loads. Bone remodelling removes micro-damage accumulating in the bone matrix during repeated mechanical loading. In this book chapter, we will focus on the implications of bone cell interactions in the control of bone remodelling and in the development of its disorders.

The renewal of bone during bone remodelling is achieved by a sequence of bone resorption and bone formation (Figure 2). This process both enables the repair of microcracks in the bone matrix and the balance between resorption and formation can be modulated for mineral homeostasis. The main cell types participating in bone remodelling are osteoblasts (bone forming cells), osteoclasts (bone resorbing cells) and osteocytes (load-sensing cells). During remodelling, these cell types are spatially and temporally organised in functional structures called basic multicellular units (BMUs). The remodelling sequence operated by a BMU follows well-defined phases. It starts with an activation phase consisting of recruitment of precursor cells. This is followed by a resorption phase characterised by active osteoclasts (OC_as) removing bone matrix and a formation phase where active osteoblasts (OB_as) lay down osteoid which then becomes mineralised to form new bone matrix. After completion of the local renewal of the bone matrix, bone cells become quiescent (resting phase) until a new remodelling event is initiated.

Figure 1.

Schematic representation of whole body calcium regulation and local bone tissue regulation via bone remodelling (figure from Peterson and Riggs 2010, reproduced with permission).

Over the past decades, a large number of regulatory factors produced by hormonal glands (such as parathyroid- and thyroid glands), lymphocytes, bone cells involved in sensing the mechanical environment (such as osteocytes and lining cells), and even tumour cells, have been shown to influence the bone remodelling sequence. These regulatory factors are essential components of the cell-cell signalling network between the bone-resorbing and bone-forming cells, and will be briefly reviewed in Section 2. Many bone disorders such as osteoporosis, Paget’s disease and cancer-related bone diseases are related to disruptions in the biochemical or cellular components of this signalling network. These disruptions lead to imbalances between bone resorption and bone formation in the BMU remodelling sequence, and/or to changes in bone turnover as expressed by the activation frequency of BMUs (Riggs and Parfitt 2005).

Figure 2.

Bone remodelling sequence executed by different bone cell types within Basic Multicellular Units (BMUs).

Bone diseases are often multifactorial and can exhibit high inter-individual variability. It is of paramount importance to characterise and treat these diseases in a patient-specific way. Different patients exhibit different temporal evolutions of bone mass and bone cell populations. These differences can be used to define a patient-specific “disease signature”, for example by the quantitative characterisation of both bone resorption rate and bone formation rate (Figure 3). The quantitative definition of such a signature is of prime importance for mathematical modelling. Mathematical modelling in biology in the next decades is expected to develop into clinical tools to help predict the evolution of a disease in an individual and to find its optimum treatment regimes. In this contribution, a mathematical model of bone cell interactions during bone remodelling is presented. This model is applied to simulate a catabolic bone disease (e.g. osteoporosis) and to investigate various treatment strategies.

A major challenge in bone biology is to understand the spatio-temporal mechanisms of action of cell-cell signalling and the interdependence between the bone structure and the bone cells actively modifying this structure at different time and length scales. Cell-cell and cell-structural interactions are complex, forming various feedback loops that are essential for bone homeostasis. Disruption of cell-cell and cell-structure feedback loops can lead to bone fractures (Chavassieux, Seeman et al. 2007). In osteoporosis for example, an important signalling pathway (i.e., the RANK-RANKL-OPG pathway – see Section 2) is strongly activated, leading to an excess of resorption over formation during remodelling (Vega, Maalouf et al. 2007). Growth factors released from the bone matrix during bone resorption further stimulate the resorptive (catabolic) pathway in a positive feedback loop, which increases remodelling activity and bone loss. To maintain bone homeostasis, formative (anabolic) pathways need to be stimulated to counterbalance the bone catabolic responses, and so to maintain skeletal integrity.

Figure 3.

Rate of bone formation [%/day] versus rate of bone resorption [%/day] diagram: utilized to represent bone disease signature; red area indicates catabolic regimes while blue area indicates anabolic regimes; points on diagonal indicate bone homeostasis and rate of bone turnover.

A systems approach to bone remodelling can help advance our current understanding of the complex system formed by the bone cells and their signalling network (Cassman, Arkin et al. 2005; Smallwood 2011). Many cellular and bio-molecular behaviours have been discovered by purely experimental research. This has led to a ‘profusion’ of biological data, summarised in detailed “component network models”. However, bone remodelling is a highly dynamic system with interdependent regulatory controls operating on multiple components. These interconnections form feed-forward and feedback controls and may lead to system behaviours that cannot be easily predicted from the behaviour of the individual components (e.g. non-linear bifurcations, instabilities, snap-through behaviours) (Strogatz 2001). Understanding how this biological system functions as a whole requires a quantitative analysis of its dynamic network behaviour by mathematical and computional modelling.

Computational approaches provide a powerful integrative methodology to investigate the complex system behaviour of bone remodelling, enabling quantitative systematic testing of various experimental and theoretical hypotheses in-silico (Pivonka and Komarova 2010; Webster and Müller 2011). As such, computational modelling can be considered an additional methodology to improve our understanding of bone biology in the same way as in-vitro cell culture experiments and in-vivo animal models are experimental methodologies. In biology, in-vitro and in-vivo models are used as controlled simplifications of the studied processes. Most experimental systems investigate physiological extremes such as maximal inhibition of specific pathways, complete ablation of a gene, or multi-fold gene over-expression. The relevance of the data obtained in these experimental manipulations for human physiology relies on the assumptions that a) the modeled processes are fundamentally similar to those occurring in humans; b) the role of the quantities of interest that are studied or manipulated is more important than the role of the quantities that are kept constant or that are not measured; and c) the effects of extreme changes reflect normal physiology.

The main drawbacks and limitations of mathematical modelling are due to the complexity of the biological systems under investigation, which requires sophisticated models. Various classes of mathematical models can be used depending on the length and time scales to investigate. As a general rule, an increase in model complexity is associated with an increase in the number of model parameters. This raises the important question of calibrating and validating a mathematical model (Babuska and Oden 2004). Many parameters introduced in mathematical models have not yet been measured or are out of experimental reach using currently available methods. For these reasons, these mathematical models are often qualitative. The development of more quantitative models will rely on the collaborative efforts of experimental biologists and applied mathematicians/bioengineers.

This book chapter is organized as follows: In Section 2, a brief overview of bone biology is presented, with emphasis on the biochemical regulation of bone remodelling, bone diseases and currently available therapeutic interventions. In Section 3, we first present a general approach to modelling a system of cell interactions mediated by signalling molecules based on receptor-ligand binding reactions. This approach is then used to present a comprehensive bone cell population model for bone remodelling. In Section 4, we apply this model to simulate a catabolic bone disease (osteoporosis) and investigate different in-silico treatment strategies. In particular, we apply the receptor-ligand binding reaction scheme to include the effect of the drug denosumab, a monoclonal antibody to RANKL. Finally, we provide a summary and conclusions in Section 5.

2. Bone biology background

Bone has multiple functions within the body with the most important being to provide support to muscles, ligaments, tendons and joints to enable movement. For this function the bone needs to maintain structure, strength, and rigidity without compromising lightness. Secondly bone provides a readily accessible store of calcium to support calcium homeostasis. Other roles include providing a protected environment for bone marrow and support for a haematopoietic stem cell niche (Calvi, Adams et al. 2003), and acting as an endocrine organ regulating energy metabolism, perhaps through osteocalcin (Clemens and Karsenty 2011).

To achieve the above functions a number of cellular processes are required. For the structural demands, there is a need for a control system for the detection of bone strain and microdamage, and signalling systems for an appropriate reparative response, which would involve the directed induction of bone resorption to remove damaged bone and new bone formation to replace it (Kidd, Stephens et al. 2010). To maintain the structural integrity of bone, signalling processes are required for the coordination of bone resorption and bone formation. In fact, on the inner endosteal surfaces of bone and within Haversian systems within cortical bone, there is strong spatial and temporal coupling of bone resorption and bone formation. Bone resorption is initially induced and bone formation follows such that, by the end of this remodeling cycle, an approximate balance is achieved between the amounts of bone removed and replaced (Jaworski 1984). In contrast on periosteal bone surfaces where the main function is modelling bone, bone resorption tends to be associated with growth related shaping of the bone, and intermittent cycles of bone formation contribute to a gradual expansion in total bone diameter without the requirement for prior bone resorption.

Bone contributes also to calcium and phosphate homeostasis as an accessible internal store for these essential minerals. To fulfill the body’s calcium homeostatic demands for a tightly controlled blood level of calcium ions, three organs respond to systemic hormones to provide complementary mechanisms for increasing or decreasing blood calcium levels. The loss of calcium from the kidney in urine can be altered by changing the rate of calcium reabsorption from the renal tubules. The proportion of calcium and phosphate absorbed from food can be changed, and the calcium and phosphate flux in and out of bone can be altered. Parathyroid hormone (PTH) is secreted by the parathyroid gland when blood calcium levels are low and acts on the kidney to increase calcium reabsorption from the distal tubules to decrease loss of calcium in the urine. Simultaneously it decreases phosphate reabsorption resulting in increased urinary phosphate loss. PTH also acts on bone to increase bone resorption by osteoclasts and may also induce a calcium and phosphate flux from labile calcium phosphate within bone and in particular that deposited around osteocytes (Parfitt 1976; Talmage and Mobley 2008). PTH also acts to induce activation of vitamin D in the kidney to produce 1,25 dihydroxyvitamin D which acts primarily as a hormone to increase calcium and phosphate absorption from the gut, and can increase bone resorption by osteoclasts to increase calcium and phosphate release from bone. Calcitonin is produced by the thyroid gland in response to elevated blood calcium levels and acts directly on osteoclasts to inhibit their activity thus reducing calcium and phosphate release from bone by osteoclasts though this effect is reversible (Chambers 1982). It also may be important in protecting the skeleton during prolonged bone stress due to pregnancy and lactation (Hurwitz 1996; Woodrow, Sharpe et al. 2006).

As conflicting structural and homeostatic signals may be concurrently present, there is a requirement for the ability to integrate multiple signals such that responses are tailored to the particular local bone environment; for example less structurally important bone should be targeted for calcium mobilization during times of high calcium demand as observed as preservation of trabecular rods seen in the maternal skeleton during lactation (Liu, Ardeshirpour et al. 2012).

To achieve these diverse requirements there are three highly specialized bone cells present in bone:

1. Osteocytes reside within the bone matrix and interact with bone and communicate to adjacent osteocytes via a dense and interconnecting canaliculae network containing osteocyte cell processes. Osteocytes orchestrate the detection and cellular responses to microdamage, fracture and changing bone strain and are the bone cells present in by far the largest numbers, comprising 90% of all bone cells (Lanyon 1993).

2. Osteoblasts are cells of mesenchymal origin which line the surface of bone and are capable of bone formation by progressive bone matrix deposition and subsequent mineralization. Cells of this lineage have receptors for PTH and 1,25 dihydroxyvitamin D and a number of local regulatory factors (Rodan and Martin 1981). After completion of their role in bone formation, the final fate for these cells is either cell death or further differentiation to osteocytes, which occurs concurrently with encasement within bone matrix.

3. Osteoclasts are large multinucleated cells of haematopoietic origin able to resorb bone by adhering to bone surfaces, secreting acid to demineralise the bone and proteolytic enzymes to break down the collagenous bone matrix. However, they are typically unable to respond directly to pro-resorptive hormones and require the presence of osteoblast lineage cells to locally regulate their differentiation and activity.

Communication between osteocytes, osteoblasts and osteoclasts enables a spatially and temporally coordinated and directed response through the integration of multiple catabolic and anabolic signals leading to skeletal responses to both physiological demands and pathological challenges.

2.1. Control of bone resorption

Cells of the osteoblast lineage maintain a pool of osteoclast precursors through the constitutive secretion of M-CSF which is a differentiation and survival factor for early osteoclast precursors which express the m-CSF receptor c-fms (Wiktor-Jedrzejczak, Bartocci et al. 1990). This pool of cells is then available for recruitment to the osteoclast population when needed to provide capability for either bone modeling (sculpting of bone shape), or remodeling (renewal of bone) or for initiating bone resorption to release calcium.

2.1.1. RANK-RANKL-OPG pathway

The regulation of bone resorption is through integration of multiple pro- and anti-resorptive stimuli with convergence of output into a dominant mediating pathway. Cells of the osteoblast lineage integrate multiple pro-and anti-resorptive signals, including hormonal, mechanical and pathological signals, with output through the changing balance of their expression of a cytokine, receptor activator of NFkappaB ligand (RANKL) (Yasuda et al 1998, Lacey et al 1998), and its inhibitor, osteoprotegerin (OPG) (Simonet, Lacey et al. 1997; Tsuda, Goto et al. 1997). RANKL is expressed as either a membrane bound cytokine or released in a soluble form by cells of the osteoblast lineage. Osteocytes, osteoblasts and osteoblast precursors can express RANKL and there is at present some controversy regarding which predominates. Selective knockdown of RANKL in osteocytes produces an apparently moderately severe osteopetrosis in mice indicating a significant role for osteocytes as a source for RANKL (Nakashima, Hayashi et al. 2011). RANKL binds to its receptor RANK (Nakagawa, Kinosaki et al. 1998) on the surface of osteoclast precursors inducing a number of intracellular signalling pathways to drive differentiation to an osteoclast phenotype, activating osteoclastic bone resorption and increasing osteoclast survival. OPG is a decoy receptor for RANKL, and if secreted in excess by cells of the osteoblast lineage, will bind RANKL and prevent its association with RANK, thereby inhibiting osteoclast differentiation and activity and promoting osteoclast apoptosis. Thus excess expression of RANKL relative to OPG will promote bone resorption and an excess of OPG relative to RANKL will inhibit bone resorption (http://www.rankligand.com/).

Typically, pro-resorptive hormones such as PTH and calcitriol will increase RANKL expression and decrease OPG expression within bone leading to increased bone resorption, while bone protective agents such as estradiol and testosterone tend to increase the expression of OPG relative to RANKL thus reducing bone resorption (Horwood, Elliott et al. 1998; Michael, Härkönen et al. 2005). Further, PTHrP, IL-1 and TNFα are examples of local regulatory factors which increase RANKL expression by cells of the osteoblast lineage during pathological bone loss related to cancer (PTHrP) or inflammation (IL-1 and TNFα). Local regulatory factors such as Mechanical signals, in contrast to hormones, are locally expressed and alter RANKL and OPG expression in a spatially restricted manner. Reduced bone strain or the presence of microfractures or fatigue damage, increase expression of RANKL relative to OPG, while increased bone strain tends to decrease the expression of RANKL, relative to OPG expression, by osteocytes (Kulkarni, Bakker et al. 2010). The pro-resorptive hormones and other local regulatory factors increase RANKL by binding to their own specific receptors in osteoblasts, activating the specific signalling pathways associated with these receptors, leading, through transcriptional activation of the RANKL promoter, to increased RANKL expression.

While RANKL/RANK signalling is necessary and dominant in the regulation of bone resorption, other modulating signalling molecules can signal directly to osteoclasts to magnify or diminish their responses to RANKL. Examples are inflammatory cytokines which can magnify osteoclast responses (Lam, Takeshita et al. 2000) and the systemic hormone calcitonin which is a potent but reversible direct inhibitor of osteoclast activity via activation of the calcitonin receptor on mature osteoclasts (Chambers 1982). Other sources of RANKL are also possible in bone but these tend to contribute to the regulation of bone resorption during disease. Activated t-cells in particular secrete RANKL and can induce bone resorption during infection or chronic inflammation (Teng, Nguyen et al. 2000).

2.1.2. Current therapeutic interventions of bone resorption – anti-resorptive drugs

There are a number of treatments currently clinically available and in development for inhibiting bone resorption. These target the osteoclast with different treatments inhibiting osteoclast differentiation, osteoclast function, osteoclast survival or a combination of these. The bisphosphonate family of drugs have been available for many years and act by binding tightly to bone mineral, being internalized by osteoclasts by endocytosis during bone resorption where they act to inactivate osteoclast function and to induce osteoclast apoptosis. Bisphosphonate treated bone is characterized by reduced numbers of osteoclasts, most of which are withdrawn from the bone surface, and many of which are undergoing apoptosis (programmed cell death). The later nitrogen containing bisphosphonates such as zoledronate are highly potent and long acting inhibitors of osteoclast function, with dosing frequency of up to one year being effective (Kavanagh, Guo et al. 2006).

Calcitonin has also been available for many years and acts by reversibly inactivating osteoclasts, however it lacks the sustained activity to be effective where the pathological stimuli for bone resorption are present such as when metastatic cancer foci are present in bone, though it is used as a treatment for postmenopausal osteoporosis or when rapid suppression of bone resorption is required (Mazzuoli, Passeri et al. 1986 ).

A recently approved treatment for osteoporosis and metastatic bone disease is the RANKL/RANK signalling inhibitor denosumab, which is a monoclonal antibody to RANKL which mimics the action of OPG by binding RANKL to prevent its association with RANK. This protein acts to prevent osteoclast differentiation as RANKL signalling is an essential requirement for osteoclast formation. It also acts to reduce osteoclast survival as RANKL is an important survival factor for osteoclasts. Histologically bone treated with denosumab is characterized by a marked reduction in osteoclast numbers (McClung, Lewiecki et al. 2006).

There are a number of treatments in development which target osteoclast function. Inhibitors of the protease cathepsin K inhibit the ability of osteoclasts to break down bone collagen and thus inhibit their function. Bone of animals treated with cathepsin K inhibitors are characterized by the presence of many osteoclasts, and shallow resorption lacunae lined by a layer of demineralised matrix reflecting the ability of osteoclasts to demineralise bone but not to break down the matrix (Eastell, Nagase et al. 2011).

Table 1.

Regulation of the relative expression of RANKL and OPG by a selection of systemic and local factors (from (Horwood, Elliott et al. 1998), (Tsukii, Shima et al. 1998), (Lee and Lorenzo 1999), (Hofbauer, Dunstan et al. 1998), (Brändström, Jonsson et al. 1998), (Hofbauer, Gori et al. 1999), (Hofbauer, Khosla et al. 1999), (Hofbauer, Lacey et al. 1999), (Vidal, Sjögren et al. 1998), (Michael, Härkönen et al. 2005))

2.2. Control of bone formation

In contrast to the regulation of osteoclasts, regulation of osteoblasts and bone formation is more complex with evidence of extensive signal duplication and redundancy. To control bone formation there are a large number of growth factors with a range of actions on cell behaviour. There are several classes of growth factors with a range of potent actions on osteoblasts and their precursors, and each class typically having multiple members and receptors, and a significant number of signalling inhibitors/enhancers. The regulatory functions of these can be divided into the following:

1. Inducing proliferation of the osteoblast precursor population.

2. Directing lineage commitment of pluripotent mesenchymal stem cells into the osteoblast lineage.

3. Inducing differentiation of osteoblast precursors into mature osteoblasts.

4. Regulation of osteoblast activity including both collagen deposition and matrix mineralization.

5. Directing further differentiation of osteoblasts to osteocytes with concurrent encasement in the bone matrix.

Regulation of bone formation occurs through the sequential action of growth factors favouring these functions in a signalling cascade characterized by feed forward and feed back paracrine and autocrine signalling.

2.2.1. Proliferative agents

Proliferative agents act as priming agents to increase the population of precursors available differentiate into functioning osteoblasts. Important proliferative growth factors include transforming growth factor beta (TGF-ß) (Dallas, Rosser et al. 2002), fibroblast growth factors (FGF), particularly FGF-1 and FGF-2 (Dunstan, Boyce et al. 1999), and platelet derived growth factors (PDGF) particularly PDGF-BB (Caplan and Correa 2011). These factors induce proliferation of osteoblast precursors in bone but tend to inhibit osteoblast differentiation. Thus their role is important in the early stages of osteoblast generation when the critical requirement is to generate a sufficiently large pool of osteoblast precursors from the relatively small number of mesenchymal stem cells present. Levels of these agents are thus increased transiently in the bone environment during bone resorption and after fracture, but progressive differentiation of precursors through to fully mature osteoblasts requires reduction in their levels.

The source of these factors in bone varies depending on the site and the nature of the stimulus for bone formation. Bone injury leads to infiltrates in the bone of inflammatory cells which are the source of PDGF-BB, which in turn is essential to expand the mesenchymal precursor pool for effective fracture repair. In contrast, in endosteal bone in which bone formation is required to balance bone removed by osteoclasts recruited for bone remodeling, a primary source of proliferative signals is the bone matrix itself. Osteoblasts, when forming bone, produce growth factors that are sequestered within the bone matrix. In the case of TGF-β this is as an inactive pro-form of the protein. However osteoclasts activate and release TGF-β during bone resorption, enabling it to be available to induce the proliferation of osteoblast precursors (Dallas, Rosser et al. 2002).

FGF’s which are highly proliferative, and IGF’s and BMP’s which have limited proliferative activity but are stimulators of osteoblast differentiation are also sequestered in bone with potential for release during bone resorption.

2.2.2. Lineage commitment factors

The wnt family of growth factors and their respective receptors, co-factors and inhibitors provide one source of regulation of lineage commitment of mesenchymal stem cells and the population of osteoblast precursors expanded by the proliferative agents above. There is evidence that mature osteoblasts can produce and secrete Wnt’s and these can act on early mesenchymal precursors to induce commitment to differentiation in the osteoblast lineage. Determining the details of wnt action has been complicated by the presence of 19 wnt ligands, 10 wnt receptors not to mention multiple wnt receptor co-factors and inhibitors (Kawano and Kypta 2003). BMP’s which are produced by osteoblast lineage cells are sequestered in bone matrix may also contribute to lineage commitment for osteoblast precursors (Abe, Yamamoto et al. 2000), though this is difficult to separate from their differentiation activity. During the inflammation such as occurs during fracture repair, transient exposure to TNF-α may influence commitment of mesenchymal stem cell to the osteoblast lineage (Lu, Wang et al. 2012).

2.2.3. Differentiation factors

The BMP’s, particularly BMP2, BMP4 and BMP7, are particularly potent differentiation agents with only limited proliferative activity. IGF1 acts to increase osteoblast activity and to increase survival and is also found in bone matrix. The activity of IGF’s is modulated by various binding proteins while that of BMP’s can be modulated by specific inhibitors such as noggin, chordin and gremlin. They are sourced in bone from osteoblastic cells as autocrine and paracrine factors, and are released from bone matrix during bone resorption. They have an important role directing the differentiation of osteoblast precursors into mature bone forming osteoblasts (Gazzerro and Canalis 2006). This involves migration of precursors to the bone surface, polarization of the cell, activation of expression of collagen type 1 and other bone matrix proteins and the development of extensive intracellular machinery for protein manufacture and secretion. Parathyroid hormone related protein (PTHrP) is a paracrine factor produced predominantly by early osteoblast precursors which acts to induce osteoblast differentiation and survival and also acts as a survival factor for osteocytes (Martin 2005; Michael, Härkönen et al. 2005). PTHrP can also mimic PTH actions to increase RANKL expression and thus also increase bone resorption. When PTHrP is produced in excessive amounts by cancer cells, it can cause profound bone loss (Suva, Winslow et al. 1987). Regulation of further differentiation of osteoblasts to osteocytes is not currently well understood, but presumably is also closely regulated by paracrine or autocrine differentiation signals.

Table 2.

Families of growth factors active in osteoblast regulation (from (Li, Hassan et al. 2008), (Rodda and McMahon 2006), (Zhou, Mak et al. 2008), (Abe, Yamamoto et al. 2000), (Tanaka, Ogasa et al. 1999), (Kawaguchi, Pilbeam et al. 1995), (Dallas, Rosser et al. 2002; Martin 2005; Miao, He et al. 2005))

Regulation of the differentiation of periosteal osteoblasts is likely to differ with regards to growth factors involved compared to those for osteoblasts in Haversian systems or on endosteal surfaces as the formation of periosteal osteoblasts is not preceded by osteoclastic bone resorption and thus growth factors stored in the bone matrix are not released to prime the osteoblast precursor population.

2.2.4. Current therapeutic agents to stimulate bone formation – anabolic drugs

There are a number of treatments that have been found to promote bone formation but it has been found essential that treatments preserve bone material properties and structure in order for improvement in bone strength. Fluoride is a potent anabolic factor used clinically for many years treating osteoporosis but was found to cause the formation of bone in which bone matrix was disorganized and poorly mineralized, and trabeculae were lacking in normal architecture. Its use significantly increased bone mineral density but was actually associated with increased fracture rates (Gutteridge, Price et al. 1990; Gutteridge, Stewart et al. 2002). Currently the most effective anabolic factor is parathyroid hormone (teriparatide) which fulfils the requirement for a factor which, at the doses used, induces the formation of new bone which retains both its normal lamellar structure and trabecular architecture and is treatment is associated with decreased risk of fracture (Paschalis, Glass et al. 2005). Interestingly, PTH is catabolic if continually present, but anabolic when given intermittently (by daily injection) (Onyia, Helvering et al. 2005). Strontium has been claimed to be a modestly anabolic factor when given orally as strontium ranelate and is approved for the treatment of osteoporosis due to it fracture prevention effects (Reginster 2002). However the mechanism of this action is not well determined.

Anabolic treatments in development include a sclerostin neutralizing antibody which in early trials promoted bone formation and increased bone mineral density (Padhi et al 2011). Presumably this antibody treatment acts by reducing the osteocyte mediated down-regulation of bone formation through the expression of the Wnt signalling inhibitor sclerostin. Interestingly one of the actions of intermittent PTH is to reduce osteocyte production of sclerostin and so PTH and the sclerostin neutralizing antibody may work in part by similar stimulatory actions on Wnt signalling (Keller and Kneissel 2005).

3. Mathematical model of bone cell interactions

The above review of some of the most important signalling pathways involved in the coupling of bone resorption and bone formation in bone remodelling indicates the complexity of the communication network between bone cells. While rapid advancement of experimental testing techniques generates large amounts of information on bone remodelling, a challenging problem is to integrate this experimental data into an understanding of system behaviour and to connect observations made at different time and length scales. In the following, an overview of mathematical models employed to integrate this different experimental information into a systems understanding of bone regulation is provided. The complexity of biological systems such as bone is characterized by the respective spatial and temporal scales involved. Due to the large variation of spatio-temporal scales involved in bone remodelling only subsets of these scales are generally investigated depending on the biological questions addressed. Mathematical modelling has been applied at the following scales:

On the subcellular scale one may be interested in signal transduction mechanisms such the NF-κB signalling pathway that becomes activated once RANKL binds to its receptor RANK expressed on osteoclastic cells and so triggers gene expression and consequently osteoclastic activity and cell survival. NF-κB signalling plays a pivotal role in the pathogenesis of osteolytic bone disorders. The NF-κB signalling pathway is regulated to maintain bone homeostasis by cytokines such as RANKL, TNF-α and IL-1. NF-κB activation involves stimulus-induced degradation of its inhibitor IkB, which allows for its translocation to the nucleus. Several mathematical models have been developed aimed at understanding how NF-κB translocation and IkB association/dissociation rate constants keep the majority of NF-κB in an inactive state in resting cells (Kearns, Basak et al. 2006).

At the cell and single BMU scale one is concerned with modelling different cell behaviour such as resorption and formation properties of bone cells and respective features of the BMU such as resorption speed, number, type and distribution of cells in the cutting and closing cone, blood vessel, etc.. A number of mathematical models have been developed to investigate single BMU behavior using either discrete (agent-based) approaches (Cacciagrano, Corradini et al. 2010; Buenzli, Jeon et al. 2012) or continuous approaches based on partial differential equations (Ryser, Komarova et al. 2010; Buenzli, Pivonka et al. 2011).

On an even larger scale in the following referred to as tissue scale one may want to characterize how bone mass changes with time during osteoporosis and other bone diseases. Such data can readily be obtained using micro-computed tomography (micro-CT) which relates measured quantities such as bone volume, bone surface area, trabecular numbers and spacing to a certain representative volume of bone tissue (such as 2x2x2 mm³). Such as a volume generally contains several BMUs and through serial sectioning histology the number of bone cells in that volume can be estimated. In mathematical terms the size of the chosen volume is crucial for defining average quantities (both mechanical and biological). A volume size where a continuous mathematical description is possible is generally denoted as representative volume element (RVE). The characterisation of an RVE requires a separation of scales between microscopic constitutive elements or components and the scale of interest at which properties are represent as continuous quantities. For example estimating the local mechanical stiffness of bone requires that the RVE contains sufficiently many pores and bone matrix, yet is small enough to represent spatial inhomogeneities (such as due to local bone morphology) and to be considered infinitesimal at this scale. In this context, different mathematical models of bone cell dynamics have been proposed by Komarova et al. (Komarova, Smith et al. 2003), Lemaire et al. (Lemaire, Tobin et al. 2004) and Pivonka et al (Pivonka, Zimak et al. 2008).

On the whole organ scale, i.e., the skeletal scale one may be interested in the action of muscles on bones during locomotion. These types of models are commonly referred to as musculoskeletal models and allow analyzing the forces produced in bones, cartilage and ligaments due to contraction and extension of muscles together with establishing in-vivo loading conditions (Shelburne and Pandy 1997; Shim, Hunter et al. 2011).

On the largest scale, i.e., the whole organism scale one may be concerned with understanding mineral homeostasis. As has been pointed out in Section 2 bone plays an important role as a reservoir of calcium and phosphate. Several mathematical models have been proposed to describe calcium and phosphorous homeostasis considering different hormonal regulators such as PTH and calcitonin (Raposo, Sobrinho et al. 2002; Peterson and Riggs 2010).

The following applications are concerned with temporal aspects of bone cell interactions using a representative volume element containing several BMUs. Using this scale of representation we investigate changes in bone matrix (volume fraction) and bone cell numbers over time both under pathological conditions and using different therapeutic interventions.

3.1. Modelling cell-cell interactions and cell responses

Receptor-ligand. binding, cell signalling and cell response

Bone balance and bone turnover during remodelling are critically dependent on the coordination between osteoclasts and osteoblasts. This coordination is carried out by cell–cell interactions mediated by several signalling molecules summarized in Section 2. In this section, we present a generic approach to modelling such cell–cell communication. This approach is based on fundamental chemical reaction principles of material balance and mass action kinetics, and can be applied to other biological systems (Lauffenburger and Linderman 1996; Lemaire, Tobin et al. 2004; Pivonka, Zimak et al. 2008).

Extracellular ligands (e.g. systemic hormones, autocrine and paracrine factors, growth factors etc.) modulate cell behaviour by binding to specific receptors on cells and activating intracellular signalling pathways. Intracellular signalling mechanisms can be complex and may include a cascade of interconnected downstream pathways involving protein trafficking, translocation to the nucleus, gene transcription factors, protein synthesis, etc. These intracellular mechanisms lead to an overall cell response, such as cell differentiation, proliferation, apoptosis or the expression of signalling molecules or receptors (see Figure 4). Both extracellular and intracellular signalling mechanisms are governed by the chemical reaction principles of material balance and mass action kinetics (Lauffenburger and Linderman 1996). To focus on the modelling of cell-cell interactions mediated by extracellular signalling molecules, the approach followed in this book chapter will treat extracellular and intracellular signalling differently. Extracellular signalling will be modelled explicitly by considering receptor–ligand binding reactions governed by mass action kinetics. By contrast, intracellular signalling will be modelled phenomenologically by assuming that it leads to an overall relation between the "input signal" perceived by a cell via extracellular receptor-ligand binding and the cell response (Figure 4).

Figure 4.

Modelling of cell responses using receptor ligand binding reactions; phenomenological response functions can either promote (π^act) or repress (π^rep) a certain cell response; relationship between input signal and output signal can be linear (α=1) or non-linear (α≠1).

Extracellular signalling: competitive binding of receptors and ligands

Receptors enable a cell to sense its local biochemical environment. If a ligand binds to a cell surface receptor and particularly if it becomes endocytosed (e.g. activating intracellular signalling mechanisms), that ligand becomes unavailable to signal to the other cells. We refer to this effect as "competitive binding". Competitive binding influences the amount of ligand binding to a single cell in presence of other cells, and so influences the "input signal" that a cell may perceive from its micro-environment. An advantage of considering explicit extracellular receptor–ligand binding reactions in a mathematical model is to account for such competitive binding in a systematic way. Competitive binding may occur either amongst identical receptors (e.g. RANKL binding RANK on several osteoclasts), or amongst different kinds of receptors binding the same ligand (e.g. RANKL binding either RANK or OPG). In the following, we first consider a simple receptor–ligand chemical reaction to illustrate the concept of competitive binding amongst identical receptors. We will then present a more general situation applicable to the RANK–RANKL–OPG signalling system, where there may be competitive binding between several types of RANKL (e.g. expressed on osteoblasts, on osteocytes, or in soluble form) and several types of receptors, (e.g. RANK, OPG and denosumab).

We first consider a single ligand L, assumed to be produced at rate P_L and degraded at rate D_L. This ligand may bind to a single receptor R, produced at rate P_R and degraded at rate D_R (production and degradation may include the generation and apoptosis of cells expressing the ligand or receptor). After binding, the receptor and ligand form a bound complex, which we denote by$$\stackrel{⌢}{RL}$$. We consider that this complex may either unbind (reverse reaction) or may be endocytosed, in which case both the receptor and the ligand become unavailable for further extracellular reactions. This case effectively corresponds to a disappearance of the bound complex from the extracellular environment. The chemical reaction flows corresponding to this receptor-ligand binding situation can thus be summarised as:

$$\underset{D_R\downarrow }{\overset{P_R\downarrow }{{\displaystyle R}}}+\underset{D_L\downarrow }{\overset{P_L\downarrow }{{\displaystyle L}}}\underset{k_{RL}^r}{\overset{k_{RL}^f}{\rightleftharpoons }}\stackrel{⌢}{RL}\underset{\text{(endocytosis)}}{\overset{D_{RL}}{\to }}\varnothing$$E1

where $k_{RL}^f,k_{RL}^r>0$are the forward and reverse binding reaction rates and $D_{RL}$is the rate at which the bound complex is degraded from the extracellular environment by endocytosis. Assuming first order reaction rates (i.e. rates proportional to the amount of reactant) and the law of mass action, the rate equations corresponding to this receptor–ligand binding situation are:

$$\begin{array}{l}{\partial }_t\stackrel{⌢}{RL}=k_{RL}^{f}R\cdot L-k_{RL}^r\stackrel{⌢}{RL}-D_{RL}\stackrel{⌢}{RL}\\ {\partial }_{t}R=P_R-D_{R}R+k_{RL}^r\stackrel{⌢}{RL}-k_{RL}^{f}R\cdot L\\ {\partial }_{t}L=P_L-D_{L}L+k_{RL}^r\stackrel{⌢}{RL}-k_{RL}^{f}R\cdot L\end{array}$$E2

where the symbols R, L and $\stackrel{⌢}{RL}$ denote the free (unbound) receptor, free ligand and bound complex concentrations (number of molecules per unit volume, e.g. mol/L). Receptor–ligand binding reactions occur on a fast time scale compared to characteristic times of cell behaviours. Local concentrations of receptors and ligands quickly converge to a quasi-steady state (Lemaire, Tobin et al. 2004; Pivonka, Zimak et al. 2008; Buenzli, Pivonka et al. 2012). This steady-stated is determined by setting time derivatives to zero in the system of ordinary differential equations (ODEs) Eqs. (2), which leads to the following system of algebraic equations:

$$R^{st}=\frac{P_R}{D_R+\frac{D_{RL}}{k_{RL}}{L}^{st}},L^{st}=\frac{P_L}{D_L+\frac{D_{RL}}{k_{RL}}{R}^{st}},{\stackrel{⌢}{RL}}^{st}=\frac{1}{k_{RL}}{R}^{st}{L}^{st}$$E3

where

$$k_{RL}=\frac{k_{RL}^r+D_{RL}}{k_{RL}^f}$$E4

is a parameter specific to the binding reactions (1). Equation (3) indicates that competitive binding effects only occur if there is a degradation of the bound complex, i.e. if D_RL > 0. Indeed, if D_RL = 0 (no degradation of the bound complex by endocytosis), then the steady-state concentration of free receptors is independent of the amount of free ligands and the steady-state concentration of free ligands is independent of the amount of free receptors. In fact, without bound complex endocytosis, the assumed production and degradation of the receptors and ligands act as "equilibrating reservoirs". This also explains why the steady-state concentrations do not depend on the initial amount of reactants.

The above description can be generalised to a situation in which a number of different kinds of receptors R = R₁,R₂,... may bind a number of different kinds of ligands L = L₁, L₂,.... The binding reactions (1) hold for all R and L and lead to the following system of algebraic equations for the steady-state concentrations of receptors, ligands and bound complexes (dropping the steady-state superscript "st" from the notation):

$R=\frac{P_R}{D_R+{\displaystyle {\sum }_L\frac{D_{RL}}{k_{RL}}L}},L=\frac{P_L}{D_L+{\displaystyle {\sum }_R\frac{D_{RL}}{k_{RL}}R}},\stackrel{⌢}{RL}=\frac{1}{k_{RL}}R\cdot L$E5

for all R = R₁,R₂,... and L = L₁, L₂,.... As before, competitive binding effects may occur provided the bound complexes are degraded (i.e. D_RL > 0). Indeed, in that case the existence of several receptors to bind the same ligand lowers the free ligand concentration and the existence of several ligands to bind the same receptor lowers the free receptor concentration. We note here that the production and degradation rates of the ligands and receptors may depend themselves on the concentrations of ligands and receptors in the above formulas, to account for example for saturation or self-limiting mechanisms whereby production of a receptor or ligand stops if they are already available in sufficient amounts (see below the production rates of RANKL and OPG). A closed form or analytical solution of the system of equations (5) is not always possible and a numerical solution is often needed. For example, substituting the first equation into the second in (5) leads to a system of equations involving only the ligand concentrations L. If there is only one ligand type, this reduces to a single polynomial equation in L whose order is equal to the number of receptors binding that ligand plus one.

Intracellular. signalling and cell response

The raw "input signal" communicated to a cell by extracellular signalling depends on the number of receptors on the cell that are bound to a ligand (occupied receptors). It is convenient to take as input signal the fraction ${\nu }_{RL}$ of occupied receptors on the cell, equal to

${\nu }_{RL}=\frac{\stackrel{⌢}{RL}}{R_{tot}}=\frac{\stackrel{⌢}{RL}}{R+\stackrel{⌢}{RL}}=\frac{L}{k_{RL}+L},$E6

where R_tot is the concentration of both free and bound receptors. The fraction of occupied receptors on a cell only depends on the concentration of extracellular free ligand L and a single binding parameter k_RL. In absence of bound complex degradation (i.e. if D_RL=0), k_RL is the so-called dissociation binding constant (Lauffenburger and Linderman 1996).

Cells may behave in several different ways in response to extracellular signalling molecules. A signalling molecule may either "activate" or "repress" a certain cell response, such as differentiation, proliferation, apoptosis, or the expression of other signalling molecules and receptors (see e.g. Tables 1 and 2). In a mathematical model, these cell responses are associated with certain model parameters. An effective way of modelling the response of a cell to a signalling molecule is to specify a phenomenological relation between the strength of the input signal to that cell (its fraction of occupied receptors) and the model parameter representing the cell response. This relation integrates potentially complex intracellular signalling mechanisms that are not explicitly modelled, and so is of phenomenological nature. In fact, such a relation is prone to experimental determination, and has been measured experimentally in other contexts.

For example, in the context of human fibroblasts stimulated by epidermal growth factor (EGF), the mitogenic response of the fibroblasts has been shown to depend linearly on the fraction of occupied EGF receptors Lauffenburger, D. A. and J. J. Linderman (1996). Receptors: Models for Binding, Trafficking, and Signaling Oxford, Oxford University Press., Fig. 6-7, p.249.

In the following, "activation" or "repression" of a cell response by a ligand will always be represented by modulating the model parameter associated with that cell response with a function of the fraction of occupied receptors. To this effect, we introduce two classes of functions of${\nu }_{RL}$: activator functions ${\pi }_{\alpha }^{\text{act}}({\nu }_{RL})$ and repressor functions${\pi }_{\alpha }^{\text{rep}}({\nu }_{RL})$, where $\alpha$is a shape parameter (see Figure 4). These activator and repressor functions are defined as follows:

$${\pi }_{\text{act}}^{\alpha }(\nu )=\frac{\nu }{\alpha +\nu (1-\alpha )},{\pi }_{\text{rep}}^{\alpha }(\nu )=1-{\pi }_{\text{act}}^{\alpha }(\nu )=\frac{\alpha (1-\nu )}{\nu +\alpha (1-\nu )},\alpha >0.$$E7

For example, if D is a model parameter representing a rate of cell differentiation, its "activation" by a signalling molecule L will be represented by$D{\pi }_{\alpha }^{\text{act}}({\nu }_{RL})$, whereas its "repression" by L will be represented by$D{\pi }_{\alpha }^{\text{rep}}({\nu }_{RL})$. Since receptor occupancy only depends on the free ligand concentration, one can also directly write the overall relation between free ligand concentration and the strength of the cell response from Eqs. (6) and (7):

${\pi }_{\text{act}}^L\equiv {\pi }_{\text{act}}^{\alpha }\left({\nu }_{RL}(L)\right)=\frac{L}{{\kappa }_{RL}+L},{\pi }_{\text{rep}}^L\equiv {\pi }_{\text{rep}}^{\alpha }\left({\nu }_{RL}(L)\right)=\frac{1}{{\kappa }_{RL}+L},\text{where }{\kappa }_{RL}=\alpha k_{RL}$E8

These functions have the same form as the fraction of occupied receptors ${\nu }_{RL}$ and the fraction of unoccupied receptors $1-{\nu }_{RL}$ (see Eq. (6), except that the parameter ${\kappa }_{RL}$ corresponds to a rescaling of the binding parameter $k_{RL}$ by the shape parameter$\alpha$. In this way, a cell is able to respond to a same receptor stimulus with different potencies in different behaviours. In the bone cell population model presented below, we will use the activator and repressor functions in the form (8). Table 3 in the appendix lists both the receptor-ligand binding parameters k_RL and the parameters κ_RL involved in the phenomenological relations between receptor occupancy and cell response.

3.2. Bone cell population model

In this section we give a summary of the bone cell population models developed by our group (Pivonka, Zimak et al. 2008; Pivonka, Zimak et al. 2010; Buenzli, Pivonka et al. 2012). The presented model takes into account the RANK-RANKL-OPG signalling pathway between bone cells, the action of TGF-β on different types of bone cells and pre-osteoblast proliferation. A revised formulation of competitive binding reactions in the RANK-RANKL-OPG system is presented, which is used to also include the effect of the anti-resorptive drug denosumab. A schematic picture of the model is presented in Figure 5.

Changes. in porosity and bone matrix fraction due to cell activity

The activity of osteoclasts and osteoblasts leads to the removal and deposition of new bone. This activity modifies the volume fraction of bone matrix in the bone tissue (f_bm=V_bm/V_RVE), or equivalently the bone porosity (f_vas=V_pores/V_RVE=1-f_bm). Osteoblasts deposit osteoid, a collagen-rich substance which later mineralizes into new bone. Primary mineralisation of osteoid is fast, i.e., 70% of the maximum mineral density is reached within a few days in humans (Parfitt 1983). Given the much larger time spans involved in evolution of bone diseases, we assume that osteoblasts “instantaneously” deposit “fully” mineralised new bone matrix as a first approximation. We further assume that the resorption rate of bone matrix k_res by an individual active osteoclast (in volume per unit time) and the formation rate of new bone matrix deposition k_form by an individual osteoblast (in volume per unit time) are constant. The evolution of the bone matrix volume fraction is thus given by

${\partial }_t{f}_{\text{bm}}\text{= }{k}_{\text{form}}\cdot {\text{OB}}_{\text{a}}-k_{\text{res}}\cdot {\text{OC}}_{\text{a}}$E9

where OB_a and OC_a denote the density of active osteoblasts and active osteoclasts respectively (number of cells per unit volume). These densities are determined by biochemical and cellular processes such as differentiation from precursor cell types and apoptosis under the control of several signalling molecules. Such processes are governed in general by the material balance equation expressed for each cell type. The equations governing the evolution of the bone cell densities are presented in the following.

Governing. equations for bone cell densities

From the osteoblastic lineage the following cell types are taken into account: uncommitted osteoblasts (OB_u) (representing a constant pool of mesenchymal stem cells or bone marrow stromal cells), osteoblast precursor cells (OB_p), and active osteoblasts (OB_a). From the osteoclastic lineage the following cell types are taken into account: osteoclast precursor cells (OC_p), and active osteoclasts (OC_a). The material balance of each cell type (in the representative volume element under consideration) is determined by specifying their production rate (source term) and elimination rate (sink term). The production and elimination of cells considered here are schematically represented in Figure 5 as flows between the different cell types. Cell differentiation accounts both for an elimination (outflux) from the pool of precursor cells and for a production (influx) into the pool of differentiated cells. In our model, OB_us differentiate into OB_ps in presence of TGF-β, but the differentiation of OB_ps in to OB_as is repressed by TGF-β. Osteoclasts become active upon RANKL binding to their receptor RANK. The availability of RANKL is in turn determined by the concentration of OPG and of the systemic hormone PTH (see below). The material balance equations for the densities of OB_p, OB_a and OC_a thus take the form of a system of rate equations with first-order production and elimination rates (i.e. rates proportional to the densities of cells) to account for the population size:

$${\partial }_t{\text{OB}}_{\text{p}}=D_{{\text{OB}}_{\text{u}}}\cdot {\text{OB}}_{\text{u}}\cdot {\pi }_{{\text{act,OB}}_{\text{u}}}^{\text{TGF-β}}+{{\text{OB}}_{\text{p}}}_{}\cdot {\text{OB}}_{\text{p}}-D_{{\text{OB}}_{\text{p}}}\cdot {\text{OB}}_{\text{p}}\cdot {\pi }_{{\text{rep,OB}}_{\text{p}}}^{\text{TGF-β}}$$E10

$${\partial }_t{\text{OB}}_{\text{a}}=D_{{\text{OB}}_{\text{p}}}\cdot {\text{OB}}_{\text{p}}\cdot {\pi }_{{\text{rep,OB}}_{\text{p}}}^{\text{TGF-β}}-A_{{\text{OB}}_{\text{a}}}\cdot {\text{OB}}_{\text{a}}$$E11

$${\partial }_t{\text{OC}}_{\text{a}}=D_{{\text{OC}}_{\text{p}}}\cdot {\text{OC}}_{\text{p}}\cdot {\pi }_{{\text{act,OC}}_{\text{p}}}^{\text{RANKL}}-A_{{\text{OC}}_{\text{a}}}\cdot {\text{OC}}_{\text{a}}\cdot {\pi }_{{\text{act,OC}}_{\text{a}}}^{\text{TGF-β}}$$E12

where D_OBu, D_OBp and D_OCp denote the differentiation rates of uncommitted osteoblast progenitors, osteoblast precursor cells, uncommitted osteoclast progenitors and osteoclast precursor cells, respectively. ${{\text{OB}}_{\text{p}}}_{}$is the proliferation rate of osteoblast precursor cells, assumed to be self-limited by the current population of OB_ps: ${{\text{OB}}_{\text{p}}}_{}=P_{{\text{OB}}_{\text{p}}}(1-{\text{OB}}_{\text{p}}/{\text{OB}}_{\text{p}}^{\text{sat}})$(Buenzli, Pivonka et al. 2012). A_OBa and A_OCa denote the apoptosis rates of active osteoblasts and active osteoclasts. Note that OB_us and OC_ps are assumed to be constant and so are not part of the state variables. In the above equations, the signalling molecules TGF- β and RANKL influence cell differentiation and apoptosis through the activator and repressor regulatory functions introduced in Eq. (8). In particular, ${\pi }_{{\text{act,OB}}_{\text{u}}}^{\text{TGF-β}},{\pi }_{{\text{rep,OB}}_{\text{p}}}^{\text{TGF-β}},$and $${\pi }_{{\text{act,OC}}_{\text{a}}}^{\text{TGF-β}}$$ represent an activation or repression by TGF-β of osteoblast differentiation and osteoclast apoptosis, whereas$${\pi }_{{\text{act,OC}}_{\text{p}}}^{\text{RANKL}}$$ represents the activation by RANKL of osteoclast differentiation.

Figure 5.

Mathematical model of bone cell interactions taking into account different bone cell types and biochemical regulatory mechanism including the RANK-RANKL-OPG signaling pathway, proliferation of osteoblast precursor cells and action of TGF-β on bone cells.

Binding. reactions in the RANK-RANKL-OPG system

We now exemplify the general approach presented earlier for modelling competitive binding reactions between receptors and ligands in the RANK–RANKL–OPG system, including also the anti-resorptive drug denosumab. To simplify, we only consider one type of RANKL, expressed on the membrane of OB_ps. However, we consider that RANKL may bind with three kinds of receptors, i.e. RANK (bound to the membrane of OC_ps), OPG and denosumab (soluble molecules). The equations (5) governing the steady-state concentration of RANKL, RANK, OPG and denosumab are fully determined once the production rate P_RANKL of RANKL, as well as the production rate of its receptors, P_R, R=RANK, OPG, denosumab are specified. All degradation rates (D_RANKL and D_R) are assumed constant (see Table 3 in the appendix). The formulation of the equations governing the binding reactions in the RANK-RANKL-OPG system presented below differs from the formulation used in (Lemaire, Tobin et al. 2004) and (Pivonka, Zimak et al. 2008; Pivonka, Zimak et al. 2010). In these previous works, competitive binding was only partially implemented.

For example, while the concentration of RANKL was dependent on that of OPG, the concentration of OPG was not dependent on that of RANKL.

A fully consistent approach is considered here instead.

The expression of RANKL on OB_p is known to be increased by the systemic hormone PTH (see Table 1). Accordingly, we assume that PTH is activating the production rate ${\beta }_{\text{OBp}}^{\text{RANKL}}$of RANKL by OB_ps.

In Lemaire et al. 2004 and Pivonka et al. 2008 regulation by PTH was assumed to limit the so-called 'carrying capacity' of RANKL on a cell instead, but there is no biological evidence for this fact. A more direct regulatory effect of PTH is assumed here.

To account for a limited space per cell where RANKL can be expressed (i.e. a limited carrying capacity), the production rate of RANKL is also multiplied by a saturation function such that production stops when the number $N_{{\text{OB}}_{\text{p}}}^{\text{RANKL}}$of RANKL molecules on the OB_p reaches a maximum number,$N_{{\text{OB}}_{\text{p}},\max }^{\text{RANKL}}$. Finally, to account for the population size of OB_p expressing RANKL, the production rate of RANKL is multiplied by the density of OB_p cells, such that:

$$P_{\text{RANKL}}={\beta }_{\text{OBp}}^{\text{RANKL}}{\pi }_{\text{act}}^{\text{PTH}}\left(1-\frac{N_{\text{OBp}}^{\text{RANKL}}}{N_{\text{OBp,max}}^{\text{RANKL}}}\right){\text{OB}}_{\text{p}}$$E13

where the total number of RANKL ligands (free and bound) present on the membrane of an OB_p cell, $N_{\text{OBp}}^{\text{RANKL}}$, is given by

$$N_{\text{OBp}}^{\text{RANKL}}=\frac{{\text{RANKL}}^{\text{tot}}}{{\text{OB}}_{\text{p}}},{\text{where RANKL}}^{\text{tot}}=\text{RANKL}+{\displaystyle {\sum }_R\stackrel{⌢}{R\text{RANKL}}}=\text{RANKL}\left(1+{\displaystyle {\sum }_R\frac{R}{k_{R,RANKL}}}\right)$$E14

The expression of OPG by OB_as is known to be inhibited by PTH (see Table 1). Accordingly, we assume that PTH is repressing the production rate of OPG per OBa cell,${\beta }_{\text{OBa}}^{\text{OPG}}$. Following (Pivonka, Zimak et al. 2008), the endogeneous production of OPG is also assumed to be self-limited, and so the production rate of OPG is multiplied by a saturation function, such that production stops when the concentration of OPG reaches a maximum concentration OPG_max. Finally, to account for the population size of OB_a expressing OPG, the production rate of OPG is multiplied by the density of OB_a cells. Hence:

$$P_{\text{OPG}}={\beta }_{\text{OBa}}^{\text{OPG}}{\pi }_{\text{rep}}^{\text{PTH}}\left(1-\frac{\text{OPG}}{{\text{OPG}}_{\text{max}}}\right){\text{OB}}_{\text{a}}$$E15

Following Refs (Lemaire, Tobin et al. 2004; Pivonka, Zimak et al. 2008), the receptor RANK is assumed be present on OC_ps in constant numbers. Thus the production rate of RANK is proportional to that of OC_p cells:

$$P_{\text{RANK}}=N_{\text{OCp}}^{\text{RANK}}{P}_{\text{OCp}}^{}$$E16

where $N_{\text{OCp}}^{\text{RANK}}$is the number of RANK receptors per OC_p cell and P_OCp is the production rate of OC_ps. The degradation rate of RANK receptors (whether free or bound) corresponds here to the differentiation rate of OC_ps into OC_as, D_OCp (in the model, active osteoclasts are not assumed to express RANK receptors for simplicity). Therefore, one has $D_{\text{RANK}}=D_{\text{RANK,RANKL}}=D_{{\text{OC}}_{\text{p}}}$in the first equation in (5) such that only the ratio P_RANK/D_RANK occurs in that equation. With Eq. (14), this ratio is equal to:

$$\frac{P_{\text{RANK}}}{D_{\text{RANK}}}=N_{{\text{OC}}_{\text{p}}}^{\text{RANK}}{\text{OC}}_{\text{p}}$$E17

since the density of OC_ps is equal to P_OCp/D_OCp.

The production rate of denosumab $P_{\text{denosumab}}^{\text{dosing}}$ will be set as given function of time according to the chosen dosing regime of the drug. All the binding and degradation properties of denosumab are assumed identical to those of OPG. In summary, the equations governing the concentrations of RANKL, RANK, OPG and denosumab are:

$$\begin{array}{l}\text{RANKL}=\frac{P_{\text{RANKL}}}{D_{\text{RANKL}}+{\displaystyle \sum _{R=\text{RANK,OPG,denosumab}}\frac{D_{R,RANKL}}{k_{R,RANKL}}}R}\\ \text{RANK}=\frac{N_{\text{OCp}}^{\text{RANK}}{\text{OC}}_{\text{p}}}{1+\frac{1}{k_{\text{RANK, RANKL}}}\text{RANKL}}\\ \text{OPG}=\frac{P_{\text{OPG}}}{D_{\text{OPG}}+\text{}\frac{D_{\text{OPG, RANKL}}}{k_{\text{OPG, RANKL}}}\text{RANKL}}\\ \text{denosumab}=\frac{P_{\text{denosumab}}^{\text{dosing}}}{D_{\text{denosumab}}}\end{array}$$E18

where P_RANKL and P_OPG are given by Eqs. (13)–(15) and all the parameters values are listed in Table 3 in the appendix. A closed form solution for these equations is not possible. Equations (10)–(12) governing the evolution of the cell populations and Eqs (18) form a system of differential algebraical equations (DAE) that need to be solved numerically.

4. Numerical simulations of bone diseases and therapeutic interventions

In this section we outline how bone diseases and drug treatments can be investigated in a computational modelling framework such as the bone cell population model of bone remodelling presented in Section 3. We refer to these computational experiments as “in-silico disease modelling” and “in-silico drug treatment modelling”. We will then exemplify such in-silico experiments in our bone cell population model by considering a simple case of an osteoporotic condition, modelled by a decrease in the endogeneous OPG production rate. While osteoporosis is often associated with dysregulations in several pathways or components, such as increased osteocyte apoptosis due to reduced estrogen levels in post-menopausal osteoporosis, a decrease in endogeneous OPG production in our model leads to a high turnover rate with a catabolic bone imbalance, which is characterisitic of the first stage of age-related osteoporosis. Simple treatments of this disease will be first simulated by prescribing constant changes in model parameters in an attempt to restore bone volume (Section 4.1). The discussion of these simple in-silico treatments is largely based on the work by (Pivonka, Zimak et al. 2010). Finally, we will investigate a temporal in-silico treatment mimicking the anti-resorptive action of the drug denosumab with different dosing regimes (Section 4.2).

To model a disease in-silico, it is important that the computational model contains suitable parameters associated with the known pathophysiology of the disease. This emphasises the importance of including a comprehensive set of signalling molecules between bone cells when modelling osteoporosis and other disorders of bone remodelling. As has been pointed out in Section 2, many catabolic bone diseases are associated with disruptions in the RANK-RANKL-OPG pathway. However, several model parameters are involved in this pathway and so several model parameters may be modified to model such a disease. While this complicates the modelling of in-silico diseases, it also opens up the possibility for patient-specific disease modelling. The measurements of different bone properties from a patient, such as the temporal evolution of the bone matrix f_bm from micro-CT scans, bone turnover rates from histological estimates of bone cell numbers in biopsies or from serum levels of turnover markers, can be used to estimate the evolution of model parameters required to simulate the evolution of the disease in this particular patient (e.g. using optimisation algorithms). Often, however, the time evolution of such properties cannot be known with precision. In a situation where for example only two time points are known in a single bone property of a patient, one may assume a constant change in a single model parameter by lack of further information. This is the situation considered in the model of osteoporosis simulated below, where it is assumed that the bone matrix volume fraction is known only at the onset of the disease (t₀) and at one further time point after disease progression and before treatment (t₁) (see Figure 6). An in-silico therapy is then applied in a second step starting from the time point t₁ to treat the disease. The effect of this therapy is checked at a later time point (t₂).

To model realistic in-silico therapies, is it crucial to identify the specific action of a therapeutic drug on the various biochemical or cellular components of the system, and to know the pharmacokinetic properties of that drug, such as the rapidity of its clearance from the system. Similarly to a disease, a drug may affect a single component or multiple components in a biological system. Numerically simulating the effect of a drug can help understand the overall action of the drug, especially when responses of different characteristic times are involved in different affected components. However, simpler in-silico therapies, in which the effect of a constant change in a model parameter is studied in a sensitivity analysis, are also worthwhile to investigate. Such theoretical treatments can give insights into the most effective therapeutic strategies and potential drug targets.

The identification of successful treatment strategies is critical for the development of clinical applications of in-silico therapies. Bone therapeutic guidelines currently classify therapeutic agents into anabolic or anti-resorptive drugs based on their effect on bone volume (or bone mass) and bone turnover (see (Riggs and Parfitt 2005) and Section 2 for more details). Using the terminology of (Riggs and Parfitt 2005) anti-resorptive drugs (such as bisphosphonates and denosomab) increase bone volume, but lead to low bone turnover, whereas anabolic drugs (such as intermittent PTH) strongly increase bone volume while producing high bone turnover. Depending on the current bone volume and bone turnover state of a particular patient, a treatment regime combining both types of drugs may thus be advised. Optimisation algorithms could help a clinician find the "optimal" treatment regime for that patient. In Figure 6, three therapeutic treatments are schematically shown to successfully restore bone volume and bone turnover rate (bone cell numbers), however with different efficiencies. The "ideal" treatment (red curve) recovers all of the bone lost and normal cell numbers instantaneously. Clearly this is impossible, and realistic treatments (blue and green curves) may only restore the normal state after some period of time. However, different treatments will usually exhibit different time courses for bone volume and cell numbers. These differences may be used in the search for optimal treatment regimes. For example, one may want to minimise the area between the "ideal treatment" curve and the actual treatment curve, i.e. to minimise the "therapy error" (gray area in Figure 6). Weighted combinations of several criteria, for example involving turnover rate and cell numbers, can also be used to define an objective "therapy error" to minimise.

Figure 6.

Schematic representation of catabolic bone diseases and applied therapy: (a) volume fraction of bone matrix [%] vs. time [years] and (b) normalized bone cell numbers [-] vs. time [years] (t₀ = onset of bone disease; t₁ = start treatment; t₂... end treatment; ideal therapy... red curve; ‘‘example’’ therapy... blue curve; therapy error... gray shaded area between ideal therapy and typical therapy, Figure modified from Pivonka et al 2010).

4.1. Simple in-silico treatment of osteoporosis in the bone cell population model: constant changes in parameter values

In this subsection we extend some of the results presented in (Pivonka, Zimak et al. 2010) using the improved model of bone cell interactions presented in Section 3, which includes proliferation of osteoblast precursor cells (Buenzli, Pivonka et al. 2012) and competitive receptor-ligand binding reactions in the RANK-RANKL-OPG system. In (Pivonka, Zimak et al. 2010), we showed that many bone disorders that act via the RANK-RANKL-OPG system induce more pronounced osteoclast responses than osteoblast responses. This indicates that targeting the RANK-RANKL-OPG pathway is not very effective at triggering bone formative responses. Furthermore, it was shown in this study that the severity of catabolic bone diseases strongly depended on how many components of the RANK-RANKL-OPG pathway were assumed to be affected. Changes of single parameters in the RANK–RANKL–OPG system led to less severe bone loss compared to multiple concurrent changes of parameters. In a subsequent study we have also shown that bone disorders related to excessive bone formation (such as van Buchem disease and sclerosteosis) may act through a different pathway, with the Wnt pathway being a major candidate (Buenzli, Pivonka et al. 2012). Regulation of osteoblast differentiation and proliferation by Wnt were shown to be very effective mechanisms to induce bone formative responses.

As mentioned above, modifying the value of a model parameter as an in-silico therapeutic treatment does not represent a realistic drug treatment scenario. However, evaluating the sensitivity of the bone response to such changes provides insights into the general mechanism of action of a potential therapeutic agent and can help identify the most effective pathways that should be targeted for disease intervention, even if such treatments may not currently exist. In Figure 7, we apply such constant changes in model parameters to restore the bone loss induced by a simulated osteoporotic condition. This osteoporotic condition is simulated by a reduction in endogeneous OPG production (i.e. a reduction in the value of the model parameter$${\beta }_{{\text{OB}}_{\text{a}}}^{\text{OPG}}$$) assumed to develop instantaneously at time t₀=0 from a healthy state. The evolution of the bone volume fraction f_bm is followed in a representative volume element located in cortical bone with initial bone volume fraction f_bm(t₀)=95%. The reduction in OPG production is chosen such that a 5% of bone volume fraction is lost after one year, i.e. at t₁=365 days. At time t₁, different in-silico treatments are commenced by changing single model parameters with the objective to restore bone volume over a period of four years (i.e., at time t₂=(1+4) 365=1825 days). Changes in the model parameters D_OBu, P_OBp, D_OBp, A_OBa (targeting osteoblast development), in A_OCa and in the binding parameter between RANK and RANKL, k_RANK-RANKL (targeting osteoclast development) have been investigated.

The numerical results shown in Figure 7 clearly show that model parameters associated with osteoblast development have a strong potential to restore bone volume within the prescribed treatment period of four years. In fact, the numbers of OB_as and OC_as are increased in these situations (not shown) while bone balance is positive. The combination of high bone turnover rate and positive bone balance enables a quick restoration of bone volume. According to the classification by (Riggs and Parfitt 2005), these model parameters can therefore be characterised as having a pro-anabolic potential.

On the other hand, parameters associated with osteoclast development (including parameters associated with the RANK-RANKL-OPG pathway) have a weak potential to restore bone volume. Increasing osteoclast apoptosis A_OCa (which can be associated with the action of a bisphosphonate) or reducing the binding affinity between RANK and RANKL via the parameter k_RANK-RANKL fail to restore bone volume within the prescribed treatment period of four years. In fact, the numbers of OB_as and OC_as are significantly reduced in these situations. Thus, although bone balance is positive, bone turnover rate is too low for this positive bone balance to build up bone volume fast enough. According to the classification by (Riggs and Parfitt 2005), these parameters can be characterised as having an anti-resorptive potential. The initial strong formative response is due to an immediate reduction of active osteoclasts without corresponding reduction in active osteoblasts. After a period of about 60 days, this initial reduction in active osteoclasts leads to a marked depletion of TGF-β and so to decreased osteoblastogenesis. This depletes RANKL which further enhances the initial reduction of active osteoclasts and leads to a low bone turnover state in the second part of the bone response.

Figure 7.

Single constant parameter therapies for treatment of catabolic bone disease due to constant reduction in OPG production rate leading to 5% bone loss over one year: (a) bone matrix volume fraction [%] vs. time [days] and (b) normalized cell numbers [-] vs. time [days] (D_OBu=differentiation of uncommitted osteoblasts, P_OBp=proliferation of osteoblast precursors, A_OBa=apoptosis of active osteoblasts, A_OCa=apoptosis of active osteoclasts, _{kD,RANK-ARNKL}=dissociation binding constant RANK-RANKL, D_OBp=differentiation of osteoblast precursors).

According to the definition of “therapy error” presented in Figure 6, the most effective treatment is achieved via manipulation of the differentiation rate of osteoblast precursor cells (D_OBp). However, the ability of D_OBp to restore bone volume depends strongly on the proliferative potential of pre-osteoblasts, i.e. on the ratio P_OBp/D_OBu (Buenzli, Pivonka et al. 2012). At low ratios P_OBp/D_OBu, an increase in D_OBp depletes the pool of pre-osteoblasts enough to lead to a low bone turnover state unable to restore bone volume within the four years treatment period (data not shown). At higher ratios P_OBp/D_OBu, this depletion occurs at higher values of D_OBp and bone volume can still be restored within the treatment period as shown in Figure 7. The results of Figure 7 suggest that an alternative efficient treatment minimising therapy error in the evolution of bone volume may consist of first using an anti-resorptive treatment during the first few months, followed by a pro-anabolic treatment.

4.2. In-silico treatment of osteoporosis by denosumab in the bone cell population model

Finally, we investigate a therapeutic strategy that mimicks the effect of the anti-resorptive drug denosumab on the simulated osteoporotic condition of Section 4.1. Although a realistic model of denosumab treatment should take into account pharmacokinetic properties, here, we simply assume that denosumab is produced in the bone compartment according to its dosing regime and that it is eliminated at a constant rate (see Eq (18)). The explicit consideration of denosumab as a new biochemical component of the system associated with its own dynamics enables to investigate temporal dosing regimes consistently. Moreover, the explicit consideration of denosumab enables competitive binding effects between denosumab, RANK and OPG to bind RANKL to be fully accounted for.

Three different dosing regimes of denosumab are shown in Figure 8, all commencing at time t₁=365 days. In Figure 8a and 8b, a continuous infusion is assumed, i.e. the production rate P_denosumab of denosumab is constant. The influence of the infusion rate ("dose", in concentration per unit time) on the evolution of the bone volume fraction is shown in Figure 8a. The influence of the infusion rate for the bone volume fraction reached after four years of treatment is shown in Figure 8b. As one may expect, the lowest infusion rate (3 pM/day) has the least effect to restore bone loss. Interestingly, a medium infusion rate (7 pM/day) is more efficient to restore bone loss than lower or higher infusion rates, an effect that would be difficult to predict without computational modelling. The similarity of the evolution of the bone volume fraction in Figure 8a and those seen in Figure 7 when modifying anti-resorptive parameters is due to the fact that denosumab acts via the RANK-RANKL-OPG pathway.

In Figure 8c, a single injection of denosumab is simulated. This situation is modelled by specifying the production rate P_denosumab(t) as a function having a single peak at time t₁=365 days. The height of the peak is taken as a measure of the dose (in concentration per unit time) and the width of the peak at half height is taken to be two days. The dosing rate is seen to influence both the maximum bone volume fraction attained and the final gain in bone volume compared to the untreated case (dotted line). Higher denosumab doses lead to larger bone gain and reduce bone loss for a longer period, although a saturation of the dose dependence is quickly reached from about 100 pM/day upwards.

Finally, dosing regimes made of multiple injections spaced by regular time intervals are shown in Figure 8d. Each injection of these dosing regimes is simulated by the case 100 pM/day of Figure 8c. The dosing regime with a small administration interval (every 5 days) is comparable to constant denosumab infusion. Administration of denosumab every 10 days still has the potential to inhibit significant bone loss. However, an administration interval of 30 days at the same dose is not sufficient to efficiently inhibit bone loss.

Figure 8.

Single time dependent parameter therapy mimicking pharmacokinetic action of denosumab for treatment of catabolic bone disease due to constant reduction in OPG production rate leading to 5% bone loss over one year: (a) bone matrix f_bm [%] vs. time [days] for constant denosumab doses, (b) bone matrix f_bm [%] vs. denosumab infusion rate [pM/day], (c) bone matrix f_bm [%] vs. time [days] for singe injections of denosumab and (d) bone matrix f_bm [%] vs. time [days] for multiple injections of denosumab at a dose of 100 pM/day.

5. Summary and conclusions

In this contribution we presented a mathematical model describing bone remodelling taking into account complex bone cell interactions. A general approach to modelling cell-cell signalling and cell behaviour was proposed. This approach considers explicit extracellular signalling mechanisms and receptor-ligand binding reactions governed by the law of mass action and material balance, while intracellular signalling mechanisms are assumed to be captured by phenomenological response functions. We highlighted the fact that different patients may exhibit different bone disease patterns which can be captured in a so-called bone disease signature diagram. Within the presented framework patient specificity can be accounted for via adjustments of model parameters.

As an application of the model, we simulated an osteoporotic condition associated with a disruption in the RANK-RANKL-OPG pathway and considered various 'in-silico' therapeutic treatments. The notion of "therapy error" was introduced to identify the most efficient ways to restore normal bone mass and bone turnover from a diseased condition. In these therapeutic interventions, we could identify two groups of model parameters: (i) anti-resorptive parameters characterised by low bone turnover and (ii) pro-anabolic parameters characterised by high bone turnover. Anti-resorptive parameters corresponded to parameters of the RANK-RANKL-OPG pathway and parameters associated to osteoclast developement. Pro-anabolic parameters corresponded to parameters associated with osteoblast development.

Finally, temporal treatment regimes were investigated using a simple model mimicking the action of the drug denosumab. Competitive binding effects between denosumab, RANK and OPG to bind RANKL could be fully accounted for by the specification of explicit receptor-ligand binding reactions in the RANK-RANKL-OPG system. For constant administration regimes, there is an intermediate optimal dose leading to maximum bone gain. For the single and multiple injection regimes bone gain was stronger with higher doses or higher administration frequencies, however, the dose and administration frequency dependences of the bone response were seen to saturate quickly.

Future advancement of bone biology research will strongly rely on how well experimental and theoretical groups are able to communicate and collaborate with each other. A prerequisite for such collaborations is the mutual understanding of research methodologies employed. We hope that this chapter will provide some guidance on how theoretical tools such as mathematical modeling can be used in biomedical research in general, and particularly in bone biology research.

Appendix

Table 3.

Parameters of the model

Acknowledgement

The authors are grateful to Prof T.J. Martin (St. Vincent’s Institute of Medical Research, University of Melbourne) for his feedback in preparing this manuscript and to Dr S. Scheiner (Vienna University of Technology) for fruitful discussions. The work of A/Prof. Pivonka is supported by the Australian Research Council (ARC) for discovery project funding (DP0988427).