Vertebral fractures are common in the elderly, but efforts to reduce their incidence have been hampered by incomplete understanding of the failure processes that are involved. This study's goal was to elucidate failure processes in the lumbar vertebra and to assess the accuracy of quantitative computed tomography (QCT)-based finite element (FE) simulations of these processes. Following QCT scanning, spine segments (n = 27) consisting of L1 with adjacent intervertebral disks and neighboring endplates of T12 and L2 were compressed axially in a stepwise manner. A microcomputed tomography scan was performed at each loading step. The resulting time-lapse series of images was analyzed using digital volume correlation (DVC) to quantify deformations throughout the vertebral body. While some diversity among vertebrae was observed on how these deformations progressed, common features were large strains that developed progressively in the superior third and, concomitantly, in the midtransverse plane, in a manner that was associated with spatial variations in microstructural parameters such as connectivity density. Results of FE simulations corresponded qualitatively to the measured failure patterns when boundary conditions were derived from DVC displacements at the endplate. However, quantitative correspondence was often poor, particularly when boundary conditions were simplified to uniform compressive loading. These findings suggest that variations in trabecular microstructure are one cause of the differences in failure patterns among vertebrae and that both the lack of incorporation of these variations into QCT-based FE models and the oversimplification of boundary conditions limit the accuracy of these models in simulating vertebral failure.

Introduction

Fractures of the vertebra are the most common type of fracture in elderly persons [1] and are associated with the elevated risk of hip fractures and additional vertebral fractures [24]. In terms of health care expenditures, loss of productivity, and mortality, the burden of these fractures is large: for example, women suffering vertebral fractures have a 20% greater risk of death in the subsequent 5 year period [5]. Bone mineral density (BMD), an estimate of the average density in the vertebra, is presently the cornerstone of clinical estimates of individual's risk of vertebral fracture yet is recognized as insufficient [6]. Considering that the vertebra is a complex, load-bearing structure, the search for better predictors of when vertebrae fracture (e.g., the force at which fracture will occur) would benefit from establishing how they fracture (i.e., the failure process). Defining where failure initiates, where it propagates, and what factors control these processes of initiation and propagation would identify characteristics of the vertebra beyond BMD that influence bone strength and fracture risk. These insights would feed the development of better predictors of vertebral fracture.

Although numerous prior studies have examined whether various structural features of the vertebra correlate with its stiffness and strength, only a small subset has examined how these features associate with the failure process. Features such as curvature of the anterior cortex, cortical thickness, endplate curvature, trabecular microstructure, and heterogeneity in the distribution of bone density throughout the centrum have been found to correlate with vertebral strength [712] or fracture severity [13] independent of average bone density. In regard to the failure process, studies have used sequences of radiographs or computed tomography scans acquired during mechanical testing of vertebrae to identify how the vertebra deforms as it reaches and exceeds the ultimate point. These studies have indicated that the endplate deflection is a principal feature of the failure process, and that the locations of large deflection and collapse of the endplate are associated with endplate porosity and the underlying trabecular microstructure [14]. Taken together, these studies provide evidence that the local microstructure of the vertebra influences how and where collapse of the vertebra initiates and progresses.

Failure processes in vertebrae have also been investigated using finite element (FE) modeling [1518]. However, the conclusions from these studies are varied, and the modeling approaches taken in these studies have not been fully validated for examining failure processes in bone. The need for validation is perhaps greatest for “patient-specific” FE models, i.e., models build from a quantitative computed tomography (QCT) scan of an individual's spine, because of the demand for using these models in a wide range of applications, including evaluating osteoporosis treatments [18], fracture risk [16,17], sex-related differences in bone strength [19], and surgical interventions such as vertebroplasty [20]. While QCT-based FE models have provided good predictions of vertebral strength, stiffness [16,2123], and fracture risk [18,24], their performance in replicating the failure process in vertebrae is less well established. Early work assessing this aspect of accuracy was limited to two-dimensional studies [2527], examination of only surface strains [25,2830], or qualitative comparisons of regions of large deformation [31]. A recent quantitative, three-dimensional comparison in thoracic vertebrae found errors as high as 280% in the FE predictions of local displacements [32]. Experimental studies of failure processes in the spine are, thus, needed both to develop an understanding of these processes and to gauge the fidelity of QCT-based FE models in capturing the true biomechanical behavior of the vertebra.

The overall goal of this study was to elucidate failure processes in the lumbar vertebra and to assess the accuracy of QCT-based FE simulations of these processes. Our specific objectives were (1) to quantify the strains occurring throughout the entire vertebral body during compression loading up to and beyond the yield point, (2) to determine how the failure patterns (i.e., the locations of large strains at the yield point) associate with local microstructure, and (3) to quantify accuracy of QCT-based FE predictions of the failure patterns.

Methods

Specimen Preparation.

Spine segments (n = 27) consisting of L1 with adjacent intervertebral disks and neighboring endplates of T12 and L2 were harvested from fresh-frozen human spines (age range: 41–91 years, mean±standard deviation: 80.5±10.4 years, 13 female and 14 male). The posterior elements were removed in order to allow the spine segment to fit within the mechanical loading device, as described later. The top and bottom endplates of each specimen were potted in circular dishes filled with 2–4 mm of polymethylmethacrylate (Fig. 1(a)). The spine segments were kept hydrated at all times and were wrapped in saline-soaked gauze at −20 °C when not in use.

