Simulating the Adaptive Remodeling of Bone in Response to Mechanical Influences using ANSYS

19 Pages·2004·2.35 MB·English

by  Uwe Winter

Simulating the Adaptive Remodeling of Bone in
Response to Mechanical Influences using ANSYS
Uwe Winter
Kx Simulation Technologies, Inc., Cincinnati, OH
Abstract
Bone remodeling refers to the self-optimizing behavior of bone tissue in response to a mechanical stimulus,
a natural phenomenon producing strong structures with minimal weight. This paper describes BORA, an
iterative finite element method incorporated into ANSYS which simulates the biological mechanism of
bone remodeling. The equations relating adaptive changes in bone density to the internal stresses and
strains are presented together with the methods used for implementing the density evolution equations into
ANSYS. Advanced features such as an adaptive time stepping scheme, solution tracking, and results
visualization utilities were developed with the goal of making BORA an efficient engineering tool that is
convenient to use.
Introduction
Living bone continually reshapes itself as it adapts to varying loading conditions. Changes in the stress and
strain distribution in bone tissue are known to stimulate cell activity, causing bone tissue to increase its
density in response to high load levels, while decreased load levels will trigger the resorption of bone
material. Examples of such changes in bone density are the weakening of bones as observed after
prolonged immobilization or exposure to a reduced gravitational field, whereas athletes will experience an
increase in bone mass and stronger bones as a result of intensive exercise programs.
The process which regulates this relationship between mechanical stimuli and bone density is called “bone
remodeling” and was first noted by Roux, who termed it “direct functional self-shaping” more than 100
years ago.
Understanding and being able to predict this phenomenon is of particular interest to the orthopedic
industry. Whenever an implant is inserted into the body, existing bone has to be removed for the implant to
take its place. This alters the load path and the strain distribution for the bone tissue in the vicinity of the
implant, causing a redistribution of bone mass at the implant-bone interface. While the deposition of
higher density bone material near the implant is desirable for good fixation, localized bone resorption may
result in the loosening of the device, resulting in pain for the patient.
Ongoing research efforts by some of the key contributors in the field of bone remodeling (refer to the
References section), have provided us with various mathematical models for relating bone density
distribution to the internal loads in the bone tissue, with the earliest bone remodeling laws presented by
Wolff and dating back to 1892. While contemporary researchers still disagree on whether stress, or strain,
or strain energy density provides the stimulus for bone remodeling, and to what extent load history and
initial conditions affect the response, they all agree on the basic premise that bone has mechanoreceptors
and that the trabecular architecture adapts itself to mechanical influences.
This paper discusses BORA (BOne Remodeling Algorithm), a finite element (FE) method based program
developed to simulate the density evolution of bone tissue and predict bone morphology with ANSYS.
BORA is the intermediate result of a multi-phase development project commissioned by DePuy
Orthopaedics (a Johnson & Johnson company based in Warsaw, Indiana) as part of their continuing drive to
develop new engineering and simulation technologies and integrate them into their existing suite of CAE
tools. It was the goal of the initial phases of this project to implement the theory of bone remodeling into
ANSYS and provide the capability of predicting how bone responds to the altered load distribution after
arthroplasty. Ultimately, BORA will allow engineers to compare different prosthesis designs and, given a
patient’s lifestyle, medical history, and specific bone structure, determine the ideal implant which leads to a
minimal amount of bone resorption and optimal trabecular bone adaptation.
The density evolution equations selected for implementation into the BORA program are based on the bone
remodeling concepts and mathematical equations presented by Weinans et al. (1992), who chose the strain
energy density as the source of the mechanical stimulus.
Theory
It is assumed that bone material has sensors, which can detect a mechanical stimulus, and, depending on the
magnitude of the stimulus, cause local bone adaptations, i.e. a change in density over a specific time
increment. Furthermore, it has been shown that there is a relationship between density and the mechanical
properties of bone.
Given a mathematical representation of these two basic concepts, the resulting bone remodeling process is
readily implemented into ANSYS, using the apparent density (porosity) as the characterization of the
internal morphology.
Density Evolution Equations
Assuming that the mechanical stimulus at a particular location (x,y,z) can be determined from the stress
tensor σ, the strain tensor ε, and the apparent density ρ, we can write the rate of change of the apparent
density as an objective function
dρ
=F(σ,ε,ρ) , (1)
dt
where the density ranges between zero and a maximum value ρ for cortical bone. The system is
cb
considered in equilibrium when the objective function F reaches zero. A generic objective function,
specified in terms of discrete locations (finite elements) and discrete time steps as in the case of linear
iterative FE methods, is proposed with equation (2) below.
∆ρ (2)
=B⋅(S−k)
∆t
B is an empirical remodeling constant, S is the mechanical stimulus and k is the reference stimulus, also
frequently referred to as the attractor state. All calculations are performed on a per element basis, i.e. there
is one sensor point per element, so the mechanical stimulus can be calculated from the structural results of
the FE model. As outlined by Weinans, the strain energy density will be used as the measure for the
mechanical stimulus which can be expressed per element as
U
S = e , (3)
e vol ⋅ρ
e e
where U is the strain energy per element, vol is the element volume, and ρ is the element density, with the
e e e
subscript e as the element number index. In other words, the mechanical stimulus is defined as the strain
energy density per unit mass (referred to as SED in this text), and the magnitude and the effect of the
stimulus is determined by comparing the SED distribution in the bone material to the reference stimulus k.
Once the mechanical stimulus is extracted from the finite element analysis according to equation (3), the
density increment for each element can be calculated with a given time increment using equation (2), and
the element density can be updated to ρ = ρ+∆ρ.
new
As the density changes for each element, the stiffness of the continuum needs to be updated with each
iteration as well. The structural properties of bone are related to density according to the following
relationship:
E=c⋅ργ, (4)
where E is the modulus of elasticity, the constant c and the exponent γ are remodeling constants. At the
end of each iteration, the mechanical properties for each element can be calculated according to equation
(4).
The calculation steps described above represent one remodeling iteration. In every subsequent iteration,
the mechanical stimulus and its proximity to the reference value will be evaluated for each element, where
the overall structural response of the bone volume is based on the biomechanically adapted stiffness of the
individual elements from the previous iteration. This iterative process continues until all bone mass is in
remodeling equilibrium (as discussed in the Convergence Checking section).
Challenges and Simplifications
While the equations and the iterative solution procedure presented are quite trivial, there are several issues
which come to mind when given the assignment of implementing them into a finite element program, and
specifically into ANSYS: How much of the standard ANSYS functionality can be used for preprocessing,
solution, and post-processing of the bone remodeling phenomenon? What mesh density should be used?
Which constitutive model in ANSYS is suitable for modeling bone? What methods are available in
ANSYS to update the mechanical properties for each element throughout the solution? Is there a better
way than trial and error to find a suitable reference stimulus value for each remodeling simulation? What
time step size is appropriate? This will be an iterative approach to a nonlinear problem – at what point
should the solution be terminated or considered converged? How can the user monitor this nonlinear
solution and determine whether the model response is “well-behaved”? Once we have results, how can we
verify and validate them? How can we plot the key result from the remodeling simulation: density?
These questions are addressed in more detail below.
Bone Tissue Properties
There are two main types of bone tissue: cancellous (“spongy”) and cortical (“compact”) bone, as
illustrated in Figure 1 showing a section through a tibia. As the alternate designations imply, the main
difference between the two is how tightly the tissue is packed together. While cortical bone is similar to a
solid mass, cancellous bone is less dense and consists of trabeculae (plates and bars) forming an irregular
network of cavities. The trabeculae are of particular interest to the bone remodeling simulation because of
their ability to rearrange themselves in response to changing loading conditions.

