A coupled mechano-biochemical model for bone adaptation

Abstract

Bone remodelling is a fundamental biological process that controls bone microrepair, adaptation to environmental loads and calcium regulation among other important processes. It is not surprising that bone remodelling has been subject of intensive both experimental and theoretical research. In particular, many mathematical models have been developed in the last decades focusing in particular aspects of this complicated phenomenon where mechanics, biochemistry and cell processes strongly interact. In this paper, we present a new model that combines most of these essential aspects in bone remodelling with especial focus on the effect of the mechanical environment into the biochemical control of bone adaptation mainly associated to the well known RANKL-RANK-OPG pathway. The predicted results show a good correspondence with experimental and clinical findings. For example, our results indicate that trabecular bone is more severely affected both in disuse and disease than cortical bone what has been observed in osteoporotic bones. In future, the methodology proposed would help to new therapeutic strategies following the evolution of bone tissue distribution in osteoporotic patients.

Appendices

Appendix A: Relation between daily strain history $\xi$ and the mechanical parameter $\mathcal {D}_{\alpha }$

The mechanical stimulus used in the biochemical model was

$$\begin{aligned} \mathcal {D}_{\alpha }=l_{v\alpha }d^{(1)}, \end{aligned}$$

(23)

where $d^{(1)}$ is the rate of volume variation and $l_{v\alpha }$ are given in Table 8. In the mechano-biochemical model presented here, we generalise this relation by including influence of shear, see (6) and (7). The relation between the newly used form of mechanical stimulation $\xi$ and the previous one $\mathcal {D}_{\alpha }$ is the following:

$$\begin{aligned} \mathcal {D}_{\alpha }=A l_{v \alpha } \xi = \frac{\mathcal {D}_{\alpha ,ref}}{\xi _{ref}} \xi , \end{aligned}$$

where the value of $A$ was determined from the reference healthy states, see the end of Sect. 2.2. The values of $\mathcal {D}_{\alpha ,ref}$ are calculated from (23) where the value of $d^{(1)}$ corresponds to the maximal strain rate found in a healthy mechanical stimulation, see Sect. 2.2 and the original model formulation (Klika and Maršík 2010):

$$\begin{aligned} \mathcal {D}_{1,ref}&= -4.58,\\ \mathcal {D}_{2,ref}&= 1.16\cdot 10^{-4},\\ \mathcal {D}_{3,ref}&= -1.31\cdot 10^{-2},\\ \mathcal {D}_{4,ref}&= -5.22\cdot 10^{-3},\\ \mathcal {D}_{5,ref}&= 2.04\cdot 10^{-4}. \end{aligned}$$

Appendix B: Biochemical model

The model presented here is an updated version of (Klika and Maršík 2010; Klika et al. 2010; Klika and Maršík 2011a, b). The core of the biochemical part of the presented model is characterised by the following interactions:

$$\begin{aligned} RR + MCELL&\mathop {\rightarrow }\limits ^{\mathrm{k}_{+1}}&MNOC + N_4,\nonumber \\ MNOC + Old\_B&\mathop {\rightarrow }\limits ^{\mathrm{k}_{+2}}&LF + N_7,\nonumber \\ LF + osteoprogenitor&\mathop {\rightarrow }\limits ^{\mathrm{k}_{+3}}&OB + N_{10},\nonumber \\ OB + N_{11}&\mathop {\rightarrow }\limits ^{\mathrm{k}_{+4}}&Osteoid + N_{13},\nonumber \\ N_{14} + Osteoid&\mathop {\rightarrow }\limits ^{\mathrm{k}_{+5}}&New\_B + N_{16}, \end{aligned}$$

(24)

where $MNOC$ are multi-nucleated osteoclasts being formed by fusion from their precursors $MCELL$ when appropriate receptor-ligand (cell-to-cell) interaction takes place ($RR$ representing RANKL-RANK interaction mediated by the RANKL-RANK-OPG chain), $Old\_B$ denotes old bone that might be resorbed, the released local factors from bone are denoted by $LF$ which activate $osteoprogenitor$ cells that proliferate and differentiate into osteoblasts $OB$, which secrete $osteoid$ (non-mineralised bone - organic part of bone tissue). The ossification of $osteoid$ into new bone tissue $New\_B$ happens after addition of appropriate substances, $N_{14}$, takes place. $N_4, N_7, N_{10}, N_{13}, N_{16}$ are waste products. For a more detailed description see (Klika and Maršík 2010).

