Energy Based Simulation of Trabecular Bone Fracture Healing Using Finite Element and Fuzzy Logic

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.


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.

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]: where The input of controller and mechanical stimulus in healing process was strain energy density: 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.

 
T , , , , , 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)).
where,   is strain energy density and the components i ij ,     represent normal strains, shear strains, normal stresses and shear stresses respectively.  u x 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.

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.
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 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,   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 3 ) and modulus of elasticity (MPa) is [12]: where a  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]: For the first time, a new parameter was introduced in our model as a remodeling index (based on modulus of elasticity): where refj is the ideal value of elastic modulus for element j which is defined in the main program; for example refj 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   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.

  u x
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.

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  achieved (Figure 5(c)).Applying diagonal force resulted in diagonal trabeculae (Figure 6).Finally, the elements in trabeculae reached to the elastic modulus of 9 GPa.

Discussion
We proposed a model to simulate trabecular bone fracture healing and remodeling, using finite element analysis and fuzzy logic.This modeling was mainly focused on mechanical factors incorporating physiological fuzzy rules.Increasing the fracture gap size had no considerable effect on final structure of trabecular bone, but increased the solution time to reach a stable state.It was shown that decreasing the fracture gap size, decreased the solution time.Similar to previous model [2], loading conditions were simplified.Applying vertical force on trabeculae resulted in vertical trabeculae and applying diagonal force resulted in diagonal trabeculae.Therefore, trabecular structure is affected by loading direction.Eventually, in biological situations where multi-directional loads are present, the trabecular bone structure is more complex [2].
The biological processes consisting of intramembranous ossification, cartilage formation, endochondral ossification and remodeling were simulated using a fuzzy logic controller.Fuzzy logic is advantageous, since one can include biological factors/variations into the biomechanical modeling of the healing process [9].
In previous models of diaphyseal fracture healing, stress, strain, fluid velocity and strain energy were used as mechanical stimuli to healing.Shefelbine et al. [2] used hydrostatic strain and octahedral shear strain as mechanical stimuli.In our model strain energy density was used.Weinans et al. [14] proved that models using strain energy or stress as the remodeling stimulus create a positive feedback system.Due to better performance and stability of FEA solution, we used quadratic elements in the model.Strain energy density,   U x , has two main advantages: 1) it is easier to deal with since it is a scalar value [9], and 2) using strain energy density means that the effects of all components of strains and stresses are considered in modeling the healing process.
An improvement of the model has been achieved by introducing a new parameter as a remodeling index, S. Figure 7 shows the final result for the model without considering S. It is obvious that only bone mass was changed and there is no indication of the directional changes of trabecular structure.This parameter is a mechanical factor and is related to apparent density of tissue.The parameter S is independent of the geometry and the loading conditions, and using it in fuzzy rules helps toward more realistic simulation of remodeling process of trabecular bone.
Many experiments on skeletal repair have been per-  formed that included wide variety of factors, such as: biological, mechanical, hormonal, sex, age, etc. [15].In our model we focused on mechanical factors and neglected vascularity, growth factors and other biological events.n future investigations, the effects of biological events and fluid flow should be taken into considerations to further improve the validity of the model.

Figure 1 .
Figure 1.Finite element model of fracture gap with loading and boundary conditions.(a) pressure loading, (b) diagonal loading.

Figure 2 .
Figure 2. Membership function of strain energy density.

Figure 3 .
Figure 3. Membership function of S.

Figure 4 .
Figure5shows 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 (Figure5(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

Figure 7 .
Figure 7. Final result for the model without considering S.

Table 1 . Material properties of tissues.
and max bone neighbor is not low THEN increase bone, and rule #3: IF cartilage is low and