- Email: [email protected]
Incorporating in vivo fall assessments in the simulation of femoral fractures with finite element models
Medical Engineering and Physics 37 (2015) 593–598
Contents lists available at ScienceDirect
Medical Engineering and Physics journal homepage: www.elsevier.com/locate/medengphy
Technical note
Incorporating in vivo fall assessments in the simulation of femoral fractures with finite element models A.M. van der Zijden a,b,∗, D. Janssen a, N. Verdonschot a,c, B.E. Groen b, B. Nienhuis b, V. Weerdesteyn b,d, E. Tanck a a
Orthopaedic Research Laboratory, Radboud University Medical Center, Nijmegen, The Netherlands Research Department, Sint Maartenskliniek, Nijmegen, The Netherlands Laboratory for Biomechanical Engineering, University of Twente, Enschede, The Netherlands d Department of Rehabilitation, Radboud University Medical Center, Nijmegen, The Netherlands b c
a r t i c l e
i n f o
Article history: Received 2 December 2013 Revised 28 February 2015 Accepted 21 March 2015
Keywords: Sideways falls Finite element model Femoral fracture load Loading configuration Osteoporosis
a b s t r a c t Femoral fractures are a major health issue. Most experimental and finite element (FE) fracture studies use polymethylmethacrylate cups on the greater trochanter (GT) to simulate fall impact loads. However, in vivo fall studies showed that the femur is loaded distally from the GT. Our objective was to incorporate in vivo fall data in FE models to determine the effects of loading position and direction, and size of simulated impact site on the fracture load and fracture type for a healthy and an osteoporotic femur. Twelve sets of loading position and angles were applied through ‘near point loads’ on the models. Additional simulations were performed with ‘cup loads’ on the GT, similar to the literature. The results showed no significant difference between fracture loads from simulations with near point loads distally from the GT and those with cup loads on the GT. However, simulated fracture types differed, as near point loads distally from the GT generally resulted in various neck fractures, whilst cup load simulations predicted superior neck and trochanteric fractures only. This study showed that incorporating in vivo fall assessments in FE models by loading the models distally from the GT results in prediction of realistic fracture loads and fracture types. © 2015 IPEM. Published by Elsevier Ltd. All rights reserved.
1. Introduction Hip fractures are a major health issue in the elderly and osteoporotic patients worldwide [1], with high morbidity and mortality rates [2,3]. Assessing the hip fracture risk is challenging, as many risk factors are involved [3,4]. The most important risk factors for hip fractures are sideways falling and diminished bone quality [3–6]. From a biomechanical perspective, the hip fracture risk can be defined as the ratio of the impact load and the bone strength [7]. Many studies examined the fracture behavior of femora by performing mechanical tests with cadaveric femora or computed tomography (CT) based finite element (FE) analyses [8–10]. In such studies, the point of application and direction of the impact load (loading configuration) are important aspects for assessing the fracture load. In 1957, Backman [11] was the first to compare the outcomes of in vivo fracture experiments to clinical fractures, providing insights ∗ Corresponding author at: 547 ORL, PO Box 9101, 6500 HB Nijmegen, The Netherlands. Tel.: +31 024 361 7461. E-mail address: [email protected] (A.M. van der Zijden). URL: http://www.biomechanics.nl, http://www.umcn.nl (A.M. van der Zijden)
http://dx.doi.org/10.1016/j.medengphy.2015.03.006 1350-4533/© 2015 IPEM. Published by Elsevier Ltd. All rights reserved.
into the loading configuration on the greater trochanter (GT) in sideways falls. Lotz and Hayes [12] converted Backman’s results to an experimental loading configuration, loading the femur at the femoral head, with the GT embedded in a polymethylmethacrylate (PMMA) cup. The femur was oriented with a 10° angle between the femoral shaft and horizontal (frontal plane) and the femoral neck 15° internally rotated (transversal plane) [12–14]. In FE studies, the trochanteric region is usually fully restrained (mimicking the PMMA cup) to prevent local crushing [15]. This ‘default’ loading configuration has been used for simulating sideways falls in a number of studies [15–19]. In vivo fall experiments showed that impact direction, point of application, and impact load distribution, in sideways falls differ between subjects and fall strategies [20]. In addition, differences in local bone quality between impact sites and subjects should also be considered, for example healthy versus osteoporotic bones [21]. Several experimental and FE studies examined the effect of loading angle (in frontal and transversal planes) on the fracture load and pattern, showing that applying a load more posteriorly from the GT results in a decrease in fracture load of up to 24% [9,10,16,22]. However, little is known about the effects of distributing the load over a relatively large
594
A.M. van der Zijden et al. / Medical Engineering and Physics 37 (2015) 593–598
(PMMA cup load) versus a small simulated impact site (point load) and of fully constraining the trochanteric region. The pre-formed shape of the cup and its fixed direction of displacement restrict the bone from (rotational) movements, affecting the internal stress–strain distribution and simulated fracture load. Therefore, applying the ‘default’ PMMA cup load may introduce artifacts for fracture risk assessment. Furthermore, recent in vivo sideways fall data [23,24] showed that it is likely that the femur is loaded more distally from the GT than commonly assumed. FE fracture simulations with cup loads on the GT have been well validated in relation to experimental loading conditions. However, past FE studies have not attempted to incorporate impact sites placed more distally from the GT. Therefore, it remains questionable whether these fracture experiments represent real-life falls in terms of applied loading configuration. Incorporating in vivo fall assessments in fracture simulations would be more realistic and can be considered as the next step to improve the prediction of realistic fracture loads and types in real-life falls. The objective of the present study was to utilize data of in vivo fall assessments in FE models, to determine the effect of varying loading configurations (direction, point of application and size of impact site) on the fracture load and fracture type of a healthy and an osteoporotic femoral bone. We hypothesized that FE models with loading configurations derived from in vivo fall data will predict realistic fracture loads and types, converging toward more realistic FE simulation of real-life falls. 2. Methods 2.1. FE model Quantitative CT scans (ACQsim, Philips, Eindhoven, The Netherlands) were derived from two fresh-frozen human cadaveric proximal femora: one healthy femur (male, 81 years, t-score 0.0) and one osteoporotic femur (female, 81 years, t-score 2.8, [25]). 3D surface meshes were derived from the QCT data (Mimics 11.0, Materialise NV, Leuven, Belgium). Femoral geometric properties of the healthy and the osteoporotic bone model, such as femoral head offset (40.0 mm and 44.3 mm) and inclination angle (127° and 127°), were within normal ranges (47.0 ± 7.2 mm and 122.9 ± 7.6°, respectively) [26]. The surface meshes were converted into solid meshes (edge length of 2 mm) (Patran 2005r2, MSC Software Corporation, Santa Ana, USA). The healthy bone model consisted of 119,441 cortical and 155,101 trabecular tetrahedral elements (54,901 nodes) and the osteoporotic bone model consisted of 81,221 cortical and 135,820 trabecular elements (44,552 nodes). Non-linear isotropic bone mechanical properties and post failure material behavior of the elements were calculated with conversion algorithms of Keyak [27] using the Von Mises yield criterion. Material properties of each element were derived from the QCT slices, with a mean cortical ash density of 0.657 ± 0.32 g/cm3 for the healthy and 0.704 ± 0.31 g/cm3 for the osteoporotic bone model. Four consecutive phases of bone material behavior were defined: an elastic, an ideal plastic, a strain softening and, after failure, an indefinite ideal plastic phase [27,28]. These specific bone models were also part of a study of Derikx et al. [28] in which FE simulations were compared to outcomes of experiments with cadaveric bones loaded to fracture in a single-limb stance loading configuration. The fracture loads predicted by the FE simulations strongly correlated with the experimental fracture loads (R2 = 0.90). 2.2. Fall loading configurations Fall loading configurations were derived from previous in vivo fall experiments [24], in which 12 subjects performed ten standardized sideways falls from kneeling height, representing a natural fall arrest strategy. We used the observed loading variations in these falls as well
Fig. 1. The FE models were oriented based on the loading configurations and virtual identification of anatomical landmarks: the GT (grey dot), the hip joint center (white dot), femoral neck axis and femoral shaft axis (dashed lines). (A) The loading angle in the frontal plane was set to a default value (α = 10°). Four points of application of the impact load were applied (y = 0, 21, 49 and 77 mm). (B) The loading angle in the transversal plane (β = 2°, 15° and 28°).
as those used in the literature to define 12 different loading configurations (three sets of loading angles and four points of application) for the FE models as described below (Fig. 1A and B). The loading angle in the frontal plane (α ) as measured in the fall experiments [24] agreed well with the loading angle of 10° (angle between femoral shaft and horizontal) as used by Lotz and Hayes [12]. As the observed range in this angle was small [24], the α -loading angle was set to 10° for all simulations (Fig. 1A). For the transversal plane, the loading angles (β ) from the in vivo fall data were expected to be unreliable in terms of their absolute values, due to the difficulty of assessing internal and external rotations of the hip joint. Therefore, an internal rotation of 15° was chosen based on literature values (Fig. 1A) [11,12]. The range of the β -loading angle, was derived from the 13° standard deviation as measured in vivo (β = 2°;15°;28°) (Fig. 1A), which was comparable to literature values (range of 15°) [9,17,22]. The main difference between the fall loading configuration from in vivo fall experiments and the literature configuration involves the point of application of the impact load. The position of the impact load relative to the GT position (y) was varied within the range of the in vivo fall data (49 ± 28 mm) [24], leading to three points of application: y = 21, 49 and 77 mm distally from the GT. For comparison
A.M. van der Zijden et al. / Medical Engineering and Physics 37 (2015) 593–598
595
Fig. 2. Boundary conditions applied to the FE models. The distal part of the femur was connected to two attachment points in the knee joint by springs. Proximally, frictionless contact was defined between the femoral head and the acetabulum of a rigid and fully restrained hemi pelvis. Selected displacement nodes were displaced in z-direction with 0.2 mm per increment, simulating a (near) point load.
with the literature, an additional point of application was defined: y = 0 mm, resembling sideways falls directly onto the GT (Fig. 1B). Using Matlab (Matlab r2009b, The MathWorks Inc., Natick, USA), the FE models were oriented based on virtual identification of anatomical landmarks (hip and knee joint centers and GT). The femur’s coordinate system was oriented with the z-axis pointing in the direction of the impact load (Fig. 1A and B). Cross-sectional slices of the FE model were made parallel to the transversal plane (xz-plane), at the y-positions of the points of application. For these slices, the node closest to the horizontal was defined as the ‘load application’ node for the particular fall configuration (Fig. 1A). 2.3. Boundary conditions To simulate anatomical constraints at the femoral head, a rigid hemi-pelvis was included in the FE models. All nodes of the pelvis were fully restrained, thereby preventing displacements of the femoral head. Frictionless contact was defined between the acetabulum and the femoral head, allowing rotation of the femoral segment around the hip joint center. The distal part of the femur was represented by springs (E = 200 GPa), which were connected to two fixed attachment points in the knee joint (Fig. 2). Limited freedom of rotation around the knee joint in the frontal plane was allowed, similar to the experimental set-up [13,18]. The fracture simulations were performed with MSC.Marc 2007r1 (MSC Software Corporation, Santa Ana, USA). The fall impact load was applied to the femur by displacing surface nodes at the various
impact sites in the direction of the force vector (z-axis) by 0.2 mm per increment. For this purpose, the load application node was displaced in the z-direction, together with surface nodes of the elements surrounding the selected node (typically 7 nodes), referred to as ‘near point load’. The material properties of the elements surrounding the displacement nodes were set to ideal elastic to prevent plasticity artifacts (Fig. 3A). For comparison with the literature, we also included simulations with the boundary conditions as used in experimental studies [12,18,22] and FE studies from the literature [10,15– 17,19]. To mimic the PMMA cups, generally used in femoral fracture experiments, the number of nodes selected for displacement was increased, (Fig. 3B). These ‘cup load’ simulations (y = 0) were performed for all three loading angles in the transversal plane (β = 2°;15°;28°). Hence, a total of 15 simulations were performed for each FE model. 2.4. Fracture load Structural fracture of our FE models was assumed to occur when the maximum total reaction force was reached, i.e. the first peak in the force–displacement (FD) curve of the simulation, similar to previous work [28,29]. In the current study, FD-curves did not always show an evident peak as in the fall loading configuration failed bone elements tended to accumulate similarly to the crumple zone of a car. Therefore, the instant of fracture was defined by assessing the volume of failed cortical elements, in line with previous studies [16,17,19]. The threshold volume was based on the sensitivity study with simulations
Fig. 3. Material properties applied to the elements of the FE models. The surface elements of the femoral head (contact area with acetabulum) were set to ideal elastic. (A) For the simulations with ‘near point’ loads, a small number of nodes (7) were displaced, with a small number of surrounding surface elements set to ideal elastic. (B) For the simulations with ‘cup’ loads on the GT, a larger number of displacement nodes (330) were used and more elements were set to ideal elastic, mimicking the use of a PMMA cup.
596
A.M. van der Zijden et al. / Medical Engineering and Physics 37 (2015) 593–598
Fig. 4. Calculated fracture loads for the healthy (A) and the osteoporotic (B) bone model. Loading configurations derived from the literature and those derived from in vivo data are separated by the vertical dashed lines. (∗) Significant differences between point loads on GT (y = 0) and both cup loads on GT (y = 0) and point loads distally from GT (y = 21, 49 and 77 mm) (p < 0.05) (loading configurations), (#) Significant differences between β loading angle (p < 0.05).
that did show an evident peak (18 out of 30 simulations in total). Fracture loads were estimated by calculating the total reaction force in the increments exceeding the threshold volumes of 100, 250, 500, 750 and 1000 mm3 . These fracture loads were compared with peak forces in the FD-curves, showing a sufficient estimation of fracture load (98 ± 8%) and correlation (R2 > 0.95) at a threshold volume of 500 mm3 . 2.5. Statistical analysis Statistical tests were conducted for the healthy and osteoporotic bone models separately. We first compared the fracture loads between the simulations based on in vivo fall data (y = 21; 49; 77 mm), the point loads on the GT (y = 0) and the cup loads on the GT (y = 0) using one-way analysis of variance (ANOVA) with post-hoc Tukey’s test (α = 0.05). For the simulations based on in vivo fall data, we also assessed the effects of points of application (y = 21; 49; 77 mm) and loading angles (β = 2°; 15°; 28°) using the same statistical methods. Finally, differences between the two bone models were also assessed with one-way ANOVA. 2.6. Fracture initiation To assess the predicted fracture type, we determined the location of the first failing cortical surface element in the bone model. Locations of fracture initiation were categorized as fracture of impact site (local crushing), trochanteric region and femoral neck [16–18,30].
3. Results Fig. 4 shows the fracture loads for all simulations from this study. For the loading configurations derived from the in vivo fall data (y = 21; 49; 77 mm), the mean fracture load was 5124 ± 1003 N for the healthy and 2373 ± 317 N for the osteoporotic bone model. For the simulations with cup loads on the GT (y = 0), mean fracture loads were 4633 ± 634 N (healthy) and 2557 ± 139 N (osteoporotic), whereas for the point load simulations (y = 0) these were 1697 ± 161 N (healthy) and 1145 ± 234 N (osteoporotic). Differences were significant for both the healthy (F(2,14) = 18.0; p < 0.001) and the osteoporotic (F(2,14) = 25.0; p < 0.001) bone model, with the simulations with point loads on the GT (y = 0) yielding lower fracture loads than the other simulations. Increasing the loading angle in the transversal plane (β ) decreased the fracture loads by 34% (F(2,8) = 17.6; p = 0.003) and 19% (F(2,8) = 13.6; p = 0.006), for the healthy and osteoporotic bone model respectively. Post-hoc Tukey’s tests showed that for the healthy bone model, fracture loads were significantly lower in the β = 28° simulations compared to the simulations with β = 2° and 15°. For the osteoporotic bone model, fracture loads from the β = 2° simulations were significantly higher than the fracture loads in the β = 15° and 28° simulations. Applying the load more distally from the GT had no significant effect on the fracture load for both the healthy (F(2,8) = 0.3; p = 0.764) and osteoporotic bone model (F(2,8) = 0.5; p = 0.642). Predicted fracture types are shown in Table 1 and examples of the fracture types are shown in Fig. 5. When loaded with a near point
Table 1 Fracture types for the healthy and osteoporotic bone models: location of fracture initiation (failing of first cortical surface element). Bone model
β loading angle
Position of load application (mm) y = 0 (cup load)
y = 0 (point load)
y = 21 (point load)
y = 49 (point load)
y = 77 (point load)
Healthy
2° 15° 28°
Superior neck Trochanteric (post.-sup.) Trochanteric (post.-sup.)
Impact site Impact site Impact site
Impact site Ant.-sup. neck Ant.-sup. neck
Inferior neck Anterior neck Ant.-sup. neck
Inferior neck Ant.-sup. neck Ant.-sup. neck
Osteoporotic
2° 15° 28°
Superior neck Superior neck Trochanteric (post.-sup.)
Impact site Impact site Impact site
Impact site Superior neck Impact site
Superior neck Ant.-sup. neck Inferior neck
Ant.-inf. neck Ant.-sup. neck Impact site
A.M. van der Zijden et al. / Medical Engineering and Physics 37 (2015) 593–598
597
Fig. 5. Predicted fracture types: (A) impact site (local crushing); (B) neck; and (C) trochanteric fractures. Phases of bone material behavior (plastic deformity) are indicated in grey values: from white (elastic phase) to black (failure phase).
load on the GT (y = 0), local crushing occurred at impact site, i.e. the GT region. When a cup load was applied on the GT, fracture occurred either in the superior neck region in case the β -loading angle was small, or in the trochanteric region when the β -loading angles were larger. In most simulations derived from the in vivo fall data (y = 21; 49; 77 mm), fractures occurred in the neck region and were often initiated at the anterior–superior side of the neck (Table 1). In simulations with the osteoporotic bone model, local crushing at impact site occurred more often than in simulations with the healthy bone model (Table 1). Fracture loads were 50% lower (F(1,29) = 24.5; p < 0.001) in the osteoporotic bone model compared to the healthy bone model (Fig. 4). 4. Discussion This study provides new insights in the importance of loading configurations (direction, point of application and size of impact site) by incorporating results from in vivo fall assessments in FE simulations of femoral fractures. Our findings corroborate our hypothesis that incorporating in vivo fall assessments in FE models, results in the prediction of realistic fracture loads and types. No significant differences were found between fracture loads from these simulations (near point load distally from GT) and those with configurations derived from the literature (cup load on GT). As for fracture types, superior neck and trochanteric fractures occurred in the cup load simulations and various neck fractures (superior, inferior and anterior) were found in simulations based on in vivo fall data. Loading the bone models with a near point load distally from the GT results in realistic fracture loads and fracture types, similar to clinical and experimental fractures [17,18,30,31]. However, a point load directly on the GT did not, as only local crushing was observed for these simulations. Increasing the number of displacement nodes (cup load on GT) prevented local crushing at the impact site. However, because of the pre-formed shape of the cup, applying a cup load more distally will restrict the bone model from (small) rotational movements during the simulation. Therefore, simulating the use of a PMMA cup is likely to affect the internal stress–strain distribution and simulated fracture. These unintended restraints can be avoided by selecting a smaller number of displacement nodes at the simulated impact sites, i.e. near point loads distally from the GT. Allowing the femoral bone model to settle between the pelvic acetabulum and the
impact load, together with the loading configurations being derived from in vivo fall data [24], improved the accuracy of the fracture simulations from an anatomical and kinematic perspective. As nearly no cortical elements failed in the distal load areas (<5), we believe that applying a distributed load instead, would not have changed the outcomes of these simulations. Experimental validation of these FE models should be performed to confirm the assumptions made in this study. To prevent bone crushing at impact sites, future fracture simulations can be further improved by implementing a rigid floor surface impacting a femur covered by a soft tissue layer [32]. Furthermore, within the simulations with point loads distally from the GT, varying the position of the point of application had no effect on the fracture load or fracture type. Increasing the loading angle in the transversal plane (β ) decreased the fracture loads by 34% and 19%, for the healthy and osteoporotic bone model respectively, which is consistent with previous reports in the literature [9,22]. The fracture simulations were validated only indirectly under stance phase loading conditions, which resulted in adequate prediction by the FE models (n = 10, R2 = 0.90) [28]. Keyak et al. [10] applied the same method for fracture simulations with stance and fall loading configurations, both resulting in good correlations with mechanical tests. For the purpose of the present study, our methods were deemed sufficiently accurate, as simulated fractures were in agreement with literature findings [17,18] and clinical practice [30,31]. The definition of the instant of fracture is a matter of controversy in the literature. Keyak defined fracture as the instant when Von Mises stresses exceed initial element strength for 15 contiguous non-surface elements [10,15]. Bessho and Koivumäki [16,17,19] defined fracture of femur models when at least one cortical (shell) surface element had failed. In our study, we assumed structural fracture when the maximum total reaction force was reached, similar to the definition used in fracture experiments [10,18,19,28]. We used a volumetric measure of failed cortical elements (500 mm3 ) to estimate this increment, as FD-curves did not always show an evident peak. A sensitivity study on the threshold volume showed that differences between fracture loads derived with various threshold volumes (500; 750; 1000 mm3 ) were negligible. Hence, the conclusions of our study do not depend heavily on the specific definition of this volumetric measure. Only one healthy and one osteoporotic bone model were used for the simulations. Differences in geometrical properties may affect the internal load distribution of the femoral bone and thereby the fracture loads and types. As geometric properties of our bone models were
598
A.M. van der Zijden et al. / Medical Engineering and Physics 37 (2015) 593–598
within ranges of the normal elderly population [26], we assumed that these two models were good general representatives of healthy and osteoporotic bones. The subjects of the in vivo experiment (trained athletes) differed from the target group. However, for safety reasons, elderly and osteoporotic people were not included in these fall experiments. We applied loading configurations based on the averages and wide ranges found in this dataset and the literature. Therefore, we expected that the loading configurations were suitable for a general population, enabling to assess trends and relative effects. Other limitations were related to the mechanical behavior of elements, as we did not include asymmetric yielding (for example the Drucker Pragner criterion [19,29]), strain rate dependent behavior [33], or anisotropic behavior [34]. Furthermore, simulation of a realistic fall impact load should involve applying a dynamic impact load instead of a quasi-static load. Further limitations concern the absence of fat and muscle tissue and the oversimplification of the pelvis model. More realistic mechanical behavior of elements, soft tissue and dynamic impact loads must be implemented in future FE models and validated against biomechanical fracture and musculoskeletal modeling studies. In summary, this study showed the potential of incorporating the loading configurations derived from in vivo fall assessments in FE simulations, applying fall impact loads distally to the GT. The fracture loads were comparable between loading configurations derived from in vivo fall data and those used in the literature, whilst from a kinematic and anatomic perspective, the loading configurations were more realistic. Simulated fracture types varied between loading configurations and differed between and within the healthy and osteoporotic bone model, reflecting its dependency on the simulated loading configuration. In conclusion, our results reveal insights in the importance of loading configurations, providing a next step in FE simulations of real-life fall-related hip fractures. Funding Nothing to declare: no external funding was received for the study. Ethical approval The bone models were part of a previous study [28], in which the models were derived from CT scans of two cadaveric femora, obtained from the Department of Anatomy with institutional approval. Fall loading configurations were derived from in vivo fall data of a previous study [24], in which each subject signed an informed consent form prior to participation. The experiment protocol was approved by the Ethical Board of the region Arnhem–Nijmegen. Conflict of interest statement The authors declare that they have no conflict of interest. References [1] Kanis JA, Oden A, McCloskey E, Johansson H, Wahl DA, Cooper C. A systematic review of hip fracture incidence and probability of fracture worldwide. Osteoporosis Int 2012;23:S308. [2] Kannus P, Parkkari J, Koskinen S, Niemi S, Palvanen M, Jarvinen M, et al. Fall-induced injuries and deaths among older adults. JAMA-J Am Med Assoc 1999;281(20):1895–9. [3] Cummings SR, Melton LJ. Epidemiology and outcomes of osteoporotic fractures. Lancet 2002;359(9319):1761–7. [4] DargentMolina P, Favier F, Grandjean H, Baudoin C, Schott AM, Hausherr E, et al. Fall-related factors and risk of hip fracture: the EPIDOS prospective study. Lancet 1996;348(9024):416. [5] Jarvinen TLN, Sievanen H, Khan KM, Heinonen A, Kannus P. Shifting the focus in fracture prevention from osteoporosis to falls. Brit Med J 2008;336(7636):124–6.
[6] Cummings-Vaughn LA, Gammack JK. Falls, osteoporosis, and hip fractures. Med Clin N Am 2011;95(3):495–506. [7] Hayes WC, Myers ER, Robinovitch SN, VandenKroonenberg A, Courtney AC, McMahon TA. Etiology and prevention of age-related hip fractures. Bone 1996;18(1):S77–86. [8] Greenspan SL, Myers ER, Maitland LA, Resnick NM, Hayes WC. Fall severity and bone-mineral density as risk-factors for hip fracture in ambulatory elderly. JAMA-J Am Med Assoc 1994;271(2):128–33. [9] Ford CM, Keaveny TM, Hayes WC. The effect of impact direction on the structural capacity of the proximal femur during falls. J Bone Miner Res 1996;11(3):377–83. [10] 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(4):539–44. [11] Backman S. The proximal end of the femur: investigations with special reference to the etiology of femoral neck fractures; anatomical studies; roentgen projections; theoretical stress calculations; experimental production of fractures. Acta Radiol. Suppl 1957(146):1–166. [12] Lotz JC, Hayes WC. The use of quantitative computed-tomography to estimate risk of fracture of the hip from falls. J Bone Joint Surg Am 1990;72A(5):689–700. [13] Courtney AC, Wachtel EF, Myers ER, Hayes WC. Effects of loading rate on strength of the proximal femur. Calcified Tissue Int 1994;55(1):53–8. [14] Keyak JH, Rossi SA, Jones KA, Skinner HB. Prediction of femoral fracture load using automated finite element modeling. J Biomech 1998;31(2):125–33. [15] Keyak JH, Rossi SA, Jones KA, Les CM, Skinner HB. Prediction of fracture location in the proximal femur using finite element models. Med Eng Phys 2001;23(9):657–64. [16] Bessho M, Ohnishi I, Okazaki H, Sato W, Kominami H, Matsunaga S, et al. Prediction of the strength and fracture location of the femoral neck by CT-based finiteelement method: a preliminary study on patients with hip fracture. J Orthop Sci 2004;9(6):545–50. [17] 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(2):226–31. [18] de Bakker PM, Manske SL, Ebacher V, Oxland TR, Cripton PA, Guy P. During sideways falls proximal femur fractures initiate in the superolateral cortex: evidence from high-speed video of simulated fractures. J Biomech 2009;42(12):1917–25. [19] Koivumaki JE, 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(4):824–9. [20] Robinovitch SN, Inkster L, Maurer J, Warnick B. Strategies for avoiding hip impact during sideways falls. J Bone Miner Res 2003;18(7):1267–73. [21] Tsouknidas A, Anagnostidis K, Maliaris G, Michailidis N. Fracture risk in the femoral hip region: a finite element analysis supported experimental approach. J Biomech 2012;45(11):1959–64. [22] Pinilla TP, Boardman KC, Bouxsein ML, Myers ER, Hayes WC. Impact direction from a fall influences the failure load of the proximal femur as much as age-related bone loss. Calcif Tissue Int 1996;58(4):231–5. [23] Choi WJ, Hoffer JA, Robinovitch SN. Effect of hip protectors, falling angle and body mass index on pressure distribution over the hip during simulated falls. Clin Biomech 2010;25(1):63–9. [24] van der Zijden AM, Groen BE, Tanck E, Nienhuis B, Verdonschot N, Weerdesteyn V. Can martial arts techniques reduce fall severity? An in vivo study of femoral loading configurations in sideways falls. J Biomech 2012;45(9):1650–5. [25] World Health Organization. Prevention and management of osteoporosis. World Health Organ Tech Rep Ser 2003;921:1–164. [26] Rubin PJ, Leyvraz PF, Aubaniac JM, Argenson JN, Esteve P, Deroguin B. The morphology of the proximal femur – a 3-dimensional radiographic analysis. J Bone Joint Surg Br 1992;74(1):28–32. [27] Keyak JH, Kaneko TS, Tehranzadeh J, Skinner HB. Predicting proximal femoral strength using structural engineering models. Clin Orthop Rel Res 2005;(437):219–28. [28] Derikx LC, van Aken JB, Janssen D, Snyers A, van der Linden YM, et al. The assessment of the risk of fracture in femora with metastatic lesions: comparing case-specific finite element analyses with predictions by clinical experts. J Bone Joint Surg Br 2012;94(8):1135–42. [29] Derikx LC, Vis R, Meinders T, Verdonschot N, Tanck E. Implementation of asymmetric yielding in case-specific finite element models improves the prediction of femoral fractures. Comput Meth Biomech 2011;14(2):183–93. [30] Zuckerman JD, Rosenberg AD. Hip fracture – reply. N Engl J Med 1996;335(26):1995–6. [31] Marks R, Allegrante JP, MacKenzie CR, Lane JM. Hip fractures among the elderly: causes, consequences and control. Ageing Res Rev 2003;2(1):57–93. [32] Laing AC, Robinovitch SN. The force attenuation provided by hip protectors depends on impact velocity, pelvic size, and soft tissue stiffness. J Biomech Eng-Trans ASME 2008;130(6):061005. [33] Hansen U, Zioupos P, Simpson R, Currey JD, Hynd D. The effect of strain rate on the mechanical properties of human cortical bone. J Biomech Eng-Trans ASME 2008;130(1):011011. [34] Keyak JH, Lee IY, Nath DS, Skinner HB. Postfailure compressive behavior of tibial trabecular bone – three anatomic directions. J Biomed Mater Res 1996;31(3):373–8.
Contents lists available at ScienceDirect
Medical Engineering and Physics journal homepage: www.elsevier.com/locate/medengphy
Technical note
Incorporating in vivo fall assessments in the simulation of femoral fractures with finite element models A.M. van der Zijden a,b,∗, D. Janssen a, N. Verdonschot a,c, B.E. Groen b, B. Nienhuis b, V. Weerdesteyn b,d, E. Tanck a a
Orthopaedic Research Laboratory, Radboud University Medical Center, Nijmegen, The Netherlands Research Department, Sint Maartenskliniek, Nijmegen, The Netherlands Laboratory for Biomechanical Engineering, University of Twente, Enschede, The Netherlands d Department of Rehabilitation, Radboud University Medical Center, Nijmegen, The Netherlands b c
a r t i c l e
i n f o
Article history: Received 2 December 2013 Revised 28 February 2015 Accepted 21 March 2015
Keywords: Sideways falls Finite element model Femoral fracture load Loading configuration Osteoporosis
a b s t r a c t Femoral fractures are a major health issue. Most experimental and finite element (FE) fracture studies use polymethylmethacrylate cups on the greater trochanter (GT) to simulate fall impact loads. However, in vivo fall studies showed that the femur is loaded distally from the GT. Our objective was to incorporate in vivo fall data in FE models to determine the effects of loading position and direction, and size of simulated impact site on the fracture load and fracture type for a healthy and an osteoporotic femur. Twelve sets of loading position and angles were applied through ‘near point loads’ on the models. Additional simulations were performed with ‘cup loads’ on the GT, similar to the literature. The results showed no significant difference between fracture loads from simulations with near point loads distally from the GT and those with cup loads on the GT. However, simulated fracture types differed, as near point loads distally from the GT generally resulted in various neck fractures, whilst cup load simulations predicted superior neck and trochanteric fractures only. This study showed that incorporating in vivo fall assessments in FE models by loading the models distally from the GT results in prediction of realistic fracture loads and fracture types. © 2015 IPEM. Published by Elsevier Ltd. All rights reserved.
1. Introduction Hip fractures are a major health issue in the elderly and osteoporotic patients worldwide [1], with high morbidity and mortality rates [2,3]. Assessing the hip fracture risk is challenging, as many risk factors are involved [3,4]. The most important risk factors for hip fractures are sideways falling and diminished bone quality [3–6]. From a biomechanical perspective, the hip fracture risk can be defined as the ratio of the impact load and the bone strength [7]. Many studies examined the fracture behavior of femora by performing mechanical tests with cadaveric femora or computed tomography (CT) based finite element (FE) analyses [8–10]. In such studies, the point of application and direction of the impact load (loading configuration) are important aspects for assessing the fracture load. In 1957, Backman [11] was the first to compare the outcomes of in vivo fracture experiments to clinical fractures, providing insights ∗ Corresponding author at: 547 ORL, PO Box 9101, 6500 HB Nijmegen, The Netherlands. Tel.: +31 024 361 7461. E-mail address: [email protected] (A.M. van der Zijden). URL: http://www.biomechanics.nl, http://www.umcn.nl (A.M. van der Zijden)
http://dx.doi.org/10.1016/j.medengphy.2015.03.006 1350-4533/© 2015 IPEM. Published by Elsevier Ltd. All rights reserved.
into the loading configuration on the greater trochanter (GT) in sideways falls. Lotz and Hayes [12] converted Backman’s results to an experimental loading configuration, loading the femur at the femoral head, with the GT embedded in a polymethylmethacrylate (PMMA) cup. The femur was oriented with a 10° angle between the femoral shaft and horizontal (frontal plane) and the femoral neck 15° internally rotated (transversal plane) [12–14]. In FE studies, the trochanteric region is usually fully restrained (mimicking the PMMA cup) to prevent local crushing [15]. This ‘default’ loading configuration has been used for simulating sideways falls in a number of studies [15–19]. In vivo fall experiments showed that impact direction, point of application, and impact load distribution, in sideways falls differ between subjects and fall strategies [20]. In addition, differences in local bone quality between impact sites and subjects should also be considered, for example healthy versus osteoporotic bones [21]. Several experimental and FE studies examined the effect of loading angle (in frontal and transversal planes) on the fracture load and pattern, showing that applying a load more posteriorly from the GT results in a decrease in fracture load of up to 24% [9,10,16,22]. However, little is known about the effects of distributing the load over a relatively large
594
A.M. van der Zijden et al. / Medical Engineering and Physics 37 (2015) 593–598
(PMMA cup load) versus a small simulated impact site (point load) and of fully constraining the trochanteric region. The pre-formed shape of the cup and its fixed direction of displacement restrict the bone from (rotational) movements, affecting the internal stress–strain distribution and simulated fracture load. Therefore, applying the ‘default’ PMMA cup load may introduce artifacts for fracture risk assessment. Furthermore, recent in vivo sideways fall data [23,24] showed that it is likely that the femur is loaded more distally from the GT than commonly assumed. FE fracture simulations with cup loads on the GT have been well validated in relation to experimental loading conditions. However, past FE studies have not attempted to incorporate impact sites placed more distally from the GT. Therefore, it remains questionable whether these fracture experiments represent real-life falls in terms of applied loading configuration. Incorporating in vivo fall assessments in fracture simulations would be more realistic and can be considered as the next step to improve the prediction of realistic fracture loads and types in real-life falls. The objective of the present study was to utilize data of in vivo fall assessments in FE models, to determine the effect of varying loading configurations (direction, point of application and size of impact site) on the fracture load and fracture type of a healthy and an osteoporotic femoral bone. We hypothesized that FE models with loading configurations derived from in vivo fall data will predict realistic fracture loads and types, converging toward more realistic FE simulation of real-life falls. 2. Methods 2.1. FE model Quantitative CT scans (ACQsim, Philips, Eindhoven, The Netherlands) were derived from two fresh-frozen human cadaveric proximal femora: one healthy femur (male, 81 years, t-score 0.0) and one osteoporotic femur (female, 81 years, t-score 2.8, [25]). 3D surface meshes were derived from the QCT data (Mimics 11.0, Materialise NV, Leuven, Belgium). Femoral geometric properties of the healthy and the osteoporotic bone model, such as femoral head offset (40.0 mm and 44.3 mm) and inclination angle (127° and 127°), were within normal ranges (47.0 ± 7.2 mm and 122.9 ± 7.6°, respectively) [26]. The surface meshes were converted into solid meshes (edge length of 2 mm) (Patran 2005r2, MSC Software Corporation, Santa Ana, USA). The healthy bone model consisted of 119,441 cortical and 155,101 trabecular tetrahedral elements (54,901 nodes) and the osteoporotic bone model consisted of 81,221 cortical and 135,820 trabecular elements (44,552 nodes). Non-linear isotropic bone mechanical properties and post failure material behavior of the elements were calculated with conversion algorithms of Keyak [27] using the Von Mises yield criterion. Material properties of each element were derived from the QCT slices, with a mean cortical ash density of 0.657 ± 0.32 g/cm3 for the healthy and 0.704 ± 0.31 g/cm3 for the osteoporotic bone model. Four consecutive phases of bone material behavior were defined: an elastic, an ideal plastic, a strain softening and, after failure, an indefinite ideal plastic phase [27,28]. These specific bone models were also part of a study of Derikx et al. [28] in which FE simulations were compared to outcomes of experiments with cadaveric bones loaded to fracture in a single-limb stance loading configuration. The fracture loads predicted by the FE simulations strongly correlated with the experimental fracture loads (R2 = 0.90). 2.2. Fall loading configurations Fall loading configurations were derived from previous in vivo fall experiments [24], in which 12 subjects performed ten standardized sideways falls from kneeling height, representing a natural fall arrest strategy. We used the observed loading variations in these falls as well
Fig. 1. The FE models were oriented based on the loading configurations and virtual identification of anatomical landmarks: the GT (grey dot), the hip joint center (white dot), femoral neck axis and femoral shaft axis (dashed lines). (A) The loading angle in the frontal plane was set to a default value (α = 10°). Four points of application of the impact load were applied (y = 0, 21, 49 and 77 mm). (B) The loading angle in the transversal plane (β = 2°, 15° and 28°).
as those used in the literature to define 12 different loading configurations (three sets of loading angles and four points of application) for the FE models as described below (Fig. 1A and B). The loading angle in the frontal plane (α ) as measured in the fall experiments [24] agreed well with the loading angle of 10° (angle between femoral shaft and horizontal) as used by Lotz and Hayes [12]. As the observed range in this angle was small [24], the α -loading angle was set to 10° for all simulations (Fig. 1A). For the transversal plane, the loading angles (β ) from the in vivo fall data were expected to be unreliable in terms of their absolute values, due to the difficulty of assessing internal and external rotations of the hip joint. Therefore, an internal rotation of 15° was chosen based on literature values (Fig. 1A) [11,12]. The range of the β -loading angle, was derived from the 13° standard deviation as measured in vivo (β = 2°;15°;28°) (Fig. 1A), which was comparable to literature values (range of 15°) [9,17,22]. The main difference between the fall loading configuration from in vivo fall experiments and the literature configuration involves the point of application of the impact load. The position of the impact load relative to the GT position (y) was varied within the range of the in vivo fall data (49 ± 28 mm) [24], leading to three points of application: y = 21, 49 and 77 mm distally from the GT. For comparison
A.M. van der Zijden et al. / Medical Engineering and Physics 37 (2015) 593–598
595
Fig. 2. Boundary conditions applied to the FE models. The distal part of the femur was connected to two attachment points in the knee joint by springs. Proximally, frictionless contact was defined between the femoral head and the acetabulum of a rigid and fully restrained hemi pelvis. Selected displacement nodes were displaced in z-direction with 0.2 mm per increment, simulating a (near) point load.
with the literature, an additional point of application was defined: y = 0 mm, resembling sideways falls directly onto the GT (Fig. 1B). Using Matlab (Matlab r2009b, The MathWorks Inc., Natick, USA), the FE models were oriented based on virtual identification of anatomical landmarks (hip and knee joint centers and GT). The femur’s coordinate system was oriented with the z-axis pointing in the direction of the impact load (Fig. 1A and B). Cross-sectional slices of the FE model were made parallel to the transversal plane (xz-plane), at the y-positions of the points of application. For these slices, the node closest to the horizontal was defined as the ‘load application’ node for the particular fall configuration (Fig. 1A). 2.3. Boundary conditions To simulate anatomical constraints at the femoral head, a rigid hemi-pelvis was included in the FE models. All nodes of the pelvis were fully restrained, thereby preventing displacements of the femoral head. Frictionless contact was defined between the acetabulum and the femoral head, allowing rotation of the femoral segment around the hip joint center. The distal part of the femur was represented by springs (E = 200 GPa), which were connected to two fixed attachment points in the knee joint (Fig. 2). Limited freedom of rotation around the knee joint in the frontal plane was allowed, similar to the experimental set-up [13,18]. The fracture simulations were performed with MSC.Marc 2007r1 (MSC Software Corporation, Santa Ana, USA). The fall impact load was applied to the femur by displacing surface nodes at the various
impact sites in the direction of the force vector (z-axis) by 0.2 mm per increment. For this purpose, the load application node was displaced in the z-direction, together with surface nodes of the elements surrounding the selected node (typically 7 nodes), referred to as ‘near point load’. The material properties of the elements surrounding the displacement nodes were set to ideal elastic to prevent plasticity artifacts (Fig. 3A). For comparison with the literature, we also included simulations with the boundary conditions as used in experimental studies [12,18,22] and FE studies from the literature [10,15– 17,19]. To mimic the PMMA cups, generally used in femoral fracture experiments, the number of nodes selected for displacement was increased, (Fig. 3B). These ‘cup load’ simulations (y = 0) were performed for all three loading angles in the transversal plane (β = 2°;15°;28°). Hence, a total of 15 simulations were performed for each FE model. 2.4. Fracture load Structural fracture of our FE models was assumed to occur when the maximum total reaction force was reached, i.e. the first peak in the force–displacement (FD) curve of the simulation, similar to previous work [28,29]. In the current study, FD-curves did not always show an evident peak as in the fall loading configuration failed bone elements tended to accumulate similarly to the crumple zone of a car. Therefore, the instant of fracture was defined by assessing the volume of failed cortical elements, in line with previous studies [16,17,19]. The threshold volume was based on the sensitivity study with simulations
Fig. 3. Material properties applied to the elements of the FE models. The surface elements of the femoral head (contact area with acetabulum) were set to ideal elastic. (A) For the simulations with ‘near point’ loads, a small number of nodes (7) were displaced, with a small number of surrounding surface elements set to ideal elastic. (B) For the simulations with ‘cup’ loads on the GT, a larger number of displacement nodes (330) were used and more elements were set to ideal elastic, mimicking the use of a PMMA cup.
596
A.M. van der Zijden et al. / Medical Engineering and Physics 37 (2015) 593–598
Fig. 4. Calculated fracture loads for the healthy (A) and the osteoporotic (B) bone model. Loading configurations derived from the literature and those derived from in vivo data are separated by the vertical dashed lines. (∗) Significant differences between point loads on GT (y = 0) and both cup loads on GT (y = 0) and point loads distally from GT (y = 21, 49 and 77 mm) (p < 0.05) (loading configurations), (#) Significant differences between β loading angle (p < 0.05).
that did show an evident peak (18 out of 30 simulations in total). Fracture loads were estimated by calculating the total reaction force in the increments exceeding the threshold volumes of 100, 250, 500, 750 and 1000 mm3 . These fracture loads were compared with peak forces in the FD-curves, showing a sufficient estimation of fracture load (98 ± 8%) and correlation (R2 > 0.95) at a threshold volume of 500 mm3 . 2.5. Statistical analysis Statistical tests were conducted for the healthy and osteoporotic bone models separately. We first compared the fracture loads between the simulations based on in vivo fall data (y = 21; 49; 77 mm), the point loads on the GT (y = 0) and the cup loads on the GT (y = 0) using one-way analysis of variance (ANOVA) with post-hoc Tukey’s test (α = 0.05). For the simulations based on in vivo fall data, we also assessed the effects of points of application (y = 21; 49; 77 mm) and loading angles (β = 2°; 15°; 28°) using the same statistical methods. Finally, differences between the two bone models were also assessed with one-way ANOVA. 2.6. Fracture initiation To assess the predicted fracture type, we determined the location of the first failing cortical surface element in the bone model. Locations of fracture initiation were categorized as fracture of impact site (local crushing), trochanteric region and femoral neck [16–18,30].
3. Results Fig. 4 shows the fracture loads for all simulations from this study. For the loading configurations derived from the in vivo fall data (y = 21; 49; 77 mm), the mean fracture load was 5124 ± 1003 N for the healthy and 2373 ± 317 N for the osteoporotic bone model. For the simulations with cup loads on the GT (y = 0), mean fracture loads were 4633 ± 634 N (healthy) and 2557 ± 139 N (osteoporotic), whereas for the point load simulations (y = 0) these were 1697 ± 161 N (healthy) and 1145 ± 234 N (osteoporotic). Differences were significant for both the healthy (F(2,14) = 18.0; p < 0.001) and the osteoporotic (F(2,14) = 25.0; p < 0.001) bone model, with the simulations with point loads on the GT (y = 0) yielding lower fracture loads than the other simulations. Increasing the loading angle in the transversal plane (β ) decreased the fracture loads by 34% (F(2,8) = 17.6; p = 0.003) and 19% (F(2,8) = 13.6; p = 0.006), for the healthy and osteoporotic bone model respectively. Post-hoc Tukey’s tests showed that for the healthy bone model, fracture loads were significantly lower in the β = 28° simulations compared to the simulations with β = 2° and 15°. For the osteoporotic bone model, fracture loads from the β = 2° simulations were significantly higher than the fracture loads in the β = 15° and 28° simulations. Applying the load more distally from the GT had no significant effect on the fracture load for both the healthy (F(2,8) = 0.3; p = 0.764) and osteoporotic bone model (F(2,8) = 0.5; p = 0.642). Predicted fracture types are shown in Table 1 and examples of the fracture types are shown in Fig. 5. When loaded with a near point
Table 1 Fracture types for the healthy and osteoporotic bone models: location of fracture initiation (failing of first cortical surface element). Bone model
β loading angle
Position of load application (mm) y = 0 (cup load)
y = 0 (point load)
y = 21 (point load)
y = 49 (point load)
y = 77 (point load)
Healthy
2° 15° 28°
Superior neck Trochanteric (post.-sup.) Trochanteric (post.-sup.)
Impact site Impact site Impact site
Impact site Ant.-sup. neck Ant.-sup. neck
Inferior neck Anterior neck Ant.-sup. neck
Inferior neck Ant.-sup. neck Ant.-sup. neck
Osteoporotic
2° 15° 28°
Superior neck Superior neck Trochanteric (post.-sup.)
Impact site Impact site Impact site
Impact site Superior neck Impact site
Superior neck Ant.-sup. neck Inferior neck
Ant.-inf. neck Ant.-sup. neck Impact site
A.M. van der Zijden et al. / Medical Engineering and Physics 37 (2015) 593–598
597
Fig. 5. Predicted fracture types: (A) impact site (local crushing); (B) neck; and (C) trochanteric fractures. Phases of bone material behavior (plastic deformity) are indicated in grey values: from white (elastic phase) to black (failure phase).
load on the GT (y = 0), local crushing occurred at impact site, i.e. the GT region. When a cup load was applied on the GT, fracture occurred either in the superior neck region in case the β -loading angle was small, or in the trochanteric region when the β -loading angles were larger. In most simulations derived from the in vivo fall data (y = 21; 49; 77 mm), fractures occurred in the neck region and were often initiated at the anterior–superior side of the neck (Table 1). In simulations with the osteoporotic bone model, local crushing at impact site occurred more often than in simulations with the healthy bone model (Table 1). Fracture loads were 50% lower (F(1,29) = 24.5; p < 0.001) in the osteoporotic bone model compared to the healthy bone model (Fig. 4). 4. Discussion This study provides new insights in the importance of loading configurations (direction, point of application and size of impact site) by incorporating results from in vivo fall assessments in FE simulations of femoral fractures. Our findings corroborate our hypothesis that incorporating in vivo fall assessments in FE models, results in the prediction of realistic fracture loads and types. No significant differences were found between fracture loads from these simulations (near point load distally from GT) and those with configurations derived from the literature (cup load on GT). As for fracture types, superior neck and trochanteric fractures occurred in the cup load simulations and various neck fractures (superior, inferior and anterior) were found in simulations based on in vivo fall data. Loading the bone models with a near point load distally from the GT results in realistic fracture loads and fracture types, similar to clinical and experimental fractures [17,18,30,31]. However, a point load directly on the GT did not, as only local crushing was observed for these simulations. Increasing the number of displacement nodes (cup load on GT) prevented local crushing at the impact site. However, because of the pre-formed shape of the cup, applying a cup load more distally will restrict the bone model from (small) rotational movements during the simulation. Therefore, simulating the use of a PMMA cup is likely to affect the internal stress–strain distribution and simulated fracture. These unintended restraints can be avoided by selecting a smaller number of displacement nodes at the simulated impact sites, i.e. near point loads distally from the GT. Allowing the femoral bone model to settle between the pelvic acetabulum and the
impact load, together with the loading configurations being derived from in vivo fall data [24], improved the accuracy of the fracture simulations from an anatomical and kinematic perspective. As nearly no cortical elements failed in the distal load areas (<5), we believe that applying a distributed load instead, would not have changed the outcomes of these simulations. Experimental validation of these FE models should be performed to confirm the assumptions made in this study. To prevent bone crushing at impact sites, future fracture simulations can be further improved by implementing a rigid floor surface impacting a femur covered by a soft tissue layer [32]. Furthermore, within the simulations with point loads distally from the GT, varying the position of the point of application had no effect on the fracture load or fracture type. Increasing the loading angle in the transversal plane (β ) decreased the fracture loads by 34% and 19%, for the healthy and osteoporotic bone model respectively, which is consistent with previous reports in the literature [9,22]. The fracture simulations were validated only indirectly under stance phase loading conditions, which resulted in adequate prediction by the FE models (n = 10, R2 = 0.90) [28]. Keyak et al. [10] applied the same method for fracture simulations with stance and fall loading configurations, both resulting in good correlations with mechanical tests. For the purpose of the present study, our methods were deemed sufficiently accurate, as simulated fractures were in agreement with literature findings [17,18] and clinical practice [30,31]. The definition of the instant of fracture is a matter of controversy in the literature. Keyak defined fracture as the instant when Von Mises stresses exceed initial element strength for 15 contiguous non-surface elements [10,15]. Bessho and Koivumäki [16,17,19] defined fracture of femur models when at least one cortical (shell) surface element had failed. In our study, we assumed structural fracture when the maximum total reaction force was reached, similar to the definition used in fracture experiments [10,18,19,28]. We used a volumetric measure of failed cortical elements (500 mm3 ) to estimate this increment, as FD-curves did not always show an evident peak. A sensitivity study on the threshold volume showed that differences between fracture loads derived with various threshold volumes (500; 750; 1000 mm3 ) were negligible. Hence, the conclusions of our study do not depend heavily on the specific definition of this volumetric measure. Only one healthy and one osteoporotic bone model were used for the simulations. Differences in geometrical properties may affect the internal load distribution of the femoral bone and thereby the fracture loads and types. As geometric properties of our bone models were
598
A.M. van der Zijden et al. / Medical Engineering and Physics 37 (2015) 593–598
within ranges of the normal elderly population [26], we assumed that these two models were good general representatives of healthy and osteoporotic bones. The subjects of the in vivo experiment (trained athletes) differed from the target group. However, for safety reasons, elderly and osteoporotic people were not included in these fall experiments. We applied loading configurations based on the averages and wide ranges found in this dataset and the literature. Therefore, we expected that the loading configurations were suitable for a general population, enabling to assess trends and relative effects. Other limitations were related to the mechanical behavior of elements, as we did not include asymmetric yielding (for example the Drucker Pragner criterion [19,29]), strain rate dependent behavior [33], or anisotropic behavior [34]. Furthermore, simulation of a realistic fall impact load should involve applying a dynamic impact load instead of a quasi-static load. Further limitations concern the absence of fat and muscle tissue and the oversimplification of the pelvis model. More realistic mechanical behavior of elements, soft tissue and dynamic impact loads must be implemented in future FE models and validated against biomechanical fracture and musculoskeletal modeling studies. In summary, this study showed the potential of incorporating the loading configurations derived from in vivo fall assessments in FE simulations, applying fall impact loads distally to the GT. The fracture loads were comparable between loading configurations derived from in vivo fall data and those used in the literature, whilst from a kinematic and anatomic perspective, the loading configurations were more realistic. Simulated fracture types varied between loading configurations and differed between and within the healthy and osteoporotic bone model, reflecting its dependency on the simulated loading configuration. In conclusion, our results reveal insights in the importance of loading configurations, providing a next step in FE simulations of real-life fall-related hip fractures. Funding Nothing to declare: no external funding was received for the study. Ethical approval The bone models were part of a previous study [28], in which the models were derived from CT scans of two cadaveric femora, obtained from the Department of Anatomy with institutional approval. Fall loading configurations were derived from in vivo fall data of a previous study [24], in which each subject signed an informed consent form prior to participation. The experiment protocol was approved by the Ethical Board of the region Arnhem–Nijmegen. Conflict of interest statement The authors declare that they have no conflict of interest. References [1] Kanis JA, Oden A, McCloskey E, Johansson H, Wahl DA, Cooper C. A systematic review of hip fracture incidence and probability of fracture worldwide. Osteoporosis Int 2012;23:S308. [2] Kannus P, Parkkari J, Koskinen S, Niemi S, Palvanen M, Jarvinen M, et al. Fall-induced injuries and deaths among older adults. JAMA-J Am Med Assoc 1999;281(20):1895–9. [3] Cummings SR, Melton LJ. Epidemiology and outcomes of osteoporotic fractures. Lancet 2002;359(9319):1761–7. [4] DargentMolina P, Favier F, Grandjean H, Baudoin C, Schott AM, Hausherr E, et al. Fall-related factors and risk of hip fracture: the EPIDOS prospective study. Lancet 1996;348(9024):416. [5] Jarvinen TLN, Sievanen H, Khan KM, Heinonen A, Kannus P. Shifting the focus in fracture prevention from osteoporosis to falls. Brit Med J 2008;336(7636):124–6.
[6] Cummings-Vaughn LA, Gammack JK. Falls, osteoporosis, and hip fractures. Med Clin N Am 2011;95(3):495–506. [7] Hayes WC, Myers ER, Robinovitch SN, VandenKroonenberg A, Courtney AC, McMahon TA. Etiology and prevention of age-related hip fractures. Bone 1996;18(1):S77–86. [8] Greenspan SL, Myers ER, Maitland LA, Resnick NM, Hayes WC. Fall severity and bone-mineral density as risk-factors for hip fracture in ambulatory elderly. JAMA-J Am Med Assoc 1994;271(2):128–33. [9] Ford CM, Keaveny TM, Hayes WC. The effect of impact direction on the structural capacity of the proximal femur during falls. J Bone Miner Res 1996;11(3):377–83. [10] 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(4):539–44. [11] Backman S. The proximal end of the femur: investigations with special reference to the etiology of femoral neck fractures; anatomical studies; roentgen projections; theoretical stress calculations; experimental production of fractures. Acta Radiol. Suppl 1957(146):1–166. [12] Lotz JC, Hayes WC. The use of quantitative computed-tomography to estimate risk of fracture of the hip from falls. J Bone Joint Surg Am 1990;72A(5):689–700. [13] Courtney AC, Wachtel EF, Myers ER, Hayes WC. Effects of loading rate on strength of the proximal femur. Calcified Tissue Int 1994;55(1):53–8. [14] Keyak JH, Rossi SA, Jones KA, Skinner HB. Prediction of femoral fracture load using automated finite element modeling. J Biomech 1998;31(2):125–33. [15] Keyak JH, Rossi SA, Jones KA, Les CM, Skinner HB. Prediction of fracture location in the proximal femur using finite element models. Med Eng Phys 2001;23(9):657–64. [16] Bessho M, Ohnishi I, Okazaki H, Sato W, Kominami H, Matsunaga S, et al. Prediction of the strength and fracture location of the femoral neck by CT-based finiteelement method: a preliminary study on patients with hip fracture. J Orthop Sci 2004;9(6):545–50. [17] 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(2):226–31. [18] de Bakker PM, Manske SL, Ebacher V, Oxland TR, Cripton PA, Guy P. During sideways falls proximal femur fractures initiate in the superolateral cortex: evidence from high-speed video of simulated fractures. J Biomech 2009;42(12):1917–25. [19] Koivumaki JE, 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(4):824–9. [20] Robinovitch SN, Inkster L, Maurer J, Warnick B. Strategies for avoiding hip impact during sideways falls. J Bone Miner Res 2003;18(7):1267–73. [21] Tsouknidas A, Anagnostidis K, Maliaris G, Michailidis N. Fracture risk in the femoral hip region: a finite element analysis supported experimental approach. J Biomech 2012;45(11):1959–64. [22] Pinilla TP, Boardman KC, Bouxsein ML, Myers ER, Hayes WC. Impact direction from a fall influences the failure load of the proximal femur as much as age-related bone loss. Calcif Tissue Int 1996;58(4):231–5. [23] Choi WJ, Hoffer JA, Robinovitch SN. Effect of hip protectors, falling angle and body mass index on pressure distribution over the hip during simulated falls. Clin Biomech 2010;25(1):63–9. [24] van der Zijden AM, Groen BE, Tanck E, Nienhuis B, Verdonschot N, Weerdesteyn V. Can martial arts techniques reduce fall severity? An in vivo study of femoral loading configurations in sideways falls. J Biomech 2012;45(9):1650–5. [25] World Health Organization. Prevention and management of osteoporosis. World Health Organ Tech Rep Ser 2003;921:1–164. [26] Rubin PJ, Leyvraz PF, Aubaniac JM, Argenson JN, Esteve P, Deroguin B. The morphology of the proximal femur – a 3-dimensional radiographic analysis. J Bone Joint Surg Br 1992;74(1):28–32. [27] Keyak JH, Kaneko TS, Tehranzadeh J, Skinner HB. Predicting proximal femoral strength using structural engineering models. Clin Orthop Rel Res 2005;(437):219–28. [28] Derikx LC, van Aken JB, Janssen D, Snyers A, van der Linden YM, et al. The assessment of the risk of fracture in femora with metastatic lesions: comparing case-specific finite element analyses with predictions by clinical experts. J Bone Joint Surg Br 2012;94(8):1135–42. [29] Derikx LC, Vis R, Meinders T, Verdonschot N, Tanck E. Implementation of asymmetric yielding in case-specific finite element models improves the prediction of femoral fractures. Comput Meth Biomech 2011;14(2):183–93. [30] Zuckerman JD, Rosenberg AD. Hip fracture – reply. N Engl J Med 1996;335(26):1995–6. [31] Marks R, Allegrante JP, MacKenzie CR, Lane JM. Hip fractures among the elderly: causes, consequences and control. Ageing Res Rev 2003;2(1):57–93. [32] Laing AC, Robinovitch SN. The force attenuation provided by hip protectors depends on impact velocity, pelvic size, and soft tissue stiffness. J Biomech Eng-Trans ASME 2008;130(6):061005. [33] Hansen U, Zioupos P, Simpson R, Currey JD, Hynd D. The effect of strain rate on the mechanical properties of human cortical bone. J Biomech Eng-Trans ASME 2008;130(1):011011. [34] Keyak JH, Lee IY, Nath DS, Skinner HB. Postfailure compressive behavior of tibial trabecular bone – three anatomic directions. J Biomed Mater Res 1996;31(3):373–8.