Inhomogeneous material property assignment and orientation definition of transverse isotropy of femur

The finite element method has been increasingly adopted to study the biomechanical behavior of biologic structures. Once the finite element mesh has been generated from CT data set, the assignment of bone tissue’s material properties to each element is a fundamental step in the generation of individualized or subject-specific finite element models. The aim of this work is to simulate the inhomogeneous and anisotropic material properties of femur using the finite element method. A program is developed to read a CT data set as well as the finite element mesh generated from it, and to assign to each element of the mesh the material properties derived from the bone tissue density at the element location. Moreover, for cancellous bone in femoral neck and cortical bone in femoral stem, the principal orientations of transverse isotropy were defined based on the trabecular structures and the haversian system respectively.


INTRODUCTION
The determination of the mechanical stresses in human bones is of great importance in both research and clinical practice.Since the stresses in bones cannot be measured non-invasively in vivo, an effective way to estimate them is through the finite element (FE) analysis which is widely used in academic research and clinical applications, such as the theory of bone remodeling [1], the design of prosthesis [2] and the evaluation of facture risk [3].In early period the methods used to get bone geometry and mechanical properties were inaccurate and sometimes highly invasive and destructive.It is well known that CT images can provide fairly accurate quantitative information on bone geometry based on high contrast between the bone tissue and the soft tissue around [4,5].Moreover, it has been demonstrated that CT numbers are almost linearly correlated with apparent density of biologic tissues [6].Good empirical relationships have been established experimentally between density and mechanical properties of bone tissues [7,8].
In early studies, if a generic bone model was constructed, then the mechanical properties of the different bone tissues were usually derived from average values reported in published experimental studies [9].For an individualized or subject-specific model, however, mechanical properties should be derived from CT data.It has been shown that the stress distribution of a bone is strongly related to the mechanical properties distribution in the bone tissue [10], and the mechanical properties of bone have been showed to depend on the subject, anatomic location, orientation, biological processes and time.Hence, it is important to find an effective method to properly map the material properties derived from CT data into finite element models.
The CT data can be regarded as a three-dimensional scalar field (related to the tissue density) sampled over a regular grid.If the finite element mesh is generated starting from the same data, the mesh and the density distribution are perfectly registered in space.The only problem is how to account for this inhomogeneous distribution of material properties into the FE model.Many approaches were proposed in literature to perform this task [11,12,13,14,15,16,17,18,19].However, these algorithms only simulated the inhomogeneity of bone, and the isotropic material property assignment was adopted without considering the orientation of material.Since the bone material is widely recognized as being anisotropic rather than isotropic [8,20], the FE simulation of inhomogeneity and isotropy cannot reflect the actual characteristics of bone structure.The aim of this study is to simulate the inhomogeneous and anisotropic material properties of femur with finite element method.

CT Data
The CT dataset of a man's femur is obtained from the public database which is created by VAKHUM project (http://www.ulb.ac.be/project/vakhum/index.html).The use of the data is free for academic purposes.The CT data are in standard DICOM formats.The slice thickness is 1mm in the epiphysis and 3mm in the diaphysis.

Finite Element Mesh
The finite element mesh of the right femur (Figure 1) generated from the corresponding CT dataset above is also obtained from the VAKHUM project.It is in a Patran Neutral file format.The mesh is made of linear hexahedral elements and is generated using the HEXAR (Cray Research, USA) automatic mesh generator that implements a grid-based meshing algorithm.The model mesh is spatially registered with the CT dataset.The complete finite element mesh consisted of 9,294 nodes and 7,934 elements.

Calculation of the Average CT Number ( HU )
For each element of the mesh, an average HU value is calculated with a numerical integration as follows [17]: where indicates the volume of the element n, (x,y,z) are the co-ordinates in the CT reference system, (r,s,t) are the local co-ordinates in the element reference system, and J represents the Jacobian of the transformation.The integrals in (1) are evaluated numerically, and the order of the numerical integration can be chosen by us.The value of HU(x,y,z) in a generic point of the CT domain is determined by a tri-linear interpolation between the eight adjacent grid points' values.n V