Figure 1. Cross-section through tibia showing morphological features of bone
Bone tissue is visco-elastic and its properties are anisotropic. Its stiffness varies for different bones in the
body and is also dependent on the age and health of a person. While some of the researchers in the field of
bone remodeling focus their efforts on taking into account these complex aspects, it was demonstrated by
several authors that an isotropic material model is quite suitable for predicting realistic bone morphology.
The simulation of bone with BORA is based on the following assumptions: bone tissue is a linear isotropic
material and its density ranges between zero (fully resorbed bone, corresponding to a void) and a maximum
value of 1.74 gm/cm3, a typical value for human cortical bone. The exponent γ in equation (4) above is
assigned a value of 2.71. The remodeling constant c = 3790 MPa/(gm·cm3) is based on the paper by
Weinans. To avoid singularities in some of the calculations to be performed, the minimum density value
used by BORA is 0.01 gm/cm3. Given equation (4), the Young’s modulus of bone will range between 0.01
MPa and 17,000 MPa. A Poisson’s ratio of 0.3 is the default with BORA.
Model Validation
With the basic coding for BORA completed, the first step was to check for blatant APDL programming
errors in BORA. The bone remodeling process represents a non-typical FEA application and we cannot
rely on experience and intuition when attempting to validate results. The paper by Weinans is a useful
source in this respect because he studied remodeling results using a simple 2-D plate model with a ramped
edge load. We can have a much greater level of confidence in the implementation of the basic remodeling
algorithm if the density distribution results presented in the Weinans paper can be reproduced.
As a second step, a 2-D model of a femur was selected as a benchmark. Typical muscle and joint loads
were applied and all bone material was assigned the same initial density value. Based on the evolution
equations presented, the remodeling algorithm was expected to alter the initially uniform density
distribution in the bone model and produce the typical femur morphology (a hollow shaft with a cortical
layer and the trabecular architecture in the vicinity of the neck).
With enough confidence established, we could then venture into simulating a three-dimensional femur
model and eventually include an implant.
Density as a DOF in a Nonlinear System
General purpose FEA codes use nodal displacements as the primary degree of freedom (DOF) in a
structural analysis, strains and stresses are the results typically calculated, and material properties remain
constant throughout the solution. This is in conflict with the density evolutions to be simulated, where it
becomes necessary that density not be a fixed material property but act as an additional DOF for the
ANSYS solution. Density is also the main result quantity of interest, which dictates the structural response
of the finite element model because of its relationship with the modulus of elasticity. Furthermore, the
density DOF has an imposed upper and lower bound.
ANSYS uses the Newton-Raphson method to solve nonlinear systems, comparing the internal forces
resulting from the internal stresses to the external forces, and iteratively updating the displacements
accordingly until equilibrium is established. The nonlinear system of remodeling bone tissue requires that
we evaluate each element’s proximity to a reference stimulus and repeatedly update the element density
based on the equations presented until remodeling equilibrium is reached. This cannot immediately be
done in ANSYS without some coercion.
Rather than modifying the existing nonlinear ANSYS solution methods and developing custom element
types and material models, it appeared more convenient to devise a separate algorithm which repeatedly
goes through the basic sequential steps of a static linear solution, post-processing results, checking for
remodeling equilibrium, modifying structural properties, and solving again.

