Assessment of factors influencing finite element vertebral model predictions

This study aimed to establish model construction and configuration procedures for future vertebral finite element analysis by studying convergence, sensitivity, and accuracy behaviors of semiautomatically generated models and comparing the results with manually generated models. During a previous study, six porcine vertebral bodies were imaged using a microcomputed tomography scanner and tested in axial compression to establish their stiffness and failure strength. Finite element models were built using a manual meshing method. In this study, the experimental agreement of those models was compared with that of semiautomatically generated models of the same six vertebrae. Both manually and semiautomatically generated models were assigned gray-scale-based, element-specific material properties. The convergence of the semiautomatically generated models was analyzed for the complete models along with material property and architecture control cases. A sensitivity study was also undertaken to test the reaction of the models to changes in material property values, architecture, and boundary conditions. In control cases, the element-specific material properties reduce the convergence of the models in comparison to homogeneous models. However, the full vertebral models showed strong convergence characteristics. The sensitivity study revealed a significant reaction to changes in architecture, boundary conditions, and load position, while the sensitivity to changes in material property values was proportional. The semiautomatically generated models produced stiffness and strength predictions of similar accuracy to the manually generated models with much shorter image segmentation and meshing times. Semiautomatic methods can provide a more rapid alternative to manual mesh generation techniques and produce vertebral models of similar accuracy. The representation of the boundary conditions, load position, and surrounding environment is crucial to the accurate prediction of the vertebral response. At present, an element size of 2x2x2 mm(3) appears sufficient since the error at this size is dominated by factors, such as the load position, which will not be improved by increasing the mesh resolution. Higher resolution meshes may be appropriate in the future as models are made more sophisticated and computational processing time is reduced.     

Factors affecting subject-specific finite element models of implant-fitted rat bone specimens: critical analysis of a technical protocol

The authors propose a protocol to derive finite element (FE) models from micro computer tomography scans of implanted rat bone. A semi-automatic procedure allows segmenting the images using specimen-specific bone mineral density (BMD) thresholds. An open-source FE model generator processes the segmented images to a quality tetrahedral mesh. The material properties assigned to each element are integrated from the BMD field. Piecewise, threshold-dependent density–elasticity relationships are implemented to limit the effects of metal artefacts. A detailed sensitivity study highlights the coherence of the generated models and quantifies the influence of the modelling parameters on the results. Two applications of the protocol are proposed. The stiffness of bare and implanted rat tibiae specimens is predicted by simulating three-point bending and inter-implant displacement, respectively. Results are compared with experimental tests. The mean value and the variability between the specimens are well captured in both tests.     

Assessment of the effect of mesh density on the material property discretisation within QCT based FE models: a practical example using the implanted proximal tibia

A three-dimensional, quantitative computed tomography based finite element model of a proximal implanted tibia was analysed in order to assess the effect of mesh density on material property discretisation and the resulting influence on the predicted stress distribution. The mesh was refined on the contact surfaces (matched meshes) with element sizes of 3, 2, 1.4, 1 and 0.8 mm. The same loading conditions were used in all models (bi-condylar load: 60% medial, 40% lateral). Significant variations were observed in the modulus distributions between the coarsest and finest mesh densities. Poor discretisation of the material properties also resulted in poor correlations of the stresses and risk ratios between the coarsest and finest meshes. Little difference in Young's modulus, von Mises stress and risk ratio distributions were observed between the three finest models; hence, it was concluded that for this particular case an element size of 1.4 mm on the contact surfaces was enough to properly describe the stiffness, stress and risk ratio distributions within the bone. Poor convergence of the material property distribution occurred when the element size was significantly larger than the pixel size of the source CT data. It was concluded that unless there is convergence in the Young's modulus distribution, convergence of the stress field or of other parameters of interest will not occur either.     

Localized damage in vertebral bone is most detrimental in regions of high strain energy density