The RANKL-RANK-OPG pathway is considered as follows:

$$\begin{aligned}&RANKL + RANK~\mathop {\rightleftarrows }\limits ^{\mathrm{k}_{\pm 1}} RR, \end{aligned}$$

(25)

$$\begin{aligned}&\quad RANKL + OPG~\mathop {\rightleftarrows }\limits ^{\mathrm{k}_{\pm 2}} RO_{\mathrm{inactive}}, \end{aligned}$$

(26)

Finally, the effect of estradiol, PTH, and NO is translated through this pathway:

$$\begin{aligned} Estradiol + OPG_{\mathrm{prod}}+Subst&\mathop {\rightleftarrows }\limits ^{\mathrm{k}_{\pm 1}}&OPG+OPG_{\mathrm{prod}},\end{aligned}$$

(27)

$$\begin{aligned} PTH + RANKL_{\mathrm{prod}}+Subst&\mathop {\rightleftarrows }\limits ^{\mathrm{k}_{\pm 1}}&RANKL+RANKL_{\mathrm{prod}},\end{aligned}$$

(28)

$$\begin{aligned} NO + RANKL~&\mathop {\rightleftarrows }\limits ^{\mathrm{k}_{\pm 1}}&Remaining\_prod,\end{aligned}$$

(29)

$$\begin{aligned} NO + OPG_{\mathrm{prod}}+Subst&\mathop {\rightleftarrows }\limits ^{\mathrm{k}_{\pm 2}}&OPG+OPG_{\mathrm{prod}}, \end{aligned}$$

(30)

where $OPG_{\mathrm{prod}},~RANKL_{\mathrm{prod}}$ represent the group of cells that expresses OPG, RANKL and a mixture of substances needed for osteoprotegerin, RANKL production is noted as $Subst$, respectively and $Remaining\_prod$ represents inactivated RANKL.

The dimensionless system of differential equations describing the kinetics of interactions of the core of the biochemical part of the presented model model (24) together with the influence of mechanical stimulus $\mathcal {D}_{\alpha }(\xi )$ (for the relationship between $\mathcal {D}_{\alpha }$ and $\xi$ see appendix 5) is provided here

$$\begin{aligned} \frac{\mathrm{d}n_{MCELL}}{\mathrm{d}\tau }&= -\delta _1(\beta _1+n_{MCELL})n_{MCELL}+\mathcal {J}_3\nonumber \\&+\mathcal {J}_{New\_B}-\mathcal {D}_1(\xi )\\ \frac{\mathrm{d}n_{Old\_B}}{\mathrm{d}\tau }&= -(\beta _3-n_{MCELL}+n_{Old\_B})n_{Old\_B}-\mathcal {D}_2(\xi )\nonumber \\&+\mathcal {J}_{New\_B} \\ \frac{\mathrm{d}n_{OB}}{\mathrm{d}\tau }&= \delta _3(\beta _6-n_{Old\_B}-(n_{OB}+n_{Osteoid}+n_{New\_B}))\nonumber \\&\times (\beta _8-(n_{OB}+n_{Osteoid}+n_{New\_B}))\\&-\delta _4(\beta _{11}-(n_{Osteoid}+n_{New\_B}))n_{OB}+\mathcal {D}_3(\xi )-\mathcal {D}_4(\xi ) \nonumber \\ \frac{\mathrm{d}n_{Osteoid}}{\mathrm{d}\tau }&= \delta _4(\beta _{11}-(n_{Osteoid}+n_{New\_B}))n_{OB}\nonumber \\&-\delta _5 (\beta _{14}-n_{New\_B})n_{Osteoid}+\mathcal {D}_4(\xi )-\mathcal {D}_5(\xi ) \\ \frac{\mathrm{d}n_{New\_B}}{\mathrm{d}\tau }&= \delta _5 (\beta _{14} -n_{New\_B})n_{Osteoid}-\mathcal {J}_{New\_B}+\mathcal {D}_5(\xi ), \end{aligned}$$