1 Realistic values for γ have been found to lie in the range of 2.0 and 3.0 (Rice et al., 1988; Hodgskinson
and Currey, 1989).
Loading
The load environment for bones is complex with force components resulting from joint contact, muscles,
tendons, and soft tissue, and the load magnitudes and directions vary constantly as a person moves about.
The close relationship between the bone remodeling response and the load environment suggests that great
care must be taken with the proposed femur benchmark to accurately apply all femoral loads similar to the
tissue force components identified qualitatively with the left image in Figure 2.
Similar remodeling simulations have shown that realistic bone morphology may be predicted by reducing
the complex load environment and the load history to two main components: the hip joint load and the
abductor load as shown on the right in Figure 2. The load values and directions indicated correspond to a
static load case representing the stance phase of gait for a person weighing 80 kg.
For the structural analysis in ANSYS, the muscle and joint load are applied as tractions (directional
pressures) using SURF154 elements (with KEYOPT(4)=1) at the locations indicated with the purple patches
in Figure 2.
Hip Joint Contact
Fx= 785 N
Fy= -589 N
Fz= 2158 N
XXX
YYY
ZZZ Abductor Muscle
Fx= -305 N
Fy= 0
Fz= -1138 N
Fixed End

Figure 2. Schematic of femoral loads2producing surface strains (left) and simplified load
environment for BORA (right)
Reference Stimulus
Equation (2) illustrated the significance of the reference stimulus k as a key parameter for the adaptive
remodeling response. It controls whether an element experiences a positive or a negative density increment
for a given SED state (over-stimulated elements will lead to an increase in bone mass, under-stimulated
elements are driven towards becoming a void). Clearly, the value k for the FE simulation must be chosen
carefully.

