Feasible Simulation of Diseases Related to Bone Remodelling and of Their Treatment

1. Introduction

A great advance in understanding mechanosensing and mechanotransduction of bone tissue has occurred in the past years. However, a clear answer is yet to come. There is plentiful evidence that most cells in human body are able to sense their mechanical environment, including osteoblasts and osteocytes. Li et al. found that marrow stromal cells change their proliferation rate and gene expression patterns in response to mechanical stimulation (Li et al, 2004). Ehrlich and Lanyon mention that osteocytes produce significantly higher levels of PGE₂ (prostaglandin E₂) and PGI₂ (prostacylin) than osteoblasts and in vivo inhibition of prostaglandins prevents bone adaptation in response to mechanical strains (Ehrlich and Lanyon, 2002). Further, nitric oxide NO is also a mediator of mechanically induced bone formation. It is released, similarly, as prostanoids, in higher levels after exposure to physiological levels of mechanical loading. Moreover, Rubin discovered that dynamic loading down-regulates osteoclastic formation (Rubin et al, 2000) (more precisely, the strained bone cell downregulates its expression of RANKL (Rubin et al, 2004)), which suggests that mechanical forces play an important role in bone adaptation process.

Candidate mechanoreceptors within a cell are stretch-activated channels, integrins (membrane spanning proteins that couple the cell to its extracellular environment), connexins (membrane spanning proteins that form channels that allow the direct exchange of small molecules with adjacent cells - including intercellular communication via gap junctions), and membrane structure (Rubin, 2006). Nowadays, there is also a completely different possible explanation for transformation mechanical signal into a biochemical function. Valle et al. explain in their review that only a proper mechanical loading may lead to exposure of binding sites and thus enabling further biochemical processes to proceed (Valle, 2007). Most probably multiple mechanosensors are involved in receiving mechanical signals. Moreover, there are other studies showing the importance of many other mechanisms beside those mentioned above (Lemaire et al, 2004).

On tissue level, it is clear that mechanical loading is important in the bone remodelling process. Heřt described an interaction between the mechanical stimulation of local cells and the bone adaptation process in the 1970s (Heřt et al, 1972), Frost observed the same behaviour in his clinical praxis and summed it up in his „Utah paradigm“ (Frost, 2004). More recently, the fundamental importance of dynamic loading was accepted. Comparison of the static versus the dynamic loading effects on bone remodelling is given in a nice and inspectional review by Ehrlich and Lanyon (2002). Further, Robling provides experimental results that confirm the essential importance of the dynamic loading (Robling, 2006). It is worth mentioning that Heřt referred to this fact in his observations more than 35 years ago (Heřt et al, 1972).

From the phenomenological relations for the rate of chemical reaction based on classical irreversible thermodynamics including coupling with mechanical processes, it was shown that, although tissues are exposed to all variety of mechanical factors: straining, shear, pressure, and even dynamic electric fields, the volume variation rate is the most important mechanical stimulus driving the processes in them (Klika, 2010). However, we believe that shear rate might be important for triggering the process. (for details see (Klika, 2010; Klika and Maršík, 2009)). Thus, in the presented manuscript, the mechanical stimulation of bone remodelling was assumed to be proportional to rate of volume change. As it will be seen from the features and results of the model, it seems to capture the response of bone to changes of its mechanical environment on tissue level. Probably, the other ways of mechanosensing are controlling the triggering of bone remodelling process in a given loci. We realize that biochemical reactions are initiated and influenced primarily by genetic effects and then by external biomechanical effects (stress changes). The aim of the presented thermodynamic model is to combine biological and biomechanical factors whereas currently available models of bone remodelling focus only on one of these factors which is actually the reason why this model was developed. It should provide an estimate of the effects of increased physical activity on quality of bone even in several disease states. Such a model may also reflect changes in remodelling behaviour resulting from pathological changes to the bone metabolism or from hip joint replacement and also may help for better assessment of the risk of osteoporosis-related fractures (Lindsay, 2003). Preliminary version of the mentioned approach has been published by our team in the past but not until this considerably improved version was the model applicable to praxis (both qualitative and quantitative results) (Klika and Maršík, 2010; Klika et al., 2010) which is here being extended.

1.1. Available models of bone remodeling

With the development of computer-aided strategies and based on the knowledge of bone geometry, applied forces, and elastic properties of the tissue, it may be possible to calculate the mechanical stress transfer inside the bone (Finite Elements analysis or FE analysis). The change of stresses is followed by a change in internal bone density distribution. This allows to formulate mathematical models that can be used to study functional adaptation quantitatively and furthermore, to create the bone density distribution patterns (Beaupré et al, 1990; Carter, 1987). Such mathematical models have been built in the past. Since they calculate just mechanical transmission inside the bone and not considering cell-biologic factors of bone physiology, they just partially correspond to the reality seen in living organisms. Basically, there are essentially two groups of models for bone remodelling. One assumes that the mechanical loading is the dominant effect, almost to the exclusion of other factors, and treatment of biochemical effects are included in parameter with no physical interpretation (e.g. Beaupré et al, 1990; Carter, 1987; Huiskes et al, 1987; Ruimerman et al, 2005; Turner et al, 1997). The results or predictions of these models yield the correct density distribution patterns in physiological cases. However, they have a limited ability to simulate disease. The second group, the biochemical models, consider control mechanisms of bone adaptation in great detail, but with limited possibilities for including mechanical effects that are known to be essential (Komarova et al, 2003; Lemaire et al, 2004). We realize that biochemical reactions are initiated and influenced primarily by genetic effects and then by external biomechanical effects (stress changes). Our thermodynamic model enables to combine biological and biomechanical factors (Klika and Maršík, 2010). Such a model may also reflect changes in remodelling behaviour resulting from pathological changes to the bone metabolism or from hip joint replacement. However, it is a model and thus it is a great simplification of the complex process of bone remodelling. In this chapter, a more detailed description of biochemical control mechanisms will be added to the mentioned model (Klika and Maršík, 2010; Klika et al, 2010) which in turn leads to possibility to study several concrete bone related diseases using this model.

1.2. Simulation of diseases and of their treatment

In our previous work, the influence of mechanical stimulation on (chemical) interactions in general was studied and it was shown how to comprise this effect into a model of studied biochemical processes (Klika and Maršík, 2009; Klika, 2010). These findings were used to describe the bone remodelling phenomenon (Klika and Maršík, 2010; Klika et al, 2010; Maršík et al, 2010). In this chapter, an extension of the mentioned bone remodelling model (influences of concrete biochemical factors) will be presented where the essential significance of dynamic loading will still be apparent.

Firstly, fundamental control factors will be mentioned. As was mentioned in the introduction, the RANKL-RANK-OPG pathway is essential in the bone remodelling control. Osteoprotegerin (OPG) inhibits binding of ligand RANKL to receptor RANK and thus prevents osteoclastogenesis. Since osteoclasts are the only resorbing agents in bone, osteoprotegerin “protects bone” (osteo-protege). Further, one of the major problems connected to bone remodelling is a rapid bone loss after menopause that affects a significant portion of women after 50 years of age. Menopause is linked to a rapid decrease in estrogen levels. And because estrogen significantly affects bone density, it would be beneficial to be able to simulate the influence of estrogen levels on the bone remodelling process. Similarly, the parathyriod hormone PTH, tumour growth factor TGF-β1, and nitric oxide NO play a significant role during the bone adaptation process.

