Parametric convergence sensitivity and validation of a finite element model of the human lumbar spine

The primary objective of this study was to generate a finite element model of the human lumbar spine (L1–L5), verify mesh convergence for each tissue constituent and perform an extensive validation using both kinematic/kinetic and stress/strain data. Mesh refinement was accomplished via convergence of strain energy density (SED) predictions for each spinal tissue. The converged model was validated based on range of motion, intradiscal pressure, facet force transmission, anterolateral cortical bone strain and anterior longitudinal ligament deformation predictions. Changes in mesh resolution had the biggest impact on SED predictions under axial rotation loading. Nonlinearity of the moment-rotation curves was accurately simulated and the model predictions on the aforementioned parameters were in good agreement with experimental data. The validated and converged model will be utilised to study the effects of degeneration on the lumbar spine biomechanics, as well as to investigate the mechanical underpinning of the contemporary treatment strategies.     

Lumbar spine finite element model for healthy subjects: development and validation

Finite element (FE) method is a proven powerful and efficient tool to study the biomechanics of the human lumbar spine. However, due to the large inter-subject variability of geometries and material properties in human lumbar spines, concerns existed on the accuracy and predictive power of one single deterministic FE model with one set of spinal geometry and material properties. It was confirmed that the combined predictions (median or mean value) of several distinct FE models can be used as an improved prediction of behavior of human lumbar spine under identical loading and boundary conditions. In light of this improved prediction, five FE models (L1-L5 spinal levels) of the human lumbar spine were developed based on five healthy living subjects with identical modeling method. The five models were extensively validated through experimental and computational results in the literature. Mesh convergence and material sensitivity analysis were also conducted. We have shown that the results from the five FE models developed in this paper were consistent with the experimental data and simulation results from the existing literature. The validated modeling method introduced in this study can be used in modeling dysfunctional lumber spines such as disc degeneration and scoliosis in future work.     

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.     

C4–C5 segment finite element model development,validation, and load-sharing investigation

Detailed cervical spine models are necessary to better understand cervical spine response to loading, improve our understanding of injury mechanisms, and specifically for predicting occupant response and injury in auto crash scenarios. The focus of this study was to develop a C4–C5 finite element model with accurate representations of each tissue within the segment. This model incorporates more than double the number of elements of existing models, required for accurate prediction of response. The most advanced material data available were then incorporated using appropriate nonlinear constitutive models to provide accurate predictions of response at physiological levels of loading. This tissue-scale segment model was validated against a wide variety of experimental data including different modes of loading (axial rotation, flexion, extension, lateral bending, and translation), and different load levels. In general, the predicted response of the model was within the single standard deviation response corridors for both low and high load levels. Importantly, this model demonstrates that appropriate refinement of the finite element mesh, representation at the tissue level, and sufficiently detailed material properties and constitutive models provide excellent response predictions without calibration of the model to experimental data. Load sharing between the disc, ligaments, and facet joints was investigated for various modes of loading, and the dominant load-bearing structure was found to correlate with typical anatomical injury sites for these modes of loading. The C4–C5 model forms the basis for the development of a full cervical spine model. Future studies will focus on tissue-level injury prediction and dynamic response.     

Are coupled rotations in the lumbar spine largely due to the osseo-ligamentous anatomy?--a modeling study

Prior studies have found that primary rotations in the lumbar spine are accompanied by coupled out-of-plane rotations. However, it is not clear whether these accompanying rotations are primarily due to passive (discs, ligaments and facet joints) or active (muscles) spinal anatomy. The aim of this study was to use a finite element (FE) model of the lumbar spine to predict three-dimensional coupled rotations between the lumbar vertebrae, due to passive spinal structures alone. The FE model was subjected to physiologically observed whole lumbar spine rotations about in vivo centres of rotation. Model predictions were validated by comparison of intra-discal pressures and primary rotations with in vivo measurements and these showed close agreement. Predicted coupled rotations matched in vivo measurements for all primary motions except lateral bending. We suggest that coupled rotations accompanying primary motions in the sagittal (flexion/extension) and transverse (axial rotation) planes are primarily due to passive spinal structures. For lateral bending the muscles most likely play a key role in the coupled rotation of the spine.     

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.     

Finite element lumbar spine facet contact parameter predictions are affected by the cartilage thickness distribution and initial joint gap size