where $\beta _{i}$ is the sum of normalised initial molar concentrations of relevant substances, $\delta _{\alpha }$ relates to the interaction rate, $D_{\alpha }(\xi )$ describes the influence of dynamic loading on reactions, $n_{i}$ is the normalised concentration of the $i$-th substance, $\mathcal {J}_3$ is a normalised outflow of $MNOC$ (i.e., $MNOC$ apoptosis) and $\mathcal {J}_{New\_B}$ is the normalised outflow of new bone $New\_B$ into old bone $Old\_B$.

By solving these kinetic equations, the time evolution of $[\mathrm{MCELL}], [\mathrm{Old\_B}], [\mathrm{OB}], [\mathrm{Osteoid}]$, and $[\mathrm{New\_B}]$ concentrations can be obtained: the relation between concentrations and their dimensionless form is:

$$\begin{aligned} n_{N_i}=[\mathrm{N_i}] \frac{1}{k_{+2} [\mathrm{Bo}]^2}, \end{aligned}$$

where $k_{+2}=6 \cdot 10^7 \frac{l}{mol\cdot s}$ and $[\mathrm{Bo}]=2.2 \frac{10^{-4}}{N_A} \frac{mol}{\mu m^3}$ with $N_A$ being the Avogardo’s number. In this model, it is assumed that the density of bone tissue is proportional to the number of osteocytes in a given volume, or that the molar concentration of osteocytes in tissue is constant; the amount of bone tissue in a given volume then determines the local bone density. Thus the quantity $[\mathrm{Old\_B}]$ is (by definition) equal to the concentration of osteocytes $[\mathrm{OCy}]$ which, with respect to above mentioned assumption, gives us the concentration of bone tissue $Old\_B$:

$$\begin{aligned}{}[\mathrm{Old\_B}]=\frac{2.3}{0.81} n_{Old\_B} \left[ \frac{kg}{m^3}\right] \end{aligned}$$

which directly relates the maximal normalised $n_{Old\_B}$ to maximal apparent bone tissue density $[\mathrm{Old\_B}]$ found in the human femur. The kinetics of other interactions is derived from the Law of Mass Action or can be found in the referred literature. Here, to provide a stand-alone model, the remaining differential equations of the biochemical model describing the considered kinetics of RANKL-RANK-OPG pathway are recapitulated.

RANKL-RANK-OPG pathway The law of mass action gives for the RANKL-RANK-OPG pathway (25) the following system of differential equations where we denote the dimensionless concentration of RANKL and OPG as $n_{RKL}$ and $n_{OPG}$, respectively:

$$\begin{aligned} \frac{\mathrm{d}n_{RKL}}{\mathrm{d}\tau }&= -n_{RKL}\left( \beta _{\mathrm{RK}}^{\mathrm{RRO}}+n_{RKL}-n_{OPG}\right) \nonumber \\&+\delta _{-1}^{\mathrm{RRO}}\left( \beta _{\mathrm{RR}}^{\mathrm{RRO}}-n_{RKL}+n_{OPG}\right) \nonumber \\&-{}\delta _{+2}^{\mathrm{RRO}}n_{RKL}n_{OPG}+\delta _{-2}^{\mathrm{RRO}}\left( \beta _{\mathrm{RO}}^{\mathrm{RRO}}-n_{OPG}\right) ,\\ \frac{\mathrm{d}n_{OPG}}{\mathrm{d}\tau }&= -\delta _{+2}^{\mathrm{RRO}}n_{RKL}n_{OPG}+\delta _{-2}^{\mathrm{RRO}}\left( \beta _{\mathrm{RO}}^{\mathrm{RRO}}-n_{OPG}\right) ,\nonumber \end{aligned}$$

