Assessment of Hip Fracture Risk Using Cross-Section Strain Energy Determined by QCT-Based Finite Element Modeling

Accurate assessment of hip fracture risk is very important to prevent hip fracture and to monitor the effect of a treatment. A subject-specific QCT-based finite element model was constructed to assess hip fracture risk at the critical locations of femur during the single-leg stance and the sideways fall. The aim of this study was to improve the prediction of hip fracture risk by introducing a novel failure criterion to more accurately describe bone failure mechanism. Hip fracture risk index was defined using cross-section strain energy, which is able to integrate information of stresses, strains, and material properties affecting bone failure. It was found that the femoral neck and the intertrochanteric region have higher fracture risk than other parts of the femur, probably owing to the larger content of cancellous bone in these regions. The study results also suggested that women are more prone to hip fracture than men. The findings in this study have a good agreement with those clinical observations reported in the literature. The proposed hip fracture risk index based on strain energy has the potential of more accurate assessment of hip fracture risk. However, experimental validation should be conducted before its clinical applications.


Introduction
Two of the major determinants of proximal femur fractures among the elderly are osteoporosis and sideways fall. The most common serious injury associated with the fall of an elderly person is hip fracture. Furthermore, hip fracture is associated with an up to 20% chance of death, a 25% chance of long term institutionalization, and less than a 50% chance of full recovery [1]. The total number of hip fractures in men and women in 1990 was estimated to be 338,000 and 917,000, respectively, over the world [2]. Assuming no change in the age-and sex-specific incidence, the number of hip fractures is estimated to approximately double to 2.6 million by the year 2025 and 4.5 million by the year 2050 over the world [2]. By the increasing trend in hip fractures because of the aging of the population, the worldwide annual costs of hip fractures in the year 2050 have been estimated to be $131.5 billion [3]. So, hip fracture risk should be assessed in individuals who are at a risk for an osteoporotic hip fracture to provide proper plans to prevent future probable fractures.
Statistical models and bone mineral density (BMD) captured by Dual-Energy X-ray Absorptiometry (DXA) are used for in vivo osteoporotic fracture risk assessment [4,5]. However, their accuracy in assessment of fractures is limited; most osteoporotic fractures actually occur with BMD measurements that are above the conventional osteoporotic threshold [6]. Fracture Risk Assessment Tool (FRAX) is a tool to evaluate an individual's fracture probability in the next 10 years, adopted by the WHO in 2008 [7]. FRAX does not take into account fall-induced impact force that is critically important in the hip fracture risk assessment [8,9]. The main limitations of the FRAX include the following: it is a statistical model and fracture risk is not consistent within 10 years with some of the treatment results [10]. Hip structure analysis (HSA) program is now commercially available and is used to automatically assess the geometric and structural parameters of the femur. Whereas HSA is based on a beam model, the deformation of bone is oversimplified, especially in the femoral neck and the intertrochanteric region where osteoporotic femoral fractures most often occur; therefore  the deformation is too complicated to be described as a beam model [11]. Finite element models constructed from quantitative computed tomography (QCT) are very helpful in assessing hip fracture risk, as they are based on well-established biomechanical principles and theories. In QCT-based finite element modeling, choosing a proper failure criterion is very important to accurate assessment of hip fracture risk. Commonly used bone failure criteria include von Mises stress and strain criteria [12][13][14][15], maximum principle stress and strain criteria [16][17][18][19], maximum shear stress criterion [20], and maximum distortion energy criterion [20]. The femur consists of inhomogeneous (porous) cancellous bone and nearly homogenous cortical bone, and their failure mechanisms are different due to their different microstructure. Failure mechanism of the cancellous bone is mostly in the form of spicule (trabecula) buckling, and the failure of denser cancellous bone and the cortical bone is mostly characterized by local cracking [21,22]. Cortical and cancellous bone have different properties. Cortical bone is usually more brittle [23], while cancellous bone is more ductile. The current failure criteria are not convenient in describing their properties. Therefore, the total strain energy, which integrates all information of stress, strain, and material properties, is a better option [21]. Mirzaei et al. [21,24] predicted failure strength and failure patterns of human proximal femur and human vertebrae using the strain energy criterion with a QCT-based FE model. Their predictions of the failure loads and failure locations were in a good agreement with the experimental findings. The strain energy criterion is widely used in fracture analysis of engineering materials. It is usually used in crack problems [25][26][27], composite laminates [28,29], and bone cement analysis [30]. Therefore, computation of hip fracture risk index (FRI) over a region of interest (ROI) based on the strain energy criterion theoretically should be more accurate for assessing hip fracture risk. To the best of our knowledge, there are currently no published studies that use the strain energy criterion for the hip fracture risk assessment. The objective of this study is to improve the hip fracture risk assessment procedure previously developed by Luo et al. [11] by introducing the strain energy criterion.