Calibration of the CT Dataset and
Simulation of the Inhomogeneity CT numbers are dependent on many factors related to the specific exam.It is assumed that the relationship between CT numbers and apparent density is linear.A calibration phantom with bone-equivalent (solution of hydroxyapatite) insertions of different densities, embedded in a water-equivalent resin-based plastic (The European Spine Phantom, [21]) was used to obtain the parameters of the linear regression.The calibration equation is then: where n  is the average density assigned to the element n of the mesh, n HU is the average CT number and  and  are the coefficients provided by calibration.
Referenced values for calibration are selected for approximate calibration from [16]: Radiographic and apparent density of water (0 HU, 1 g/cm 3 ); Average radiographic density in the cortical region and the apparent density value for cortical bone (1840 HU, 1.73 g/cm 3 ).
With the steps above, each element of the mesh is assigned different apparent density.Therefore, this procedure makes the simulation of inhomogeneity of femoral material come true.

The Density-Modulus Relationship and the
Simulation of Anisotropy Various empirical models of the relationship between Young's modulus and the apparent bone density can be found in the literatures [7,8].The relations are generally reported: where n E is the average Young's modulus assigned to the element n of the mesh, n  is its apparent density and a, b and c are the coefficients.In most of these works [11,12,13,14,15,16,17,18,19], the isotropic material property assignment was adopted without considering the direction of material property as its simplicity.Many previous studies had clearly demonstrated the anisotropic behavior of bone, and the determination of bone material properties as a function of direction is the essential requirement [22].The axial direction is defined according to the Haversian osteons of the cortical bone and according to the spatial main direction of the trabecular structures of cancellous bone.The relationship between bone material properties and apparent density is given by [8]: coincide with the structure of bone.
As we know, bone structure is customarily recognized as confirming to 'wolff's law' which is essentially the observation that bone changes its external shape and internal architecture in response to stresses acting on it.Thus, the structure of bone (or material orientation) strongly coincides with the principal stress track.The directions of the femoral neck and stem can be approximately regarded as stress track in femur.
 Cortical bone: Femoral bone can be recognized as transversely isotropic material.According to the cortical bone structure in femoral stem and cancellous bone structure in femoral neck, the principal material orientation of cancellous bone is defined by the direction of the trabecular structures and the principal material orientation of cortical bone by the direction of the haversian system.

RESULTS
where E is the Young's modulus (MPa), G the shear modulus (MPa),  is the apparent density (g/cm 3 ),  the Poisson's ratio.All parameters are set in elemental Cartesian coordinate system of Patran (MSC, USA).

Material Properties Distribution
This material assignment procedure produces 165 different material cards.The distribution of the material properties in femur are shown in Figure 2. The maximum and minimum values for apparent density and elastic modulus are listed in Table 1.The maximum values are corresponding to the Mat-1 and the minimum to Mat-165 as a result of the definition in the program.The numbers of elements for each material property card are showed in Figure 3.
This procedure may, theoretically, lead to a different material card (Mat) for each element of the mesh, which may result in computational problems with those codes that can handle a limited number of materials.We can reduce the number of materials by choosing a E  threshold.Then the maximum computed value of the elastic modulus, is assigned to the material Mat-1.All the elements with are assigned the material Mat-1.Mat-2 is characterized by the maximum E of the remaining elements and so on until the whole set of the model material is defined.In this paper, the threshold

The FE Simulation of Transversely Isotropic Material Property
Figure 2 only showed the inhomogeneous distribution of material properties.After separating the femoral stem and neck, the principal material orientations for elements in neck and stem are defined respectively.Figure 4 illustrates the vector form of the principal orientations in femoral neck.Figure 5 presents the vector form of the principal orientations in femoral stem.
A program was developed to read the original data including CT data, finite element mesh data and some parameters and then to implement the procedure described above.In this way, each element is assigned different density, elastic modulus and Poisson's ratio.At last, a finite element mesh provided with the assigned material properties is written out in a Patran Neutral file format.

