ABSTRACT

The trabecular bone fracture healing differs from diaphyseal fracture healing, in which trabecular bone heals based on intramembraneous ossification. The process includes a small callus formation, then woven bone forms, it follows by remodeling process to form regular trabecular bone. The objective of this study was to present an energy based model to simulate bone formation and remodeling during trabecular bone fracture healing. This modeling mainly focused on the mechanical factors. The model distinguishes three basic type of tissue: bone, cartilage and soft tissue. In order to determine tissue differentiation a fuzzy controller was proposed. An algorithm was developed to link the fuzzy logic controller to a finite element model (FEM) of trabecular bone. In general, finite element analysis provides input for fuzzy controller. Based on the input data, the fuzzy system selects the type of tissue to build. Strain energy density was used as the mechanical stimulus and a new parameter was incorporated in to the healing process as the remodeling index.

1. Introduction

Trabecular bone or cancellous bone is found in cuboidal bones, flat bones, and at the ends of long bones. Its porosity is between 50% and 95% and the pores are filled with bone marrow. Bone marrow is a tissue composed of blood vessels, nerves and various types of cells. The trabecular bone matrix consists of plates and struts called trabeculae [1] and trabecular bone fracture is due to breaking these trabeculae.

The trabecular bone fracture healing differs from diaphyseal fracture healing, in that trabecular bone heals in the form of intramembraneous ossification. In intramembraneous ossification bone is formed directly on the existing bone. The first stage of healing begins by resorption of necrotic tissue and formation of soft tissue in the fracture gap. In the next stage, osteoblasts start their activity around the existing trabeculae. After woven bone formation, bone cells start the remodeling process. During remodeling the woven bone, which can be considered isotropic and homogeneous, is replaced by anisotropic and inhomogeneous trabecular bone [2].

Several biomechanical models have been proposed to describe the influence of mechanical stimuli on the fracture healing process, particularly on diaphyseal fracture healings. Claes et al. [3] (1999) developed a theory which explains the relation between the local stress and strain. Bailon-Plaza et al. [4] (2001) presented a two dimensional mathematical model of bone healing, they assumed that concentration of growth factors regulate tissue differentiation. Lacroix et al. [5] (2002) used a biphasic poroelastic finite element model; they assumed that shear strain and fluid flow are the tissue differentiation stimuli. Bailon-Plaza et al. [6] (2003) incorporated the effects of mechanical stimulation on cell differentiation into their previous mathematical model of growth factors, which regulates fracture healing. Gomez-Benito et al. [7] (2005) proposed a mathematical model to simulate the effect of mechanical stimuli on the cellular processes, during healing, which were implemented in a finite element code as combination of three coupled analysis stages: a biphasic, a diffusion and a thermoelastic steps. Isaksson et al. [8] (2006) used deviatoric strain, pore pressure and fluid velocity as mechanical stimuli in healing process and they suggested that the deviatoric component of strain may be the most significant mechanical parameter to guide tissue differentiation during indirect fracture healing. Also fuzzy logic has been used in combination with finite element models with strain energy density [9]

and dilatational and distortional strain components [10] as mechanical stimuli to simulate tissue differentiation for the compact bone.For the first time Shefelbine et al. [2] (2005) developed a model which can simulate fracture healing and remodeling process in trabecular bone by combining the fuzzy logic and the FEM. They used hydrostatic strain and octahedral shear strain as mechanical stimuli for healing in trabecular bone.

The objective of our investigation was to present an energy based model to simulate bone formation and remodeling during trabecular bone fracture healing. This modeling mainly focused on the mechanical factors and newly introduced remodeling index.

2. Methods

2.1. FEA Model

In our modeling approach, finite element analysis (FEA) was widely used. FEA is a general technique for obtaining numerical solution to the mathematical equations governing problems of interest. Using this method, one can transform the partial differential equations into a set of linear algebraic equations which can be solved by digital computers. The essence of the FEA is to break large complex structures into smaller interconnected components of “finite size” called “elements”. Each element has a function which is assumed to satisfy the required differential equations over the volume of the element. Since the differential equations are only solved over the volume of the element, this leads to a “piecewise approximation” of the actual response of the domain. In FEA modeling of structural mechanics, governing differential equations are differential equations of equilibrium, and field variable (unknown) is displacement. Solution of FEA results in finding displacement field as the primary unknown, which finally leads to determination of secondary unknowns, such as strains and stresses.

Our FEA model consists of 100 two-dimensional quadratic elements and the fracture gap size was considered 360 µm with thetrabeculae spicules on both sides (Figure 1). In trabecular bone matrix, trabeculae have thickness of about 200 µm with a variable arrangement [1]. A fracture gap in trabecular bone was modeled using a micro-scale simple square (600*600 µm). This 2D model was similar to previous 3D model which was proposed by Shefelbine et al. [2] (2005). The fracture gap was filled with soft tissue to represent bone marrow and hematoma. To understand the effect of gap size on trabecular bone fracture healing, models with gap size of 120 µm and 600 µm were studied too.