with $n_{RR}(t)=\beta _{\mathrm{RR}}^{\mathrm{RRO}}-n_{RKL}(t)+n_{OPG}(t)$ and parameter values are given in Appendix 8. The amount of RANKL-RANK bonds ($RR$) mediate the outcome effect of the whole RANKL-RANK-OPG pathway on bone remodelling. The amount of $RR$ bonds is included in the $\beta _1$ parameter. The differential equations are in dimensionless form and to provide prediction of changes in real serum levels they need to be compared with the physiological value of RANKL that was used for obtaining the dimensionless form. In total, the initial conditions are given by relationships

$$\begin{aligned} n_{RKL_0}&=\frac{[\mathrm{RKL}]}{[\mathrm{RKL}]_{st}}\underbrace{n_{RKL,st}}_{1},\\ n_{OPG_0}&=\frac{[\mathrm{OPG}]}{[\mathrm{OPG}]_{st}}\underbrace{n_{OPG,st}}_{0.2333}, \end{aligned}$$

and standard serum levels in humans were estimated to be $[\mathrm{RKL}]_{st}=0.84\frac{pmol}{l}=46.2\frac{pg}{ml}$ and $[\mathrm{OPG}]_{st}=1.8\frac{pmol}{l}=36\frac{pg}{ml}$. The normalised standard OPG concentration $n_{OPG,st}$ is determined from molar masses of RANKL and OPG, see (Klika et al. 2010). Finally, the interconnection with the core model is

$$\begin{aligned} \beta _1=(1.41/0.79) n_{RR}\left( \tau _{7days}^{RRO}\right) - 0.81, \end{aligned}$$

giving $\beta _1=0.6$ for standard serum values of RANKL and OPG. For more details see (Klika et al. 2010).

As was mentioned, many of the factors having a significant influence on RANKL-RANK-OPG signalling pathway, and thus on bone quality, are translated through this same pathway. Below we will add impacts of PTH, estradiol and nitric oxide. Here, we model their impact by modifying the initial values of RANKL and OPG in this RANKL-RANK-OPG model as predicted from the following sub-models. The final initial conditions for RANKL and OPG are determined by the sum of all predicted deflections from their standard values, namely

$$\begin{aligned} n_{RKL_0}&= n_{RKL,st}+\left( n_{RKL_0}^{RRO,PTH}-n_{RKL,st}\right) + \left( n_{RKL_0}^{RRO,NO}-n_{RKL,st}\right) ,\qquad \end{aligned}$$

(31)

$$\begin{aligned} n_{OPG_0}&= n_{OPG,st}+\left( n_{OPG_0}^{RRO,estr}-n_{OPG,st}\right) + \left( n_{OPG_0}^{RRO,NO}-n_{OPG,st}\right) ,\qquad \end{aligned}$$

(32)

where the terms in brackets denote differences from standard initial values caused by the control mechanism noted in the upper indices.

Further, as each following sub-model has its own normalisation there are two scaling factors present in each one: (1) scaling of used in-vitro data by $C_{in-vitro}$ (it is unclear if this scaling finds an equivalent amount, in terms of its effects on bone remodelling, of in-vivo serum levels to that of used in-vitro ones but it seems to capture their impact on bone remodelling as reported in literature), (2) scaling to match the different normalisation used— matching is such that standard value of a given factor are met in both normalisations by $C_{norm}$, see below.

PTH The impact of PTH on RANKL-RANK-OPG pathway (28) can be described as follows:

$$\begin{aligned} \frac{\mathrm{d}{\tilde{n}_{PTH}}}{\mathrm{d}\tau }=-{\tilde{n}_{PTH}}\left( \beta _{\mathrm{Substr}}^{\mathrm{PTH}}+{\tilde{n}_{PTH}}\right) +\delta _{-1}^{\mathrm{PTH}}\left( \beta _{\mathrm{RKL}}^{\mathrm{PTH}}-{\tilde{n}_{PTH}}\right) , \end{aligned}$$

which affects the amount of RANKL available for making receptor-ligand bonds. This equation for normalised in-vitro concentration $\tilde{n}_{PTH}$ was fitted with in-vitro dose-dependent data (no other available). The rescaling of in-vitro data is included in ${\tilde{n}_{PTH}}=C_{in-vitro}^{PTH} n_{PTH}$, $C_{in-vitro}^{PTH}=0.296$ and thus the initial condition needed for solving the above differential equation is

