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.     

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.     

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.     

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.     

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.     

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?     

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.     

Subject-specific finite element models of long bones: An in vitro evaluation of the overall accuracy

The determination of the mechanical stresses induced in human bones is of great importance in both research and clinical practice. Since the stresses in bones cannot be measured non-invasively in vivo, the only way to estimate them is through subject-specific finite element modelling. Several methods exist for the automatic generation of these models from CT data, but before bringing them in the clinical practice it is necessary to assess their accuracy in the predictions of the bone stresses. Particular attention should be paid to those regions, like the epiphyseal and metaphyseal parts of long bones, where the automatic methods are typically less accurate. Aim of the present study was to implement a general procedure to automatically generate subject-specific finite element models of bones from CT data and estimate the accuracy of this general procedure by applying it to one real femur. This femur was tested in vitro under five different loading scenarios and the results of these tests were used to verify how the adoption of a simplified two-material homogeneous model would change the accuracy with respect to the density-based inhomogeneous one, with special attention paid to the epiphyseal and metaphyseal proximal regions of the bone. The results showed that the density-based inhomogeneous model predicts with a very good accuracy the measured stresses (R(2)=0.91, RMSE=8.6%, peak error=27%), while the two-material model was less accurate (R(2)=0.89, RMSE=9.6%, peak error=35%). The results showed that it is possible to automatically generate accurate finite element models of bones from CT data and that the strategy of material properties mapping has a significant influence on its accuracy.     

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.     

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.     

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?     

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.     

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.     

The Role of the Size and Location of the Tumors and of the Vertebral Anatomy in Determining the Structural Stability of the Metastatically Involved Spine: a Finite Element Study

Vertebral fractures associated with the loss of structural integrity of neoplastic vertebrae are common, and determined to the deterioration of the bone quality in the lesion area. The prediction of the fracture risk in metastatically involved spines can guide in deciding if preventive solutions, such as medical prophylaxis, bracing, or surgery are indicated for the patient. In this study, finite element models of 22 thoracolumbar vertebrae were built based on CT scans of three spines, covering a wide spectrum of possible clinical scenarios in terms of age, bone quality and degenerative features, taking into account the local material properties of bone tissue. Simulations were performed in order to investigate the effect of the size and location of the tumoral lesion, the bone quality and the vertebral level in determining the structural stability of the neoplastic vertebrae. Tumors with random size and positions were added to the models, for a total of 660 simulations in which a compressive load was simulated. Results highlighted the fundamental role of the tumor size, whereas the other parameters had a lower, but non-negligible impact on the axial collapse of the vertebra, the vertebral bulge in the transverse plane and the canal narrowing under the application of the load. All the considered parameters are radiologically measurable, and can therefore be translated in a straightforward way to the clinical practice to support decisions about preventive treatment of metastatic fractures.     

NACOB presentation to ASB Young Scientist Award: Postdoctoral. The impact of boundary conditions and mesh size on the accuracy of cancellous bone tissue modulus determination using large-scale finite-element modeling. North American Congress on Biomechanics.

The apparent properties of cancellous bone are determined by a combination of both hard tissue properties and microstructural organization. A method is desired to extract the underlying hard tissue properties from simple mechanical tests, free from the complications of microstructure. It has been suggested that microCT voxel-based large-scale finite element models could be employed to accomplish this goal (van Rietbergen et al., 1995, Journal of Biomechanics, 28, 69-81). This approach has recently been implemented and it is becoming increasingly popular as finite element models increase in size and sophistication (Fyhrie et al., 1997, Proceedings of the 43rd Annual Meeting of the Orthopaedic Research Society, San Francisco, CA, p. 815; van Rietbergen et al., 1997, Proceedings of the 43rd Annual Meeting of the Orthopaedic Research Society, San Francisco, CA, p. 62). However, no direct quantitative measurements of the accuracy of this method applied to porous structures such as cancellous bone have been made. This project demonstrates the feasibility of this approach by quantifying its best-case accuracy in determining the trabecular hard tissue modulus of analogues fabricated of a material with known material properties determined independently by direct testing. In addition we were able to assess the impact of mesh size and boundary conditions on accuracy. We found that the assumption of a frictionless boundary condition in the parallel plate compression loading configuration was a significant source of error that could be overcome with the use of rigid end-caps similar to those used by Keaveny et al. (1997 Journal of Orthopaedic Research, 15(1), 101-110). In conclusion, we found that this approach is an effective method for determining the average trabecular hard tissue properties of human cancellous bone with an expected practical accuracy level better than 5%.     

A new approach for assigning bone material properties from CT images into finite element models

Generation of subject-specific finite element (FE) models from computed tomography (CT) datasets is of significance for application of the FE analysis to bone structures. A great challenge that remains is the automatic assignment of bone material properties from CT Hounsfield Units into finite element models. This paper proposes a new assignment approach, in which material properties are directly assigned to each integration point. Instead of modifying the dataset of FE models, the proposed approach divides the assignment procedure into two steps: generating the data file of the image intensity of a bone in a MATLAB program and reading the file into ABAQUS via user subroutines. Its accuracy has been validated by assigning the density of a bone phantom into a FE model. The proposed approach has been applied to the FE model of a sheep tibia and its applicability tested on a variety of element types. The proposed assignment approach is simple and illustrative. It can be easily modified to fit users’ situations.     

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.