Finite Element Analysis of Bone Remodeling - Implementation of a Remodeling Algorithm in MATLAB and ANSYS

Abstract

English: The process of adaptive bone remodeling can be described mathemati- cally and simulated in a computer model, integrated with the finite element method. The main focus of this thesis is the implementation of a bone remod- eling algorithm in MATLAB and ANSYS on the basis of FEM. The strain energy density is used as mechanical stimulus. The cortical and trabecular bone are described as continuous materials with variable density. This thesis can be divided into four main parts. The first part is due to the ma- terial properties of cortical and trabecular bone. The second part is about the remodeling theory and gives an historical review of the developed numerical approaches up to now. The implementation of the remodeling algorithm in ANSYS and MATLAB as well as its validation is topic of part three. In last main part, the algorithm is applied to a 2D FE-model of a human proximal femur.

Deutsch: Der Vorgang des adaptiven Knochenumbaus kann mathematisch beschrieben und mit Hilfe der Finiten Elemente Methode an einem Comput- ermodell simuliert werden. Das Hauptaugenmerk dieser Arbeit liegt in der Implementierung eines Algorithmus in MATLAB und ANSYS, mit dem der Knochenumbau auf Basis der FEM simuliert werden kann. Die Dehnungsen- ergiedichte dient hierbei als mechanischer Stimulus. Die Spongiosa und Ko- rtikalis werden dabei als kontinuierliches Material mit variierender Dichte beschrieben.

Die Arbeit kann in vier Hauptteile gegliedert werden. Im ersten Teil wer- den die Materialeigenschaften von kortikalen und spongiösen Knochen be- trachtet. Der zweite Teil geht über die Theorie des Knochenumbaus und gibt einen historischen Rückblick über bisher entwickelte numerische An- sätze. Die Implementierung des Algorithmus in ANSYS und MATLAB, sowie dessen Validierung ist Thema des dritten Teils. Im letzten Teil wird der Algorithmus an einem 2D FE-Modell eines menschlichen proximalen Femurs angewendet.

Keywords Bone remodeling, finite element method, strain energy density, mechan- ical stimulus

ii

Preface

This Master thesis project was carried out between January 2005 and February 2006 at the Lehrstuhl für Statik at the Technische Universität at München in coopera- tion with the Labor für Werkstoffkunde und Metallographie (LWM) at the Fach- hochschule Regensburg.

First of all I like to thank the first supervisor of this thesis, Sebastian Dendorfer, for his excellent support in every respect. His deep knowledge regarding computational mechanics and biomechanics in general and especially regarding the remodeling theory of bone proofed vital for this thesis.

I would also like to thank all my colleagues of the LWM, especially Prof. Dr. Joachim Hammer, for providing the necessary software and hardware, and for giv- ing me the opportunity to work at the LWM.

And finally I would like to thank Dr. habil. Manfred Bischoff, who also supervised this thesis from the Lehrstuhl für Statik.

Martin Groß 1 , München March 2006

this thesis or the implementation

iii

iv

Contents

Contents v

List of Figures viii

List of Tables xi

1 Introduction 2

1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Aim . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 Bone as Material 4

2.1 Composition of Bone . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.2 Difference between Cortical and Cancellous Bone . . . . . . . . . . 5

2.3 Material Properties of Bone Tissue . . . . . . . . . . . . . . . . . . 8

2.3.1 Cortical Bone Tissue . . . . . . . . . . . . . . . . . . . . . 8

2.3.2 Cancellous Bone Tissue . . . . . . . . . . . . . . . . . . . 10

3 Bone Remodeling - From Nature to Model 14

3.1 Remodeling Theory . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3.1.1 Difference between Modeling and Remodeling . . . . . . . 14

3.1.2 The Remodeling Process . . . . . . . . . . . . . . . . . . . 15

. . . . . . . . . . . . . . . . . . . 16

Remodeling Cycle Duration . . . . . . . . . . . . . . . . . 20

CONTENTS

3.1.3 Types of Remodeling . . . . . . . . . . . . . . . . . . . . . 21