Materials and Methods
The proposed methodology for assessment of hip fracture risk in the critical regions of femur using the strain energy criterion determined from QCT-based finite element model is shown in Figure 1. The procedure is explained in detail in the following.

QCT-Based Finite Element Model
The purpose of this study is to accurately assess hip fracture risk, so, a 3D finite element model of subject's femur is required to achieve it. The 3D model can be constructed from the subject's femur QCT images. QCT slices are produced using multiple scanners with a set of proper acquisition and reconstruction parameters (Figure 2(a)). Slice thickness of 1 mm is commonly used. The scanned QCT images are stored in the format of Digital Imaging and Communications in Medicine (DICOM), which can be used for the construction of a 3D FE model. A proper segmentation is done to separate the femur for constructing the 3D model. Each voxel in the QCT-scan has an intensity (or grey scale) that is expressed as Hounsfield Unit (HU), which is correlated to bone density [31,32]. Threshold value   Table 1.

Generation of Finite Element Mesh.
In the first step, the geometrical model of the femur is generated from clinical QCT images using Mimics (Materialise, Leuven, Belgium). QCT images (in DICOM format) are imported to Mimics for segmentation (Figure 2(a)) and construction of 3D geometric model of the femur (Figure 2(b)). With the 3D geometric model, a FE mesh is generated using the 3matic module in Mimics (Figure 2(c)). The 4-node linear tetrahedral element SOLID72 in ANSYS was used in this study. To investigate model convergence, FE models with different maximum element edge lengths were created. For each FE model, displacement was calculated under the same loading and boundary conditions. The maximum element edge length that produced converged finite element solutions was obtained and used in all the rest of FE simulations. The number of elements for each case is assigned based on the maximum element edge length that provided converged FE solution, meaning that the number of elements is different from case to case.

Assignment of Material Properties.
To construct a more faithful FE model, bone material properties are considered inhomogeneous and isotropic in this study. Information on the inhomogeneous isotropic mechanical properties of the bone can be derived from the CT data using a mathematical relationship between the CT numbers and the mechanical properties of bone. The following empirical equation was used to determine bone ash density ( ash ) from HU number [33,34]: Equations (2) and (3), derived by Keller [35], were, respectively, used to assign Young's modulus ( ) and yield stress ( ) according to the bone ash density: = 116 2.03 ash (MPa) .
A constant Poisson's ratio (] = 0.4) was considered as suggested in the literature [13,36,37]. To assign material properties, elements are grouped into a number of discrete material bins using Mimics (Materialise, Leuven, Belgium), to approximately represent the continuous distribution of the inhomogeneous bone mechanical properties. To determine the maximum number of material bins, convergence study was performed. Models with different material bins were created for convergence study. For each FE model, displacement was calculated under the same loading and boundary conditions. The maximum number of material bins that generated converged finite element solutions was obtained. Figure 2(d) shows the isotropic inhomogeneous distribution of material properties.