2 taken from Duda et al., (1998)
The “correct” value for the reference stimulus k is highly dependent on the structure and the load
environment, and realistic bone material architecture will be produced for only a very specific range of
values. The suitability of a selected value for the attractor state can only be determined by inspecting the
results from a remodeling simulation.
Choosing too large a k value will leave all elements under-stimulated, producing a volume with all bone
completely resorbed. In contrast to a value which is too small, in which case all elements may assume the
maximum density because they are over-stimulated.
A better method than the ever-frustrating “HI-LO trial and error” has to be found to give the user at least a
starting point from which the magnitude of the reference stimulus may be tweaked to obtain realistic bone
density distributions.
Solution Efficiency
All solution procedures described in the related papers studied for the development of this FEM approach
use a fixed time step for each remodeling iteration. The density increment for each element, ∆ρ, is given
e
with equation (5).
∆ρ =∆t⋅B⋅(S −k) (5)
e e
Because ∆t in this equation applies to all elements and for every iteration, the fixed time increment
approach forces a choice between accuracy and efficiency. Too large a time increment will cause the
elements to “overshoot” their targets: the density extrema and the attractor state. Using a small time
increment would be an obvious solution, but at the cost of solution efficiency.
BORA is intended for use as a simulation tool allowing engineers to solve complex bone shapes with at
least 100,000 elements within a few hours. With this goal in mind, an alternate approach which combines
accuracy and efficiency was developed for BORA by adapting the time step increment to the convergence
behavior of the model and its remodeling history as explained in section Adaptive Time Stepping.
Convergence
It is desirable to allow the user to control the accuracy of the solution, in other words, a convergence
criterion and convergence tolerance for the remodeling equilibrium has to be established. We also want to
provide the user with some feedback about the convergence status of the model and to give an option to
“cleanly” terminate the solution if necessary.
The following sections detail how these features were implemented into BORA.
Implementation – Basic Concepts
BORA vs. ANSYS
Any ANSYS database using 2-D or 3-D structural elements may readily be used for a bone remodeling
simulation, as long as the FE model has some load applied and is fully constraint. No re-compiled version
of ANSYS is required to run BORA. BORA is merely a set of input decks with APDL code, and the
complete BORA functionality (solution, post-processing, other utilities) is accessed from the toolbar.
Beyond these similarities, there are several unique characteristics of BORA where it differs greatly from
the general use of ANSYS.
• There is no need to define a Young’s modulus, but all elements must have density defined.
• A remodeling convergence history graph is displayed during solution.
• BORA results quantities are unique (no stresses or displacement - BORA uses density as a result).
• Results are stored in arrays and will not readily be accessible using the general /POST1 commands.
Remodeling Equations in ANSYS
The basic steps performed by BORA are illustrated in ANSYS terms with the flowchart in Appendix A.
The substitutions required to solve the bone remodeling equations are easily obtained from a structural
analysis in /POST1: strain energy may be extracted via ETABLEs, element volume and density are
available via *GETs, which is all the data required to calculate S (or SED), the strain energy density per
e
unit mass for each element.
Calculations and data manipulations are performed using a combination of *VOPER, *VPUT, and element
table operations and sorting. Further discussions of the programming for BORA are limited to the key
ANSYS commands highlighted in the flowchart.
Convergence Checking
Expressed in simple terms, a remodeling solution is converged when all elements have reached one of three
states:
• The element density has reached the minimum value, i.e. all bone tissue has fully resorbed and a void
has formed at the element location.
• The element has reached the maximum density value, i.e. the element has turned into cortical bone.
• The elemental strain energy density is equal to the reference stimulus.
To allow the solution to continue until either of these three conditions is met for ALL elements is desirable
but computationally an unrealistic expectation. As for all nonlinear and iterative solutions, we have to
allow a certain amount of error or deviation from these three idealized element states. This admissible
amount of error is defined through the use of two key convergence parameters in BORA.
The ANSYS parameter ERR_K controls how close to the attractor state K an element has to be before it can
be considered converged. It defines a proximity band with a width of ERR_K*K centered around the
reference stimulus. Elements with SED values in the range of ERR_K/2-K and ERR_K/2+K will be treated
as if they had reached the attractor state. The parameter CNV_CNV allows the user to control the
percentage of elements expected to meet either of the three converged states.
But, as with any nonlinear iterative solution, a converged solution does not necessarily imply an accurate
solution. We still need to establish that the predicted bone morphology reflects an adequate and realistic
density distribution.
Element Status
An element is considered “active” as long it has not reached the minimum or maximum density level. Only
active elements are considered for the convergence metrics and the decision logic used for the adaptive
time stepping scheme introduced below. Active elements which do not fall within the convergence
tolerance band near the attractor state are assigned a status of “remodeling”.
Implementation – Enhancements and Automation
Considering BORA’s intended use as an engineering tool, several enhancements and convenience features
were added to the program.
Adaptive Time Stepping
Because the time increment is an arbitrary factor used to calculate the density increment and the
remodeling process described disregards any transient effects, it is just as feasible to use a fixed density
increment instead of a preset time increment. In other words, we will allow BORA to adjust the time
increment to suit the goal of solution efficiency such that a desirable density increment is maintained
throughout the solution according to the following equation:
∆ρ (6)
∆t = fixed
m 1 m  U 
B⋅ ∑ e −k
m  vol ⋅ρ 
e=1 e e 
The strain energy density per unit mass S from equation (2) was substituted with the average value for a
group of m elements and solved for ∆t. As this averaged strain energy density changes throughout the
solution, we now have a variable time step size ∆t which is a function of a conveniently selected density
m
increment ∆ρ and a representative group of m elements.
fixed
But which m elements should dictate the calculations in equation (6)? Mainly driven by observations of the
convergence history of a large number of test model, the program employs two different strategies for
calculating a “weighted” time increment based on two groups of elements: over-stimulated and under-
stimulated elements. For each group, the average SED is calculated and the difference to the reference
stimulus determined which yields the corresponding average time increment for each group. With ∆t for
hi
the over-stimulated group and ∆t for the other, BORA calculates a weighted time increment ∆t' as shown.
lo
∆t ⋅m +∆t ⋅m
∆t′= hi hi lo lo (7)
m +m
hi lo
The time stepping algorithm is called “adaptive” because the before-mentioned weighted time increment
method is used until the SED of a given percentage of elements (a value of 80% proved most successful)
fall within a proximity band near the reference stimulus. From that point on, a more conservative time step
calculation method based on the maximum density increment for any element takes over, as shown with
equation (8).
∆ρ (8)
∆t′′= fixed
 U 