It was hypothesized that damage to bone tissue would be most detrimental to the structural integrity of the vertebral body if it occurred in regions with high strain energy density, and not necessarily in regions of high or low trabecular bone apparent density, or in a particular anatomic location. The reduction in stiffness due to localized damage was computed in 16 finite element models of 10-mm-thick human vertebral sections. Statistical analyses were performed to determine which characteristic at the damage location--strain energy density, apparent density, or anatomic location--best predicted the corresponding stiffness reduction. There was a strong positive correlation between regional strain energy density and structural stiffness reduction in all 16 vertebral sections for damage in the trabecular centrum (p < 0.05, r2 = 0.43-0.93). By contrast, regional apparent density showed a significant negative correlation to stiffness reduction in only four of the sixteen bones (p < 0.05, r2 = 0.47-0.58). While damage in different anatomic locations did lead to different reductions in stiffness (p < 0.0001, ANOVA), no single location was consistently the most critical location for damage. Thus, knowledge of the characteristics of bone that determine strain energy density distributions can provide an understanding of how damage reduces whole bone mechanical properties. A patient-specific finite element model displaying a map of strain energy density can help optimize surgical planning and reinforcement of bone in individuals with high fracture risk.     

A patient-specific computer tomography-based finite element methodology to calculate the six dimensional stiffness matrix of human vertebral bodies

In most finite element (FE) studies of vertebral bodies, axial compression is the loading mode of choice to investigate structural properties, but this might not adequately reflect the various loads to which the spine is subjected during daily activities or the increased fracture risk associated with shearing or bending loads. This work aims at proposing a patient-specific computer tomography (CT)-based methodology, using the currently most advanced, clinically applicable finite element approach to perform a structural investigation of the vertebral body by calculation of its full six dimensional (6D) stiffness matrix. FE models were created from voxel images after smoothing of the peripheral voxels and extrusion of a cortical shell, with material laws describing heterogeneous, anisotropic elasticity for trabecular bone, isotropic elasticity for the cortex based on experimental data. Validated against experimental axial stiffness, these models were loaded in the six canonical modes and their 6D stiffness matrix calculated. Results show that, on average, the major vertebral rigidities correlated well or excellently with the axial rigidity but that weaker correlations were observed for the minor coupling rigidities and for the image-based density measurements. This suggests that axial rigidity is representative of the overall stiffness of the vertebral body and that finite element analysis brings more insight in vertebral fragility than densitometric approaches. Finally, this extended patient-specific FE methodology provides a more complete quantification of structural properties for clinical studies at the spine.     

Effects of suppression of bone turnover on cortical and trabecular load sharing in the canine vertebral body

The relative biomechanical effects of antiresorptive treatment on cortical thickness vs. trabecular bone microarchitecture in the spine are not well understood. To address this, T-10 vertebral bodies were analyzed from skeletally mature female beagle dogs that had been treated with oral saline (n=8 control) or a high dose of oral risedronate (0.5 mg/kg/day, n=9 RIS-suppressed) for 1 year. Two linearly elastic finite element models (36-μm voxel size) were generated for each vertebral body—a whole-vertebra model and a trabecular-compartment model—and subjected to uniform compressive loading. Tissue-level material properties were kept constant to isolate the effects of changes in microstructure alone. Suppression of bone turnover resulted in increased stiffness of the whole vertebra (20.9%, p=0.02) and the trabecular compartment (26.0%, p=0.01), while the computed stiffness of the cortical shell (difference between whole-vertebra and trabecular-compartment stiffnesses, 11.7%, p=0.15) was statistically unaltered. Regression analyses indicated subtle but significant changes in the relative structural roles of the cortical shell and the trabecular compartment. Despite higher average cortical shell thickness in RIS-suppressed vertebrae (23.1%, p=0.002), the maximum load taken by the shell for a given value of shell mass fraction was lower (p=0.005) for the RIS-suppressed group. Taken together, our results suggest that—in this canine model—the overall changes in the compressive stiffness of the vertebral body due to suppression of bone turnover were attributable more to the changes in the trabecular compartment than in the cortical shell. Such biomechanical studies provide an unique insight into higher-scale effects such as the biomechanical responses of the whole vertebra.     

Assessment of the Effect of Mesh Density on the Material Property Discretisation Within QCT Based FE Models: A Practical Example Using the Implanted Proximal Tibia