DISCUSSION
Every element of the finite element model has its own element coordinate system used to determine the material's principal axe of this element.As the element coordinate systems vary from element to element, the principal material orientations need to be defined so as to As is well-known, it is of significance to explore the biomechanical behavior of bone.FE analysis has been adopted by many researchers based on CT data that can provide useful information on the geometrical topology     and material properties of bone.However, how to construct accurate FE models that reflect the real material properties of bone remains a problem.This work is aimed to simulate the inhomogeneity and anisotropy of the femur in order to make the FE material model more accurate.
As the relationship between density and Young's modulus is power law, different results will be produced when firstly averaging on each element the HU field and then calculating the element Young's modulus, or, on the contrary, averaging Young's modulus derived from HU value.It has been demonstrated that the latter can achieve more accurate results [18].In the present work, CT number over each element volume is averaged first and then the Young's modulus is calculated with the experiential relationships.However, it has no influences on the FE simulation of anisotropy.
Different algorithms are adopted to assign material properties into elements of the mesh model [17].Most of these algorithms are based on the assumption that the variability of the CT numbers within the volume of each element can be neglected.All differences rely on how the average value of each element is computed.
As we know, inner of the femur is marrow cavity that contains marrow or air.However, the mesh model in this study does not take into account the endosteal surface.Thus, the elements that belong to marrow cavity are assigned density and elastic modulus as well.Some simple methods are introduced to assign to each element the value by averaging the CT numbers of the CT sampling grids which are nearest to corresponding nodes of element or by averaging the CT numbers of the eight CT sampling points surrounding the element centroid [10,15].These methods may produce inaccurate results when the element size is significantly larger than the spacing of the CT sampling grid.FE simulation of complete anisotropy of femoral bone is impossible at present.Thus, femoral material is simplified as being transversely isotropic in this study.It has been demonstrated that the structure of femur is highly variable, especially to cancellous bone.Thus, a clear and exact definition of the principal axes of transversely isotropy is impossible.In this study, we only separated the femoral neck and stem.Then, the principal orientations of neck were defined on the basis of the direction of trabecular structure and the principal orientations of stem on the basis of the direction of harversian system.As the structure of femur (or material orientation) coincides with the track of principal stress, the orientation definition based on pass of stress is reasonable.
A second approach that can avoid those shortcomings determines all CT sampling points that fall inside the element volume and assigns to the element the average of these values [16].However, this method may not produce satisfactory results when the elements are of comparable dimension or smaller than the voxel size.
An advanced method adopted in our study estimates the average CT numbers by numerical integration over the element volume [17].In this case, the dimension of the finite element is not influent on the accuracy of the average CT number estimation since the calculation of CT numbers takes into account the spatial distribution of the CT number in the eight surrounding CT lattice vertices.
Muscle, as well as bone, is one of the most important parts of the biologic body.Musculoskeletal loading influences the stresses and strains within the human bone and thereby affects the processes of bone modeling and remodeling.Therefore, three-dimensional FEA for muscle is also of great importance.Assignment of the material properties to muscle FE mesh using the method in this paper may encounter some difficulties as the muscle belong to soft tissue that do not have steady shape.Thus, further research may be focused on the assignment of muscle material property.
Although the inhomogeneous and transversely isotropic material properties simulated in this work are theoretically close to that of real femur, experimental validation need be performed in the future.

Figure 2 .
Figure 2. Right femur with different material properties mapped on it, observed from different angles of view.

Figure 3 .
Figure 3.The number of elements for each material property card.

Figure 4 .
Figure 4.The principal material orientations in femoral neck showed in vector form (the purple arrows).

Figure 5 .
Figure 5.The principal material orientations in femoral stem showed in vector form (the purple arrows).

Table 1 .
Material properties of femur (the unit for density is g/cm 3 , and for elastic modulus is MPa).