PTH causes a release of calcium from the bone matrix and induces MNOC differentiation from precursor cells, estrogen has complex effects with final outcome in decreasing bone resorption by MNOC, calcitonin decreases levels of blood calcium by inhibiting MNOC function, and osteocalcin inhibits mineralisation (Sikavitsas et al, 2001). The discovery of the RANKL-RANK-OPG pathway enabled a more detailed study of the control mechanisms of bone remodelling. Robling et al. states that all PTH, PGE (prostaglandin), IL (interleukin), and vitamin D are “translated” by corresponding cells (osteoblasts) into RANKL levels (Robling et al, 2006). Further, nitric oxide NO is known to be a strong inhibitor of bone resorption and recently it has been known that it works in part by suppressing the expression of RANKL and, moreover, by promoting the expression of OPG (Robling et al, 2006). Both these effects eventually lead to a decrease of numbers of active osteoclasts MNOC, which in turn causes decrease of bone resorption. Kong and Penninger mention that the OPG expression is induced by estrogen (Kong and Penninger, 2000). Boyle et al. add that OPG production by osteoblasts is based on anabolic stimulation from TGF-β or estrogen (Boyle et al, 2003). Martin also deals with the question how hormones and cytokines influence contact-dependent regulation of MNOC by osteoblasts. He summaries results from the carried out experiments (mainly in vitro) that PTH, IL-11, and vitamin D (1.25(OH)₂D₃ more precisely) promotes RANKL formation, which in turn increases osteoclastogenesis (Martin, 2004).

RANKL-RANK-OPG pathway mediates many of these above mention biochemical factors. Moreover, RANKL levels also reflect microcrack density. Hence, it is essential to incorporate this pathway into our model. The connection will be enabled through the amount RANKL-RANK bonds that are one of the components of developed model, noted as RR, see (Klika and Maršík, 2010).

2. Methods

Bone remodelling is a very complex mechanism. It is necessary to make reasonable simplification to identify the essential control mechanisms. The used model covers the accepted crucial features (Klika and Maršík, 2010).

In this manuscript, several schemes for treating cellular interactions (as interactions among mixtures of chemical substances) are proposed, which are commonly used in similar problems (Keener and Sneyd, 2009; Lemaire et al, 2004). In our model, individual components of bone are substrates and products of biochemical interactions, whose rates are determined by the difference in chemical potentials. The affinity and chemical potential are (in simpler cases) proportional to the logarithm of the molar concentration of the involved substances (deGroot and Mazur, 1962). Here we actually use a modified version of this law of mass action with the additional effect of dynamic loading on the interaction rates (Klika, 2010; Klika and Maršík, 2009). The law of mass action is frequently used in modelling biological processes (Keener and Sneyd, 2009). However, it is only a simplified model that can be justified by agreement between calculated data and clinical observations. From the submitted results, it can be seen that features of the presented model are consistent with clinical and experimental data.

2.1. Incorporation of RANKL-RANK-OPG pathway into the bone remodelling model

A new model for RANKL-RANK-OPG chain kinetics has been formulated and added (Klika et al., 2010) to the mentioned model of bone remodelling (Klika and Maršík, 2010). This section is being RANKL is a ligand molecule and binds to RANK forming a bond, here noted as RR and its molar concentration as$[\text{RR}]$, between osteoblasts and precursors of osteoclasts. Osteoblasts also secrete a decoy receptor osteoprotegerin OPG that binds with high affinity to RANKL and thus prevents the needed connection between osteoblasts and osteoclastic precursors.

The reaction scheme of interaction of the mentioned molecules can be described as follows:

$$\begin{array}{ll}RANKL+RANK\hfill & \underleftrightarrow{{\text{k}}_{\pm 1}}RR\hfill \\ RANKL+OPG\hfill & \underleftrightarrow{{\text{k}}_{\pm 2}}R{O}_{\text{inactive}}\hfill \end{array}$$E1

where $R{O}_{\text{inactive}}$represents the bond between the decoy OPG and ligand RANKL. Using the law of mass action (Klika and Maršík, 2009) we may infer kinetics of the above mentioned interactions. Only the simplification, when assuming a relation between forward and backward reaction rates$k_{+i}\gg k_{-i}$, is not applicable here. We get

$$\begin{array}{l}\frac{\text{d}{n}_{RANKL}}{\text{d}\tau }=-n_{RANKL}({\beta }_{\text{RK}}^{\text{RRO}}+n_{RANKL}-n_{OPG})+\hfill \\ +{\delta }_{-1}^{\text{RRO}}({\beta }_{\text{RR}}^{\text{RRO}}-n_{RANKL}+n_{OPG})-\hfill \\ -{\delta }_{+2}^{\text{RRO}}{n}_{RANKL}{n}_{OPG}+{\delta }_{-2}^{\text{RRO}}({\beta }_{\text{RO}}^{\text{RRO}}-n_{OPG}),\hfill \\ \frac{\text{d}{n}_{OPG}}{\text{d}\tau }=-{\delta }_{+2}^{\text{RRO}}{n}_{RANKL}{n}_{OPG}+{\delta }_{-2}^{\text{RRO}}({\beta }_{\text{RO}}^{\text{RRO}}-n_{OPG}),\hfill \end{array}$$E2

where

$$\begin{array}{l}{\delta }_{-1}^{\text{RRO}}=\frac{k_{-1}}{k_{+1}[{\text{RANKL}}_{\text{stand}}]},\hfill \\ {\delta }_{-2}^{\text{RRO}}=\frac{k_{-2}}{k_{+1}[{\text{RANKL}}_{\text{stand}}]},\hfill \\ {\delta }_{+2}^{\text{RRO}}=\frac{k_{+2}}{k_{+1}},\hfill \\ {\beta }_{\text{RO}}^{\text{RRO}}=\frac{C_{\text{RO}}}{[{\text{RANKL}}_{\text{stand}}]}=\frac{[{\text{RO}}_0]+[{\text{OPG}}_0]}{[{\text{RANKL}}_{\text{stand}}]},\hfill \\ {\beta }_{\text{RR}}^{\text{RRO}}=\frac{C_{\text{RR}}}{{\text{RANKL}}_{\text{stand}}}=\frac{[{\text{RR}}_0]+[{\text{RANKL}}_0]-[{\text{OPG}}_0]}{[{\text{RANKL}}_{\text{stand}}]}\text{.}\hfill \end{array}$$E3

Again $k_{\pm i}$are reaction rate coefficients, ${\delta }_i$are interaction rates, and ${\beta }_j^{\text{RRO}}$represents the normalised initial molar concentrations of corresponding substances, denoted with index 0 and finally $[{\text{RANKL}}_{\text{stand}}]$represents standard serum level of RANKL used for normalisation of molar concentrations of substance$i$,$n_i$. All the parameters have evidently a physical interpretation and are measurable. However, hardly any such in vivo data for humans is available. Fortunately, the recent progress in the understanding of bone remodelling control enabled in vitro studies of individual factors.

Quinn et al. studied the influence of RANKL and OPG concentration on a number of osteoclasts (more precisely, TRAP positive multinucleated osteoclasts) in a dose-dependent way (Quinn et al, 2001). We would like to use this data to determine the above mentioned parameters of the RANKL-RANK-OPG model. Because the carried out experiments are studying effects of RANKL and OPG separately, the reaction scheme (1) may be splitted into two separate reactions for parameter setting. This is convenient because the kinetics of a single biochemical reaction can be described using a single differential equation (in this case non-linear). Moreover, both normalised differential equations corresponding to these two reactions can be written in the same form:

