MULTI-SCALE MODELING OF THE HUMAN VERTEBRAL BODY: COMPARISON OF MICRO-CT BASED HIGHRESOLUTION AND CONTINUUM-LEVEL MODELS
SENTHIL K. ESWARAN AARON J. FIELDS PREM NAGARATHNAM TONY M. KEAVENY Orthopaedic Biomechanics Laboratory, Department of Mechanical Engineering, University of California, Berkeley, CA 94720, USA

The overall goal of this study was to assess the mechanistic fidelity of continuum-level finite element models of the vertebral body, which represent a promising tool for understanding and predicting clinical fracture risk. Two finite element (FE) models were generated from micro-CT scans of each of 13 T9 vertebral bodies -- a micro-FE model at 60-micron resolution and a coarsened, continuum-level model at 0.96-mm resolution. Two previously-reported continuum-level modulus-density relationships for human vertebral bone were parametrically varied to investigate their effects on model fidelity using the micro-CT models as a gold standard. We found that the modulusdensity relation, particularly that assigned to the peripheral bone, substantially altered the regression coefficients, but not the degree of correlation between continuum and micro-FE predictions of whole-vertebral stiffness. The major load paths through the vertebrae compared well between the continuum-level and micro-FE models (vonMises distribution), but the distributions of minimum principal strain were notably different. We conclude that continuum-level models provide robust measures of whole-vertebral behavior, describe well the load transfer paths through the vertebra, but provide strain distributions that are markedly different than the volume-averaged micro-scale strains. Appreciation of these multi-scale differences should improve interpretation of results from these sorts of continuum models and may improve their clinical utility.

1.

Introduction

Patient-specific voxel-type continuum-level finite element models of the v e r t e b r a l body, that are generated automatically from clinical-resolution quantitative computed tomography (CT) scans, represent a promising alternative t o standard densitometric techniques in understanding and predicting osteoporotic fracture risk [1, 2]. Due to the averaging process during clinical CT imaging, the microstructure of the trabeculae, which has a characteristic dimension on the order of 150 microns, is homogenized into a density-varying continuum of 1 mm-sized image voxels. The thin cortical shell, having a thickness of less than half a millimeter, is also averaged, and indeed there is blurring at the bone edges due to partial volume averaging effects. As a result, it is unclear if there remains a correspondence in the ability of these continuum-


level models to properly capture the underlying physics of the bone failure. Understanding the relation between the behavior of these continuum-level m o d e l s and more realistic micro-mechanical models that capture the microstructure is important since it may lead to identification of means to improve the fidelity of the continuum level models and improve interpretation of the outcomes of such models. While a previous study that performed multi-scale comparison of the continuum-level and micro-FE models for the human femora provided much insight [3], differences in the specific modeling details as well as the mechanics between the femur and the vertebral body confounds the extrapolation of their results to the vertebral body. Due to difficulties in explicitly characterizing the continuum behavior of the cortical shell and associated limitations in using a generic model of the cortical shell [4-6], continuum-level models typically ignore explicit modeling of the cortical shell [1, 2, 7] and apply experimental mechanical property-density relation -- developed for vertebral trabecular bone but not the cortical shell and adjacent bone -- to all bone elements. The implications of using different material property-density relationships are unclear, and in particular, the validity of extrapolating these relations beyond their reported density range. Silva et al. [8] found that strain contours from continuum-level models were better indicators of experimental fracture locations than stress contours, but that study was limited to vertebral sections in which load transfer patterns through the vertebral body may have been over-simplified. Thus, there remains an incomplete understanding of ability of continuum-level models to predict stress or strain distributions within the vertebral body. The overall goal of this study was to assess the mechanistic fidelity of continuum-level models of the whole vertebra -- and some of the assumptions used in creating those models -- by using a multi-scale approach to compare its predictions to that of more mechanistic micro-FE models. Our first objective was to determine the effect of different modulus-density relations on the agreement between the continuum-level and micro-FE model predictions of vertebral stiffness, and our second objective was to compare the stress and strain distributions of the multi-scale models. This study extends the use of the multiscale modeling approach, previously applied to human femora [3], to a cohort of elderly whole vertebrae and is unique in its focus on the micromechanics of the multi-scale models. 2. Methods

