Bone Growth and Remodeling

Overview

In FEBio the theory of bone growth and remodeling is inspired by the work of Huiskes, Weinans and Mullender. It uses the framework of reactive multiphasic theory. In this framework, trabecular bone is modeled as a porous solid whose referential apparent density $\rho_r^s$ (mass of bone in current configuration, per mixture volume in reference configuration),

$\rho_r^s = \frac{dm^s}{dV_r}$

evolves over time due to a reaction that adds or removes bone mass in response to loading conditions. The units of $\rho_r^s$ are [mass]/[volume]. The reaction may involve any number of constituents, but the simplest embodiement of bone growth and remodeling only models the solid bone mass explicitly, whereas other constituents, such as cells that synthesize or resorb bone, and soluble nutrients or waste products, are implicitly assumed to exist in the mixture.

$\text{cells} + \text{nutrients} \stackrel{\text{loading}}{\to} \text{cells} + \text{solid bone}$

The solid bone constituent is modeled using a solid-bound molecule (sbm). The above reaction is modeled using using a mass-action-forward reaction with a Huiskes reaction rate, to produce the mass balance equation

$\hat{\rho}_{r}^{s}=\frac{d\rho_{r}^{s}}{dt}=B\left(\frac{\Psi_{r}}{\rho_{r}^{s}}-\psi_{0}\right)$

In this expression, $\hat{\rho}_{r}^{s}$ is the mass density supply resulting from the reactive process (units of [mass]/[volume][time]), $\Psi_{r}$ is the strain energy density of the trabecular bone (units of [energy]/[volume]), and $\psi_{0}$ is the specific strain energy threshold above which bone mass is added, and below which bone mass is resorbed (units of [energy]/[mass]). The parameter $B$ represents the rate at which the reaction proceeds, in units of [mass][time]/[length]$^5$.

The strain energy density $\Psi_{r}$ of trabecular bone may be calculated by assuming that this material is isotropic linear elastic. Thus,

$\Psi_{r}=\frac{E_{Y}}{2\left(1+\nu\right)}\left[\frac{\nu}{1-2\nu}\left(\text{tr}\boldsymbol{\varepsilon}\right)^{2}+\text{tr}\boldsymbol{\varepsilon}^{2}\right]$

where $E_Y$ is Young’s modulus, $\nu$ is Poisson’s ratio, and $\boldsymbol{\varepsilon}$ is the infinitesimal strain tensor. Based on the work of Currey and Rice et al., it is assumed that $E_Y$ varies with $\rho_r^s$ according to

$E_{Y}=E_{0}\left(\frac{\rho_{r}^{s}}{\rho_{0}}\right)^{\gamma}$

where $E_0$ is the modulus in the limit when $\rho_{0}=\rho_{r}^{s}$, and $\gamma$ is a power exponent derived from experimental measurements. In FEBio this constitutive model is called Carter-Hayes due to their similar model for bone strength.

Multiphasic models normally expect to solve for the interstitial fluid pressure within the mixture. However, the bone remodeling problems illustrated in this case study are performed assuming that the fluid pressure is negligible. Therefore, boundary conditions are prescribed that fix all the pressure degrees of freedom (setting them to zero) to avoid performing extraneous calculations. Therefore, the values of the hydraulic permeability and osmotic coefficient adopted in the examples below are arbitrary.

The models described in the two examples below can be downloaded from the FEBioStudio Repository.

Example: 1D Bar Remodeling

Consider a one-dimensional bar along the $x-$direction ($0 \le x \le h$), with uniform cross section, subjected to a uniform body force $f_0$ along $x$ (units of [force]/[volume]), and a prescribed normal traction $\sigma_0$ at the rightmost end ($x = h $).

Because of the 1D nature of this analysis we may set Poisson’s ratio to $\nu=0$, implying that $\sigma_{xx}=E_Y \varepsilon_{xx}$. The normal stress $\sigma_{xx}$ along the bar reduces to $\sigma_{xx}=\sigma_0 + f_0 \left(h – x\right)$. The resulting strain energy density may be evaluated as $\Psi_r=\frac{1}{2}\sigma_{xx}\varepsilon_{xx}=\frac{\sigma_{xx}^2}{2 E_Y}$. At steady state (when $\frac{d\rho_{r}^{s}}{dt}=0$), the solution to this remodeling problem can be obtained analytically from $\frac{\Psi_{r}}{\rho_{r}^{s}}=\psi_{0}$, thus