Finite Element Analysis
Using ANSYS. The finite element model of femur with the assigned material properties was output from Mimics and then imported to ANSYS for finite element analysis. For a precise assessment of hip fracture risk during the single-leg stance and the sideways fall, loading and boundary conditions simulating the single-leg stance and the sideways fall configurations are required in the FE model. To simulate the single-leg stance configuration, 2.5 times the patient's body weight was applied as a distributed load on the femoral head in direction of femoral shaft axis [38] and femur was fixed at the distal end [13,39]; see Figure 2(e). Consider where is the subject's body weight in Newton (N). To simulate sideways fall, the distal end of femur was completely fixed and the surface of femoral head was fixed in the loading direction ( Figure 2(f)) [40,41]. The impact force during the sideways fall acting on the greater trochanter perpendicularly (Figure 2(f)) is given by [38,42]: where ℎ is the height of the subject in centimeter (cm). Loading and boundary conditions on the greater trochanter, the femoral head, and the distal end of femur were applied to a group of nodes using APDL codes (Figures 2(e) and 2(f)). After importing the QCT-based FE model and applying the loading and boundary conditions, finite element analysis was performed and finite element solutions were obtained. In all the analysis, the nodal displacements, stresses, and strains were obtained for each subject.

Detection of the Three Critical Cross-Sections on the Femur.
Hip fractures usually occur at one of the anatomical locations: the femoral neck, the intertrochanter, and the subtrochanter as illustrated in Figure 3. According to clinical observations, 49 percent of hip fractures are intertrochanteric, 37 percent are at femoral neck, and 14 percent are subtrochanteric [43]. Therefore, the smallest femoral neck cross-section (SFN CS), the intertrochanteric cross-section (IntT CS), and the subtrochanteric cross-section (SubT CS) are the three critical cross-sections of femur that usually have the highest fracture risk ( Figure 3). To determine the smallest femoral neck cross-section and the intertrochanteric cross-section, neckshaft angle is needed. The neck-shaft angle is the angle  between the femoral neck axis and the femoral shaft axis. This angle traditionally is measured on conventional radiography images, or using 2D images projected from CT/MRI data. In spite of their popularity, these methods are based on oversimplification of the real 3D anatomy and may lead to large errors due to the inaccuracies in selection of the measurement plane [44][45][46]. In this study, the neck-shaft angle was measured using a 3D measurement technique based on fitting functions. In this technique, the shapes of particular parts of the femur are approximated using geometric entities such as circle, cylinder, and sphere, which are well fitted to the actual anatomy, and the geometrical relationships between these entities are obtained to estimate the neck-shaft angle.
First, a sphere was fitted to the femoral head, to obtain the position of the joint's centre of rotation, which is also the femoral head centre. Then, the femoral neck axis and the femoral shaft axis are identified by applying the "fit ruled surface direction" function on the femoral neck and shaft. All fitting functions were applied using the 3-matic module in Mimics. The neck-shaft angle was also measured by 3matic module of Mimics ( Figure 4). With the femoral neckshaft angle, the intertrochanteric cross-section and the smallest femoral neck cross-section were found using in-house computer codes. The smallest femoral neck cross-section is chosen as the cross-section with the smallest area in the neck region and the intertrochanteric cross-section is chosen as the cross-section that has the largest area in the intertrochanteric region [47]. By using APDL codes, perpendicular planes on the femoral neck axis were determined and then areas of the cross-sections were calculated. The planes with the smallest and the largest areas were chosen, respectively, as the smallest femoral neck cross-section and the intertrochanteric crosssection. The subtrochanteric cross-section is considered five centimeters below the lesser trochanter [48] (Figure 3).