$$\begin{aligned} {\tilde{n}_{PTH_0}}=C_{in-vitro}^{PTH} \frac{[\mathrm{PTH}]_0}{[\mathrm{PTH}]_{st}}, \end{aligned}$$

(33)

where $[\mathrm{PTH}]$ represents serum levels of PTH, typically in $\frac{pg}{ml}$. The concentration of RANKL after action of PTH is ($n_{PTH}(\tau )$ being the scaled solution of the above differential equation):

$$\begin{aligned} n_{RKL_0}^{RRO,PTH}&= C_{norm}^{PTH}\left( \beta _{\mathrm{RKL}}^{\mathrm{PTH}}-n_{PTH}\left( \tau _{10days}^{PTH}\right) \right) \\ {}&= C_{norm}^{PTH}\left( n_{RKL_0}^{PTH}+n_{PTH_0}-n_{PTH}\left( \tau _{10days}^{PTH}\right) \right) , \end{aligned}$$

where $n_{RKL_0}^{PTH}$ represents the initial concentration of RANKL before interacting with PTH according to (28) and $C_{norm}^{PTH}=1.0$ is the second scaling to match the different normalisations (see above). The amount of available RANKL for forming bonds with its receptor RANK is denoted as $n_{RKL_0}^{RRO,PTH}$ and is used as an input into the RANKL-RANK-OPG model taking into account the influence of PTH on RANKL concentration, see (31)–(32). Parameter values can be found again in Appendix 8.

Estradiol The idea behind adding an impact of estradiol on RANKL-RANK-OPG pathway (27) is the same as in PTH, only that estradiol affects OPG instead of RANKL (over-tilde has the same meaning as in PTH—rescaled values of estradiol):

$$\begin{aligned} \frac{\mathrm{d}{\tilde{n}_{estr}}}{\mathrm{d}\tau }=-{\tilde{n}_{estr}}(\beta _{\mathrm{Substr}}^{\mathrm{estr}}+{\tilde{n}_{estr}})+ \delta _{-1}^{\mathrm{estr}}\left( \beta _{\mathrm{OPG}}^{\mathrm{estr}}-{\tilde{n}_{estr}}\right) , \end{aligned}$$

and the initial condition is obtained from serum levels analogously to that one in PTH, Eq. (33), with $C_{in-vitro}^{estr}=1.19$. The concentration of OPG after action of estradiol is

$$\begin{aligned} n_{OPG_0}^{RRO,estr}&= C_{norm}^{estr}\left( \beta _{\mathrm{OPG}}^{\mathrm{estr}}-n_{estr}\left( \tau _{\mathrm{24h}}^{\mathrm{estr}}\right) \right) \\ {}&= C_{norm}^{estr}\left( n_{OPG_0}^{estr}+n_{estr_0}-n_{estr}\left( \tau _{\mathrm{24h}}^{\mathrm{estr}}\right) \right) , \end{aligned}$$

with $C_{norm}^{estr}=1.04$. Again the amount of available OPG for blocking bond formation is denoted as $n_{OPG_0}^{RRO,estr}$ and is used as an input into the RANKL-RANK-OPG model taking into account the influence of estradiol on OPG concentration, see (31)-(32). Parameter values can be found again in Appendix 8.

