|Home | About | Journals | Submit | Contact Us | Français|
Skeletal fractures associated with bone mass loss are a major clinical problem and economic burden, and lead to significant morbidity and mortality in the aging population. Clinical image based measures of bone mass show only moderate correlative strength with bone strength. However, engineering models derived from clinical image data predict bone strength with significantly greater accuracy. Currently, image-based finite element (FE) models are time consuming to construct and are non parametric. The goal of this study was to develop a parametric proximal femur FE model based on a statistical shape and density model (SSDM) derived from clinical image data. A small number of independent SSDM parameters described the shape and bone density distribution of a set of cadaver femurs and captured the variability affecting proximal femur FE strength predictions. Finally, a 3D FE model of an ”unknown” femur was reconstructed from the SSDM with an average spatial error of 0.016 mm and an average bone density error of 0.037 g/cm3.
The problem of increased risk of skeletal fractures due to bone mass loss in aging or disease is a major clinical problem leading to estimated health care costs of nearly $17 billion in the US (Burge et al. 2007, Kanis et al. 2000). Not withstanding the economic burden, non-vertebral fractures are a significant cause of morbidity and mortality in the aging population (Johnell et al. 2004, Kanis et al. 2003, Melton 2003). Thus, concerted efforts are needed to not only identify those at risk for bone fractures, but also to identify treatment strategies that can maintain the health of the skeleton with age.
Bone mineral density (BMD) measurements are widely used to determine bone strength, especially in women, and can account for up to 70% of bone strength. Attempts to directly predict fracture risk based on bone mass have been made using x-ray, ultrasound, dual energy x-ray absorptiometry (DXA), computed tomography (CT), high-resolution CT, and magnetic resonance imaging (MRI) techniques (Bauer et al. 2007, Frost et al. 2002, Genant et al. 1996, Gnudi et al. 1998, Gonnelli et al. 2005, Jiang et al. 2000, Lecouvet et al. 1997, Majumdar et al. 1999, Nicholson et al. 1997, Oden et al. 1999, Rehman et al. 2002, Ross 1997, Stewart et al. 2006, Sweeney et al. 2002, Tabensky et al. 1996, Taneichi et al. 1997). While some correlations of fracture risk with the various imaging and biochemical data have been demonstrated, they remain nonspecific and have low sensitivity (Kanis 2002a, Kanis 2002b).
Engineering models based on clinical imaging data greatly improve estimates of bone strength (Cody et al. 1999, Keyak 2001, Keyak et al. 1990, Lotz et al. 1991, Lotz et al. 1995, Testi et al. 2002). Automated finite element (FE) meshing of Quantitative CT (QCT) data has been used to predict proximal femur strength and demonstrated accuracy similar to DXA-based methods (Keyak, Meagher, Skinner and Mote 1990, Keyak et al. 1998) using linear FEA analyses, and greater accuracyusing non-linear material models (Keyak 2001). More recent applications of FE modeling to bone strength prediction incorporate smooth surface representations of the proximal femur and non-linear material models to predict failure load (Bessho et al. 2007, Bessho et al. 2004, Taddei et al. 2006) demonstrating good correlation with results from experiments.
Current applications of finite element methods are limited in that each model represents a single bone and, in order to construct a model of a similar bone from a different patient, the complete model construction procedure must be repeated; only small gains in efficiency are made when generating models of multiple bones. Furthermore, the finite element models used to date have not included parametric descriptions of bone geometry and density distribution that would allow full model definition in terms of a relatively small number of parameters. The use of parametric modeling methods would greatly reduce the time required to construct a model of an individual’s femur and, by preserving model fidelity, would not negatively impact the predictive accuracy of the model. In addition, parametric models can be used to investigate the effects of variations in specific traits on bone strength and thus could help identify those individuals at greatest risk of bone fracture.
Statistical shape modeling methods have been used to describe variability in the morphology of a population of a structure in terms of a random field representation (Cootes et al. 1994, Kaus et al. 2003, Lorenz and Krahnstöver 2000, Rueckert et al. 2003) and current applications include image segmentation, registration, object recognition, and disease diagnosis (Benameur et al. 2005, Dornaika and Ahlberg 2006, Ferrarini et al. 2006, Koikkalainen et al. 2007, Rao et al. 2006, Rueckert, Frangi and Schnabel 2003, Shan et al. 2006). Statistical shape models capture the variability of biological structures by projecting a high dimensional representation of the structure onto a lower dimensional subspace of possible shapes constructed from a population of training shapes.
The goal of this study was to implement a statistical shape modeling based method to develop a parametric finite element representation of the proximal femur shape and bone density distribution of a group of cadaveric femurs. We used this modeling methodology to identify the major components (independent combinations) of bone geometry and density that explain the majority of geometry and bone density variation within this group of cadaver femurs and we show how each independent component independently affects the predicted strength of the proximal femur. Finally, we show how this parametric modeling methodology can be used to reconstruct the three-dimensional shape and bone density distribution of a previously unseen femur.
A training set of seven fresh-frozen human female right femurs (69.9 ± 8.8 years old) (Anatomical Board of the State of Texas, San Antonio, TX; Anatomical Gifts Registry, Hanover, MD) were submerged in distilled water along with a density calibration standard (Mindways, Austin, TX) and scanned using a clinical helical CT scanner (Lightspeed Ultra, GE Healthcare, Chalfont St. Giles, UK). CT data was reconstructed with 0.488 × 0.488 × 1.25 mm voxels. An iterative thresholding scheme was used to extract the bones from the imaging data and triangulated surfaces were defined to describe the outer cortical boundary of each femur (Amira v. 4.0, Mercury Computing, Chelmsford, MA). One femur surface was arbitrarily chosen as the template surface and iterative closest point analyses were performed to register the remaining surfaces to the template surface (MATLAB R2006b, The Mathworks, Inc., Natick, MA).
A three-dimensional (3D) finite element mesh was defined (TrueGrid, v. 2.3, XYZ Scientific Applications, Inc., Livermore, CA) and projected onto the template femur surface (Figure 1). A set of landmarks was defined by identifying the 47 specific attachment points of the computational FE mesh on the template femur surface. These landmarks were mapped onto the remaining femur surfaces using a modified closest point transform (Mauch 2000). Parametric descriptions of the femur surfaces were developed by mapping each triangulated surfaces to a unit sphere (Figure 2) (Brechbuhler et al. 1995). Point-to-point correspondence was achieved between landmarks on all femurs using a hierarchical optimization approach (Davies et al. 2002) by repositioning landmarks such that landmark variance within the training set was minimized (MATLAB) (Figure 3) (Davies, Twining, Cootes, Waterton and Taylor 2002). The FE mesh was then tied to each set of optimized landmarks and projected onto the corresponding femur surfaces to generate a total of 7 individual 3D FE meshes with mesh-to-mesh correspondence across all femurs.
QCT image intensity values corresponding to each node in each FE mesh were converted to apparent ash density using a linear relationship between known K2HPO4 densities in the calibration phantom and corresponding QCT image intensities (Taddei et al. 2004) along with a relationship between K2HPO4 density and apparent ash density (Les et al. 1994).
The finite element mesh and corresponding nodal bone density values are described by a shape and density parameter vector as
where vj(xyz) are the three dimensional coordinates of the nodes in the finite element mesh, vjd is the corresponding bone density at that node, j=1,…, number of nodes in the finite element mesh, and i = 1…n=7 denote each femur in the training set (Davies, Twining, Cootes, Waterton and Taylor 2002). The mean finite element model of all femurs in the training set is then defined as:
and the correlation between individual FE meshes is given by the empirical covariance matrix (Rajamani et al. 2004).
Principal component analysis (PCA) of S results in a set of eigenvalues (λi) and eigenvectors (qi) (Rueckert, Frangi and Schnabel 2003) that are the principal directions of a shape and bone density space centered at . Each eigenvalue defines the variance of the finite element mesh and corresponding bone density values from the mean in the direction of the corresponding eigenvector. The proportion of the total variance described by each eigenvector is equal to its corresponding eigenvalue normalized by the total variance (sum of all eigenvalues). The finite element mesh with associated bone density for each bone in the training set is now described as a statistical shape and density model (SSDM):
where pv is a vector containing coordinates and apparent bone density value for all nodes in the FE model, m is the number of eigenvalues, λj, , and deviation from the average femur, , was determined as the sum of the products of a set of scalar weighting factors, cj, and SSDM standard deviations, , along the qj (eigenvector) direction (Bredbenner et al. 2007)
Multiple finite element models were developed using Eq. 4 and by setting m = mM, the number of major eigenvalues (i.e. the eigenmodes that describe greater than 75% of the total training set variability).. FE models were developed for the mean geometry and bone density distribution as well as one each investigating the effects of ±1 standard deviation (SD) for each major independent eigenmode
For each FE model, apparent ash density was determined for each element as the mean of included nodal ash densities and binned into a total of 20 apparent ash density levels. Elastic- plastic material behavior with linear hardening was assumed, and isotropic elastic and hardening moduli and ultimate stress values were determined as functions of ash density for all bone elements using empirical relationships (Keller 1994). Fall loading was simulated by fixing the distal and greater trochanter nodes in each FE model and a vertical displacement was prescribed to a flexible cup contacting the medial side of the femoral head (Figure 4). Finite element models were solved using LS-DYNA (v. 970, LSTC, Livermore, CA).
To relate the principal components of the SSDM to physical descriptors of bone geometry, five measures of proximal femur geometry were quantified to determine femoral neck axis length, femoral neck diameter, neck-shaft angle, femoral head diameter, and calcar femoral cortex width (Pulkkinen et al. 2006, Pulkkinen et al. 2004) (Figure 5). These shape descriptors were determined for the mean FE model and mean ±1 standard deviation (SD) for the each of the first three independent eigenmodes of the SSDM.
Individual finite element models of each femur in the training set were reconstructed from the SSDM using (Cootes, Hill, Taylor and Haslam 1994):
where is the mean shape and bone density model, Q contains the shape model eigenvectors qi
for n-1 of non-zero shape model eigenvalues and
The ability of the statistical shape model to reconstruct an unknown or previously unseen proximal femur was quantified using a leave one out analysis (Cootes, Hill, Taylor and Haslam 1994). In this analysis, the shape model was reformulated using all but one femur in the training set and the resulting eigenmodes were then used to reconstruct the remaining “unknown” femur using a modified version of Eq. 5:
where pu is the unknown femur shape and density vector, ′ is the mean shape and density vector of the training set leaving one femur out, and b′ are shape parameters.. The SSDM parameter values, b′, , that minimize the square root of the sum of the squares of the differences between the original shape and density vector of the unknown femur and the SSMD model predicted shape and density vector of the unknown femur (Eq. 8) were determined using non-linear optimization (Mathcad, v.13, PTC, Needham, MA). Error between the original finite element shape and density vector and the SSDM predicted shape and density vector was defined as the difference between corresponding shape and density vector elements of the original vector and the SSDM predicted vector.
Principal component analysis of the statistical shape and density model demonstrated that 78% of the variation in the set of 7 human proximal femur finite element models was captured by the first three independent modes of variation (eigenmodes) (Figure 6).
A +1 SD variation from the mean shape and density distribution due to the first eigenmode resulted in an increase in proximal femur neck axis length, a decrease in neck-shaft angle, increases in both the femoral neck and femoral head diameter, and a decrease in calcar femoral cortex width whereas a −1 SD perturbation resulted in the opposite effect (Figure 7). Shape and density variations due to the second eigenmode had the opposite effect of the first eigenmode. Shape changes of the proximal femur due to the third eignmode for the neck axis length and neck-shaft angle generally followed those changes due to the second mode of variation and in most measures, the changes were of slightly lower magnitude (Figure 8).
Variations from the mean shape and density distribution described by the major eigenmodes of the SSDM resulted in significant changes in the simulated fall loading response of the proximal femur (Figure 9). The mean shape and density model resulted in a predicted maximum force of 1590 N. The shape and density variations due to the first eigenmode resulted in an increase in predicted peak force to 1996 N for a mean −1 SD variation and a peak force of 1727 N for a mean + 1 SD variation. Shape and bone density changes due to the second eigenmode has the greatest effect on the model predicated peak force ranging from 4693 N for the mean −1 SD perturbation to 524 N for the mean + 1 SD perturbation.
The SSDM was used to reconstruct an individual femur contained within the training set with high accuracy. Using all principal components, the maximum spatial error between the SSDM derived FE model and the original finite element model of the femur was 0.009 mm with an average error of 9.0E-5 mm and the maximum bone density error was 0.0031 g/cm3 with an average error of 9.4E-6 g/cm3. Using a decreasing number of eigenmodes to reconstruct the femur finite element model, the spatial error increased to a maximum of less than 0.5 mm with an average error of 5.9E-3 mm and the bone density error increased to a maximum of 2 g/cm3 with an average of 0.013 g/cm3 when only one eigenmode was used (Figure 10). The shape and bone density distribution of a previously unseen femur was approximated using the SSDM reformulated using all but one femur FE model (defined as the “unknown” femur). The maximum spatial error between the original FE model of the unknown femur and the shape model approximated femur was 2.813 mm with an average error of 0.016 mm (Figure 11a). The maximum bone density error was 0.796 g/cm3 with an average error of 0.037 g/cm3 (Figure 11b).
A parametric finite element modeling procedure based on statistical shape modeling was developed to model complex shape and bone density variations of a group of cadaver proximal femurs. The statistical shape and density model (SSDM) compactly and efficiently describes the mean and variability in the proximal femur geometry and bone density distribution for the given set of cadaveric femurs using a reduced variable space. The SSDM was used to efficiently re-construct high-fidelity finite element models of individual femurs within the training set with a high degree of accuracy. The parametric finite element modeling procedure was also used to construct a finite element model of an “unknown” femur by varying the SSDM independent parameters to minimize the difference between the SSDM model approximated femur and the actual femur model. The traditional statistical shape modelling approach allows investigations related to geometry variation in a range of musculoskeletal applications including investigating osteoarthritis risk (Bredbenner et al. 2010)) and guiding implant design (Fitzpatrick et al. 2008). In the present study, the SSDM is constructed using a finite element model of each individual femur in the training set, differing from the traditional application of statistical shape modeling in which image data alone is used to construct the shape model. By using a 3D finite element model derived from clinical image data for each femur in the training set, the resulting SSDM can be used to directly construct high fidelity engineering models of the proximal femur; these models can then be used to evaluate the effect of specific variations in bone shape and bone density on predicted bone strength. Since the variability in shape and density distribution of the proximal femur is represented by a set of uncorrelated, independent variables defined by the eigenmodes of the statistical shape and density model, specific combinations of eigenmodes can be used to reconstruct finite element models of an individual femurs and the resulting model used to assess bone strength.
A limitation of the current SSDM is the small number of femurs in the training set. If the femurs in the SSDM training set do not adequately represent shape and density variations in the population of interest, the resulting SSDM will not adequately reconstruct unknown femurs that significantly differ from those within the training set with a high level of accuracy. To overcome this limitation, a larger number of femurs should be included in the training set. The existing SSDM can be updated by reformulating the model to include additional femurs along with the existing set.
The finite element model representing the mean − 1 SD variation in the second eigenmode resulted in an increase in neck angle length, decrease in neck-shaft angle, increase in neck diameter, and an increase in cortex width (Figure 6), all of which have been shown to decrease fracture risk (Center et al. 1998, Crabtree et al. 2002, El-Kaissi et al. 2005, Gnudi et al. 1999, Patron et al. 2006, Pulkkinen, Eckstein, Lochmuller, Kuhn and Jamsa 2006, Pulkkinen, Partanen, Jalovaara and Jamsa 2004, Specker and Binkley 2005). This independent mode of variation (eigenmode) resulted in the greatest predicted maximum load in a simulated fall configuration (Figure 7). Thus, this modeling approach directly relates variations in bone geometry and bone density to bone strength through a high fidelity predictive engineering model. The advantage of the approach described herein is that a physics based predictive finite element model is used to directly assess bone strength, rather than a statistical correlation based on epidemiological data of fracture incidents, bone density measures, and geometry in a given population.
In an investigation employing similar FE-based statistical shape modelling methods to investigate fracture risk, a SSDM was created using QCT data for 21 human cadaver femurs (although the pointwise elastic modulus distribution was used in the formulation, rather than the pointwise BMD distribution used in the present study) (Bryan et al. 2009). The multivariate parameter space of the resulting shape and modulus model was sampled using a Monte Carlo approach and a set of 1000 “virtual femurs” were created within the variability bounds of the shape and modulus model. This interesting study predicted that 28 of the 1000 FE models were at highest risk of fracture due to fall loading and that the overall percentage of cortical bone in the femur played a large role in differentiating between failed and non-failed femurs. However, it would seem that, as in the present study, variability in the multivariate parameter space is limited due to the creation of the SSDM from a small set of femurs (n=21), meaning that the 1000 virtual FE models simulate a group of individuals with relatively similar femur geometry and BMD distribution. Furthermore, it is unclear how the choice of an elastic material model, rather than a more realistic nonlinear material model, would affect the FE results during simulation of overload fracture.
In summary, a statistical shape and density modeling approach was implemented to define a parametric representation of proximal femur geometry and bone density distribution. The SSDM was used to describe the major modes of variation of proximal femur geometry and bone density distribution within the training set of cadaver femurs using a small number of independent modes. The SSDM derived finite element models demonstrated the effect of the major modes of variation in shape and bone density on predicted bone strength. By determining the optimal combination of eigenmodes, a finite element model of a previously unseen femur was constructed. The application of this method to reconstruct the three dimensional shape and bone density distribution of the proximal femur from ubiquitous clinical data such as DXA scans has the potential to significantly improve the assessment of an individual’s risk of bone fracture.
This research was supported by SwRI internal research project R9541 and NIH/NIAMS research grant AR052013.
Daniel P. Nicolella, Southwest Research Institute, 6220 Culebra Road, San Antonio, TX 78238-5166, Voice: (210) 522-3222, Fax: (210) 522-6965.
Todd L. Bredbenner, Southwest Research Institute, 6220 Culebra Road, San Antonio, TX 78238-5166, Voice: (210) 522-3565, Fax: (210) 522-6965.