Journal of Biomechanics 48 (2015) 479–488

Contents lists available at ScienceDirect

Journal of Biomechanics journal homepage: www.elsevier.com/locate/jbiomech www.JBiomech.com

Nonlinear dynamics of the human lumbar intervertebral disc Giacomo Marini a,n, Gerd Huber b, Klaus Püschel c, Stephen J. Ferguson a a b c

Institute for Biomechanics, ETH Zurich, Hönggerbergring 64, HPP-O14, 8093 Zurich, Switzerland Institute of Biomechanics, TUHH Hamburg University of Technology, Denickestraße 15, K, 21073 Hamburg, Germany Institute for Forensic Medicine, University Medical Center Hamburg—Eppendorf, Butenfeld 34, D-22529 Hamburg, Germany

art ic l e i nf o

a b s t r a c t

Article history: Accepted 3 December 2014

Systems with a quasi-static response similar to the axial response of the intervertebral disc (i.e. progressive stiffening) often present complex dynamics, characterized by peculiar nonlinearities in the frequency response. However, such characteristics have not been reported for the dynamic response of the disc. The accurate understanding of disc dynamics is essential to investigate the unclear correlation between whole body vibration and low back pain. The present study investigated the dynamic response of the disc, including its potential nonlinear response, over a range of loading conditions. Human lumbar discs were tested by applying a static preload to the top and a sinusoidal displacement at the bottom of the disc. The frequency of the stimuli was set to increase linearly from a low frequency to a high frequency limit and back down. In general, the response showed nonlinear and asymmetric characteristics. For each test, the disc had different response in the frequency-increasing compared to the frequency-decreasing sweep. In particular, the system presented abrupt changes of the oscillation amplitude at specific frequencies, which differed between the two sweeps. This behaviour indicates that the system oscillation has a different equilibrium condition depending on the path followed by the stimuli. Preload and amplitude of the oscillation directly influenced the disc response by changing the nonlinear dynamics and frequency of the jump-phenomenon. These results show that the characterization of the dynamic response of physiological systems should be readdressed to determine potential nonlinearities. Their direct effect on the system function should be further investigated. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Intervertebral disc Nonlinear dynamic Softening Hardening Jump-phenomenon

1. Introduction The mechanical response of the intervertebral disc (IVD) is due to the interaction, composition and structural organization of its components at different scales. This provides unique load bearing abilities to the disc and a flexibility that is often referred to as nonlinear, i.e. progressive stiffening (White and Panjabi, 1978). In particular, the uniaxial quasi-static response resembles the mechanics of stiffening systems, where the stiffness increases with the deformation. Several studies have investigated the disc response for different loading scenario. In particular, past studies were focused on the quasi-static and low frequency (o10 Hz) response to pure and combined cyclic loading (Costi et al., 2008; Koeller et al., 1986; Walsh and Lotz, 2004). However, such loading conditions are only partially representative of the in-situ loads that the IVD normally experiences. During occupational whole body vibration (WBV), a combined static and dynamic load acts on the spine. The internal static load is determined by the body mass, posture and the load bearing due to the task-specific payload or

n

Corresponding author. Tel.: þ 41446334440; fax: þ41446331453. E-mail address: [email protected] (G. Marini).

http://dx.doi.org/10.1016/j.jbiomech.2014.12.006 0021-9290/& 2014 Elsevier Ltd. All rights reserved.

equipment. The dynamic load is caused by the high frequency stimuli of the body due to machine vibrations or schocks from the environment (Seidel et al., 1998). Low back pain disorders are often related to WBV, with a direct influence on the IVD. Longterm exposure was seen to be detrimental to normal disc metabolism (Bovenzi and Hulschof, 1999; Urban et al., 2004), whereas short term exposure to high frequency and large amplitude hydrostatic stress was reported to be beneficial for protein synthesis and reduction of protein degradation (Kasra et al., 2003). Although the guidelines proposed by the International Organization for Standardization (ISO-2631-1, 1997; ISO-2631-5, 2005) establish threshold values for daily vibration dose, the precise relationship between vibration dose values and risk of injuries is unclear (Mansfield, 2005). Further, comparison of evaluations and assessments methods, suggested by these standards, reveals discrepancies in the chosen evaluating parameters, and proposed daily exposure limits (Griffin, 1998). Since in-vivo studies evaluate the risk based on the individual perceived parameters (i.e. level of discomfort, perception of vibration) (Thamsuwan et al., 2013), potential vibration-induced injuries could be not perceived if they do not involve nervous receptors. Studies have suggested fatigue failure of vertebral endplates by WBV as cause of subsequent degenerative changes of the lumbar spine, and in particular of the

G. Marini et al. / Journal of Biomechanics 48 (2015) 479–488

480