B⋅max e −k
 
vol ⋅ρ 
e e
Program-Suggested Reference Stimulus
The goal of the remodeling process is (S–k)=0 for all active elements. To provide the remodeling
simulation with a reasonable starting value for k, BORA calculates a suggested value k during the first
auto
iteration per equation (9) below, which will be used by the program unless the user inputs a specific value.
k is based on the ad hoc assertion that the average SED per element at the state of remodeling
auto
equilibrium is reduced to a percentage of the initial SED value (c = 0.4 works well).
k
1 n
k =c ∑S (9)
auto k n e
e=1
Minimum DOF requirement
In addition to the convergence tolerance controls discussed earlier, the remodeling progress is monitored to
catch a solution which has only minor density changes in the active elements. Independent of the element
convergence status, the solution will terminate when the average change in density for all active elements is
less than a set percentage of the maximum density ρ (suggested values are in the neighborhood of
cb
0.02·ρ )
cb
Solution Monitoring
There are four ways in which BORA provides the user with feedback about the progress of a remodeling
simulation: a status bar, a graph displaying the remodeling convergence statistics, a parameter file which is
updated for each iteration, and a convergence history file.
The status bar displayed during a BORA run (Figure 3) shows the percent complete of the remodeling
simulation based on the maximum number of iterations specified. The status bar also provides an option to
cleanly abort the solution and to restart it at a later point.