After removal of posterior elements, thirteen T9 cadaveric vertebral bodies (88 ± 4 years, n=3 male, n=10 female) were micro-CT scanned at 30-mm voxel size


(SCANCO Medical AG, Brüttisellen, Switzerland). Segmented images were obtained using a global threshold value. Two finite element models were generated for each vertebra ­ a micro-FE model with a 60-micron element size and a coarsened, continuum-level model at 0.96 mm resolution. The micro-FE model was compartmentalized into trabecular bone, cortical shell and cortical endplates using previously validated algorithms [9, 10]. All bone tissue in the micro-FE models was assigned a Young's modulus of 18.5 GPa and Poisson's ratio of 0.3 [11]. A tissue density of 2.05 g/cm3 [12] was used to convert the volume fraction of each continuum element into apparent density. Two different modulus-apparent density relations reported in the literature for human vertebral trabecular bone [12, 13] were modified to generate a total of four parametric variations (Table 1, Fig. 1). To study the effect of extrapolating the modulus-density relations beyond the experimental data, the two reported relations [12, 13] (Models A, B) were modified within or beyond the experimental density range (Models C, D). The number of continuum-level elements (expressed as a percentage of the total number of elements in the model) beyond the experimental density range was calculated for each vertebral body. Transversely isotropic properties were assumed for all continuum-level elements [2]. A layer of polymethylmethacrylate (PMMA) was added to the superior and inferior endplates in order to provide uniform loading conditions across all vertebrae. The resulting micro-FE models contained up to 80 million elements. Using custom code [14] and a maximum of 640 processors in parallel, a uniform compressive strain of 1% was applied in the superiorinferior direction to each model.
Table 1. Different modulus-density relations analyzed in this study. Model A B C D Modulus (MPa) ­ Apparent Density (g/cm3 ) Mapping

E zz = 2130r - 97.1 [13] E zz = 4730r1.56 [12] Model B for r £ 0.38 , Model A for r > 0.38 M  odel A for r £ 0.38 , Model B for r > 0.38

  

 


7

Percent of Total Elements

6 5 4 3 2 1 0 0.0

12 10 8

Experimental range

6 4 2

0.5

1.0

1.5

2.0

0

Apparent Density (g/cc)
Figure 1. Modulus-density relations (Models A, B) for vertebral trabecular bone and a typical density distribution extrapolated to the entire density range. See Table 1 for equations for Models A and B.

Whole vertebral stiffnesses of the continuum-models (Models A­D) were correlated with that of the micro-FE model. The difference in the stiffness predictions between models A and B which use the modulus-density relations reported in the literature was correlated with the micro-FE stiffness in order to statistically test whether the choice of the modulus-density relation affected the slope and/or intercept of the correlation between continuum-level and micro-FE stiffness predictions. The stiffness predictions of models C and D were compared in a similar way to test whether differences in the modulus-density relations within (Model C vs. A) or beyond (Model D vs. A) the experimental density range had a greater effect on the regression coefficients. The minimum principal strain and von-Mises stress distributions were compared between the multi-scale models. The results from the micro-FE models were averaged over each continuum-level element to facilitate this comparison (Fig. 2). Moreover, the distribution within each finite element model was compartmentalized into three groups: 0­50th percentile, 51­75th p e r c e n t i l e , and 76th­100th percentile. Quantitative comparison of the correspondence between the continuum-level elements identified as being at high stress/strain (76th­100th percentile) and the averaged results from the micro-FE models was performed. Preliminary investigation showed that the strain distribution of the continuum-level models may largely be governed by the axial stiffness of each transverse layer, i.e., SEi Ai summed over i = 1, n continuumlevel elements in that transverse slice. For each vertebral body and for both the

Apparent Modulus (GPa)

Model A Model B Density Distribution

14


continuum-level and micro-FE models, the number of elements belonging to the high strain category (76th-100th percentile) in each transverse layer was plotted against the axial stiffness of that layer. The transverse layers that included the endplates were not included as part of this regression since the curvature of the endplates confounds measurement of axial stiffness.

