A multiscale mechanobiological model of bone remodelling predicts site-specific bone loss in the femur during osteoporosis and mechanical disuse

Abstract

We propose a multiscale mechanobiological model of bone remodelling to investigate the site-specific evolution of bone volume fraction across the midshaft of a femur. The model includes hormonal regulation and biochemical coupling of bone cell populations, the influence of the microstructure on bone turnover rate, and mechanical adaptation of the tissue. Both microscopic and tissue-scale stress/strain states of the tissue are calculated from macroscopic loads by a combination of beam theory and micromechanical homogenisation. This model is applied to simulate the spatio-temporal evolution of a human midshaft femur scan subjected to two deregulating circumstances: (i) osteoporosis and (ii) mechanical disuse. Both simulated deregulations led to endocortical bone loss, cortical wall thinning and expansion of the medullary cavity, in accordance with experimental findings. Our model suggests that these observations are attributable to a large extent to the influence of the microstructure on bone turnover rate. Mechanical adaptation is found to help preserve intracortical bone matrix near the periosteum. Moreover, it leads to non-uniform cortical wall thickness due to the asymmetry of macroscopic loads introduced by the bending moment. The effect of mechanical adaptation near the endosteum can be greatly affected by whether the mechanical stimulus includes stress concentration effects or not.

Appendices

Appendix 1: Complements on the model description

1.1 Differentiation rates and signalling molecules in the cell populations model

In Sect. 2.3, we presented the simplified equations of the bone cells population model. Here are the developments of these equations.

$$\begin{aligned}&\mathcal {D}_{\text {OC}_{\mathrm{u}}}\big (\text {MCSF}, \text {RANKL}(\varPsi ,\text {PTH})\big )\nonumber \\&\quad = D_{\text {OC}_{\mathrm{u}}} \pi ^\text {act}\big (\tfrac{\text {MCSF}}{k^{\text {MCSF}}_{\text {OC}_{\text {u}}}}\big ) \pi ^\text {act}\big (\tfrac{\text {RANKL}}{k^{\text {RANKL}}_{\text {OC}_{\text {u}}}} \big ),\nonumber \\&\mathcal {D}_{\text {OC}_{\mathrm{p}}}\big (\text {RANKL}(\varPsi ,\text {PTH})\big ) = D_{\text {OC}_{\mathrm{p}}} \pi ^\text {act}\big (\tfrac{\text {RANKL}}{k^{\text {RANKL}}_{\text {OC}_{\text {p}}}}\big ),\nonumber \\&\mathcal {A}_{\text {OC}_{\mathrm{a}}}(\text {TGF} \upbeta ) = A_{\text {OC}_{\mathrm{a}}}\pi ^\text {act}\big (\tfrac{\text {TGF} \upbeta }{k^{\text {TGF} \upbeta }_{\text {OC}_{\text {a}}}}\big ),\nonumber \\&\mathcal {D}_{\text {OB}_{\mathrm{u}}}(\text {TGF} \upbeta ) = D_{\text {OB}_{\mathrm{u}}}\pi ^\text {act}\big (\tfrac{\text {TGF} \upbeta }{k^{\text {TGF} \upbeta }_{\text {OB}_{\text {u}}}}\big ),\nonumber \\&\mathcal {D}_{\text {OB}_{\mathrm{p}}}(\text {TGF} \upbeta ) = D_{\text {OB}_{\mathrm{p}}}\pi ^{\text {rep}}\big (\tfrac{\text {TGF} \upbeta }{k^{\text {TGF} \upbeta }_{\text {OB}_{\text {p}}}}\big ). \end{aligned}$$

(23)

In those equations, several signalling molecules play a role: $\text {TGF} \upbeta$, RANK, RANKL, OPG, MCSF and PTH. The concentrations of these molecules follow the principles of mass action kinetics of receptor–ligand reactions. Due to the separation of scale between the cells differentiation and apoptosis rates and the receptor–ligand binding reactions, we solve them in a quasi-steady-state hypothesis:

$$\begin{aligned}&{\text {PTH}}({\varvec{r}}, t) = {\left\{ \begin{array}{ll} {\text {P}}_{\mathrm{PTH}}, &{}\quad \text {without deregulation} \\ {\text {P}}_{\mathrm{PTH}}^{\mathrm{OP}}, &{}\quad \text {with simulated OP} \end{array}\right. }, \end{aligned}$$

(24)

$$\begin{aligned}&\text {TGF} \upbeta ({\varvec{r}}, t) =\frac{{P^{\mathrm{ext}}_{{\mathrm{TGF}} \upbeta }} + n^{\mathrm{bone}}_{\mathrm{TGF} \upbeta } ~k_\text {res}~{\text {OC}}_{\mathrm{a}}}{D_{\mathrm{TGF} \upbeta }}, \end{aligned}$$

(25)

$$\begin{aligned}&{\text {RANK}}({\varvec{r}}, t) = N^{\mathrm{RANK}}_{{\mathrm{OC}}_{\mathrm{p}}}\, \text {OC}_{\mathrm{p}}, \end{aligned}$$

(26)

$$\begin{aligned}&\text {OPG}({\varvec{r}}, t) = \frac{{P_{\text {OPG}} + \beta ^{\text {OPG}}_{\text {OB}_{\text {a}}}\, \text {OB}_{\mathrm{a}}\, \pi ^{\text {rep}}\!\Big (\frac{\text {PTH}}{k^{\text {PTH}}_{\text {OB}}}\Big )}}{\beta ^{\text {OPG}}_{\text {OB}_{\text {a}}}\,\text {OB}_{\mathrm{a}}\,\pi ^{\text {rep}}\!\Big (\frac{\text {PTH}}{k^{\text {PTH}}_{\text {OB}}}\Big )/\text {OPG}_\text {sat} + D_{\text {OPG}}}, \end{aligned}$$

(27)

(28)

Table 1 Nomenclature

1.2 Model parameter values

Table 1 lists all the parameters of the model.

Most model parameters were calibrated and validated in previous studies. The stiffness tensors $\mathbbm {c}^{\text {micro}}_\text {bm}$ and $\mathbbm {c}^{\text {micro}}_\text {vas}$ arising in the micromechanical theory were measured by mechanical experiments (Ashman et al. 1984; Murdock 1996); see also (Fritsch and Hellmich 2007). The secretory rate $\text {k}_{\mathrm{form}}$(volume of bone matrix formed per osteoblast per day) was measured in animal models (Marotti and Zallone 1976; Jones 1974). The value used here is based on osteonal data from dogs, rescaled to match human osteon dimensions (Buenzli et al. 2014). The resorption rate $k_\text {res}$ (volume of bone matrix resorbed per osteoclast per day) was estimated based on the progression speed of a cortical BMU in bone and the average number of $\text {OC}_{\mathrm{a}}$ per BMU (Buenzli et al. 2014).

The $\text {OB}_{\mathrm{a}}$ and $\text {OC}_{\mathrm{a}}$ cell densities are calibrated in our model to give rise to site-specific turnover rates in steady state in agreement with experimental measurements (see Sect. 2.3). The parameters $\lambda$ and $\kappa$ related to the strength of mechanical adaptation and the parameter $\text {P}^{\text {OP}}_{\text {PTH}}$ related to the rate of bone loss during osteoporosis were calibrated in Sect. 2.6 based on experimental data.

The densities and dynamical behaviours of the precursor cell types ($\text {OB}_{\mathrm{p}}$, $\text {OB}_{\mathrm{u}}$, $\text {OC}_{\mathrm{p}}$ and $\text {OC}_{\mathrm{u}}$) and of the biochemical signalling molecules ($\text {TGF} \upbeta$, PTH, RANK, RANKL, OPG and MCSF) depend on all the remaining parameters, which represent binding reaction rates, differentiation rates, proliferation rates and apoptosis rates (see Table 1). These parameters were validated qualitatively in previous studies (Pivonka et al. 2008; Buenzli et al. 2012; Pivonka et al. 2013), i.e. they give rise to physiologically expected qualitative behaviours of the model. A quantitative calibration of these parameters to physiological values of precursor cell densities and molecule concentrations was not done due to the lack of experimental data. However, such a quantitative calibration is not an issue since it can be done without altering the time evolution of the model so long as the densities of the active cells ($\text {OC}_{\mathrm{a}}$, $\text {OB}_{\mathrm{a}}$) remain properly calibrated in Sect. 2.3. Indeed, in our simulations, both the signalling molecules and precursor cell densities ($\text {OC}_{\mathrm{u}}$, $\text {OC}_{\mathrm{p}}$, $\text {OB}_{\mathrm{u}}$, $\text {OB}_{\mathrm{p}}$) are effectively in a quasi-steady at all times, the former by model assumption (see Sect. 6.1) and the latter due to the slow evolution of bone adaptation (years) compared to the faster bone cell dynamics (days, weeks). To increase the value of density or concentration of a component by a factor $\alpha$, one may increase its production rate by the factor $\alpha$ or reduce its elimination rate by the factor $1/\alpha$. To maintain the system’s behaviour with this new value of density or concentration, it suffices to adjust the sensitivity of interactions with this component with the relevant parameters (such as binding reaction rate parameters). For example, to double the concentration of $\text {TGF} \upbeta$ without affecting how it signals $\text {OC}_{\mathrm{a}}$, $\text {OB}_{\mathrm{u}}$ and $\text {OB}_{\mathrm{p}}$, one may double $n^\text {bone}_{\text {TGF} \upbeta }$ or halve $D_{\text {TGF} \upbeta }$ (see Eq. (25)) and then double $k^{\text {TGF} \upbeta }_{\text {OC}_{\mathrm{a}}}$, $k^{\text {TGF} \upbeta }_{\text {OB}_{\mathrm{u}}}$ and $k^{\text {TGF} \upbeta }_{\text {OB}_{\mathrm{p}}}$ so that the ratios $\text {TGF} \upbeta /k^{\text {TGF} \upbeta }_{\text {OC}_{\mathrm{a}}}$, $\text {TGF} \upbeta /k^{\text {TGF} \upbeta }_{\text {OB}_{\mathrm{u}}}$ and $\text {TGF} \upbeta /k^{\text {TGF} \upbeta }_{\text {OB}_{\mathrm{p}}}$ in Eq. (23) have the same values. A similar approach was used in Buenzli et al. (2014) to calibrate the cell populations within a single BMU. This approach is also used in Sect. 6.3 below.

1.3 Recalibration of the model

Since $\text {OC}_{\mathrm{u}}$ and $\text {OB}_{\mathrm{u}}$ vary with ${f_\text {bm}}$ so as to retrieve experimentally valid turnover rates, some other parameters required modification compared with previous versions of the cell population model in which $\text {OC}_{\mathrm{u}}$ and $\text {OB}_{\mathrm{u}}$ were constant and uncalibrated (Buenzli et al. 2012; Pivonka et al. 2013).

By comparing the cell densities between this model and the previously published one (Buenzli et al. 2012), we can determine scaling coefficients which allows a systematic calibration of $\pi ^\text {act}\big (\tfrac{\text {TGF} \upbeta }{k^{\text {TGF} \upbeta }_{\text {OC}_{\mathrm{a}}}}\big )$ and $\pi ^\text {act}\big (\tfrac{\text {RANKL}}{k^{\text {RANKL}}_{\text {OC}}}\big )$. Indeed, these functions depend on the active and precursor cell densities. In the original models, the constants in these functions were calibrated such as to obtain a strong biochemical feedback response. Maintaining this strong biochemical response is the aim of this re-calibration.

The calibration is realised at ${f_\text {bm}}$ = 0.90. Both the turnover rate value and the values of $k_\text {res}$ and $\text {k}_{\mathrm{form}}$ have been changed according to the literature. Hence, by isolating $\text {OB}_{\mathrm{a}}$ and $\text {OC}_{\mathrm{a}}$ in the two new constraints of the steady state, the active osteoblast and active osteoclast read:

$$\begin{aligned} \text {OC}_{\text {a}}^{\text {new}}= & {} \frac{\chi _\text {BV}^{\text {new}}}{k_{\text {res}}^{\text {new}}} = \beta \cdot \text {OC}_{\text {a}} \end{aligned}$$

(29)

$$\begin{aligned} \text {OB}_{\text {a}}^{\text {new}}= & {} \frac{k_{\text {res}}^{\text {new}} \cdot \text {OC}_{\text {a}}^{\text {new}}}{k_{\text {form}}^{\text {new}}} = \gamma \cdot \text {OB}_{\text {a}} \end{aligned}$$

(30)

if $\delta$ is the coefficient of proportionality between the new bone turnover rate and the previous one, $\beta = \delta \cdot k_\text {res}/ k_{\text {res}}^{\text {new}}$ and $\gamma = \delta \cdot \text {k}_{\mathrm{form}}/ k_{\text {form}}^{\text {new}}$. These coefficients are introduced in the determination of $\text {TGF} \upbeta$ and $\text {OPG}$. Previously, $\text {TGF} \upbeta$ was (Buenzli et al. 2012):

$$\begin{aligned} \text {TGF} \upbeta = \frac{P^{\text {ext}}_{\text {TGF} \upbeta } + n^{\text {bone}}_{\text {TGF} \upbeta } ~k_\text {res}~\text {OC}_{\text {a}}}{D_{\text {TGF}\upbeta }} \end{aligned}$$

(31)

The new one becomes:

$$\begin{aligned} {\text {TGF} \upbeta }^{\text {new}} = \frac{P^{\text {ext}}_{\text {TGF} \upbeta } + n^{\text {bone}}_{\text {TGF}\upbeta } ~k_{\text {res}}^{\text {new}} ~\text {OC}_{\text {a}}^{\text {new}} \cdot \delta ^{-1}}{D_{\text {TGF}\upbeta }} \end{aligned}$$

(32)

The same manipulation is realised on the determination of OPG. The factor $\beta ^{\text {OPG}}_{\text {OB}_{\text {a}}}\, \text {OB}_{\text {a}}$ in Eq. (27) becomes $\beta ^{\text {OPG}}_{\text {OB}_{\text {a}}}\, \text {OB}_{\text {a}}^{\text {new}} \gamma ^{-1}$.

Appendix 2: Update frequency of mechanical state in the numerical algorithm

In our model, to simulate osteoporosis and the change of porosity with time, we need to solve the temporal equations of the bone cell populations model, Eqs. (14)–(17) and Eq. (12). Those equations via the mechanical feedback are correlated with the spatial Eqs. (2) and (42). Knowing the porosity distribution is required to determine the stress and strain distributions. Hence, we have a semi-coupled algorithm (Fig. 2).

Fig. 11

The evolution of the bone matrix volume fraction for different time steps. Note that this RVE is in the intracortical region which undergoes first resorption and then formation due to the redistribution of the mechanical loads

However, due to the separation of time scale, we can decompose the problem into two parts. Indeed, it takes more time for the microstructure to change significantly enough to influence the bone cell populations model. Therefore, we solve the bone cell populations model for a duration $\varDelta$t, assuming the mechanical feedback to be constant in this time interval. Then, we recalculate the stress and strain distributions based on the new porosity distribution, and this becomes the new mechanical feedback.

A sensitivity analysis of the solution in the time step $\varDelta$t of evolution of cell densities and bone matrix volume fraction is required. For very small time steps ($\varDelta$t $\le$ 1 day), one would expect the algorithm to converge to the exact solution. On the other hand for very large time steps ($\varDelta$t $\ge$ 5 years), a large deviation from the exact solution is expected. Figure 11 shows the evolution of the bone matrix volume fraction for one selected RVE (y = 0, z = $-10$ mm) in the cross section. These simulations show that time steps of $\varDelta$t = 250 days, 1 and 2 years lead to very similar evolution of the bone matrix volume fraction. On the other hand, $\varDelta$t = 5 years and 10 years lead to strong deviations from the smaller time increments. For all the simulations of 40 years of osteoporosis, we used a time step of 2 years.

Appendix 3: Generalised Beam theory for inhomogeneous materials

In the following, we represent the governing equations using a Cartesian (x, y, z) coordinate system. The x-axis represents the beam axis and coincides with the direction of the vascular pores (i.e. Haversian systems). The y and z coordinates describe a material point in the cross section (Fig. 1c). The origin of the system is known as the normal force centre: NC. Since our cross sections are inhomogeneous all the quantities, including the stiffness, are functions of y and z.

First, based on the constitutive relation: Hooke’s law, we determine the strain and stress relation:

$$\begin{aligned} {\varvec{\sigma }}^\text {tissue}(y, z, t) = \mathbb {C}^\text {tissue}(y, z, t) : {\varvec{\varepsilon }}^\text {tissue}(y, z, t) \end{aligned}$$

(33)

where ${\varvec{\sigma }}^\text {tissue}(y, z, t)$ and ${\varvec{\varepsilon }}^\text {tissue}(y, z, t)$ are the “tissue” stress and strain and $\mathbb {C}^\text {tissue}(y, z, t)$ the tissue stiffness matrix. The stiffness matrix is determined at the tissue scale, and the explanation is presented in Sect. 2.1.

Based on the Bernoulli hypothesis, the strain distribution appears to be a plane and remains plane even after deformation. This is why we can decompose the strain by introducing three constants: $\varepsilon _1$, $\kappa _3$ and $\kappa _2$.

$$\begin{aligned} \varepsilon ^\text {tissue}_{xx}(y,z,t) = \varepsilon _1(t) - \kappa _3(t) y + \kappa _2(t) z. \end{aligned}$$

(34)

By introducing this relation into Hooke’s law, we obtain:

$$\begin{aligned} \sigma ^\text {tissue}_{xx}(y, z, t) = \mathbb {C}^\text {tissue}_{xx}(y, z, t)~(\varepsilon _1(t) - \kappa _3(t) y + \kappa _2(t) z) \end{aligned}$$

(35)

Because we assume the shear force to be null, the stress tensor is reduced to one component: $\sigma ^\text {tissue}_{xx}(y, z, t)$. And with Bernoulli hypothesis, the strain tensor contains only one component. Hence, the stiffness matrix can be replaced by the component $\mathbb {C}^\text {tissue}_{xx}(y, z, t)$.

Here we can see that if we determine the strain constants, we would know the stress distribution. The mechanical loadings, the inputs of this model, allow us to determine the strain. Indeed, the cross section is supposed to be under a normal force: ${\varvec{N}}$ and a bending moment ${\varvec{M}}$ here divided in two bending moments: $M_y$ and $M_z$, such as $M ~\hat{{\varvec{m}}} = M_y ~\hat{{\varvec{y}}} + M_z ~\hat{{\varvec{z}}}$. By definition of the stress, we have the relations:

$$\begin{aligned}&N= \int \sigma ^\text {tissue}_{xx}(y, z, t) \text {d}A \end{aligned}$$

(36)

$$\begin{aligned}&M_y = \int z \cdot \sigma ^\text {tissue}_{xx}(y, z, t) \text {d}A \end{aligned}$$

(37)

$$\begin{aligned}&M_z = - \int y \cdot \sigma ^\text {tissue}_{xx}(y, z, t) \text {d}A \end{aligned}$$

(38)

By introducing the static moments of first and second order: EA, $ES_{y}$, $ES_{z}$, $EI_{yy}$, $EI_{zz}$, $EI_{yz}$, the equations become the following constitutive relation:

$$\begin{aligned} \begin{bmatrix} N \\ M_y \\ M_z \end{bmatrix} = \begin{bmatrix} EA&ES_{y}&ES_{z} \\ ES_{y}&EI_{yy}&EI_{yz} \\ ES_{z}&EI_{yz}&EI_{zz} \end{bmatrix} \begin{bmatrix} \varepsilon _1 \\ \kappa _2 \\ \kappa _3 \end{bmatrix} \end{aligned}$$

(39)

where:

$$\begin{aligned}&EA = \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \text {d}A\\&ES_{y} = \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \cdot y \text {d}A\\&EI_{yy} = \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \cdot y^{2} \text {d}A\\&ES_{z} = \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \cdot z \text {d}A\\&EI_{zy} = \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \cdot yz \text {d}A = EI_{yz}\\&EI_{zz} = \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \cdot z^{2} \text {d}A \end{aligned}$$

If we chose the origin of the coordinate system at the normal centre (NC) of the cross section, the coupling terms between extension and bending vanish since they become null by definition of the NC:

$$\begin{aligned} ES_{y}= & {} \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \cdot y \text {d}A = 0 \end{aligned}$$

(40)

$$\begin{aligned} ES_{z}= & {} \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \cdot z \text {d}A = 0 \end{aligned}$$

(41)

The constitutive relation can be simplified as:

$$\begin{aligned} \begin{bmatrix} N \\ M_{y} \\ M_{z} \end{bmatrix} = \begin{bmatrix} EA&0&0 \\ 0&EI_{yy}&EI_{yz} \\ 0&EI_{yz}&EI_{zz} \end{bmatrix} \begin{bmatrix} \varepsilon _1 \\ \kappa _2 \\ \kappa _3 \end{bmatrix} \end{aligned}$$

(42)

1.1 Determination of the normal force centre NC

The special location of the origin of the coordinate system for which the coupling terms ($ES_{y}$ and $ES_{z}$) between extension and bending become zero is by definition called the normal force centre NC. The result of this definition is that an axial force ${\varvec{N}}$ which acts at the NC only causes straining and no bending. The coupling terms are also referred to as weighted static moments or weighted first order moments. To find the position of the coordinate system for which the coupling terms become zero requires a tool.

Assume a temporary coordinate system: $\overline{y} - \overline{z}$ from which the porosity distribution is known. The shift in origin of this coordinate system with respect to the y–z coordinate system through the unknown NC is denoted with $\overline{y}_{NC}$ and $\overline{z}_{NC}$. The temporary coordinate system can be expressed in terms of the y–z coordinate system as:

$$\begin{aligned} \overline{y} = y + \overline{y}_{NC} ~~~~~~~~ \overline{z} = z + \overline{z}_{NC} \end{aligned}$$

Hence:

$$\begin{aligned} ES_{\overline{y}}&= \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \cdot \overline{y} \mathrm{d}A \nonumber \\&= \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \cdot y \mathrm{d}A + \overline{y}_{NC} \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \text {d}A \nonumber \\&= ES_{y} + EA \cdot \overline{y}_{NC} \end{aligned}$$

(43)

$$\begin{aligned} ES_{\overline{z}}&= \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \cdot \overline{z} \mathrm{d}A \nonumber \\&= \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \cdot z \mathrm{d}A + \overline{z}_{NC} \int \mathbb {C}^\text {tissue}_{xx}(y, z, t) \text {d}A \nonumber \\&= ES_{z} + EA \cdot \overline{z}_{NC} \end{aligned}$$

(44)

By definition, $ES_{y}$ and $ES_{z}$ with respect to the y–z coordinate system are zero from which the unknown position of the NC with respect to the known position of the $\overline{y}$– $\overline{z}$ coordinate system can be found:

$$\begin{aligned} \overline{y}_{NC}= & {} \frac{ES_{\overline{y}}}{EA} \end{aligned}$$

(45)

$$\begin{aligned} \overline{z}_{NC}= & {} \frac{ES_{\overline{z}}}{EA} \end{aligned}$$

(46)

To conclude, here is the step-by-step methodology we are using to find the stress and strain distribution in the cross section:

• 1. Localise the normal centre (NC) by computing the integrations: EA, $ES_{y}$, $ES_{z}$.

• 2. Compute the integrations: $EI_{yy}$, $EI_{zz}$ and $EI_{yz}$.

• 3. Determine the cross-sectional forces: ${\varvec{N}}$, $M_{y}$ and $M_{z}$.

• 4. Calculate the cross-sectional deformations: $\varepsilon _1$, $\kappa _2$ and $\kappa _3$ based on Eq. (42).

• 5. Find the strain distribution based on the kinematic relation, Eq. (34). Here it is important to remember to use the coordinate centred in NC.

• 6. Find the stress distribution based on Hooke’s law, Eq.  (33).

The initial cross section is extracted from a microradiograph, as it is explained in Sect. 2.4; the mechanical loading is not symmetrical. Hence, the position of the NC is changing. This is why we need to localise it after each step.