QCT Imaging.

Each specimen was scanned using a GE Lightspeed VCT system (GE Healthcare, Waukesha, WI; 120 kV, 240 mA, 0.32 mm in-plane resolution, 0.625 mm slice thickness, bone kernel reconstruction) (Fig. 1(b)). Each specimen was held in a radiolucent fixture and immersed in degassed water. The scan axis was aligned with the superior-inferior direction of the specimen. A three-component, calcium hydroxyapatite phantom of densities 150, 75, and 0 mg/cm3 (Image Analysis, Columbia, KY) was included in the scan field of view and was used to convert intensity values into mineral densities.

Mechanical Testing and Microcomputed Tomography (μCT) Imaging.

Each specimen was placed in a custom-designed, radiolucent, loading device (Fig. 1(c)). A 22 kN load cell (LLB450, Futek Advanced Sensor Technology, Irvine, CA) was placed in the loading device beneath the specimen to record the axial force applied at each load increment. The cup was filled with 60% saline and 40% of 25% ethyl alcohol to maintain hydration of the specimen while also slowing decomposition. After ten cycles of preconditioning to 400 N, the vertebra was scanned via μCT (μCT80, Scanco Medical, Brüttisellen, Switzerland) at a resolution 37 μm/voxel. The settings for voltage, current, and integration time were 70 kVp, 114 mA, and 300 ms, respectively. The specimen was then compressed axially in a stepwise fashion by turning the top screw cap (one step = 1 mm applied at 0.25 mm/s). Following a 20 min hold period to allow for the load to equilibrate, the specimen was imaged at each step using the same settings as the initial scan. Stepwise loading was continued until the ultimate force was reached (Fig. 1(c)).

Quantification of Trabecular Microstructure.

Using the pretest μCT scans, bone volume fraction (BV/TV), apparent density (ρapp), trabecular separation (Tb.Sp*), trabecular number (Tb.N*), connectivity density (ConnD), degree of anisotropy (DA), and structural model index (SMI) were calculated (Scanco Medical) for 5-mm cubic regions located throughout the centrum, in order to quantify the spatial variation of microstructure within and among coronal, transverse, and sagittal planes (Fig. 2(a)) [33]. Measurements of trabecular thickness (Tb.Th*) were not used in this study, because the ratio of the measured values of Tb.Th* to voxel size did not exceed the minimum recommended by the system manufacturer.

Quantification of Strains Throughout the Vertebra.

Strains occurring throughout each L1 vertebra during mechanical testing were measured experimentally using a previously published method [34]. Briefly, this method entails using image registration (Scanco Medical) to align the stacks of μCT images from different load increments, creating a mask and an accompanying mesh that subdivides the vertebral body into hexahedral regions approximately 5 mm on a side (IA-FEM, The University of Iowa, Iowa City, IA) (Fig. 2(b)), and then computing the displacements and strain occurring within each region at each load increment via a custom digital volume correlation (DVC) technique [35]. This technique works by tracking the movement of the nodes of the mesh from one load increment to the next; this movement can be tracked because the porous, irregular microstructure of the vertebral trabecular bone and cortex bestows each region with a unique grayscale appearance, or “texture,” in the μCT images. The displacement field was assumed to vary linearly within each region, and the green strain was computed for each region from the displacement field. The measurement error for the displacements was estimated as 15 μm, which was the mean displacement error obtained from using the DVC technique on a pair of repeated μCT scans of a vertebra, one of which was shifted artificially by 74 μm (two voxels). The corresponding measurement error for strain was 0.0007.

The DVC measurements of displacement were also used to calculate the stiffness of the L1 vertebra. First, the average axial displacement of each L1 endplate was computed by averaging the axial displacements of the nodes located along the endplate. Then, the change in vertebral height (i.e., the total displacement applied only to the L1 vertebral body) was calculated as the difference between the average displacement of the superior and inferior L1 endplates at the yield point. The yield point was defined as the increment after which there was a change in the slope of the load–displacement curve (indicated by “C” on the load–displacement curves in Figs. 3 and 4). Vertebral stiffness was calculated as the ratio of the total applied force to the change in vertebral height.

FE Analyses.

Finite element models were generated from the QCT images following coarsening of the images to 1.24 × 1.24 × 1.25 mm3. The outer boundary of the vertebral body was defined using a semi-automated segmentation technique (Amira 5.4, Visage Imaging, Inc., San Diego, CA) (Fig. 2(c)). Each voxel was converted into a hexahedral finite element. The element size was chosen based on a previous study where numerical convergence was examined [22,36,37]. Linear, transversely isotropic elastic properties were assigned based on the local value of bone mineral density [36]. Any negative values of elastic modulus produced by this assignment, which can happen in regions with large amounts of fatty tissue, were set to 0.0001 MPa [22].