Figure 2. Flowchart describing the procedure by which regions of high-stress were compared between micro-FE and continuum-level models. A similar approach was used for minimum principal strain.

3.

Results

Stiffness predictions of all parametric variations of the continuum-level model were highly correlated with micro-FE stiffness (Fig. 3, R 2=0.89­0.93). Choice of the modulus-density relation had a substantial, statistically significant effect on the regression coefficients between the predicted continuum-level FE stiffness and micro-FE stiffness but did not appreciably affect the R2 values. Despite the majority of elements being within the experimental density range (73.2±5.8%), differences in the modulus-density relations beyond the reported density range (Model D vs. A) had a greater effect (p<0.0001) on the regression with the micro-FE stiffness than did any differences associated with changes to the modulus-density relation for bone within the experimental density range (Model C vs. A). The elements with density values beyond the experimental density range were typically at the periphery of the vertebra and included the cortical shell and endplates. While the major load paths compared well between the continuum-level and micro-FE models (as seen from the von-Mises distributions), the strain distribution in the continuum-level model was notably different from that of the micro-FE model (Fig. 4). Comparison of regions with high von-Mises stress


showed that the correspondence between the continuum-level and micro-FE models was similar in the trabecular centrum and cortical shell (Fig. 5).
Model A: y = 3.50 + 0.35x Model B: y = 1.55 + 0.63x Model C: y = 5.04 + 0.36x Model D: y = -0.08 + 0.61x
2 R = 0.90 R2= 0.92 R2= 0.89 R2= 0.93

40

Continuum-level FE Stiffness (kN/mm)

35 30 25 20 15 10 15 20

25

30

35

40

45

50

Micro-FE Stiffness (kN/mm)
Figure 3. Stiffness predictions of all parametric variations of the continuum-level model were highly correlated with micro-FE stiffness. Material property changes to the bone beyond the experimental density range had a greater effect on the regression coefficients than changes to the bone within the experimental range. See Table 1 for descriptions for Models A-D.

Figure 4. Typical distribution of von-Mises stress and minimum principal strain at a mid-coronal slice (0.96 mm thick) of a vertebral body. Middle and right-hand panels show volume-averaged measures. While the continuum-level models captured the major load paths for this loading mode (von-Mises distribution), the distribution of minimum principal strain of the continuum-level model was notably different compared to that of the micro-FE model.


Model A Model B

100

Variable 1: Total Agreement (%) Variable 2: Agreement within highly-stressed trabecular elements (%) Variable 3: Agreement within highly-stressed shell elements (%)

Agreement between micro-FE and continuum-level highly-stressed elements (%)

80 60 40 20 0

p<0.0001

p<0.001

p=0.37

Variable 3 Variable 1 Variable 2 Figure 5. Agreement of von-Mises and minimum principal strain distributions (75th-100th percentiles) between the micro-FE and continuum-level models. See Table 1 for descriptions for models A and B.

In general, the transverse layers in the continuum-model that had lower axial stiffness had greater number of elements experiencing high strains whereas this trend was not observed in the micro-FE models. In the continuum model, the trend for transverse slices with lower axial stiffness to have greater number of elements in the high-strain category was observed in 11 specimens (p<0.05), while two specimens showed no significant correlations. On the contrary, in the micro-FE models, there was only 1 specimen that had a significant negative correlation (p<0.05) between axial stiffness and number of elements in the highstrain category; 5 specimens had a significant positive correlation (p<0.05); while 7 specimens had no significant correlation. 4. Discussion

The goal of this study was to investigate the mechanistic fidelity of the continuum-level models by comparing its predictions with that of a more detailed micro-FE model. The absolute predictions of stiffness in the continuum-level models were more sensitive to the material properties assigned to the peripheral bone as compared to the trabecular centrum. The degree of correlation between stiffness predictions of the continuum-level and micro-FE m o d e l s was insensitive to the choice of the modulus-density relation, presumably because the spatial distribution of density was captured in all parametric variations. The slope of the stiffness correlation was substantially lower than an ideal slope (Y=X) for all parametric variations, indicating that the