IVD (Huber et al., 2010; Sandover, 1998; Schwarze et al., 1998). Therefore, the objective characterization of IVD dynamics is extremely relevant to reduce the risks of degenerative processes triggered by exposure to high frequency loading. However, these studies which have investigated the IVD dynamics provided different outcome with respect to the frequency range of the vertical mode (Izambert et al., 2003; Kasra et al., 1992). This could be due to nonlinearities of the system oscillation which normally affects the dynamic of systems with similar nonlinear quasi-static response. The aim of this experimental study was to investigate the nonlinear dynamic response of human IVDs. The influence of the amplitude of the stimuli and the preload on the oscillation was a specific focus of the study.

2.3. Setup dynamic test A test fixture was designed to duplicate a base excitation model, where the upper part of the body was simulated by a dead weight, referred to as the preload (m), and a cyclic displacement (y) was applied at the bottom of the specimen (Fig. 1a–c). The preload was applied by weights (m/4) fastened to a four-arm crossbar fixed to the top embedding cup. Four pillars structurally coupled the motion of the bottom embedding cup to the top mounted actuator of a servo-hydraulic machine (MTSs Bionix, USA). To limit the vibrational degree of freedom to the longitudinal direction, a guiding cylinder was mounted on the top embedding cup, and its oscillation was constrained by eight Teflon holders fixed on the pillars. Silicone lubricant was sprayed on the contact surfaces to reduce friction. Accelerometers (4508 B 004, Brüel & Kjaer, Nærum, Denmark) .. recorded the acceleration (sampling frequency 11 kHz) of the disc (x) and the frame (ÿ). A camera system (VICONs, Oxford, United Kingdom) tracked the displacement of markers on the disc (two markers) and on the frame (one marker). The volume of view was adjusted to achieve 400 fps.

2. Material and methods

2.4. Test protocol

2.1. Cadaveric material and its screening

Before the dynamic test, a quasi-static compression-tension test was performed to precondition the specimen and to obtain the nonlinear quasi-static response (five cycles, compression 1000 N, tension 200 N, displacement rate 0.09 mm/s). The setup for the test is illustrated in Fig. 1d. According to the base excitation model, the dynamic test was performed by applying a sinusoidal displacement to the lower embedding cup (y). The displacement amplitude (yi) was kept constant and the system was driven by a sinusoidal sweep function:

Twelve male human cadaver lumbar spines (T12-S1) were harvested at autopsy (Institute for Forensic Medicine, University Medical Center Hamburg-Eppendorf, Hamburg, Germany). They were sealed in plastic bags and stored at  20 1C. Frozen specimens were scanned by computed tomography scanner (Brilliance 16, Philips Healthcare, Hamburg, Germany, settings: 120 kVp; 1 mm slices). The CT-images of the tested specimens were segmented to evaluate the area of the inferior and superior disc boundary, and the volume of the disc. The disc height was calculated by dividing the volume by the mean boundary area. Mean gross morphology data, age, and classification of each specimen are reported in Table 1. 2.2. Specimen Before testing, each spine was thawed overnight at room temperature. On the day of testing muscle tissue and posterior elements were removed. The anterior and posterior longitudinal ligaments were kept intact. Superior and inferior vertebra were accurately clean to expose the cortical bone. To improve fixation between the vertebrae and the mould, screws were inserted into the vertebrae. The specimen was then embedded cranially and caudally in two steel cups with polymethylmethacrylate (Technovits 7100, Heraus Kulzer GmbH, Wehrheim, Germany). The specimen position was adjusted to align the disc’s longitudinal axis to the vertical direction of the test fixture. IVDs were kept moist by spraying with Ringer solution and by wrapping with wet gauze. After the mechanical testing, the specimens were cut in the mid sagittal plane, photographed, and degeneration degree graded (Thompson et al., 1990). The disc was then scraped off the vertebra to reach the endplates and macroscopically assess them for fracture.

y ¼ yi sin ð2 π f ðt Þ t Þ The frequency was set to increase linearly from a low frequency to a high frequency boundary (Increasing Frequency Sweep, IFS) and back down (Decreasing Frequency Sweep, DFS). The frequency sweep rate was set to 0.1 Hz/s, after verifying that it was low enough to have a transient steady response of the system. The frequency range was defined by the analysis of the last cycle of the quasi-static test, which was force-averaged. For small stimuli the dynamic of nonlinear system can be approximated with linear model, hence, assuming a linear undamped model, the ideal resonant frequency would be given by: sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kpreload 1 f id ¼ 2π preload The initial frequency range was defined by f id 76 Hz (value defined during the pilot study). The influence of the preload was investigated by repeating the dynamic test with different weights, from 16 kg to 56 kg in steps of 10 kg. For each preload, different yi were used: 0.1 mm, 0.1 mm, 0.15 mm, 0.2 mm 0.15 mm, 0.1 mm, 0.1 mm. After the test with 56 kg-preload, a test with 16 kg-preload and yi of 0.1 mm was repeated three times. Tests were coded [preload]kg–[yi]mm, e.g. 36 kg–0.1 mm. Three minutes were allowed for creep to occur after the application of each load increment.

Table 1 Characteristics of the tested specimens with the classification based on the comparison of the quasi-static test before and after the dynamic tests (specimen status, I ¼intact, F ¼failed, FE¼ failed-endplate), and on the dynamic response exhibited in the test sequence (H ¼hardening, S¼ softening). No.

