Concurrent consideration of cortical and cancellous bone within continuum bone remodelling

The kinematic description is characterised by the deformation map $\boldsymbol{\varphi}$ that maps the referential placement $\mathbf{X}$ of a physical particle in the material (reference) configuration $\mathcal{B}_{0}$ at time $t_{0}$ to its actual placement $\mathbf{x}$ in the spatial (actual) configuration $\mathcal{B}_{t}\subset\mathbb{E}^{3}$ at time $t\in\mathbb{R}_{+}$

The referential gradient $\nabla_{0}$ of the nonlinear deformation map denotes the deformation gradient $\mathbf{F}$, it maps from the material tangent space $T\mathcal{B}_{0}$ to the spatial tangent space $T\mathcal{B}_{t}$

The determinant of the deformation gradient is the Jacobian $J=\det\mathbf{F}>0$ and relates the volume elements between the two configurations. The right and left Cauchy-Green tensor $\mathbf{C}=\mathbf{F}^{T}\cdot\mathbf{F}$ and $\mathbf{b}=\mathbf{F}\cdot\mathbf{F}^{T}$ are introduced as characteristic strain measures.

In the following, the expression $\dot{\overline{\{\circ\}}}=\partial_{t}\{\circ\}|_{X}$ is used for the material time derivative of a quantity $\{\circ\}$ at fixed material placement $\mathbf{X}$. The presented formulation refers completely to the material configuration. Therefore, $\nabla_{0}\{\circ\}$ and $\mbox{Div}\{\circ\}$ denote the (referential) gradient and divergence of any field $\{\circ\}$ with respect to the material placement $\mathbf{X}$. A density quantity $\{\circ\}$ may be expressed per unit volume $\{\circ\}_{0}$ or per unit mass $\{\circ\}$, both related via the current referential nominal density $\rho_{0}$ as $\{\circ\}_{0}=\rho_{0}\{\circ\}$.

The balance of linear momentum is given by

with $\mathbf{P}=\mathbf{P}(\nabla_{0}\boldsymbol{\varphi},\rho_{0},...)$ the Piola stress tensor. Since the time scale of inertial effects is significantly smaller than that of the remodelling process, the inertial term is neglected herein (Frost 1987). The loads applied during locomotion typically correspond to a multiple of gravity. Therefore, body forces are usually not considered in bone remodelling as they are not the dominant stimulus (Jacobs et al. 1995).

The balance of mass is given by

where $R_{0}$ is the scalar-value mass source term. As the actual nominal density $\rho_{t}={\rho_{0}}/{J}$ is time-dependent also due to the coupling with the Jacobian, it is useful to capture the evolution of the referential nominal density $\rho_{0}$. The relative change in (referential) nominal density $\tilde{\rho_{0}}$ is defined in terms of the initial nominal density $\rho_{0}^{\star}$ as

Therefore, density gain results for values $\tilde{\rho_{0}}>0$ whereas density loss is obtained for $-1<\tilde{\rho_{0}}<0$.

We distinguish between the homogenised macro scale with nominal (time-evolving, referential) density $\rho_{0}=\rho_{0}(t)\in(0,\rho_{0}^{s}]$ occupying the (time-independent, referential) volume element $\,\mbox{d}V$ and the heterogeneous micro scale with true (time-independent, referential) density $\rho_{0}^{s}$ of the solid skeleton material occupying the (time-evolving, referential) partial volume element $\,\mbox{d}V^{s}=\,\mbox{d}V^{s}(t)\in(0,\,\mbox{d}V]$, see Fig. 1. The (time-evolving, referential) volume fraction of solid skeleton material thus follows as

The (referential) volume fraction $\nu_{0}(t)$ indicates how much space within a (referential) volume element $\,\mbox{d}V$ at the macro scale is occupied by solid skeleton material, thus $1-\nu_{0}(t)$ denotes the (referential) volume fraction of the (empty) pore space.

At time $t=0$ the corresponding initial values for the (referential) volume fraction $\nu_{0}(0)$, the (referential) partial volume element $\,\mbox{d}V^{s}(0)$ and the nominal (referential) density $\rho_{0}(0)$ are related via

