- Email: [email protected]
Analysis of strength and failure pattern of human proximal femur using quantitative computed tomography (QCT)-based finite element method
Bone 64 (2014) 108–114
Contents lists available at ScienceDirect
Bone journal homepage: www.elsevier.com/locate/bone
Original Full Length Article
Analysis of strength and failure pattern of human proximal femur using quantitative computed tomography (QCT)-based finite element method Majid Mirzaei ⁎, Maziyar Keshavarzian, Vahid Naeini Department of Mechanical Engineering, Tarbiat Modares University, P.O. Box: 14115-143, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 14 December 2013 Revised 10 March 2014 Accepted 4 April 2014 Available online 13 April 2014 Edited by Richard Eastell Keywords: Bone Fracture Imaging Noninvasive methods Strain energy density
a b s t r a c t This paper presents a novel method for fast and reliable prediction of the failure strength of human proximal femur, using the quantitative computed tomography (QCT)-based linear finite element analysis (FEA). Ten fresh frozen human femora (age: 34 ± 16) were QCT-scanned and the pertinent 3D voxel-based finite element models were constructed. A specially-designed holding frame was used to define and maintain a unique geometrical reference system for both FEA and in-vitro mechanical testing. The analyses and tests were carried out at 8 different loading orientations. A new scheme was developed for assortment of the element risk factor (defined as the ratio of the strain energy density to the yield strain energy for each element) and implemented for the prediction of the failure strength. The predicted and observed failure patterns were in correspondence, and the FEA predictions of the failure loads were in very good agreement with the experimental results (R2 = 0.86, slope = 0.96, p b 0.01). The average computational time was 5 min (on a regular desktop personal computer) for an average element number of 197,000. Noting that the run-time for a similar nonlinear model is about 8 h, it was concluded that the proposed linear scheme is overwhelmingly efficient in terms of computational costs. Thus, it can efficiently be used to predict the femoral failure strength with the same accuracy of similar nonlinear models. © 2014 Elsevier Inc. All rights reserved.
Introduction Hip fractures cause significant numbers of disability and mortality worldwide. The incidence of these fractures varies globally, e.g., the city of Shiraz in Iran is reported to have the highest incidence of hip fractures in Asia [1]. Previous studies estimated that the worldwide number of hip fractures was 1.3 million in 1990 and predicted an increase between 7.3 and 21.3 million by 2050 [1–3]. Thus, many researchers have tried to devise new methods for noninvasive assessment of the hip fracture risk, among which the QCT-based finite element analysis (FEA) has been recognized as a very reliable tool [4–12]. However, there are some limitations that should be resolved before this method can be clinically feasible. One major problem is the high computational time and cost for implementation of nonlinear FEA which has been shown to have the required accuracy [4,6–9,13]. A typical QCT-based nonlinear FEA of human proximal femur requires at least 8–10 h of runtime, even by implementing advanced computational techniques [4,7,9,13]. On the other hand, the linear FE analyses are very much faster, quite easier to apply, and can provide comparable predictions of the failure location and pattern for the proximal femur and vertebra [13,14]. However, they require further data processing schemes for predictions of the failure load that may give less accurate results in certain cases. ⁎ Corresponding author. Fax: +98 21 82884909. E-mail addresses: [email protected], [email protected] (M. Mirzaei).
http://dx.doi.org/10.1016/j.bone.2014.04.007 8756-3282/© 2014 Elsevier Inc. All rights reserved.
For instance, Keyak et al. proposed a linear FEA method for the prediction of the failure load in the sideways fall and stance configurations and argued that failure of a certain number of contiguous non-surface elements can represent the failure of the whole proximal femur [15]. On the other hand, Nishiyama et al. defined a constant yield strain level and suggested that the failure of a constant percentage of elements would represent the failure of proximal femur in the sideways fall configuration [16]. Nevertheless, we believe that the percentage or the number of critical elements should naturally be different for various specimens and loading orientations. Moreover, our previous studies on human vertebrae and proximal femur have shown favorable results in the prediction of the failure patterns using a strain energy density measure within the framework of linear QCT-based FEA [13,14]. Thus, this study was aimed at development and experimental verification of a new method for fast and accurate prediction of the failure strength of human proximal femur through specific assortments and usage of high strain energy elements of the linear QCT-based FEA results. Materials and methods Sample preparation Ten fresh-frozen human femora, from 9 cadavers (5 male, 4 female, average age: 34 ± 16), were used in this study (see Table 1) (specimens 9 and 10 were a contralateral pair from a single cadaver).
M. Mirzaei et al. / Bone 64 (2014) 108–114
plane (the femoral coronal plane) using three different views, i.e., the medial, the axial top, and the axial bottom views as follows (see Fig. 1(c)):
Table 1 Age, sex, side, and loading orientation of specimens. Number
1 2 3 4 5 6 7 8 9 10
Age
45 50 24 25 44 35 37 40 20 20
Sex
M M F F M F F M M M
109
Side Right/left
Orientation α
β
L R R L L R L R R L
10 10 10 15 15 0 20 −20 0 0
0 0 0 30 0 −28 0 0 −20 20
The initial treatment of the specimens included standard excision, freezing, and bagging procedures by the Iranian Tissue Bank (ITB). Moreover, the death causes announced by ITB were the three categories of coronary failure, cerebral death, and fatal trauma. Nevertheless, the plain radiographs of the samples were examined to ascertain the lack of metastatic diseases, pathologic defects, or insufficiency fractures. Definition and implementation of the femoral reference system The orientations of mechanical loadings on proximal femurs are usually adjusted with respect to the coronal plane, which is defined as a plane containing the femoral cervical and shaft axes. Nevertheless, certain ambiguities exist about recognition of these axes on femoral specimens, because the complex geometries of femoral parts lack any axisymmetric feature. Moreover, the lack of a unique definition of these axes for both FEA and mechanical testing can be troublesome. We used the reference system proposed by Mirzaei et al. [13], in which 3 distinct points are marked on the femoral surface to establish a reference
• The center point of the femoral head from the medial view (point 1). • The center point of the femoral head from the axial top view (point 2). • The center point of the femoral shaft from the axial bottom view (point 3). Despite the fact that the femoral head is not quite spherical and the femoral shaft is not quite cylindrical, the three points can be marked without ambiguity by circumscription of the related views in quadrilaterals as depicted in Fig. 1(c). In practice, a specially-designed holding frame was used to designate these points on the samples as depicted in Fig. 1(b). The frame has two side plates, which move simultaneously in opposite directions using special connecting rods so that, upon contact with the femoral head, the median plane of the head coincides with the frame plane. The frame also has a two-pin slider for alignment of the femoral shaft and three screws which are used to specify both the frame plane and the femoral coronal plane. The role of the forth screw is just to add extra fixity. This reference frame remains attached to the sample during the QCT process and is naturally transferred to the 3D and the FE models. On the other hand, the role of the Plexiglas disk attached to the lower end of the distal femur is to preserve and transfer the reference system to the mechanical testing. The disk has three screws which are used to bring it to the desired location and orientation with respect to both the femoral shaft and the holding frame. In its final position, the disk plane is parallel to the lower side of the holding frame, its axis of rotation passes through point 3, and the direction of its long screw lies in the defined coronal plane. Thus, when the femur is detached from the holding frame, the disk preserves the location and orientation of the defined coronal plane.
Fig. 1. (a) A proximal femur specimen. (b) The specimen inside the holding frame. (c) Three different views of a proximal femur and the three points identified on the femoral surface for definition of the coronal plane. (d) The imaging and calibration phantoms inside the scanner. (e) The CT image of a section of the specimen and calibration rods. (f) The spherical coordinate system used in this study. (g) Conventional coordinate system for loading of femur; α (the angle between the applied load and the sagittal plane) and β (the angle between the load and the coronal plane). (h) The 3D solid model of the specimen. (i) The FE model of the specimen.
110
M. Mirzaei et al. / Bone 64 (2014) 108–114
QCT scanning and image processing Each framed specimen was placed inside a Plexiglas container filled with water, designed as an imaging phantom for QCT studies (Fig. 1(d)). The QCT scans were carried out using a clinical scanner (SIEMENSSOMATOM 64, 140 kV, 80 mAs, 0.5 × 0.5 mm/pixel resolution, and 1 mm slice thickness). Using a calibration phantom (MINDWAYS Software, Inc., San Francisco, CA) (Figs. 1(d,e)), gray scale values were mapped to K2HPO4 equivalent density (ρKHP) using five tubes with reference densities and the Hounsfield Units (HUs) were calibrated. The segmentation of the bone hard tissue from the surroundings was performed for each slice, using a homemade computer code that implements the MATLAB's Image Processing Toolbox [14]. The raw QCT images (in DICOM format) were converted into binary format and a combination of user-defined threshold limit with an edgefollowing scheme was used for generation of the hard-tissue contours and elimination of the soft tissue.
Finite element analysis A complementary homemade MATLAB code was used to build the FE models by conversion of each voxel into an 8-noded brick element. In brief, this code treats each voxel as a basis for construction of a single brick element and prepares an input file (node and element file) which is compatible with the format of the ANSYS software (ANSYS. Inc., Canonsburg, PA) (see Fig. 1(i)). Calibrated CT scan data [17] (ρash = 1.22ρKHP + 0.0526) were used to calculate the ash density (ρash, g/cm3). The bone tissue was assumed as an inhomogeneous isotropic material with linearly elastic properties for linear analyses and linearly elastic-perfectly plastic behavior for nonlinear analyses. Young's modulus of elasticity (E, MPa) and yield strength (S, MPa) values were assigned to each element using the empirical relationships based on calibrated ash density for the trabecular and the cortical bone [18–20]; E = 33,900ρash2.20 for ρash ≤ 0.27 (trabecular bone), E = 10,200ρash2.01 for ρash ≥ 0.6 (cortical bone), E = 5,307ρash + 469 for 0.27 b ρash b 0.6 (transition), S = 137ρash1.88 for ρash b 0.317 (trabecular bone), S = 114ρash1.72 for ρash ≥ 0.317 (cortical bone). The three reference points were identified from CT scans and the coronal plane was recognized and used for the application of the load and boundary conditions. However, definition of a general coordinate system for application of mechanical loads on femur is another important issue in the FEA of femoral samples. Conventionally, two angles are defined and used for loading of the femur: α (the angle between the applied load and the sagittal plane) and β (the angle between the
load and the coronal plane) (see Fig. 1(g)). In this study, we used a spherical coordinate system on the femoral head as depicted in Fig. 1(f). The following relationship exists between the conventional α–β and the new ϕ–θ systems: tan θ ¼ tan β= tan α; tan ϕ ¼ tan α= cos θ:
ð1Þ
The boundary conditions in the finite element model were consistent with the experiments, i.e., the load was equally distributed among the nodes of femoral head (which were in contact with the 30-mm steel cap) and the nodes on the clamped cross section of the shaft were fully restrained (Fig. 2). All elements with the modulus below 5 MPa were assigned a low modulus of 0.01 MPa [20]. The elements loaded on the femoral head were identified and assigned an elastic modulus of 20 GPa and a yield strength of 200 MPa to prevent element distortion [5]. The Poisson's ratio was assumed to be 0.4 [20]. The models were loaded at 8 different orientations (for replication of physiologic and non-physiological loadings) to test the validity of models under various loading conditions (see Table 1). The linear and nonlinear FE models were analyzed using ANSYS software (ANSYS. Inc., Canonsburg, PA) on a desktop PC (2 × 2.7 GHz Quad-Core INTEL, 4 GB RAM). The FE models consisted of 197,000 brick elements (on average) and the average computation time was 5 min for the linear and 8 h for the nonlinear models. It should be mentioned that in voxel-based FE, the initial element size is dictated by the voxel size, which in our case was 0.5 × 0.5 × 1 mm3. In order to enhance the elemental aspect ratio and also reduce the computation time, the element size was changed to 1 × 1 × 1 mm3. The comparison between the analysis results of the two different mesh sizes showed that the difference was less than 0.2%. A simple homemade computer code was used to sort the elemental risk factor (RF) by computing the ratio of the strain energy density to the yield strain energy density for each element. The locus of the elements with the upmost RFs (critical elements) was considered as the failure initiation site. The development of the failure pattern was simulated by increasing the percentage of the screened critical elements [13,14] (Fig. 3). In the strength analysis procedure, a new scheme was developed and implemented for assortment of the elements based on their RF. The average of risk factors of all the elements (excluding those with an elastic modulus of 0.01 and those related to boundary conditions) was designated as RFav. On the other hand, the load modification factor (LMF) was defined as the ratio of the applied force in the FEA to the experimental fracture load, i.e., LMF = FFEA / Fexp. Initially, we used the holdout method, in which the specimens were randomly divided into the training and the test datasets. In the training dataset, a percentage of elements (with uppermost RFs) that would
Fig. 2. (a) Schematic of the mechanical testing fixture. (b) A specimen with steel sleeve. (c) The mechanical gripping fixture mounted inside the testing machine.
M. Mirzaei et al. / Bone 64 (2014) 108–114
111
Fig. 3. Left: The analyses results for specimens 4, 5, and 9 (the pictures a, b, and c show the development of the failure pattern predicted by linear-FEA), along with the experimental results. Right: Linear-FE predictions of local failure initiation for 4 specimens loaded at stance configuration along with the results of similar case studies reported in the literature, Bessho et al. [8], Keyak and Falkinstein [18], and Cody et al. [12].
have an average RF equal to LMF were selected for each model. These percentages were plotted against RFav and a linear regression was performed to obtain a relating expression for these two parameters. This expression was then used for the prediction of the fracture load in the test dataset through a reverse procedure, i.e., the RFav was calculated for each model and fed into the obtained linear relationship (diagram (a) in Fig. 4) to determine the percentage of the critical elements. The average risk factor of these critical elements (LMF), along with the applied FE load, was then used to predict the experimental failure load. The above procedure was repeated in two other ways for cross validation. In the first approach, the whole data set was used for training and the resulting percentage-RFav expression was used for the prediction of the failure load of every single specimen. The second approach was carried out on a leave-one-out basis, in which the percentage-RFav expressions were obtained for 10 sets of 9-specimen groups and were used for the prediction of the failure load for the excluded specimen from each group. Mechanical testing Mechanical tests were performed using a gripping fixture specifically designed for testing specimens in the stance configuration (Fig. 2). This fixture is capable of working within the range of − 30° b α b 30°
and − 30° b β b 30°, which covers almost all the physiologic and the non-physiologic axial loadings. The following procedure was used to grip and align the specimen in the appropriate direction [13]. 1. The distal diaphysis (length of 85 mm) was placed within a steel sleeve filled with PMMA. 2. A steel pin was inserted through the whole assembly. 3. The grooves on the outer surface of the sleeve were matched with the three sliding faces of the gripping fixture and thereby the slip rotation was prevented without damaging the sample. 4. The coronal planes of the specimens were aligned with the fixture reference. A mechanical testing machine (INSTRON Corporation, Canton, MA) was used to perform the experimental tests. In order to create a uniform load distribution, steel caps of different sizes were used based on the femoral head diameter. Initially, a preload of 100N was applied and the application of the main load was sufficiently delayed to ensure the stabilization of load and displacement readings. The continuation of loading until fracture was performed at a displacement rate of 1 mm/ min. The vertical deflection of the base plate of the fixture was measured and subtracted from the readings of the total displacement. The
112
M. Mirzaei et al. / Bone 64 (2014) 108–114
Fig. 4. The correlations between: (a) the percentage of the critical elements and the average RF of all elements for the training dataset, (b) linear FE predictions and experimental loads obtained for the test dataset, (c) the entire experimental loads and the linear FE loads predicted by the method of Nishiyama et al. [16], (d) the percentage of the critical elements and the average RF of all elements for the entire dataset, (e) experimental and linear FE loads obtained for the entire dataset, (f) experimental and linear FE loads (cross validation using the leave-one-out method), (g) experimental and nonlinear FE loads, and (h) nonlinear and linear FE loads.
femoral stiffness was defined as the slope of initial linear section of the obtained load–displacement data, and the maximum load during each test was considered as the failure strength. Results Failure patterns Fig. 3 shows the predicted and experimental failure patterns for three specimens (numbers 4, 5, and 9) under different loading orientations (see Table 1). It should be noted that in this figure only the screened critical elements are shown for clarity. For the specimen 4, the initiation of a trochanteric local failure was distinguished. The occurrence of a similar failure pattern at the same location (in the form of a distinct crack) was clearly visible on the experimental specimen (see Fig. 3). It should be mentioned that the fracture shown in the picture is in fact the result of damage initiation adjacent to the lower trochanter. Further development of this damage zone is in the form of crack growth towards the distal part. Nevertheless, the proposed linear method, which is based on the locus of the elements with upmost risk factors, can only predict the location of damage initiation and limited growth.
In fact, the prediction of the crack propagation path beyond this region certainly requires extra modeling efforts based on fracture mechanics approach. The initiation and limited-growth areas are shown by arrows in the picture. For the specimen 5, which experienced a totally different initiation and propagation pattern, a distinct subcapital damaged zone was observed. The same failure pattern was predicted by the linear analysis as shown in Fig. 3. Moreover, the intertrochanteric fracture of specimen 9 was quite comparable with the predicted failure pattern. Finally, the right column of Fig. 3 shows the linear FE results obtained for 4 specimens loaded at stance configuration of α = 20°, β = 0°, along with the results of similar case studies reported in the literature. It is clear that for this loading orientation the specimens of this study show subcapital failure patterns, which agree with the experimental and computational results reported by other researchers. Failure loads Fig. 4(a) shows a very good correlation (R2 = 0.86, p b 0.01) between the percentage of the critical elements and the average RF of all elements for the training dataset. The application of the obtained linear expression in the prediction of the experimental loads of the test dataset
M. Mirzaei et al. / Bone 64 (2014) 108–114
was very successful (R2 = 0.88, p b 0.01) (as depicted in Fig. 4(b)). Also note that the slope of the correlation between the predicted and experimental results (0.96) is very close to 1, which indicates a sound correlation. Fig. 4(c) shows that the application of the method proposed by Nishiyama et al. [16] to the entire set of experimental loads, resulted in a weaker correlation (R2 = 0.80, p b 0.01) with the slope of 0.59. The correlations shown in Figs. 4(d) and (e), which have been obtained for the entire dataset, can be considered more comprehensive and appropriate for the prediction of the femoral strength using the proposed method. Fig. 4(f) shows the results of cross validation using the leaveone-out method, which is another indication of the accuracy and robustness of the method. The correlation between the experimental loads and nonlinear FE predictions is shown in Fig. 4(g). The resulting correlation (R2 = 0.85, slope = 0.79, intercept = 20.5 N, for 10 specimens) can be compared with the results of a very comprehensive study reported by Koivumäki et al. (R2 = 0.86, slope = 0.92, intercept = 258 N, for 61 specimens) [9]. Also, it should be noted that our correlation was for nonlinear analyses of different loading orientations in the stance configuration, but their results were reported for elasto-plastic analyses of a sideways fall configuration. Finally, Fig. 4(h) presents a side-to-side nonlinear versus linear comparison.
Discussion Many investigations have been carried out on non-invasive assessment of the hip fracture risk using the QCT-based finite element analysis (FEA). The valuable advantages of this method are the ability of development of very accurate image-based 3D models of the femoral geometry, plus a pointwise description of BMD-based material properties. In general, the results have shown that nonlinear FEA can predict the femoral strength with very good accuracy. However, several researchers have faced certain difficulties on the practical aspects of this method. For instance, Koivumäki et al. carried out bi-linear elasto-plastic FEA to predict the fracture load of proximal femur in the sideways fall configuration and mentioned that the method was time consuming and not practical for clinical use [9]. Bessho et al. implemented bi-linear elastoplastic models for the prediction of the failure load in the stance configuration and reported average computational times of 10 and 30 h for 3 and 2 mm mesh sizes, respectively [7]. On the other hand, some previous attempts on the usage of the QCTbased linear analysis for strength prediction have not shown very favorable results (R2 = 0.72 and slope = 1.66 for stance, and a nonlinear correlation for fall configuration [20]). A more successful method was recently reported by Nishayama et al. who suggested that failure of the whole proximal femur can occur when 7% of the voxels exceed 0.9% strain. Accordingly, they calculated the failure load by performing fast linear analysis in the fall configuration (R2 = 0.81, slope = 0.68) [16]. Nevertheless, the purpose of the current study was to develop an efficient and reliable linear method for the prediction of the failure strength of the proximal femur, with the same accuracy of comparable nonlinear models. In practice, bone fractures occur beyond the range of linear loading and different parts of the nonlinear portion of the overall load–displacement diagram are representative of various failure mechanisms. However, it should be noted that our aim is to predict a single load (ultimate strength) on the nonlinear diagram without construction of the nonlinear structural behavior. We reason that the overall nonlinear behavior of the load–displacement diagram prior to the ultimate strength is the manifestation of the local failure of high-risk regions (not the failure of all material points). Thus, in our method, the regions with the uppermost strain energy densities are identified and assorted as high-risk or critical regions using the linear analysis data. Accordingly, the overall ultimate strength is considered as the load magnitude at which the strain energy of the critical regions equals the yield strain
113
energy. We believe that this idea is logical and realistic, so the method can give a safe estimate of the failure load. In practice, this method requires a linear analysis of a QCT-voxel based FE model at an arbitrary load level, which can be lower or higher than the actual load that is to be predicted. Next, the elements are sorted based on their risk factor (RF), defined as the ratio of the strain energy density to the yield strain energy for each element. The former is a direct outcome of the FE postprocessor and the latter is calculated for each element by additional BMD-based information about the material yield properties (see the expressions in the Finite element analysis Section). Accordingly, the average RF of all the elements (except for the elements related to the boundary conditions) is calculated and fed in to the proposed linear relationship (percentage = 46.33 × RFav − 0.68) to determine the percentage of the critical elements. The average risk factor of these critical elements is the load modification factor (LMF), which is the required correction factor for the applied FE load to predict the femoral strength (Fexp = FFEA / LMF). We reason that the strain energy is a product of the stress and strain tensors, so it can be a very good representative of both force and deformation intensities. We also argue that the number of the critical elements not only depends on the loading orientation, but is also related to the distribution of the weak points throughout the femur. For instance, it is quite reasonable that for the cases that the applied loading orientation and/or the particular distribution of material property can cause several favorable yielding or failure sites, the average RF of all the elements (RFav) can be relatively high. Therefore, the number of critical elements that should be screened for the failure load prediction must naturally be higher. On the other hand, for the cases that loading condition causes extreme stress–strain gradients and/or there might be relatively fewer weak material spots throughout the bone, the number of favorable failure sites (also the RFav) would be less and the critical elements should naturally be limited. It should be noted that the above argument is naturally reflected in the proposed linear expression. The very good agreement between the predicted and experimental failure loads (R2 = 0.86, slope = 0.96) shows that the suggested method is reliable for the prediction of the failure load. We also found a very good agreement between the predicted and actual failure patterns using the strain energy density measure. The ability of this measure in the prediction of the failure pattern can be attributed to the micro-mechanism of local failures in the bone tissue. The micromechanism of failure of porous trabecular tissue is mostly in the form of spicule (trabecula) buckling, and the failure mechanism of denser trabecular and nearly-homogenous cortical tissues is mostly of the form of local cracking. In either case (buckling or cracking), the strain energy density can be considered as a viable damage controlling parameter from a solid mechanics point of view [14]. At this point, it is beneficial to discuss some specific issues related to precision, efficiency, and applicability of the proposed method. We believe that increasing the accuracy and versatility of our image processing stage by incorporation of new techniques (like the proposed methods reported by Hangartner and Short [23]) can enhance the performance of the method. Varghese et al. [24] showed that the implementation of these delicate segmentation techniques can enhance the accuracy of QCT-based FE models and give precise predictions of surface strains (so that they are even comparable with the strain gage measurements). However, since the focus of our present study is on fast prediction of the global failure strength, the overall accuracy of our method (in its current form) might be less sensitive to such delicate enhancements of the segmentation process. Another aspect is the usage of the same properties for the compressive and tensile mechanical behavior of the bone tissue. Although this might be a source of error, it did not result in a significant deviation of FEA results from experiments. Moreover, the lack of a thoroughly validated failure theory for bone, and the findings of Keyak et al. [21,22], indicates that accounting for this difference will not necessarily increase the accuracy of finite element simulations.
114
M. Mirzaei et al. / Bone 64 (2014) 108–114
Although the focus of the current study was on determination of the ultimate failure strength, we believe that by implementation of more sophisticated assortments of critical elements, it might be possible to predict other characteristic loads, like the onset of the overall nonlinear behavior using linear analysis. Finally, we speculate that the proposed method can be applied to the whole femur, or bones other than the femur, because it includes a set of systematic steps that naturally leads to the construction of a reliable computational model of the bone and implements a series of logical steps of data processing which is independent of the type of bone specimens. Acknowledgments The authors wish to thank Dr. Mohsen Sadeghi for his continuous help during several stages of this work. The valuable help of Dr. Mahdavi and Dr. Khodadadi from the Iranian Tissue Bank (ITB, Tehran University of Medical Sciences) in providing Cadaveric samples is highly appreciated. The QCT scans were carried out using the facilities of the Noor Clinic with the help and support of Dr. Akhlaghpour. The valuable help of Mr. Kargar (Materials Laboratory of Tarbiat Modares University) with the mechanical tests is highly appreciated. This study was funded by Tarbiat Modares University. Conflicts of interest statement There is no conflict of interests regarding this manuscript. References [1] Soveid M, Serati AR, Masoompoor M. Incidence of hip fracture in Shiraz, Iran. Osteoporos Int 2005;16:1412–6. [2] Ballane G, Cauley JA, Arabi A, Fuleihan GE-H. Chapter 27 — geographic variability in hip and vertebral fractures. In: Marcus R, Feldman D, Dempster DW, Luckey M, Cauley JA, editors. Osteoporos, 4th ed., vol. 45. San Diego: Academic Press; 2013. p. 623–44. [3] Gullberg B, Johnell O, Kanis JA. World-wide projections for hip fracture. Osteoporos Int 1997;7:407–13. [4] Keyak J. Improved prediction of proximal femoral fracture load using nonlinear finite element models. Med Eng Phys 2001;23:165–73. [5] Keyak JH, Kaneko TS, Tehranzadeh J, Skinner HB. Predicting proximal femoral strength using structural engineering models. Clin Orthop Relat Res 2005;437:219–28. [6] Dragomir-Daescu D, Op Den Buijs J, McEligot S, Dai Y, Entwistle RC, Salas C, et al. Robust QCT/FEA models of proximal femur stiffness and fracture load during a sideways fall on the hip. Ann Biomed Eng 2011;39:742–55.
[7] Bessho M, Ohnishi I, Matsuyama J, Matsumoto T, Imai K, Nakamura K. Prediction of strength and strain of the proximal femur by a CT-based finite element method. J Biomech 2007;40:1745–53. [8] Bessho M, Ohnishi I, Matsumoto T, Ohashi S, Matsuyama J, Tobita K, et al. Prediction of proximal femur strength using a CT-based nonlinear finite element method: differences in predicted fracture load and site with changing load and boundary conditions. Bone 2009;45:226–31. [9] Koivumäki JEM, Thevenot J, Pulkkinen P, Kuhn V, Link TM, Eckstein F, et al. Ct-based finite element models can be used to estimate experimentally measured failure loads in the proximal femur. Bone 2012;50:824–9. [10] Schileo E, Taddei F, Cristofolini L, Viceconti M. Subject-specific finite element models implementing a maximum principal strain criterion are able to estimate failure risk and fracture location on human femurs tested in vitro. J Biomech 2008;41:356–67. [11] Keaveny TM, McClung MR, Wan X, Kopperdahl DL, Mitlak BH, Krohn K. Femoral strength in osteoporotic women treated with teriparatide or alendronate. Bone 2012;50:165–70. [12] Cody DD, Gross GJ, Hou FJ, Spencer HJ, Goldstein SA, Fyhrie DP. Femoral strength is better predicted by finite element models than QCT and DXA. J Biomech 1999;32: 1013–20. [13] Mirzaei M, Samiezadeh S, Khodadadi A, Ghazavi MR. Finite element prediction and experimental verification of the failure pattern of proximal femur using quantitative computed tomography images. Proceedings of the International Conference on Biomechanics and Biomedical Engineering, Copenhagen, Denmark; 2012. p. 111–7. [14] Mirzaei M, Zeinali A, Razmjoo A, Nazemi M. On prediction of the strength levels and failure patterns of human vertebrae using quantitative computed tomography (QCT)-based finite element method. J Biomech 2009;42:1584–91. [15] Keyak JH, Skinner HB, Fleming JA. Effect of force direction on femoral fracture load for two types of loading conditions. J Orthop Res 2001;19:539–44. [16] Nishiyama KK, Gilchrist S, Guy P, Cripton P, Boyd SK. Proximal femur bone strength estimated by a computationally fast finite element analysis in a sideways fall configuration. J Biomech 2013;46:1231–6. [17] Les CM, Keyak JH, Stover SM, Taylor KT, Kaneps AJ. Estimation of material properties in the equine metacarpus with use of quantitative computed tomography. J Orthop Res 1994;12:822–33. [18] Keyak JH, Falkinstein Y. Comparison of in situ and in vitro CT scan-based finite element model predictions of proximal femoral fracture load. Med Eng Phys 2003;25: 781–7. [19] Keller TS. Predicting the compressive mechanical behavior of bone. J Biomech 1994;27:1159–68. [20] Keyak JH, Rossi SA, Jones KA, Skinner HB. Prediction of femoral fracture load using automated finite element modeling. J Biomech 1998;31:125–33. [21] Les CM, Keyak JH, Stover SM, Taylor KT. Development and validation of a series of three-dimensional finite element models of the equine metacarpus. J Biomech 1997;30:737–42. [22] Keyak JH, Rossi SA. Prediction of femoral fracture load using finite element models: an examination of stress- and strain-based failure theories. J Biomech 2000;33: 209–14. [23] Hangartner TN, Short DF. Accurate quantification of width and density of bone structures by computed tomography. Med Phys 2007;10:3777–84. [24] Varghese B, Short D, Penmesta R, Goswami T, Hangartner T. Computed-tomographybased finite-element models of long bones can accurately capture strain response to bending and torsion. J Biomech 2011;44:1374–9.
Contents lists available at ScienceDirect
Bone journal homepage: www.elsevier.com/locate/bone
Original Full Length Article
Analysis of strength and failure pattern of human proximal femur using quantitative computed tomography (QCT)-based finite element method Majid Mirzaei ⁎, Maziyar Keshavarzian, Vahid Naeini Department of Mechanical Engineering, Tarbiat Modares University, P.O. Box: 14115-143, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 14 December 2013 Revised 10 March 2014 Accepted 4 April 2014 Available online 13 April 2014 Edited by Richard Eastell Keywords: Bone Fracture Imaging Noninvasive methods Strain energy density
a b s t r a c t This paper presents a novel method for fast and reliable prediction of the failure strength of human proximal femur, using the quantitative computed tomography (QCT)-based linear finite element analysis (FEA). Ten fresh frozen human femora (age: 34 ± 16) were QCT-scanned and the pertinent 3D voxel-based finite element models were constructed. A specially-designed holding frame was used to define and maintain a unique geometrical reference system for both FEA and in-vitro mechanical testing. The analyses and tests were carried out at 8 different loading orientations. A new scheme was developed for assortment of the element risk factor (defined as the ratio of the strain energy density to the yield strain energy for each element) and implemented for the prediction of the failure strength. The predicted and observed failure patterns were in correspondence, and the FEA predictions of the failure loads were in very good agreement with the experimental results (R2 = 0.86, slope = 0.96, p b 0.01). The average computational time was 5 min (on a regular desktop personal computer) for an average element number of 197,000. Noting that the run-time for a similar nonlinear model is about 8 h, it was concluded that the proposed linear scheme is overwhelmingly efficient in terms of computational costs. Thus, it can efficiently be used to predict the femoral failure strength with the same accuracy of similar nonlinear models. © 2014 Elsevier Inc. All rights reserved.
Introduction Hip fractures cause significant numbers of disability and mortality worldwide. The incidence of these fractures varies globally, e.g., the city of Shiraz in Iran is reported to have the highest incidence of hip fractures in Asia [1]. Previous studies estimated that the worldwide number of hip fractures was 1.3 million in 1990 and predicted an increase between 7.3 and 21.3 million by 2050 [1–3]. Thus, many researchers have tried to devise new methods for noninvasive assessment of the hip fracture risk, among which the QCT-based finite element analysis (FEA) has been recognized as a very reliable tool [4–12]. However, there are some limitations that should be resolved before this method can be clinically feasible. One major problem is the high computational time and cost for implementation of nonlinear FEA which has been shown to have the required accuracy [4,6–9,13]. A typical QCT-based nonlinear FEA of human proximal femur requires at least 8–10 h of runtime, even by implementing advanced computational techniques [4,7,9,13]. On the other hand, the linear FE analyses are very much faster, quite easier to apply, and can provide comparable predictions of the failure location and pattern for the proximal femur and vertebra [13,14]. However, they require further data processing schemes for predictions of the failure load that may give less accurate results in certain cases. ⁎ Corresponding author. Fax: +98 21 82884909. E-mail addresses: [email protected], [email protected] (M. Mirzaei).
http://dx.doi.org/10.1016/j.bone.2014.04.007 8756-3282/© 2014 Elsevier Inc. All rights reserved.
For instance, Keyak et al. proposed a linear FEA method for the prediction of the failure load in the sideways fall and stance configurations and argued that failure of a certain number of contiguous non-surface elements can represent the failure of the whole proximal femur [15]. On the other hand, Nishiyama et al. defined a constant yield strain level and suggested that the failure of a constant percentage of elements would represent the failure of proximal femur in the sideways fall configuration [16]. Nevertheless, we believe that the percentage or the number of critical elements should naturally be different for various specimens and loading orientations. Moreover, our previous studies on human vertebrae and proximal femur have shown favorable results in the prediction of the failure patterns using a strain energy density measure within the framework of linear QCT-based FEA [13,14]. Thus, this study was aimed at development and experimental verification of a new method for fast and accurate prediction of the failure strength of human proximal femur through specific assortments and usage of high strain energy elements of the linear QCT-based FEA results. Materials and methods Sample preparation Ten fresh-frozen human femora, from 9 cadavers (5 male, 4 female, average age: 34 ± 16), were used in this study (see Table 1) (specimens 9 and 10 were a contralateral pair from a single cadaver).
M. Mirzaei et al. / Bone 64 (2014) 108–114
plane (the femoral coronal plane) using three different views, i.e., the medial, the axial top, and the axial bottom views as follows (see Fig. 1(c)):
Table 1 Age, sex, side, and loading orientation of specimens. Number
1 2 3 4 5 6 7 8 9 10
Age
45 50 24 25 44 35 37 40 20 20
Sex
M M F F M F F M M M
109
Side Right/left
Orientation α
β
L R R L L R L R R L
10 10 10 15 15 0 20 −20 0 0
0 0 0 30 0 −28 0 0 −20 20
The initial treatment of the specimens included standard excision, freezing, and bagging procedures by the Iranian Tissue Bank (ITB). Moreover, the death causes announced by ITB were the three categories of coronary failure, cerebral death, and fatal trauma. Nevertheless, the plain radiographs of the samples were examined to ascertain the lack of metastatic diseases, pathologic defects, or insufficiency fractures. Definition and implementation of the femoral reference system The orientations of mechanical loadings on proximal femurs are usually adjusted with respect to the coronal plane, which is defined as a plane containing the femoral cervical and shaft axes. Nevertheless, certain ambiguities exist about recognition of these axes on femoral specimens, because the complex geometries of femoral parts lack any axisymmetric feature. Moreover, the lack of a unique definition of these axes for both FEA and mechanical testing can be troublesome. We used the reference system proposed by Mirzaei et al. [13], in which 3 distinct points are marked on the femoral surface to establish a reference
• The center point of the femoral head from the medial view (point 1). • The center point of the femoral head from the axial top view (point 2). • The center point of the femoral shaft from the axial bottom view (point 3). Despite the fact that the femoral head is not quite spherical and the femoral shaft is not quite cylindrical, the three points can be marked without ambiguity by circumscription of the related views in quadrilaterals as depicted in Fig. 1(c). In practice, a specially-designed holding frame was used to designate these points on the samples as depicted in Fig. 1(b). The frame has two side plates, which move simultaneously in opposite directions using special connecting rods so that, upon contact with the femoral head, the median plane of the head coincides with the frame plane. The frame also has a two-pin slider for alignment of the femoral shaft and three screws which are used to specify both the frame plane and the femoral coronal plane. The role of the forth screw is just to add extra fixity. This reference frame remains attached to the sample during the QCT process and is naturally transferred to the 3D and the FE models. On the other hand, the role of the Plexiglas disk attached to the lower end of the distal femur is to preserve and transfer the reference system to the mechanical testing. The disk has three screws which are used to bring it to the desired location and orientation with respect to both the femoral shaft and the holding frame. In its final position, the disk plane is parallel to the lower side of the holding frame, its axis of rotation passes through point 3, and the direction of its long screw lies in the defined coronal plane. Thus, when the femur is detached from the holding frame, the disk preserves the location and orientation of the defined coronal plane.
Fig. 1. (a) A proximal femur specimen. (b) The specimen inside the holding frame. (c) Three different views of a proximal femur and the three points identified on the femoral surface for definition of the coronal plane. (d) The imaging and calibration phantoms inside the scanner. (e) The CT image of a section of the specimen and calibration rods. (f) The spherical coordinate system used in this study. (g) Conventional coordinate system for loading of femur; α (the angle between the applied load and the sagittal plane) and β (the angle between the load and the coronal plane). (h) The 3D solid model of the specimen. (i) The FE model of the specimen.
110
M. Mirzaei et al. / Bone 64 (2014) 108–114
QCT scanning and image processing Each framed specimen was placed inside a Plexiglas container filled with water, designed as an imaging phantom for QCT studies (Fig. 1(d)). The QCT scans were carried out using a clinical scanner (SIEMENSSOMATOM 64, 140 kV, 80 mAs, 0.5 × 0.5 mm/pixel resolution, and 1 mm slice thickness). Using a calibration phantom (MINDWAYS Software, Inc., San Francisco, CA) (Figs. 1(d,e)), gray scale values were mapped to K2HPO4 equivalent density (ρKHP) using five tubes with reference densities and the Hounsfield Units (HUs) were calibrated. The segmentation of the bone hard tissue from the surroundings was performed for each slice, using a homemade computer code that implements the MATLAB's Image Processing Toolbox [14]. The raw QCT images (in DICOM format) were converted into binary format and a combination of user-defined threshold limit with an edgefollowing scheme was used for generation of the hard-tissue contours and elimination of the soft tissue.
Finite element analysis A complementary homemade MATLAB code was used to build the FE models by conversion of each voxel into an 8-noded brick element. In brief, this code treats each voxel as a basis for construction of a single brick element and prepares an input file (node and element file) which is compatible with the format of the ANSYS software (ANSYS. Inc., Canonsburg, PA) (see Fig. 1(i)). Calibrated CT scan data [17] (ρash = 1.22ρKHP + 0.0526) were used to calculate the ash density (ρash, g/cm3). The bone tissue was assumed as an inhomogeneous isotropic material with linearly elastic properties for linear analyses and linearly elastic-perfectly plastic behavior for nonlinear analyses. Young's modulus of elasticity (E, MPa) and yield strength (S, MPa) values were assigned to each element using the empirical relationships based on calibrated ash density for the trabecular and the cortical bone [18–20]; E = 33,900ρash2.20 for ρash ≤ 0.27 (trabecular bone), E = 10,200ρash2.01 for ρash ≥ 0.6 (cortical bone), E = 5,307ρash + 469 for 0.27 b ρash b 0.6 (transition), S = 137ρash1.88 for ρash b 0.317 (trabecular bone), S = 114ρash1.72 for ρash ≥ 0.317 (cortical bone). The three reference points were identified from CT scans and the coronal plane was recognized and used for the application of the load and boundary conditions. However, definition of a general coordinate system for application of mechanical loads on femur is another important issue in the FEA of femoral samples. Conventionally, two angles are defined and used for loading of the femur: α (the angle between the applied load and the sagittal plane) and β (the angle between the
load and the coronal plane) (see Fig. 1(g)). In this study, we used a spherical coordinate system on the femoral head as depicted in Fig. 1(f). The following relationship exists between the conventional α–β and the new ϕ–θ systems: tan θ ¼ tan β= tan α; tan ϕ ¼ tan α= cos θ:
ð1Þ
The boundary conditions in the finite element model were consistent with the experiments, i.e., the load was equally distributed among the nodes of femoral head (which were in contact with the 30-mm steel cap) and the nodes on the clamped cross section of the shaft were fully restrained (Fig. 2). All elements with the modulus below 5 MPa were assigned a low modulus of 0.01 MPa [20]. The elements loaded on the femoral head were identified and assigned an elastic modulus of 20 GPa and a yield strength of 200 MPa to prevent element distortion [5]. The Poisson's ratio was assumed to be 0.4 [20]. The models were loaded at 8 different orientations (for replication of physiologic and non-physiological loadings) to test the validity of models under various loading conditions (see Table 1). The linear and nonlinear FE models were analyzed using ANSYS software (ANSYS. Inc., Canonsburg, PA) on a desktop PC (2 × 2.7 GHz Quad-Core INTEL, 4 GB RAM). The FE models consisted of 197,000 brick elements (on average) and the average computation time was 5 min for the linear and 8 h for the nonlinear models. It should be mentioned that in voxel-based FE, the initial element size is dictated by the voxel size, which in our case was 0.5 × 0.5 × 1 mm3. In order to enhance the elemental aspect ratio and also reduce the computation time, the element size was changed to 1 × 1 × 1 mm3. The comparison between the analysis results of the two different mesh sizes showed that the difference was less than 0.2%. A simple homemade computer code was used to sort the elemental risk factor (RF) by computing the ratio of the strain energy density to the yield strain energy density for each element. The locus of the elements with the upmost RFs (critical elements) was considered as the failure initiation site. The development of the failure pattern was simulated by increasing the percentage of the screened critical elements [13,14] (Fig. 3). In the strength analysis procedure, a new scheme was developed and implemented for assortment of the elements based on their RF. The average of risk factors of all the elements (excluding those with an elastic modulus of 0.01 and those related to boundary conditions) was designated as RFav. On the other hand, the load modification factor (LMF) was defined as the ratio of the applied force in the FEA to the experimental fracture load, i.e., LMF = FFEA / Fexp. Initially, we used the holdout method, in which the specimens were randomly divided into the training and the test datasets. In the training dataset, a percentage of elements (with uppermost RFs) that would
Fig. 2. (a) Schematic of the mechanical testing fixture. (b) A specimen with steel sleeve. (c) The mechanical gripping fixture mounted inside the testing machine.
M. Mirzaei et al. / Bone 64 (2014) 108–114
111
Fig. 3. Left: The analyses results for specimens 4, 5, and 9 (the pictures a, b, and c show the development of the failure pattern predicted by linear-FEA), along with the experimental results. Right: Linear-FE predictions of local failure initiation for 4 specimens loaded at stance configuration along with the results of similar case studies reported in the literature, Bessho et al. [8], Keyak and Falkinstein [18], and Cody et al. [12].
have an average RF equal to LMF were selected for each model. These percentages were plotted against RFav and a linear regression was performed to obtain a relating expression for these two parameters. This expression was then used for the prediction of the fracture load in the test dataset through a reverse procedure, i.e., the RFav was calculated for each model and fed into the obtained linear relationship (diagram (a) in Fig. 4) to determine the percentage of the critical elements. The average risk factor of these critical elements (LMF), along with the applied FE load, was then used to predict the experimental failure load. The above procedure was repeated in two other ways for cross validation. In the first approach, the whole data set was used for training and the resulting percentage-RFav expression was used for the prediction of the failure load of every single specimen. The second approach was carried out on a leave-one-out basis, in which the percentage-RFav expressions were obtained for 10 sets of 9-specimen groups and were used for the prediction of the failure load for the excluded specimen from each group. Mechanical testing Mechanical tests were performed using a gripping fixture specifically designed for testing specimens in the stance configuration (Fig. 2). This fixture is capable of working within the range of − 30° b α b 30°
and − 30° b β b 30°, which covers almost all the physiologic and the non-physiologic axial loadings. The following procedure was used to grip and align the specimen in the appropriate direction [13]. 1. The distal diaphysis (length of 85 mm) was placed within a steel sleeve filled with PMMA. 2. A steel pin was inserted through the whole assembly. 3. The grooves on the outer surface of the sleeve were matched with the three sliding faces of the gripping fixture and thereby the slip rotation was prevented without damaging the sample. 4. The coronal planes of the specimens were aligned with the fixture reference. A mechanical testing machine (INSTRON Corporation, Canton, MA) was used to perform the experimental tests. In order to create a uniform load distribution, steel caps of different sizes were used based on the femoral head diameter. Initially, a preload of 100N was applied and the application of the main load was sufficiently delayed to ensure the stabilization of load and displacement readings. The continuation of loading until fracture was performed at a displacement rate of 1 mm/ min. The vertical deflection of the base plate of the fixture was measured and subtracted from the readings of the total displacement. The
112
M. Mirzaei et al. / Bone 64 (2014) 108–114
Fig. 4. The correlations between: (a) the percentage of the critical elements and the average RF of all elements for the training dataset, (b) linear FE predictions and experimental loads obtained for the test dataset, (c) the entire experimental loads and the linear FE loads predicted by the method of Nishiyama et al. [16], (d) the percentage of the critical elements and the average RF of all elements for the entire dataset, (e) experimental and linear FE loads obtained for the entire dataset, (f) experimental and linear FE loads (cross validation using the leave-one-out method), (g) experimental and nonlinear FE loads, and (h) nonlinear and linear FE loads.
femoral stiffness was defined as the slope of initial linear section of the obtained load–displacement data, and the maximum load during each test was considered as the failure strength. Results Failure patterns Fig. 3 shows the predicted and experimental failure patterns for three specimens (numbers 4, 5, and 9) under different loading orientations (see Table 1). It should be noted that in this figure only the screened critical elements are shown for clarity. For the specimen 4, the initiation of a trochanteric local failure was distinguished. The occurrence of a similar failure pattern at the same location (in the form of a distinct crack) was clearly visible on the experimental specimen (see Fig. 3). It should be mentioned that the fracture shown in the picture is in fact the result of damage initiation adjacent to the lower trochanter. Further development of this damage zone is in the form of crack growth towards the distal part. Nevertheless, the proposed linear method, which is based on the locus of the elements with upmost risk factors, can only predict the location of damage initiation and limited growth.
In fact, the prediction of the crack propagation path beyond this region certainly requires extra modeling efforts based on fracture mechanics approach. The initiation and limited-growth areas are shown by arrows in the picture. For the specimen 5, which experienced a totally different initiation and propagation pattern, a distinct subcapital damaged zone was observed. The same failure pattern was predicted by the linear analysis as shown in Fig. 3. Moreover, the intertrochanteric fracture of specimen 9 was quite comparable with the predicted failure pattern. Finally, the right column of Fig. 3 shows the linear FE results obtained for 4 specimens loaded at stance configuration of α = 20°, β = 0°, along with the results of similar case studies reported in the literature. It is clear that for this loading orientation the specimens of this study show subcapital failure patterns, which agree with the experimental and computational results reported by other researchers. Failure loads Fig. 4(a) shows a very good correlation (R2 = 0.86, p b 0.01) between the percentage of the critical elements and the average RF of all elements for the training dataset. The application of the obtained linear expression in the prediction of the experimental loads of the test dataset
M. Mirzaei et al. / Bone 64 (2014) 108–114
was very successful (R2 = 0.88, p b 0.01) (as depicted in Fig. 4(b)). Also note that the slope of the correlation between the predicted and experimental results (0.96) is very close to 1, which indicates a sound correlation. Fig. 4(c) shows that the application of the method proposed by Nishiyama et al. [16] to the entire set of experimental loads, resulted in a weaker correlation (R2 = 0.80, p b 0.01) with the slope of 0.59. The correlations shown in Figs. 4(d) and (e), which have been obtained for the entire dataset, can be considered more comprehensive and appropriate for the prediction of the femoral strength using the proposed method. Fig. 4(f) shows the results of cross validation using the leaveone-out method, which is another indication of the accuracy and robustness of the method. The correlation between the experimental loads and nonlinear FE predictions is shown in Fig. 4(g). The resulting correlation (R2 = 0.85, slope = 0.79, intercept = 20.5 N, for 10 specimens) can be compared with the results of a very comprehensive study reported by Koivumäki et al. (R2 = 0.86, slope = 0.92, intercept = 258 N, for 61 specimens) [9]. Also, it should be noted that our correlation was for nonlinear analyses of different loading orientations in the stance configuration, but their results were reported for elasto-plastic analyses of a sideways fall configuration. Finally, Fig. 4(h) presents a side-to-side nonlinear versus linear comparison.
Discussion Many investigations have been carried out on non-invasive assessment of the hip fracture risk using the QCT-based finite element analysis (FEA). The valuable advantages of this method are the ability of development of very accurate image-based 3D models of the femoral geometry, plus a pointwise description of BMD-based material properties. In general, the results have shown that nonlinear FEA can predict the femoral strength with very good accuracy. However, several researchers have faced certain difficulties on the practical aspects of this method. For instance, Koivumäki et al. carried out bi-linear elasto-plastic FEA to predict the fracture load of proximal femur in the sideways fall configuration and mentioned that the method was time consuming and not practical for clinical use [9]. Bessho et al. implemented bi-linear elastoplastic models for the prediction of the failure load in the stance configuration and reported average computational times of 10 and 30 h for 3 and 2 mm mesh sizes, respectively [7]. On the other hand, some previous attempts on the usage of the QCTbased linear analysis for strength prediction have not shown very favorable results (R2 = 0.72 and slope = 1.66 for stance, and a nonlinear correlation for fall configuration [20]). A more successful method was recently reported by Nishayama et al. who suggested that failure of the whole proximal femur can occur when 7% of the voxels exceed 0.9% strain. Accordingly, they calculated the failure load by performing fast linear analysis in the fall configuration (R2 = 0.81, slope = 0.68) [16]. Nevertheless, the purpose of the current study was to develop an efficient and reliable linear method for the prediction of the failure strength of the proximal femur, with the same accuracy of comparable nonlinear models. In practice, bone fractures occur beyond the range of linear loading and different parts of the nonlinear portion of the overall load–displacement diagram are representative of various failure mechanisms. However, it should be noted that our aim is to predict a single load (ultimate strength) on the nonlinear diagram without construction of the nonlinear structural behavior. We reason that the overall nonlinear behavior of the load–displacement diagram prior to the ultimate strength is the manifestation of the local failure of high-risk regions (not the failure of all material points). Thus, in our method, the regions with the uppermost strain energy densities are identified and assorted as high-risk or critical regions using the linear analysis data. Accordingly, the overall ultimate strength is considered as the load magnitude at which the strain energy of the critical regions equals the yield strain
113
energy. We believe that this idea is logical and realistic, so the method can give a safe estimate of the failure load. In practice, this method requires a linear analysis of a QCT-voxel based FE model at an arbitrary load level, which can be lower or higher than the actual load that is to be predicted. Next, the elements are sorted based on their risk factor (RF), defined as the ratio of the strain energy density to the yield strain energy for each element. The former is a direct outcome of the FE postprocessor and the latter is calculated for each element by additional BMD-based information about the material yield properties (see the expressions in the Finite element analysis Section). Accordingly, the average RF of all the elements (except for the elements related to the boundary conditions) is calculated and fed in to the proposed linear relationship (percentage = 46.33 × RFav − 0.68) to determine the percentage of the critical elements. The average risk factor of these critical elements is the load modification factor (LMF), which is the required correction factor for the applied FE load to predict the femoral strength (Fexp = FFEA / LMF). We reason that the strain energy is a product of the stress and strain tensors, so it can be a very good representative of both force and deformation intensities. We also argue that the number of the critical elements not only depends on the loading orientation, but is also related to the distribution of the weak points throughout the femur. For instance, it is quite reasonable that for the cases that the applied loading orientation and/or the particular distribution of material property can cause several favorable yielding or failure sites, the average RF of all the elements (RFav) can be relatively high. Therefore, the number of critical elements that should be screened for the failure load prediction must naturally be higher. On the other hand, for the cases that loading condition causes extreme stress–strain gradients and/or there might be relatively fewer weak material spots throughout the bone, the number of favorable failure sites (also the RFav) would be less and the critical elements should naturally be limited. It should be noted that the above argument is naturally reflected in the proposed linear expression. The very good agreement between the predicted and experimental failure loads (R2 = 0.86, slope = 0.96) shows that the suggested method is reliable for the prediction of the failure load. We also found a very good agreement between the predicted and actual failure patterns using the strain energy density measure. The ability of this measure in the prediction of the failure pattern can be attributed to the micro-mechanism of local failures in the bone tissue. The micromechanism of failure of porous trabecular tissue is mostly in the form of spicule (trabecula) buckling, and the failure mechanism of denser trabecular and nearly-homogenous cortical tissues is mostly of the form of local cracking. In either case (buckling or cracking), the strain energy density can be considered as a viable damage controlling parameter from a solid mechanics point of view [14]. At this point, it is beneficial to discuss some specific issues related to precision, efficiency, and applicability of the proposed method. We believe that increasing the accuracy and versatility of our image processing stage by incorporation of new techniques (like the proposed methods reported by Hangartner and Short [23]) can enhance the performance of the method. Varghese et al. [24] showed that the implementation of these delicate segmentation techniques can enhance the accuracy of QCT-based FE models and give precise predictions of surface strains (so that they are even comparable with the strain gage measurements). However, since the focus of our present study is on fast prediction of the global failure strength, the overall accuracy of our method (in its current form) might be less sensitive to such delicate enhancements of the segmentation process. Another aspect is the usage of the same properties for the compressive and tensile mechanical behavior of the bone tissue. Although this might be a source of error, it did not result in a significant deviation of FEA results from experiments. Moreover, the lack of a thoroughly validated failure theory for bone, and the findings of Keyak et al. [21,22], indicates that accounting for this difference will not necessarily increase the accuracy of finite element simulations.
114
M. Mirzaei et al. / Bone 64 (2014) 108–114
Although the focus of the current study was on determination of the ultimate failure strength, we believe that by implementation of more sophisticated assortments of critical elements, it might be possible to predict other characteristic loads, like the onset of the overall nonlinear behavior using linear analysis. Finally, we speculate that the proposed method can be applied to the whole femur, or bones other than the femur, because it includes a set of systematic steps that naturally leads to the construction of a reliable computational model of the bone and implements a series of logical steps of data processing which is independent of the type of bone specimens. Acknowledgments The authors wish to thank Dr. Mohsen Sadeghi for his continuous help during several stages of this work. The valuable help of Dr. Mahdavi and Dr. Khodadadi from the Iranian Tissue Bank (ITB, Tehran University of Medical Sciences) in providing Cadaveric samples is highly appreciated. The QCT scans were carried out using the facilities of the Noor Clinic with the help and support of Dr. Akhlaghpour. The valuable help of Mr. Kargar (Materials Laboratory of Tarbiat Modares University) with the mechanical tests is highly appreciated. This study was funded by Tarbiat Modares University. Conflicts of interest statement There is no conflict of interests regarding this manuscript. References [1] Soveid M, Serati AR, Masoompoor M. Incidence of hip fracture in Shiraz, Iran. Osteoporos Int 2005;16:1412–6. [2] Ballane G, Cauley JA, Arabi A, Fuleihan GE-H. Chapter 27 — geographic variability in hip and vertebral fractures. In: Marcus R, Feldman D, Dempster DW, Luckey M, Cauley JA, editors. Osteoporos, 4th ed., vol. 45. San Diego: Academic Press; 2013. p. 623–44. [3] Gullberg B, Johnell O, Kanis JA. World-wide projections for hip fracture. Osteoporos Int 1997;7:407–13. [4] Keyak J. Improved prediction of proximal femoral fracture load using nonlinear finite element models. Med Eng Phys 2001;23:165–73. [5] Keyak JH, Kaneko TS, Tehranzadeh J, Skinner HB. Predicting proximal femoral strength using structural engineering models. Clin Orthop Relat Res 2005;437:219–28. [6] Dragomir-Daescu D, Op Den Buijs J, McEligot S, Dai Y, Entwistle RC, Salas C, et al. Robust QCT/FEA models of proximal femur stiffness and fracture load during a sideways fall on the hip. Ann Biomed Eng 2011;39:742–55.
[7] Bessho M, Ohnishi I, Matsuyama J, Matsumoto T, Imai K, Nakamura K. Prediction of strength and strain of the proximal femur by a CT-based finite element method. J Biomech 2007;40:1745–53. [8] Bessho M, Ohnishi I, Matsumoto T, Ohashi S, Matsuyama J, Tobita K, et al. Prediction of proximal femur strength using a CT-based nonlinear finite element method: differences in predicted fracture load and site with changing load and boundary conditions. Bone 2009;45:226–31. [9] Koivumäki JEM, Thevenot J, Pulkkinen P, Kuhn V, Link TM, Eckstein F, et al. Ct-based finite element models can be used to estimate experimentally measured failure loads in the proximal femur. Bone 2012;50:824–9. [10] Schileo E, Taddei F, Cristofolini L, Viceconti M. Subject-specific finite element models implementing a maximum principal strain criterion are able to estimate failure risk and fracture location on human femurs tested in vitro. J Biomech 2008;41:356–67. [11] Keaveny TM, McClung MR, Wan X, Kopperdahl DL, Mitlak BH, Krohn K. Femoral strength in osteoporotic women treated with teriparatide or alendronate. Bone 2012;50:165–70. [12] Cody DD, Gross GJ, Hou FJ, Spencer HJ, Goldstein SA, Fyhrie DP. Femoral strength is better predicted by finite element models than QCT and DXA. J Biomech 1999;32: 1013–20. [13] Mirzaei M, Samiezadeh S, Khodadadi A, Ghazavi MR. Finite element prediction and experimental verification of the failure pattern of proximal femur using quantitative computed tomography images. Proceedings of the International Conference on Biomechanics and Biomedical Engineering, Copenhagen, Denmark; 2012. p. 111–7. [14] Mirzaei M, Zeinali A, Razmjoo A, Nazemi M. On prediction of the strength levels and failure patterns of human vertebrae using quantitative computed tomography (QCT)-based finite element method. J Biomech 2009;42:1584–91. [15] Keyak JH, Skinner HB, Fleming JA. Effect of force direction on femoral fracture load for two types of loading conditions. J Orthop Res 2001;19:539–44. [16] Nishiyama KK, Gilchrist S, Guy P, Cripton P, Boyd SK. Proximal femur bone strength estimated by a computationally fast finite element analysis in a sideways fall configuration. J Biomech 2013;46:1231–6. [17] Les CM, Keyak JH, Stover SM, Taylor KT, Kaneps AJ. Estimation of material properties in the equine metacarpus with use of quantitative computed tomography. J Orthop Res 1994;12:822–33. [18] Keyak JH, Falkinstein Y. Comparison of in situ and in vitro CT scan-based finite element model predictions of proximal femoral fracture load. Med Eng Phys 2003;25: 781–7. [19] Keller TS. Predicting the compressive mechanical behavior of bone. J Biomech 1994;27:1159–68. [20] Keyak JH, Rossi SA, Jones KA, Skinner HB. Prediction of femoral fracture load using automated finite element modeling. J Biomech 1998;31:125–33. [21] Les CM, Keyak JH, Stover SM, Taylor KT. Development and validation of a series of three-dimensional finite element models of the equine metacarpus. J Biomech 1997;30:737–42. [22] Keyak JH, Rossi SA. Prediction of femoral fracture load using finite element models: an examination of stress- and strain-based failure theories. J Biomech 2000;33: 209–14. [23] Hangartner TN, Short DF. Accurate quantification of width and density of bone structures by computed tomography. Med Phys 2007;10:3777–84. [24] Varghese B, Short D, Penmesta R, Goswami T, Hangartner T. Computed-tomographybased finite-element models of long bones can accurately capture strain response to bending and torsion. J Biomech 2011;44:1374–9.