Current finite element modeling techniques utilize geometrically inaccurate cartilage distribution representations in the lumbar spine. We hypothesize that this shortcoming severely limits the predictive fidelity of these simulations. Specifically, it is unclear how these anatomically inaccurate cartilage representations alter range of motion and facet contact predictions. In the current study, cadaveric vertebrae were serially sectioned, and images were taken of each slice in order to identify the osteochondral interface and the articulating surface. A series of custom-written algorithms were utilized in order to quantify each facet joint's three-dimensional cartilage distribution using a previously developed methodology. These vertebrae-dependent thickness cartilage distributions were implemented on an L1 through L5 lumbar spine finite element model. Moments were applied in three principal planes of motion, and range of motion and facet contact predictions from the variable thickness and constant thickness distribution models were determined. Initial facet gap thickness dimensions were also parameterized. The data indicate that the mean and maximum cartilage thickness increased inferiorly from L1 to L5, with an overall mean thickness value of 0.57 mm. Cartilage distribution and initial facet joint gap thickness had little influence on the lumbar range of motion in any direction, whereas the mean contact pressure, total contact force, and total contact area predictions were altered considerably. The data indicate that range of motion predictions alone are insufficient to establish model validation intended to predict mechanical contact parameters. These data also emphasize the need for the careful consideration of the initial facet joint gap thickness with respect to the spinal condition being studied.     

Biomechanical implications of lumbar spinal ligament transection

Many lumbar spine surgeries either intentionally or inadvertently damage or transect spinal ligaments. The purpose of this work was to quantify the previously unknown biomechanical consequences of isolated spinal ligament transection on the remaining spinal ligaments (stress transfer), vertebrae (bone remodelling stimulus) and intervertebral discs (disc pressure) of the lumbar spine. A finite element model of the full lumbar spine was developed and validated against experimental data and tested in the primary modes of spinal motion in the intact condition. Once a ligament was removed, stress increased in the remaining spinal ligaments and changes occurred in vertebral strain energy, but disc pressure remained similar. All major biomechanical changes occurred at the same spinal level as the transected ligament, with minor changes at adjacent levels. This work demonstrates that iatrogenic damage to spinal ligaments disturbs the load sharing within the spinal ligament network and may induce significant clinically relevant changes in the spinal motion segment.     

Three-dimensional geometrical and mechanical modelling of the lumbar spine.

The main objective of this study is to design a three-dimensional geometrical and mechanical finite element model of the lumbar spine. The model's geometry is constructed using six parameters per vertebra. These parameters are digitized from two X-rays (anterio-posterior and lateral), thus yielding an individualized model which can be arrived at from the radiographs of a tested specimen. This procedure makes the model validation easier, as geometry is generally a factor of dispersion in experimental results. The geometrical reconstruction, in the form of a finite elements mesh, was effected for the whole lumbar spine. The global coherence of the model was verified.     

Application of homogenization theory to the study of trabecular bone mechanics.

It is generally accepted that the strength and stiffness of trabecular bone is strongly affected by trabecular microstructure. It has also been hypothesized that stress induced adaptation of trabecular bone is affected by trabecular tissue level stress and/or strain. At this time, however, there is no generally accepted (or easily accomplished) technique for predicting the effect of microstructure on trabecular bone apparent stiffness and strength or estimating tissue level stress or strain. In this paper, a recently developed mechanics theory specifically designed to analyze microstructured materials, called the homogenization theory, is presented and applied to analyze trabecular bone mechanics. Using the homogenization theory it is possible to perform microstructural and continuum analyses separately and then combine them in a systematic manner. Stiffness predictions from two different microstructural models of trabecular bone show reasonable agreement with experimental results, depending on metaphyseal region, (R2 greater than 0.5 for proximal humerus specimens, R2 less than 0.5 for distal femur and proximal tibia specimens). Estimates of both microstructural strain energy density (SED) and apparent SED show that there are large differences (up to 30 times) between apparent SED (as calculated by standard continuum finite element analyses) and the maximum microstructural or tissue SED. Furthermore, a strut and spherical void microstructure gave very different estimates of maximum tissue SED for the same bone volume fraction (BV/TV). The estimates from the spherical void microstructure are between 2 and 20 times greater than the strut microstructure at 10-20% BV/TV.     

Quantitative comparison of ligament formulation and pre-strain in finite element analysis of the human lumbar spine