and are denoted the initial (referential) volume fraction $\nu_{0}^{\star}:=\nu_{0}(0)$, the initial (referential) partial volume $\,\mbox{d}V^{s\star}:=\,\mbox{d}V^{s}(0)$ and the initial nominal (referential) density $\rho_{0}^{\star}:=\rho_{0}(0)$, respectively. Typically, $\nu_{0}^{\star}$ is significantly higher for cortical bone as compared to cancellous bone.

Furthermore, we may distinguish between the nominal relative (referential) density and its true counterpart

when phenomenologically modelling the homogenised mechanical properties of bone based on its skeleton structure at the micro scale, i.e. its microstructure. In the following, the present contribution aims to rationalise various modelling options stemming from this distinction.

Bone possesses a microstructure characterised by the spatially heterogeneous distribution of solid skeleton and pore space. The resulting structural bone properties of the heterogeneous solid skeleton at the micro scale may then be homogenised into nominal bone properties of the effective material at the macro scale.

A simplistic, however efficient and effective homogenisation rule – typically employed for cellular materials – relates the true material parameters of the solid skeleton at the micro scale to the corresponding effective, nominal parameters of the homogenised material at the macro scale via a power $N$ of the initial (referential) volume fraction $\nu_{0}^{\star}$. Let for example $\lambda^{s}$ and $\mu^{s}$ denote the Lamé parameters of the solid skeleton material. Then, the homogenised, effective bone material is characterised by the nominal material parameters

Typically, the concrete values assumed for the homogenisation exponent $N$ capture different microstructures and homogenisation approaches. In this contribution, either the most elementary volumetric mixture rule with $N=1$ or a relation with $N=n$ shall be considered. The latter is particularly motivated by the proposal of (Gibson 2005) postulating a dependence of the nominal Young’s modulus of a homogenised (linear elastic) porous material on the Young’s modulus of the solid skeleton material via a power of the true relative density.

In what follows we assume the concept of strain equivalence as also employed in damage mechanics. To elucidate this aspect more concretely, we employ, without loss of generality, a (referential) volume-specific Neo-Hookean true free energy density for the solid skeleton material

Note that strain equivalence postulates identical ’strain’ at the macro and the micro scale. Thus, based on the nominal material parameters, the corresponding (referential) volume- specific initial nominal free energy follows from its true counterpart as

Concretely, the (referential) volume-specific Neo-Hookean initial nominal free energy reads in terms of the nominal Lamé parameters as

The initial microstructure evolves with time due to remodelling into the current microstructure, thus necessitating to express the (referential) volume-specific current nominal free energy density $\Psi_{0}$ in terms of its initial counterpart.

To this end, again motivated by (Gibson 2005), when taking the microstructure of cellular materials into account, the (referential) volume-specific current nominal free energy density $\Psi_{0}$ is determined by weighting its initial nominal counterpart $\Psi_{0}^{\star}$ with a power $n$ of the nominal relative density

Then the Piola stress entering the equilibrium equation derives from the (referential) volume-specific current nominal free energy density $\Psi_{0}$ as

and expands concretely for the Neo-Hookean case with the initial Piola stress

Incorporating the (referential) volume-specific Neo-Hookean true free energy density $\Psi_{0}^{s}$ eventually renders $\Psi_{0}$ in terms of the power $n$ of the true relative density and the power $N-n$ of the initial (referential) volume fraction

The exponent $n$ is typically set to $n=2$ for trabecular bone with low nominal density and to $n>2$ for denser cancellous bone, compare to (Rho et al. 1993) and (Gibson 2005).

Setting finally the homogenisation exponent $N$ in particular to $N=n$ results in an expression for the (referential) volume-specific current nominal free energy density $\Psi_{0}$ in terms of its true counterpart and the true relative density

Correspondingly, the Piola stress expands in this case as