Nitric oxide Similarly, the influence of nitric oxide on RANKL-RANK-OPG pathway (29)–(30) can be included. It affects both RANKL and OPG concentrations. It was assumed that the effects of NO on OPG and RANKL can be separated (Klika and Maršík 2011a) leading to the same differential equations as in previous cases for estradiol and PTH. Due to unavailability of in-vivo data that could quantify NO effects on bone quality we used studies where the amount of NO was regulated by NO donor intake in milligrams per day (typically nitroglycerin). A non-linear scaling in the input data was required to obtain a fit with observed outcome effects on BMD. It is unclear if this scaling finds an equivalent amount (in terms of its effects on bone remodelling) of in-vivo NO serum levels to that of donor but it seems to capture the impact of NO on bone remodelling. Concretely, the scaling transforms were logarithmic ${\tilde{n}_{NO,OPG}}=C_{in-vitro}^{NO,OPG}\left( 1+\ln \left( [\mathrm{NO^{OPG}}]/[\mathrm{NOdonor}]_{stand}\right) \right)$ and ${\tilde{n}_{NO^{RKL}}}=C_{in-vitro}^{NO,RKL}\left( 1+\ln \left( [\mathrm{NO^{RKL}}]/[\mathrm{NOdonor}]_{stand}\right) \right)$ with normalisation constants $C_{in-vitro}^{NO,OPG}=1.01$, $C_{in-vitro}^{NO,RKL}=0.036$ and with $[\mathrm{NOdonor}]_{stand}=0.044\frac{mg}{kg\times day}$ being the standard or reference amount of NO donor that does not cause any change in BMD (Wimalawansa 2007).

Finally, the differences from standard values of OPG and RANKL caused by NO are

$$\begin{aligned} n_{RKL_0}^{RRO,NO}&= C_{norm}^{NO,RKL}\left( \beta _{\mathrm{RKL}}^{\mathrm{NO,RKL}}-n_{NO,RKL}\left( \tau _{\mathrm{24h}}^{\mathrm{NO,RKL}}\right) \right) \\&= C_{norm}^{NO,RKL}\left( n_{OPG_0}^{NO,RKL}+n_{NO,RKL_0}-n_{NO,RKL}\left( \tau _{\mathrm{24h}}^{\mathrm{NO,RKL}}\right) \right) ,\\ n_{OPG_0}^{RRO,NO}&= C_{norm}^{NO,OPG}\left( \beta _{\mathrm{OPG}}^{\mathrm{NO,OPG}}-n_{NO,OPG}\left( \tau _{\mathrm{24h}}^{\mathrm{NO,OPG}}\right) \right) \\&= C_{norm}^{NO,OPG}\left( n_{OPG_0}^{NO,OPG}+n_{NO,OPG_0}-n_{NO,OPG}\left( \tau _{\mathrm{24h}}^{\mathrm{NO,OPG}}\right) \right) , \end{aligned}$$

where the two normalisation constants are $C_{norm}^{NO,RKL}=1.02$, $C_{norm}^{NO,OPG}=2.16$. Again the amount of available OPG for blocking bond formation is denoted as $n_{OPG_0}^{RRO,NO},~n_{RKL_0}^{RRO,NO}$ and is used as an input into the RANKL-RANK-OPG model taking into account the influence of nitric oxide on this essential control pathway, see (31)–(32). Parameter values can be found in Appendix 8.

Implementation of the whole model in Abaqus as a user subroutine is enclosed in the supplementary material.

Appendix C: Stationary solution

The biochemical model, see Sect. 6, has exactly one positive stationary solution:

$$\begin{aligned} \overline{n_{MCELL}}&= 1/2\Bigg (-\beta _1+\sqrt{\beta _1^2+4\frac{-\mathcal {D}_1(\xi )+\mathcal {J}_3+\mathcal {J}_{14}}{\delta _1}}\Bigg )\nonumber \\ \overline{n_{Old\_B}}&= 1/2\Bigg (\bigg (-\beta _3+\overline{n_{MCELL}}\bigg )\nonumber \\&+\sqrt{\bigg (-\beta _3+\overline{n_{MCELL}}\bigg )^2+4(\mathcal {J}_{14}-\mathcal {D}_2(\xi ))}\Bigg )\nonumber \\ Sum_1&= 1/2\Bigg (\bigg (\beta _8+\beta _6-\overline{n_{Old\_B}}\bigg )\\&{}-\sqrt{\bigg (\beta _8-(\beta _6-\overline{n_{Old\_B}})\bigg )^2+4\frac{\mathcal {J}_{14}-\mathcal {D}_3(\xi )}{\delta _3}}\Bigg ) \nonumber \\ Sum_2&= 1/2\Bigg (\bigg (\beta _{11}+Sum_1\bigg )-\sqrt{\bigg (\beta _{11}-Sum_1\bigg )^2+4\frac{\mathcal {J}_{14}-\mathcal {D}_4(\xi )}{\delta _4}}\Bigg ) \nonumber \\ \overline{New\_B}&= 1/2\Bigg (\bigg (\beta _{14}+Sum_2\bigg )-\sqrt{\bigg (\beta _{14}-Sum_2\bigg )^2+4\frac{\mathcal {J}_{14}-\mathcal {D}_5(\xi )}{\delta _5}}\Bigg ) \nonumber \\ \overline{n_{Osteoid}}&= Sum_2-\overline{New\_B}\nonumber \\ \overline{n_{OB}}&= Sum_1-Sum_2.\nonumber \end{aligned}$$

