Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Biomech. Author manuscript; available in PMC 2017 March 21.
Published in final edited form as:
PMCID: PMC4801716

A database of lumbar spinal mechanical behavior for validation of spinal analytical models


Data from two experimental studies with eight specimens each of spinal motion segments and/or intervertebral discs are presented in a form that can be used for comparison with finite element model predictions. The data include the effect of compressive preload (0, 250 and 500 N) with quasistatic cyclic loading (0.0115 Hz) and the effect of loading frequency (1, 0.1, 0.01 and 0.001 Hz) with a physiological compressive preload (mean 642 N). Specimens were tested with displacements in each of six degrees of freedom (three translations and three rotations) about defined anatomical axes. The three forces and three moments in the corresponding axis system were recorded during each test. Linearized stiffness matrices were calculated that could be used in multi-segmental biomechanical models of the spine and these matrices were analyzed to determine whether off-diagonal terms and symmetry assumptions should be included.

These databases of lumbar spinal mechanical behavior under physiological conditions quantify behaviors that should be present in finite element model simulations. The addition of more specimens to identify sources of variability associated with physical dimensions, degeneration, and other variables would be beneficial. Supplementary data provide the recorded data and Matlab® codes for reading files. Linearized stiffness matrices derived from the tests at different preloads revealed few significant unexpected off-diagonal terms and little evidence of significant matrix asymmetry.

Keywords: Spine, Lumbar Spine, Stiffness, Preload, Frequency-dependence, Database


Analytical models of the biomechanical function of the spine and spinal motion segments together with experimental studies provide understanding of in vivo spinal loading as well as the biomechanics of the intervertebral discs and motion segments. Finite element models have been used to gain understanding of elastic load-deformation behavior [Dreischarf et al. 2014, Weisse et al. 2012], time-dependent behavior [Castro et al. 2014, Galbusera et al., 2011, Schmidt et al. 2013, Schroeder et al. 2010] and insights into damage accumulation [Qasim et al. 2012], degeneration and transport of nutrients and metabolites [Natarajan, Williams and Andersson (2004), Schmidt et al., 2013]. Models are also used to predict the behavior of surgical implants [Zhang and Teo 2008], and in multi-segmental spine models to estimate in vivo spinal and muscular forces with each motion segment represented by a stiffness matrix or equivalent beam [Gardner-Morse and Stokes 2004]. While motion segment and disc models are becoming increasingly sophisticated by including complex elastic formulations, creep, viscoelasticity and swelling behavior, there is a shortage of experimental data for validation of the complex nonlinear, time-dependent six-degree of freedom behavior of the spine. Although tissue properties have been extensively studied [e.g. Cloyd et al. 2007, O’Connell, Sen and Elliott 2012., Wagnac et al., 2011], the structural behaviors of discs and motion segments including the neural arch and facet joints are not well documented.

The intervertebral disc is a complex avascular structure having nonlinear behavior in six degrees of freedom with stiffness that increases with axial compressive load [Gardner-Morse and Stokes 2004], and with time dependent response (creep, hysteresis, etc.). The disc is commonly modeled as a biphasic tissue with fluid and solid phases, also sometimes including fixed charges responsible for swelling behavior and retention of fluid. The nonlinear elastic behavior is thought to result primarily from the cable-like nature of collagen fibers, and its time-dependent behavior resulting from both fluid-flow effects and solid phase viscoelasticity [Costi et al., 2008].

Motion segment behavior in multi-segmental analyses of spinal and trunk loading and stability can be represented efficiently by use of a continuum formulation such as a linearized stiffness matrix or ‘equivalent’ beam [Gardner-Morse and Stokes 2004]. Most of the available motion segment data are limited to axial compression stiffness and creep, and elastic behavior in the three principal angular degrees of freedom. The full linear stiffness matrix as required for equivalent structure representation comprises both diagonal and off-diagonal terms and matrix symmetry. Goel et al. [1987] noted that certain forces or moments are associated with at least two displacements (in a flexibility experiment) and identified them as ‘primary’ and ‘secondary’. For example, an anterior force would generally produce anterior shear (primary) as well as flexion (secondary) motion. Secondary motion is often also identified as ‘coupling’ behavior. Also, any misalignment of the applied force relative to structural symmetry axes of the specimen would produce additional out-of-plane motion. Here in the context of a stiffness experiment (applied displacements) we identify primary terms as the diagonal terms of the stiffness matrix, and secondary (expected) off-diagonal terms; also additional off-diagonal terms resulting from structural asymmetry or testing axis misalignment. In experimentally determined stiffness matrices, translational stiffness should be independent of the center of the axis system, but rotational stiffness terms would be axis-center dependent. (The inverse is the case for a flexibility matrix).

