Hadi S. Hosseiní institute lor Surgical Technology and Biomectianics, university ot Bern, Stauttaciierstr. 78, Bern CH-3014, Switzerland e-maii: [email protected]
Allison L Clouthier Institute tor Surgicai Technoiogy and Biomectianics, iJniversity ol Bern, Stauttactierstr. 78, BernCt1-30t4, Switzeriand
Philippe K. Zysset Institute for Surgicai Technoiogy and Biomechanics, University ot Bern, Stauttacherstr. 78, Bern CH-3014, Switzerland
Experimental Validation of Finite Element Analysis of Human Vertebral Collapse Under Large Compressive Strains Osteoporosis-related vertebral fi-actures represent a major health problem in elderly populations. Such fi-actures can offen only be diagnosed after a substantial deformation history of the vertebral body. Therefore, it remains a challenge for clinicians to distinguish between stable and progressive potentially harmful fractures. Accordingly, novel criteria for selection of the appropriate conservative or surgical treatment are urgently needed. Comptiter tomography-based finite element analysis is an increasingly accepted method to predict the quasi-static vertebral strength and to follow up this .•small strain property longitudinally in time. A recent development in constitutive modeling allows us to simtdate strain localization and densification in trabecular bone under large compressive strains withotit mesh dependence. The aim of this work was to validate this recently developed constitutive model of trabecular bone for the prediction of strain localization and dettsification in the human vertebral body subjected to large compressive deformation. A custom-made stepwise loading device mounted in a high resolution peripheral computer tomogtaphy system was used to describe the progressive collapse of 13 human vertebrae utider axial compression. Continuum firtite element analyses of the 13 compression tests were realized and the zones of high volumetric strain were compared with the experiments. A fair qualitative correspondence of the strain localization zone betv.'een the experiment and finite element analysis was achieved in 9 out of 13 tests and significant correlations of the volumetric strains were obtained throughotu the range of applied axial compression. Interestingly, the stepwise propagating localization zones in trabecular bone converged to the buckling locations in the cortical shell. While the adopted continuum finite element approach still suffers from several limitations, these encouraging preliminary results towardsthe prediction of extended vertebral collapse may help in assessing fracture stability in future work. [DOI: 10.1115/1.4026409] Keywords: vertebral failure, fracture pattern, stepwise loading, densification, large deformation
1
Introduction
Vertebral compression fractures are the most common fractures of the spine [1]. These fractures greatly reduce the quality of life of the elderly and have a profound impact on health care budgets [2,3]. Vertebral compression fractures lead to a reduction in height along its cranio-caudal direction. It results in trabecular bone densification in moderate to severe vertebral deformities since the bone does not usually break into fragments in this case [1]. This type of spine injury occurs mainly as a result of osteoporosis, a skeletal disease which reduces bone mass and increases bone fragility [4]. Osteoporotic vertebral compression fractures range from mild to severe deformities and may cause severe pain, progressively reducing the ability of the patient to perform daily activities [5]. Accordingly, the gradual increase in bone fragility may cause vertebral collapse due to nontraumatic loading conditions. Due to the absence of trauma and the fact that there is no clear onset for the occurrence of vertebral collapse, it may stay silent and undetected [6,7]. The clinical treatment of such fractures demands knowledge of the history of bone failure and severity of deformation. It is important for clinicians to determine Contributed by the Bioengineering Division of ASME for publication in the JOURNAL OF BIOMECHANICAL ENGINEERING. Manuscript received July 9, 2013; final manuscript received December 29, 2013; accepted manuscript posted January 2. 2014; published online March 24, 2014. Assoc. Editor: Pasquale Vena.
Journal of Biomechanical Engineering
first if the fracture is stable or unstable. A stable fracture will not involve critical motions due to physiologic loading and is rarely associated with neurologic hazards [5]. The prediction of deformation patterns in mild to severe vertebral deformities can help advance the accuracy of the diagnosis of such fractures. Therefore, it is important to understand the biomechanical behavior of bone in compression and develop sensitive and reliable diagnostic tools to predict both strength and deformation patterns beyond the ultimate stress of osteoporotic vertebral bodies. High-resolution peripheral quantitative computed tomography (HR-pQCT) and stepwise loading techniques refine our understanding of the local and overall deformation of the bone and provide detailed information about the failure and collapse of the vertebral body during the deformation. This effectively enables the investigation of the mechanical behavior of bone and its failure pattem in compression. The deformation and failure of trabecular bone have been previously investigated by simultaneously coupling the stepwise loading technique and time-lapsed micro computed tomography (,i¡CT) imaging of bone, which is termed image guided failure assessment (IGFA) [8,9]. In a more recent study, the elastic deformation of the endplate under increasing compressive load was measured in a nondestructive manner by applying stepwise loading on functional spine units inside an HR-pQCT machine [10]. Three-dimensional digital image correlation techniques employed in these studies have measured the
Copyright © 2014 by ASiUIE
APRIL2014, Vol. 136 / 041006-1
deformation of bone from the sequence of the ^CT images. These methods, however, have only been used in small deformation analysis of bone [11,12]. It has been reported that the CT-based homogenized finite element (hFE) analysis efficiently enables the prediction of the overall mechanical properties and the localization pattem of the vertebral body in compression [13-17]. This approach is based on the homogenization of the bone micro-architec:ture over local volumetric regions of interest. The bone volume fraction (bone volume over total volume, or BV/TV) accounts for the heterogeneity, while the orientation of the trabeculae (fabric anisotropy) takes into account the anisotropy of the bone micro-architecture. These two parameters are obtained from the CT images and govern the bone's mechanical response [18,19]. Current CT-based hFE models of bone have been mainly validated for the small strain deformations while the prediction of trabecular bone densification and vertebral fracture in large comprí^ssion has been rarely investigated. These hFE models, however, have the potential to realize the large deformation of the bone, as shown in a recent study [20]. Their accuracy in predicting trabecular bone failure and localization in the vertebral body is highly influenced by factors such as the resolution of the CT scans, generation of the hFE models from the CT scans, material ccmstants, applied boundary conditions, exactness of the inner and outer cortical shell geometries, and mesh dependency of the; model due to strain-softening, among others [14-16,21]. The main purpose of the current work was to validate our recent nonlocal damage-plastic model of bone and to experimentally and numerically investigate the deformation pattems seen in mild to severe vertebral deformities. Moreover, we analyzed the trabecular bone localization and densification and the buckling of the cortical shell in vertebral collapse. To this end, a novel experimental setup to monitor the three-dimensional (3D) deformation pattem of the vertebral body up to large compression was developed. The current experimental setup can be used to evaluate the exactness of the FE models in predicting the strength of the vertebral body in small deformation and the trabecular bone densification in moderate to large compression. The specific aims included: (1) designing an experimental setup foi" a better understanding of the trabecular bone failure and collapsie of the cortical shell as a highly coupled structure, (2) modeling the deformation of the vertebral body up to large compression using the CT-based hFE simulations, (3) predicting the localization p.ïttem of the trabecular bone by measuring the volumetric straiti, (4) predicting the stability of vertebral collapse in mild to seven; vertebral compression fractures, and (5) investigating the accuracy of these predictions based on the experiments by quantifying the local evolution of the BV/TV during compression.
2
Methods
2.1 Specimen Preparation. Thirteen cadaveric human vertebral bodies were taken from five thoracolumbar spines, obtained from the Medical University of Vienna, Austria after approval of the ethics commission. The spines were received from three female (ages: 59, 76, and 84) and two male (ages: 75 and 84) donors. None of the patients had cancer reported as a cause of death. The reason for excluding the vertebral specimens affected by pathologies such as cancer in the current stuidy is due to the fact that bone tumors may add more complexity and uncertainty into both the experimental and numerical analyses of vertebral failure. The material properties of lytic or blastic tumors may be quite different than the healthy bone tissue and we; do not yet have proper experimentally validated constants for tumor tissues. In contrast, no study suggests that the tissue properties of osteoporotic bones differ from those of healthy bones [22]. The specimens were stored at —20°C prior to and after preparation. Each vertebra was cleaned from the soft tissues and the postericr elements were sectioned at the pedicles (see Fig. l(a)). Polymethylmethalcrylate 041006-2 / Vol. 136, APRIL 2014
Fig. 1 (a) A vertebral body cleaned from the soft tissues and the posterior elements, and (b) the PMMA end-cap, upon which the shape of the endplate is printed
(PMMA) was used to support the endplates of the vertebra during the mechanical testing. The PMMA provides well-controlled boundary conditions during the compression, which makes it desirable for experiments and numerical simulations. These boundary conditions resemble highly degenerated intervertebral disks where the central pressure of the nucleus pulposus is lost, which is often the case in the elderly [23]. The usual approach to bind the vertebra to the PMMA is to embed the specimen with parallel lOmm-thick PMMA disks with 4 mm thickness above their superior and inferior endplates [24]. This procedure, however, is not appropriate for large deformations, since part of the cortical shell close to the endplates might be embedded, resulting in stiff boundaries preventing the realistic deformations of the vertebra under the endplates. Additionally, the exact reproducibility of the FE mesh for the PMMA becomes more challenging. Therefore, it is important to avoid such a problem in large compression. In the current study, we used an altemative method to avoid this issue: the PMMA components were first mixed and the liquid mixture was then cast into a cylindrical mold with a height of 15 mm and an inner diameter of 65 mm in order to produce the cranial and caudal end-caps for the vertebra. Approximately 90 seconds after casting, the vertebra was gently pressed on the surface of the solidifying PMMA to imprint the shape of the endplate (see Fig. ]{b)). Prior to printing. Vaseline was spread on the endplate and the inner surface of the mold in order to ensure easy removal of the vertebra from the end-cap and the end-cap from the mold after solidification. UHU plus schnellfest 2-K-Epoxidkleber was used to glue the imprinted PMMA end-caps to the vertebral endplates. This procedure replicates more realistic degenerated intervertebral disks compared to a vertebral body which is, up to 6 mm, embedded into the PMMA; moreover, the accurate FE mesh for the PMMA can be easily generated by extruding the nodes of the endplate.
Thrust Bearing Relief Valve Assembly Specimen \ Load Screw
Loading Platens
Attachment to CT
Fig. 2 Compression device used to apply the axial load to the vertebral body. The device is mounted on the XtremeCT machine after each loading step. Transactions of the ASME
2.2 Stépwise Loading and HR-pQCT Imaging. The experimental setup consists of compression tests and HR-pQCT imaging in a coupled fashion. Two compression devices (see Fig. 2) similar to a previous work [10] were designed and manufactured at the Institute for Surgical Technology and Biomechanics (ISTB), Faculty of Medicine, University of Bern, Switzerland. The outer cylinder of the compression device was made of polyoxymethylene-copolymer, which is stiff enough to bear, the tensile loads during the compression and is radiolucent to avoid imaging artifacts. A miniature load cell (Honeywell model 31; capacity, 3000 lb) was integrated into each compression device to allow for continuous measurement of the specimen's reaction force during the loading and relaxation period. The load cell was isolated using a sealing membrane since the specimen was immersed in 0.9% saline solution during the test. The specimens were placed into the loading chamber filled with the saline solution for 30 min to ensure full hydration prior to testing and the wet conditions were maintained during the whole test. The compressive loading was applied via a loading screw (1 mm pitch) and the approximate axial translation of the screw in each loading step was controlled by measuring the rotation angle of the screw using 24 perforated index holes on the upper plate of the compression device next to the head of the loading screw. The exact axial displacement undergone by specimen was then measured using the registered HR-pQCT images with a voxel size of 82 ^m by counting the change in the number of the voxels between the inferior and superior endplates and multiplying by the voxel size. To avoid imaging artifacts as a result of stress-relaxation, which leads to the local infinitesimal deformations of the bone stmcture, the loading device was mounted into the /iCT scanner approximately 25 min after loading. During the relaxation period of one loading device, the other one was placed into the ^iCT scanner for imaging. The approach of continuous imaging and relaxation using two loading devices enabled the testing of two specimens per testing day, up to large deformations (14 steps per specimen). The relaxed force at the end of each loading step was then selected as the corresponding force to the measured displacement of the step. 2.3 FE Modeling. A nonlocal implicit gradient-enhanced damage-plastic model for trabecular bone was recently developed and implemented into the FE software ABAQUS using a user element subroutine [25]. The basic theory of the constitutive model used in this work was explained in detail in Ref. [16]. To further improve the plastic yield criterion, a generalized quadric yield surface for bone recently introduced in Ref. [26] has been adopted in the current model. Moreover, to circumvent the mesh sensitivity of the model to strain-softening phenomenon, a nonlocal gradient approach similar to that used in Refs. [27,28] was used. The model was tested against the experimental stepwise loading results of 16 human vertebral trabecular bone biopsies in a previous work to study the mesh sensitivity of the constitutive law to strain-softening. Each biopsy was meshed using three different mesh sizes and was simulated using both local and nonlocal damage-plastic models. The force-displacement curves and the damage evolution were compared between the different mesh sizes for each sample simulated using both local and nonlocal damage-plastic models. While the force-displacement curves and the damage pattems simulated by the local model varied between the different mesh sizes due to strain-softening, the nonlocal simulations showed no mesh dependency [25], In the current study, the nonlocal damage-plastic model was used to simulate the mechanical behavior of^ human vertebral bodies in large compression. The vertebral bodies were meshed by 10-node quadratic tetrahedron elements with four integration points per element. The characteristic element length varied from 8 to 10 mm for different specimens and between the cortex and the trabecular core of each specimen. The number of nodes varied from 5441 to 13,214 depending on the size of the specimen. The algorithm used to generate the hFE models based on the Journal of Biomechanical Engineering
BVAV
Fig. 3 (a) The /iCT reconstruction of a vertebral body, and (b) the hFE mesh showing the distribution of the homogenized BV/TV for each element
high-resolution CT images is explained in Ref. [29]. The image processing and the generation of the hFE meshes were performed accordingly, using an in-house software inheriting the same algorithm. The BV/TV and fabric anisotropy were mapped from the ßCT images to the FE mesh. To account for the variation of the BV/TV within the element, the BV/TV values were separately measured at each integration point. The 10-node quadratic tetrahedron elements carry four integration points; therefore, assigning a distinct BV/TV to each integration point accounts for the small variation of the bone heterogeneity within the element. This may result in a more sensitive update of the mechanical response since the gradient of the BV/TV within a continuum element may be high, especially in osteoporotic trabecular bones. This improvement also enables the use of coarser continuum elements and stays as accurate by accounting for bone density variation within the element compared to the case that the BV/TV is only acquired for the center of gravity of the element. The BV/TV values were calculated based on the bone mineral density (BMD) values taken from the image data using an empirical equation relating the two values [30] (see Fig. 3). The fabric anisotropy was calculated using the mean intercept length method separately for each element. Since the variation of the fabric within the element is less critical compared to that of BV/TV, the fabric was calculated at the center of gravity of each continuum element [31]. Therefore, all four integration points of every element possess the same fabric tensor. The FE meshes for the PMMA end-caps were generated by extruding the nodes located on the surface of the endplate in the cranio-caudal direction. The bottom surface of the caudal end-cap was fixed in all directions and the top surface of the cranial one was subjected to displacement-controlled uniaxial compression to replicate the experimental conditions. The simulations were performed for the 13 specimens using the nonlocal damage-plastic model. The deformed FE meshes corresponding to the same displacement levels of the experimental loading steps were extracted from the output database of each simulation. The deformed FE meshes were then used for the numerical analysis of the vertebral collapse. The material constants used in the nonlocal model are given in Tables 1 to 5. Here, SQ, (IQ, CTQ, a^. To, 7,0- and y^o are in MPa, the nonlocal length scale /„i is in mm, and the rest of the material constants are dimensionless. All of the material constants denoted by a '0' subscript identify the material properties of the bone tissue, i.e., when BV/TV = 1. Table 1 summarizes the orthotropic elasticity model parameters. Here, ËO, HQ, and uo are the elastic modulus, shear modulus, and Poisson's ratio of the bone tissue and k and /
Table 1 Co (MPa)
1189.6
Fabric-based elasticity material constants Ho (MPa) 249.2
i/o 0.181
k
/
0.972
0.820
APRIL2014, Vol. 136 / 041006-3
Table 2 criterion (7ñ (MPa)
Material constants of the fabric-based
Table 4
strength
(T+ (MPa)
lo (MPa)
Co
p
q
27.79
16.03
0.1
1.290
0.593
660.0 37.24
iVIateriai parameters of the densification i'l
)>(MPa)
2.928
65.0
Tabie 5 Tabie 3 iVIateriai functions
constants
for
hardening
and
damage
' • / .
2.77
100.0
0.75
-0.1
/„I (mm)
m
0.5
7.0
are the power terms of the BV/TV and fabric anisotropy, respectively [32,33]. The values of EQ and /IQ were taken from Ref. [34]; however, they were scaled by a factor of 0.4 to take into account the stress relaxation during stepwise loading. Ttie material constants of the volume fraction and fabric-based strength criterion are shown in Table 2. Here, CTQ , a'^, and TO are tne uniaxial compressive, uniaxial tensile, and shear strength properties of the bone tissue, respectively. In Table 2, Co is an interaction parameter determining the shape of the yield surface and p and q are the power terms of the BV/TV and fabric anisotropy, respectively [26]. The strength parameters of the recently proposed quadric yield surface were taken from Ref. [26], which are based on the strength data of Refs. [35,36]. Accordingly, a^, a'^, and to were scaled by a factor of 0.7 in order to account for stress relaxation. The least squares method was used to determine the two scale factors of the material parameters to best fit the predicted apparent stiffness and the strength values of the model to the ones experimentally measured. Table 3 displays the constants of the post-yield hardening and damage functions. Here, K^ax is the prescribed cumulated plastic strain at the onset of softening, s determines the rate of hardening, Dc is the critical damage value, and a is the rate of damage. Table 4 shows the density-based constants of the densification whose constitutive equations were developed in a previous study [20]. Here, 7,0 and ;•; are non-negative real values of the linear part and y^o ™d /^ are non-negative real values cif the power term and t is the exponent of the power term of the densification. In Table 4,fi^iis a negative real value accounting for the volumetric threshold of the densification. The nonlocal parameters are given in Table 5. Here, m denotes the nonlocal weight and /„i is the introduced nonlocal length scale. The nonlocal weight m has been chosen to be larger than one to ensure the well-posedness of the nonlocal fonnulation [27,28]. The nonlocal length scale /„i has been measured as the average trabecular bone space of the vertebral bodies. 2.4 Qualitative Analysis of Trabecular Bone Densification. The densification of trabecular bone under progressive compression occurs due to pore closure. Densification fihenomenon may result from the increasing contact forces as the trabecular structure collapses under compression. Therefore, bone densification is triggered when local volumetric changes of the trabecular bone reaches a certain value and may exponentially grow in large compaction of the bone [20]. In continuum mechanics, the change in volume is calculated by the determinant of the deformation gradient tensor det F (1) where dv is the current volume of the continuum and dV is its original volume. In the current large strain simulations, an updated Lagrangian formulation is employed, in which the strain 041006-4 / Vol. 136, APRIL 2014
6.0
Ei
Nonlocal material parameters
1.1 0.02
'
0.12
BV/TV
0.40
0.12
Fig. 4 The BV/TV of the transverse siices caicuiated from the (a) /JCT image, and (b) FE mesh in their intact state
measure is the rate of deformation D [37]. It can be mathematically demonstrated that the trace of the deformation rate is linked to the rate of volumetric change [38] trD = -
(2)
The constitutive equations of the densification, as explained in Ref. [20], employs the incremental integration of Eq. (2) for the change in volume (3)
where In/ has been measured at each integration point and their predicted values have been demonstrated next to the ^¿CT images of the deformed specimens; see Figs. 5-8. To facilitate the detection of the collapsed zones compared to the original vertebra, the CT image of the intact specimen has been provided.
2.5 Quantitative Analysis of Trabecular Bone Densification. The tracking of the trabecular bone collapse between large loading steps using 3D digital image correlation was perceived both as a challenge and an overkill for the planned assessment of densification. The evolution of bone mass during the deformation can be alternatively tracked instead of the displacement field. In fact, the change in the BV/TV during deformation correlates with the change in the total volume (TV) of the region of interest (ROI), assuming that the bone volume of the ROI is preserved. This approach enables the quantitative analysis of bone localization for large deformations in a robust manner. To quantify the trabecular bone densiñcation, each specimen was divided into 10 equally distant transverse slices along its craniocaudal direction (see Fig. 4). For this purpose, first the vertebral body height of the registered ^iCT image was measured based on the axial coordinates of the most cranial and caudal bone voxels appearing on the superior and inferior endplates. The height was then divided by 10, the axial coordinates of each slice were Transactions of the ASME
35« •
—Simulation • Expenment
2800 • z 2100 • S 1400 •
^
•
'
.
•
700 • 0 •
3
3SD0 -
6 9 Dispiacement [mm]
[F76.TI|
12
IS
—Simulation
• Expertment 2800 •
• 52100 • iS 1400 • 700 • 0 •
3
6 9 Displacement [mm]
3SO0 -1
12
simulation Experiment
2800
£2100 • I 1400 • 700 • 0 6 9 Displacement {mm| 3500 •
12
—Simulation
|F76.Tl|
2800 .
A
|2100 •
•
• *
£ 1400 • 700 •
1
0¡ li
$
$
9
12
1
Displacement [mm] 3500 •
—Simulation • Experiment
2S00 -
g 2100 • fi 1400 700 • D •
0.15
-1.0
3
6 9 Displacement [mm]
12
Fig. 5 Three-dimensional representation of a vertebral body collapse up to very large compression. Trabecular bone densification is represented by the volumetric strain Iny. The corresponding applied compressive deformation is distinguished by the red circle on the force-displacement curves. obtained, and the BV/TV of the slice was measured based on its BMD value. A similar procedure ^iß, used to measure the height of the smoothly meshed vertebra by finding the corresponding endplate nodes of the FE mesh. The upper and lower coordinates of each slice in the axial direction were determined and the elements whose centers locate inside the slice were collected. The bone volume (BV) of each slice was then calculated by summing up the BV of the collected elements. Assuming that the BV of Journal of Biomechanical Engineering
the elements is conserved during the compression, the numerical BV/TV of each slice in each loading step may be calculated as (BV/TV)^, = ( B V ) Q / ( T V ) J
(4)
where (BV)Q is the sum of the bone volumes of the slice in the undeformed configuration and (TV)J is the sum of the element volumes of the slice at the loading step n. Having the numerical APRIL2014, Vol. 136 / 041006-5
Intact
step
Fig. 6 Qualitative analysis of the vertebral body collapse in five specimens. Trabecuiar bone densification is represented by the volumetric strain In/. The corresponding applied compressive deformation is distinguished by the red circle on the forcedisplacement curves.
hFE models. To prevent the possible geometrical errors arising from this issue, the top slice of the superior and the bottom slice of the inferior endplates were excluded from the quantitative analysis. The coefficient of determination (R-) and the concordance (5) correlation coefficient (p,.) [39] were calculated based on eight slices for all 13 specimens. Here, R^ was plotted next to p,. for nine distinct ranges of appawhere (BV/TV)', is the BV/TV of slice s at step n and (BV/TV)Q is the BV/TV of slice s of the intact specimen. Thus, rent compressive strain (up to 68%), as shown in Fig. 9, in order A(BV/TV) may be introduced as a measure for trabecular bone to analyze how well the numerical and the experimental data fit localization and densification in moderate to large compressive the regression line and to measure the agreement between the two strains. Therefore, in order to evaluate the accuracy of our hFE variables during the compression, respectively. simulations in predicting the local pattem of trabecular bone collapse and densification, A(BV/TV) of the transverse slices obtained from the ^CT images and the hFE simulations was 3 Results tracked during the deformation. Due to the fact that the exactness A representative vertebra has been chosen to qualitatively illusof the FE mesh of the endplate geometry is consti^ained by the ele- trate the three-dimensional (3D) progressive collapse of the bone ment size, very fine geometrical details cannot bs captured in our using both the experimental and numerical outcomes (see Fig. 5).
and the experimental BV/TV of all of the loading steps at hand, the change in the volumetric bone density of eacli transverse slice s at loading step n may be measured as
041006-6 / Vol. 136, APRIL 2014
Transactions of the ASIVIE
Intact -?«*%
Step • um
—Sfittutatiöo
IQQO
§00 50O
-
40O
•
0 •
ß
Î
2 3 Í Dl^iaeemem (mm)
5
—Smultliim
z wa • / \
Ï
*
•
300 -
tt (trtm)
0.15
.-1,0
Fig. 7 Qualitative analysis of the vertebral body collapse in four specimens. Trabecular bone densification is represented by the volumetric strain Iny. The corresponding applied compressive deformation is distinguished by the red circle on the forcedisplacement curves.
Part of the vertebral body has been removed from the image to demonstrate the trabecular core failure on the sagittal and coronal planes. Trabecular bone densification can be distinguished by the increasing gray regions of the collapsed trabeculae. To facilitate the distinction between the localized region(s) from the rest of the specimen, the undeformed state of the specimen has been provided in the corresponding figures. The force-displacement curves include the numerical response (continuous curve) and the experimental data (solid points). The corresponding displacement level of the deformed specimen at each loading step is shown by a red solid circle for clarity. The reaction force of each step has been taken from the recorded force-time data during loading referred to the end of the stressrelaxation period of the step. To quantify the deviation of the measured force from the numerical predictions, the root mean square error (RMSE) has been calculated for the force data and has been normalized with respect to the ultimate force values of the numerical data. Due to the stepwise loading nature of the experiments, the accurate experimental ultimate force values were not obtained with certainty; therefore, the normalization of the error has been performed with respect to the numerical values. The normalized RMSE ranged from 8% to 62%.
Journal of Biomechanicai Engineering
The deformed slices of the 13 specimens are illustrated in Figs. 6-8, in which In / has been chosen as a qualitative criterion in order to demonstrate the localization and the densification of the trabecular bone. Similarly, the red solid circles on the forcedisplacement curves correspond to the deformation level of each collapsing vertebra in these figures. It was observed that the progressive collapse of the vertebral body during compression is characterized by the defiection of the endplate, buckling of the cortical thickness, and the densification of the trabecular core; see Figs. 5-8. The stepwise loading results showed that the hFE simulations correctly predicted the localization pattem of the trabecular bone and the buckling of the cortical thickness in several specimens; see Figs. 6 and 7. However, the simulations failed in the accurate prediction of the vertebral body deformation in a few samples, separately shown in Fig. 8. The statistics (see Fig. 9) showed that the model starts to predict the A(BV/TV) in apparent compressive strains; larger than 4% (end of strain-softening). The accuracy of the predictions progressively improves in higher compressions (up to 60% apparent compressive strain). The systematic analysis of A(BV/TV) as a measure for bone densification showed that this quantity does not provide enough APRIL2014, Vol. 136 / 041006-7
Intact
Step • 4^^ -
MM-Tl
—Hmtiîatiw.
AWQ • 3SKI _K)Ô6 S2SOÜ • |2£ffiO1503 '
iwa •
y «
ísa • í
0.15 4
Allison L Clouthier Institute tor Surgicai Technoiogy and Biomectianics, iJniversity ol Bern, Stauttactierstr. 78, BernCt1-30t4, Switzeriand
Philippe K. Zysset Institute for Surgicai Technoiogy and Biomechanics, University ot Bern, Stauttacherstr. 78, Bern CH-3014, Switzerland
Experimental Validation of Finite Element Analysis of Human Vertebral Collapse Under Large Compressive Strains Osteoporosis-related vertebral fi-actures represent a major health problem in elderly populations. Such fi-actures can offen only be diagnosed after a substantial deformation history of the vertebral body. Therefore, it remains a challenge for clinicians to distinguish between stable and progressive potentially harmful fractures. Accordingly, novel criteria for selection of the appropriate conservative or surgical treatment are urgently needed. Comptiter tomography-based finite element analysis is an increasingly accepted method to predict the quasi-static vertebral strength and to follow up this .•small strain property longitudinally in time. A recent development in constitutive modeling allows us to simtdate strain localization and densification in trabecular bone under large compressive strains withotit mesh dependence. The aim of this work was to validate this recently developed constitutive model of trabecular bone for the prediction of strain localization and dettsification in the human vertebral body subjected to large compressive deformation. A custom-made stepwise loading device mounted in a high resolution peripheral computer tomogtaphy system was used to describe the progressive collapse of 13 human vertebrae utider axial compression. Continuum firtite element analyses of the 13 compression tests were realized and the zones of high volumetric strain were compared with the experiments. A fair qualitative correspondence of the strain localization zone betv.'een the experiment and finite element analysis was achieved in 9 out of 13 tests and significant correlations of the volumetric strains were obtained throughotu the range of applied axial compression. Interestingly, the stepwise propagating localization zones in trabecular bone converged to the buckling locations in the cortical shell. While the adopted continuum finite element approach still suffers from several limitations, these encouraging preliminary results towardsthe prediction of extended vertebral collapse may help in assessing fracture stability in future work. [DOI: 10.1115/1.4026409] Keywords: vertebral failure, fracture pattern, stepwise loading, densification, large deformation
1
Introduction
Vertebral compression fractures are the most common fractures of the spine [1]. These fractures greatly reduce the quality of life of the elderly and have a profound impact on health care budgets [2,3]. Vertebral compression fractures lead to a reduction in height along its cranio-caudal direction. It results in trabecular bone densification in moderate to severe vertebral deformities since the bone does not usually break into fragments in this case [1]. This type of spine injury occurs mainly as a result of osteoporosis, a skeletal disease which reduces bone mass and increases bone fragility [4]. Osteoporotic vertebral compression fractures range from mild to severe deformities and may cause severe pain, progressively reducing the ability of the patient to perform daily activities [5]. Accordingly, the gradual increase in bone fragility may cause vertebral collapse due to nontraumatic loading conditions. Due to the absence of trauma and the fact that there is no clear onset for the occurrence of vertebral collapse, it may stay silent and undetected [6,7]. The clinical treatment of such fractures demands knowledge of the history of bone failure and severity of deformation. It is important for clinicians to determine Contributed by the Bioengineering Division of ASME for publication in the JOURNAL OF BIOMECHANICAL ENGINEERING. Manuscript received July 9, 2013; final manuscript received December 29, 2013; accepted manuscript posted January 2. 2014; published online March 24, 2014. Assoc. Editor: Pasquale Vena.
Journal of Biomechanical Engineering
first if the fracture is stable or unstable. A stable fracture will not involve critical motions due to physiologic loading and is rarely associated with neurologic hazards [5]. The prediction of deformation patterns in mild to severe vertebral deformities can help advance the accuracy of the diagnosis of such fractures. Therefore, it is important to understand the biomechanical behavior of bone in compression and develop sensitive and reliable diagnostic tools to predict both strength and deformation patterns beyond the ultimate stress of osteoporotic vertebral bodies. High-resolution peripheral quantitative computed tomography (HR-pQCT) and stepwise loading techniques refine our understanding of the local and overall deformation of the bone and provide detailed information about the failure and collapse of the vertebral body during the deformation. This effectively enables the investigation of the mechanical behavior of bone and its failure pattem in compression. The deformation and failure of trabecular bone have been previously investigated by simultaneously coupling the stepwise loading technique and time-lapsed micro computed tomography (,i¡CT) imaging of bone, which is termed image guided failure assessment (IGFA) [8,9]. In a more recent study, the elastic deformation of the endplate under increasing compressive load was measured in a nondestructive manner by applying stepwise loading on functional spine units inside an HR-pQCT machine [10]. Three-dimensional digital image correlation techniques employed in these studies have measured the
Copyright © 2014 by ASiUIE
APRIL2014, Vol. 136 / 041006-1
deformation of bone from the sequence of the ^CT images. These methods, however, have only been used in small deformation analysis of bone [11,12]. It has been reported that the CT-based homogenized finite element (hFE) analysis efficiently enables the prediction of the overall mechanical properties and the localization pattem of the vertebral body in compression [13-17]. This approach is based on the homogenization of the bone micro-architec:ture over local volumetric regions of interest. The bone volume fraction (bone volume over total volume, or BV/TV) accounts for the heterogeneity, while the orientation of the trabeculae (fabric anisotropy) takes into account the anisotropy of the bone micro-architecture. These two parameters are obtained from the CT images and govern the bone's mechanical response [18,19]. Current CT-based hFE models of bone have been mainly validated for the small strain deformations while the prediction of trabecular bone densification and vertebral fracture in large comprí^ssion has been rarely investigated. These hFE models, however, have the potential to realize the large deformation of the bone, as shown in a recent study [20]. Their accuracy in predicting trabecular bone failure and localization in the vertebral body is highly influenced by factors such as the resolution of the CT scans, generation of the hFE models from the CT scans, material ccmstants, applied boundary conditions, exactness of the inner and outer cortical shell geometries, and mesh dependency of the; model due to strain-softening, among others [14-16,21]. The main purpose of the current work was to validate our recent nonlocal damage-plastic model of bone and to experimentally and numerically investigate the deformation pattems seen in mild to severe vertebral deformities. Moreover, we analyzed the trabecular bone localization and densification and the buckling of the cortical shell in vertebral collapse. To this end, a novel experimental setup to monitor the three-dimensional (3D) deformation pattem of the vertebral body up to large compression was developed. The current experimental setup can be used to evaluate the exactness of the FE models in predicting the strength of the vertebral body in small deformation and the trabecular bone densification in moderate to large compression. The specific aims included: (1) designing an experimental setup foi" a better understanding of the trabecular bone failure and collapsie of the cortical shell as a highly coupled structure, (2) modeling the deformation of the vertebral body up to large compression using the CT-based hFE simulations, (3) predicting the localization p.ïttem of the trabecular bone by measuring the volumetric straiti, (4) predicting the stability of vertebral collapse in mild to seven; vertebral compression fractures, and (5) investigating the accuracy of these predictions based on the experiments by quantifying the local evolution of the BV/TV during compression.
2
Methods
2.1 Specimen Preparation. Thirteen cadaveric human vertebral bodies were taken from five thoracolumbar spines, obtained from the Medical University of Vienna, Austria after approval of the ethics commission. The spines were received from three female (ages: 59, 76, and 84) and two male (ages: 75 and 84) donors. None of the patients had cancer reported as a cause of death. The reason for excluding the vertebral specimens affected by pathologies such as cancer in the current stuidy is due to the fact that bone tumors may add more complexity and uncertainty into both the experimental and numerical analyses of vertebral failure. The material properties of lytic or blastic tumors may be quite different than the healthy bone tissue and we; do not yet have proper experimentally validated constants for tumor tissues. In contrast, no study suggests that the tissue properties of osteoporotic bones differ from those of healthy bones [22]. The specimens were stored at —20°C prior to and after preparation. Each vertebra was cleaned from the soft tissues and the postericr elements were sectioned at the pedicles (see Fig. l(a)). Polymethylmethalcrylate 041006-2 / Vol. 136, APRIL 2014
Fig. 1 (a) A vertebral body cleaned from the soft tissues and the posterior elements, and (b) the PMMA end-cap, upon which the shape of the endplate is printed
(PMMA) was used to support the endplates of the vertebra during the mechanical testing. The PMMA provides well-controlled boundary conditions during the compression, which makes it desirable for experiments and numerical simulations. These boundary conditions resemble highly degenerated intervertebral disks where the central pressure of the nucleus pulposus is lost, which is often the case in the elderly [23]. The usual approach to bind the vertebra to the PMMA is to embed the specimen with parallel lOmm-thick PMMA disks with 4 mm thickness above their superior and inferior endplates [24]. This procedure, however, is not appropriate for large deformations, since part of the cortical shell close to the endplates might be embedded, resulting in stiff boundaries preventing the realistic deformations of the vertebra under the endplates. Additionally, the exact reproducibility of the FE mesh for the PMMA becomes more challenging. Therefore, it is important to avoid such a problem in large compression. In the current study, we used an altemative method to avoid this issue: the PMMA components were first mixed and the liquid mixture was then cast into a cylindrical mold with a height of 15 mm and an inner diameter of 65 mm in order to produce the cranial and caudal end-caps for the vertebra. Approximately 90 seconds after casting, the vertebra was gently pressed on the surface of the solidifying PMMA to imprint the shape of the endplate (see Fig. ]{b)). Prior to printing. Vaseline was spread on the endplate and the inner surface of the mold in order to ensure easy removal of the vertebra from the end-cap and the end-cap from the mold after solidification. UHU plus schnellfest 2-K-Epoxidkleber was used to glue the imprinted PMMA end-caps to the vertebral endplates. This procedure replicates more realistic degenerated intervertebral disks compared to a vertebral body which is, up to 6 mm, embedded into the PMMA; moreover, the accurate FE mesh for the PMMA can be easily generated by extruding the nodes of the endplate.
Thrust Bearing Relief Valve Assembly Specimen \ Load Screw
Loading Platens
Attachment to CT
Fig. 2 Compression device used to apply the axial load to the vertebral body. The device is mounted on the XtremeCT machine after each loading step. Transactions of the ASME
2.2 Stépwise Loading and HR-pQCT Imaging. The experimental setup consists of compression tests and HR-pQCT imaging in a coupled fashion. Two compression devices (see Fig. 2) similar to a previous work [10] were designed and manufactured at the Institute for Surgical Technology and Biomechanics (ISTB), Faculty of Medicine, University of Bern, Switzerland. The outer cylinder of the compression device was made of polyoxymethylene-copolymer, which is stiff enough to bear, the tensile loads during the compression and is radiolucent to avoid imaging artifacts. A miniature load cell (Honeywell model 31; capacity, 3000 lb) was integrated into each compression device to allow for continuous measurement of the specimen's reaction force during the loading and relaxation period. The load cell was isolated using a sealing membrane since the specimen was immersed in 0.9% saline solution during the test. The specimens were placed into the loading chamber filled with the saline solution for 30 min to ensure full hydration prior to testing and the wet conditions were maintained during the whole test. The compressive loading was applied via a loading screw (1 mm pitch) and the approximate axial translation of the screw in each loading step was controlled by measuring the rotation angle of the screw using 24 perforated index holes on the upper plate of the compression device next to the head of the loading screw. The exact axial displacement undergone by specimen was then measured using the registered HR-pQCT images with a voxel size of 82 ^m by counting the change in the number of the voxels between the inferior and superior endplates and multiplying by the voxel size. To avoid imaging artifacts as a result of stress-relaxation, which leads to the local infinitesimal deformations of the bone stmcture, the loading device was mounted into the /iCT scanner approximately 25 min after loading. During the relaxation period of one loading device, the other one was placed into the ^iCT scanner for imaging. The approach of continuous imaging and relaxation using two loading devices enabled the testing of two specimens per testing day, up to large deformations (14 steps per specimen). The relaxed force at the end of each loading step was then selected as the corresponding force to the measured displacement of the step. 2.3 FE Modeling. A nonlocal implicit gradient-enhanced damage-plastic model for trabecular bone was recently developed and implemented into the FE software ABAQUS using a user element subroutine [25]. The basic theory of the constitutive model used in this work was explained in detail in Ref. [16]. To further improve the plastic yield criterion, a generalized quadric yield surface for bone recently introduced in Ref. [26] has been adopted in the current model. Moreover, to circumvent the mesh sensitivity of the model to strain-softening phenomenon, a nonlocal gradient approach similar to that used in Refs. [27,28] was used. The model was tested against the experimental stepwise loading results of 16 human vertebral trabecular bone biopsies in a previous work to study the mesh sensitivity of the constitutive law to strain-softening. Each biopsy was meshed using three different mesh sizes and was simulated using both local and nonlocal damage-plastic models. The force-displacement curves and the damage evolution were compared between the different mesh sizes for each sample simulated using both local and nonlocal damage-plastic models. While the force-displacement curves and the damage pattems simulated by the local model varied between the different mesh sizes due to strain-softening, the nonlocal simulations showed no mesh dependency [25], In the current study, the nonlocal damage-plastic model was used to simulate the mechanical behavior of^ human vertebral bodies in large compression. The vertebral bodies were meshed by 10-node quadratic tetrahedron elements with four integration points per element. The characteristic element length varied from 8 to 10 mm for different specimens and between the cortex and the trabecular core of each specimen. The number of nodes varied from 5441 to 13,214 depending on the size of the specimen. The algorithm used to generate the hFE models based on the Journal of Biomechanical Engineering
BVAV
Fig. 3 (a) The /iCT reconstruction of a vertebral body, and (b) the hFE mesh showing the distribution of the homogenized BV/TV for each element
high-resolution CT images is explained in Ref. [29]. The image processing and the generation of the hFE meshes were performed accordingly, using an in-house software inheriting the same algorithm. The BV/TV and fabric anisotropy were mapped from the ßCT images to the FE mesh. To account for the variation of the BV/TV within the element, the BV/TV values were separately measured at each integration point. The 10-node quadratic tetrahedron elements carry four integration points; therefore, assigning a distinct BV/TV to each integration point accounts for the small variation of the bone heterogeneity within the element. This may result in a more sensitive update of the mechanical response since the gradient of the BV/TV within a continuum element may be high, especially in osteoporotic trabecular bones. This improvement also enables the use of coarser continuum elements and stays as accurate by accounting for bone density variation within the element compared to the case that the BV/TV is only acquired for the center of gravity of the element. The BV/TV values were calculated based on the bone mineral density (BMD) values taken from the image data using an empirical equation relating the two values [30] (see Fig. 3). The fabric anisotropy was calculated using the mean intercept length method separately for each element. Since the variation of the fabric within the element is less critical compared to that of BV/TV, the fabric was calculated at the center of gravity of each continuum element [31]. Therefore, all four integration points of every element possess the same fabric tensor. The FE meshes for the PMMA end-caps were generated by extruding the nodes located on the surface of the endplate in the cranio-caudal direction. The bottom surface of the caudal end-cap was fixed in all directions and the top surface of the cranial one was subjected to displacement-controlled uniaxial compression to replicate the experimental conditions. The simulations were performed for the 13 specimens using the nonlocal damage-plastic model. The deformed FE meshes corresponding to the same displacement levels of the experimental loading steps were extracted from the output database of each simulation. The deformed FE meshes were then used for the numerical analysis of the vertebral collapse. The material constants used in the nonlocal model are given in Tables 1 to 5. Here, SQ, (IQ, CTQ, a^. To, 7,0- and y^o are in MPa, the nonlocal length scale /„i is in mm, and the rest of the material constants are dimensionless. All of the material constants denoted by a '0' subscript identify the material properties of the bone tissue, i.e., when BV/TV = 1. Table 1 summarizes the orthotropic elasticity model parameters. Here, ËO, HQ, and uo are the elastic modulus, shear modulus, and Poisson's ratio of the bone tissue and k and /
Table 1 Co (MPa)
1189.6
Fabric-based elasticity material constants Ho (MPa) 249.2
i/o 0.181
k
/
0.972
0.820
APRIL2014, Vol. 136 / 041006-3
Table 2 criterion (7ñ (MPa)
Material constants of the fabric-based
Table 4
strength
(T+ (MPa)
lo (MPa)
Co
p
q
27.79
16.03
0.1
1.290
0.593
660.0 37.24
iVIateriai parameters of the densification i'l
)>(MPa)
2.928
65.0
Tabie 5 Tabie 3 iVIateriai functions
constants
for
hardening
and
damage
' • / .
2.77
100.0
0.75
-0.1
/„I (mm)
m
0.5
7.0
are the power terms of the BV/TV and fabric anisotropy, respectively [32,33]. The values of EQ and /IQ were taken from Ref. [34]; however, they were scaled by a factor of 0.4 to take into account the stress relaxation during stepwise loading. Ttie material constants of the volume fraction and fabric-based strength criterion are shown in Table 2. Here, CTQ , a'^, and TO are tne uniaxial compressive, uniaxial tensile, and shear strength properties of the bone tissue, respectively. In Table 2, Co is an interaction parameter determining the shape of the yield surface and p and q are the power terms of the BV/TV and fabric anisotropy, respectively [26]. The strength parameters of the recently proposed quadric yield surface were taken from Ref. [26], which are based on the strength data of Refs. [35,36]. Accordingly, a^, a'^, and to were scaled by a factor of 0.7 in order to account for stress relaxation. The least squares method was used to determine the two scale factors of the material parameters to best fit the predicted apparent stiffness and the strength values of the model to the ones experimentally measured. Table 3 displays the constants of the post-yield hardening and damage functions. Here, K^ax is the prescribed cumulated plastic strain at the onset of softening, s determines the rate of hardening, Dc is the critical damage value, and a is the rate of damage. Table 4 shows the density-based constants of the densification whose constitutive equations were developed in a previous study [20]. Here, 7,0 and ;•; are non-negative real values of the linear part and y^o ™d /^ are non-negative real values cif the power term and t is the exponent of the power term of the densification. In Table 4,fi^iis a negative real value accounting for the volumetric threshold of the densification. The nonlocal parameters are given in Table 5. Here, m denotes the nonlocal weight and /„i is the introduced nonlocal length scale. The nonlocal weight m has been chosen to be larger than one to ensure the well-posedness of the nonlocal fonnulation [27,28]. The nonlocal length scale /„i has been measured as the average trabecular bone space of the vertebral bodies. 2.4 Qualitative Analysis of Trabecular Bone Densification. The densification of trabecular bone under progressive compression occurs due to pore closure. Densification fihenomenon may result from the increasing contact forces as the trabecular structure collapses under compression. Therefore, bone densification is triggered when local volumetric changes of the trabecular bone reaches a certain value and may exponentially grow in large compaction of the bone [20]. In continuum mechanics, the change in volume is calculated by the determinant of the deformation gradient tensor det F (1) where dv is the current volume of the continuum and dV is its original volume. In the current large strain simulations, an updated Lagrangian formulation is employed, in which the strain 041006-4 / Vol. 136, APRIL 2014
6.0
Ei
Nonlocal material parameters
1.1 0.02
'
0.12
BV/TV
0.40
0.12
Fig. 4 The BV/TV of the transverse siices caicuiated from the (a) /JCT image, and (b) FE mesh in their intact state
measure is the rate of deformation D [37]. It can be mathematically demonstrated that the trace of the deformation rate is linked to the rate of volumetric change [38] trD = -
(2)
The constitutive equations of the densification, as explained in Ref. [20], employs the incremental integration of Eq. (2) for the change in volume (3)
where In/ has been measured at each integration point and their predicted values have been demonstrated next to the ^¿CT images of the deformed specimens; see Figs. 5-8. To facilitate the detection of the collapsed zones compared to the original vertebra, the CT image of the intact specimen has been provided.
2.5 Quantitative Analysis of Trabecular Bone Densification. The tracking of the trabecular bone collapse between large loading steps using 3D digital image correlation was perceived both as a challenge and an overkill for the planned assessment of densification. The evolution of bone mass during the deformation can be alternatively tracked instead of the displacement field. In fact, the change in the BV/TV during deformation correlates with the change in the total volume (TV) of the region of interest (ROI), assuming that the bone volume of the ROI is preserved. This approach enables the quantitative analysis of bone localization for large deformations in a robust manner. To quantify the trabecular bone densiñcation, each specimen was divided into 10 equally distant transverse slices along its craniocaudal direction (see Fig. 4). For this purpose, first the vertebral body height of the registered ^iCT image was measured based on the axial coordinates of the most cranial and caudal bone voxels appearing on the superior and inferior endplates. The height was then divided by 10, the axial coordinates of each slice were Transactions of the ASME
35« •
—Simulation • Expenment
2800 • z 2100 • S 1400 •
^
•
'
.
•
700 • 0 •
3
3SD0 -
6 9 Dispiacement [mm]
[F76.TI|
12
IS
—Simulation
• Expertment 2800 •
• 52100 • iS 1400 • 700 • 0 •
3
6 9 Displacement [mm]
3SO0 -1
12
simulation Experiment
2800
£2100 • I 1400 • 700 • 0 6 9 Displacement {mm| 3500 •
12
—Simulation
|F76.Tl|
2800 .
A
|2100 •
•
• *
£ 1400 • 700 •
1
0¡ li
$
$
9
12
1
Displacement [mm] 3500 •
—Simulation • Experiment
2S00 -
g 2100 • fi 1400 700 • D •
0.15
-1.0
3
6 9 Displacement [mm]
12
Fig. 5 Three-dimensional representation of a vertebral body collapse up to very large compression. Trabecular bone densification is represented by the volumetric strain Iny. The corresponding applied compressive deformation is distinguished by the red circle on the force-displacement curves. obtained, and the BV/TV of the slice was measured based on its BMD value. A similar procedure ^iß, used to measure the height of the smoothly meshed vertebra by finding the corresponding endplate nodes of the FE mesh. The upper and lower coordinates of each slice in the axial direction were determined and the elements whose centers locate inside the slice were collected. The bone volume (BV) of each slice was then calculated by summing up the BV of the collected elements. Assuming that the BV of Journal of Biomechanical Engineering
the elements is conserved during the compression, the numerical BV/TV of each slice in each loading step may be calculated as (BV/TV)^, = ( B V ) Q / ( T V ) J
(4)
where (BV)Q is the sum of the bone volumes of the slice in the undeformed configuration and (TV)J is the sum of the element volumes of the slice at the loading step n. Having the numerical APRIL2014, Vol. 136 / 041006-5
Intact
step
Fig. 6 Qualitative analysis of the vertebral body collapse in five specimens. Trabecuiar bone densification is represented by the volumetric strain In/. The corresponding applied compressive deformation is distinguished by the red circle on the forcedisplacement curves.
hFE models. To prevent the possible geometrical errors arising from this issue, the top slice of the superior and the bottom slice of the inferior endplates were excluded from the quantitative analysis. The coefficient of determination (R-) and the concordance (5) correlation coefficient (p,.) [39] were calculated based on eight slices for all 13 specimens. Here, R^ was plotted next to p,. for nine distinct ranges of appawhere (BV/TV)', is the BV/TV of slice s at step n and (BV/TV)Q is the BV/TV of slice s of the intact specimen. Thus, rent compressive strain (up to 68%), as shown in Fig. 9, in order A(BV/TV) may be introduced as a measure for trabecular bone to analyze how well the numerical and the experimental data fit localization and densification in moderate to large compressive the regression line and to measure the agreement between the two strains. Therefore, in order to evaluate the accuracy of our hFE variables during the compression, respectively. simulations in predicting the local pattem of trabecular bone collapse and densification, A(BV/TV) of the transverse slices obtained from the ^CT images and the hFE simulations was 3 Results tracked during the deformation. Due to the fact that the exactness A representative vertebra has been chosen to qualitatively illusof the FE mesh of the endplate geometry is consti^ained by the ele- trate the three-dimensional (3D) progressive collapse of the bone ment size, very fine geometrical details cannot bs captured in our using both the experimental and numerical outcomes (see Fig. 5).
and the experimental BV/TV of all of the loading steps at hand, the change in the volumetric bone density of eacli transverse slice s at loading step n may be measured as
041006-6 / Vol. 136, APRIL 2014
Transactions of the ASIVIE
Intact -?«*%
Step • um
—Sfittutatiöo
IQQO
§00 50O
-
40O
•
0 •
ß
Î
2 3 Í Dl^iaeemem (mm)
5
—Smultliim
z wa • / \
Ï
*
•
300 -
tt (trtm)
0.15
.-1,0
Fig. 7 Qualitative analysis of the vertebral body collapse in four specimens. Trabecular bone densification is represented by the volumetric strain Iny. The corresponding applied compressive deformation is distinguished by the red circle on the forcedisplacement curves.
Part of the vertebral body has been removed from the image to demonstrate the trabecular core failure on the sagittal and coronal planes. Trabecular bone densification can be distinguished by the increasing gray regions of the collapsed trabeculae. To facilitate the distinction between the localized region(s) from the rest of the specimen, the undeformed state of the specimen has been provided in the corresponding figures. The force-displacement curves include the numerical response (continuous curve) and the experimental data (solid points). The corresponding displacement level of the deformed specimen at each loading step is shown by a red solid circle for clarity. The reaction force of each step has been taken from the recorded force-time data during loading referred to the end of the stressrelaxation period of the step. To quantify the deviation of the measured force from the numerical predictions, the root mean square error (RMSE) has been calculated for the force data and has been normalized with respect to the ultimate force values of the numerical data. Due to the stepwise loading nature of the experiments, the accurate experimental ultimate force values were not obtained with certainty; therefore, the normalization of the error has been performed with respect to the numerical values. The normalized RMSE ranged from 8% to 62%.
Journal of Biomechanicai Engineering
The deformed slices of the 13 specimens are illustrated in Figs. 6-8, in which In / has been chosen as a qualitative criterion in order to demonstrate the localization and the densification of the trabecular bone. Similarly, the red solid circles on the forcedisplacement curves correspond to the deformation level of each collapsing vertebra in these figures. It was observed that the progressive collapse of the vertebral body during compression is characterized by the defiection of the endplate, buckling of the cortical thickness, and the densification of the trabecular core; see Figs. 5-8. The stepwise loading results showed that the hFE simulations correctly predicted the localization pattem of the trabecular bone and the buckling of the cortical thickness in several specimens; see Figs. 6 and 7. However, the simulations failed in the accurate prediction of the vertebral body deformation in a few samples, separately shown in Fig. 8. The statistics (see Fig. 9) showed that the model starts to predict the A(BV/TV) in apparent compressive strains; larger than 4% (end of strain-softening). The accuracy of the predictions progressively improves in higher compressions (up to 60% apparent compressive strain). The systematic analysis of A(BV/TV) as a measure for bone densification showed that this quantity does not provide enough APRIL2014, Vol. 136 / 041006-7
Intact
Step • 4^^ -
MM-Tl
—Hmtiîatiw.
AWQ • 3SKI _K)Ô6 S2SOÜ • |2£ffiO1503 '
iwa •
y «
ísa • í
0.15 4