$$\dot{x}=-A{x}^2-Bx+C,A>\mathrm{0,}C>0$$E4

where $A=\mathrm{1,}B={\beta }_{\text{RK}}+{\delta }_{-1},C={\delta }_{-1}{\beta }_{\text{RR}}$for the RANKL reaction and $A={\delta }_{+2},B={\delta }_{+2}{\beta }_{\text{RANKL}}+{\delta }_{-2},C={\delta }_{-2}{\beta }_{\text{RO}}$for OPG reaction. The normalised form is also useful because it decreases the number of unknown parameters. The differential equation (4) has the following solution for positive constants A, C and for initial value$x_0$:

$$\begin{array}{l}x(\tau )={\left[\frac{2A}{\sqrt{B^2+4AC}}\left(1+\frac{1+\frac{2A}{\sqrt{B^2+4AC}}\left(x_0+\frac{B}{2A}\right)}{1-\frac{2A}{\sqrt{B^2+4AC}}\left(x_0+\frac{B}{2A}\right)}{e}^{\sqrt{B^2+4AC}\tau }\right)\right]}^{-1}\cdot \hfill \\ \cdot \left[\left(1-\frac{B}{\sqrt{B^2+4AC}}\right)\frac{1+\frac{2A}{\sqrt{B^2+4AC}}\left(x_0+\frac{B}{2A}\right)}{1-\frac{2A}{\sqrt{B^2+4AC}}\left(x_0+\frac{B}{2A}\right)}{e}^{\sqrt{B^2+4AC}\tau }-\frac{B}{\sqrt{B^2+4AC}}-1\right]\text{.}\hfill \end{array}$$E5

Because we know the analytic form of function describing the kinetics of RANKL (and OPG), we may use the least square method for determination of the unknown parameters according to the carried out experiments. Data from the Quinn et al. in vitro experiment relates RANKL (and OPG) concentration to MNOC concentration (the number of osteoclasts per well). The mentioned reaction scheme (1) of RANKL-RANK-OPG interaction has an output product denoted as RR. Thus, to be able to use the mentioned data from Quinn et al., we need to relate RANKL-RANK bonds ($[\text{RR}]$) to the number of osteoclasts ($[\text{MNOC}]$). To get a precise prediction of this relationship from the presented model we would also need to know the analytical solution of the system of ODEs that describe the bone remodelling process (Klika and Maršík, 2010), which is not possible. On the other hand, the interaction that describes the relation between RANKL-RANK bonds and MNOC concentration is the first one in our bone remodelling scheme (Klika and Maršík, 2010) and it will be assumed that the number of formed and active osteoclasts is proportional to the RR concentration. It means that it was assumed that in vitro, where no remodellation occurs, the formation of osteoclasts may be described by:

$$RR\rightleftarrows MNOC\text{.}$$E6

This assumption will be used just for purposes of parameter setting and from final results it will be possible to see if this simplification was too great or not.

The next issue we have to deal with is finding a possible relation between in vitro and in vivo data. In vivo ones are more or less unavailable, especially in such a detail that is needed for parameter setting. Further, determination of standard serum levels of OPG and RANKL is needed. The problem is that in most cases in vitro concentrations have to be much higher to reach a similar effect as in vivo. Moreover, no such relation may exist. It will be assumed that there is a correspondence among these two approaches and that it is linear, i.e. in vivo data can be gained from in vitro after appropriate scaling of concentrations.

The search for standard serum levels of osteoprotegerin and RANKL was not simple. Studies differ greatly in the presented values. Kawasaki et al states that the standard level of osteoprotegerin is $250\frac{pg}{\mu l}$(Kawasaki et al, 2006) and Moschen et al. mention $800\frac{pg}{\mu l}$(Moschen et al, 2005). Further, Eghbali-Fatourechi et al. determined OPG serum levels to be $2.05\frac{pmol}{l}$(Eghbali-Fatourechi et al, 2003). The probable cause of these discrepancies lies in differently used techniques of gaining osteoprotegerin and measuring its concentration. Kawasaki et al. measured the amount of RANKL in gingival crevicular fluid, Moschen et al. performed collonic explant cultures from biopsies and consequently measured RANKL and OPG levels using an ELISA kit, and Eghbali-Fatourechi used a different cell preparation technique followed by measurement with an ELISA kit. One of the manufacturers of the ELISA kit for assessment OPG levels cites several studies on OPG levels in humans and also submits results from their own research (Elisa OPG kit, 2006). At least all these measurements are carried out by the same measurement technique and are comparable. Therefore, we set standard OPG and RANKL levels according to data that are there referred to - $[{\text{RANKL}}_{\text{stand}}]=0.84\frac{pmol}{l}=55\cdot 0.84\frac{pg}{ml}=46.2\frac{pg}{ml}$and $[{\text{OPG}}_{\text{stand}}]=1.8\frac{pmol}{l}=20\cdot 1.8\frac{pg}{ml}=36\frac{pg}{ml}$in serum (Kudlacek et al, 2003), where the knowledge of molecular weights $M{W}_{\text{RANKL}}={5510}^3,M{W}_{\text{OPG}}={2010}^3$was used (Elisa OPG kit, 2006; RANKL product data sheet, 2008). Now it is needed to find a reasonable relation with in vitro data from Quinn et al that will be used for the least squares method for parameter estimation. The following consideration will be used: the physiological range of levels of OPG and RANKL will be found and consequently related to studied effective in vitro range by Quinn. OPG serum levels found in human are $12$-$138\frac{pg}{ml}=0.6$-$6.9\frac{pmol}{l}$and RANKL serum levels are 0-$250\frac{pg}{ml}=0$-$4.55\frac{pmol}{l}$with standard values of $0.84\frac{pmol}{l}$for RANKL and $1.8\frac{pmol}{l}$for OPG, respectively. When we relate these values to the in vitro ranges of RANKL 0-$500\frac{ng}{ml}$and of OPG 0-$30\frac{ng}{ml}$, we get the in vitro equivalents for standard values:$[{\text{RANKL}}_{\text{invitrostand}}]=92.3\frac{ng}{ml}$,$[{\text{OPG}}_{\text{invitrostand}}]=7.83\frac{ng}{ml}$.

A list of parameters that will be determined by least squares from the RANKL experiment are the following:

$${\delta }_{-1}^{\text{RRO}},{\tau }_{7\text{days}}^{\text{RRO}},n_{R{K}_0},n_{R{R}_0}$$E7

where ${\tau }_{7\text{days}}^{\text{RRO}}$is the dimensionless time that corresponds to 7 days. Before the parameter setting is by curve fitting (least square method) is carried out, it is reasonable to have at least some estimation of parameter values. Because the normalisation was done by division with term $k_{+1}{[{\text{RANKL}}_{\text{stand}}]}^2$and from (3), we get:

$$\begin{array}{l}{\tau }_{7\text{days}}^{\text{RRO}}=t{k}_{+1}{C}_{\text{RR}}{|}_{t=7days}={610}^5{10}^7{10}^{-12}\approx {10}^0,\hfill \\ n_{R{R}_0}\approx {10}^0,\hfill \\ n_{R{K}_0}\approx {10}^0,\hfill \\ {\delta }_{-1^{\text{RRO}}}=\frac{k_{-1}}{k_{+1}[{\text{RANKL}}_{\text{stand}}]}\approx k_{-1}{10}^5,\hfill \end{array}$$E8