The purpose of this paper is to provide data for use in validation of analytical models or other analyses of the lumbar spine. Data were recorded from two sets of human lumbar motion segments (without posterior elements for the second set) in six degrees of freedom, under different axial preloads and with different loading frequencies. The data have previously been used to report stiffening by preload and representation of the motion segment as an ‘equivalent’ beam [Gardner-Morse and Stokes 2004] and in a study of frequency-dependent apparent stiffness and hysteresis behavior under cyclic loading [Costi et al., 2008] with testing frequencies covering the range from quasistatic to walking.

Experimental data were also examined to identify significant off-diagonal terms in the derived linearized stiffness matrices and to determine their degree of symmetry. These matrices are symmetrical in a linear elastic structure, but experimentally have been found to be asymmetrical [Holsgrove et al., 2015] in porcine motion segments. Previously published experimentally derived stiffness matrices using this database [Gardner-Morse and Stokes 2004] assumed symmetry.


Data for the six-degrees of freedom behavior of human cadaveric motion segments were recorded using a ‘hexapod’ (Stewart platform) apparatus [Stokes et al., 2002]. This computer-controlled apparatus has six linear actuators, six displacement transducers and a six-degree-of freedom load cell. It was programmed to impose each of the three principal displacements and three principal rotations in the motion segment’s local axis system while the three principal forces and moments were recorded. Human motion segments were dissected from available human spines (thus a sample of convenience) that had been stored at −80 °C. Each specimen was radiographed and no evidence of anatomical abnormality or gross degeneration was observed. Some osteophytes were observed on the older specimens.

In the first ‘preload’ test series [Gardner-Morse and Stokes 2004], there were eight lumbar motion segments (L2-3 and L4-5 from each of four human females, aged 17, 21, 52 and 58 years). See Table 1. These specimens were tested in physiological saline at 4 °C. Each specimen was tested with axial compressive preloads of 0, 250 and 500 N (stress approximately 0, 0.15 and 0.3 MPa, comparable with the lower range of physiological values [Wilke et al., 1999]) and was allowed to equilibrate with each preload for at least three hours before each load-displacement test. The six tests (three pure translations and three pure rotations) were sequentially performed with four sawtooth-waveform (ramp-loading) cycles of 87 s (1.15 x 10−2 Hz) in each displacement direction. This loading rate was the slowest possible representing quasi-static conditions, compatible with an acceptable total testing time per specimen (~80 hours). The applied displacements and resulting forces were recorded at 1 Hz. The displacements were cycles of +/−0.5 mm in anterior–posterior and lateral displacements, +/−0.35 mm axial displacement, +/−1.5 degrees lateral bend rotation and +/−1 degrees flexion-extension and axial rotations. After testing each intact specimen, the facets and ligaments (posterior elements) were removed and the tests were repeated.

Table 1
Specimen details for database of lumbar spinal mechanical behavior.

In the second ‘frequency’ tests series [Costi et al., 2008], eight vertebra-disc-vertebra motion segments (posterior elements removed) from seven human lumbar spines (five males, two females, mean [SD] age: 41 [18] years, range: 16–58 years; levels: L1/L2 n=2, L2/L3 n=2, L3/L4 n=3, L4/L5 n=1; weight: 84 [19] kg) were tested. According to Thompson’s criteria [Thompson et al. 1987] for disc degeneration grade modified for transverse sections, 2 discs were grade 1, 5 were grade 3, and 1 was grade 5. See Table 1. Specimens were loaded with a sine waveform at each of four frequencies (0.001, 0.01, 0.1, and 1 Hz, with displacements and forces recorded at 0.2, 2, 32 and 128 Hz respectively) after equilibration overnight under a preload based on estimated area to produce 0.4 MPa to simulate in vivo static loading conditions [Costi et al., 2008]. The displacements amplitudes were +/−0.6 mm in anterior–posterior and lateral displacements, +/−0.25 mm axial displacement, +/−3 degrees flexion-extension rotation and +/−2 degrees lateral bend rotation and axial rotation. These specimens were tested immersed in phosphate buffered saline with protease inhibitors at 37 °C, the solution having a pH of 6.8. Total testing time per specimen was ~52 hours.