A three-dimensional, quantitative computed tomography based finite element model of a proximal implanted tibia was analysed in order to assess the effect of mesh density on material property discretisation and the resulting influence on the predicted stress distribution. The mesh was refined on the contact surfaces (matched meshes) with element sizes of 3, 2, 1.4, 1 and 0.8 mm. The same loading conditions were used in all models (bi-condylar load: 60% medial, 40% lateral). Significant variations were observed in the modulus distributions between the coarsest and finest mesh densities. Poor discretisation of the material properties also resulted in poor correlations of the stresses and risk ratios between the coarsest and finest meshes. Little difference in Young's modulus, von Mises stress and risk ratio distributions were observed between the three finest models; hence, it was concluded that for this particular case an element size of 1.4 mm on the contact surfaces was enough to properly describe the stiffness, stress and risk ratio distributions within the bone. Poor convergence of the material property distribution occurred when the element size was significantly larger than the pixel size of the source CT data. It was concluded that unless there is convergence in the Young's modulus distribution, convergence of the stress field or of other parameters of interest will not occur either.     

Personalised Mechanical Properties of Scoliotic Vertebrae Determined In Vivo Using Tomodensitometry

This in vivo study investigated the mechanical properties of apical scoliotic vertebrae using computed tomography (CT) and finite element (FE) meshing. CT examination was performed on seven scoliotic girls. FE meshing of each vertebral body allowed automatic mapping of the CT scan and the visualisation of the bone density distribution. Centroids and mass centres were compared to analyse the mechanical properties distribution. Compared to the centroid, the mass centre migrated into the concavity of the curvature. The three vertebrae of a same patient had the same body migration behaviour because they were located at the curvature apex. This observation was verified in the coronal plane, but not in the sagittal plane. These results represent new data over few geometrical analyses of scoliotic vertebrae. Same in vivo personalisation of mechanical properties should be performed on intervertebral discs or ligaments to personalise stiffness properties of the spine for the biomechanical modelling of human torso. Moreover, do this mechanical deformation of scoliotic vertebrae, that appears before the vertebral wedging, could be a predictive tool in scoliosis treatment?     

Discordance between Prevalent Vertebral Fracture and Vertebral Strength Estimated by the Finite Element Method Based on Quantitative Computed Tomography in Patients with Type 2 Diabetes Mellitus

Background

Bone fragility is increased in patients with type 2 diabetes mellitus (T2DM), but a useful method to estimate bone fragility in T2DM patients is lacking because bone mineral density alone is not sufficient to assess the risk of fracture. This study investigated the association between prevalent vertebral fractures (VFs) and the vertebral strength index estimated by the quantitative computed tomography-based nonlinear finite element method (QCT-based nonlinear FEM) using multi-detector computed tomography (MDCT) for clinical practice use.

Research Design and Methods

A cross-sectional observational study was conducted on 54 postmenopausal women and 92 men over 50 years of age, all of whom had T2DM. The vertebral strength index was compared in patients with and without VFs confirmed by spinal radiographs. A standard FEM procedure was performed with the application of known parameters for the bone material properties obtained from nondiabetic subjects.

Results

A total of 20 women (37.0%) and 39 men (42.4%) with VFs were identified. The vertebral strength index was significantly higher in the men than in the women (P<0.01). Multiple regression analysis demonstrated that the vertebral strength index was significantly and positively correlated with the spinal bone mineral density (BMD) and inversely associated with age in both genders. There were no significant differences in the parameters, including the vertebral strength index, between patients with and without VFs. Logistic regression analysis adjusted for age, spine BMD, BMI, HbA1c, and duration of T2DM did not indicate a significant relationship between the vertebral strength index and the presence of VFs.

Conclusion

The vertebral strength index calculated by QCT-based nonlinear FEM using material property parameters obtained from nondiabetic subjects, whose risk of fracture is lower than that of T2DM patients, was not significantly associated with bone fragility in patients with T2DM. This discordance may indirectly suggest that patients with T2DM have deteriorated bone material compared with nondiabetic subjects, a potential cause of bone fragility in T2DM patients.     