tissue-level modulus assumed for the micro-FE models (E=18.5 GPa) may be too high (indeed, subsequent independent tests in our laboratory have shown the effective tissue modulus of the trabecular bone to be closer to 10 GPa.) The major load paths compared well between the continuum-level and micro-FE models (von-Mises distribution) since the high-modulus continuum elements carry significant load and thus, dominate the stress distribution. The largely transverse strain bands observed in the continuum-level models were consistent with a springs-in-series analogy wherein the highest strains were observed in the most compliant layers of the vertebral body. The volume-averaged strain distribution in the micro-FE model, however, was not consistent with this analogy and was similar to that of the stress distribution of micro-FE models, indicative of an axial tissue-level stress and strain state. Taken together, these results indicate continuum-level models of human vertebrae capture the stress distributions better that the strain distributions, and that whole-vertebra metrics such as stiffness are highly correlated across the scales of modeling. A major strength of this study was its focus on both the apparent-level and tissue-level predictions of the multi-scale models. The micro-CT models, using 60-micron-sized voxel elements, had sufficient geometric resolution to model the physics associated with bending and deformation of individual trabeculae. These models also included a detailed geometric description of the thin cortical shell. The substantial porosity between individual trabeculae is also included explicitly in these models. By contrast, the continuum models had no such porosity, the cortical shell was averaged with surrounding material, and the physics associated with deformation of individual trabeculae was absent. Coarsening of the micro-CT based vertebral body models allowed us to simulate the behavior of clinically-relevant continuum-level models and to directly relate the mechanics of the two scales. The results presented here represent the bestcase performance of the continuum-level models since errors due to such factors as image streaks or low signal-to-noise ratio -- both of which may be present in clinical QCT-scans -- have not been simulated in this study. Parametric variations of the modulus-density relationships in the continuum models enabled us to extract the relative roles of the peripheral bone and the trabecular c e n t r u m in determining the accuracy of the continuum-level stiffness predictions. Despite these strengths, there are a few notable limitations. First, the boundary conditions in these models were simplistic compared to physiological loading through a non-linearly and inhomogeneous elastic disc. This approach was justified at this juncture since these analyses represent a necessary first step in understanding the mechanics of the multi-scale models under controlled loading conditions without any confounding effects of the complex disc.


Second, linear analyses were performed in all models and hence, failure distribution or strength predictions were not compared between the multi-scale models. The results from the comparison of the elastic stress/strain distribution need not necessarily hold for failure distributions since there may be load redistribution around fractured trabeculae in the real bone. However, since the locations of high strain in the micro-FE model are likely representative of regions of initial failure [15], these results may well represent measures of a yield strength. This study provides contrasting results to that of a previous study on multiscale analysis of human femora. Verhulp et al. [3] found that both the stress and strain distributions of the continuum-level model compared well with that of the micro-FE model of a healthy and an osteoporotic femur during fall loading conditions. However, the micro-mechanics of the femur under fall loading conditions are much different than those for the spine under compression. Deformation mechanisms depend on the bone volume fraction, role of the cortical shell, and orientation of the trabecular bone with respect to the loading axis -- all different between hip and spine. Thus, the relation between multiscale strains is expected to depend on such factors as anatomic site and loading conditions and in that sense the results are entirely consistent. The results from this study provide unique insight into the mechanics of the multi-scale models. Consistent with previous studies which have shown that the peripheral, dense bone plays a large role in vertebral stiffness [9, 10], our results demonstrated that absolute predictions of stiffness were more sensitive to the material properties assigned to the bone around the periphery as compared to the bone within the trabecular centrum. This finding suggests that improved characterization of the continuum-level material properties of the peripheral bone may improve apparent-level predictions of strength at the continuum-level model. Using the current modeling approach of the vertebral body at the c o n t i n u u m scale, our results indicate that the volume-averaged stress distribution predicted by the continuum model in the elastic range may be more accurate compared to the strain distribution. However, the discrepancy for the l a t t e r is likely due to the volume-averaging process since the results demonstrated how this process allows high strains to occur in the continuum level model in regions of very low apparent density (and consequently, in transverse layers with low axial stiffness), whereas in the micro-mechanical model the strain in such regions were zero because of the absence of bone. Interesting, regions of high apparent strain appear to agree well with directly observed regions of failure in thick sections of vertebrae loaded to failure in the cadaver laboratory [8]. Thus, it remains to be seen if -- despite the disagreement between continuum-level and micro-level volume-averaged strain distributions