Figure 3. BORA solution status bar with option to terminate

Figure 4. Convergence history plot for BORA simulation
Throughout the BORA run, a set of convergence metrics is graphed and updated with each iteration
showing vital information about the convergence behavior of the remodeling solution in progress. A
typical convergence history for a "well-behaved" simulation is shown with Figure 4 below. The
significance of each graph is discussed below. All metrics are expressed as percentages and are graphed as
a function of the remodeling iteration number.
e_minmax tracks how many elements have reached either the minimum density (fully resorbed) or
maximum density (cortical bone) relative to all bone elements in the model. e_cnvrg displays a measure of
the ratio between converged elements (at the specified reference stimulus proximity) and all active
elements (elements which are NOT at either density extreme). e_remodl shows the percentage of elements
that are remodeling for each iteration.
sat_levl is not really a convergence measure but a means of monitoring the average normalized element
density level, the bone density saturation level.
Whenever a significant change in the convergence history occurs, an entry is added to a convergence
history file (shown in Figure 5). This feature is intended to help the user better understand the model
response and to assist in debugging undesirable remodeling solutions.
Figure 5. Convergence history file created during BORA solution
Option to Restart
With all general model and element-specific information (SED, ρ, etc.) stored in arrays, the solution can be
stopped and restarted from a continuously updated restart file which contains all parameters and arrays.
The model state immediately before the abort can be recreated by resuming the original database and then
restoring all parameters and filling all arrays. This restart procedure is automated and available from the
toolbar.
BORA Utilities
Frequently performed tasks are automated in the form of macros which are accessed through the BORA
toolbar. A few automated tasks of note are the saturation level plot, the remodeling status plot, and the
BORA report in HTML format.
Other post-processing options available from the toolbar include the graphing of the history for various
remodeling metrics such as the change in total bone mass, the change in the overall stiffness, the maximum
and average density increment per element, and the average proximity to the reference stimulus. The
convergence history displayed during solution may be recalled and displayed as well.
Visualizing Results
The results of interest from a remodeling solution are the density distribution and the strain energy density
per unit mass relative to the reference stimulus, none of which are available directly from an ANSYS
solution. BORA uses *VPUT to move all relevant data from arrays to ETABLEs, which can then
conveniently be sorted, grouped, and displayed as contour plots. Again, the post-processing toolbar
provides access to macros which perform all data manipulation and plot commands for viewing the
remodeling results.
Two of the automated results plots in BORA are illustrated by means of a familiar structure, a short
cantilevered beam with a localized pressure load applied (refer to Figure 6). However, this beam is unique
because it is made of bone and free to adapt itself to the loading environment.

Figure 6. Half-symmetry model of short cantilevered bone beam with pressure load