Hip Fracture Risk Index Definition Using the Strain Energy
Criterion. Based on the previous discussion of the bone failure mechanism and microstructure, the strain energy criterion is theoretically more favourable for hip fracture risk assessment. The strain energy at the three critical cross-sections of femur induced by the applied forces was computed using in-house developed MATLAB codes and the data extracted by APDL codes from the obtained finite element solutions. The plane boundaries of the three critical cross-sections, extracted from the finite element mesh, were imported to MATALB to generate a 2D mesh for calculating the cross-section strain energy. Figure 5 shows the generated triangle elements over the smallest femoral neck cross-section, the intertrochanteric cross-section, and the subtrochanteric cross-section. The strain energy at the three critical cross-sections induced by the applied forces is the sum of strain energy in all the triangle elements; that is, where is the strain energy in Element induced by the applied forces and is the number of triangle elements created over the concerned cross-section. Gaussian integration method was used to calculate the strain energy in elements. Integration points in each triangle element were determined using in-house MATLAB codes. By using Gaussian integration method, the strain energy of Element induced by the applied forces is calculated as wherêis the strain energy density of Element ;̂is the strain energy density at the integration point of in Element ; is the weight at the integration point; | | is the determinant of the Jacobean matrix of the triangle element; and is the number of integration points over the triangle element. The strain energy density at an integration point ( ) was determined from the finite element solutions obtained by the QCT-based FE model; that is, where { } is the displacement vector consisting of displacements at the nodes of Element ; matrix [ ] consists of the derivatives of shape functions of the element; and [ ] is the material property matrix. Consider where Poisson's ratio is constant (] = 0.4) and Young's modulus is a function of the bone density as given in (2). For each integration point, its Young's modulus is determined by the bone density at the point. The maximum allowable strain energy over a critical cross-section of the femur was also computed from the obtained finite element solutions using in-house MATLAB codes. The maximum allowable strain energy (or the yield strain energy) over a critical cross-section is obtained as where is the yield strain energy in Element and is the number of triangle elements over the concerned crosssection. The Gaussian integration method was also used to calculate the maximum allowable strain energy in each triangle element. The maximum allowable strain energy that a triangle element ( ) can sustain is given by wherêis the yield strain energy density in Element ; is the number of integration points; and̂is the yield strain energy density at integration point and is calculated aŝ where and are, respectively, the yield stress and Young's modulus at the integration point. Both of them are function of bone density, as given in (2) and (3).
Hip fracture risk index ( ) at a critical cross-section is defined as where and are, respectively, obtained from (6) and (11).

Element Size in Femur Finite Element Analysis.
The convergence of finite element solutions in a representative case is shown in Figure 6. The convergence study showed that the finite element displacements converged with the maximum element edge length smaller than 8 mm. Therefore, in construction of the rest of femur FE models, the maximum element edge length was set to 8 mm.  neck cross-section was monitored under the same loading and boundary conditions. Data on the displacement were compared among the FE models with different material bins. The results of the convergence study showed that the displacement did not change significantly with the number of material bins larger than 50. Therefore, in the assignment of material properties for all the cases, 50 discrete material bins were considered.

Element Size in Calculating Cross-Section Strain Energy.
Convergence study was also performed to determine the element size used in integrating cross-section strain energy, as it affects the calculated fracture risk index (FRI). The FRI at the concerned cross-section was calculated with different maximum element edge lengths. The results are plotted in Figure 7. The FRI did not change significantly with the maximum element edge length smaller than 5 mm. Therefore, the maximum element edge length was set to 5 mm in calculating cross-section strain energy.

The Number of Integration Points in Calculating Cross-Section Strain
Energy. The effect of the number of integration points on the calculated FRI was investigated. FRI at the smallest femoral neck cross-section was computed for 5 clinical cases with different number of integration points. The relative errors between FRIs obtained with 3 and 7 integration points are shown in Table 2. As it can be seen, the errors are not significant. Therefore, the 3-point integration rule was used in this study to reduce computational time.