Two sets of boundary conditions were applied separately (Fig. 5). In the first, the DVC-measured values of endplate displacement at the yield point were applied to the superior and inferior endplates in the model (“experimentally matched loading”). In the second, the superior endplate in the model was uniformly displaced downward by an amount corresponding to the experimentally measured change in vertebral height at the yield point, while the inferior surface was fixed in all three directions (“uniform compression loading”). The second set of boundary conditions was, thus, a simplified, idealized version of the first, and also one that would be feasible for patient-specific FE modeling in the present day.

For comparison to the experimentally measured values of vertebral stiffness, the FE-computed stiffness of the L1 vertebra was calculated using the sum of the reaction forces in the axial direction divided by the average displacement applied at the endplates.

Statistical Analyses.

Digital volume correlation calculated strains (minimum principal strain and maximum shear strain) for each of the hexahedral volumes and microstructural properties for each of the 5 mm cubes were grouped by location within each of the three anatomic planes (Fig. 2(a)). The dependence of each of the outcome variables (e.g., a type of strain or microstructural property) on anatomic location was determined by carrying out a repeated-measure analyses of variance for each of the three anatomic groupings and, when necessary, was followed by Wilcoxon signed-rank test for post hoc pairwise comparisons (JMP 10.0, SAS Institute Inc., Cary, NC). Subsequently, to examine whether the regions of high strain occurred preferentially in locations with certain microstructural characteristics (e.g., low BV/TV), analyses of covariance were used to examine the dependence of minimum principal strain and maximum shear strain on a given microstructural property (the covariate), anatomic location (defined by transverse plane or coronal plane), and sex.

The FE-computed and experimentally measured axial displacements were compared point-by-point throughout each L1 vertebral body using linear regression and paired t-tests (JMP 10.0, SAS Institute Inc., Cary, NC). In these analyses, each measurement of displacement was paired with the average of the FE nodal displacements within the 5 mm region at the corresponding location. This region-averaging was necessary because the finite elements were smaller than the DVC regions. The effect of boundary conditions on the accuracy of the FE-computed displacements was determined by first computing, for each specimen, the median percent difference (i.e., median error) between measured and FE-computed displacements for each set of boundary conditions, and then comparing via a paired t-test, the median errors for experimentally matched versus uniform loading boundary conditions. The FE-computed and experimentally measured values of vertebral stiffness were compared using linear regression and paired t-tests.

A significance level of 0.05 was used for all statistical analyses and Bonferroni correction was used when applicable.

One sample was excluded from this study due to an error occurring during mechanical testing (male, age 82). All of the remaining samples (n = 26) were used for examining failure patterns, because the strains occurring in the vertebra exceeded the error threshold (0.0007). However, for the comparison of measured and FE-computed displacements, an additional nine samples were excluded (leaving n = 17) because the former did not exceed the error threshold (15 μm). For the comparison of measured and computed vertebral stiffnesses, only an additional seven samples were excluded (leaving n = 19), because in these seven samples the average endplate displacement did not exceed the error threshold.

Results

Progression of Vertebral Failure.

The progression of deformation throughout the vertebra fell into three general categories (Figs. 3 and 4, Supplemental Figure 1 which is available under the “Supplemental Data” tab for this paper on the Digital Collection). In the first category (12 of the 26 specimens; six male, six female), large compressive and shear strains initiated in the inferior/midtransverse cortex and periphery of the trabecular centrum and then propagated superiorly and posteriorly through the middle of the centrum and into the posterior-superior regions (Figs. 3(a) and 4(a)). In the second category (10 of the 26 specimens; five male, five female), the largest compressive and shear strains occurred initially in the central superior endplate and then became concentrated in the anterior and posterior cortices (Figs. 3(b) and 4(b)). In the third category (four of the 26 specimens; two male, two female), the large compressive strains were first seen at the superior, anterior region of the vertebra. The locations of large compressive and shear strain then propagated inward toward the middle of the centrum, and in some specimens, intensified in the anterior aspect (Figs. 3(c) and 4(c)).

Despite the differences in the progression of deformation among the three categories, however, there were some common features among them. Up to the yield point, the large strains were localized within only a small to moderately sized region of the vertebra, whereas immediately following the yield point, there was a marked increase in the total size of the regions with large strain. Averaged over all vertebrae, strains at the yield point were highest and lowest in the superior and inferior thirds of the centrum, respectively (p < 0.001; Figs. 6(a) and 6(b)). These regional differences mirrored those in Tb.N* and Tb.Sp* but not any of the other parameters of trabecular microstructure (Figs. 6(c)6(f)). The magnitudes of maximum shear strain and minimum principal strain at the yield point were correlated with one another for the majority of the specimens (n = 23, R2 = 0.70 ± 0.20, p ≤ 0.035; for the remaining n = 3, R2 = 0.08 ± 0.05, p > 0.081) (Fig. 7(a)). The direction of the minimum principal strain tended to be oblique to the loading axis (superior-inferior direction; Fig. 7(b)).

Associations Between Vertebral Failure Patterns and Spatial Distributions of Trabecular Microstructure.