which for the considered parameter values fulfil all the assumptions of Poincare-Ljapunov theorem and thus the stationary solution is an asymptotically stable fixed point.

Appendix D: List of parameter values

$$\begin{aligned}&\text{ initial } \text{ values } \text{ of } \text{ damage, } v_b, \alpha \text{(ash } \text{ fraction) }\\&\quad d_0=0.0, v_b^0=0.6, \alpha _0=0.69\\&\mathbf{Biochemical}\,\, \mathbf{part}\,\, \mathbf{of}\,\, \mathbf{the}\,\, \mathbf{model}\\&\quad [\mathrm{PTH}]_{st}=34.0\frac{pg}{ml}, [\mathrm{estr}]_{st}=50.0\frac{pg}{ml},[\mathrm{NO_{donor}}]_{st}=0.044\frac{mg}{kg\cdot day},\\&\qquad [\mathrm{RKL}]_{st}=46.2\frac{pg}{ml},[\mathrm{OPG}]_{st}=36.0\frac{pg}{ml}\\&\text{ Estradiol } \text{ submodel }\\&\quad \tau _{24h}^{estr}=26.166,\delta _1^{estr}=0.1448, nSUBST_0^{estr}=0.018\\&\text{ NO-OPG } \text{ submodel }\\&\quad \tau _{24h}^{NO,OPG}=0.00177,\delta _1^{NO,OPG}=476.7, nSUBST_0^{NO,OPG}=7.67\\&\text{ NO-RANKL } \text{ submodel }\\&\quad \tau _{24h}^{NO,RKL}=50.0, \delta _1^{NO,RKL}=0.8, nREMAIN_0^{NO,RKL}=1E-5\\&\text{ PTH } \text{ submodel }\\&\quad \tau _{10days}^{PTH}=1687.3, \delta _1^{PTH}=1.763, nSUBST_0^{PTH}=2.6,\\&\quad nRKL_0^{PTH}=1.31\\&\text{ RANKL-RANK-OPG } \text{ submodel }\\&\quad \tau _{7days}^{RRO}=4.64,\delta _1^{RRO}=4.92E-6, \delta _{2+}^{RRO}=12.96,\\&\quad \delta _{2-}^{RRO}=5.86E-19,RKL_{Stand}=1.0,OPG_{Stand}=0.2333\\&\text{ main } \text{ model }\\&\quad \tau _{Day}=1 \text{ day },\tau _{BR}=315 \text{ days }\\&\quad MCELL_0=0.8,osteoprog_0=0.0,Osteoid_0=0.0,\\&\quad \delta _1=4.0,\delta _3=2.0,\delta _4=0.142857,\delta _5=0.142857,\\&\quad \beta _6=1.05,\beta _8=0.1,J_{Old\_B}=2.6E-4,J_3=4/25 J_{Old\_B},\\&\quad \beta _{11}=1.0,\beta _3=0.0005,\\&\quad \beta _{14}=0.05,\mathcal {D}_1^{ref}=-4.58,\mathcal {D}_2^{ref}=1.16E-4,\mathcal {D}_3^{ref}=-1.31E-2,\\&\quad \mathcal {D}_4^{ref}=-5.22E-3,\mathcal {D}_5^{ref}=2.04E-4 \end{aligned}$$