Setting $N=n$ tacitly assumes that the geometry of the microstructure is not distorted by the deformation (for example a circular pore remains a circular pore). In contrast, setting $N\not=n$ conceptually allows to distinguish between the initial and the (deformation-dependent) current microstructure, thus in general $n=n(\mathbf{F})$ with $N=n(\mathbf{I})$.

This contribution shall explore the consequences of two limiting cases by setting $N$ to either $N=1$ (volumetric mixture rule) or to $N=n$ (cellular material like homogenisation rule). More sophisticated scenarios with $n=n(\mathbf{F})$ where $N=n(\mathbf{I})$ shall not be analysed here.

To model changes in the nominal (referential) density $\rho_{0}$ incurred by mechanical stimuli, the volume-localised mass balance $\dot{\rho}_{0}=R_{0}$ is equipped with a (referential) source term $R_{0}$. Following (Harrigan and Hamilton 1993), the (referential) mass source $R_{0}$ is expressed by the (referential) volume-specific current nominal free energy density $\Psi_{0}$ and the nominal relative density as

Here, the parameter $c$ controls the speed of density evolution and has dimension of time divided by length squared. The attractor stimulus $\Psi_{0}^{a}$ defines the biological equilibrium (homeostasis) to which the system is driving at constant (mechanical) stimulus. The exponent $m$ is a material parameter particularising the mechanical stimulus (Harrigan and Hamilton 1993) in terms of the nominal relative density.

The (referential) mass source may alternatively be expressed by the (referential) volume-specific initial nominal free energy density $\Psi_{0}^{\star}$ or, likewise, by its true counterpart $\Psi_{0}^{s}$ as

We may deduce from this representation that the exponent $m$ is also crucial for the stability of the density evolution: at homeostasis ($R_{0}=0$) $\rho_{0}>\rho_{0}^{\star}$ for $\Psi_{0}^{\star}>\Psi_{0}^{a}$ only if $m-n>0$. In conclusion, the exponent $m$ should be chosen as $m>n$ (Harrigan and Hamilton 1994), (Kuhl and Steinmann 2003).

Setting finally the homogenisation exponent $N$ in particular to $N=n$ results in alternative expressions for the (referential) mass source

To differentiate between cancellous and cortical bone tissue, the effect of the initial volume fraction shall be examined. For this purpose, a simple unit square sample is studied, which is loaded in tension by a force function that increases stepwise as illustrated in Figure 2.

The reference parameters are chosen similar to (Kuhl and Steinmann 2003): true elastic modulus $E^{s}=1.0$, Poisson’s ratio $\nu^{s}=0.0$, attractor stimulus $\Psi_{0}^{a}=1.0$, parameter $m=3$ and parameter $n=2$. The homogenisation exponent is chosen as $N=n$ as motivated by the proposal of (Gibson 2005).

A true density of the solid material $\rho_{0}^{s}=2$, similar to (Gibson 2005), is assumed for both bone microstructure types in the following. It is known from measurements that the true density of trabeculae is slightly smaller and the one of cortical bone fibrils slightly larger, see (An and Draughn 2000). However, as the objective of this analysis is to point out the general trend in the evolution of the different tissue microstructures, the assumption of equal true density is sufficient in this context. An initial nominal density of $\rho_{0}^{\star}=0.5$ leading to an initial volume fraction of $\nu_{0}^{\star}={\rho_{0}^{\star}}/{\rho_{0}^{s}}={0.5}/{2}=0.25$ is chosen as representative of cancellous bone. A comparatively high value of $\nu_{0}^{\star}={\rho_{0}^{\star}}/{\rho_{0}^{s}}={1.8}/{2}=0.9$ is set for the volume fraction representing compact bone.

As shown in Figure 3, a small initial volume fraction (gray line) leads to rapid bone growth while a high value (black line) causes no density change in the first load step. Only after increasing the load in the second step, and thus increasing the mechanical stimulus, an adaptation of bone density can be observed.

Cortical bone in the model therefore requires a significantly higher stimulus and grows less in relative terms. But despite its rapid growth, the absolute density of cancellous bone still remains at a lower level.