The distributions of strain occurring throughout the vertebra at the yield point were associated with several features of the trabecular microstructure. Low values of ConnD and Tb.N*, and high values of Tb.Sp*, were associated with large magnitudes of compressive and shear strains regardless of whether the regions were grouped according to transverse or coronal planes (p < 0.006) (Fig. 6). High values of DA were associated with large compressive and shear strains, and low values of SMI were associated with large compressive strains, only when grouped according to coronal plane (p < 0.004) (See Supplemental Figure 2 which is available under the “Supplemental Data” tab for this paper). Further, there was an interaction between SMI and coronal plane (p = 0.008) in that the association between SMI and compressive strain was found in the anterior and middle (p < 0.032), but not posterior (p = 0.734), coronal planes. No associations were found between BV/TV and the strains in any of the anatomical planes (p > 0.093).

The inclusion of sex as a factor in the analyses of covariances that examined associations between the distributions of strains at the yield point and the distributions of trabecular microstructure allowed us to examine the possibility of differences in failure patterns between male and female vertebrae. There was a significant interaction between sex and transverse plane (p < 0.001): in the female vertebrae, compressive and shear strains were higher in the superior versus middle transverse plane and higher in the middle versus inferior transverse plane; in contrast, in the male vertebrae, compressive and shear strains were lowest in the inferior transverse plane but were not different between superior and middle transverse planes. Shear strains were higher in female versus male vertebrae in the superior transverse plane (p = 0.050). No other differences in compressive or shear strains, or in microstructure, were found between sexes (p > 0.078) (See Supplementary Material which is available under the “Supplemental Data” tab for this paper).

Accuracy of Failure Patterns Predicted by QCT-Based FE Analyses.

The failure patterns in the vertebrae were more accurately predicted by the QCT-based FE models that used experimentally matched loading than those that used uniform compression loading. Qualitative comparison of the predicted deformation fields indicated that the simulations using uniform compression loading did not reproduce as much of the spatial heterogeneity in displacement and strain (Fig. 8). R2 values for regressions of FE-computed versus experimentally measured displacements improved from 0.000 to 0.180 (p = 0.013–0.946) for uniform compressive loading to 0.175–0.842 (p < 0.001) for experimentally matched loading (Figs. 9(a) and 9(b)). In addition, the median displacement error (the median percent difference between FE-computed and measured displacements) was higher for the uniform compressive loading (range: 19.7–239.8%) versus experimentally matched loading (range: 8.5–48%) (p = 0.003; Fig. 9(c)). However, the latter set of boundary conditions still resulted in inaccuracies in the displacement predictions: paired t-tests comparing the FE-computed and measured displacements for each sample and each set of boundary conditions indicated that, for most of the samples, the computed displacements were different from their measured counterparts (experimentally matched: p < 0.043 for all except two samples, for which p > 0.351; uniform compressive: p < 0.042 for all except four samples, for which p > 0.193; Fig. 9(c)).

Experimentally Measured Stiffness Versus FE Model Prediction.

For both sets of boundary conditions, predictions of vertebral stiffness were different from (p < 0.043) and not correlated with (R2 < 0.085, p > 0.23) the measured values of stiffness (Table 1).

Discussion

The goal of this study was to quantify failure processes in the lumbar vertebra and to assess how well current implementations of QCT-based FE models of the vertebra simulate these processes. The progression of deformation exhibited by the vertebrae in this study fell into three general categories, which differed from one another primarily in the location(s) at which the largest compressive and shear strains developed up to the yield point. The analyses of trabecular microstructure indicated that the failure patterns in the vertebra at yield were associated with properties such as connectivity density and trabecular number, thus suggesting that microstructural variations are one cause of the differences in failure patterns among vertebrae. Overall, the FE simulations, which account for variations in local density but not any other features of the trabecular microstructure, did not accurately replicate the failure patterns observed experimentally, particularly when the simulations used uniform boundary conditions. Taken together, these results indicate that the QCT-based FE models may not correctly capture the structural behavior of the vertebra at the onset of failure because they do not incorporate specimen-specific information on trabecular microstructure and because of oversimplification of the boundary conditions.

The primary strength of this study is the use of an experimental method for measuring deformations throughout the vertebra over the course of the mechanical test, such that we could quantify how the vertebra deforms as failure initiates and progresses, how these deformations associate with intraspecimen variations in microstructure, and how accurately the QCT-based FE models simulate these deformations. This approach revealed that, although different categories of deformation progression were noted up until the yield point, these differences waned over subsequent loading increments (Figs. 3 and 4) in that all three categories exhibited deformation fields characterized by higher strains superiorly than inferiorly, and comparatively high compressive strains at the anterior side of the midtransverse plane. These results indicate that different vertebrae reach yield via different paths, even though their deformations postyield may be more similar. The evolution of the deformation following yield may explain why the microstructural features that correlated with the failure patterns at yield were not all the same as those found to correlate with the residual strains existing throughout the vertebra following compressive loading past the ultimate point [11]. Both studies found similar results regarding ConnD, Tb.Sp*, and Tb.N*; however, DA was negatively correlated with shear residual strains but positively correlated with shear and compressive strains at yield. It is also important to note that the experimental method for measuring deformations provided the experimentally matched boundary conditions so that we could specifically examine how the simplification to uniform compression boundary conditions affected the accuracy of the FE simulations. This examination demonstrated that this simplification does have a large impact on the fidelity of the simulations.