For the preload tests, the motion segment coordinate system was centered in the superior vertebral body with the Z-axis joining the centers of the vertebrae and the X and Y coordinates aligned with the anatomical planes as identified on biplanar radiographs [Stokes, et al., 2002] with +X anterior, +Y left and +Z superior. For the frequency tests, the coordinate system was similar but centered in the intervertebral disc. Thus the displacements/rotations were applied and the forces were measured at the center of the disc. The vertebral body center convention was appropriate for generating a stiffness matrix for continuum finite element analyses, whereas the center of the disc convention was appropriate for rigid body kinematics with viscoelastic joints with linear combinations of springs and dashpots (e.g. Maxwell, Kelvin–Voigt and/or standard linear solid models).

For each preload condition in the preload tests, and for the data obtained with and without posterior element, all 36 terms of a stiffness matrix were obtained by least-squares best fit (Matlab® ‘backslash’ operator) to experimental data. The resulting stiffness matrices were then analyzed statistically (by two-tailed t-tests) to identify which terms were significantly different from zero and whether pairs of off-diagonal terms expected to be equal by symmetry were equal. Non-zero off-diagonal terms would result from inherent coupling [Goel et al. 1987] and any experimental axis misalignment relative to the segment’s structural axis, while nonlinear behavior would also result in matrix asymmetry. For the data from the frequency tests, the primary terms of the stiffness matrix were obtained by linear regression analysis of recorded displacement-force and rotation-moment data.


The dataset, together with sample Matlab® programs for reading the files are contained in Supplementary Data. The programs generate the graphical output for the ‘preload’ data in Figure 1 and for the ‘frequency’ data in Figure 2. Analyses of the ‘preload’ dataset showed that off-diagonal (coupling) terms were significantly different from zero for all preload conditions and for specimens with and without posterior elements for lateral shear/axial rotation, for anterior shear/flexion-extension and for lateral shear/lateral bending. See Table 2. Anterior shear/flexion-extension and lateral shear/lateral bending coupling terms are ‘primary’ coupling terms according to Goel et al. [1987]. Also, the coupling terms for lateral bend/axial rotation were significantly different from zero in all cases except zero preload without posterior elements. The stiffness matrix for a beam with a rigid anterior-posterior offset also has these four coupling terms plus an additional coupling term between axial displacements and flexion-extension. Beyond these four coupling terms, there were only three cases in which both of the complementary coupling terms were both significantly different from zero. However, there was statistically significant evidence of matrix asymmetry (complementary stiffness terms not equal) in 16 of the 35 cases in which at least one of the coupling terms were significant. For the four significant coupling terms, only nine of the 23 significant cases (39.1%) were not symmetrical and the average asymmetry (calculated as the mean of the absolute values of differences between two values, divided by their average and expressed as a percentage) was 12.5%. See Table 2. Thus, overall, the linearized stiffness matrices included the expected off- diagonal coupling terms, and there was little evidence of additional coupling terms, nor of any substantial degree of matrix asymmetry for the range of input displacements in these tests. Load- displacement behavior was almost linear, except for axial displacement (Dz-direction).

Figure 1
Illustrative graphs of ‘preload’ data generated by Matlab® code plot_preload.m in the Supplemental data. The graphs for each degree of freedom show the mean and standard deviation (error bars) across specimens of displacement-increasing ...
Figure 2
Illustrative graphs of ‘frequency’ data generated by Matlab® code plot_freq.m in the Supplemental data. Anterior-posterior force Fx was recorded over time at each of four frequencies for all eight specimens under anterior-posterior ...
Table 2
Off-diagonal values in stiffness matrices (in units of N, mm and radians) obtained from the preload test series (mean of 8 motion segments).

The primary stiffnesses (diagonal terms in the stiffness matrices) were all somewhat smaller than the diagonal stiffness values in the frequency test series for all six degrees of freedom. See Table 3. The differences were probably due to a higher preload (mean 642 N) and possibly due to sex differences (Table 1), differences in the location of the axes center (for rotational degrees of freedom) and/or different dimensional or tissue properties.

Table 3
Mean ± standard deviation (SD) of diagonal stiffness terms of disc specimens without posterior elements.