Due to its small surface area in relation to the volume, cortical bone adapts to the load conditions much slower and less than spongy bone. Considering the mass source term in Equation 19, it appears that there are three parameters which can influence the behaviour of bone adaptation. By reducing the parameter $c$, the speed of density evolution can be slowed down, but after a sufficiently long time the result always converges towards the same saturation level. Therefore, a parameter adjustment of the exponent $m$ and the attractor stimulus $\Psi_{0}^{a}$ shall be considered for cortical bone with $\nu_{0}^{\star}=0.9$ and $N=n=2$.

In order to satisfy the stability criterion, the exponent $m$ must be selected greater than $n$, see (Harrigan and Hamilton 1993). Thus, $m=3$, $m=2.5$ and $m=10$ are investigated as depicted in Figure 4.

As can be seen, reducing the exponent leads to an increase in density growth (dashed line). In contrast, choosing a large value for the exponent $m$ tends to slow down the density evolution (dotted line).

One possibility to simulate the significantly reduced density evolution of compact bone compared to spongy bone would therefore be setting a relatively high exponent $m$.

A variation of the attractor stimulus is analysed in Figure 5.

The dotted line in Figure 5 based on a relatively high attractor stimulus displays low density results. Bone growth can not be achieved until the third load level. Reducing the attractor stimulus (dashed line) results in an overall increased gain in nominal density.

Since cortical bone should neither grow nor absorb faster than cancellous bone, in the following, the same attractor stimulus will be assumed for both microstructure types.

To analyse the local distribution of cortical and cancellous bone and their interaction, a two-dimensional proximal femur head is examined.

Three typical load cases representing the midstance phase of gait, extreme abduction and adduction are applied to the model as reported in (Carter and Beauprè 2010). In the calculations bi-linear element expansions with a 2x2 Gauss quadrature rule are used. The mesh, see Figure 6 (left), consists of 13317 elements and 13602 nodes. The parameters are set as follows: Poisson’s ratio $\nu^{s}=0.2$, attractor stimulus $\Psi_{0}^{a}=0.01$, exponent $n=2$ and true solid material density $\rho_{0}^{s}=2.0$.

An initial nominal cortical density of $\rho_{0}^{\star}=1.825$, leading to an initial volume fraction of $\nu_{0}^{\star}={\rho_{0}^{\star}}/{\rho_{0}^{s}}={1.825}/{2}=0.9125$ is set for cortical bone (elements at the external boundary), and $\rho_{0}^{\star}=0.55$ with $\nu_{0}^{\star}={\rho_{0}^{\star}}/{\rho_{0}^{s}}={0.55}/{2}=0.275$ for cancellous bone. To avoid an abrupt unrealistic change between elements of different microstructure types, a small transition zone between both regions is created using a sigmoid function, see Figure 6 (middle).

Four different models are studied and compared in the following, all using the homogenisation exponent $N=n=2$. Selected parameters and thus the difference between the four models are illustrated in Table 3.3.

The first model represents the full spongy bone as reported in previous literature with an initial nominal density of $\rho_{0}^{\star}=0.55$. The second model differentiates between cortical and cancellous bone using the initial volume fraction or rather the initial nominal density distribution as depicted in Figure 6 (middle).

With a true elastic modulus of $E^{s}=10000$, which is similarly documented for solid trabecular material by (Rho et al. 1993) and (Yamada et al. 2014), and based on the quadratic relation $N=n=2$ between stiffness and nominal density according to (Gibson 2005), this results in a nominal stiffness of $E^{\star}=756.25$ in the spongy region at an initial nominal density of $\rho_{0}^{\star}=0.55$, which roughly reflects the magnitude of stiffness as reported in (Wirtz et al. 2000) and which is also utilised in the full spongy model. The maximum nominal elastic modulus of cortical bone is approximately $E^{\star}=8327$ and thus significantly higher, see Figure 6 (right).

Motivated by the results of section 3.2, the third and fourth model include a variation of the exponent $m$. While the exponent is always set to $m=3$ in the first and second model, the others include an increase within the cortical bone area to a maximum of $m=10$ in the third and $m=18$ in the fourth model. The distribution of the exponent $m$ follows the respective volume fraction which indicates the type of bone microstructure. A slight transition is implemented between the two bone types using an hyperbolic tangent function, see Figure 7.