Age (yr)

Disc level

Disc grade

Volume (mm  3)

Height (mm)

Status

Dynamic response

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

21 21 28 31 31 34 34 35 35 39 39 47 47 56 56 61 63 63 67 67 69 69

L2L3 L4L5 L1L2 L1L2 L4L5 L2L3 L4L5 L2L3 L4L5 L2L3 L4L5 L2L3 L4L5 L2L3 L4L5 L4L5 L2L3 L4L5 L2L3 L4L5 L2L3 L4L5

1 1 1 1 1 3 3 1 1 2 2 2 2 3 3 2 3 3 2 2 3 3

12,921 16,261 14,491 14,504 16,565 9,794 11,367 14,319 17,208 15,052 16,717 22,739 20,974 13,284 9,839 12,095 18,926 20,185 14,099 17,407 10,399 10,184

8.5 10.7 9.8 9.9 11.1 6.4 7.1 8.7 9.6 10.5 10.8 10.8 10.2 8.7 6.4 7.4 9.6 10.2 8.5 9.6 7.9 7.1

F I I F I F I F I F I FE FE FE F I FE I FE FE FE F

H þS S H þS H þS H þS H þS H þS H þS H þS H þS S H þS S H þS H þS H þS H þS S H þS S H þS H þS

G. Marini et al. / Journal of Biomechanics 48 (2015) 479–488

481

Fig. 1. Mechanical system analogue for the characterization of the disc dynamic through a base-excitation model (a) and developed testing fixture (b and c). The excitation (y) was transmitted from the MTS-actuator to the base of the disc through a four-pillar cage. The upper part of the body was loaded by a mass (m) distributed over four dead weights (m/4). These were connected to the disc by a four-arm crossbar. The centre of mass of this construct lay approximately at the height of the specimen. A numerical eigenfrequency analysis of the crossbar’s structure showed that the first mode of oscillation lies above 210 Hz. The oscillation of the mass (x) perpendicular to the excitation direction – due to any imperfections – was reduced by a guiding cylinder whose movement was radially restrained by low-contact-area Teflon blocks. Acceleration of the parts (ÿ, x€ ) was recorded by accelerometers placed at the center of the crossbar and at the MTS flange (see exploded view, c). The latter was chosen to avoid contamination of the sensors after checking that the signal was similar to the acceleration measured on the low embedding cup. The input and output displacement were recorded through a camera system (VICON). The markers were placed on two arms of the crossbar, to check for any wobbling, and on the MTS-flange near the accelerometer. To perform the quasi-static test, the embedded specimen was mounted on a sliding table linked on a six degree of freedom load cell (d). The top embedding cup was connected to the pillars by a coupling flange.