Parametric finite element analysis of vertebral bodies affected by tumors

The vertebral column is the most frequent site of metastatic involvement of the skeleton. Due to the proximity to the spinal cord, from 5% to 10% of all cancer patients develop neurologic manifestations. As a consequence, fracture risk prediction has significant clinical importance. In this study, we model the metastatically involved vertebra so as to parametrically investigate the effects of tumor size, material properties and compressive loading rate on vertebral strength. A two-dimensional axisymmetric finite element model of a spinal motion segment consisting of the first lumbar vertebral body (no posterior elements) and adjacent intervertebral disc was developed to allow the inclusion of a centrally located tumor in the vertebral body. After evaluating elastic, mixed, and poroelastic formulations, we concluded that the poroelastic representation was most suitable for modeling the metastatically involved vertebra's response to compressive load. Maximum principal strains were used to localize regions of potential vertebral trabecular bone failure. Radial and axial vertebral body displacements were used as relative indicators of spinal canal encroachment and endplate failure. Increased tumor size and loading rate, and reduced trabecular bone density all elevated axial and radial displacements and maximum tensile strains. The results of this parametric study suggest that vertebral tumor size and bone density contribute significantly to a patients risk for vertebral fracture and should be incorporated in clinical assessment paradigms.     

Simulation of the behaviour of the L1 vertebra for different material properties and loading conditions

Three-dimensional finite element models of the thoracolumbar junction (T12–L2) and isolated L1 vertebra were developed to investigate the role of material properties and loading conditions on vertebral stresses and strains to predict fracture risk. The geometry of the vertebrae was obtained from computed tomography images. The isolated vertebra model included an L1 vertebra loaded through polymethylmethacrylate plates located at the top and bottom of the vertebra, and the segment model included T12 to L2 vertebrae and seven ligaments, fibrous intervertebral discs and facet joints. Each model was examined with both homogeneous and spatially varying bone tissue properties. Stresses and strains were compared for uniform compression and flexion. Including material heterogeneity remarkably reduced the stiffness of the isolated L1 vertebra and increased the magnitudes of the minimum principal strains and stresses in the mid-transverse section. The stress and strain distributions further changed when physiological loading was applied to the L1 vertebra. In the segment models, including heterogeneous material properties increased the magnitude of the minimum principal strain by 158% in the centre of the mid-transverse section. Overall, the inclusion of heterogeneity and physiological loading increased the magnitude of the strains up to 346% in flexion and 273% in compression.     

Finite element analysis of vertebral body mechanics with a nonlinear microstructural model for the trabecular core.

In this study, a finite element model of a vertebral body was used to study the load-bearing role of the two components (shell and core) under compression. The model of the vertebral body has the characteristic kidney shape transverse cross section with concave lateral surfaces and flat superior and inferior surfaces. A nonlinear unit cell based foam model was used for the trabecular core, where nonlinearity was introduced as coupled elastoplastic beam behavior of individual trabeculae. The advantage of the foam model is that architecture and material properties are separated, thus facilitating studies of the effects of architecture on the apparent behavior. Age-related changes in the trabecular architecture were considered in order to address the effects of osteoporosis on the load-sharing behavior. Stiffness changes with age (architecture and porosity changes) for the trabecular bone model were shown to follow trends in published experimental results. Elastic analyses showed that the relative contribution of the shell to the load-bearing ability of the vertebra decreases with increasing age and lateral wall curvature. Elasto-plastic (non-linear) analyses showed that failure regions were concentrated in the upper posterior region of the vertebra in both the shell and core components. The ultimate load of the vertebral body model varied from 2800 N to 5600 N, depending on age (architecture and porosity of the trabecular core) and shell thickness. The model predictions lie within the range of experimental results. The results provide an understanding of the relative role of the core and shell in vertebral body mechanics and shed light on the yield and post-yield behavior of the vertebral body.     

A patient-specific finite element methodology to predict damage accumulation in vertebral bodies under axial compression, sagittal flexion and combined loads