The results for the evolution of the nominal density are displayed for each model at different points in time in Figure 8.

Comparing the first two rows, it becomes clear that the pattern of density evolution has changed remarkably. The density growth in the femoral head and neck is significantly less intense when taking cortical bone into consideration. Density loss increases and slightly shifts in the intertrochanter region, which leads to an obvious change in density pattern. In areas with subliminal stress like the medial femur head, cortical bone decomposes slightly and at the same time ensures a less distinct density reduction in the adjacent spongy region.

Increasing the exponent $m$ in cortical bone (see row 3,4) suppresses its resorption and cortical bone remains, although weaker and thinner. Density loss of cancellous bone in the intertrochanter region is reduced and an increase in density in the femoral neck, albeit small, can be observed compared to the second model.

In summary, it appears that cortical bone acts as a kind of stiffer shell. Therefore, the interior of the femur head bears less load and cancellous bone does not adapt its density as much as in the first, full spongy model.

Finally, the alternative volumetric mixture rule with $N=1$ shall be examined. We first re-analyse the simple square example as described in section 3.1 with a variation in the initial volume fraction and homogenisation exponent $N=1$, whereby all the other parameters remain unchanged. The temporal evolution of the nominal density is shown in Figure 9.

Comparing the results of $N=1$ with those of section 3.1 in Figure 3, it can be seen that density growth is reduced for an initial volume fraction of $\nu_{0}^{\star}=0.25$ (gray lines), whereas a reduction in density growth for an initial volume fraction of $\nu_{0}^{\star}=0.9$, which even leads to a small density loss in the first load step, is hardly visible (black lines).

The reason for this difference lies in the modified relationship between stiffness and nominal density from quadratic to linear. Whereas the initial nominal stiffness has quadrupled ($1/0.25$) at an initial volume fraction of $\nu_{0}^{\star}=0.25$, the initial volume fraction of $\nu_{0}^{\star}=0.9$ yields only a factor of $1.111=1/0.9$.

Keeping this in mind, we again consider the proximal femur model, similar to section 3.3. Setting the true elastic modulus to $E^{s}=10000$ yields an initial nominal stiffness of $E^{\star}=2750$ for spongy bone with $\nu_{0}^{\star}=0.275$ via the linear relationship, which is far too high. Reducing the Young’s modulus of the solid material to $E^{s}=2750$ results in the same $E^{\star}$ for cancellous bone as in section 3.3, see Figure 10. However, the stiffness for cortical bone becomes much too small.

Nevertheless, Figure 11 shows the results of the four different models including the full spongy model, the distinction between cancellous and cortical bone with $m=3$, and two variations of the exponent $m$ for cortical bone, all calculated with a true elastic modulus of $E^{s}=2750$ and an homogenisation exponent of $N=1$. Since the true Young’s modulus has been properly adjusted and thus the nominal stiffness as well as all the other parameters in the full spongy model correspond to those of section 3.3, the result is the same as can be seen in the first rows of Figure 11 and 8.

When comparing the different homogenisation exponents for the second model, we note the significantly higher density growth in the area of the femoral neck, but also an increased bone loss in the area of the medial femoral head. Furthermore, the pattern of density evolution more closely resembles that of the full spongy model, especially in the area of the greater trochanter. This effect is enhanced when increasing the exponent $m$ for cortical bone, see Figure 11 row 3 and 4. In this scenario, cortical bone decomposes much more in subliminally loaded regions with $N=1$, especially with a small exponent $m$.

Overall, the application of the volumetric mixture rule with a homogenisation exponent of $N=1$ leads to a nominal density distribution that is more similar to the full-spongy model. However, ensuring a stiffness-density ratio for cancellous bone as motivated by literature will not simultaneously achieve a realistic ratio in the cortical area. A further possible option would therefore be a separate definition of an individual true Young’s modulus for cortical bone, but will not be examined in this survey.