|Home | About | Journals | Submit | Contact Us | Français|
Osteoporotic vertebral body fractures are an increasing clinical problem among the aging population. Specimen-specific finite element models, derived from quantitative computed tomography (QCT) have the potential to more accurately predict failure loads in the vertebra. Additionally, the use of extended finite element modeling (X-FEM) allows for a detailed analysis of crack initiation and propagation in various materials. Our aim was to study the feasibility of QCT/X-FEM analysis to predict fracture properties of vertebral bodies. Three cadaveric specimens were obtained and the L3 vertebrae were excised. The vertebrae were CT scanned to develop computational models and mechanically tested in compression to measure failure load, stiffness and to observe crack location. One vertebra was used for calibration of the material properties from experimental results and CT grey scale values. The two additional specimens were used to assess the model prediction. The resulting QCT/X-FEM model of the specimen used for calibration had 2% and 4% errors in stiffness and failure load, respectively, compared with the experiment. The predicted failure loads of the additional two vertebrae were larger by about 41-44% when compared to the measured values, while the stiffness differed by 129% and 40%. The predicted fracture patterns matched fairly well with the visually observed experimental cracks. Our feasibility study indicated that the QCT/X-FEM method used to predict vertebral compression fractures is a promising tool to consider in future applications for improving vertebral fracture risk prediction in the elderly.
Osteoporotic vertebral fractures are an increasing problem among the elderly population and are thought to result from normal daily activities due to a decrease in bone mass and altered trabecular architecture. Once a fracture develops there is an increased risk for subsequent fractures which are often associated with significant physiological and behavioral function limitations [23,38]. Dual energy x-ray absorptiometry (DXA) and quantitative computed tomography (QCT) are the two most widely-used diagnostic methods for measuring bone mineral density (BMD) and predict which patients are at a risk for vertebral fracture [2,22,27,30,37,38]. However, a significant number of individuals with osteoporotic fractures have BMD values above the threshold for osteoporosis . Furthermore, DXA does not account for 3D vertebral geometry, regional variability in bone density nor identifies the site or region at risk of fracture .
Quantitative computed tomography-based finite element models (QCT/FEM) account for vertebral geometry, architecture and heterogeneous distribution of BMD within the bone . The extended finite element method (X-FEM) first introduced by Belytschko and Blackard  allows for the analysis of crack initiation without the need for repeated re-meshing or explicit geometric modeling of the discontinuity during crack propagation. Some studies have implemented X-FEM analysis to study dental, ceramics and brittle materials [3,26,33,43], cracks upon impact on windshields , the effects and characteristics of cracks on 2D structures , hip fracture pattern and repair  or cortical bone damage . However, X-FEM has not been previously used to study complex structures such as cadaveric vertebral compression fractures. There is an urgent need to better identify high-risk individuals for treatment and this requires application of novel approaches. Thus, the overall purpose of this study was to create a vertebral compression model to better predict failure load and location. Specifically, our objectives were to develop a specimen-specific QCT/X-FEM model of the lumbar vertebra (L3) implementing density dependent-linear elastic properties derived from a CT calibration phantom, and a user-defined compression failure criterion. Assessing fracture risk by using a patient specific CT-based X-FEM model could help therapies to be tailored to the specific fracture risk, location and pattern present in the individual.
Three fresh frozen cadaveric specimens were obtained (Mayo Clinic Foundation, Anatomy Department) after internal approval from the bio-specimen committee. Radiographs of the torso were acquired and visually screened by a trained clinician to ensure that no history of spine trauma or pathologies affecting bone quality (ex: fractures or cancer) other than postmenopausal osteoporosis was present in the spine. Areal bone mineral density (aBMD in g/cm2) of the thoracolumbar spinal segment (T12-L4) was determined prior to dissection of the specimens, using the General Electric Lunar Prodigy DXA equipment. Measurements were obtained in the anterior-posterior condition and a degree of osteoporosis (T-score) was assigned.
The L3 vertebrae were excised from the specimens and all soft tissue and posterior elements removed. To restrain the specimens for QCT imaging scanning and experimental testing, both upper and lower regions were embedded in polymethylmethacrylate (PMMA), as previously described . The vertebral bodies were wrapped, kept moist with physiologic saline solution and stored at −20°C until scanning and testing were performed.
Before scanning, the specimens were thawed for 24 hours. A calibration phantom (Midways Inc., San Francisco, CA) was placed beneath the samples during the scanning process. The phantom contains five rods of reference materials allowing Hounsfield units (HU) to be converted to equivalent K2HPO4 density (ρash, g/cm3). The vertebrae and phantom were imaged in air at room temperature using a Siemens Somatom Definition scanner (Siemens, Malvern, PA) and the scanner was operated at 120 kVp, 450 mAs. The images were reconstructed using a sharp kernel (U70) and 0.4mm isotropic voxels were obtained.
The vertebrae were fractured using a Mini Bionix 858 servohydraulic test machine (MTS, Eden Prairie, MN). Testing was performed at room temperature and the specimens were compressed at a rate of 5 mm/min [8,11,12,15,29] to 25% decrease of vertebral height . The potted specimens were loaded in compression between two aluminum platens to allow for uniform vertical displacements . The 25% reduction in vertebral body height was chosen as it represents a grade 1 mild compression fracture as described by Genant et al. . Force and displacement data were collected at 102.4 Hz. Compressive failure load was defined as the peak force in the load-displacement curve and stiffness was calculated in the linear portion of the curve between 30-70% peak-force, as previously described .
The QCT/X-FEM process is summarized in Fig. 1. In order to account for discontinuities across the crack tip and along the faces, X-FEM incorporates enrichment terms, extra degrees of freedom at selected nodes, in the form :
where, I represents the set of nodes in the mesh; Ni(x) is the standard finite element shape function associated with node i, and ui is the classical nodal displacement at node i. The first additional term involves the Heaviside “jump or signed” function (H(x)), with its corresponding degrees of freedom bj, which determines whether the nodes are above or below the crack when the element is bisected:
This Heaviside function introduces eventual discontinuities across the crack faces. The second additional term (Fl(x)) models the displacement around the crack-tip with its respective degrees of freedom :
QCT DICOM images of the vertebrae were imported into Mimics 14.12 (Materialise Inc., Ann Arbor, MI). Thresholding and segmentation routines were used to create a solid mask of the vertebral bodies. An initial automatic segmentation was performed using a threshold of 400 HU; each slice was then manually edited to include the entire cortical and trabecular regions and to exclude any soft tissue at the surface of the vertebrae . A 3D volume model was generated from the segmented trabecular and cortical bone. Hexahedral, 0.8mm isotropic, elements were created from the 3D models and material properties assigned to each voxel. Hounsfield units in the QCT scans were converted to apparent mineral density (ρapp) values using a linear regression model obtained from the calibration phantom used in the imaging process. A stiffness sensitivity analysis evaluating the robustness of the QCT/X-FEM model to changes in allocation of discrete material properties was performed on one specimen as previously described . Briefly, a custom Matlab™ (The Mathworks, Inc. Natick, MA, USA) script was implemented to group the CT grey values of the voxels into a discrete number of bins that approximated the continuous density distribution. To find an appropriate number of discrete material property bins, the HU values were grouped in 8, 18, 42 and 50 equal sized bins. The first bin was assigned a density of ρash = 0.01 g/cm3 to avoid conversion of negative HU values to unphysical negative densities [9,14]. The mean HU number of each finite element was averaged from the values of the contained voxels. Isotropic Young's modulus in each element (Eq. 1) was calculated as a function of apparent density (ρapp) [24,35] and based on the ratio ρapp/ρash of 0.6 at the center of each bin [14,20,25]:
The average, minimum and maximum Young's modulus in the models was also analyzed. The maximum Young's modulus was estimated from the material bins that contained any elements with densities larger than 1 g/cm3. Poisson's ratio for each element was set to 0.4, as used in previous studies [7,22,24,28] and yield strain (εy) was related to ρash as follows (Eq. 2) :
The constant value of 0.0065 in Eq. (2) was determined by a trial and error optimization procedure to improve the agreement between the experimentally measured and QCT/X-FEM predicted failure load in the same specimen used for calibration of the material . Two additional specimens were used to assess the model. In order to simulate failure of the vertebrae, a user-defined failure criterion, f (Eq. 3), was implemented to analyze crack initiation and propagation based on yield strain (εy), ρash and Minimum Normal Principal Strain (MINNPE):
Where, εz is the compressive strain in the vertical, z, direction; and εfail is the strain at which the bone tissue fails in compression.
The compressive fracture criterion implemented allowed the elements in the model that reached a MINNPE that exceeded the yield strain to fail and the crack to propagate perpendicular to the loading direction (Fig. 2). A power-law energy based damage evolution law (GIC (mode I), GIIC (mode II) and GIIIC (mode III)) governed crack propagation. In order to allow for crack to propagate once the element failed these rates were set to 0.0001 (J/mm2), assuming brittle cortical and trabecular bone. In this case, the crack propagated in the XY plane once the compressive strain reached the criterion in equation 3. The failure criterion defined in equation 3 was implemented using the FORTRAN user subroutine UDMGINI in conjunction with X-FEM using ABAQUS v.6.14 (Dassault Systèmes Simulia Corp., Providence, RI) standard for crack propagation modeling.
Models with density-dependent linear isotropic element properties were run in ABAQUS using boundary conditions that mimicked the experimental testing. The inferior vertebral surface nodes embedded in PMMA were constrained in all directions. The superior surface nodes embedded in PMMA were constrained in the mediolateral/anteroposterior translations and all three axes of rotation while allowing a vertical displacement . To mimic loading of the platen on the vertebrae, reference nodes were visually placed respectively at the superior and inferior endplates of the specimens. Superior and inferior surface nodes from the vertebrae and the reference nodes were selected to form rigid bodies and create a kinematic coupling (Fig 1). A vertical displacement to the superior reference node was applied and compression similar to the experimental testing was modeled.
In order to prepare for the X-FEM analysis, enrichment regions were identified in the models. To avoid subjectivity in the selection of the regions, and to allow for cortical and trabecular bone to fail independently, two regions were chosen based on the anatomical characteristics of the vertebrae. One enrichment region, defined as cortical bone, consisted of the outer surface elements of the vertebral bodies. A second region consisted of the remaining elements within the vertebrae, namely trabecular bone. Elements built-in the boundary conditions were excluded from these selections (Fig 1). The compressive failure criterion based on equations (2) and (3) was applied to both enrichment regions. With element failure, mimicking bone damage and fracture, a non-contact interface between the newly formed surfaces allowed for crack propagation during compression of the vertebrae. Predicted failure load, measured as the reaction force at the inferior reference node of the models, was considered as the peak force before the first large drop in the load-displacement curve. Stiffness was calculated from the initial linear portion of the load-displacement curve.
A summary of the specimens, vertebrae and DXA outcomes is shown in Table 1. One specimen had a T-score of 0.9 (normal) in the thoracolumbar spinal segment, while its individual L3 vertebra presented a T-score of +2.5 (normal, aBMD: 1.527 g/cm2). The other two specimens had a T-score of −2.6 (osteoporotic) and −2.9 (osteoporotic) in the thoracolumbar spine, while their individual L3 vertebra presented a T-score of −1.9 (osteopenic, aBMD: 0.980 g/cm2) and −1.9 (osteopenic, aBMD: 0.975 g/cm2), respectively.
Predicted stiffness of the calibrated model was sensitive to the number of materials assigned through the binning process when compared to the experimental value. A 21% difference was observed when using 8 bins, 6% for 18 bins but less than 1% for 42 or 50 bins, leading us to use of 42 material bins in the QCT/X-FEM modeling process. Load-displacement curves for the experimental and modeling procedures are shown in Fig.3. The experimental measured and predicted compressive failure loads and stiffness outcomes are presented in Table 2. Material property distribution (average, minimum and maximum Young's modulus) is also described. QCT/X-FEM of specimen 2, used for model calibration, showed good agreement with the experimental stiffness (2% error) and failure load (4% error). The predicted stiffness of specimens 1 and 3, used for model assessment, varied to a greater extent compared to the measured values. The model of specimen 1over predicted the stiffness and failure load by 129% and 41%, respectively, while the model of specimen 3 under predicted the stiffness by 40% and over predicted the failure load by 44%. Fig. 4 shows the experimental and predicted crack location and failure pattern for the specimens. Predicted fracture patterns matched fairly well with the visually observed experimental cracks, with the vertebrae failing closer to the superior vertebral rim. For specimen 2, the predicted and experimental cracks closely matched. For specimen 1, the experimental crack was observed at ~60% of the vertebral measured from the bottom, while the predicted crack was placed at ~80% of the height. For specimen 3, the experimental crack was at ~80% of the height compared with ~95% of the height for the predicted failure location. In all three cases the cortical region of the vertebra was predicted to fail first, with a subsequent propagation of the crack through the cortex and internally into the trabecular bone structure.
Segmentation and meshing of the vertebrae required 1-2 days of operator's time and the QCT/X-FEM technique involved only an additional selection of the enrichment regions and minor pre-processing modifications associated with the user-defined failure criterion and material assignment. The models had a computing time of about 8 hrs. when using a workstation with Intel Xeon L5520 (8 CPUs) with 24 GB RAM memory.
This study introduced a finite element model of vertebral compression fracture analyzed with the extended finite element method. We created specimen-specific QCT/X-FEM models of the lumbar vertebra (L3) implementing density dependent-linear elastic properties derived from computed tomography and applied a user-defined density dependent- compressive strain failure criterion to study vertebral failure load and crack location. The QCT/X-FEM method allowed for global vertebral failure load analysis as well as local evidence of the origin of the crack and its propagation. One vertebra (specimen 2) from the study was used to calibrate the material property assignment strategies and study the performance of our QCT/X-FEM model. Because yield strains differ between sites , an original strain-based equation for femoral fracture  was further optimized to better match the failure characteristic of the vertebra. These model parameters were then used to assess the QCT/X-FEM prediction of stiffness, failure load and failure location in a separate set of two vertebrae.
The predicted failure loads of the vertebrae used for model assessment were larger by about 41-44% when compared to the measured values, while the predicted stiffness differed to a greater extent. It should be noted that these values are highly affected by the density-linear elastic material property relations . Predicted values were obtained by applying previously published empirical equations which might have increased our estimated error. Experimental and QCT/X-FEM loading curves observed in Fig. 3 show the characteristic response of the models compared to the measured outcomes. There is a large drop in the predicted loads due to failure of the elements. This response is due to the non-contact interface between the newly formed surfaces of the models. Although X-FEM can be used to model propagation with contact interactions, the contact force between the crack generated from this type of compression loading profile and failure mechanism will prevent the crack from propagating.
While others have looked at bone fractures using post-processing methods, the element deletion , cohesive element approach , or the reduced-stiffness methods, these techniques do not model the discontinuities presented by the propagating crack. In the element deletion method, discontinuities associated with a cracked element are modeled by modifying the constitutive relations and setting the failed element stresses and response in the model to zero . Similarly, the stiffness reduction method decreases the element's stiffness upon reaching a failure criterion [18,21]. An element with reduced stiffness can still withstand forces and distort the analysis. Also, these methods do not always accurately represent the fracture location, pattern or propagation when using a user-defined failure criterion. Furthermore, element size has a significant effect in the outcomes of these models as the failed elements still sustain structural forces and can present large distortions due to their reduced stiffness, eventually leading to convergence issues. To compensate for these shortcomings and still predict accurately the failure loads, would require developing a finer mesh with smaller elements which is more computationally expensive. The QCT/X-FEM analysis overcomes these problems by handling the discontinuities of the crack independent of element size.
The efficiency of strain patterns and distributions as a predictor of vertebral bone failure has previously been described . Although the failure criterion was tuned in one vertebra to better match vertebral strain values, QCT/X-FEM analysis in the other two specimens did not involve any additional input from the user to predict the experimental failure patterns and loads. Ali et al. implemented the X-FEM technique to predict femoral fracture in specimen-specific models, observing a crack and failure in the cortical bone region which was close to the experimental outcome . The study showed predicted errors in failure loads and stiffness of 69% and 40%, respectively. The femur models used in  presented tensile-dominated loading and failure, while the failure criterion used in our study was based on compressive damage. Similarly, Feerick et al. implemented the X-FEM technique to predict cortical bone damage and crack propagation using an anisotropic model during screw pullout. The specimen-specific QCT/X-FEM analysis performed in this study resulted in failure patterns which were similar to the observed experimental failure characteristics of the vertebrae. More importantly, in our current analysis the assumption of a starter crack was not used, thus the origin and initiation of the crack was not forced to start at a specific region, it rather developed naturally. A study by Fradet et al. predicted vertebral fractures on a spine segment under varying kinematic conditions . The study reported that under compression loading, fracture initiated either on the endplates or in the cortical wall of the vertebra, similar to the crack origination of our study.
Vertebral fractures are classified as mild (Grade 1: 20-25% deformity), moderate (Grade 2: 25-40%) and severe (Grade 3: 40% or greater deformation) . Vertebral fractures at a young age are most often encountered in men due to trauma or injury. However, at older ages, although still common in men, the prevalence shifts towards a postmenopausal population , with an initial mild fracture, similar to the fracture model used in this study, representing many of the fractures associated with early osteoporotic deformities . Vertebral failure load might not be the only parameter to consider in fracture risk prediction. As depicted in Fig. 3, a small sudden drop in the linear curve indicates small cracks of bone failure occurring at load levels below the ultimate failure load of the vertebra. This crack formation might also indicate the initial local point of failure and its propagation pattern. This could be of importance to clinicians when targeting patient-specific therapeutics, for example in the case of vertebroplasty, in which localized interventions could be applied to a specific region of the weak vertebrae instead of the whole region. Fang et al. reported failure strength of PMMA-treated vertebra to be ~1.5 times larger compared to the intact control . However, the PMMA-induced stiffness mismatch between the treated and intact vertebrae will induce adjacent disc degeneration and a vertebral fracture cascade. Compared to traditional cement augmentation techniques, local reinforcement may provide the treated vertebrae with a closer to normal structural behavior.
Yield strain and Young's modulus values in our models depended on the measured HU values from the CT images. The trabecular compartment, filled with fat and marrow content, lowered the average, locally measured grey value and ash density of bone tissue, thus leading to smaller Young's modulus and larger yield strains, as described by the power law equations. In contrast, dense bone in the cortex, was represented by a larger modulus and smaller yield strains. Partial volume effects can highly affect the grey values of the voxels in the images. Although the average Young's modulus values of the three specimens were similar, their maximum values showed a larger variability. Voxels at the boundary of the vertebrae will contain both, cortical and trabecular bone, thus resulting in a mixed value of such voxel. Similarly, although care was taken to not include voxels outside the vertebral cortex, partial volume effects related to averaging of the CT voxel values between the cortex and the surrounding air could also bring these values down from those corresponding to only cortical tissue.
Our study has several limitations. First, X-FEM analysis on ABAQUS v.6.14 allowed only a single crack per failed element and the direction of the crack could not change by more than 90 degrees inside the element. Second, our user-defined failure criterion permitted for cohesive-type crack propagation allowing the crack to move and propagate in only one plane. Third, we are presenting data of three specimens and future studies with a larger sample size should be considered to further validate the process over a wider range of bone density distributions (normal/osteopenic/osteoporotic). The stiffness was under-predicted in one case and over-predicted in another specimen. This is related to the material equations of choice which will highly impact the results. Future studies with a larger sample size will also allow for optimization of the fracture criterion based on vertebral yield strain, aid a material properties analysis to find the best-fit elastic modulus estimation equation and propose a more robust estimation of non-linear material properties to capture the post-yield behavior of the vertebrae. A mesh size sensitivity study was not performed to find a mesh independent of crack propagation. However, the element size is small enough for prediction of a steady crack trajectory while at the same time not being computationally expensive. While enriching the entire vertebral body increases the computation time, enriching only a smaller local area might fail to capture fractures originating in other regions. Furthermore, with computers and imaging techniques constantly evolving we expect the run time of the models to not be the limiting factor in a clinical setting. One of the models (specimen 3) experienced convergence issues at large displacements due to distortions of some elements. However this occurred beyond the limits of structural failure of the vertebra. While this study shows the crack pattern and stiffness of the specimens using a linear approach, plastic deformation could be further modeled by modifying the energy release rates during failure, capturing the non-linear behavior of the loading characteristics and avoiding an abrupt decrease in the loading profile. Finally, the models allowed for tracking of crack initiation and propagation. However, we only observed the final location of the experimental process and did not record the crack origination and propagation during testing. Although experimental fracture and cracks will originate before a 25% reduction in height, it was not possible to observe this failure location during the experiment without involving imaging techniques such as micro computed tomography, which is clinically unfeasible. For this reason, the vertebrae were compressed to a 25% reduction in height to ascertain final crack location. Future studies should quantitatively measure predicted crack location and propagation in the vertebrae and compare to the experimental measured cracks.
In summary, we have introduced the first specimen-specific QCT/X-FEM model of the human vertebra to predict failure load, stiffness and fracture location under compression loading. Fracture patterns are subject-specific as these depend on bone density, quality and geometrical characteristics which are not captured by generic models. While this is a first attempt to use QCT/X-FEM in the prediction of osteoporotic vertebral compression fractures, the model could be adapted to other loading conditions and could become a valuable tool to consider in future applications improving fracture risk prediction in our aging population.
Funding: This project was supported by OREF (Orthopedic Research and Education Foundation), SRS (Scoliosis Research Society) and internal funding from the Mayo Foundation. The authors would like to acknowledge the Opus CT Imaging Resource of Mayo Clinic (NIH construction grant RR018898) for CT imaging of the vertebra. The study sponsors had no role in the study design, collection, analysis or interpretation of data.
Compliance with Ethical Standards
Conflict of Interest: The authors declare that they have no conflict of interest.