Due to the inherent limitations of DXA, assessment of the biomechanical properties of vertebral bodies relies increasingly on CT-based finite element (FE) models, but these often use simplistic material behaviour and/or single loading cases. In this study, we applied a novel constitutive law for bone elasticity, plasticity and damage to FE models created from coarsened pQCT images of human vertebrae, and compared vertebral stiffness, strength and damage accumulation for axial compression, anterior flexion and a combination of these two cases. FE axial stiffness and strength correlated with experiments and were linearly related to flexion properties. In all loading modes, damage localised preferentially in the trabecular compartment. Damage for the combined loading was higher than cumulated damage produced by individual compression and flexion. In conclusion, this FE method predicts stiffness and strength of vertebral bodies from CT images with clinical resolution and provides insight into damage accumulation in various loading modes.     

A nonlocal constitutive model for trabecular bone softening in compression

Using the three-dimensional morphological data provided by computed tomography, finite element (FE) models can be generated and used to compute the stiffness and strength of whole bones. Three-dimensional constitutive laws capturing the main features of bone mechanical behavior can be developed and implemented into FE software to enable simulations on complex bone structures. For this purpose, a constitutive law is proposed, which captures the compressive behavior of trabecular bone as a porous material with accumulation of irreversible strain and loss of stiffness beyond its yield point and softening beyond its ultimate point. To account for these features, a constitutive law based on damage coupled with hardening anisotropic elastoplasticity is formulated using density and fabric-based tensors. To prevent mesh dependence of the solution, a nonlocal averaging technique is adopted. The law has been implemented into a FE software and some simple simulations are first presented to illustrate its behavior. Finally, examples dealing with compression of vertebral bodies clearly show the impact of softening on the localization of the inelastic process.     

Personalised mechanical properties of scoliotic vertebrae determined in vivo using tomodensitometry

This in vivo study investigated the mechanical properties of apical scoliotic vertebrae using computed tomography (CT) and finite element (FE) meshing. CT examination was performed on seven scoliotic girls. FE meshing of each vertebral body allowed automatic mapping of the CT scan and the visualisation of the bone density distribution. Centroids and mass centres were compared to analyse the mechanical properties distribution. Compared to the centroid, the mass centre migrated into the concavity of the curvature. The three vertebrae of a same patient had the same body migration behaviour because they were located at the curvature apex. This observation was verified in the coronal plane, but not in the sagittal plane. These results represent new data over few geometrical analyses of scoliotic vertebrae. Same in vivo personalisation of mechanical properties should be performed on intervertebral discs or ligaments to personalise stiffness properties of the spine for the biomechanical modelling of human torso. Moreover, do this mechanical deformation of scoliotic vertebrae, that appears before the vertebral wedging, could be a predictive tool in scoliosis treatment?     

Quantitative computed tomography-based finite element models of the human lumbar vertebral body: effect of element size on stiffness,damage, and fracture strength predictions

This study investigated the numerical convergence characteristics of specimen-specific "voxel-based" finite element models of 14 excised human cadaveric lumbar vertebral bodies (age: 37-87; M = 6, F = 8) that were generated automatically from clinical-type CT scans. With eventual clinical applications in mind, the ability of the model stiffness to predict the experimentally measured compressive fracture strength of the vertebral bodies was also assessed. The stiffness of "low"-resolution models (3 x 3 x 3 mm element size) was on average only 4% greater (p = 0.03) than for "high"-resolution models (1 x 1 x 1.5 mm) despite interspecimen variations that varied over four-fold. Damage predictions using low- vs high-resolution models were significantly different (p = 0.01) at loads corresponding to an overall strain of 0.5%. Both the high (r2 = 0.94) and low (r2 = 0.92) resolution model stiffness values were highly correlated with the experimentally measured ultimate strength values. Because vertebral stiffness variations in the population are much greater than those that arise from differences in voxel size, these results indicate that imaging resolution is not critical in cross-sectional studies of this parameter. However, longitudinal studies that seek to track more subtle changes in stiffness over time should account for the small but highly significant effects of voxel size. These results also demonstrate that an automated voxel-based finite element modeling technique may provide an excellent noninvasive assessment of vertebral strength.     

Hexahedral meshing of subject-specific anatomic structures using mapped building blocks