There are also limitations to this study. First, in the FE models, the bone tissue was modeled as linear elastic and transversely isotropic, with the elastic constants determined purely by the local value of BMD and fixed values from the literature. Further, no distinction other than that contained in the spatial variation in BMD was made among the trabecular centrum, cortical shell, and endplate. Although more sophisticated models are feasible using high-resolution peripheral QCT (HRpQCT) or microcomputed tomography (μCT), these imaging methods are not permissible for in vivo imaging of the human spine, and thus, their use was not well aligned with our goal to assess the accuracy of clinically feasible patient-specific FE models. Second, some imprecision in the matching of DVC-measured and FE-computed displacements arose from the region averaging of the latter due to the lower spatial resolution of the experimental displacement fields. As is typical for image correlation methods, the magnitudes of the errors decreased with increasing subregion size. However, as the subregion size increases, the spatial resolution of the strain measurements decreases. We, thus, chose the smallest subregion size that provided acceptable levels of accuracy and precision. Third, the posterior elements of the vertebra were removed prior to mechanical testing (and QCT imaging), rendering the loading somewhat nonphysiologic. Removal of the posterior elements is a common procedure for biomechanical testing of vertebrae, at times because of limitations on the size of test samples [21,38,39], as was the case herein due to the size constraints imposed by the μCT scanner. Finally, given the quasi-static nature of the loading protocol, the results of this study relate more directly to vertebral compression fractures caused by prolonged loading and may not apply to fractures due to short-term or dynamic loading conditions.

The data obtained in this study advance current understanding of failure processes in the vertebra by quantifying the deformations present earlier in the failure process than has been examined previously. Prior studies have focused on the ultimate point or subsequent points during the failure process [14,16,17,28,32]. These studies have found that the loss of load-carrying capacity of thoracic and lumbar vertebrae coincides with large deflection of the endplate, that the superior endplate is much more often the endplate that collapses, and that the locations of large endplate deflection are associated with many microstructural features of the trabecular bone immediately underlying the endplate. These studies have also found that well past the ultimate point, large strains occurred in the midtransverse plane and near the endplates. What is unclear from these published data, however, is whether the endplate deflection occurs due to high strains immediately under the endplate or deeper into the centrum, and whether those high strains develop progressively from early stages of loading or occur more suddenly as the ultimate point nears. The results of the present study show that, on average, the comparatively high strains in the superior third of the vertebra originate at the yield point (Figs. 6(a) and 6(b)) or earlier (Figs. 3 and 4). Although strains are also often high in some locations of the midtransverse plane prior to and at the yield point, both shear and compressive strains tend to be highest more superiorly. Thus, by piecing together the progression of deformations presented in this study with those reported previously for later stages of loading, it is apparent that collapse of the superior endplate has its origins in comparatively high compressive and shear strains that develop in the superior third of the vertebra by the yield point if not earlier. Moreover, associations between strain magnitudes and trabecular microstructure exist both at the yield and ultimate points. These findings emphasize an integral role of the microstructure underlying the superior endplate in the deformation and failure behavior of the lumbar vertebra.

Despite evidence from a number of studies that QCT-based FE models of the vertebra can provide good predictions of vertebral stiffness and strength, the QCT-based models constructed in this study did not accurately predict vertebral deformation at the yield point, particularly for the boundary conditions representing uniform compression. This latter result is consistent with an earlier study that assessed the accuracy of QCT-based FE predictions of the deformation of thoracic vertebrae immediately following the ultimate point [32]. Collectively, these results emphasize that without boundary conditions that account for the nonuniform loading of the endplate by the intervertebral disks and further improvements to the constitutive representation of the bone tissue in the models, QCT-based FE models do not accurately represent the true structural behavior of the vertebra even though they might successfully predict the overall stiffness and strength. It is further noteworthy that, in contrast to prior studies, the FE-computed values of vertebral stiffness in the present study were not correlated with the measured values (Table 1), even though the latter are in the range of previously published values [16,22]. We believe that the discrepancy with prior work arose because we loaded the vertebrae via their adjacent intervertebral disks, whereas prior studies have removed the disks and embedded the endplates in potting material prior to applying compression, thus ensuring more uniform loading across the endplate. Even though the experimentally matched boundary conditions mimicked the presence of the disks, the nonuniformity in the loading across the endplate could exacerbate errors caused by inaccuracies in the material properties used for the bone tissue.

In summary, the findings of this study indicate that, while some diversity among samples was observed in the way that deformations develop and progress throughout the lumbar vertebra, common features of the failure processes on average were large strains that develop progressively in the superior third of the vertebra and, concomitantly, in the midtransverse plane, in a manner that is associated with spatial variations in microstructure. FE simulations of models built from QCT scans of the vertebrae tested in this study provided a reasonable qualitative correspondence to the observed failure patterns when the boundary conditions were derived from the measured displacements at the endplate. However, the quantitative correspondence was poor for many vertebrae, particularly when the boundary conditions were simplified to uniform compressive loading. These results suggest that the development of clinically feasible strategies to improve the estimation of physiological boundary conditions and microstructural features of the vertebra are likely to yield major improvement in the accuracy with which QCT-based FE models replicate failure processes in the vertebra. These insights may lead to more sensitive and specific indicators of vertebral fracture risk.