Osteonal Remodeling . . . . . . . . . . . . . . . . . . . . . 21 Trabecular Remodeling . . . . . . . . . . . . . . . . . . . . 22 Endosteal and Periosteal Remodeling . . . . . . . . . . . . 22 3.2 Theories from the Beginning up to now . . . . . . . . . . . . . . . 24 3.2.1 Mechanically Excited Bone Adaption Theories (1865 -1920) 24 3.2.2 Bone Adaptation: General Relationships of Mechanics to

Bone Physiology (1920 - 1970) . . . . . . . . . . . . . . . 27 3.2.3 Bone Adaptation: Experimental Study of Mechanically

Mediated Bone (1970 - 1984) . . . . . . . . . . . . . . . . 28 3.2.4 Theories of Bone Adaptation: Numerical Simulations

(1985 to present) . . . . . . . . . . . . . . . . . . . . . . . 31

4 Simulation of Remodeling 35

4.1 Implementation of Optimization Algorithm . . . . . . . . . . . . . 35 4.1.1 Reference System . . . . . . . . . . . . . . . . . . . . . . . 36

Material Properties . . . . . . . . . . . . . . . . . . . . . . 38 Load . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.1.2 Optimization Part . . . . . . . . . . . . . . . . . . . . . . . 39

Stimulus . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 Adaptation Functions . . . . . . . . . . . . . . . . . . . . . 41 Convergence Criteria . . . . . . . . . . . . . . . . . . . . . 43 4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.2.1 Reference System 1 . . . . . . . . . . . . . . . . . . . . . . 44 4.2.2 Reference System 2 . . . . . . . . . . . . . . . . . . . . . . 46 4.2.3 Convergence behavior . . . . . . . . . . . . . . . . . . . . 47

vi

CONTENTS

5 Applications 48 5.1 Modeling of Human Proximal Femur . . . . . . . . . . . . . . . . . 49 5.1.1 Load Definitions . . . . . . . . . . . . . . . . . . . . . . . 50 5.1.2 Initial Configuration and Boundary Conditions . . . . . . . 52 5.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.2.1 Proximal Femur with Stepwise Adaptation Function . . . . 55

Initial Homogeneous Density Distribution . . . . . . . . . . 55 Initial Stochastic Density Distribution . . . . . . . . . . . . 58 5.2.2 Proximal Femur with Linear Adaptation Function . . . . . . 61

Initial Homogeneous Density Distribution . . . . . . . . . . 61 Initial Stochastic Density Distribution . . . . . . . . . . . . 62 5.2.3 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . 63

6 Discussion 65

6.1 Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 6.2 Details of Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 65 6.2.1 Building the Model . . . . . . . . . . . . . . . . . . . . . . 65 6.2.2 Development of Bone Structure . . . . . . . . . . . . . . . 66 6.2.3 Convergence Behavior . . . . . . . . . . . . . . . . . . . . 68 6.3 Conclusion and Outlook . . . . . . . . . . . . . . . . . . . . . . . 68

Index 70

Bibliography 72

A Implementation in MATLAB and ANSYS A

A.1 MATLAB main file run.m . . . . . . . . . . . . . . . . . . . . . . A A.2 ANSYS macro . . . . . . . . . . . . . . . . . . . . . . . . . . . . E A.3 MATLAB function dr.m . . . . . . . . . . . . . . . . . . . . . . . . K

B Basic Anatomic Terminology L

vii

List of Figures

2.1 Sketch of some important features of typical long bone, from [93] . 6

2.2 SEM micrograph of ground trabecular vertebra bone, from [76] . . . 6

cube) from [4] . 6

2.4 The influence of loading rate on the tensile strength and modulus of

cortical bone from [92]. . . . . . . . . . . . . . . . . . . . . . . . . 8

2.5 Typical stress strain curves for trabecular bone of different densities,

from [43]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.6 Stress against strain of trabecular bone specimen under compression

from [76]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

against density

From [85, 14, 2, 57, 64, 79]. . . . . . . . . . . . . . . . . . . . . 13