These data provide information for use in validation of finite element models. By including the recordings in all six degrees of freedom for each principal displacement, they provide complete values of the nonlinear stiffness and hysteresis at any recorded loading condition. A finite element model that reproduces this behavior could then provide information about the state of strain at any location. The stiffening with preload effect is thought to result from the stiffness difference in axial compression (much less stiffness in tension than in compression), probably due to the cable-like properties of collagen fibers in the annulus, together with some large displacement (geometrical) effects (e.g. disc bulge).

Finite element models generally predict that stiffness depends on the state of disc degeneration [Galbusera et al. 2011, Natarajan, Williams and Andersson 2006], supported by some experimental evidence [Zirbel et al. 2013] and on disc dimensions, disc thickness and cross sectional area properties [Chagnon et al. 2010, Meijer et al. 2010, Natarajan and Andersson, 1999, Niemeyer, Wilke and Schmidt 2012]. However, no significant correlations between recorded behavior and dimensions and degeneration state were found by Gardner-Morse and Stokes [2004] and Berkson et al. [1979].

For the relatively small displacements in these experiments, there were few significant unexpected off-diagonal terms and little evidence of significant matrix asymmetry. Assuming testing misalignments are random, statistical tests provided a method to detect significant off-diagonal stiffness terms and significant asymmetries. In comparing the linearized stiffness of the two groups of specimens (Table 3), differences were expected for rotational degrees of freedom since the axis center for the preload tests was at the center of the upper vertebral body while it was at the disc center for the frequency data. Stiffness values in finite element analyses can be transformed to different axis systems by means of a rigid offset.

The variability in the stiffness of these specimens indicates the need for future work to include a greater number of specimens to identify the sources of the variability (age, gender, physical dimensions and state of degeneration, etc.) These relationships in turn should be compared with the predictions made by finite element modeling.

Supplementary Material


Data were collected with support from NIH grants R01 44119 and R01 049370


Conflict of interest statement