Fig. 2. Gain and phase of the transmissibility function vs. frequency of the control test performed with a rubber damper with amplitude of the stimuli (yi) of 0.5 mm and preload of 36 kg (a) and 46 kg (b). Assuming an undamped model, the analysis of the quasi-static test suggested a resonant frequency for the axial mode of 9.4 Hz for 36 kg preload and 8.6 Hz for 46 kg preload. The plot distinguishes between Increasing Frequency Sweep (blue) and Decreasing Frequency Sweep (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 2.5. Data evaluation

The crossbar displacement was computed by averaging the signal of the two optical markers. Acceleration and displacement signals were filtered by a low-pass

filter (Butterworth, passband ripple 0.1 dB, stopband ripple 200 dB, passband frequency 140 Hz, stopband frequency 200 Hz) to remove the noise. The filtering induced a smoothing of the crest of about 4.5%. The peak values (maximum and minimum within one cycle) of the relative displacement between the crossbar and base (x–y), with respect to the reference state identified by the specific preload,

482

G. Marini et al. / Journal of Biomechanics 48 (2015) 479–488

were selected and plotted against the input frequency, distinguishing between IFS and DFS. To evaluate the intensity of the response with respect to the specific frequency, the transmissibility function was computed, T(ω). This is the ratio between the Fourier transform of the input and the Fourier transform of the output. In particular, the acceleration was used because of its higher sampling frequency, T ðωÞ ¼ X€ ðωÞ=Ÿ ðωÞ. Of each preload-sequence, only the second (0.1 mm), third (0.15 mm) and fourth (0.2 mm) test were considered in the analysis. The first test (0.1 mm) was used to check the effect of the mass change. The tests with descending amplitude were used to investigate changes during consecutive tests.

2.6. Statistical analysis The data groups showed a normal distribution (Shapiro–Wilk) and equal variance (Levene and Brown–Forsythe) after a logarithmic transformation was applied. Since the same specimen was tested sequentially with different loads, the repeated measures Student’s t-test and ANOVA were used. A post-hoc analysis was performed to investigate patterns between groups with specific preload and amplitude of stimuli. The assumed significant level (0.05) was corrected for multiple comparisons (Bonferroni and Tukey-Kramer).

Fig. 3. Specimen classification based on the force–displacement curves of the quasi static elongation–compression test for a failed specimen, no. 20 (a), and for an intact specimen, no. 5 (b). Intact specimens did not have significant differences in compression response, whereas the failed specimens presented an increase hysteresis cycle and larger maximum displacement in compression and tension after compared to before the test. Plots distinguish between the response before (green) and after (blue) the dynamic test. By the analysis of the mean change of specimen height normalized with respect to the initial disc height (Table 1) over the dynamic test, the specimens were further sorted in: intact specimens, failed specimens without endplates fracture and failed specimens with endplates fracture (c). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

G. Marini et al. / Journal of Biomechanics 48 (2015) 479–488

483

2.7. Control study

3.2. Height loss

The reliability of the measurement system and testing procedure were verified by performing control tests on a rubber specimen. Although the system exhibited a nonlinear quasi-static response, major nonlinearities were not evident in the system dynamics (Fig. 2). The resonant frequency for the longitudinal mode indicated by the test with a preload of 36 kg (9.9 Hz) and 46 kg (9.0 Hz) differed by circa þ 5.4% from the frequency obtained by the analysis of the quasi-static response assuming an undamped linear model (9.36 Hz and 8.59 Hz).

The sudden failure of the specimens was evaluated by analysing the change of specimens’ height recorded by the camera system which was normalized by the initial disc heights. Each specimen had a height reduction due to the dynamic test (Fig. 3c). However, the intact group had a significantly lower height loss compared to the failed group (p ¼0.005). Half of the failed specimens presented a macroscopic endplate fracture which was recognizable by a significant sudden drop of the specimen height for consecutive tests (p ¼0.0055). The absolute change in specimen height significantly correlated (Spearman, R¼ 0.81, po 0.01) with the sum of the absolute change between the maximum displacement (Δe þ Δc) reached before and after the dynamic test (Fig. 3a and b).

3. Results 3.1. Quasi-static response Depending on the specimen, the quasi-static response before and after the dynamic test differed. Two-thirds of the specimens presented an increased area of the hysteresis cycle (p o0.001) and larger maximum displacement in extension (p o0.04) and compression (p o0.001) compared to before the dynamic test (Fig. 3a). The remaining discs maintained a similar response in compression (p ¼0.65) but they had a softer response in extension (p ¼0.033) after the dynamic test (Fig. 3b). We will refer to the latter set as intact and to the former as failed.

3.3. Dynamic response The dynamic response of the IVD showed nonlinear asymmetric characteristics. In general, the system response was different in the IFS compared to the DFS with abrupt changes of oscillation at different frequencies (p o0.002) and with different

Fig. 4. Peaks of the relative displacement between the crossbar and base (x–y) vs. input frequency for the test with 16 kg preload and amplitude of the forced displacement of 0.1 mm (a) and 0.15 mm (b) for an L4L5 disc specimen. The plot distinguishes between Increasing Frequency Sweep (blue) and Decreasing Frequency Sweep (red). The .. plots on the left are the analytical response of undamped nonlinear models described by mxþ kxþ μx3 ¼f0sin(2πft) (Thomson, 1981). The frequency axis was normalized with respect to the natural frequency fn ¼sqrt(k/m)/2π. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

484

G. Marini et al. / Journal of Biomechanics 48 (2015) 479–488

initial and final amplitude (Fig. 4). Such a jump-phenomenon (JP) is characteristic for the dynamic behaviour of nonlinear systems, where the oscillation changes without following a continuous path. Two different nonlinear behaviours were seen in the disc response. A pure softening response (SR), where decreasing the frequency of the stimuli, within the frequency range of the instabilities, is associated with an increase of the oscillation amplitude (Fig. 4a); and a hardening-like response (HLR), where

the oscillation amplitude grows when increasing the frequency of the stimuli (Fig. 4b). Five specimens had only a SR in all the tests (Table 1), whereas the other discs, depending on the amplitude of the stimuli and applied preload, had either a SR or a HLR (Fig. 5). After failure of the endplate, the specimens exhibited only a SR. The tests were allocated to two groups: intact endplate tests (IET) and failed endplate tests (FET). FET comprises the tests after failure of the endplate occurs. The IET contains the remaining tests, since

Fig. 5. Peaks of the relative displacement between the crossbar and base (x–y) vs. input frequency vs. preload at forced displacement (yi) of 0.1 mm (a), 0.15 mm (b) and 0.2 mm (c) for specimen no. 5. The plots distinguish between Increasing Frequency Sweep (blue) and Decreasing Frequency Sweep (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

G. Marini et al. / Journal of Biomechanics 48 (2015) 479–488

485

Table 2 Boundary frequencies of the jump-phenomenon and peak of the transmissibility magnitude, |T|, for the intact endplate tests (IET) and failed endplate tests (FET). Data were grouped based on the dynamic response exhibited in the test sequence (HL¼ hardening-like, S¼ softening), and on the sweep (DFS=Decreasing Frequency Sweep, IFS=Increasing Frequency Sweep). yi (mm)

IET

S

HL

FET

S

Peak of the transmissibility magnitude, |T|(a.u.) Preload (kg)

Boundary frequencies of the jump-phenomenon (Hz) Preload (kg) 16

26

36

46

56

16

26

36

46

56

IFS

0.1 0.15 0.2

40.92 (3.56) 36.08 (2.10) 34.60 (4.82)

34.07 (2.31) 31.46 (1.90) 29.02 (1.60)

30.81 (1.96) 28.88 (1.60) 26.98 (1.50)

29.31 (1.63) 27.53 (1.45) 26.25 (1.05)

28.62 (1.29) 27.02 (1.07) 25.83 (0.97)

4.44 (0.98) 4.24 (1.09) 3.75 (0.92)

5.48 (0.69) 4.56 (0.63) 3.65 (0.40)

5.52 (0.41) 4.56 (0.94) 3.73 (0.49)

5.40 (0.36) 4.25 (0.26) 4.00 (1.56)

5.13 (0.36) 3.99 (0.16) 3.60 (0.30)

DFS

0.1 0.15 0.2

38.92 (3.92) 33.44 (1.65) 30.68 (3.66)

30.14 (2.67) 27.30 (1.99) 24.71 (0.86)

27.03 (1.72) 24.45 (1.14) 22.93 (0.63)

25.46 (1.08) 23.50 (0.91) 22.28 (0.66)

24.98 (0.86) 23.25 (0.71) 22.24 (0.58)

4.69 (1.01) 4.54 (1.31) 3.84 (0.84)

6.57 (1.38) 6.12 (1.10) 5.52 (1.16)

7.14 (1.29) 7.06 (1.13) 6.59 (1.02)

7.09 (0.68) 6.83 (0.96) 6.63 (0.80)

6.72 (0.61) 6.26 (0.60) 6.06 (0.79)

IFS

0.1 0.15 0.2

37.79 38.71 (2.27) 37.25 (2.89)

30.40-33.19 32.27 (2.23) 30.99 (2.24)

n.p. 28.52-31.73 28.33 (1.27)

n.p. 26.53 25.83-28.47

n.p. n.p. n.p.

7.48 6.96 (2.01) 7.22 (2.17)

6.89-10.68 9.23 (3.20) 7.74 (2.36)

n.p. 4.66-13.72 7.11 (3.01)

n.p. 4.49 3.59-9.98

n.p. n.p. n.p.

DFS

0.1 0.15 0.2

34.25 35.38 (5.14) 32.64 (4.84)

28.65 (1.29) 27.08 (2.97) 26.03 (2.88)

n.p. 24.41-27.89 23.60 (1.51)

n.p. 23.58 21.62-24.47

n.p. n.p. n.p.

8.56 6.72 (2.01) 6.39 (2.17)

9.73 -10.68 8.60 (3.20) 7.28 (2.36)

n.p. 7.73-13.72 8.24 (3.01)

n.p. 8.13 7.85-9.98

n.p. n.p. n.p.

IFS

0.1 0.15 0.2

n.p. n.p. n.p.

30.2 29.13 25.79

26.69-29.87 25.01-27.53 23.85-25.96

27.12-31.92 22.67-24.06 21.95 (0.78)

28.69 (3.55) 26.24 (2.72) 24.81 (1.98)

n.p. n.p. n.p.

6.19 6.2 3.72

4.74-4.88 3.94-4.11 3.30-3.87

4.58-4.71 3.63-3.99 3.62 (0.90)

4.68 (0.57) 4.25 (0.86) 4.06 (0.96)

DFS

0.1 0.15 0.2

n.p. n.p. n.p.

28.99 26.29 23.23-24.6

24.61-26.46 22.52-24.16 21.33-22.72

24.63-26.9 25.37-29.46 25.42 (2.02)

25.84 (1.99) 23.28 (0.76) 21.98 (0.64)

n.p. n.p. n.p.

6.61 6.2 5.26-5.28

5.82-6.73 5.61-6.63 5.30-6.13

5.93-6.32 5.38-6.26 5.61 (0.59)

5.71 (0.59) 5.63 (0.53) 5.34 (0.59)

If the sample size was larger than three, mean and standard deviation () were reported. If the sample size was larger than one, maximum and minimum values were reported. n.p. (not present) was used in the case no test exhibited that specific dynamic response.

Fig. 6. Mean (symbol) and standard deviation (error bar) of the jump-phenomenon frequency vs. amplitude of the forced displacement (yi) for the applied preload in case of softening response (a) or hardening-like response (b), and vs. applied preload for yi of 0.1 mm, 0.15 mm and 0.2 mm in case of softening response (c) or hardening-like response (d). The plots distinguish between Increasing Frequency Sweep (blue) and Decreasing Frequency Sweep (red) and regard the Intact Endplates Tests. The black dashed line is the frequency computed by assuming a linear undamped system with stiffness such that its natural frequency is the average of the IFS and DFS mean frequencies for a preload of 16 kg. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the responses of the intact specimens and failed specimens without endplate failure did not differ significantly, i.e. JP-frequencies (p o0.012) and |T| (p o0.044). Furthermore, the difference between HLRs and SRs were investigated considering only the IETs.

In each test, the amplitude of the oscillation was larger in tension than in compression (p o0.01). The acceleration had an opposite trend (p o0.01) with larger values in compression compared to tension. The FET had lower |T| in the DFS compared to the IET (po 0.02) whereas in the IFS no significant difference

486

G. Marini et al. / Journal of Biomechanics 48 (2015) 479–488

was visible (p4 0.21). HLR had larger |T| compared to the SRs in both sweeps (unpaired test, p o0.007). In most specimens, the HLR-JP-frequencies did not differ significantly from the SR-JP-frequencies (unpaired t-test, p40.09) with the exception of the case 16 kg–0.15 mm. The FET had significant lower frequencies than the IET in the IFS for the cases with preload of 26 kg and 36 kg and yi 40.1 mm (po0.03). The range of IET-JP-frequency resulted between 22.24 Hz and 40.92 Hz whereas for the FET the range was between 21.1 Hz and 30 Hz (Table 2). A two factor repeated measures ANOVA did not highlight any interaction between preload and yi for both the JP-frequencies (p 40.06) and the |T| (p 40.08) for the IETs. However, both factors had a significant main effect (p o0.001) on the JP-frequencies and |T| with the exception of yi in the DFS for |T|. The factor’s main effect was investigated separately. 3.4. Amplitude of the stimuli Given a preload, increasing the amplitude of the stimuli induces a transition of the disc dynamic from a SR to a HLR (e.g. from 36–0.1 mm to 36–0.15 mm, Fig. 5). Increasing yi reduces the JP-frequency range of SRs (p o0.021, Fig. 6a and b and Table 2). The |T| decreases significantly with increasing yi for preload larger than 26 kg (p o0.01, Table 2). 3.5. Preload For a given amplitude of the forced displacement, increasing the preload induces a transition from HLR to SR (e.g. from 36–0.15 mm to 46–0.15 mm, Fig. 5). Increasing the preload mostly decreases the JP-frequencies (Fig. 6c and d and Table 2) (p o0.001). The |T| was maximum with preload of 26 kg in the IFS and 36 kg in the DFS (p o0.01, Table 2).

4. Discussion The presented results demonstrate the nonlinear dynamic response of the human lumbar IVD. Depending on the preload and amplitude of the stimulus, the specimen showed a dynamic behaviour similar to pure softening or hardening-like system (Fig. 4). Hardening systems have stiffness which increases with increasing deformation, whereas for softening systems the stiffness decreases with increasing deformation (Virgin, 2000). The performed quasi-static tests would classify the IVDs as a purely stiffening system with respect to the neutral condition (0 N force, 0 mm displacement). However, this is not the normal condition of the IVD. Assuming a certain preload is applied, the system will be in a different initial force equilibrium, with a higher initial stiffness compared to the zero-preload state. During oscillation, when the disc passes into compression, there would be an increase in stiffness. In extension, the system would have a decrease of stiffness until the neutral zone (area of minimal internal resistance, Panjabi, 1992) is passed, and then re-stiffening occurs. When both responses were seen for a given preload, the system had a SR at low amplitude of the stimuli (yi), which passed to a HLR when the amplitude was increased. Assuming the increase of the neutral zone between the two tests negligible, this entails that increasing the energy of the stimuli provided the system enough energy to neglect the initial decrease of stiffness in elongation. Increasing the preload requires more energy in the system to exhibit a HLR, hence this behaviour was seen only at large amplitudes. One could speculate that the different frequency between the JP in the IFS and in the DFS could be due to a change in stiffness

Fig. 7. Peaks of the displacement amplitude of the crossbar (x) vs input frequency for a C2C3 specimen (female, 65 yr) tested with two different loading protocols: (i) solid line) Increasing Frequency Sweep (IFS) followed by the Decreasing Frequency Sweep (DFS); (ii) dotted line) DFS was followed by the IFS. Preload was 16 kg and amplitude of the forced stimuli (yi) 0.15 mm. The plots distinguish between IFS (blue) and DFS (red). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

