Journal of Biomechanics 47 (2014) 3272–3278
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
Mapping anisotropy of the proximal femur for enhanced image based finite element analysis William S. Enns-Bray a,b,c, Jan S. Owoc a,b,c, Kyle K. Nishiyama a,b,c, Steven K. Boyd a,b,c,d,n a
Schulich School of Engineering, University of Calgary, Calgary, AB, Canada T2N 4Z6 Roger Jackson Centre for Health and Wellness Research, University of Calgary, Calgary, AB, Canada T2N 4Z6 c McCaig Institute for Bone and Joint Health, University of Calgary, Calgary, AB, Canada T2N 4Z6 d Department of Radiology, Faculty of Medicine, University of Calgary, Calgary, AB, Canada T2N 4Z6 b
art ic l e i nf o
a b s t r a c t
Article history: Accepted 18 August 2014
Finite element (FE) models of bone derived from quantitative computed tomography (QCT) rely on realistic material properties to accurately predict bone strength. QCT cannot resolve bone microarchitecture, therefore QCT-based FE models lack the anisotropy apparent within the underlying bone tissue. This study proposes a method for mapping femoral anisotropy using high-resolution peripheral quantitative computed tomography (HR-pQCT) scans of human cadaver specimens. Femur HR-pQCT images were sub-divided into numerous overlapping cubic sub-volumes and the local anisotropy was quantified using a ‘direct-mechanics’ method. The resulting directionality reflected all the major stress lines visible within the trabecular lattice, and provided a realistic estimate of the alignment of Harvesian systems within the cortical compartment. QCT-based FE models of the proximal femur were constructed with isotropic and anisotropic material properties, with directionality interpolated from the map of anisotropy. Models were loaded in a sideways fall configuration and the resulting whole bone stiffness was compared to experimental stiffness and ultimate strength. Anisotropic models were consistently less stiff, but no statistically significant differences in correlation were observed between material models against experimental data. The mean difference in whole bone stiffness between model types was approximately 26%, suggesting that anisotropy can still effect considerable change in the mechanics of proximal femur models. The under prediction of whole bone stiffness in anisotropic models suggests that the orthotropic elastic constants require further investigation. The ability to map mechanical anisotropy from high-resolution images and interpolate information into clinical-resolution models will allow testing of new anisotropic material mapping strategies. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Bone Proximal femur Finite element Anisotropy Computed tomography
1. Introduction Image-based finite element (FE) analysis of bone provides a noninvasive estimate of bone strength. FE models derived from clinical quantitative computed tomography (QCT) have been used to estimate bone mechanics, including stiffness and failure load, in physiological and non-physiological loading configurations, showing good correlation with experimental data (Dragomir-Daescu et al., 2011; Schileo et al., 2008; Keyak et al., 1998). QCT-based FE models have demonstrated superior prediction of femoral fracture compared to dual-energy x-ray absorptiometry (DXA) in both cadaver (Cody et al., 1999) and clinical studies (Orwoll et al., 2009). Patient-specific FE n Corresponding author at: McCaig Institute for Bone and Joint Health, University of Calgary, 3280 Hospital Drive NW, Room HRIC 3AC64, Calgary, AB, Canada T2N 4Z6. Tel.: þ 1 403 220 3664; fax: þ 1 403 210 9573. E-mail address: [email protected] (S.K. Boyd).
http://dx.doi.org/10.1016/j.jbiomech.2014.08.020 0021-9290/& 2014 Elsevier Ltd. All rights reserved.
models are also becoming increasingly utilized for the assessment of therapeutic interventions (Keaveny et al., 2012; Lewiecki et al., 2009). QCT-based FE models typically form a continuum of elements based on 3D image voxels that describe the bone’s geometry and density distribution. Several empirical relationships exist, summarized by Helgason et al. (2008), relating element density to isotropic tissue stiffness. However this approach does not incorporate the anisotropic behavior of the underlying microarchitecture contained within elements representing the cancellous bone compartment. The heterogeneous trabecular architecture has a varying degree of anisotropy (DA) with minimal correlation between DA and bone volume fraction (BVF) (Charlebois and Zysset, 2010; Goulet et al., 1994). Studies looking at clinically relevant loading configurations, such as a side-ways fall, have been able to explain approximately 85–90% of the variance in experimental stiffness using FE models constructed with isotropic material properties (Dragomir-Daescu et al., 2011; Koivumäki et al., 2012; Keyak et al., 1998). It is possible that the remaining variance in
W.S. Enns-Bray et al. / Journal of Biomechanics 47 (2014) 3272–3278
whole bone stiffness is partly related to unaccounted anisotropy of the microarchitecture, which is difficult to resolve in the proximal femur at clinical CT resolutions. Incorporating anisotropic material properties into image-based FE models of the proximal femur requires properly oriented principal stiffness directions for each element, and the DA that reflects the mechanics of the microarchitecture. The orthotropic DA is intrinsically described by the nine independent orthotropic elastic constants, which have thus far been derived from modal analysis (Taylor et al., 2002), image-based fabric tensors (Zysset and Curnier, 1995), or image-based FE analysis (Yang et al., 1999). Directionality on the other hand, has been measured and applied to FE models of the proximal femur in a variety of ways. Early studies applied orthotropic material properties to their models with little account for alignment between principal stiffness direction and orientation of the underlying microarchitecture, but the resulting models showed very little change compared to isotropic models (Peng et al., 2006). Meanwhile, studies accounting for directionality only in predictable skeletal regions, such as the femoral neck and shaft, found a 13% difference in maximum von Mises stress and 15% difference in nodal displacement in certain regions of the femoral neck (Yang et al., 2010). Recent studies have turned to computational methods in order to develop a robust measurement of bone directionality. Stiffness directions have been approximated as the principal stress directions calculated on a per element basis after simulating complex physiological loads during gait, including joint contact and muscle attachment forces (San Antonio et al., 2012). Alternatively, directionality can be described with fabric tensor derived from mean intercept length (MIL) analysis of micro-CT images (Harrigan and Mann, 1984). In bone mechanics, fabric provides a description of local anisotropy (Cowin, 1985), which can be combined with the element density to determine the compliance tensor for each element (Zysset and Curnier, 1995). Anisotropic femur models constructed using information from numerous local fabric measurements, has shown improved whole-bone stiffness prediction compared to isotropic models relative to a micro-CT gold standard for two specimens (Marangalou et al., 2012). Methods have also been developed for assessing structural anisotropy directly from clinical resolution images using gradient structure tensors (Tabor et al., 2013). While this technique has demonstrated good directionality agreement with fabric tensors in tissue samples from vertebral (Wolfram et al., 2009) and femoral (Kersh et al., 2013) specimens, the resulting DA measured in these studies has only shown weak correlation with fabric tensors thus far. This study proposes a new method for quantifying and mapping anisotropy of the proximal femur using the direct mechanics (DM) method (van Rietbergen et al., 1996). This approach provides a direct computational assessment of mechanical anisotropy by deriving directionality and DA of trabecular bone from FE-based stiffness tensors. This provides a robust measurement of mechanical anisotropy, which served as the original gold standard to validate other anisotropy measurements, such as MIL (Odgaard et al., 1997). The goal of this research was to develop a DM mapping methodology followed by interpolation of anisotropic data into clinical resolution, QCT-based models of the proximal femur. We hypothesize that the addition of anisotropic information will lead to a significant change in FE model behavior and improve correlation with experimental data, when compared to isotropic models.
2. Methods Seven cadaveric proximal femur specimens were fresh frozen, with all soft tissue removed, and stored in saline solution prior to scanning in DXA, QCT, and High Resolution-peripheral Quantitative Computed Tomography (HR-pQCT). The heterogeneous sample set included both healthy and osteoporotic femurs, with a
3273
Table 1 Descriptive characteristics of femur specimens. Specimen
Age
Side
Weight (kg)
Height (cm)
aBMD (g/cm2)
T-score
1 2 3 4 5 6 7
67 67 83 55 55 95 95
L R L L R R L
59 59 NA 75 75 50 50
158 158 NA 165 165 158 158
0.828 0.792 0.717 0.523 0.230 0.515 NA
0.93 1.23 3.43 5.83 4.74 3.5 NA
wide range of reported T-scores as determined from the DXA scans (Table 1). After QCT scanning (GE Discovery CT750HD, GE Healthcare; 120 kVp, 60 mAs, 512 512 matrix size) periosteal surfaces were identified with a semi-automatic contouring method (Stradwin 4.3; Cambridge, UK) (Treece et al., 2010). QCT image data was then calibrated into ash density (ρash, g/cm3) using a calibration phantom included within each scan, and in-house software developed with the Visualization Toolkit (VTK 5.10; Kitware Inc.; Clifton Park, NY). The original image voxels were resliced into 0.625 mm isometric voxels using cubic interpolation (McErlain et al., 2012). HR-pQCT scanning used the default human in vivo scanning protocol, and images were reconstructed with 82 mm isometric voxel size (XtremeCT, Scanco Medical, Brütisellen, Switzerland; 60 kVp, 1000 μA, 250 ms integration time, 750 projections). As per standard HR-pQCT analysis, a Gaussian filter was applied with sigma of 1.2 and a kernel size of 1.0 to smooth the grey-scale data to remove the inherent noise in the image data. The bone mineral phase was then extracted with a global threshold of 136.5 mg HA/cm3, which sufficiently preserved the internal trabecular architecture. An outline of the entire pipeline for mapping anisotropic material properties from HR-pQCT images, and subsequent interpolation to QCT-based FE models is provided in Fig. 1. The segmented HR-pQCT images were sub-divided over a grid into numerous overlapping cubes, with a side length of 4.92 mm (Image Processing Language v5.15, Scanco Medical, Brütisellen, Switzerland). The overlapping volumes effectively increased the sampling rate of the mapping process, where each cube overlaps 50% of the volume of the six adjacent cubes. The sample volume parameters were determined based on sensitivity testing described later in this section. Each cube provided a single local anisotropy measurement using the DM method (Van Rietbergen et al., 1996). Briefly, cubes were converted into FE models using direct voxel-to-element conversion (FAIM v6.0, Numerics88 Solutions, Calgary, Canada). Six uniaxial strains were applied to each cube model, three in compression and three in symmetric shear, to define a fourth-order stiffness tensor fully characterizing the apparent mechanical anisotropy of each cube volume. The stiffness tensor for each cube was rotated into optimal orthotropic orientation using a cost-function that minimized the non-orthotropic terms of the stiffness tensor (van Rietbergen et al., 1996). The optimized rotation matrix contains the orientation of the three orthogonal planes of symmetry, which correspond to the principal directions of strength at each measurement site. The vectors describing principal stiffness directions at each cube location formed the framework for an interpolation volume generated using radial basis functions (RBF) according to global position within the femur. A multilayer RBF algorithm was implemented with a compact (Gaussian) basis function (ALGLIB v3.7.0, ALGLIB Project; Nizhny Novgorod, Russia; LambdaV¼100, Nlayers¼ 5, Rbase 5.75 mm). The RBF provided a means to interpolate site-specific anisotropy measures (based on the cubes of HR-pQCT data) onto the QCT model elements for the entire proximal femur. QCT-based FE models were constructed by directly converting isometric, calibrated QCT voxels into 8-node hexahedral elements, which were assigned stiffnesses based on each element’s scalar density (Table 2). Two models were generated for each specimen with either isotropic or anisotropic material properties. Elements representing cancellous bone were assigned a principal stiffness based on the experimentally derived relationship described by Keller (1994), while the remaining orthotropic elastic constants were based on computationally derived relationships developed by Yang et al. (1999), similar to the recently published study by San Antonio et al. (2012). Cortical bone was assigned transversely isotropic material properties based on observations by Dong and Guo (2004). The stiffness matrix of each element in the QCT model was rotated to align with the underlying microarchitecture determined from the corresponding map of anisotropy from HRpQCT. Subsequently, boundary conditions were applied to replicate a sideways fall configuration which was performed experimentally on each bone and are described in detail elsewhere (Nishiyama et al., 2013). Simulating the experimental conditions in the FE models, the femoral head, shaft, and greater trochanter were embedded in PMMA caps with a specific depth measured in situ. Elements representing PMMA were assigned an isotropic modulus of 2500 MPa and a Poisson’s ratio of 0.3 (Lewis, 1997). The head was fixed in the loading direction with unconstrained transverse motion, the femoral shaft was completely fixed, and a 1 mm compressive displacement was applied to the surface nodes of the trochanter PMMA caps. FE models were solved (FAIM v6.0, Numerics88 Solutions; Calgary Canada) on a desktop workstation (Ubuntu v12.04, Canonical Group Ltd, London, UK; 3.1 GHz
3274
W.S. Enns-Bray et al. / Journal of Biomechanics 47 (2014) 3272–3278
HR-pQCT scan (82µm)
QCT scan (625µm)
Smoothing & t hreshold
Semi auto segmentation
Sub division into cubes
Density calibration
Direct mechanics
Finite element meshing
RBFs generated from stiffness directions
Density based orthotropic properties
Rigid image registration between HR-pQCT and QCT
Element by element interpolation of stiffness directions Fig. 1. Outline of methodology for mapping anisotropy of HR-pQCT images, and applying orthotropic material properties to QCT-based FE models. The map of anisotropy and QCT-based FE models are combined by first registering their global coordinates then interpolating the principal directions using RBFs.
Table 2 Density–elasticity relationships used to derive orthotropic and isotropic elastic constants as a function of the element ash density ρ (g/cm2) (Yang et al., 1999). Elastic constant
E3a E2 E1 G12 G23 G31 υ12 υ23 υ31
Orthotropic models
Isotropic model
Cancellous bone
Cortical bone
Cortical and cancellous bone
10,500ρ2.29 E3 0.76ρ0.09 E3 0.47ρ0.12 E3 0.26ρ0.24 E3 0.29ρ0.17 E3 0.45ρ0.18 0.27ρ 0.09 0.14ρ 0.16 0.14ρ 0.07
10,500ρ2.29 0.57 E3 0.57 E3 0.29 E3 0.2 E3 0.2 E3 0.37 0.3 0.3
10,500ρ2.29 E3 E3 E3/2(1 þ υ12) G12 G12 0.3 υ12 υ12
a The density–elasticity relationship for the primary stiffness direction was based on the observations by Keller (1994).
Intel Xeon, 4 threads, 256 GB RAM). The resulting whole bone stiffness predicted by the linear models was derived from the total reaction force calculated at the surface nodes of the trochanter divided by the applied displacement. This stiffness was also compared to the experimental ultimate strength, defined as the peak reaction force, as some evidence suggests that linear elastic stiffness estimated by FE can be a good predictor of compressive strength in bone (MacNeil and Boyd, 2008). Whole bone stiffness estimated by both isotropic and anisotropic femur models was compared to experimental data using linear regression. Significance of the correlation was determined by comparing the slope of the regression line to zero, and the Fisher transform was used to convert the correlation coefficient into z-values to detect
significant differences between material models. Bland–Altman plots provided a basis to compare the agreement between the two methods, which also describes the interchangeability between the isotropic and anisotropic modeling techniques (Bland and Altman, 1986). A sensitivity analysis was performed on sample volumes from the femoral head, neck, and intertrochantic regions, to examine the process variables most susceptible to user bias. Specifically, the segmentation threshold of HR-pQCT images and sample volume used for DM analysis were varied to determine the resulting variance on the measured directionality. Sample volume side-lengths and the global threshold were varied from 3–10 mm, and 0–300 mg HA/cm3, respectively. The threshold range includes extreme values intended to represent poor segmentation parameters. An additional sensitivity test used a range of cube overlap values (0, 15, 33, and 50%) to effectively alter the sampling rate to test the sensitivity of whole-bone stiffness to the number of local anisotropy measurements.
3. Results The processes of subdivision required approximately 5 days of computation time per bone using a single thread on an OpenVMS machine running IPL (Scanco Medical). This lengthy processing time can be attributed to the lack of parallelization of operations within the IPL framework. Subdivision of each femur specimen with 50% sub-volume overlap resulted in approximately 1.3 104 cubes. DM analysis of an entire set of cube required approximately 24 h, using three GPU cores. Comparing a frontal slice of the map of anisotropy to the cutaway of a segmented HR-pQCT image of the same specimen revealed good qualitative agreement between the measured directionality and the visible orientation of trabecular microarchitecture (Fig. 2). Observations in the sagittal and transverse plane reveal similar agreement in directionality (Fig. 3). A pattern of decreased DA was also observed in stiffer (i.e. denser) regions of trabecular bone. Although the DM approach worked well for most sub-volumes, it was problematic in the cortical bone region (i.e. bordering the periosteal surface), particularly when sub-volumes extracted from the periosteal surface were not cubic. However, in most cases, the resulting principal directions from DM measurements of cortical bone yielded a realistic directionality that reflected the alignment of the Harvesian systems. Anisotropic QCT-based FE models of the proximal femur contained approximately 10 M elements, and were solved on average in 189763 min. Plots of the elemental von Mises equivalent stress show similar stress distributions, but anisotropic models displayed overall lower magnitude of von Mises stress (Fig. 4). Whole bone stiffness predicted by isotropic and anisotropic models demonstrated good correlation with experimental whole bone stiffness ðR2ISO ¼ 0:78; R2Ani ¼ 0:79; p o0:01Þ, but poorly with experimental ultimate strength ðR2ISO ¼ 0:36; R2Ani ¼ 0:35; p o 0:01Þ (Fig. 5A and B). No statistically significant difference was detected between isotropic and anisotropic correlation coefficients, with two-tailed p-values of 0.97 and 0.99 for experimental stiffness and strength, respectively. Bland–Altman plots, on the other hand, revealed a considerable mean percent difference in stiffness (26.178.2%) between model types (Fig. 5C). Interestingly, these plots also depict a clear error bias that can be modeled as a negative linear relationship between percent difference in stiffness and average model stiffness (R2 ¼0.84). Sensitivity analyses (data not shown) revealed minimal variance in principal direction measured by DM for segmentation thresholds between 120 and 180 mg HA/cm3. Higher thresholds resulted in perforation of thin trabeculae in low density regions, leading to sudden changes in directionality. In practice, the user should recognize these improper thresholds due to the obvious lack of agreement with the reconstructed image data. Increasing the sample volume resulted in a change in principal direction that reached an asymptote representing the average directionality, as a large sample region of bone microarchitecture is included in the measurement volume. Cube lengths below 5 mm resulted in highly variable directionality due to insufficient architecture to satisfy continuum material assumptions,
W.S. Enns-Bray et al. / Journal of Biomechanics 47 (2014) 3272–3278
1
5
2
3
4
Fig. 2. General orientation of the trabecular network within the proximal femur (left), with major directions outlined describing the primary compression (1) and tension lines (2), secondary compressive (3) and tension (4) lines, and the greater trochanter lines (5). Corresponding slice of the HR-pQCT based map of anisotropy (right), with arrow glyphs represent the orientation of principal stiffness measured by direct mechanics analysis of each cube. Gaps in the image represent areas with insufficient bone architecture to complete DM analysis, which only occurs in regions of extremely low bone volume.
particularly in the low density femurs where trabecular architecture was already diminished. Lastly, the resulting whole-bone stiffness was insensitive to the anisotropy sampling rate, with a difference of 1.6% in FE predicted stiffness between models with maximum (50%) and minimum (0%) sample volume overlap in the mapping procedure.
4. Discussion This study implemented a novel approach for quantitatively mapping the anisotropy of the proximal femur, and interpolating directionality from this map into the material properties of
3275
QCT-based FE models. This was achieved by subdividing HRpQCT images into numerous cubic sub-volumes for local anisotropy measurements using DM analysis. Qualitatively, the directionality visualized by this mapping methodology displayed good alignment with the general trabecular alignment observable from the HR-pQCT scans of the proximal femurs, suggesting that this is an effective method for mapping realistic 3D principal stiffness directions throughout the femur volume. Sensitivity analyses also demonstrated minimal variation in directionality due to the segmentation threshold, sub-volume dimensions, and sampling frequency. These findings suggest that the segmentation and subdivision could feasibly be automated, and that the sampling frequency could actually be reduced (i.e. reduce the overall number of sample volumes) which would significantly hasten the mapping process. Comparing the directionality in this study to the fabric measures of Marangalou et al. (2012) shows good agreement within the femoral head, neck, and intertrochantic region. In fact, we also experimented with directionality based on fabric tensors using the MIL method and found excellent agreement with the DM results (data not shown). Together, these findings support the potential use of fabric measures to determine directionality, which are less computationally demanding than DM measurements. In the femoral head and greater trochanter we also observed a pronounced medial-lateral directionality that is notably lacking in non-image-based anisotropy maps (San Antonio et al., 2012). This pattern represents the primary tension line within the proximal femur, which is the trabecular network corresponding to compressive strength in a sideways fall scenario. The absence of medial-lateral directionality observed by San Antonio and colleagues is likely due to the boundary conditions used to derive principal stress directions which define the principal stiffness directions. Thus architecture-based methods such as DM or MIL appear to provide more accurate measurements of trabecular directionality in regions with complex microarchitecture such as the proximal femur. Consistent with other studies (San Antonio et al., 2012; Yang et al., 2010; Baca et al., 2008), we implemented different material properties to elements representing cortical and cancellous compartments, and distinguished them using a fixed apparent density threshold of 1.2 g/cm3 representing the density of fully mineralized bone tissue. Although other studies have used isotropic material properties for cortical bone (Marangalou et al., 2012), we chose to apply orthotropic properties to cortical elements, since the measured directionality in cortical regions showed good agreement with the Harvesian system directionality revealed by Baca et al. (2008). However, using a density-based threshold for identifying cortical elements could only consistently identify cortical bone in the femoral shaft, excluding the thinner cortices in the femoral head, neck, and trochanters, which are likely important under sideways fall loading configurations. It is challenging to extract cortices from CT data at clinical resolutions (Prevrhal et al., 2003; Hangartner and Gilsanz, 1996), and future studies should consider image-based techniques (Treece et al., 2012; Pahr and Zysset, 2009) for distinguishing compartments. We found no statistically significant differences between the mechanical properties of organ-level anisotropic and isotropic femur models. This suggests that accurate directionality combined with existing orthotropic density–elasticity relationships does not adequately capture additional variation in experimental bone strength, or that the variation in experimental data masks the small improvements realized by incorporating anisotropy. There was a consistent underestimation of experimental stiffness within anisotropic models, which given the similar correlation and slopes, suggests that the orthotropic elastic constants were not appropriate for this anisotropic material mapping strategy. Similar to
3276
W.S. Enns-Bray et al. / Journal of Biomechanics 47 (2014) 3272–3278
Fig. 3. Trabecular orientation and corresponding principal directions for sagittal (top) and transverse (bottom) planes.
another study (San Antonio et al., 2012), these constants were expressed as a fraction of the density–elasticity relationship proposed by Keller (1994) based on element density. This could explain both the consistent underestimation of stiffness as well as the error bias observed in the Bland–Altman plots. The bias signifies that differences between isotropic and anisotropic models increases as the whole-bone stiffness decreases, which can be explained by
the increased DA assigned to lower density bone. To overcome this problem, an alternative approach could use orthotropic elastic constants scaled by eigenvalues from normalized tensors as suggested by Zysset and Curnier (1995). Our mapping methodology provides an effective assessment of structural directionality, yet there are some important limitations. In addition to realistic stiffness directions, quantifying and
W.S. Enns-Bray et al. / Journal of Biomechanics 47 (2014) 3272–3278
FEA Predicted Stiffness (N/mm)
Fig. 4. von Mises equivalent stress field for orthotropic (left) and isotropic (right) models after 1 mm compressive displacement was applied to the PMMA cap on the greater trochanter.
Anisotropic
Isotropic
3000 2500
y = 1.23x-161 R² = 0.783
2000 1500 1000 500
y = 1.18x-498 R² = 0.792
0 500
1000
1500
2000
2500
% Difference in Stiffness (Iso-Aniso)
Experimental Stiffness (N/mm)
based on the lack of improvement compared to experimental data, we cannot recommend these orthotropic density–elasticity relationships to future studies of anisotropic organ-level FE modeling. Other similar studies have implemented a DA based on modal analysis (Taylor et al., 2002), or fabric tensors (Marangalou et al., 2012), but these femur models have not yet been experimentally validated. Alternatively, gradient structure tensors have been developed in attempt to measure patient-specific anisotropy information from clinical resolution scans, however the accuracy of the resulting DA remains a concern (Tabor et al., 2013). For future studies we will consider extracting additional information, including directionality and DA, from high-resolution images to fully realize the impact of anisotropy in QCT-based FE of the proximal femur. The advantage of the DM method in this instance is a direct evaluation of orthotropic stiffness constants throughout the bone, which could be feasibly extracted from HR-pQCT data and interpolated into QCT-based models. In conclusion, we are now able to map the mechanical anisotropy of a whole bone, which can be represented with a RBF function for easy integration into new and emerging constitutive FE models. While we initially observed little improvement by the addition of these complex material properties, it remains possible that the full potential of incorporating anisotropy may not be fully realized using the orthotropic constants we have applied. By developing a pipeline for mapping anisotropy from highresolution CT images and interpolating to clinical-resolution FE models, we are well positioned to test different anisotropic material mapping strategies and explore FE models incorporating
FEA Predicted Stiffness (N/mm)
applying bone anisotropy to QCT-based models requires an accurate DA (i.e. ratio between stiffness constants in orthogonal planes), which is a major challenge to determine. This study implemented a DA based on QCT voxel density similar to previous study (San Antonio et al., 2012). The motivation for this decision was to initially minimize the amount of information extracted from high resolution images, thus only directionality was interpolated from the map of anisotropy. The goal was to test these elastic constants with the best directional orientation possible, but
3277
Anisotropic
Isotropic
3000 y = 0.697x + 426 R² = 0.355
2500 2000 1500 1000
y = 0.658x + 95.9 R² = 0.350
500 0 500
1500
2500
3500
Experimental Ultimate Strength (N)
45 40 35
R² = 0.8383
30 25 20 15 10 5 0 0
1000
2000
3000
Average Stiffness Between Model Types (N/mm) Fig. 5. Whole bone stiffness predicted by isotropic and anisotropic QCT-based FE models compared to experimental stiffness (A) and ultimate strength (B). Linear regression lines are shown for each material model along with regression equation and coefficient of determination. A Bland–Altman plot depicts the percent difference in whole bone stiffness between model types as a function of the average stiffness between isotropic and anisotropic models (C). Horizontal lines indicate the mean difference and 95% confidence interval, and an apparent error bias is highlighted with a negative linear relationship.
3278
W.S. Enns-Bray et al. / Journal of Biomechanics 47 (2014) 3272–3278
new yield criterion to improve the predictive capabilities of FE models of the proximal femur. Conflict of interest statement None of the authors have a conflict of interest related to this work. Acknowledgements We would like to thank Mr. Eric Nodwell for technical programming assistance and development of the FAIM software for this project, and Dr. Yves Pauchard for his assistance in developing image registration software. This work was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC). References Baca, V., Horak, Z., Mikulenka, P., Dzupa, V., 2008. Comparison of an inhomogeneous orthotropic and isotropic material models used for FE analyses. Med. Eng. Phys. 30, 924–930. Bland, J.M., Altman, D.G., 1986. Statistical methods for assessing agreement between two methods of clinical measurement. Lancet 327, 308–310. Charlebois, M. Pretterklieber, Zysset, M., 2010. The role of fabric in the large strain compressive behavior of human trabecular bone. J. Biomech. Eng. 132 (1210061-9). Cody, D.D. Gross, Hou, G.J., Spencer, F.J., Goldstein, H.J., Fyhrie, D.P., S.A., 1999. Femoral strength is better predicted by finite element models than QCT and DXA. J. Biomech. 32, 1013–1020. Cowin, S.C., 1985. The relationship between the elasticity tensor and the fabric tensor. Mech. Mater. 4, 137–147. Dong, X.N., Guo, X.E., 2004. The dependence of transversely isotropic elasticity of human femoral cortical bone on porosity. J. Biomech. 37, 1281–1287. Dragomir-Daescu, D., Op Den Buijs, J., McEligot, S., Dai, Y., Entwistle, R.C., Salas, C., Melton, L.J., Bennet, K.E., Khosla, S., Amin, S., 2011. Robust QCT/FEA models of proximal femur stiffness and fracture load during a sideways fall on the hip. Ann. Biomed. Eng. 39, 742–755. Goulet, R.W., Goldstein, S.A., Ciarelli, M.J., Kuhn, J.L., Brown, M.B., Feldkamp, L.A., 1994. The relationship between structural and orthogonal compressive properties of trabecular bone. J. Biomech. 27, 375–389. Hangartner, T.N., Gilsanz, V., 1996. Evaluation of cortical bone by computed tomography. J. Bone Miner. Res. 11, 1518–1525. Harrigan., T.P., Mann., R.W., 1984. Characterization of microstructural anisotropy in orthotropic materials using a second rank tensor. J. Mater. Sci. 19 (761-761). Helgason, B., Perilli, E., Schileo, E., Fulvia, T., Sigurour, B., Viceconti, M., 2008. Mathematical relationships between bone density and mechanical properties: a literature review. Clin. Biomech. 23, 135–146. Keaveny, T.M., McClung, M.R., Wan, X., Kopperdahl, D.L., Mitlak, B.H., Krohn, K., 2012. Femoral strength in osteoporotic women treated with teriparatide or alendronate. J. Biomech. 50, 165–170. Keller, T.S., 1994. Predicting the compressive mechanical behavior of bone. J. Biomech. 27, 1159–1168. Kersh, M.E., Zysset, P.K., Pahr, D.H., Wolfram, U., Larsson, D., Pandy, M.G., 2013. Measurement of structural anisotropy in femoral trabecular bone using clinicalresolution CT images. J. Biomech. 46 (15), 2659–2666. Keyak, J.H., Rossi, S.A., Jones, K.A., Skinner, H.B., 1998. Prediction of femoral fracture load using automated finite element modeling. J. Biomech. 31, 125–133. Koivumäki, J.E.M., Thevenot, J., Pulkkinen, P., Kuhn, V., Link, T.M., Eckstein, F., Jämsä, T., 2012. Ct-based finite element models can be used to estimate experimentally measured failure loads in the proximal femur. Bone 50, 824–829.
Lewiecki, E.M., Keaveny, T.M., Kopperdahl, D.L., Genant, H.K., Engelke, K., Fuerst, T., Kivitz, A., Davies, R.Y., Fitzpatrick, L.A., 2009. Once-monthly oral ibandronate improves biomechanical determinants of bone strength in women with postmenopausal osteoporosis. J. Clin. Endocrinol. Metab. 94, 171–180. Lewis, G., 1997. Properties of acrylic bone cement: state of the art review. J. Biomed. Mater. Res. 38, 155–182. MacNeil, J.A., Boyd, S.K., 2008. Bone strength at the distal radius can be estimated from high-resolution peripheral quantitative computed tomography and the finite element method. Bone 42, 1203–1213. Marangalou, J., Ito, K., van Rietbergen, B., 2012. A new approach to determine the accuracy of morphology–elasticity relationships in continuum FE analyses of human proximal femur. J. Biomech. 45, 2884–2892. McErlain, D.D., Nishiyama, K.K., Sandino, C., Boyd, S.K., 2012.Evaluation of the Effect of CT Image Resolution on Voxel-based, Subject Specific FEA Models. [Abstract] European Society for Biomechanics. Nishiyama, K., Gilchrist, S., Guy, P., Cripton, P., Boyd, S., 2013. Proximal femur bone strength estimated by a computationally fast finite element analysis in a sideways fall configuration. J. Biomech. 46, 1231–1236. Odgaard, A., Kabel, J., van Rietbergen, B., Dalstra, M., Huiskes, R., 1997. Fabric and elastic principal directions of cancellous bone are closely related. J. Biomech. 30, 487–495. Orwoll, E.S., Marshall, L.M., Nielson, C.M., Cummings, S.R., Lapidus, J., Cauley, J.A., Ensrud, K., Lane, N., Hoffmann, P.R., Kopperdahl, D.L., Keaveny, T.M., 2009. Finite element analysis of the proximal femur and fracture risk in older men. J. Bone Miner. Res. 24, 475–483. Pahr, D.H., Zysset, P.K., 2009. From high-resolution CT data to finite element models: development of an integrated modular framework. Comput. Methods Biomech. Biomed. Eng 12, 45–57. Peng, L., Bai, J., Zeng, X., Zhou, Y., 2006. Comparison of isotropic and orthotropic material property assignments on femoral finite element models under two loading conditions. Med. Eng. Phys. 28, 227–233. Prevrhal, S., Fox, J.C., Shepherd, J.A., Genant, H.K., 2003. Accuracy of CT-based thickness measurement of thin structures: modeling of limited spatial resolution in all three dimensions. Med. Phys. 30, 1–8. San Antonio, T., Ciaccia, M., Müller-Karger, C., Casanova, E., 2012. Orientation of orthotropic material properties in a femur FE model: a method based on the principal stresses directions. Med. Eng. Phys. 34, 914–919. Schileo, E., Taddei, F., Cristofolini, L., Viceconti, M., 2008. 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. 41, 356–367. Tabor, Z., Petryniak, R., Latala, Z., Konopka, T., 2013. The potential of multi-slice computed tomography based quantification of the structural anisotropy of vertebral trabecular bone. Med. Eng. Phys. 35, 7–15. Taylor, W.R., Roland, E., Ploeg, H., Hertig, D., Klabunde, R., Warner, M.D., Hobatho, M.C., Rkotomanana, L., Clift, S.E., 2002. Determination of orthotropic bone elastic constants using FEA and modal analysis. J. Biomech. 35, 767–773. Treece, G.M., Poole, K.E.S., Gee, A.H., 2012. Imaging the femoral cortex: thickness, density and mass from clinical CT. Med. Image Anal. 16, 952–965. Treece, G.M., Gee, A.H., Mayhew, P.M., Poole, K.E., 2010. High resolution cortical bone thickness measurement from clinical CT data. Med. Image Anal. 14, 276–290. Van Rietbergen, B., Odgaard, A., Kabel, J., Huiskes, R., 1996. Direct mechanics assessment of elastic symmetries and properties of trabecular bone architecture. J. Biomech. 29, 1653–1657. Wolfram, U., Schmitz, B., Heuer, F., Reinehr, M., Wilke, H.J., 2009. Vertebral trabecular main direction can be determined from clinical CT datasets using the gradient structure tensor and not the inertia tensor—a case study. J. Biomech. 42, 1390–1396. Yang, G., Kabel, J., Van Rietbergen, B., Odgaard, A., Huiskes, R., Cowin, S., 1999. The anisotropic Hooke’s law for cancellous bone and wood. J. Elast. 53, 125–146. Yang, H., Ma, X., Guo, T., 2010. Some factors that affect the comparison between isotropic and orthotropic inhomogeneous finite element material models of femur. Med. Eng. Phys. 32, 553–560. Zysset, P.K., Curnier, A., 1995. An alternative model for anisotropic elasticity based on fabric tensors. Mech. Mater. 21, 243–250.
Contents lists available at ScienceDirect
Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com
Mapping anisotropy of the proximal femur for enhanced image based finite element analysis William S. Enns-Bray a,b,c, Jan S. Owoc a,b,c, Kyle K. Nishiyama a,b,c, Steven K. Boyd a,b,c,d,n a
Schulich School of Engineering, University of Calgary, Calgary, AB, Canada T2N 4Z6 Roger Jackson Centre for Health and Wellness Research, University of Calgary, Calgary, AB, Canada T2N 4Z6 c McCaig Institute for Bone and Joint Health, University of Calgary, Calgary, AB, Canada T2N 4Z6 d Department of Radiology, Faculty of Medicine, University of Calgary, Calgary, AB, Canada T2N 4Z6 b
art ic l e i nf o
a b s t r a c t
Article history: Accepted 18 August 2014
Finite element (FE) models of bone derived from quantitative computed tomography (QCT) rely on realistic material properties to accurately predict bone strength. QCT cannot resolve bone microarchitecture, therefore QCT-based FE models lack the anisotropy apparent within the underlying bone tissue. This study proposes a method for mapping femoral anisotropy using high-resolution peripheral quantitative computed tomography (HR-pQCT) scans of human cadaver specimens. Femur HR-pQCT images were sub-divided into numerous overlapping cubic sub-volumes and the local anisotropy was quantified using a ‘direct-mechanics’ method. The resulting directionality reflected all the major stress lines visible within the trabecular lattice, and provided a realistic estimate of the alignment of Harvesian systems within the cortical compartment. QCT-based FE models of the proximal femur were constructed with isotropic and anisotropic material properties, with directionality interpolated from the map of anisotropy. Models were loaded in a sideways fall configuration and the resulting whole bone stiffness was compared to experimental stiffness and ultimate strength. Anisotropic models were consistently less stiff, but no statistically significant differences in correlation were observed between material models against experimental data. The mean difference in whole bone stiffness between model types was approximately 26%, suggesting that anisotropy can still effect considerable change in the mechanics of proximal femur models. The under prediction of whole bone stiffness in anisotropic models suggests that the orthotropic elastic constants require further investigation. The ability to map mechanical anisotropy from high-resolution images and interpolate information into clinical-resolution models will allow testing of new anisotropic material mapping strategies. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Bone Proximal femur Finite element Anisotropy Computed tomography
1. Introduction Image-based finite element (FE) analysis of bone provides a noninvasive estimate of bone strength. FE models derived from clinical quantitative computed tomography (QCT) have been used to estimate bone mechanics, including stiffness and failure load, in physiological and non-physiological loading configurations, showing good correlation with experimental data (Dragomir-Daescu et al., 2011; Schileo et al., 2008; Keyak et al., 1998). QCT-based FE models have demonstrated superior prediction of femoral fracture compared to dual-energy x-ray absorptiometry (DXA) in both cadaver (Cody et al., 1999) and clinical studies (Orwoll et al., 2009). Patient-specific FE n Corresponding author at: McCaig Institute for Bone and Joint Health, University of Calgary, 3280 Hospital Drive NW, Room HRIC 3AC64, Calgary, AB, Canada T2N 4Z6. Tel.: þ 1 403 220 3664; fax: þ 1 403 210 9573. E-mail address: [email protected] (S.K. Boyd).
http://dx.doi.org/10.1016/j.jbiomech.2014.08.020 0021-9290/& 2014 Elsevier Ltd. All rights reserved.
models are also becoming increasingly utilized for the assessment of therapeutic interventions (Keaveny et al., 2012; Lewiecki et al., 2009). QCT-based FE models typically form a continuum of elements based on 3D image voxels that describe the bone’s geometry and density distribution. Several empirical relationships exist, summarized by Helgason et al. (2008), relating element density to isotropic tissue stiffness. However this approach does not incorporate the anisotropic behavior of the underlying microarchitecture contained within elements representing the cancellous bone compartment. The heterogeneous trabecular architecture has a varying degree of anisotropy (DA) with minimal correlation between DA and bone volume fraction (BVF) (Charlebois and Zysset, 2010; Goulet et al., 1994). Studies looking at clinically relevant loading configurations, such as a side-ways fall, have been able to explain approximately 85–90% of the variance in experimental stiffness using FE models constructed with isotropic material properties (Dragomir-Daescu et al., 2011; Koivumäki et al., 2012; Keyak et al., 1998). It is possible that the remaining variance in
W.S. Enns-Bray et al. / Journal of Biomechanics 47 (2014) 3272–3278
whole bone stiffness is partly related to unaccounted anisotropy of the microarchitecture, which is difficult to resolve in the proximal femur at clinical CT resolutions. Incorporating anisotropic material properties into image-based FE models of the proximal femur requires properly oriented principal stiffness directions for each element, and the DA that reflects the mechanics of the microarchitecture. The orthotropic DA is intrinsically described by the nine independent orthotropic elastic constants, which have thus far been derived from modal analysis (Taylor et al., 2002), image-based fabric tensors (Zysset and Curnier, 1995), or image-based FE analysis (Yang et al., 1999). Directionality on the other hand, has been measured and applied to FE models of the proximal femur in a variety of ways. Early studies applied orthotropic material properties to their models with little account for alignment between principal stiffness direction and orientation of the underlying microarchitecture, but the resulting models showed very little change compared to isotropic models (Peng et al., 2006). Meanwhile, studies accounting for directionality only in predictable skeletal regions, such as the femoral neck and shaft, found a 13% difference in maximum von Mises stress and 15% difference in nodal displacement in certain regions of the femoral neck (Yang et al., 2010). Recent studies have turned to computational methods in order to develop a robust measurement of bone directionality. Stiffness directions have been approximated as the principal stress directions calculated on a per element basis after simulating complex physiological loads during gait, including joint contact and muscle attachment forces (San Antonio et al., 2012). Alternatively, directionality can be described with fabric tensor derived from mean intercept length (MIL) analysis of micro-CT images (Harrigan and Mann, 1984). In bone mechanics, fabric provides a description of local anisotropy (Cowin, 1985), which can be combined with the element density to determine the compliance tensor for each element (Zysset and Curnier, 1995). Anisotropic femur models constructed using information from numerous local fabric measurements, has shown improved whole-bone stiffness prediction compared to isotropic models relative to a micro-CT gold standard for two specimens (Marangalou et al., 2012). Methods have also been developed for assessing structural anisotropy directly from clinical resolution images using gradient structure tensors (Tabor et al., 2013). While this technique has demonstrated good directionality agreement with fabric tensors in tissue samples from vertebral (Wolfram et al., 2009) and femoral (Kersh et al., 2013) specimens, the resulting DA measured in these studies has only shown weak correlation with fabric tensors thus far. This study proposes a new method for quantifying and mapping anisotropy of the proximal femur using the direct mechanics (DM) method (van Rietbergen et al., 1996). This approach provides a direct computational assessment of mechanical anisotropy by deriving directionality and DA of trabecular bone from FE-based stiffness tensors. This provides a robust measurement of mechanical anisotropy, which served as the original gold standard to validate other anisotropy measurements, such as MIL (Odgaard et al., 1997). The goal of this research was to develop a DM mapping methodology followed by interpolation of anisotropic data into clinical resolution, QCT-based models of the proximal femur. We hypothesize that the addition of anisotropic information will lead to a significant change in FE model behavior and improve correlation with experimental data, when compared to isotropic models.
2. Methods Seven cadaveric proximal femur specimens were fresh frozen, with all soft tissue removed, and stored in saline solution prior to scanning in DXA, QCT, and High Resolution-peripheral Quantitative Computed Tomography (HR-pQCT). The heterogeneous sample set included both healthy and osteoporotic femurs, with a
3273
Table 1 Descriptive characteristics of femur specimens. Specimen
Age
Side
Weight (kg)
Height (cm)
aBMD (g/cm2)
T-score
1 2 3 4 5 6 7
67 67 83 55 55 95 95
L R L L R R L
59 59 NA 75 75 50 50
158 158 NA 165 165 158 158
0.828 0.792 0.717 0.523 0.230 0.515 NA
0.93 1.23 3.43 5.83 4.74 3.5 NA
wide range of reported T-scores as determined from the DXA scans (Table 1). After QCT scanning (GE Discovery CT750HD, GE Healthcare; 120 kVp, 60 mAs, 512 512 matrix size) periosteal surfaces were identified with a semi-automatic contouring method (Stradwin 4.3; Cambridge, UK) (Treece et al., 2010). QCT image data was then calibrated into ash density (ρash, g/cm3) using a calibration phantom included within each scan, and in-house software developed with the Visualization Toolkit (VTK 5.10; Kitware Inc.; Clifton Park, NY). The original image voxels were resliced into 0.625 mm isometric voxels using cubic interpolation (McErlain et al., 2012). HR-pQCT scanning used the default human in vivo scanning protocol, and images were reconstructed with 82 mm isometric voxel size (XtremeCT, Scanco Medical, Brütisellen, Switzerland; 60 kVp, 1000 μA, 250 ms integration time, 750 projections). As per standard HR-pQCT analysis, a Gaussian filter was applied with sigma of 1.2 and a kernel size of 1.0 to smooth the grey-scale data to remove the inherent noise in the image data. The bone mineral phase was then extracted with a global threshold of 136.5 mg HA/cm3, which sufficiently preserved the internal trabecular architecture. An outline of the entire pipeline for mapping anisotropic material properties from HR-pQCT images, and subsequent interpolation to QCT-based FE models is provided in Fig. 1. The segmented HR-pQCT images were sub-divided over a grid into numerous overlapping cubes, with a side length of 4.92 mm (Image Processing Language v5.15, Scanco Medical, Brütisellen, Switzerland). The overlapping volumes effectively increased the sampling rate of the mapping process, where each cube overlaps 50% of the volume of the six adjacent cubes. The sample volume parameters were determined based on sensitivity testing described later in this section. Each cube provided a single local anisotropy measurement using the DM method (Van Rietbergen et al., 1996). Briefly, cubes were converted into FE models using direct voxel-to-element conversion (FAIM v6.0, Numerics88 Solutions, Calgary, Canada). Six uniaxial strains were applied to each cube model, three in compression and three in symmetric shear, to define a fourth-order stiffness tensor fully characterizing the apparent mechanical anisotropy of each cube volume. The stiffness tensor for each cube was rotated into optimal orthotropic orientation using a cost-function that minimized the non-orthotropic terms of the stiffness tensor (van Rietbergen et al., 1996). The optimized rotation matrix contains the orientation of the three orthogonal planes of symmetry, which correspond to the principal directions of strength at each measurement site. The vectors describing principal stiffness directions at each cube location formed the framework for an interpolation volume generated using radial basis functions (RBF) according to global position within the femur. A multilayer RBF algorithm was implemented with a compact (Gaussian) basis function (ALGLIB v3.7.0, ALGLIB Project; Nizhny Novgorod, Russia; LambdaV¼100, Nlayers¼ 5, Rbase 5.75 mm). The RBF provided a means to interpolate site-specific anisotropy measures (based on the cubes of HR-pQCT data) onto the QCT model elements for the entire proximal femur. QCT-based FE models were constructed by directly converting isometric, calibrated QCT voxels into 8-node hexahedral elements, which were assigned stiffnesses based on each element’s scalar density (Table 2). Two models were generated for each specimen with either isotropic or anisotropic material properties. Elements representing cancellous bone were assigned a principal stiffness based on the experimentally derived relationship described by Keller (1994), while the remaining orthotropic elastic constants were based on computationally derived relationships developed by Yang et al. (1999), similar to the recently published study by San Antonio et al. (2012). Cortical bone was assigned transversely isotropic material properties based on observations by Dong and Guo (2004). The stiffness matrix of each element in the QCT model was rotated to align with the underlying microarchitecture determined from the corresponding map of anisotropy from HRpQCT. Subsequently, boundary conditions were applied to replicate a sideways fall configuration which was performed experimentally on each bone and are described in detail elsewhere (Nishiyama et al., 2013). Simulating the experimental conditions in the FE models, the femoral head, shaft, and greater trochanter were embedded in PMMA caps with a specific depth measured in situ. Elements representing PMMA were assigned an isotropic modulus of 2500 MPa and a Poisson’s ratio of 0.3 (Lewis, 1997). The head was fixed in the loading direction with unconstrained transverse motion, the femoral shaft was completely fixed, and a 1 mm compressive displacement was applied to the surface nodes of the trochanter PMMA caps. FE models were solved (FAIM v6.0, Numerics88 Solutions; Calgary Canada) on a desktop workstation (Ubuntu v12.04, Canonical Group Ltd, London, UK; 3.1 GHz
3274
W.S. Enns-Bray et al. / Journal of Biomechanics 47 (2014) 3272–3278
HR-pQCT scan (82µm)
QCT scan (625µm)
Smoothing & t hreshold
Semi auto segmentation
Sub division into cubes
Density calibration
Direct mechanics
Finite element meshing
RBFs generated from stiffness directions
Density based orthotropic properties
Rigid image registration between HR-pQCT and QCT
Element by element interpolation of stiffness directions Fig. 1. Outline of methodology for mapping anisotropy of HR-pQCT images, and applying orthotropic material properties to QCT-based FE models. The map of anisotropy and QCT-based FE models are combined by first registering their global coordinates then interpolating the principal directions using RBFs.
Table 2 Density–elasticity relationships used to derive orthotropic and isotropic elastic constants as a function of the element ash density ρ (g/cm2) (Yang et al., 1999). Elastic constant
E3a E2 E1 G12 G23 G31 υ12 υ23 υ31
Orthotropic models
Isotropic model
Cancellous bone
Cortical bone
Cortical and cancellous bone
10,500ρ2.29 E3 0.76ρ0.09 E3 0.47ρ0.12 E3 0.26ρ0.24 E3 0.29ρ0.17 E3 0.45ρ0.18 0.27ρ 0.09 0.14ρ 0.16 0.14ρ 0.07
10,500ρ2.29 0.57 E3 0.57 E3 0.29 E3 0.2 E3 0.2 E3 0.37 0.3 0.3
10,500ρ2.29 E3 E3 E3/2(1 þ υ12) G12 G12 0.3 υ12 υ12
a The density–elasticity relationship for the primary stiffness direction was based on the observations by Keller (1994).
Intel Xeon, 4 threads, 256 GB RAM). The resulting whole bone stiffness predicted by the linear models was derived from the total reaction force calculated at the surface nodes of the trochanter divided by the applied displacement. This stiffness was also compared to the experimental ultimate strength, defined as the peak reaction force, as some evidence suggests that linear elastic stiffness estimated by FE can be a good predictor of compressive strength in bone (MacNeil and Boyd, 2008). Whole bone stiffness estimated by both isotropic and anisotropic femur models was compared to experimental data using linear regression. Significance of the correlation was determined by comparing the slope of the regression line to zero, and the Fisher transform was used to convert the correlation coefficient into z-values to detect
significant differences between material models. Bland–Altman plots provided a basis to compare the agreement between the two methods, which also describes the interchangeability between the isotropic and anisotropic modeling techniques (Bland and Altman, 1986). A sensitivity analysis was performed on sample volumes from the femoral head, neck, and intertrochantic regions, to examine the process variables most susceptible to user bias. Specifically, the segmentation threshold of HR-pQCT images and sample volume used for DM analysis were varied to determine the resulting variance on the measured directionality. Sample volume side-lengths and the global threshold were varied from 3–10 mm, and 0–300 mg HA/cm3, respectively. The threshold range includes extreme values intended to represent poor segmentation parameters. An additional sensitivity test used a range of cube overlap values (0, 15, 33, and 50%) to effectively alter the sampling rate to test the sensitivity of whole-bone stiffness to the number of local anisotropy measurements.
3. Results The processes of subdivision required approximately 5 days of computation time per bone using a single thread on an OpenVMS machine running IPL (Scanco Medical). This lengthy processing time can be attributed to the lack of parallelization of operations within the IPL framework. Subdivision of each femur specimen with 50% sub-volume overlap resulted in approximately 1.3 104 cubes. DM analysis of an entire set of cube required approximately 24 h, using three GPU cores. Comparing a frontal slice of the map of anisotropy to the cutaway of a segmented HR-pQCT image of the same specimen revealed good qualitative agreement between the measured directionality and the visible orientation of trabecular microarchitecture (Fig. 2). Observations in the sagittal and transverse plane reveal similar agreement in directionality (Fig. 3). A pattern of decreased DA was also observed in stiffer (i.e. denser) regions of trabecular bone. Although the DM approach worked well for most sub-volumes, it was problematic in the cortical bone region (i.e. bordering the periosteal surface), particularly when sub-volumes extracted from the periosteal surface were not cubic. However, in most cases, the resulting principal directions from DM measurements of cortical bone yielded a realistic directionality that reflected the alignment of the Harvesian systems. Anisotropic QCT-based FE models of the proximal femur contained approximately 10 M elements, and were solved on average in 189763 min. Plots of the elemental von Mises equivalent stress show similar stress distributions, but anisotropic models displayed overall lower magnitude of von Mises stress (Fig. 4). Whole bone stiffness predicted by isotropic and anisotropic models demonstrated good correlation with experimental whole bone stiffness ðR2ISO ¼ 0:78; R2Ani ¼ 0:79; p o0:01Þ, but poorly with experimental ultimate strength ðR2ISO ¼ 0:36; R2Ani ¼ 0:35; p o 0:01Þ (Fig. 5A and B). No statistically significant difference was detected between isotropic and anisotropic correlation coefficients, with two-tailed p-values of 0.97 and 0.99 for experimental stiffness and strength, respectively. Bland–Altman plots, on the other hand, revealed a considerable mean percent difference in stiffness (26.178.2%) between model types (Fig. 5C). Interestingly, these plots also depict a clear error bias that can be modeled as a negative linear relationship between percent difference in stiffness and average model stiffness (R2 ¼0.84). Sensitivity analyses (data not shown) revealed minimal variance in principal direction measured by DM for segmentation thresholds between 120 and 180 mg HA/cm3. Higher thresholds resulted in perforation of thin trabeculae in low density regions, leading to sudden changes in directionality. In practice, the user should recognize these improper thresholds due to the obvious lack of agreement with the reconstructed image data. Increasing the sample volume resulted in a change in principal direction that reached an asymptote representing the average directionality, as a large sample region of bone microarchitecture is included in the measurement volume. Cube lengths below 5 mm resulted in highly variable directionality due to insufficient architecture to satisfy continuum material assumptions,
W.S. Enns-Bray et al. / Journal of Biomechanics 47 (2014) 3272–3278
1
5
2
3
4
Fig. 2. General orientation of the trabecular network within the proximal femur (left), with major directions outlined describing the primary compression (1) and tension lines (2), secondary compressive (3) and tension (4) lines, and the greater trochanter lines (5). Corresponding slice of the HR-pQCT based map of anisotropy (right), with arrow glyphs represent the orientation of principal stiffness measured by direct mechanics analysis of each cube. Gaps in the image represent areas with insufficient bone architecture to complete DM analysis, which only occurs in regions of extremely low bone volume.
particularly in the low density femurs where trabecular architecture was already diminished. Lastly, the resulting whole-bone stiffness was insensitive to the anisotropy sampling rate, with a difference of 1.6% in FE predicted stiffness between models with maximum (50%) and minimum (0%) sample volume overlap in the mapping procedure.
4. Discussion This study implemented a novel approach for quantitatively mapping the anisotropy of the proximal femur, and interpolating directionality from this map into the material properties of
3275
QCT-based FE models. This was achieved by subdividing HRpQCT images into numerous cubic sub-volumes for local anisotropy measurements using DM analysis. Qualitatively, the directionality visualized by this mapping methodology displayed good alignment with the general trabecular alignment observable from the HR-pQCT scans of the proximal femurs, suggesting that this is an effective method for mapping realistic 3D principal stiffness directions throughout the femur volume. Sensitivity analyses also demonstrated minimal variation in directionality due to the segmentation threshold, sub-volume dimensions, and sampling frequency. These findings suggest that the segmentation and subdivision could feasibly be automated, and that the sampling frequency could actually be reduced (i.e. reduce the overall number of sample volumes) which would significantly hasten the mapping process. Comparing the directionality in this study to the fabric measures of Marangalou et al. (2012) shows good agreement within the femoral head, neck, and intertrochantic region. In fact, we also experimented with directionality based on fabric tensors using the MIL method and found excellent agreement with the DM results (data not shown). Together, these findings support the potential use of fabric measures to determine directionality, which are less computationally demanding than DM measurements. In the femoral head and greater trochanter we also observed a pronounced medial-lateral directionality that is notably lacking in non-image-based anisotropy maps (San Antonio et al., 2012). This pattern represents the primary tension line within the proximal femur, which is the trabecular network corresponding to compressive strength in a sideways fall scenario. The absence of medial-lateral directionality observed by San Antonio and colleagues is likely due to the boundary conditions used to derive principal stress directions which define the principal stiffness directions. Thus architecture-based methods such as DM or MIL appear to provide more accurate measurements of trabecular directionality in regions with complex microarchitecture such as the proximal femur. Consistent with other studies (San Antonio et al., 2012; Yang et al., 2010; Baca et al., 2008), we implemented different material properties to elements representing cortical and cancellous compartments, and distinguished them using a fixed apparent density threshold of 1.2 g/cm3 representing the density of fully mineralized bone tissue. Although other studies have used isotropic material properties for cortical bone (Marangalou et al., 2012), we chose to apply orthotropic properties to cortical elements, since the measured directionality in cortical regions showed good agreement with the Harvesian system directionality revealed by Baca et al. (2008). However, using a density-based threshold for identifying cortical elements could only consistently identify cortical bone in the femoral shaft, excluding the thinner cortices in the femoral head, neck, and trochanters, which are likely important under sideways fall loading configurations. It is challenging to extract cortices from CT data at clinical resolutions (Prevrhal et al., 2003; Hangartner and Gilsanz, 1996), and future studies should consider image-based techniques (Treece et al., 2012; Pahr and Zysset, 2009) for distinguishing compartments. We found no statistically significant differences between the mechanical properties of organ-level anisotropic and isotropic femur models. This suggests that accurate directionality combined with existing orthotropic density–elasticity relationships does not adequately capture additional variation in experimental bone strength, or that the variation in experimental data masks the small improvements realized by incorporating anisotropy. There was a consistent underestimation of experimental stiffness within anisotropic models, which given the similar correlation and slopes, suggests that the orthotropic elastic constants were not appropriate for this anisotropic material mapping strategy. Similar to
3276
W.S. Enns-Bray et al. / Journal of Biomechanics 47 (2014) 3272–3278
Fig. 3. Trabecular orientation and corresponding principal directions for sagittal (top) and transverse (bottom) planes.
another study (San Antonio et al., 2012), these constants were expressed as a fraction of the density–elasticity relationship proposed by Keller (1994) based on element density. This could explain both the consistent underestimation of stiffness as well as the error bias observed in the Bland–Altman plots. The bias signifies that differences between isotropic and anisotropic models increases as the whole-bone stiffness decreases, which can be explained by
the increased DA assigned to lower density bone. To overcome this problem, an alternative approach could use orthotropic elastic constants scaled by eigenvalues from normalized tensors as suggested by Zysset and Curnier (1995). Our mapping methodology provides an effective assessment of structural directionality, yet there are some important limitations. In addition to realistic stiffness directions, quantifying and
W.S. Enns-Bray et al. / Journal of Biomechanics 47 (2014) 3272–3278
FEA Predicted Stiffness (N/mm)
Fig. 4. von Mises equivalent stress field for orthotropic (left) and isotropic (right) models after 1 mm compressive displacement was applied to the PMMA cap on the greater trochanter.
Anisotropic
Isotropic
3000 2500
y = 1.23x-161 R² = 0.783
2000 1500 1000 500
y = 1.18x-498 R² = 0.792
0 500
1000
1500
2000
2500
% Difference in Stiffness (Iso-Aniso)
Experimental Stiffness (N/mm)
based on the lack of improvement compared to experimental data, we cannot recommend these orthotropic density–elasticity relationships to future studies of anisotropic organ-level FE modeling. Other similar studies have implemented a DA based on modal analysis (Taylor et al., 2002), or fabric tensors (Marangalou et al., 2012), but these femur models have not yet been experimentally validated. Alternatively, gradient structure tensors have been developed in attempt to measure patient-specific anisotropy information from clinical resolution scans, however the accuracy of the resulting DA remains a concern (Tabor et al., 2013). For future studies we will consider extracting additional information, including directionality and DA, from high-resolution images to fully realize the impact of anisotropy in QCT-based FE of the proximal femur. The advantage of the DM method in this instance is a direct evaluation of orthotropic stiffness constants throughout the bone, which could be feasibly extracted from HR-pQCT data and interpolated into QCT-based models. In conclusion, we are now able to map the mechanical anisotropy of a whole bone, which can be represented with a RBF function for easy integration into new and emerging constitutive FE models. While we initially observed little improvement by the addition of these complex material properties, it remains possible that the full potential of incorporating anisotropy may not be fully realized using the orthotropic constants we have applied. By developing a pipeline for mapping anisotropy from highresolution CT images and interpolating to clinical-resolution FE models, we are well positioned to test different anisotropic material mapping strategies and explore FE models incorporating
FEA Predicted Stiffness (N/mm)
applying bone anisotropy to QCT-based models requires an accurate DA (i.e. ratio between stiffness constants in orthogonal planes), which is a major challenge to determine. This study implemented a DA based on QCT voxel density similar to previous study (San Antonio et al., 2012). The motivation for this decision was to initially minimize the amount of information extracted from high resolution images, thus only directionality was interpolated from the map of anisotropy. The goal was to test these elastic constants with the best directional orientation possible, but
3277
Anisotropic
Isotropic
3000 y = 0.697x + 426 R² = 0.355
2500 2000 1500 1000
y = 0.658x + 95.9 R² = 0.350
500 0 500
1500
2500
3500
Experimental Ultimate Strength (N)
45 40 35
R² = 0.8383
30 25 20 15 10 5 0 0
1000
2000
3000
Average Stiffness Between Model Types (N/mm) Fig. 5. Whole bone stiffness predicted by isotropic and anisotropic QCT-based FE models compared to experimental stiffness (A) and ultimate strength (B). Linear regression lines are shown for each material model along with regression equation and coefficient of determination. A Bland–Altman plot depicts the percent difference in whole bone stiffness between model types as a function of the average stiffness between isotropic and anisotropic models (C). Horizontal lines indicate the mean difference and 95% confidence interval, and an apparent error bias is highlighted with a negative linear relationship.
3278
W.S. Enns-Bray et al. / Journal of Biomechanics 47 (2014) 3272–3278
new yield criterion to improve the predictive capabilities of FE models of the proximal femur. Conflict of interest statement None of the authors have a conflict of interest related to this work. Acknowledgements We would like to thank Mr. Eric Nodwell for technical programming assistance and development of the FAIM software for this project, and Dr. Yves Pauchard for his assistance in developing image registration software. This work was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC). References Baca, V., Horak, Z., Mikulenka, P., Dzupa, V., 2008. Comparison of an inhomogeneous orthotropic and isotropic material models used for FE analyses. Med. Eng. Phys. 30, 924–930. Bland, J.M., Altman, D.G., 1986. Statistical methods for assessing agreement between two methods of clinical measurement. Lancet 327, 308–310. Charlebois, M. Pretterklieber, Zysset, M., 2010. The role of fabric in the large strain compressive behavior of human trabecular bone. J. Biomech. Eng. 132 (1210061-9). Cody, D.D. Gross, Hou, G.J., Spencer, F.J., Goldstein, H.J., Fyhrie, D.P., S.A., 1999. Femoral strength is better predicted by finite element models than QCT and DXA. J. Biomech. 32, 1013–1020. Cowin, S.C., 1985. The relationship between the elasticity tensor and the fabric tensor. Mech. Mater. 4, 137–147. Dong, X.N., Guo, X.E., 2004. The dependence of transversely isotropic elasticity of human femoral cortical bone on porosity. J. Biomech. 37, 1281–1287. Dragomir-Daescu, D., Op Den Buijs, J., McEligot, S., Dai, Y., Entwistle, R.C., Salas, C., Melton, L.J., Bennet, K.E., Khosla, S., Amin, S., 2011. Robust QCT/FEA models of proximal femur stiffness and fracture load during a sideways fall on the hip. Ann. Biomed. Eng. 39, 742–755. Goulet, R.W., Goldstein, S.A., Ciarelli, M.J., Kuhn, J.L., Brown, M.B., Feldkamp, L.A., 1994. The relationship between structural and orthogonal compressive properties of trabecular bone. J. Biomech. 27, 375–389. Hangartner, T.N., Gilsanz, V., 1996. Evaluation of cortical bone by computed tomography. J. Bone Miner. Res. 11, 1518–1525. Harrigan., T.P., Mann., R.W., 1984. Characterization of microstructural anisotropy in orthotropic materials using a second rank tensor. J. Mater. Sci. 19 (761-761). Helgason, B., Perilli, E., Schileo, E., Fulvia, T., Sigurour, B., Viceconti, M., 2008. Mathematical relationships between bone density and mechanical properties: a literature review. Clin. Biomech. 23, 135–146. Keaveny, T.M., McClung, M.R., Wan, X., Kopperdahl, D.L., Mitlak, B.H., Krohn, K., 2012. Femoral strength in osteoporotic women treated with teriparatide or alendronate. J. Biomech. 50, 165–170. Keller, T.S., 1994. Predicting the compressive mechanical behavior of bone. J. Biomech. 27, 1159–1168. Kersh, M.E., Zysset, P.K., Pahr, D.H., Wolfram, U., Larsson, D., Pandy, M.G., 2013. Measurement of structural anisotropy in femoral trabecular bone using clinicalresolution CT images. J. Biomech. 46 (15), 2659–2666. Keyak, J.H., Rossi, S.A., Jones, K.A., Skinner, H.B., 1998. Prediction of femoral fracture load using automated finite element modeling. J. Biomech. 31, 125–133. Koivumäki, J.E.M., Thevenot, J., Pulkkinen, P., Kuhn, V., Link, T.M., Eckstein, F., Jämsä, T., 2012. Ct-based finite element models can be used to estimate experimentally measured failure loads in the proximal femur. Bone 50, 824–829.
Lewiecki, E.M., Keaveny, T.M., Kopperdahl, D.L., Genant, H.K., Engelke, K., Fuerst, T., Kivitz, A., Davies, R.Y., Fitzpatrick, L.A., 2009. Once-monthly oral ibandronate improves biomechanical determinants of bone strength in women with postmenopausal osteoporosis. J. Clin. Endocrinol. Metab. 94, 171–180. Lewis, G., 1997. Properties of acrylic bone cement: state of the art review. J. Biomed. Mater. Res. 38, 155–182. MacNeil, J.A., Boyd, S.K., 2008. Bone strength at the distal radius can be estimated from high-resolution peripheral quantitative computed tomography and the finite element method. Bone 42, 1203–1213. Marangalou, J., Ito, K., van Rietbergen, B., 2012. A new approach to determine the accuracy of morphology–elasticity relationships in continuum FE analyses of human proximal femur. J. Biomech. 45, 2884–2892. McErlain, D.D., Nishiyama, K.K., Sandino, C., Boyd, S.K., 2012.Evaluation of the Effect of CT Image Resolution on Voxel-based, Subject Specific FEA Models. [Abstract] European Society for Biomechanics. Nishiyama, K., Gilchrist, S., Guy, P., Cripton, P., Boyd, S., 2013. Proximal femur bone strength estimated by a computationally fast finite element analysis in a sideways fall configuration. J. Biomech. 46, 1231–1236. Odgaard, A., Kabel, J., van Rietbergen, B., Dalstra, M., Huiskes, R., 1997. Fabric and elastic principal directions of cancellous bone are closely related. J. Biomech. 30, 487–495. Orwoll, E.S., Marshall, L.M., Nielson, C.M., Cummings, S.R., Lapidus, J., Cauley, J.A., Ensrud, K., Lane, N., Hoffmann, P.R., Kopperdahl, D.L., Keaveny, T.M., 2009. Finite element analysis of the proximal femur and fracture risk in older men. J. Bone Miner. Res. 24, 475–483. Pahr, D.H., Zysset, P.K., 2009. From high-resolution CT data to finite element models: development of an integrated modular framework. Comput. Methods Biomech. Biomed. Eng 12, 45–57. Peng, L., Bai, J., Zeng, X., Zhou, Y., 2006. Comparison of isotropic and orthotropic material property assignments on femoral finite element models under two loading conditions. Med. Eng. Phys. 28, 227–233. Prevrhal, S., Fox, J.C., Shepherd, J.A., Genant, H.K., 2003. Accuracy of CT-based thickness measurement of thin structures: modeling of limited spatial resolution in all three dimensions. Med. Phys. 30, 1–8. San Antonio, T., Ciaccia, M., Müller-Karger, C., Casanova, E., 2012. Orientation of orthotropic material properties in a femur FE model: a method based on the principal stresses directions. Med. Eng. Phys. 34, 914–919. Schileo, E., Taddei, F., Cristofolini, L., Viceconti, M., 2008. 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. 41, 356–367. Tabor, Z., Petryniak, R., Latala, Z., Konopka, T., 2013. The potential of multi-slice computed tomography based quantification of the structural anisotropy of vertebral trabecular bone. Med. Eng. Phys. 35, 7–15. Taylor, W.R., Roland, E., Ploeg, H., Hertig, D., Klabunde, R., Warner, M.D., Hobatho, M.C., Rkotomanana, L., Clift, S.E., 2002. Determination of orthotropic bone elastic constants using FEA and modal analysis. J. Biomech. 35, 767–773. Treece, G.M., Poole, K.E.S., Gee, A.H., 2012. Imaging the femoral cortex: thickness, density and mass from clinical CT. Med. Image Anal. 16, 952–965. Treece, G.M., Gee, A.H., Mayhew, P.M., Poole, K.E., 2010. High resolution cortical bone thickness measurement from clinical CT data. Med. Image Anal. 14, 276–290. Van Rietbergen, B., Odgaard, A., Kabel, J., Huiskes, R., 1996. Direct mechanics assessment of elastic symmetries and properties of trabecular bone architecture. J. Biomech. 29, 1653–1657. Wolfram, U., Schmitz, B., Heuer, F., Reinehr, M., Wilke, H.J., 2009. Vertebral trabecular main direction can be determined from clinical CT datasets using the gradient structure tensor and not the inertia tensor—a case study. J. Biomech. 42, 1390–1396. Yang, G., Kabel, J., Van Rietbergen, B., Odgaard, A., Huiskes, R., Cowin, S., 1999. The anisotropic Hooke’s law for cancellous bone and wood. J. Elast. 53, 125–146. Yang, H., Ma, X., Guo, T., 2010. Some factors that affect the comparison between isotropic and orthotropic inhomogeneous finite element material models of femur. Med. Eng. Phys. 32, 553–560. Zysset, P.K., Curnier, A., 1995. An alternative model for anisotropic elasticity based on fabric tensors. Mech. Mater. 21, 243–250.