Data has been published that quantifies the nonlinear, anisotropic material behaviour and pre-strain behaviour of the anterior longitudinal, supraspinous (SSL), and interspinous ligaments of the human lumbar spine. Additionally, data has been published on localized material properties of the SSL. These results have been incrementally incorporated into a previously validated finite element model of the human lumbar spine. Results suggest that the effects of increased ligament model fidelity on bone strain energy were moderate and the effects on disc pressure were slight, and do not justify a change in modelling strategy for most clinical applications. There were significant effects on the ligament stresses of the ligaments that were directly modified, suggesting that these phenomena should be included in FE models where ligament stresses are the desired metric.     

Intervertebral reaction force prediction using an enhanced assembly of OpenSim models

OpenSim offers a valuable approach to investigating otherwise difficult to assess yet important biomechanical parameters such as joint reaction forces. Although the range of available models in the public repository is continually increasing, there currently exists no OpenSim model for the computation of intervertebral joint reactions during flexion and lifting tasks. The current work combines and improves elements of existing models to develop an enhanced model of the upper body and lumbar spine. Models of the upper body with extremities, neck and head were combined with an improved version of a lumbar spine from the model repository. Translational motion was enabled for each lumbar vertebrae with six controllable degrees of freedom. Motion segment stiffness was implemented at lumbar levels and mass properties were assigned throughout the model. Moreover, body coordinate frames of the spine were modified to allow straightforward variation of sagittal alignment and to simplify interpretation of results. Evaluation of model predictions for level L1–L2, L3–L4 and L4–L5 in various postures of forward flexion and moderate lifting (8 kg) revealed an agreement within 10% to experimental studies and model-based computational analyses. However, in an extended posture or during lifting of heavier loads (20 kg), computed joint reactions differed substantially from reported in vivo measures using instrumented implants. We conclude that agreement between the model and available experimental data was good in view of limitations of both the model and the validation datasets. The presented model is useful in that it permits computation of realistic lumbar spine joint reaction forces during flexion and moderate lifting tasks. The model and corresponding documentation are now available in the online OpenSim repository.     

A microstructural finite element simulation of mechanically induced bone formation.

A finite element method to simulate the formation of an interconnected trabectular bone microstructure oriented with respect to applied in vivo mechanical forces is introduced and quantitatively compared to experimental data from a hydraulic bone chamber implant model. Randomly located 45 microm mineralized nodules were used as the initial condition for the model simulations to represent an early stage of intramembranous bone formation. Boundary conditions were applied consistent with the mechanical environment provided by the in vivo bone chamber model. A two-dimensional repair simulation algorithim that incorporated strain energy density (SED), SED gradient, principal strain, or principal strain gradient as the local objective criterion was utilized to simulate the formation of an oriented trabecular bone microstructure. The simulation solutions were convergent, unique, and relatively insensitive to the assumed initial distribution of mineralized nodules. Model predictions of trabecular bone morphology and anisotropy were quantitatively compared to experimental results. All simulations produced structures that qualitatively resembled oriented trabecular bone. However only simulations utilizing a gradient objective criterion yielded results quantitatively similar to in vivo observations. This simulation approach coupled with an experimental model that delivers controlled in vivo mechanical stimuli can be utilized to study the relationship between physical factors and microstructural adaptation during bone repair.     

A biomechanical model of the lumbosacral joint during lifting activities

A biomechanical model of the lumbosacral region was constructed for the purpose of systematically studying the combined stresses and strains on the local ligaments, muscles and disc tissue during sagittal plane two-handed lifting. The model was validated in two ways. The first validation was a comparison of experimental study results with model predictions. In general predictions compared very reasonably with observed values of several authors with the exception of strain predictions on the articular ligaments. Second, a sensitivity analysis was performed over a wide range of lifting tasks. The predicted stress/strain values followed anticipated patterns and were of reasonable magnitudes. On the basis of the results of the sensitivity analysis it was concluded that typical lifting tasks can lead to excessive disc compressive forces, muscle moment generation requirements, and possibly lumbodorsal fascia strains. Conversely, annulus rupture of a healthy disc due to overstrain appears very unlikely.     

Development of a model based method for investigating facet articulation