3.1 Photomicrograph of an osteonal basic multicellular unit . . . . . . . 16

3.2 Schematic sketch of an osteonal BMU. Cross-sectional view at the

bottom right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.3 The Six Phases of an Osteon’s Lifetime. a) Activation AC, b) Re-

sorption RES, c) Reversal REV, d) Formation FO, e) Mineralization MI, f) Quiescence QU. from Ott with permission [75] . . . . . . . . 19

3.4 BMU activation rate vs. age for human ribs. (From data by [34]) . . 22

3.5 Stress trajectories in curved Culmann crane (left) compared with a

schematic representation of the trabecular pattern in the proximal femur, from Wolff, 1870 . . . . . . . . . . . . . . . . . . . . . . . 24

3.6 Change of trabecular structure in post-fracture , from Wolff, 1870 . 25

3.7 Frost’s description of the different adaptive responses for the ado-

lescent and the adult skeleton. . . . . . . . . . . . . . . . . . . . . 29

LIST OF FIGURES

3.8 Density distribution in the femoral head by Fyhrie. . . . . . . . . . 33 3.9 Adaptation function according to equation (3.6) . . . . . . . . . . . 34

4.1 Simulation algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . 35 4.2 Reference systems. . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.3 Higher order 2-D element, from Zienkiewicz [94]. . . . . . . . . . 37 4.4 Young’s modulus against discretized density ρ according to equa-

tion (4.2) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.5 Defined load steps i against time steps. . . . . . . . . . . . . . . . . 40 4.6 Overview of implemented adaptation functions. . . . . . . . . . . . 41 4.7 Resulting material (density) distribution of reference system 1. Fig-

ure 4.7(o) is done with Ole Sigmud’s code from [86] . . . . . . . . 44 4.8 Resulting material distribution of reference system 2. Figure 4.8(o)

is done with Ole Sigmud’s code from [86] . . . . . . . . . . . . . . 46 4.9 Convergence plot of reference system 1. . . . . . . . . . . . . . . . 47 4.10 Convergence plot of reference system 2. . . . . . . . . . . . . . . . 47

5.1 Anatomy of the human proximal femur. From [89]. . . . . . . . . . 48 5.2 3-D femur model with section plane. . . . . . . . . . . . . . . . . . 49 5.3 2-D Finite element mesh of the proximal femur with 7124 elements. 49 5.4 Element FLUID79 from [1] . . . . . . . . . . . . . . . . . . . . . . 50 5.5 Proximal femur with muscles. . . . . . . . . . . . . . . . . . . . . 51 5.6 Hip contact force against time for human normal walking from [7]. . 51 5.7 Overview load cases for normal walking. . . . . . . . . . . . . . . . 52 5.8 Initial configurations with two different density distributions. . . . . 53 5.9 Remodeling ratio coefficient B(n)in ( g

cm 3 ) 2 MPa −1 . . . . . . . . . . 53

5.10 Remodeling progress in human proximal femur with a stepwise

adaptation function and initial homogeneous density distribution. . . 55 5.11 Comparison of v. Mises stresses at initial and converged state. . . . 56 5.12 Comparison of principal stresses at initial and converged state. . . . 57

ix

LIST OF FIGURES

5.13 Remodeling progress in human proximal femur with a stepwise

adaptation function and initial stochastic density distribution. . . . . 58 5.14 Comparison of v. Mises stresses at initial and converged state. . . . 59 5.15 Comparison of principal stresses at initial and converged state. . . . 60 5.16 Remodeling progress in human proximal femur with a linear adap-

tation function and initial homogeneous density distribution. . . . . 61 5.17 Remodeling progress in human proximal femur with a linear adap-

tation function and initial stochastic density distribution. . . . . . . 62 5.18 Sum of averaged strain energy density u ∗ . . . . . . . . . . . . . . . 63 5.19 Difference ˙

|q| in linear-linear scale. . . . . . . . . . . . . . . . . . . 64 5.20 Difference ˙