Acknowledgment

The authors thank Professor Paul Barbone, Dr. Kadin Tseng, John Gallagher, Daniel Hogan, Spencer Shore, Alexander Adams, and Gabriel McDonald.

Funding Data

  • Division of Chemical, Bioengineering, Environmental, and Transport Systems (Grant No. NSF BES 0521255).

  • National Institute of Arthritis and Musculoskeletal and Skin Diseases (Grant No. NIH R01 AR054620).

  • International Osteoporosis Foundation and Servier Research Group (EFM), BU/CIMIT Applied Healthcare Fellowship (AIH).

References

References
1.
Burge
,
R.
,
Dawson‐Hughes
,
B.
,
Solomon
,
D. H.
,
Wong
,
J. B.
,
King
,
A.
, and
Tosteson
,
A.
,
2007
, “
Incidence and Economic Burden of Osteoporosis‐Related Fractures in the United States, 2005–2025
,”
J. Bone Miner. Res.
,
22
(
3
), pp.
465
475
.
2.
Black
,
D. M.
,
Arden
,
N. K.
,
Palermo
,
L.
,
Pearson
,
J.
, and
Cummings
,
S. R.
,
1999
, “
Prevalent Vertebral Deformities Predict Hip Fractures and New Vertebral Deformities but Not Wrist Fractures
,”
J. Bone Miner. Res.
,
14
(
5
), pp.
821
828
.
3.
Cauley
,
J. A.
,
Hochberg
,
M. C.
,
Lui
,
L.-Y.
,
Palermo
,
L.
,
Ensrud
,
K. E.
,
Hillier
,
T. A.
,
Nevitt
,
M. C.
, and
Cummings
,
S. R.
,
2007
, “
Long-Term Risk of Incident Vertebral Fractures
,”
JAMA
,
298
(
23
), pp.
2761
2767
.
4.
Pongchaiyakul
,
C.
,
Nguyen
,
N. D.
,
Jones
,
G.
,
Center
,
J. R.
,
Eisman
,
J. A.
, and
Nguyen
,
T. V.
,
2005
, “
Asymptomatic Vertebral Deformity as a Major Risk Factor for Subsequent Fractures and Mortality: A Long‐Term Prospective Study
,”
J. Bone Miner. Res.
,
20
(
8
), pp.
1349
1355
.
5.
Bliuc
,
D.
,
Nguyen
,
N. D.
,
Nguyen
,
T. V.
,
Eisman
,
J. A.
, and
Center
,
J. R.
,
2013
, “
Compound Risk of High Mortality Following Osteoporotic Fracture and Refracture in Elderly Women and Men
,”
J. Bone Miner. Res.
,
28
(
11
), pp.
2317
2324
.
6.
Kanis
,
J. A.
,
Johnell
,
O.
,
Oden
,
A.
,
De Laet
,
C.
,
Jonsson
,
B.
, and
Dawson
,
A.
,
2002
, “
Ten-Year Risk of Osteoporotic Fracture and the Effect of Risk Factors on Screening Strategies
,”
Bone
,
30
(
1
), pp.
251
258
.
7.
Eswaran
,
S. K.
,
Gupta
,
A.
,
Adams
,
M. F.
, and
Keaveny
,
T. M.
,
2006
, “
Cortical and Trabecular Load Sharing in the Human Vertebral Body
,”
J. Bone Miner. Res.
,
21
(
2
), pp.
307
314
.
8.
Wegrzyn
,
J.
,
Roux
,
J. P.
,
Arlot
,
M. E.
,
Boutroy
,
S.
,
Vilayphiou
,
N.
,
Guyen
,
O.
,
Delmas
,
P. D.
,
Chapurlat
,
R.
, and
Bouxsein
,
M. L.
,
2010
, “
Role of Trabecular Microarchitecture and Its Heterogeneity Parameters in the Mechanical Behavior of Ex Vivo Human L3 Vertebrae
,”
J. Bone Miner. Res.
,
25
(
11
), pp.
2324
2331
.
9.
Yerramshetty
,
J.
,
Kim
,
D. G.
, and
Yeni
,
Y. N.
,
2009
, “
Increased Microstructural Variability Is Associated With Decreased Structural Strength but With Increased Measures of Structural Ductility in Human Vertebrae
,”
ASME J. Biomech. Eng.
,
131
(9), p.
094501
.
10.
Cody
,
D. D.
,
Goldstein
,
S. A.
,
Flynn
,
M. J.
, and
Brown
,
E. B.
,
1991
, “
Correlations Between Vertebral Regional Bone Mineral Density (rBMD) and Whole Bone Fracture Load
,”
Spine
,
16
(
2
), pp.
146
154
.
11.
Hussein
,
A. I.
, and
Morgan
,
E. F.
,
2013
, “
The Effect of Intravertebral Heterogeneity in Microstructure on Vertebral Strength and Failure Patterns
,”
Osteoporosis Int.
,
24
(
3
), pp.
979
989
.
12.
Fields
,
A. J.
,
Lee
,
G. L.
, and
Keaveny
,
T. M.
,
2010
, “
Mechanisms of Initial Endplate Failure in the Human Vertebral Body
,”
J. Biomech.
,
43
(
16
), pp.
3126
3131
.
13.
Graeff
,
C.
,
Marin
,
F.
,
Petto
,
H.
,
Kayser
,
O.
,
Reisinger
,
A.
,
Peña
,
J.
,
Zysset
,
P.
, and
Glüer
,
C.-C.
,
2013
, “
High Resolution Quantitative Computed Tomography-Based Assessment of Trabecular Microstructure and Strength Estimates by Finite-Element Analysis of the Spine, but Not DXA, Reflects Vertebral Fracture Status in Men With Glucocorticoid-Induced Osteoporosis
,”
Bone
,
52
(
2
), pp.
568
577
.
14.
Jackman
,
T. M.
,
Hussein
,
A. I.
,
Adams
,
A. M.
,
Makhnejia
,
K. K.
, and
Morgan
,
E. F.
,
2014
, “
Endplate Deflection Is a Defining Feature of Vertebral Fracture and Is Associated With Properties of the Underlying Trabecular Bone
,”
J. Orthopaedic Res.
,
32
(
7
), pp.
880
886
.
15.
Faulkner
,
K. G.
,
Cann
,
C. E.
, and
Hasegawa
,
B. H.
,
1991
, “
Effect of Bone Distribution on Vertebral Strength: Assessment With a Patient-Specific Nonlinear Finite Element Analysis
,”
Radiology
,
179
(3), pp.
669
674
.
16.
Dall'Ara
,
E.
,
Pahr
,
D.
,
Varga
,
P.
,
Kainberger
,
F.
, and
Zysset
,
P.
,
2012
, “
QCT-Based Finite Element Models Predict Human Vertebral Strength In Vitro Significantly Better Than Simulated DEXA
,”
Osteoporosis Int.
,
23
(
2
), pp.
563
572
.
17.
Dall'Ara
,
E.
,
Schmidt
,
R.
,
Pahr
,
D.
,
Varga
,
P.
,
Chevalier
,
Y.
,
Patsch
,
J.
,
Kainberger
,
F.
, and
Zysset
,
P.
,
2010
, “
A Nonlinear Finite Element Model Validation Study Based on a Novel Experimental Technique for Inducing Anterior Wedge-Shape Fractures in Human Vertebral Bodies In Vitro
,”
J. Biomech.
,
43
(
12
), pp.
2374
2380
.
18.
Keaveny
,
T. M.
,
Donley
,
D. W.
,
Hoffmann
,
P. F.
,
Mitlak
,
B. H.
,
Glass
,
E. V.
, and
San Martin
,
J. A.
,
2007
, “
Effects of Teriparatide and Alendronate on Vertebral Strength as Assessed by Finite Element Modeling of QCT Scans in Women With Osteoporosis
,”
J. Bone Miner. Res.
,
22
(
1
), pp.
149
157
.
19.
Keyak
,
J.
,
Sigurdsson
,
S.
,
Karlsdottir
,
G.
,
Oskarsdottir
,
D.
,
Sigmarsdottir
,
A.
,
Zhao
,
S.
,
Kornak
,
J.
,
Harris
,
T.
,
Sigurdsson
,
G.
, and
Jonsson
,
B.
,
2011
, “
Male–Female Differences in the Association Between Incident Hip Fracture and Proximal Femoral Strength: A Finite Element Analysis Study
,”
Bone
,
48
(
6
), pp.
1239
1245
.
20.
Rohlmann
,
A.
,
Zander
,
T.
, and
Bergmann
,
G.
,
2006
, “
Spinal Loads After Osteoporotic Vertebral Fractures Treated by Vertebroplasty or Kyphoplasty
,”
Eur. Spine J.
,
15
(
8
), pp.
1255
1264
.
21.
Buckley
,
J. M.
,
Loo
,
K.
, and
Motherway
,
J.
,
2007
, “
Comparison of Quantitative Computed Tomography-Based Measures in Predicting Vertebral Compressive Strength
,”
Bone
,
40
(
3
), pp.
767
774
.
22.
Crawford
,
R. P.
,
Cann
,
C. E.
, and
Keaveny
,
T. M.
,
2003
, “
Finite Element Models Predict In Vitro Vertebral Body Compressive Strength Better Than Quantitative Computed Tomography
,”
Bone
,
33
(
4
), pp.
744
750
.
23.
Imai
,
K.
,
Ohnishi
,
I.
,
Bessho
,
M.
, and
Nakamura
,
K.
,
2006
, “
Nonlinear Finite Element Model Predicts Vertebral Bone Strength and Fracture Site
,”
Spine
,
31
(
16
), pp.
1789
1794
.
24.
Kopperdahl
,
D. L.
,
Aspelund
,
T.
,
Hoffmann
,
P. F.
,
Sigurdsson
,
S.
,
Siggeirsdottir
,
K.
,
Harris
,
T. B.
,
Gudnason
,
V.
, and
Keaveny
,
T. M.
,
2014
, “
Assessment of Incident Spine and Hip Fractures in Women and Men Using Finite Element Analysis of CT Scans
,”
J. Bone Miner. Res.
,
29
(
3
), pp.
570
580
.
25.
Bay
,
B. K.
,
Yerby
,
S. A.
,
McLain
,
R. F.
, and
Toh
,
E.
,
1999
, “
Measurement of Strain Distributions Within Vertebral Body Sections by Texture Correlation
,”
Spine
,
24
(
1
), pp.
10
17
.
26.
Kopperdahl
,
D. L.
,
Morgan
,
E. F.
, and
Keaveny
,
T. M.
,
2002
, “
Quantitative Computed Tomography Estimates of the Mechanical Properties of Human Vertebral Trabecular Bone
,”
J. Orthopaedic Res.
,
20
(
4
), pp.
801
805
.
27.
Pollintine
,
P.
,
van Tunen
,
M. S.
,
Luo
,
J.
,
Brown
,
M. D.
,
Dolan
,
P.
, and
Adams
,
M. A.
,
2010
, “
Time-Dependent Compressive Deformation of the Ageing Spine: Relevance to Spinal Stenosis
,”
Spine
,
35
(
4
), pp.
386
394
.
28.
Yerby
,
S. A.
,
Bay
,
B. K.
,
Toh
,
E.
,
McLain
,
R. F.
, and
Drews
,
M. J.
,
1998
, “
The Effect of Boundary Conditions on Experimentally Measured Trabecular Strain in the Thoracic Spine
,”
J. Biomech.
,
31
(
10
), pp.
891
897
.
29.
Silva
,
M. J.
,
Keaveny
,
T. M.
, and
Hayes
,
W. C.
,
1998
, “
Computed Tomography-Based Finite Element Analysis Predicts Failure Loads and Fracture Patterns for Vertebral Sections
,”
J. Orthopaedic Res.
,
16
(3), pp.
300
308
.
30.
Kopperdahl
,
D. L.
,
Roberts
,
A. D.
, and
Keaveny
,
T. M.
,
1999
, “
Localized Damage in Vertebral Bone Is Most Detrimental in Regions of High Strain Energy Density
,”
ASME J. Biomech. Eng.
,
121
(6), pp.
622
628
.
31.
Hosseini
,
H. S.
,
Clouthier
,
A. L.
, and
Zysset
,
P. K.
,
2014
, “
Experimental Validation of Finite Element Analysis of Human Vertebral Collapse Under Large Compressive Strains
,”
ASME J. Biomech. Eng.
,
136
(
4
), p.
041006
.
32.
Jackman
,
T. M.
,
DelMonaco
,
A. M.
, and
Morgan
,
E. F.
,
2016
, “
Accuracy of Finite Element Analyses of CT Scans in Predictions of Vertebral Failure Patterns Under Axial Compression and Anterior Flexion
,”
J. Biomech.
,
49
(
2
), pp.
267
275
.
33.
Hussein
,
A.
,
Jackman
,
T.
,
Morgan
,
S.
,
Barest
,
G.
, and
Morgan
,
E.
,
2013
, “
The Intravertebral Distribution of Bone Density: Correspondence to Intervertebral Disc Health and Implications for Vertebral Strength
,”
Osteoporosis Int.
,
24
(
12
), pp.
3021
3030
.
34.
Hussein
,
A. I.
,
Barbone
,
P. E.
, and
Morgan
,
E. F.
,
2012
, “
Digital Volume Correlation for Study of the Mechanics of Whole Bones
,”
Procedia IUTAM
,
4
, pp.
116
125
.
35.
Liu
,
L.
, and
Morgan
,
E. F.
,
2007
, “
Accuracy and Precision of Digital Volume Correlation in Quantifying Displacements and Strains in Trabecular Bone
,”
J. Biomech.
,
40
(
15
), pp.
3516
3520
.
36.
Crawford
,
R. P.
,
Rosenberg
,
W. S.
, and
Keaveny
,
T. M.
,
2003
, “
Quantitative Computed Tomography-Based Finite Element Models of the Human Lumbar Vertebral Body: Effect of Element Size on Stiffness, Damage, and Fracture Strength Predictions
,”
ASME J. Biomech. Eng.
,
125
(
4
), pp.
434
438
.
37.
Pahr
,
D. H.
, and
Zysset
,
P. K.
,
2009
, “
A Comparison of Enhanced Continuum FE With Micro FE Models of Human Vertebral Bodies
,”
J. Biomech.
,
42
(
4
), pp.
455
462
.
38.
Kopperdahl
,
D. L.
,
Pearlman
,
J. L.
, and
Keaveny
,
T. M.
,
2000
, “
Biomechanical Consequences of an Isolated Overload on the Human Vertebral Body
,”
J. Orthopaedic Res.
,
18
(5), pp.
685
690
.
39.
Mosekilde
,
L.
,
Bentzen
,
S.
,
Ørtoft
,
G.
, and
Jørgensen
,
J.
,
1989
, “
The Predictive Value of Quantitative Computed Tomography for Vertebral Body Compressive Strength and Ash Density
,”
Bone
,
10
(
6
), pp.
465
470
.

Supplementary data