caused by creep. Initial quasi-static tests likely contributed strongly to a consistent state of the disc’s intrinsic properties. Nevertheless, to investigate this hypothesis, a specific test was performed with the standard sweep protocol (IFS-DFS). The test was then repeated with the sweep protocol inverted (DFS-IFS). Both tests exhibited similar response, with differences of frequency of circa 0.5 Hz, whereas the difference between the JP-frequencies of the IFS and DFS was of 5.7 Hz (Fig. 7). Moreover, by the comparison of test 1–2 of each preload-sequence, the frequency differed by circa 0.68 Hz (std 0.62 Hz) for 16 kg preload, which decreased for larger preload (o 0.14 Hz). These two observations would suggest that the reduction in stiffness due to the creep while the system oscillates in the nonlinear region is negligible compared to the intrinsic nonlinear dynamic. Displacement measurements from the camera system showed that each specimen lost disc height during the dynamic test. In some cases, this indicated that a failure of the specimen had occurred and correlated with the change in mechanical properties of the system (R¼ 0.81). Failure of the specimen did not always involve a visible endplate fracture. A visual inspection of the endplate alone is not suitable to assess any micro-damage of the subchondral bone, and failures which could also affect the integrity of the annulus. However, in the specimens where failure of the endplate was visible, a sudden loss of disc height was recorded (45%). After this height loss, the disc had purely a SR, which would confirm the hypothesis that the load transfer between nucleus pulposus and annulus fibrosus plays an important role in the disc dynamics by increasing the neutral zone after failure. This aspect also indicates that these peculiar dynamics are mainly related to the collagenous component of the annulus. In fact, although the mechanism of load transfer is compromised, the IVD still exhibits a nonlinear dynamic, thus implying that the residual mechanical properties offered by the structure, in particular by the collagenous component of the annulus, are sufficient to trigger these phenomena. The remaining failed specimens had a more gradual height reduction, which could correspond to a different mechanism of disc failure (i.e. progressive fissure or annulus-endplate junction

G. Marini et al. / Journal of Biomechanics 48 (2015) 479–488

487

Table 3 Range of resonance from the literature for the lumbar intervertebral disc. Study

Segment

Sex

Oscillating mass [kg]

Preload [N]

Testing method

Range of resonance [Hz]

Izambert et al. (2003)

6  L1L2 1  L2L3 1  L3L4

M, F

40 Applied by a lever arm

400 Applied by a lever arm

Discrete frequency (bath)

8.0–10.4

Kasra et al. (1992)

1  T12L1 6  L2-L3

M, F

40

30, 230, 390, 530, 680

Discrete frequency (air)

23.5–33

Present study

2  L1L2 9  L2L3 11  L4L5

M

16–56 With step of 10

160–550 With step of 100

Continuous frequency (air)

22.2–40.9

failure). Studies have proposed high rate compression loads as a mechanism of damage initiation and herniation in human lumbar disc (Wade et al., 2014). The present study showed that vibrations can induce load-shocks due to the change between two stable states of oscillation, which can in turn trigger the failure of the disc. The magnitude of this shock is dependent on loading condition of the disc (preload and stimuli) and is higher if hardening occurs (þ246%, Table 2). Whole body vibration has also been suggested as a factor to trigger disc degeneration by inducing fatigue failure of the endplates (Seidel et al., 2008), especially if bone mineral density is low. The mean age of the endplate-failed specimens (59.43 years) resulting significantly lower (p o0.02) than the intact specimens (39 years) indicates that this hypothesis could be valid and suggest care with respect to whole-body-vibrations in the workplace but also in the use of vibrating platform for training purpose in elderly (von Stengel et al., 2011). The present study was compared to the available literature related to the disc dynamics (Table 3). These studies (Kasra et al., 1992; Izambert et al., 2003) characterized the disc response at discrete frequencies, whereas our study controlled the input displacement by a continuous function, and swept from a low frequency to a high frequency and back down. The path of the stimuli has no influence on a linear system but for a nonlinear system is critical and it could be the reason any nonlinearity was not previously reported. Increasing the preload decreases the JP-frequency range, which would be expected in a linear system. However, the frequency change due to the mass increase did not follow a similar relationship (Fig. 6c and d), which was expected since the increase of preload translates into an increase of stiffness due to the nonlinear quasi-static response of the disc. For the observed nonlinear system, the principle of superposition does not hold. Furthermore, as an additional aspect of the nonlinear dynamic of the IVD, increases of the displacement amplitude decreased the frequency range of the JPs (Fig. 6a and b). For a linear system, the amplitude of the stimuli does not influence the resonant frequency of a mode. A similar characteristic was reported also in studies focused on WBV for different kind of exposure (Huang and Griffin, 2009; Matsumoto and Griffin, 2002). The presented study provided the first characterization of the nonlinear dynamics of the IVD. The response of the IVD was asymmetric and dependent on the applied preload and amplitude of the stimuli. Future development will concern the analysis and comparison of the mean point of oscillation in the nonlinear region which has been seen to be different from the stable branches.