|q| in linear-logarithmic scale. . . . . . . . . . . . . . . 64

6.1 Comparison of converged results with X-ray plot of proximal femur. 67

B.1 Anatomic planes with labels from [61] . . . . . . . . . . . . . . . . M

x

List of Tables

2.1 Strength of femoral cortical bone from [43]. Mean values from [78]. 9 2.2 Moduli of femoral cortical bone, from [43]. Mean values from [78]. 10 2.3 Mean values of modulus and ultimate strength for various anatomic

sites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

5.1 Joint and muscle forces during normal walking according to

Bergmann [7]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

B.1 Body components . . . . . . . . . . . . . . . . . . . . . . . . . . . L B.2 Terms of position . . . . . . . . . . . . . . . . . . . . . . . . . . . M

Chapter 1: Introduction

Chapter 1

Introduction

1.1 Motivation

Bone is a living material which has its main function in building the skeleton and therefore enabling locomotion and protection of the organism. It is subjected to permanent and transient loads caused by the daily active or special events like acci- dents. In contrast to inert materials from standard mechanics, this tissue is able to react adaptively to its environment. Aside from skeletal growth and fracture heal- ing, which are of temporary character, the internal bone structure is maintained and adapted continuously. This process is termed remodeling .

Within this process microdamage is removed, leading to an increase of the fatigue life of bone tissue. Furthermore, the structural adaptation to changes in the mechan- ical environment plays an important role in conjunction with implants and prosthe- ses. In fact, latest developments of such devices have been analyzed numerically in order to predict the long-term reaction of the tissue to this impact. Osteoporosis, nowadays a widespread bone disease, underlies similar concepts as the remodeling process. Due to the enormous social damage caused by such dis- eases on one side and failure of implants and prostheses on the other side, an ad- vance in the understanding and computer simulation of remodeling is of great im- portance.

2

Chapter 1: Introduction

Therefore, the interest in this phenomenon has been increasing in the last century. Especially in the last twenty years, many numerical algorithms have been developed in order to simulate this process. Although early models were capable of predicting good approaches to the real behavior, a quantitative analysis has been impossible and many biologic aspects have been neglected.

Recently, new methods have been published which take into account aspects of the microstructure and cell activities. These models provide good results for predictions of bone loss due to osteoporosis but suffer deficiencies from a mechanical point of view.

In the following, the composition of bone, the difference of cortical and trabecular bone and its mechanical properties are described in chapter 2. Afterwards the pro- cess of remodeling will be outlined in chapter 3 coming with a historical review of research work done in this field from the of the 19th century up to now. In chapter 4 the implementation of the remodeling algorithm is described and simulations for

evaluation are shown. The implemented algorithm is applied to a proximal femur FE-model with different initial configurations in chapter 5. Finally, the results are discussed in chapter 6 and the appendix provides the source code for implementa- tion in MATLAB and ANSYS.

1.2 Aim

The aim of this work is to implement an easy to use and extendable numerical algorithm which can build up the remodeling process of bone due to mechanical stimulus. Besides mechanical stimulus, a huge number of items, such as age, race, gender, possible diseases etc. play an import role in this context. These highly com- plex and mostly underdetermined boundary conditions of the optimization, how- ever, shall not be taken into account in this approach. Therefore, this thesis will be mainly focused on bone’s internal structural changes in response to a change in the mechanical environment.

3

Chapter 2: Bone as Material

Chapter 2

Bone as Material

2.1 Composition of Bone

Bone is not an uniform material, it is composed of collagen, water, hydroxyapatite mineral and small amounts of proteglycans and noncollagenous proteins . [67]

Collagen is a structural protein, that can spontaneously organize itself into strong fibers. More than a dozen types of collagen have been identified. Type I collagen is the most abundant collagen of the human body. It is present in scar tissue, the end product when tissue heals by repair. Beside in bone, it is also found in tendons, lig- aments and skin. Collagen is responsible for bone’s flexibility and tensile strength. It also provides loci for the nucleation of bone mineral crystals, which give bone rigidity and compressive strength.