reported here -- the former remains highly correlated with directly observed failure regions in whole vertebrae. Further work is also recommended to extend these studies to more physiological type loading conditions, including load transfer through a disc with and without forward flexion. Acknowledgments Funding from National Institute of Health (AR49828). Computational resources from grant UCB-266 (NPACI). Finite element analyses performed on an IBM Power4 supercomputer (Datastar, San Diego Supercomputer Center). Specimen preparation by J. Buckley and micro-CT scanning by Dr. M.K. Liebschner (Rice University). Dr. Keaveny has a financial interest in O.N. Diagnostics and both he and the company may benefit from the results of this research. References 1. Faulkner, K.G., C.E. Cann, and B.H. Hasegawa, Effect of bone d i s t r i b u t i o n on vertebral strength: assessment with patient-specific nonlinear finite element analysis. Radiology, 1991. 179: 669-74. 2. Crawford, R.P., C.E. Cann, and T.M. Keaveny, Finite element models p r e d i c t in vitro vertebral body compressive strength better than quantitative computed tomography. Bone, 2003. 33: 744-50. 3. Verhulp, E., B. van Rietbergen, and R. Huiskes, Comparison of microlevel and continuum-level voxel models of the proximal femur. Journal of Biomechanics, 2005. 9: 2951-7. 4. Liebschner, M.A., et al., Finite element modeling of the human thoracolumbar spine. Spine, 2003. 28: 559-65. 5. Cao, K.D., M.J. Grimm, and K.H. Yang, Load sharing within a human lumbar vertebral body using the finite element method. Spine, 2001. 26: E253-60. 6. Silva, M.J., T.M. Keaveny, and W.C. Hayes, Load sharing between the shell and centrum in the lumbar vertebral body. Spine, 1997. 22: 140-50. 7. Homminga, J., et al., Osteoporosis changes the amount of vertebral trabecular bone at risk of fracture but not the vertebral load distribution. Spine, 2001. 26: 1555-61. 8. Silva, M.J., T.M. Keaveny, and W.C. Hayes, Computed tomographybased finite element analysis predicts failure loads and fracture patterns for vertebral sections. Journal of Orthopaedic Research, 1998. 16: 300-8. 9. Eswaran, S.K., et al., Cortical and trabecular load sharing in the human vertebral body. J Bone Miner Res, 2006. 21: 307-14. 10. Eswaran, S.K., et al., The micromechanics of cortical shell removal in the human vertebral body. Computer Methods in Applied Mechanics and Engineering. 2006. 196: 3025-32.


11. Bevill, G., et al., Influence of bone volume fraction and architecture on computed large-deformation failure mechanisms in human trabecular bone. Bone, 2006. 39: 1218-25. 12. Morgan, E.F., H.H. Bayraktar, and T.M. Keaveny, Trabecular bone modulus-density relationships depend on anatomic site. J Biomech, 2003. 36: 897-904. 13. Kopperdahl, D.L., E.F. Morgan, and T.M. Keaveny, Quantitative computed tomography estimates of the mechanical properties of human vertebral trabecular bone. Journal of Orthopaedic Research, 2002. 20: 8015. 14. Adams, M.F., et al. Ultrascalable implicit finite element analyses in solid mechanics with over a half a billion degrees of freedom. in ACM/IEEE Proceedings of SC2004: High Performance Networking and Computing. 2004. 15. Eswaran, S.K., A. Gupta, and T.M. Keaveny, Locations of bone tissue at high risk of initial failure during compressive loading of the human vertebral body. Bone, 2007. 41: 733-9.