Reported investigations of facet articulation in the human spine have often been conducted through the insertion of pressure sensitive film into the joint space, which requires incision of the facet capsule and may alter the characteristics of interaction between the facet surfaces. Load transmission through the facet has also been measured using strain gauges bonded to the articular processes. While this method allows for preservation of the facet capsule, it requires extensive instrumentation of the spine, as well as strain-gauge calibration, and is highly sensitive to placement and location of the strain gauges. The inherently invasive nature of these techniques makes it difficult to translate them into medical practice. A method has been developed to investigate facet articulation through the application of test kinematics to a specimen-specific rigid-body model of each vertebra within a lumbar spine segment. Rigid-body models of each vertebral body were developed from CT scans of each specimen. The distances between nearest-neighboring points on each facet surface were calculated for specific time frames of each specimen's flexion/extension test. A metric describing the proportion of each facet surface within a distance (2 mm) from the neighboring surface, the contact area ratio (CAR), was calculated at each of these time frames. A statistically significant difference (p<0.037) was found in the CAR between the time frames corresponding to full flexion and full extension in every level of the lumbar spine (L1-L5) using the data obtained from the seven specimens evaluated in this study. The finding that the contact area of the facet is greater in extension than flexion corresponds to other findings in the literature, as well as the generally accepted role of the facets in extension. Thus, a biomechanical method with a sufficiently sensitive metric is presented as a means to evaluate differences in facet articulation between intact and treated or between healthy and pathologic spines.     

How does the geometry affect the internal biomechanics of a lumbar spine bi-segment finite element model? Consequences on the validation process

Numerical modelling can provide a thorough understanding of the mechanical influence on the spinal tissues and may offer explanations to mechanically linked pathologies. Such objective might be achieved only if the models are carefully validated. Sensitivity study must be performed in order to evaluate the influence of the approximations inherent to modelling. In this study, a new geometrically accurate L3-L5 lumbar spine bi-segmental finite-element model was acquired by modifying a previously existing model. The effect of changes in bone geometry, ligament fibres distribution, nucleus position and disc height was investigated in flexion and extension by comparison of the results obtained from the model before and after the geometrical update. Additional calculations were performed in axial rotation and lateral bending in order to compare the computed ranges of motion (ROM) with experimental results. It was found that the geometrical parameters affected the stress distribution and strain energy in the zygapophysial joints, the ligaments, and the intervertebral disc, changing qualitatively and quantitatively their relative role in resisting the imposed loads. The predicted ROM were generally in good agreement with the experimental results, independently of the geometrical changes. Hence, although the model update affected its internal biomechanics, no conclusions could be drawn from the experimental data about the validation of a particular geometry. Hence the validation of the lumbar spine model should be based on the relative role of its structural components and not only on its global mobility.     

Search for critical loading condition of the spine--a meta analysis of a nonlinear viscoelastic finite element model

The relative vulnerability of spinal motion segments to different loading combinations remains unknown. The meta-analysis described here using the results of a validated L2-L3 nonlinear viscoelastic finite element model was designed to investigate the critical loading and its effect on the internal mechanics of the human lumbar spine. A Box-Behnken experimental design was used to design the magnitude of seven independent variables associated with loads, rotations and velocity of motion. Subsequently, an optimization method was used to find the primary and secondary variables that influence spine mechanical output related to facet forces, disc pressure, ligament forces, annulus matrix compressive/shear stresses and anulus fibers strain. The mechanical responses with respect to the two most-relevant variables were then regressed linearly using the response surface quadratic model. Axial force and sagittal rotation were identified as the most-relevant variables for mechanical responses. The procedure developed can be used to find the critical loading for finite element models with multi input variables. The derived meta-models can be used to predict the risk associated with various loading parameters and in setting safer load limits.     

Computational modeling of three-dimensional electrodiffusion in biological systems: application to the node of Ranvier

A computational model is presented for the simulation of three-dimensional electrodiffusion of ions. Finite volume techniques were used to solve the Poisson-Nernst-Planck equation, and a dual Delaunay-Voronoi mesh was constructed to evaluate fluxes of ions, as well as resulting electric potentials. The algorithm has been validated and applied to a generalized node of Ranvier, where numerical results for computed action potentials agree well with cable model predictions for large clusters of voltage-gated ion channels. At smaller channel clusters, however, the three-dimensional electrodiffusion predictions diverge from the cable model predictions and show a broadening of the action potential, indicating a significant effect due to each channel's own local electric field. The node of Ranvier complex is an elaborate organization of membrane-bound aqueous compartments, and the model presented here represents what we believe is a significant first step in simulating electrophysiological events with combined realistic structural and physiological data.