Both authors do not have any financial and personal relationships with other people or organizations that could inappropriately influence (bias) their work.

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


  • Berkson MH, Nachemson AL, Schultz AB. Mechanical properties of human lumbar spine motion segments-Part II: responses in compression and shear; influence of gross morphology. Journal of Biomechanical Engineering. 1979;101:53–57.
  • Castro AP, Wilson W, Huyghe JM, Ito K, Alves JL. Intervertebral disc creep behavior assessment through an open source finite element solver. Journal of Biomechanics. 2014;47(1):297–301. [PubMed]
  • Chagnon A, Aubin CE, Villemure I. Biomechanical influence of disk properties on the load transfer of healthy and degenerated disks using a poroelastic finite element model. Journal of Biomechanical Engineering. 2010;132(11):111006. [PubMed]
  • Cloyd JM, Malhotra NR, Weng L, Chen W, Mauck RL, Elliott DM. Material properties in unconfined compression of human nucleus pulposus, injectable hyaluronic acid-based hydrogels and tissue engineering scaffolds. European Spine Journal. 2007;11:1892–8. [PMC free article] [PubMed]
  • Costi JJ, Stokes IA, Gardner-Morse MG, Iatridis JC. Frequency dependent behavior of the intervertebral disc in response to each of six degree of freedom dynamic loading: Solid phase and fluid phase contributions. Spine. 2008;33(16):1731–8. [PMC free article] [PubMed]
  • Dreischarf M, Zander T, Shirazi-Adl A, Puttlitz CM, Adam CJ, Chen CS, Goel VK, Kiapour A, Kim YH, Labus KM, Little JP, Park WM, Wang YH, Wilke HJ, Rohlmann A, Schmidt H. Comparison of eight published static finite element models of the intact lumbar spine: predictive power of models improves when combined together. Journal of Biomechanics. 2014;47(8):1757–66. [PubMed]
  • Galbusera F, Schmidt H, Noailly J, Malandrino A, Lacroix D, Wilke HJ, Shirazi-Adl A. Comparison of four methods to simulate swelling in poroelastic finite element models of intervertebral discs. Journal of the Mechanical Behavior of Biomedical Materials. 2011;4(7):1234–41. [PubMed]
  • Galbusera F, Schmidt H, Neidlinger-Wilke C, Gottschalk A, Wilke HJ. The mechanical response of the lumbar spine to different combinations of disc degenerative changes investigated using randomized poroelastic finite element models. European Spine Journal. 2011;20(4):563–71. [PMC free article] [PubMed]
  • Gardner-Morse MG, Stokes IAF. Structural behavior of human lumbar spinal motion segments. Journal of Biomechanics. 2004;37(2):205–121. [PubMed]
  • Goel VK. Three-dimensional motion behavior of the human spine-a question of terminology. Journal of Biomechanical Engineering. 1987;109:353–355. [PubMed]
  • Holsgrove TP, Gill HS, Miles AW, Gheduzzi S. The dynamic, six-axis stiffness matrix testing of porcine spinal specimens. Spine Journal. 2015;15(1):176–84. [PubMed]
  • Meijer GJ, Homminga J, Hekman EE, Veldhuizen AG, Verkerke GJ. The effect of three-dimensional geometrical changes during adolescent growth on the biomechanics of a spinal motion segment. Journal of Biomechanics. 2010;43(8):1590–7. [PubMed]
  • Natarajan RN, Andersson GB. The influence of lumbar disc height and cross-sectional area on the mechanical response of the disc to physiologic loading. Spine. 1999;24(18):1873–81. [PubMed]
  • Natarajan RN, Williams JR, Andersson GB. Recent advances in analytical modeling of lumbar disc degeneration. Spine. 2004;29(23):2733–41. [PubMed]
  • Natarajan RN, Williams JR, Andersson GB. Modeling changes in intervertebral disc mechanics with degeneration. Journal of Bone and Joint Surgery. 2006;88(Suppl 2):36–40. [PubMed]
  • Niemeyer F, Wilke HJ, Schmidt H. Geometry strongly influences the response of numerical models of the lumbar spine--a probabilistic finite element analysis. Journal of Biomechanics. 2012;45(8):1414–23. [PubMed]
  • O’Connell GD, Sen S, Elliott DM. Human annulus fibrosus material properties from biaxial testing and constitutive modeling are altered with degeneration. Biomechanics Modeling Mechanobiology. 2012;11(3–4):493–503. [PMC free article] [PubMed]
  • Qasim M, Natarajan RN, An HS, Andersson GB. Initiation and progression of mechanical damage in the intervertebral disc under cyclic loading using continuum damage mechanics methodology: A finite element study. Journal of Biomechanics. 2012;45(11):1934–40. [PMC free article] [PubMed]
  • Schmidt H, Galbusera F, Rohlmann A, Shirazi-Adl A. What have we learned from finite element model studies of lumbar intervertebral discs in the past four decades? Journal of Biomechanics. 2013;46(14):2342–55. [PubMed]
  • Schmidt H, Shirazi-Adl A, Galbusera F, Wilke HJ. Response analysis of the lumbar spine during regular daily activities--a finite element analysis. Journal of Biomechanics. 2010;43(10):1849–56. [PubMed]
  • Schroeder Y, Huyghe JM, van Donkelaar CC, Ito K. A biochemical/biophysical 3D FE intervertebral disc model. Biomechanics Modeling Mechanobiology. 2010;9:641–650. [PubMed]
  • Stokes IAF, Gardner-Morse M, Churchill D, Laible JP. Measurement of a spinal motion segment stiffness matrix. Journal of Biomechanics. 2002;35:517–521. [PubMed]
  • Thompson JP, Pearce RH, Schechter MT, et al. Preliminary evaluation of a scheme for grading the gross morphology of the human intervertebral disc. Spine. 1990;15:411–5. [PubMed]
  • Wagnac E, Arnoux PJ, Garo A, El-Rich M, Aubin CE. Calibration of hyperelastic material properties of the human lumbar intervertebral disc under fast dynamic compressive loads. Journal of Biomechanical Engineering. 2011;133(10):101007. [PubMed]
  • Weisse B, Aiyangar AK, Affolter Ch, Gander R, Terrasi GP, Ploeg H. Determination of the translational and rotational stiffnesses of an L4-L5 functional spinal unit using a specimen-specific finite element model. Journal of the Mechanical Behavior of Biomedical Materials. 2012;13:45–61. [PubMed]
  • Wilke HJ, Neef P, Caimi M, Hoogland T, Claes LE. New in vivo measurements of pressures in the intervertebral disc in daily life. Spine. 1999;24(8):755–62. [PubMed]
  • Zhang QH, Teo EC. Finite element application in implant research for treatment of lumbar degenerative disc disease. Medical Engineering & Physics. 2008;30(10):1246–56. [PubMed]
  • Zirbel SA, Stolworthy DK, Howell LL, Bowden AE. Intervertebral disc degeneration alters lumbar spine segmental stiffness in all modes of loading under a compressive follower load. Spine Journal. 2013;13(9):1134–47. [PubMed]