Does mechanical stimulation really protect the architecture of trabecular bone? A simulation study

Abstract

Although it is beyond doubt that mechanical stimulation is crucial to maintain bone mass, its role in preserving bone architecture is much less clear. Commonly, it is assumed that mechanics helps to conserve the trabecular network since an “accidental” thinning of a trabecula due to a resorption event would result in a local increase of load, thereby activating bone deposition there. However, considering that the thin trabecula is part of a network, it is not evident that load concentration happens locally on the weakened trabecula. The aim of this work was to clarify whether mechanical load has a protective role for preserving the trabecular network during remodeling. Trabecular bone is made dynamic by a remodeling algorithm, which results in a thickening/thinning of trabeculae with high/low strain energy density. Our simulations show that larger deviations from a regular cubic lattice result in a greater loss of trabeculae. Around lost trabeculae, the remaining trabeculae are on average thinner. More generally, thin trabeculae are more likely to have thin trabeculae in their neighborhood. The plausible consideration that a thin trabecula concentrates a higher amount of strain energy within itself is therefore only true when considering a single isolated trabecula. Mechano-regulated remodeling within a network-like architecture leads to local concentrations of thin trabeculae.

Appendices

Appendix A

Relationship between trabecular thickness and strain energy density frequency distributions

With $T$ the thickness of a single trabecula having a Young’s modulus $E$ and subjected to an axial load $F$, the strain energy density SED can be written as:

$$\begin{aligned} \hbox {SED}(T)=g\left( T \right) =\frac{8}{\pi ^{2}}\frac{F^{2}}{E}\left( {\frac{1}{T}} \right) ^{4}. \end{aligned}$$

(A1)

Considering a collection of independent trabeculae and assuming that $T$ has a Gaussian probability density function $n_T \left( T \right)$ (blue curve, Fig. 2c) with mean value $\mu$ and standard deviation $\sigma$:

$$\begin{aligned} n_T \left( T \right) =\frac{1}{\sigma \sqrt{2\pi }}\hbox {exp}\left( {-\frac{\left( {T-\mu } \right) ^{2}}{2\sigma ^{2}}} \right) , \end{aligned}$$

(A2)

the probability density function of SED is calculated as (blue curve, Fig. 2d):

$$\begin{aligned}&\!\!\!f_S \left( {\hbox {SED}} \right) =n_T \left( {g^{-1}\left( {\hbox {SED}} \right) } \right) \left| {\frac{dg^{-1}\left( {\hbox {SED}} \right) }{d(\hbox {SED})}} \right| \nonumber \\&\quad =\frac{\gamma ^{1/4}}{4\sigma \sqrt{2\pi }}\left( {\frac{1}{\hbox {SED}}} \right) ^{5/4}\hbox {exp}\left\{ {-\frac{\left[ {\left( {\frac{\gamma }{\hbox {SED}}} \right) ^{1/4}-\mu } \right] ^{2}}{2\sigma ^{2}}} \right\} ,\nonumber \\ \end{aligned}$$

(A3)

where $\gamma =\frac{8F^{2}}{\pi ^{2}E}$ .

Appendix B

1.1 Mechanical control of the remodeling of a single trabecula: recurrence relation

Considering one single trabecula loaded by a constant axial force $F$ and characterized by a Young’s modulus $E$, the relationship between the cross-sectional area $A$ and the strain energy density SED at a discrete time point $i$ is:

$$\begin{aligned} \hbox {SED}_i =\frac{F^{2}}{2E}\left( {\frac{1}{A_i }} \right) ^{2} . \end{aligned}$$

(B1)

According to the remodeling rule introduced in Fig. 2a and assuming linearity (i.e., SED remains smaller than 2 ${\hbox {SED}}_\mathrm{ref})$, the change in cross-sectional area $\Delta A$ is given by:

$$\begin{aligned} \Delta A=\Delta A_\mathrm{max} \left( {\frac{\hbox {SED}}{\hbox {SED}_\mathrm{ref} }-1} \right) . \end{aligned}$$

(B2)

Hence, by inserting (B1) into (B2), the recurrence relation for the cross-sectional area reads:

$$\begin{aligned} A_{i+1} =A_i +\Delta A=A_i +\Delta A_\mathrm{max} \left( {\frac{F^{2}}{2E\,\mathrm{SED}_\mathrm{ref} }\frac{1}{A_i^2 }-1} \right) .\nonumber \\ \end{aligned}$$

(B3)

The additional normalization by the factor $\Delta A_\mathrm{max}$ (for its definition see Method section) allows writing (B3) in a dimensionless form:

$$\begin{aligned} x_{i+1} =x_i +\hat{{\gamma }}\frac{1}{x_i^2 }-1, \end{aligned}$$

(B4)

with $x_i =\frac{A_i }{\Delta A_\mathrm{max} }$ and $\hat{{\gamma }}=\frac{F^{2}}{2E\,\mathrm{SED}_\mathrm{ref} }\frac{1}{\Delta A_\mathrm{max}^2 }$. The fixed points of a general recurrence relation $x_{i+1} =f(x_i )$ are found by setting $f\left( {x^{*}} \right) =x^{*}$ (Strogatz 2001), hence:

$$\begin{aligned} x^{*}=x^{*}+\hat{{\gamma }}\frac{1}{x^{*2}}-1, \end{aligned}$$

(B5)

resulting in

$$\begin{aligned} x^{*}=\pm \sqrt{\hat{{\gamma }}}. \end{aligned}$$

(B6)

The stability of the fix points (B6) can be then analyzed by looking at the first derivative of the recurrence relation (B4):

$$\begin{aligned} \lambda =f^{{\prime }}\left( {x^{*}} \right) =1-\frac{2}{\sqrt{\hat{{\gamma }}}}. \end{aligned}$$

(B7)

The model is unstable when ${\vert }{\lambda }{\vert }>1$, hence for $0<\hat{{\gamma }}<1$. Assuming $F$= 1.5 N, $E$= 10 GPa and ${\Delta A}_\mathrm{max} = 0.0039\hbox { mm}^{2}$, the corresponding value of ${\hbox {SED}}_\mathrm{ref}$ above which the recurrence relation gives rise to oscillations in the cross-sectional area is $7.4 \times 10^{6}\hbox { J/m}^{3}$ (Fig. 7).

Fig. 7

Time evolution of the cross-sectional area of an individual trabecula having an initial thickness of 0.14 mm and remodeled according to the recurrence relation derived in Appendix 2 (Eq. (B7)) for three different values of ${\hbox {SED}}_\mathrm{ref}$. When ${\hbox {SED}}_\mathrm{ref}$ is well below the instability threshold (e.g., ${\hbox {SED}}_\mathrm{ref }=5 \times 10^{6}\, \hbox {J/m}^{3} < 7.4 \times 10^{6}\hbox { J/m}^{3})$, the trabecula can be controlled by a mechano-regulated remodeling process. If ${\hbox {SED}}_\mathrm{ref}$ approaches the instability threshold (e.g., ${\hbox {SED}}_\mathrm{ref }= 7 \times 10^{6}\hbox { J/m}^{3})$, despite some initial oscillations, the cross-sectional area still converges to a steady state. For values of ${\hbox {SED}}_\mathrm{ref}$ above the instability threshold (e.g., ${\hbox {SED}}_\mathrm{ref }= 9 \times 10^{6} \hbox { J/m}^{3})$, the cross-sectional area oscillates without approaching a steady state