Mineral in bone consists almost entirely of hydroxyapatite crystals, Ca 10 (PO 4 ) 6 (OH) 2 . The individual crystals are rods with hexagonal symme- try, measuring about 50 × 50 × 400 angstroms (1 [µm] = 10,000 angstroms [Å]). Bone mineral is impure, containing many structural substitutions (e.g., carbonate, fluoride, citrate). These impurities are governed by the composition of body fluids and in turn affect the solubility of the bone mineral.

4

Chapter 2: Bone as Material

Ground substance of bone consists of proteglycans. In particular, decorin and biglycan are small species of proteglycans found in bone. Although the specific role of the proteglycans is not known exactly. Decorin is known to modulate col- lagen fibril assembly. Proteglycans may also act to control the location or rate of mineralization in through their calcium-binding properties.

Noncollagenous proteins contain quite a few molecules whose functions are also unclear. The most abundant noncollagenous protein is osteocalcin, which is se- creted by osteoclasts and appears to be important in the mineralization of new bone. It also is a chemoattractant for bone cells and or its serum concentration are an excellent method of noninvasively determining rates of bone turnover. Other non- collagenous proteins in bone include osteopontin and osteonectin.

Water appears in the calcified bone matrix in two different conditions, one part is free and the other part is bound to other molecules. The mineralization of osteoid (the organic portion of extracellular bone) displaces part of its water. Therefore, the water content of new bone tissue changes as it mineralizes.

2.2 Difference between Cortical and Cancellous

Bone

In principle there are two types of bone, as determined by porosity (volume fraction of soft tissues). The porosity of bone can vary step less from zero to 100%, but most bone tissues are of either very low or very high porosity and just a small part of intermediate porosity. These two types of bone tissue are referred to compact or cortical bone and trabecular or cancellous bone, respectively as shown in Figure 2.1. This combination of trabecular and cortical bone forms a sandwich-type

structure, well known in engineering for its optimal structural properties. [27]

5

Chapter 2: Bone as Material

Figure 2.1: Sketch of some important features of typical long bone, from [93]

Trabecular bone is the spongy, porous type of bone as shown in Figure 2.2 and can be found in the cuboidal bones, the flat and irregular bones, such as the sternum, pelvis and spine (vertebra) and at the end of all long bones [54]. Its porosity varies from 75% to 95%. The interconnected pores, which scale is on the order of 1 mm, are filled with marrow. The bone matrix is in the shape of plates or struts called trabeculae each with about 200 µm in diameter. The arrangement of the trabeculae is not unique. Mostly they build a randomly oriented meshwork. Sometimes they appear to be organized into orthogonal arrays.

6

Chapter 2: Bone as Material

Compact bone is the dense bone crop up in shafts of long bones and forming a cortex or shell around vertebral bodies and other cancellous bones. Its porosity varies from 5% to 10% and its pores consist of spaces categorized as follows

Haversian canals are approximately aligned to the long axis of the bone. They are about 50µm in diameter and contain nerves and capillaries. Haversian canals are named after an English physician, Clopton Havers (1691).

Volkmann’s canals are short, transverse canals connecting Haversian canals to each other and to the outside surfaces of the bone. These canals also include blood vessels and presumably nerves. They are named after Richard von Volkmann (1830-1889), a surgeon and early advocate of Lister’s antiseptic surgical methods.

Resorption cavities are temporary spaces caved by osteoclasts in the initial stage of remodeling, described in 3.1.2. These cavities are about 200µm in diameter.

It is important to keep in mind that bone is a dynamic porous structure. Its porosity may change as a result of a pathologic condition or in a normal adaptive response to a mechanical or physiologic stimulus. This leads to higher density in trabecular bone, or to lower density in compact bone. Such changes strongly affect bone’s mechanical properties.

Besides trabecular and compact bone, two further aspects how bone can be characterized should be mentioned. Namely lamellar vs. woven bone and primary and secondary bone. More details can be found in [67, 72].

7

Chapter 2: Bone as Material

2.3 Material Properties of Bone Tissue