Stress and Strain Patterns at the Three Critical Cross-
Sections. For the 10 clinical cases (5 females and 5 males, totally 20 right and left femurs), the maximum von Mises stress and strain at the three critical cross-sections of femur during both the single-leg stance and the sideways fall are shown in Figures 8 and 9. It can be observed that, during the sideways fall, the femoral neck and the intertrochanteric region received higher stresses than the subtrochanteric region (Table 4). But during the single-leg stance, the patterns in the stresses are different (Table 3); first, the differences between the stresses over the three regions are much smaller; for some cases, the stresses at the subtrochanteric region are higher than those in the other two regions (Figure 8). Strains in the three regions have similar patterns during both the single-leg stance and the sideways fall (Tables 5 and 6).

Comparison of Hip Fracture Risk at the Three Critical
Cross-Sections. For the 60 clinical cases (30 females and 30 males), hip fracture risk indices based on the strain energy    As shown in Tables 7, 8, 9, and 10 and Figures 12 and 13, the average FRI at the smallest femoral neck was higher than those at the intertrochanteric and the subtrochanteric crosssection during both the single-leg stance and sideways fall.
For the single-leg stance, FRIs for all cases at the three critical cross-sections are much lower than one (FRI ≪ 1), indicating that the possibility of hip fracture incidence in the single-leg stance was low. For the sideways fall, FRIs of 8 right femurs and 7 left femurs at the smallest femoral neck crosssection and FRIs of 7 right femurs and 5 left femurs at the intertrochanteric cross-section were higher than one (FRI > 1), meaning that there is possibility for fracture occurring in these regions; but the FRIs at the subtrochanteric crosssection in all cases were much lower than one (FRI ≪ 1), indicating that there is lower possibility of fracture in this region. Figure 14 shows the number of possible fractures at the three critical cross-sections of femur during the sideways fall, that is, the cases that have FRI larger than one (FRI > 1) at one of the three critical cross-sections.  1  3  5  7  9  11  13  15  17  19  21  23  25  27  29  31  33  35  37  39  41  43  45  47  49  51  53  55 1  3  5  7  9  11  13  15  17  19  21  23  25  27  29  31  33  35  37  39  41  43  45  47  49  51  53  55

Comparison of Hip Fracture Risk in Women and Men.
In this study, hip fracture risk for the 30 females and 30 males was assessed. The average hip fracture risk of females was generally higher than that of males. As it can be seen from Tables 11-13 and Figures 15 and 16, the average FRI at the smallest femoral neck cross-section, the intertrochanteric cross-section, and the subtrochanteric cross-section of the 30 females is higher than that of the 30 males for both the single-leg stance and the sideways fall. Figure 17 shows the number of possible hip fractures in women and men during the sideways fall for the studied cases.