Conflict of interest statement The authors have no conflict of interests to declare.

Acknowledgement Funding for this research project was provided by the European Union through a Marie Curie action (FPT7-PITN-GA-2009-238690SPINEFX). The authors would like to thank Dipl.-Ing. Kay Sellenschloh and Dipl.-Ing. Matthias Vollmer of the Institute of Biomechanics of the TUHH, and Prof. Dr. Remco Leine of the ETH Zurich, for their help and valuable suggestions. References Bovenzi, M., Hulschof, C.T.J., 1999. An updated review of epidemiologic studies on the relationship between exposure to whole-body vibration and low back pain (1986–1997). Int. Arch. Occup. Environ. Health 72, 351–365. Costi, J.J., Stokes, I.A., Gardner-Morse, M.G., Iatridis, J.C., 2008. 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 33, 1731–1738. Griffin, M.J., 1998. A comparison of standardized methods for predicting the hazards of whole-body vibration and repeated shocks. J. Sound Vib. 215, 883–914. Huang, Y., Griffin, M., 2009. Nonlinearity in apparent mass and transmissibility of the supine human body during vertical whole-body vibration. J. Sound Vib. 324, 429–452. Huber, G., Skrzypiec, D.M., Klein, A., Püschel, K., Morlock, M.M., 2010. High cycle fatigue behaviour of functional spinal units. Ind. Health 48, 550–556. ISO-2631-1, 1997. Mechanical Vibration and Shock—Evaluation of Human Exposure to Whole Body Vibration. Part 1: General Requirements. International Organization for Standardization, Switzerland. ISO-2631-5. 2005. Mechanical Vibration and Shock—Evaluation of Human Exposure to Whole Body Vibration. Part 5: Method for Evaluation of Vibration Containing Multiple Shocks. International Organization for Standardization, Switzerland. Izambert, O., Mitton, D., Thourot, M., Lavaste, F., 2003. Dynamic stiffness and damping of human intervertebral disc using axial oscillatory displacement under a free mass system. Eur. Spine J. 12, 562–566. Kasra, M., Shirazi-Adl, A., Drouin, G., 1992. Dynamics of human lumbar intervertebral joints. Experimental and finite-element investigations. Spine 17, 93–102. Kasra, M., Goel, V., Martin, J., Wang, S.-T., Choi, W., Buckwalter, J., 2003. Effect of dynamic hydrostatic pressure on rabbit intervertebral disc cells. J. Orthop. Res. 21, 597–603. Koeller, W., Muehlhaus, S., Meier, W., Hartmann, F., 1986. Biomechanical properties of human intervertebral discs subjected to axial dynamic compressioninfluence of age and degeneration. J. Biomech. 19, 807–816. Mansfield, N., 2005. Human Response to Vibration. CRC Press, Boca Raton, Florida, USA. Matsumoto, Y., Griffin, M.J., 2002. Non-linear characteristics in the dynamic responses of seated subjects exposed to vertical whole-body vibration. J. Biomech. Eng. 124, 527–532. Panjabi, M., 1992. The stabilizing system of the spine. Part II. Neutral zone and instability hypothesis. J. Spinal Disord. 5, 390–397. Sandover, J., 1998. The fatigue approach to vibration and health: is it a practical and viable way of predicting the effects on people? J. Sound Vib. 215, 699–721. Schwarze, S., Notbohm, G., Dupuis, H., Hartung, E., 1998. Dose–response relationships between whole-body vibration and lumbar disk disease-a field study on 388 drivers of different vehicles. J. Sound Vib. 215, 613–628. Seidel, H., Blüthner, R., Hinz, B., Schust, M., 1998. On the health risk of the lumbar spine due to whole-body vibration-Theoretical approach, experimental data and evaluation of whole-body vibration. J. Sound Vib. 215, 723–741. Seidel, H., Pöpplau, B.M., Morlock, M.M., Püschel, K., Huber, G., 2008. The size of lumbar vertebral endplate areas—prediction by anthropometric characteristics and significance for fatigue failure due to whole-body vibration. Int. J. Ind. Ergon. 38, 844–855. Thamsuwan, O., Blood, R.P., Ching, R.P., Boyle, L., Johnson, P.W., 2013. Whole body vibration exposures in bus drivers: a comparison between a high coach and a low-floor city bus. Int. J. Ind. Ergon. 43, 9–17.