$\rho_{r}^{s}=\left(\frac{\sigma_{xx}^{2}\rho_{0}^{\gamma}}{2E_{0}\psi_{0}}\right)^{\frac{1}{1+\gamma}}$

We solve this problem in FEBioStudio using $h=10$, $\sigma_0=-10$, $f_0=-1$, $E_0=10^4$, $\rho_0=1$, $\gamma=2$, $B=1$ and $\psi_0=0.01$.

Mesh Geometry and Boundary Conditions

1. Create a rectangular block 10 x 1 x 1
2. Mesh uniformly with 100 x 1 x 1 elements
3. Zero X,Y,Z displacements at left face ($x=0$)
4. Prescribe normal Traction {-10,0,0} on right face ($x=10$)
5. Prescribe tangential Traction {-0.25,0,0} on four lateral faces
6. Prescribe zero fluid pressure on four lateral faces

Material Definition

Solid-bound molecules can alter the referential solid volume fraction $\phi_r^s$ of a mixture according to the relation outlined in the User’s Manual. The formula in this case reduces to $\phi_r^s = \phi_0^s + \frac{\rho_r^s}{\rho_T^s}$, where $\rho_T^s$ is the true density of the remodeling solid and $\phi_0^s$ is the initial solid volume fraction contributed by other solid constituents not modeled explicitly. In this example we select $\phi_0^s = 0$ and $\rho_T^s=1000$, while we limit the range of $\rho_r^s$ to $0.01 \le \rho_r^s \le 2$, so that the remodeling process has a negligible influence on $\phi_r^s$.

1. Add a solid-bound molecule (name: bone) to the SBM table (Physics menu) with charge = 0, (true) density=1000, molar mass = 1
2. Create multiphasic material with
   • Name: Remodeling Bar
   • Solid volume fraction = 0
   • Carter-Hayes solid (E0=10000, rho0=1, gamma=2, sbm=bone
   • Constant permeability (perm=1)
   • Constant osmotic coefficient (osmotic coefficient = 1)
   • Add SBM bone to list of solid-bound molecules, set initial density=1, min density=0.01, max density=2)
3. Assign material to bar

Chemical Reaction

1. Add a chemical reaction (Physics menu)
   • Name: Solid remodeling
   • Multiphasic material: Remodeling Bar
   • Reaction: mass-action-forward
   • Fwd Rate: Huiskes reaction rate
   • Check ‘Override Calculated Vbar’ box
   • Products: Add SBM bone
2. Edit material Remodeling Bar
   • B=1, psi0=0.01, vP=1

Analysis

1. Add analysis step (Multiphasic)
2. Time stepping:
   • Time steps = 200
   • Step size = 1
   • Max step size = 1
3. Quasi-Newton method: Full Newton
4. Output plotfile requests
   • Keep defaults
   • Add referential solid volume fraction, sbm referential apparent density, strain energy density, specific strain energy
5. Globals: Set absolute temperature to non-zero value (e.g., 298 K)

Results

Plotting the specific strain energy over time shows that it evolves toward the uniform value of $\psi_0=0.01$.

A plot of $\rho_r^s$ along $x$ shows the distribution as it evolves over time.

A comparison of the solution at $t=200$ with the analytical solution given above shows an error of $0.13\%$ or less.

Example: 2D Plate Remodeling

This example follows a nearly identical approach as the 1D bar remodeling example. Create a plate using a rectangular box (10 x 10 x 1). Mesh it uniformly with 50 x 50 x 1 elements. Set zero displacements along X, Y and Z at the bottom face, and zero displacements along Z on the front and back faces (plane strain analysis). Prescribe pressure on the top face and use the <math> input format to produce a linearly increasing normal compressive traction from the left edge to the right edge.

Set the remodeling rate parameter B to 10, the reference specific strain energy psi0 to 0.0005, and run the analysis for 5000 steps, using a plot stride of 10 to avoid saving every time point to the plot file. Switch the quasi-Newton method to Broyden, with 50 maximum updates, 5 maximum reformations, and 50 optimal iterations.

Results show that a trabecular-like pattern emerges over time. However, because of the constraint on the range of $\rho_r^s$, the steady-state solution does not produce a uniform value of the specific strain energy.

Was this article helpful?

0 out of 5 stars

5 Stars		0%
4 Stars		0%
3 Stars		0%
2 Stars		0%
1 Stars		0%

How can we improve this article?

Please submit the reason for your vote so that we can improve the article.

Name

Email

Submit