Discussion
Hip fracture may occur anywhere from the articular cartilage of the hip joint to the femur shaft [48]. Not only the location of the fracture types but also the etiology differs. It was reported that women with the intertrochanteric fracture have significantly lower BMD than those with the femoral neck fracture [49][50][51]. On the other hand, the femoral neck fracture     may not be mainly attributed to low BMD but may be related to external causes such as sideways fall [52]. Femoral neck and intertrochanteric fractures are often the result of falls from standing height and impact onto the greater trochanter, particularly for the elderly. The subtrochanteric fractures, on the other hand, are typically the result of high energy impacts such as motor vehicle accidents and falls from a height [53]. The selection of bone failure criterion is challenging. In literatures, stress and strain based failure criteria such as   the von Mises stress and strain criteria and the maximum principle stress and strain criteria were commonly used to assess hip fracture risk. To the best of our knowledge, the strain energy based failure criterion has not been used yet for hip fracture risk assessment. Whereas the cancellous bone failure is in the form of buckling and deformation (strain intensity) and the cortical bone failure is related to its local cracking (stress intensity), strain energy failure criterion, which is a combination of both stress and strain intensities, is theoretically more reasonable than other failure criteria for hip fracture risk assessment. The differences between the strains in the three critical regions of femur during both the single-leg stance and the sideways fall (Figure 9, Tables  5 and 6) are much higher than the differences between   the corresponding stresses (Figure 8, Tables 3 and 4), indicating that bone failure is more sensitive to the strains because of its fragility property and the effects of strains should also be considered in bone fracture risk assessment.
Results of this study show that the femoral neck and the intertrochanteric region have higher fracture risk than the subtrochanteric region (Figures 12 and 13 and Tables 7-10), which is consistent with the fact that the femoral neck and the intertrochanteric region have a larger proportion of cancellous bone than the subtrochanteric region; and the cancellous bone is generally weaker than the cortical bone. Therefore, hip fracture is most likely to initiate first at the femoral neck and then in the intertrochanteric region or in the subtrochanteric region. For all subjects, there is very low hip fracture risk during the single-leg stance (FRI < 1). Among all the cases, during the sideways fall, 15 femurs at the femoral neck and 12 femurs at the intertrochanteric region have FRI higher than one (FRI > 1), while there is very low fracture risk at the subtrochanteric region for all the subjects (FRI < 1) (Figure 14). Our findings have a good agreement with the previous clinical observations. According to the clinical observations, the femoral neck is the most common location for a hip fracture, accounting for 45% to 53% of hip fractures; the intertrochanteric fractures account for approximately 38% to 50% of all hip fractures; and the subtrochanteric fractures are less common than the femoral neck and intertrochanteric fractures, accounting for approximately 5% to 15% of hip fractures [54,55]. Also it was observed that, in old people with age between 65 and 99, femoral neck and intertrochanteric fractures occurred with approximately the same frequency [56]. Data on 176 geriatric patients with hip fractures showed that 59% were the intertrochanteric fractures while the rest (41%) were the intracapsular neck fractures [57]. In another study by Michelson et al. [43], it was observed that 37% of hip fractures are in the femoral neck, 49% are intertrochanteric, and 14% are subtrochanteric. Therefore, the femoral neck and the intertrochanteric region are more disposed to fracture than the subtrochanteric region.
Both the fracture risk level and the potential fracture location are patient-dependent and depend on BMD and other conditions that have not been studied in this paper. We have found that women generally have lower bone strength than men and thus are exposed to higher hip fracture risk during both the single-leg stance and the sideways fall (Figures 15 and 16 and Tables 11-13). Based on the calculated FRI, in the 30 women potentially there are 12 femoral neck fractures and 7 intertrochanteric fractures, while for the 30 men there are only 3 potential femoral neck and 5 potential intertrochanteric fractures ( Figure 17). Our finding is consistent with other studies that have found that hip fracture risk in women is generally higher than in men. The results of a worldwide study by Dhanwal et al. [58] on incidence and epidemiology of hip fractures in Asia, Africa, Europe, Latin America, North America, and Oceania show that women are more disposed to the hip fracture risk than men in different countries over the world. In the study by Jacobsen et al. [59], the age-adjusted hip fracture incidence rates in white US males and females were 4.3 and 8.1 per 1,000 per year, respectively. Higher hip fracture risk in women may be related to their lower BMD.

Conclusion
A reliable methodology to assess hip fracture risk in individuals is crucially important for preventing hip fracture and initiating a treatment. The purpose of this study is to propose a more effective hip fracture risk index that is based on the strain energy failure criterion, and it is able to better describe bone failure mechanism and microstructure. The proposed fracture risk index can predict not only the fracture risk level, but also the potential fracture location. The results of this study showed that there is a very low hip fracture risk during the single-leg stance, while, during the sideways fall, there is a high fracture risk at the femoral neck and the intertrochanteric region, compared to the subtrochanteric region. Based on the results obtained from this study, women are more prone to hip fracture than men. The procedure described in this study can be implemented into computer programming and used in hip fracture prevention and monitoring of osteoporosis treatments in the elderly. The method may also help design more effective hip protectors by providing feedback information such as stresses/strains in the hip for adjusting the design parameters. The main limitation of this study is that no experiment has been conducted to validate the predicted fracture risk levels and potential fracture locations, which has been set as future work.