488

G. Marini et al. / Journal of Biomechanics 48 (2015) 479–488

Thompson, J.P., Pearce, R.H., Schechter, M.T., Adams, M.E., Tsang, I.K., Bishop, P.B., 1990. Preliminary evaluation of a scheme for grading the gross morphology of the human intervertebral disc. Spine 15, 411–415. Thomson, W.T., 1981. Theory of Vibration with Applications. Prentice-Hall, Englewood Cliffs, NJ, USA. Urban, J.P.G., Stanton, S., Fairbank, J.C.T., 2004. Nutrition of the intervertebral disc. Spine 29, 2700–2709. Virgin, L.N., 2000. Introduction to Experimental Nonlinear Dynamics. Cambridge University Press, Cambridge, UK. von Stengel, S., Kemmler, W., Bebenek, M., Engelke, K., Kalender, W.A., 2011. Effects of whole-body vibration training on different devices on bone mineral density. Med. Sci. Sports Exerc. 43, 1071–1079.

Wade, K.R., Robertson, P.A., Thambyah, A., Broom, N.D., 2014. How healthy discs herniate: a biomechanical and microstructural study investigating the combined effects of compression rate and flexion. Spine 39, 1018–1028. Walsh, A.J.L., Lotz, J.C., 2004. Biological response of the intervertebral disc to dynamic loading. J. Biomech. 37, 329–337. White, A.A.I., Panjabi, M.M., 1978. Clinical Biomechanics of the Spine. Lippincott Williams & Wilkins, Philadelphia, Pennsylvania, USA.