To determine the mechanical properties of bone tissue, small uniform specimens are loaded under well-defined conditions. In homogeneous materials, such testing conditions lead to uniform stresses throughout the specimen. The resulting defor- mation can be measured and the stress-strain relationship can be established. With several load types, material properties in tension, compression, bending and torsion can be determined [8, 10, 9, 17, 26, 25, 78, 77, 85].

2.3.1 Cortical Bone Tissue

The material properties of cortical bone are influenced by several factors. One is the rate at which the bone tissue is loaded. Rapid loading of cortical bone specimens leads to increased elastic moduli and ultimate strength as compared to specimens loaded more slowly. To quantify the rate of deformation, one can refer to the strain rate [s −1 ] to which the tissue is exposed. In normal activities, bone is subjected to strain rates generally below 0.01 s −1 . Materials like bone, are said to be time-

dependent or viscoelastic , as shown in Figure 2.4. This means that the stress-strain

Figure 2.4: The influence of loading rate on the tensile strength and modulus of corti- cal bone from [92].

8

Chapter 2: Bone as Material

characteristics and strength properties depend on the applied strain rate. According to [32], however, this rate dependency is relatively weak. The stress-strain behavior of cortical bone is also strongly dependent on the ori- entation of the bone microstructure with respect to the loading direction. Cortical bone is stronger and stiffer in longitudinal direction (direction of osteon orienta- tion) than in transverse direction. In addition, bone specimens loaded in a direction perpendicular to the osteons tend to fail in a more brittle manner, with little non- elastic deformation after yielding. Materials such as bone, for which elastic and strength properties are dependent on the direction of applied loading, are said to be anisotropic materials . Therefore, both strain rate and direction of applied load has to be specified when describing material behavior. In Table 2.1 the ultimate strength values of adult femoral cortical bone under vari- ous modes of loading, in both the longitudinal and transverse direction, are summa- rized: These indicate that the material strength of bone tissue depends on the type

Table 2.1: Strength of femoral cortical bone from [43]. Mean values from [78].

of loading as well as on the loading direction. The compressive strength is greater than the tensile strength in both longitudinal and transverse directions. Transverse specimens are weaker the longitudinal specimens in both tension and compression. The shear strength (determined by torsion tests about the longitudinal axis and re- flection shear stresses along transverse and longitudinal planes) is about one-third of the compressive strength. The modulus values for adult femoral cortical bone are shown in Table 2.2.

9

Chapter 2: Bone as Material

The longitudinal elastic modulus is about 50% greater than the transverse elas-

Table 2.2: Moduli of femoral cortical bone, from [43]. Mean values from [78].

tic modulus. The shear modulus for torsion about the longitudinal axis is about one-fifth of the longitudinal modulus.

2.3.2 Cancellous Bone Tissue

Cancellous bone is a complex material with significant heterogeneity. Its elastic and strength properties vary across anatomic sites, with aging and disease. Tra- becular bone is classified from an engineering materials perspective as a compos- ite, anisotropic, open porous cellular solid. Like many biological materials, it dis- plays viscoelastic behavior, as well as damage susceptibility during cyclic loading [54, 29].

A critical issue that distinguishes trabecular bone from many other biological tis- sues is its substantial heterogeneity, which leads to wide variations in mechanical properties. This heterogeneity results from underlying variations in volume frac- tion, architecture and tissue properties, in that order of importance. Across sites and species, mean values of modulus and strength can differ by more than an order of magnitude (Table 2.3). Substantial loss of mechanical properties also occurs with aging in humans. For example, ultimate stress is reduced by almost 7% and 11% per decade for the human proximal femur and spine, respectively, from ages 20 - 100 [68, 71, 70]. Strength does not decrease significantly until after about age 30 [71, 70]. Trabecular bone is anisotropic in both modulus and strength [71]. Com- pared with such materials as fiber-reinforced composites, the extent of anisotropy is mild, but its biomechanical significance in terms of whole bone strength or bone- implant performance remains to be quantified. With increasing porosity for exam-

10