The model distinguishes three basic type of tissue: bone, cartilage and soft tissue. All tissues were assumed to be isotropic, homogenous and linear elastic. Trabecular bone on tissue level has an apparent elastic modulus of

500 - 5000 MPa.The material properties of single trabeculaeare similar to cortical bone [2]. In our model the single trabeculae were modeled.

Due to iterative nature of the solution procedure, prior to every iteration, the element material properties were updated. Each element consists of soft tissue, cartilage and bone. Material properties were similar to those by Shefelbine et al. (Table 1). To determine the elastic modulus and poisson’s ratio for each element, mixture rule was applied [2,5,9]:

(1)

(2)

where and are volume fractions and young modulus of cartilage, bone and soft tissue respectively. These volume fractions were fuzzy controller outputs and for each element there is a balance Equation:

(3)

In initial state, the model consisted of bone and soft tissue elements.

It was assumed that 3 MPa pressure loading was applied on the top of trabeculae while the other end of the trabeculae is fixed (Figure 1(a)) and soft tissue elements were not constrained.

To investigate the influence of load direction on trabecular structure a diagonal load at a 45˚ angle was applied on the top of one of the trabeculae while the end of trabecula on its opposite side was fixed. The load magnitude was similar to previous load case (Figure 1(b)).

2.2. Fuzzy Controller

In order to model biological processes during healing (including intramembranous ossification, cartilage formation, endochondral ossification and remodeling) a fuzzy logic controller, consisting of 18 linguistic rules, 5 inputs and 2 outputs was used. The most important characteristic of the fuzzy logic is its ability to model the qualitative processes, so for the modeling of the biological processes that are not easy to model quantitativelyfuzzy logic could be a useful tool.

The fuzzy controller input variables and their possible rangesare: 1) percentage of bone in the element = {low, nedium, high}, 2) percentage of cartilage in the element, 3) strain energy density = {low, , 4) maximum percentage of bone in neighborin gelements = {low}, 5) ratio = {low} (ratio is a remodeling index) and outputs of fuzzy controller are percentage of bone and percentage of cartilage in each element. Using balance equation, the percentage of soft tissue can be calculated.

The input of controller and mechanical stimulus in healing process was strain energy density:

(4)

where, is strain energy density and the components represent normal strains, shear strains, normal stresses and shear stresses respectively. was calculated by the finite element analysis and its membership functions were similar to diagrams of Ament et al. [9] (Figure 2). The fuzzy controller uses input from the finite element analysis and decides what type of tissue to build.

In healing process, after the fracture gap is stabilized, remodeling starts. During remodeling, the trabecular bone density changes. Bone density is related to its strain energy density and it could be used for simulation of remodeling process.

In the remodeling process, density variation changes

the elastic modulus of bone. Carter and Hayes [11] found that the elastic modulus and the strength of trabecular and cortical bone are closely related to cube and square of the apparent bone density respectively and the relationship between density (gr/cm³) and modulus of elasticity (MPa) is [12]:

(5)

where is bone apparent density. Note that in this simulation was calculated in each iteration applying mixture rule, so based on Equation (5) the elastic modulus can be used instead of bone apparent density [13]:

(6)

For the first time, a new parameter was introduced in our model as a remodeling index (based on modulus of elasticity):

(7)

where is the ideal value of elastic modulus for element j which is defined in the main program; for example is defined 10 Pa for bone elements and is the elastic modulus of element j in iteration i. The ideal values of elastic modulus are given in Table 1. For being smaller than 1, bone formation occurs whereas for S being larger than 1, bone resorption occurs (Figure 3).

Fuzzy rules were a set of IF/THEN statements (Table 2). Each rule was related to a biological process and all conditions must be met in order to activate the rules.In fuzzy logic controller only one rule could be activated or if the inputs were related to overlapping areas in diagrams (Figures 2 and 3), then multiple rules could be activated. For example, rule #2: IF cartilage is low and is increased and max bone neighbor is not low THEN increase bone, and rule #3: IF cartilage is low and is physiologic and max bone neighbor is not

low THEN increase bone. If was in overlapping area of “increased” and “physiologic” membership functions and other inputs satisfy all conditions of rules 2,3, then both rules will be activated.

After determination of outputs from the fuzzy controller, the material properties of each element in the finite element model will be updated and the next iteration begins (Figure 4). This process continues until the structure will not change and the state is stable.

3. Results

Figure 5 shows all the stages of trabecular bone fracture healing, bone formation, and remodeling. In the initial step of healing, after applying load and producing mechanical stimulus, based on the values of strain energy density and maximum percentage of bone in neighboring elements, fuzzy rules were activated. For example, if the elements were near the trabeculae spicules, then the intramembranous ossification rules would be activated. If the strain energy density was very high or percentage of bone in neighboring elements was low, then the cartilage formation and subsequently the endochondral ossification rules would be activated. For each iteration elastic modulus of elements were changed. After the fracture gap was filled with bone elements (Figure 5(b)), based on the magnitude of the S, remodeling rules were activated and bone resorption began, therefore the elastic modulus decreased in some elements. The simulation continued until a stable structure and satisfied conditions were