where the value of $k_{+1}$was estimated from the parameter setting in the bone remodelling model, standard value of RANKL $[{\text{RANKL}}_{\text{stand}}]\approx 1\frac{pmol}{l}$was mentioned above, and the $k_{-1}$value may be anywhere in $({\mathrm{0,10}}^7)$but most probably lower than one.

The least square method with the used data from Quinn et al. (Quinn et al, 2001) and the analytic function as described above gives the following estimates:

$$$$E9

$${\delta }_{-1}^{\text{RRO}}={4.9210}^{-6},{\tau }_{7\text{days}}^{\text{RRO}}=\mathrm{4.64,}{n}_{R{K}_0}=\mathrm{1.037,}{n}_{R{R}_0}=0.0947$$E10

If we compare these values with their order estimation above, we see that the values are acceptable and the curve fit is as well, see Fig. 1a.

Now, we may proceed with OPG parameters. The difference is that if we use only the second reaction of RANKL-RANK-OPG reaction scheme (1), we do not know how initial OPG concentration influences the number of bonds between RANKL and RANK. However, this influence is mediated by a decrease in number of available ligands RANKL by binding with OPG. Because OPG binds with higher affinity to ligand RANKL than this ligand to its receptor RANK (otherwise the decoy effects of OPG would be very limited), it will be assumed that OPG binds to RANKL more rapidly than the competiting reaction. The reason for this is again in the need of analytic solution of differential equations that govern the kinetics of mentioned processes (it was not possible to solve the full system of two differential equations (2) so the mentioned simplification was needed; again, from the results that follow it seems reasonable). Thus, the influence of levels of osteoprotegerin on the RR concentration may be mediated by an appropriate modification of initial concentration of RANKL which in turn affects the resulting RR concentration. Schematically:

$$2^{\text{nd}}reactionin(1)\to [\text{OPG}](t)$$E11

and consequently$[{\text{RANKL}}_0]=[\text{OPG}]({\tau }_{\text{OPG}})$, which is used in

$$1^{\text{st}}reactionof(1)\to [\text{RR}][{\tau }_{7\text{days}}]$$E12

were ${\tau }_{\text{OPG}}$is a time to be determined.

Figure 1.

RANKL and OPG fitted solutions (blue curves) by least squares method to data measured (dots) by Quinn et al. (Quinn et al, 2001). Firstly,${\text{n}}_{\text{RR}}$as a function of${\text{n}}_{{\text{RANKL}}_0}$is determined and consequently${\text{n}}_{\text{RR}}$as a function of${\text{n}}_{{\text{OPG}}_0}$, created by embedding dependency of$[\text{RANKL}]$on$[\text{OPG}]$and of$[\text{RR}]$on$[\text{RANKL}]$concentration, was found.

The already determined parameters from the RANKL setting will be used and only the yet unknown will be determined, i.e.

$${\delta }_{-2}^{\text{RRO}},{\delta }_{+2}^{\text{RRO}},{\tau }_{\text{OPG}}^{\text{RRO}},n_{R{O}_0}$$E13

Again, the least squares in the case of OPG give the following estimates (based on data from Quinn and the fact that molecular weight of RANKL is $55*{10}^3$and of OPG$20*{10}^3$):

$${\delta }_{-2}^{\text{RRO}}={5.8610}^{-19},{\delta }_{+2}^{\text{RRO}}=\mathrm{12.96,}{\tau }_{\text{OPG}}^{\text{RRO}}=\mathrm{11.36,}{n}_{R{O}_0}=\mathrm{6.135.}$$E14

Also, the values are admissible and the curve fit is reasonable (the function here is much more complex because OPG concentration is firstly used to determine an initial RANKL concentration for a consecutive reaction that finally gives $[\text{RR}]$outcome), see Fig 1b.

If the mentioned results of parameter estimation are combined, all the needed values of parameters of RANKL-RANK-OPG model (3) may be inferred:

$$\begin{array}{l}{\delta }_{-1}^{\text{RRO}}=\frac{k_{-1}}{k_{+1}[{\text{RANKL}}_{\text{stand}}]}={4.9210}^{-6},\hfill \\ {\delta }_{-2}^{\text{RRO}}=\frac{k_{-2}}{k_{+1}[{\text{RANKL}}_{\text{stand}}]}={5.8610}^{-19},\hfill \\ {\delta }_{+2}^{\text{RRO}}=\frac{k_{+2}}{k_{+1}}=\mathrm{12.96,}\hfill \\ {\beta }_{\text{RO}}^{\text{RRO}}=\frac{C_{\text{RO}}}{[{\text{RANKL}}_{\text{stand}}]}=\frac{[{\text{RO}}_0]+[{\text{OPG}}_0]}{[{\text{RANKL}}_{\text{stand}}]}=6.135+n_{OP{G}_0},\hfill \\ {\beta }_{\text{RR}}^{\text{RRO}}=\frac{C_{\text{RR}}}{[{\text{RANKL}}_{\text{stand}}]}=\frac{[{\text{RR}}_0]+[{\text{RANKL}}_0]-[{\text{OPG}}_0]}{[{\text{RANKL}}_{\text{stand}}]}=0.0947+n_{RANK{L}_0}-[{\text{OPG}}_0],\hfill \\ {\beta }_{\text{RK}}^{\text{RRO}}=\frac{C_{\text{RK}}}{[{\text{RANKL}}_{\text{stand}}]}=\frac{[{\text{RK}}_0]-[{\text{RANKL}}_0]+[{\text{OPG}}_0]}{[{\text{RANKL}}_{\text{stand}}]}=1.037-n_{RANK{L}_0}+n_{OP{G}_0},\hfill \\ {\tau }_{7\text{days}}^{\text{RRO}}=\mathrm{4.64.}\hfill \end{array}$$E15

The predicted effects of RANKL and OPG serum levels on bone density
[RANKL]; [pmol/l]	[OPG]; [pmol/l]	nRR[1]	Normalised bone; Density [1]
0.84 (standard)	1.8 (standard)	0.790	100% (0.811)
4.55	1.8	1.132	76.9% (0.624)
0.1	1.8	0.13	**172.6% (1.40)
0.84	6.9	0.276	**152.9% (0.1.24)
0.84	0.6	0.892	92.5% (0.75)

Table 1.

The predicted effects of the RANKL-RANK-OPG pathway on bone density. ${\text{n}}_{\text{RR}}$is a result from the RANKL-RANK-OPG pathway model, and consequently, bone density (the number in parentheses in the last column) is predicted from the presented thermodynamic bone remodelling model based on the calculated${\text{n}}_{\text{RR}}$. The asterisk in the front of values notices that it may be necessary to intermit the treatment after a certain time: * - after a longer time, ** - after a shorter period. Simulated or predicted data by model that are bold are in accordance with data found in literature – (Kudlacek et al, 2003).

Interconnection between this RRO model and bone remodelling model is mediated by$[\text{RR}]$. The concentration of RR influences the value of parameter ${\beta }_1$in the developed thermodynamic bone remodelling model, see (Klika and Maršík, 2010). There are different normalisations used in these two mentioned models and we assume that in the case of standard values of RANKL and OPG, the parameter ${\beta }_1$should have its standard value (corresponding to ``healthy'' state). Further, the typical normalised concentration of RR in bone remodelling model is $n_{RR}\in (\mathrm{1.35,1.41})$in standard state (see (Klika and Maršík, 2010)). Thus:

$${\beta }_1=1.41/0.79n_{RR}-\mathrm{0.81,}$$E16

which gives the value ${\beta }_1=0.6$for standard values of RANKL and OPG because $n_{RR}$under these condition equals 0.79 and $n_{RR}$is a result of the interaction in RANKL-RANK-OPG pathway at time${\tau }_{7\text{days}}^{\text{RRO}}$. As can be seen, the value of $n_{RR}$influences only${\beta }_1$, i.e. it acts only as a modification of initial conditions of the bone remodelling model. However, it will be seen in the results below that it sufficiently captures the influence of the whole pathway.

The increase in ligand concentration RANKL should lead to an increase in osteoclast formation, and consequently, the decrease of bone tissue density, and conversely, osteoprotegerin OPG prevents osteoclastogenesis. Modelling of this pathway is carried out through solving kinetic equations (2) with the above mentioned parameter values (8). Consequently, the output value of $n_{RR}$is used as an input variable in the bone adaptation model – (9). Tab. 1 gives an idea of how the added RANKL-RANK-OPG pathway may influence bone density (percentual changes of $n_{RR}$are more or less in accordance with data found in Quinn et al. (Quinn et al., 2001).

2.2. Incorporation of PTH effects into the bone remodelling model

The parathyroid hormone plays an important role in bone remodelling (see introduction). Again, the discovery of the essential control RANKL-RANK-OPG pathway enabled to describe in higher detail the influence of PTH on bone adaptation process. Robling mentions that PTH effect is ``translated'' into RANKL (Robling et al, 2006), more precisely PTH promotes RANKL formation (Martin, 2004). Fukushima et al. studied in vitro responses of the RANKL expression to PTH (Fukushima et al, 2005). They clearly showed that RANKL levels are dose-dependent on PTH concentrations in vitro (in the studied range ${10}^{-9}M$-${10}^{-6}M$). Using this data, the influence of PTH into the presented bone remodelling model will be incorporated.

As was mentioned, PTH promotes RANKL expression. Thus, we may describe this fact using the following interaction:

$$PTH+RANK{L}_{\text{producers}}+Substratum\underleftrightarrow{{\text{k}}_{\pm 1}}RANKL+RANK{L}_{\text{producers}},$$E17

where $RANK{L}_{\text{producers}}$represents the group of cells that are expressing RANKL and a mixture of substances needed for ligand RANKL production is denoted as$Substratum$. Again, the differential equation describing kinetics of PTH concentration can be derived:

$$\frac{\text{d}[\text{PTH}]}{\text{d}\tau }=-[\text{PTH}]({\beta }_{\text{Substr}}^{\text{PTH}}+[\text{PTH}])+{\delta }_{-1}^{\text{PTH}}({\beta }_{\text{RANKL}}^{\text{PTH}}-[\text{PTH}]),$$E18

where

$$\begin{array}{l}{\delta }_{-1}^{\text{PTH}}=\frac{k_{-1}}{k_{+1}[{\text{RANKL}}_{\text{stand}}]},\hfill \\ {\beta }_{\text{RANKL}}^{\text{PTH}}=\frac{C_{\text{RANKL}}}{[{\text{RANKL}}_{\text{stand}}]}=\frac{[{\text{RANKL}}_0]+[{\text{PTH}}_0]}{[{\text{RANKL}}_{\text{stand}}]},\hfill \\ {\beta }_{\text{Substr}}^{\text{PTH}}=\frac{C_{\text{Substr}}}{[{\text{RANKL}}_{\text{stand}}]}=\frac{[{\text{Substr}}_0]-[{\text{PTH}}_0]}{[{\text{RANKL}}_{\text{stand}}]}\text{.}\hfill \end{array}$$E19

Figure 2.

PTH fitted solution (blue curve) by the least squares method to data measured (dots) by Fukushima et al. (Fukushima et al, 2005).

Again, this differential equation can be rewritten into (4) where$A=\mathrm{1,}B={\beta }_{\text{Substr}}^{\text{PTH}}+{\delta }_{-1}^{\text{PTH}},C={\delta }_{-1}{\beta }_{\text{RANKL}}$. Thus, the analytical function that describes the evolution of PTH concentration in time from its initial concentration is known. The mentioned in vitro data from Fukushima et al. will be used for estimation of these parameter values using the least square method. Finally, the initial concentration of PTH influences RANKL concentration after 10 days in culture:

$$[\text{RANKL}]={\beta }_{\text{RANKL}}^{\text{PTH}}-[\text{PTH}]\text{.}$$E20

The least squares method gives the following results (see Fig. 2):

$${\delta }_{-1}^{\text{PTH}}=\mathrm{0.145,}{\tau }_{10\text{days}}^{\text{PTH}}=\mathrm{26.17,}{n}_{Subst{r}_0}^{\text{PTH}}=\mathrm{0.018.}$$E21

In vitro behaviour may significantly vary from in vivo observation. Usually, the range of concentration found in humans is much narrower than examined in vitro. A relation between in vitro and in vivo study will be found. Firstly, it is needed to determine standard in vivo serum level of PTH in humans. Khosla et al. state that PTH serum levels in humans are in $0-6\frac{pmol}{l}$(Khosla et al, 1998). Further, Jorde et al. mentions that standard in vivo PTH level in non-smokers is $3.6\frac{pmol}{l}=33.84\frac{pg}{ml}$(Jorde et al, 2005), where the fact that $M{W}_{\text{PTH}}=9400\text{daltons}$was used (Potts, 2005). Now, we need to find a study that would capture PTH effects on bone remodelling. Charopoulos et al. studied effects of primary hyperparathyroidism in both calcemic and non-calcemic groups (Charopoulos et al, 2006). They mention normal PTH levels to be $37.13\frac{pg}{ml}$and in the case of the non-calcemic hyperparathyroidism group, PTH levels were significantly higher$97.09=\frac{pg}{ml}$. Moreover, these values significantly correlated to BMD (bone mineral density).

The idea is similar as in the previous case, the influence of PTH on bone density is mediated by RANKL-RANK-OPG pathway. The predicted value of RANKL based on PTH level will be used as an input into RRO model that will consequently translate this influence into the appropriate change in number of active osteoclasts.

If one uses in vitro equivalent of in vivo standard level of PTH equal to ${510}^{-9}M$and $n_{RANK{L}_0}=1.31$(again these values are opted to get the described in vivo effect of PTH concentration and that standard PTH values would lead to standard RANKL concentration and consequently also ``healthy'' bone density) we get similar results as observed in Charopoulos et al. study - see Tab. 2.

We see that in normal range of PTH $25$-$45\frac{pg}{ml}$the predicted bone density varies slightly, see Tab. 2. However, disease state such as primary hyper-parathyriodism with PTH levels 3 times higher than normal ($97.09\frac{pg}{ml}$) the BMD dropped by 7-15% in the tibia and by 15% in the lumbar vertebrae (Charopoulos et al, 2006). The presented model estimates a decrease in bone density by 12.8% after such an increase in PTH.

The predicted effects of PTH serum levels on bone density
[PTH] [pg/ml]	normalised bone density [1]
25	102,8% (0.834)
34 (standard)	100% (0.811)
45	97.0%(0.787)
97	87.3% (0.708)

Table 2.

The predicted effects of PTH serum levels on bone density. PTH influences RANKL expression, which in turn influences osteoclastogenesis. Consequently, bone density is predicted from the presented bone remodelling model based on the calculated$[\text{RR}]$. Simulated or predicted data by model that are bold are in accordance with the data found in literature – (Charopoulos et al, 2006; Jorde et al,2005).

2.3. Incorporation of NO effects into the bone remodelling model

Nitric oxide is a small molecule that can freely diffuse across cell membrane without any use of channels (and thus without any control from cell). Probably this easy access is reason for its signalling function. Also, because there is no need of passive or even active transportation across membrane it provides very fast communications. From the mentioned nature of signalling molecule NO, it is evident that it plays role in many processes in body, and concretely, it influences bone remodelling.

There has been a quite thorough investigation of NO action on bone remodelling. Before the discovery of the RANKL-RANK-OPG pathway, van't Hof already carried out an in vitro experiment showing that NO plays a significant role in bone resorption by inducing apoptosis of osteoclast progenitors and supressing osteoclast activity (van’t Hof and Ralston, 1997). NO also seems to have a significant effect in inflammation-induced osteoporosis (such as rheumatoid arthritis), where NO production is increased, and consequently, bone density decreased (Armour et al, 1997). Wimalawansa et al. studied the effect of NO donor (concretely nitroglycerin) treatment on BMD change in rats. After treatment, they had significantly lower BMD and femur weights and also the influence of NO on the trabecular and cortical bone - loss of trabecular bone volume and decrease in cortical area in cross-section were prevented by the administration of nitroglycerin (Wimalawansa et al, 1997). Van't Hof and Ralston give a very nice overview of the NO role in bone - how nitric oxide is created, in what form, mechanism of NO action on bone remodelling, effects of NO on bone resorption (activation of inducible nitric oxide synthase - iNOS - pathway is essential for IL-1 stimulated bone resorption), and mention several in vivo animal models that exemplify this behaviour (van’t Hof and Ralston, 2001). Again, the discovery of RANKL and OPG enabled a more detailed study of NO effects on bone remodelling. Robling mentions that NO suppresses the expression of RANKL and promotes the expression of OPG (Robling et al, 2006). Both these effects lead to the decrease of bone resorption. Fan et al. carried out an in vitro study where a dose-dependent response to NO in both RANKL and OPG production has been shown (Fan et al, 2004). Rahnert et al. studied the influence of 2% exercise ($20000\mu \epsilon$) with 10 cycles/min (thus approximately$7000\mu \epsilon s^{-1}$), observed that mechanical loading represses RANKL production (i.e. resorption), and that nitric oxide influences this repression (Rahnert et al, 2008).

Thus the in vitro effects of NO are rather known; however, there is lack of clinical studies with concrete effects of NO on bone quality. Fortunately, there is a group that has recently carried out a control study aimed on NO effects on bones in humans (the NOVEL clinical study - http://clinicaltrials.gov/show/NCT00043719). They performed a pilot study of NO effects on humans (surgically induced menopause women), which had a promising results. Nitroglycerin has similar maintenance effects as estrogen on BMD (Wimalawansa, 2000). Further, there is a nice review by Wimalawansa on the treatment of bone osteoporosis with NO: nitric oxide is a very reactive agent (half life is 1 second) and thus it has a very localised effects; after menopause, the circulating NO is significantly lower than before (can be rectified after hormonal therapy); NO donors have high potential as a novel therapeutic agent in bone diseases (follows from clinical studies of his group); NO seems to have a biphasic effect on bone quality - high levels of NO generation decrease BMD whereas lower levels increase BMD (from here it implies that the possible therapeutic window is $<35mg/day$- in this range the treatment is beneficial on skeleton); animal studies are confirming that NO has a biphasic effect on BMD. Further, that nitroglycerin is superior to the other 10 studied NO donors, and showing what is the in vivo dose-response to nitroglycerin (Wimalawansa, 2007).

An estimation of NO effects on bone may be carried out (there is yet no clinical study that we may refer to) that will be based on in vitro studies and on studies with NO donor delivery in humans and animals as mentioned in the previous paragraph (in the animal study with dose response to NO donors, the same levels of nitroglycerin per kilogram of body weight were distributed and the same biphasic effect with the same threshold - $\frac{35mg}{75kg}\approx 0.5\frac{mg}{kg}$of nitroglycerin per day - were observed as in humans).

For a description of NO influences on bone remodelling, the following reaction scheme will be used. It is based on the previously mentioned observation:

$$\begin{array}{ll}NO+RANKL\hfill & \underleftrightarrow{{\text{k}}_{\pm 1}}Remaining\text{\_}product,\hfill \\ NO+OP{G}_{\text{producers}}+Substratum\hfill & \underleftrightarrow{{\text{k}}_{\pm 2}}OPG+OP{G}_{\text{producers}},\hfill \end{array}$$E22

where $Remaining\text{\_}product$represents inactivated RANKL. Again, we may simply obtain differential equations describing the kinetics of these interactions. Because the influence of NO levels are relatively small (max 6% increase in BMD (Wimalawansa, 2007)), we may separate the influence of NO on RANKL and on OPG, and their collective effect can be obtained by summing the individual ones (confirmed by the numerical solution of both ODEs at once and compared with separated effects). By the same technique (least squares - Fig. 3) as in previous cases, the unknown parameters of the mentioned scheme may be determined (normalisation in the two differential equations is quite different):

$$\begin{array}{ll}\hfill & RANKL:{\delta }_{-1}^{\text{NO},\text{R}}={2.0610}^{-12},{\tau }_{24\text{h}}^{\text{NO},\text{R}}=\mathrm{1.036,}{n}_{Remaining\text{\_}produc{t}_0}=\mathrm{0.930,}\hfill \\ \hfill & OPG:{\delta }_{-1}^{\text{NO},\text{O}}=\mathrm{528.2,}{\tau }_{24\text{h}}^{\text{NO},\text{O}}={1.610}^{-3},n_{Subst{r}_0}=\mathrm{8.50.}\hfill \end{array}$$E23

Figure 3.

NO fitted solutions (blue curve) by the least squares method to data measured (dots) by Wimalawansa (Wimalawansa, 2007) - NO influences both RANKL and OPG expression.

The in vivo standard was determined from a graph in (Wimalawansa, 2007): $0.044\frac{mg}{kg}$per day of nitroglycerin - there is no BMD change observed. This dosage of nitroglycerin treatment (and consequently NO levels) in ovariectomized rats is assumed to correspond to standard in vivo levels of NO (or a more stable$N{O}_3$). Here, the previously used linear relation between in vitro and in vivo scales does not suitably characterise the observed relationship. Therefore the following logarithmic relation was proposed and served for finding a suitable $[{\text{NO}}_{\text{in\_vitro\_stand}}]$(so that NO would have the same effects on bone remodelling as observed in vivo):

$$\begin{array}{l}{[{\text{OPG}}_0]}^{\text{RRO}}=\frac{n_{OP{G}_{stand}}}{fitNOOPG([{\text{NO}}_{\text{in\_vitro\_stand}}])}\cdot \hfill \\ \cdot fitNOOPG((1+\ln ([{\text{NO}}_{\text{in\_vivo}}]/[{\text{NO}}_{\text{in\_vivo\_stand}}])[{\text{NO}}_{\text{in\_vivo\_stand}}])\hfill \end{array}$$E24

where $fitNOOPG([{\text{NO}}_{\text{in\_vitro}}])$is a function describing in vitro effect of NO concentration on OPG, $n_{OP{G}_{stand}}=\frac{[{\text{OPG}}_{\text{stand}}]}{[{\text{RANKL}}_{\text{stand}}]}$is a normalised standard OPG molar concentration, $[{\text{NO}}_{\text{in\_vivo}}]$represents in vivo NO level corresponding to served dosage of nitroglycerin, $[{\text{NO}}_{\text{in\_vivo\_stand}}]$is the standard in vivo NO level, $\frac{n_{OP{G}_{stand}}}{fitNOOPG([{\text{NO}}_{\text{in\_vitro\_stand}}])}$represents a scaling factor guaranteeing linkage between this model of NO influence and RANKl-RANK-OPG model (with different normalisation). Thus, this relation means that a logarithmically different response in vivo to NO levels than in vitro is assumed. Similarly for RANKL:

$$\begin{array}{ll}\hfill & RANKL:{[{\text{NO}}_{\text{in\_vitro\_stand}}]}^{\text{NO},\text{R}}=\mathrm{0.036,}\hfill \\ \hfill & OPG:{[{\text{NO}}_{\text{in\_vitro\_stand}}]}^{\text{NO},\text{O}}=\mathrm{1.01.}\hfill \end{array}$$E25

From results in Tab. 3, it can be seen that the simulated effects of NO on bone density are very similar to those shown in (Wimalawansa, 2007). As was mentioned, after menopause a decrease in NO levels is observed. As can be seen from Tab. 3, this may also slightly contribute to the increase in bone resorption. Moreover, there seems to be a possibility to

The predicted effects of NG treatment on bone density
NG (nitroglycerin); [pg/ml] per day	nOPG[1]	nRANKL[1]	normalised bone density [1]
0.001	0.157	1.078	92.4% (0.749)
0.02	0.21	1.016	98.3% (0.797)
0.044 (standard)	0.233	1	100% (0.811)
0.1	0.250	0.984	102.0%(0.827)
0.5	0.281	0.953	106.0% (0.860)

Table 3.

The predicted effects of nitroglycerin (NO donor) treatment on bone density. NO promotes OPG expression and suppress the RANKL production - the same behaviour can be observed here. Consequently, bone density is predicted from the presented bone remodelling model based on the calculated$[\text{RR}]$. It is worth mentioning that NO has biphasic effect on bone, and in the window here studied, i.e. NO donor nitroglycerin is served up to $0.044\frac{\text{mg}}{\text{kg}}$per day, has anabolic effects. Simulated or predicted data by model that are bold are in accordance with data found in literature – (Wimalawansa, 2007).

decrease these unwanted effects by nitroglycerin treatment, but only to some extent because when the dosage increase approximately over $0.5\frac{mg}{kg}$per day it ceases to be beneficial. On the contrary, it may be deleterious (so-called biphasic effect of NO on skeleton; it may still have a positive impact on cardiovascular system etc.).

2.4. Notes on estradiol $1.25{(\text{OH})}_2{\text{D}}_3$(vitamin D) and TGF-${\text{β}}_1$

Estradiol

Analogously to incorporation of PTH effects on bone remodelling through RANKL-RANK-OPG pathway can be carried out so it will be omitted here and can be found in (Klika et al, 2010). However, the simulation results will be shown even with changing levels of estradiol to demonstrate the usage of this model – see section 3.

$1.25{(OH)}_2{D}_3$(vitamin D)

Similarly, as in the previous sections, we wanted to incorporate also the influences of vitamin D (more precisely$1.25{(\text{OH})}_2{\text{D}}_3$) and TGF-${\beta }_1$on the bone remodelling phenomenon. However, both these factors had to be omitted (for different reasons).

Vitamin D seems to be quite a significant factor acting on bone remodelling, e.g. seasonal effects due to sun exposure (and consequently increased vitamin D production in body) are observed (Saquib et al, 2006). Tanaka et al. showed that vitamin D promotes active osteoclast formation in vitro (Tanaka et al, 1993). There are in vitro studies that confirm again the mediation of vitamin D effects by the RANKL-RANK-OPG pathway (Martin, 2004). Moreover, it seems that there isa clinical study about to start that documents the mentioned influence on RANKL in humans but not enough of information is available yet.

But the problem is that there are many other studies that do not observe any significant influence of vitamin D levels on bone tissue, e.g. (Lips et al, 2001). One possible reason for these discrepancies in the observed results may lie in problems with its measurement. Carter mentions that there are several different methods for its determination and give significantly different results (Carter et al, 2004). In the year 2007, a thorough recherche was carried out by University of Ottawa, Evidence-based Practice Center for Agency for Healthcare Research and Quality (Cranney et al, 2007). They mention that 11 out of 19 studies are confirming that there is a significant correlation between vitamin D and bone mineral density, but also the remaining 8 are stating that there is no such effect.

Due to this ambiguity in clinical studies, we do not incorporate vitamin D into the presented model. If there is more infomation in future that would clarify these discrepancies, the vitamin D effect may be added by the same procedure as was used in estradiol, PTH, and NO.

TGF-

$${\beta }_1$$E26

There is clear evidence that this growth factor plays a significant role in the control of bone remodelling. Murakami et al. showed that the OPG expression is promoted dose-dependently by TGF-${\beta }_1$(Murakami et al, 1998) and Quinn et al. further found that RANKL expression is down-regulated by actuation of this growth factor (Quinn et al, 2001). However, clinical data collected by Zhang et al. show a fairly different behaviour of effect (Zhang et al, 2009). This different behaviour may be due to biphasic effect of TGH-${\beta }_1$that was observed by Karst et al. in vitro (Karst et al, 2004): TGF-${\beta }_1$in lower concentration elevated and in higher decreased the RANKL/OPG ratio. However, still from the mentioned clinical study it follows that actuation of TGF-$\beta$on bone remodelling is not mediated only by the RANKL-RANK-OPG pathway because response to TGF-$\beta$seems to be more complex. It is known that this factor also induces $MNOC$apoptosis (Mundy et al, 1995; Murakami et al, 1998), which is described by $J_3$in the presented model. Further, TGF-${\beta }_1$attracts osteoblastic precursors by chemotaxis to the resorbed site (Mundy et al, 1995). Even though transforming growth factor $\beta$is a component of local factors $LF$in the presented model, the influence of MNOC apoptosis can be described through$J_3$, and actuation on RANKL-RANK-OPG pathway as well, still the effects on bone tissue are not in accordance with clinically observed behaviour. Thus, the influence of TGF-$\beta$on bone remodelling cannot be described by the presented model in enough detail - it is a limitation of the model.

3. Examples of predictions of bone remodelling based on the presented model

We may now simulate the response of bone remodelling to changing environment, both mechanical and biochemical. Similarly, as was described in (Maršík et al, 2010), density distribution patterns may be obtained using Finite Elements Method (FEM). The results from the previous section will be used.

To demonstrate the usage of the presented model, a prediction of density distribution in human femur was carried out (the FE model of femur consists of 25636 3D 10 node tetrahedron elements). As an initial state, a homogeneous distribution of apparent density throughout the whole bone was used (1.8 g/cm³). Since each iteration is calculated by solving differential equations representing the whole model of bone remodelling, it is actually a time evolution of bone remodelling in bone (Ansys v11.0 program was used to calculate the mechanical stimuli in each element). Real geometry was gained from a CT-scan and external forces were applied in the usual directions and magnitude - including muscle forces acting on the femur (Heller et al., 2005). This technique was used for all computed results included in this chapter. One may observe that cortical bone (the denser part of bone) and cancellous bone are formed (and maintained) as a consequence of mechanical stimuli distribution in human bone.

It was needed to use a relation between predicted bone density and mechanical properties (elastic moduli). There is a great disparity in the proposed elastic-density relationships (Helgason et al, 2008; Rice et al, 1988; Rho et al, 1993; Hodgskinson and Currey, 1992). We used a most common relationship for the human femur because the relation seems to be site-specific (Helgason et al, 2008):$\text{E}\propto {\rho }^2$, which defines through bone density the constitutive relation being considered. If a different power law is used, the pattern of density distribution still remains the same.

As was mentioned in the introduction, the most important mechanical stimulus for maintaining bone tissue is the most common daily activity - walking. Nowadays, many people have sedentary jobs and that is why they spent less hours by walking than would be appropriate. An example of a person whose walking activity is around 55% of a healthy standard (1.5 hours of walking compared to 2.75h) is shown - inappropriate loading leads to decrease in bone density, see Fig. 4. Naturally, a possible treatment would be in spending more time walking (1.25h more) or similar effect can be reached by running, which is a higher osteogenic stimuli - 30 minutes of running every other day will compensate for the disuse - see Fig. 4.

Hyperparathyroidism was chosen as an example of biochemical control of bone remodelling. As was discussed in section 2.2, during this disease, the PTH serum levels are$97\frac{pg}{ml}$, whereas a standard value is$34\frac{pg}{ml}$. This should lead to a significant decrease in bone density. An example of a person who is physically active (correct mechanical stimuli on regular daily basis, i.e. approximately 20 000 steps per day) but suffers from hyperparathyroidism (serum level of PTH is$97\frac{pg}{ml}$) is depicted on Fig. 5. We may observe an expected decrease in bone density. Because mechanical loading has an osteogenic effect even in a disease state, we may try to increase bone density (possibly insufficiently) only with an increase in mechanical stimulation, 30 minute-running per every other day was chosen. The predicted improvement is depicted in Fig.5.

During menopause, a decline in estradiol levels occur. In some women, the decrease is very dramatic (a drop bellow $5\frac{pg}{ml}$is observed, whereas a standard serum level is $40$-$60\frac{pg}{ml}$) while in some not (serum level remains above$20\frac{pg}{ml}$), see (Klika et al, 2010). Further it was observed that, together with estradiol, there is a decline in nitric oxide levels - see section 2.3. An example of a woman who is physically active (correct mechanical stimuli on regular daily basis, i.e. approximately 20 000 steps per day) but in a consequence of menopause hasdecreased serum levels of estradiol and nitric oxide (as mentioned above - estradiol:$2.5\frac{pg}{ml}$, NO level correspond to $0.02\frac{mg}{kg}$per day of nitroglycerin) is depicted in Fig. 6. The presentedmodel predicts a decrease of 8% in bone tissue density, which does not seem to beosteoporosis yet.This may be because menopause is accompanied by more effects thanthese two mentioned and also most probably because they are less physically active (may be caused by pain). If we combine the 8% decrease (Fig. 6) caused by menopause alone with another 9% decline (Fig. 4) caused by improper loading, we get a significant drop by almost 20% in the overall bone density of the femur, which can be considered as osteoporotic state. One possible treatment of bone loss connected with menopause is treated with hormone therapy (HRT). Simulation of such a treatment that increased estradiol serum levels to $20\frac{pg}{ml}$(or by $0.107\frac{mg}{kg}$per day of nitroglycerin treatment which is actually less expensive) is given in Fig. 6. Again, the importance of mechanical stimulation shown when increased physical activity (running 30 minutes every other day) increases bone density in similar fashion as HRT treatment (the same figure). And best results are reached when both effects are combined and even the original bone tissue density can be restored - Fig. 6.

4. Conclusion

A natural goal of the modelling of a process in the human body is to help in understanding its mechanisms and ideally to help in the treatment of diseases related to this phenomenon. For this reason, more detailed influences of various biochemical factors were added. Nowadays, the RANKL/RANK/OPG chain is deemed to be one of the most important biochemical controls of the bone remodelling process. The direct cellular contact of osteoclast precursor with stromal cells is needed for osteoclastogenesis. This contact is mediated by the receptor on osteoclasts and their precursor, RANK, and ligand RANKL on osteoblasts. Osteoprotegerin binds with higher affinity to RANK which inhibits the receptor-ligand interaction and as a result, it reduces osteoclastogenesis. Thus, the raise in OPG concentration results in a smaller number of resorbing osteoclasts, which leads to higher bone tissue density. The results discussed in the presented work have exactly the same behaviour. Similarly, the effects of RANKL, RANK, and estradiol were added to the mentioned model. Consequently, simulation capabilities were demonstrated on examples of diseases and their treatment. These results were partially validated by clinical studies found in literature.

However, the impression that the presented model is able to simulate the bone remodelling process in the whole complexity is not correct. It has limitations, as mentioned below, in the spatial precision of the results (i.e. actual structure of bone tissue) and also some control mechanisms cannot be included (e.g. TGF-β effects, as was mentioned at the end of section 2.4). But still, the model can be at least considered as a summary of known important factors, comprising the fundamental aspects of the currently known knowledge of the bone remodelling phenomenon, with some predictive capabilities and encouraging predictive simulations. Since the presented model is a concentration model, it cannot be used arbitrarily. The limitation is, of course, in the spatial precision of results. The minimal volume unit (finite element) should be sufficiently large to contain enough of all the

Figure 4.

Prediction of disuse effect on bone quality - simulation of insufficient loading (half of the recommended daily stimulation), treatment proposal - running (30 minutes every other day), and its simulation. Notice the change of bone mass (BM) of the whole femur.

Figure 5.

Prediction of the PTH effect on bone quality - simulation of hyperparathyroidism ($97\frac{\text{pg}}{\text{ml}}$where normal PTH levels are$34\frac{\text{pg}}{\text{ml}}$), treatment proposal - running (30 minutes every other day), and its simulation. Notice the change of bone mass (BM) of the whole femur.

Figure 6.

Prediction of the menopause effect on bone quality (estradiol levels decreased to$2.5\frac{\text{pg}}{\text{ml}}$and NO to half of its normal level), treatment proposal, and its simulation - hormonal treatment (HRT), running (30 minutes every other day). Notice the change of bone mass (BM) of the whole femur.

substances entering the reaction schemes, namely osteoclasts and osteoblasts. It surely cannot be used on the length scales of BMU where it is no longer guaranteed that any osteoclast is present. There are approximately 10⁷ BMU in a human skeleton present at any moment (Klika and Maršík, 2010) and, because bones have a total volume of 1.75l, there is 1 BMU per 0.175 mm³ on average at any moment. In other words, the presented model cannot be used for length scales smaller than $\sqrt[3]{0.175}$mm³ and we recommend that it is not used on length scales below $\sqrt[3]{0.5}$mm³$\approx$0.8 mm.

Ongoing applications of the model include simulations of the 3D geometries of the femur and vertebrae (FE models) under various conditions (both biochemical and mechanical). The preliminary results are encouraging and show the correct density distribution. Currently, we are working on bone modelling (change of shape of bone) model that would add the possibility to adapt bone to its mechanical environment as it is observed in vivo. Further, we would like to have a more detailed description of the inner structure of bone as an outcome of the model. Most probably, a homogenisation technique will be used for addressing this goal. Most importantly, a validation of the model‘s predictions (or finding limits of its application) should start in near future in cooperation with Ambulant Centre for Defects of Locomotor Aparatus in Prague.