To extend the use of computational techniques like finite element analysis to clinical settings, it would be beneficial to have the ability to generate a unique model for every subject quickly and efficiently. This work is an extension of two previously developed mapped meshing tools that utilised force and displacement control to map a template mesh to a subject-specific surface. The objective of this study was to map a template block structure, common to multiblock meshing techniques, to a subject-specific surface. The rationale is that the blocks are considerably less refined and may be readily edited after mapping, thereby yielding a mesh of high quality in less time than mapping the mesh itself. In this paper, the versatility and robustness of the method was verified by processing four data-sets. The method was found to be robust enough to cope with the variability of bony surface size, spatial position and geometry, producing building block structures (BBSs) that generated meshes comparable to those produced using BBSs that were created manually.     

Discretization error when using finite element models: Analysis and evaluation of an underestimated problem

Mesh convergence tests are often insufficiently performed in finite element analyses. There are many parameters which may have an effect on the mesh convergence behavior. The aim of this study was to identify the influence of different parameters on the mesh convergence behavior.For this purpose we used a simplified axis-symmetrical model of a single pedicle screw flank with surrounding bone to simulate a pull-out test. In parameter studies, the flank radii and the contact conditions at the bone–screw interface were varied. These parameter studies were carried out using an implicit and explicit solver. Thereby, the convergence criteria and the number of the substeps for the implicit nonlinear iteration process as well as the velocity and the material density for the explicit approach were considered.The mesh convergence behavior was influenced by varying the flank radii and the contact conditions. The implicit calculations led to a reaction force, which converged rapidly to a certain value with increasing mesh density, whereas the maximum von-Mises stress showed substantial convergence problems. The number of substeps and the convergence criteria of the iteration process strongly influenced the implicit solutions. In contrast, the maximum von-Mises stresses resulting from explicit calculations converged to a certain value after only a few refinement steps. Different pull-out velocities substantially affected the mesh convergence behavior, while the material density showed only a negligible influence.The results indicated the need to perform an appropriate mesh convergence test when using finite element methods. We were able to show that different parameters strongly influence the mesh convergence behavior and we demonstrated that convergence tests do not always lead to a satisfactory or acceptable solution.     

HR-pQCT-based homogenised finite element models provide quantitative predictions of experimental vertebral body stiffness and strength with the same accuracy as μFE models

This study validated two different high-resolution peripheral quantitative computer tomography (HR-pQCT)-based finite element (FE) approaches, enhanced homogenised continuum-level (hFE) and micro-finite element (μFE) models, by comparing them with compression test results of vertebral body sections. Thirty-five vertebral body sections were prepared by removing endplates and posterior elements, scanned with HR-pQCT and tested in compression up to failure. Linear hFE and μFE models were created from segmented and grey-level CT images, and apparent model stiffness values were compared with experimental stiffness as well as strength results. Experimental and numerical apparent elastic properties based on grey-level/segmented CT images (N=35) correlated well for μFE (r2=0.748/0.842) and hFE models (r2=0.741/0.864). Vertebral section stiffness values from the linear μFE/hFE models estimated experimental ultimate apparent strength very well (r2=0.920/0.927). Calibrated hFE models were able to predict quantitatively apparent stiffness with the same accuracy as μFE models. However, hFE models needed no back-calculation of a tissue modulus or any kind of fitting and were computationally much cheaper.     

A finite element model technique to determine the mechanical response of a lumbar spine segment under complex loads

This study presents a CT-based finite element model of the lumbar spine taking into account all function-related boundary conditions, such as anisotropy of mechanical properties, ligaments, contact elements, mesh size, etc. Through advanced mesh generation and employment of compound elements, the developed model is capable of assessing the mechanical response of the examined spine segment for complex loading conditions, thus providing valuable insight on stress development within the model and allowing the prediction of critical loading scenarios. The model was validated through a comparison of the calculated force-induced inclination/deformation and a correlation of these data to experimental values. The mechanical response of the examined functional spine segment was evaluated, and the effect of the loading scenario determined for both vertebral bodies as well as the connecting intervertebral disc.