|Home | About | Journals | Submit | Contact Us | Français|
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.
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.  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.
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  years, range: 16–58 years; levels: L1/L2 n=2, L2/L3 n=2, L3/L4 n=3, L4/L5 n=1; weight: 84  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. . 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).
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.
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  and Berkson et al. .
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.
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.