A 1-D model of the nonlinear dynamics of the human lumbar intervertebral disc

Journal of Sound and Vibration 387 (2017) 194–206

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

A 1-D model of the 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, Zurich, Switzerland Institute of Biomechanics, TUHH Hamburg University of Technology, Hamburg, Germany Institute for Forensic Medicine, University Medical Center Hamburg-Eppendorf, Hamburg, Germany

a r t i c l e in f o

abstract

Article history: Received 7 July 2015 Received in revised form 30 August 2016 Accepted 19 September 2016 Handling Editor: J. Macdonald Available online 18 October 2016

Lumped parameter models of the spine have been developed to investigate its response to whole body vibration. However, these models assume the behaviour of the intervertebral disc to be linear-elastic. Recently, the authors have reported on the nonlinear dynamic behaviour of the human lumbar intervertebral disc. This response was shown to be dependent on the applied preload and amplitude of the stimuli. However, the mechanical properties of a standard linear elastic model are not dependent on the current deformation state of the system. The aim of this study was therefore to develop a model that is able to describe the axial, nonlinear quasi-static response and to predict the nonlinear dynamic characteristics of the disc. The ability to adapt the model to an individual disc's response was a specific focus of the study, with model validation performed against prior experimental data. The influence of the numerical parameters used in the simulations was investigated. The developed model exhibited an axial quasi-static and dynamic response, which agreed well with the corresponding experiments. However, the model needs further improvement to capture additional peculiar characteristics of the system dynamics, such as the change of mean point of oscillation exhibited by the specimens when oscillating in the region of nonlinear resonance. Reference time steps were identified for specific integration scheme. The study has demonstrated that taking into account the nonlinear-elastic behaviour typical of the intervertebral disc results in a predicted system oscillation much closer to the physiological response than that provided by linear-elastic models. For dynamic analysis, the use of standard linear-elastic models should be avoided, or restricted to study cases where the amplitude of the stimuli is relatively small. & 2016 Elsevier Ltd. All rights reserved.

Keywords: intervertebral disc nonlinear dynamic softening hardening jump-phenomenon

1. Introduction In industrialized societies, the body is often exposed to high frequency loads ( > 10 Hz) due to working conditions or transmitted from the environment. Although the effect of vibrations on spinal health has been extensively studied, the direct consequences of this exposure are not completely understood. In particular, many epidemiological studies have focused on the hypothesized correlation with low-back pain (LBP) [1–4], which is one of the primary causes for hospitalization and work loss for industrial workers [5,6]. Nevertheless, the findings have been quite contradictory [7–9] and do not allow

n

Corresponding author. E-mail address: [email protected] (G. Marini).

http://dx.doi.org/10.1016/j.jsv.2016.09.021 0022-460X/& 2016 Elsevier Ltd. All rights reserved.

G. Marini et al. / Journal of Sound and Vibration 387 (2017) 194–206

195

the establishment of a conclusive dose–response relationship. The intervertebral disc (IVD), with its particular anatomical and functional characteristics, has been often associated with LBP induced by whole body vibration. Lying between two adjacent vertebrae, the disc is also in direct contact with the spinal cord and the fascicles of branching nerves (nerve roots). Therefore, any alterations affecting the IVD (e.g. herniation) are highly likely to be perceived as focal pain and/or radiating pain and numbness of the lower limbs. After the third decade of life, the IVD typically enters a physiological process of degeneration, which results in altered structural, biochemical and functional properties [10]. Although in some cases the degenerated discs are asymptomatic [11], when LBP manifests advanced disc degeneration is usually visible. Studies have shown that this phenomenon could be also triggered, or enhanced, by the failure of the cranial and caudal boundary confinement provided by the vertebrae endplates [12,13]. High frequency loads of the spine have been postulated to accelerate the physiological disc degeneration by fatigue failure of the endplates [14,15]. To investigate the spine's response to different dynamic loads, numerical models are often used [16,17]. The ability to easily modify the model parameters to cover the physiological variance and to adapt the model to a specific scenario provides the opportunity to gain insight into the dose–response relationships without performing countless experiments. However, these models require a thorough validation to rely on the numerical outcome provided. In most of these studies a linear-elastic viscous model (LEVM) was assumed to describe the IVD response [18,16,19,20]. This results in a system whose mechanical properties are not dependent on its current deformation and, in particular, does not distinguish between the tension–compression asymmetry typical of the IVD axial response [21]. The influence of this aspect on the numerical outcome is unknown but could be highly relevant, since analytical models with analogous quasi-static mechanical response present peculiar dynamics such as various stable and unstable modes of oscillation [22]. In our recent experimental study, we have investigated the IVD dynamic response, focusing on the possible presence of nonlinearities with respect to different preloads and amplitudes of the mechanical stimuli [23]. The system response was characterized by sudden changes of oscillation amplitude at specific frequencies, which were dependent on the stimulus history (increasing or decreasing frequency sweep). This behaviour is typical of the dynamic response of systems with similar nonlinear quasi-static response. Moreover, the specific response of the disc was shown to be influenced by the preload and amplitude of the oscillation. Current models are not able to describe these complex dynamics because they mostly linearize the IVD elastic response. The aim of this study, therefore, was to develop an analytical model able to correctly describe the quasi-static response of the disc, including the tension-compression asymmetry and viscoelastic response. We hypothesize that such a model would exhibit the characteristic nonlinear dynamic response of the disc. The responses predicted by the nonlinear-elastic viscous model (NLEVM) were compared to corresponding dynamic experiments on individual IVDs for validation. The NLEVM was then compared to a LEVM, which assumed a local linearization of the quasi-static response. To provide reference parameters for the use of the NLEVM, the influence of the computational approach was also investigated.

2. Methods 2.1. Experiment The previously performed experiment has been described in Marini et al. [23] and served also to define the simulation model. Briefly, this experiment duplicates a base excitation model (Fig. 1a). In this configuration, the excitation is applied at the bottom of the system (y) and transmitted to the upper part (x), which is free to move. The specific mechanical properties of the system (i.e. elastic stiffness and damping) determine the body dynamics. Since the study focuses on the axial mode of oscillation, the lateral movement of the mass was constrained by lubricated Teflon holders. Thus, the setup resulted in a one degree-of-freedom system. The input displacement was prescribed by a sinusoidal function with constant amplitude (yi) but a sweeping frequency f(t):

y = yi sin (2πf (t )t )

(1)

The frequency was set to increase linearly from a low frequency ( flow ) to a high frequency boundary ( fhigh ) and back down. Therefore, we distinguished between the Increasing Frequency Sweep (IFS) and Decreasing Frequency Sweep (DFS):

⎧ for t < tturn(IFS) ⎪ flow + frate t f (t ) = ⎨ ⎪ f ⎩ high − frate t for t > tturn(DFS)

with tturn =

fhigh − flow frate

(2)

The frequency sweep rate ( frate ) was set to 0.1 Hz/s, after verifying that it was low enough to demonstrate a transient steady response of the system. Weights (mi) were rigidly attached cranially to mimic the static preload (preloadi = g ·mi with g the gravity) induced by the upper part of the body. Before the dynamic test, a quasi-static tension–compression test was performed to precondition the specimen and to obtain the nonlinear quasi-static response (5 cycles, tension 200 N, compression 1000 N, and displacement rate 0.09 mm/s). Frequency boundaries flow and fhigh were defined by the analysis of the last cycle of the quasi-static test, which was force-

196

G. Marini et al. / Journal of Sound and Vibration 387 (2017) 194–206

Fig. 1. 1-dof base-excitation model used as reference setup for the experiment (a); typical force–displacement curve obtained by force-averaging the last cycle of the tension–compression test performed before and after the dynamic characterization (b). The point pi represents the reference position around which the system oscillates when a specific mass (mi) is applied (preloadi = g ·mi , with g the gravity).

averaged. For small stimuli the dynamic of nonlinear systems can be approximated with a linear model, hence, assuming a linear undamped model, the ideal resonant frequency ( fid ) would be given by:

fid =

1 k preload/preload 2π

(3)

The initial frequency range was defined by fid ± 6 Hz . The dynamic test was repeated with different weights from 16 kg to 56 kg, in steps of 10 kg. For each preload, three different yi were used (0.1 mm, 0.15 mm, 0.2 mm). The acceleration and displacement of the upper and lower end of the system were recorded by accelerometers and a camera system. The specimens (N ¼22, Table 1) exhibited abrupt changes of the oscillation amplitude (jump phenomenon, JP) at specific frequencies, which were different for the IFS and DFS. These are characteristic of the dynamics of nonlinear systems with quasi-static response similar to the IVD axial response. Two responses were observed in the tests: a pure softening response (SR), where the oscillation amplitude increased for decreasing frequency of the stimuli, and a hardening-like response (HLR), where the oscillation amplitude grew when increasing the frequency of the stimuli [23]. Table 1 Characteristics of the tested specimens, with the classification based on the comparison of the area of the hysteresis cycle and change in maximum displacement in tension and compression in quasi-static tensile tests performed before and after the dynamic tests (specimen status, I¼ intact, F ¼ failed, and FE ¼failed-endplate) and on the dynamic response exhibited in the test sequence (HLR ¼hardening-like response and SR ¼ softening response). The scaling factor (α) is used to describe the change of viscous properties due to the high strain rate ( > 20 Hz) . No.

Age (yrs.)

Disc level

Disc grade

Volume (mm3)

Height (mm)

Status

Dynamic response

Scaling factor α (10  3)

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 9764 11 367 14 319 17 208 15 052 16 717 22 739 20 974 13 284 9839 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

HLR þ SR SR HLR þ SR HLR þ SR HLR þ SR HLR þ SR HLR þ SR HLR þ SR HLR þ SR HLR þ SR SR HLR þ SR SR HLR þ SR HLR þ SR HLR þ SR HLR þ SR SR HLR þ SR SR HLR þ SR HLR þ SR

2.45 2.71 3.14 2.32 3.33 3.92 4.01 2.05 2.56 3.81 3.63 3.01 2.01 3.23 2.42 3.17 2.03 1.95 2.60 3.51 2.84 2.32

G. Marini et al. / Journal of Sound and Vibration 387 (2017) 194–206

197

2.2. Disc model The reference uniaxial elastic response of the disc was described by a polynomial function (Eq. (4)), which distinguished between the two tensile loading conditions: tension and compression.

⎧ 3 5 ⎪k u + k u for u > 0 mm tension 3t 5t S(u) = S0 + k u + ⎨ ⎪ 3 5 ⎩ k3cu + k5cu for u < 0 mm compression

(4)

where u = x − y − d0 , d0 and S0 are respectively a displacement and force offset, k is the stiffness at the point (d0, S0), the point of minimum stiffness, and k3t, k3c, k5t, k5c are the parameters which characterize the nonlinear response in tension and compression. This formulation allows the description of the two axial states independently, while preserving the continuity of the function at the point (d0, S0). When the dynamic response is investigated numerically, spurious oscillations can develop by the system passing between two states with different stiffness. The dissipative properties of the disc were described by an additional term, D:

D(u̇) =

1 2

( μ ( 1 − erf(τu )) + μ ( 1 + erf(τu ))) v

c

v

(5)

t

v d0v , d0

where u̇ is the relative velocity between the base and the top (x ̇ − ẏ ) and uv = x − y − is an offset parameter related to the change of the dissipation factor from tension (μt) to compression (μc), and viceversa, erf() is the error function and τ an additional parameter modelling the transition between tension and compression state of the damping factor. For further details, refer to Appendix A. The parameters were characterized in a two-step optimization procedure. In the first step, the parameters characterizing the elastic response were obtained, Eq. (4). In the next step, the fitting of the force–displacement curve (not force averaged) was repeated using Eq. (5) and the remaining parameters were determined. 2.3. Base excitation model The free body diagram of mi for the experimental setup (Fig. 1a) is described by the following equation:

mi x¨ + D(ẋ − ẏ ) + S(x − y) = 0

(6)

where D(x ̇ − ẏ ) is the dissipative term which is a function of the relative velocity between the base and mi (x ̇ − ẏ ) while S (x − y ) is the elastic term and it is related to the relative displacement between mi and the base (x − y ). Eq. (6) models a specific equilibrium condition defined by the applied mi. In the tension–compression tests, the reference force–displacement response was obtained without weights (mi ¼ 0 kg). This corresponds to the reference condition about (0 N, 0 mm) (p0, Fig. 1b). When a preload is applied (pi, Fig. 1b), the system shifts to a new equilibrium condition due to the disc compression (preloadi , disp0i ). Since the mechanical properties of the disc depend on the deformation, the stiffness at pi is different than at p0. To consider the effect of the mass change on the model response, Eq. (6) was modified as follows:

mi x¨ + D(ẋ − ẏ ) + S(x − y + disp0i ) + preloadi = 0

(7)

2.4. Numerical analysis The model (Eq. (7)) was implemented in Simulink® 8.3 (The Mathworks, Natick, Massachusetts). An explicit Dormand– Prince solver (ode8) with fixed time increment (0.0005 s) was used to numerically solve the equation. The displacement recorded during the experiments by the testing machine (MTS® Bionix, Minneapolis, Minnesota) was used as input for the model. The peak values of the mass oscillation (maximum and minimum) within one y-cycle were then plotted as functions of the frequency and compared to the experimental results [23]. Assuming mi of 36 kg, the influence of linear- vs. nonlinear-elastic response was investigated with model No. 3 (Table 1). The stiffness of the LEVM was set equal to the stiffness of the NLEVM at the reference point defined by the preloadi . A continuous sinusoidal sweep ( flow = 1 Hz, fhigh = 60 Hz ) with constant amplitude was applied as input. Both the IFS and DFS were simulated. The influence of the numerical integration on the response of the NLEVM was examined with model No. 1 (Table 1). The loading scenario replicates the preloads and yi applied in the experiments. Beginning with 0.01 s, the time step was progressively reduced by half to a minimum of 0.0001 s. The analysis was repeated with different integration schemes: Heun's method (ode2, second-order accuracy), fourth-order Runge–Kutta (ode4, fourth-order accuracy) and Dormand–Price RK8 (ode8, eighth-order accuracy). Furthermore, to investigate the influence of the numerical approach, specifically the path followed by the stimuli on the response predicted by the model, the continuous sweep was compared to a discrete excitation. Both the IFS and the DFS were investigated. In this case the frequencies of the spectrum were examined separately with steps of 0.2 Hz. For each frequency, a simulation was run with an input of 45 cycles. These were organized into three consecutive steps: (i) 15 oscillations where the input frequency was set to 1 Hz for the forward scan (IFS) or 60 Hz for the backward scan (DFS); (ii) 15

198

G. Marini et al. / Journal of Sound and Vibration 387 (2017) 194–206

Fig. 2. Comparison of the reaction force predicted by the NLEVM (dashed grey), after optimization of the parameters, with the quasi-static tension– compression tests (black) performed before (a) and after (b) the dynamic tests.

oscillations with input frequency set to the desired discrete test value to adapt the system to the new condition; and (iii) 15 oscillations where the maximum and minimum amplitude of oscillation were recorded [24]. 2.5. Data analysis Specimen condition strongly influenced the dynamic response measured in the experiments; therefore the specimens were classified into intact (I) or failed (F) groups, depending on the quasi-static response before and after the dynamic test [23]. The failed specimens were further distinguished based on the presence of a macroscopic endplate fracture, which was also recognizable by a sudden loss of the specimen height during the dynamic testing sequence (failed-endplate, FE). Depending on the normality of the data (Shapiro–Wilk), parametric tests (t-test and ANOVA) or non-parametric tests (Wilcoxon and Kruskall–Wallis) were used to compare the groups. In some cases, a logarithmic transformation was applied to ensure equal variance (Levene and Brown–Forsythe). Spearman's correlation was used to assess the model prediction.

3. Results 3.1. Quasi-static response The response predicted by the NLEVM, after fitting of the parameters for the elastic and dissipative response, is shown in Fig. 2. The results of the optimization are reported in Table 2 for the quasi-static test performed before the dynamic characterization and in Table 3 for the tests conducted at the end of the experiment. After transformation of the parameters to respect unit requirements, a further scaling of μc, μt was required for each specimen (Table 1) to consider the reduction of the viscous properties of the disc with increasing strain rate ([25]). The scaling factor was obtained through optimization with respect to the dynamic tests. In general, the model was able to correctly describe the experimental quasi-static force–displacement curves measured before (R > 0.996, p < 0.001) and after (R > 0.982, p < 0.001) the dynamic characterization. Furthermore, significant differences in the correlation coefficient were not found when comparing the model prediction for tension–compression (unpaired K–W, p > 0.064 ), or between groups (unpaired K–W, p > 0.07). Significant differences in the parameters k, k5c and k5t were found between the I and FE specimens for the tests performed before the dynamic characterization (unpaired, p < 0.06 ). After the dynamic test, the parameters of the I and FE specimens showed differences, except for f0, k5t, k3t, k3c, μt (unpaired, p < 0.029). Additionally, significant changes were found between the F and FE for k, k5c and d0v (unpaired, p < 0.042). After the dynamic testing, all the specimens showed a significant change (paired, p < 0.001) in the parameters regulating the system dissipation (μc, μt). Except for k3t, the defining parameters for the I group were not influenced significantly by the dynamic test. The F and FE groups showed significant differences in k, k3c and d0 (paired, p < 0.031) before and after dynamic testing. Furthermore, the F group showed significant differences in k3t and k5c before and after dynamic testing (paired, p < 0.032). 3.2. Dynamic response: model vs. experiment The NLEVM showed in many cases (67.7 percent) a dynamic response, HLR or SR, which was similar (R > 0.5) to the corresponding experiments especially in the DFS (mean R ¼0.66) and for preload larger than 16 kg. However, circa 13.7 percent of the simulations produced an erroneous response, HLR instead of a SR. The model of the specimens that exhibited only a SR in the experiments consistently predicted only this kind of response. Fig. 3 shows the response predicted by the

G. Marini et al. / Journal of Sound and Vibration 387 (2017) 194–206

199

Table 2 Model parameters describing the elastic and viscous response of the quasi-static tension–compression test performed before the dynamic characterization. Data were grouped according to specimen status after the dynamic characterization (specimen status, I¼ intact, F¼ failed, and FE¼ failed-endplate). Status

No.

Elastic parameters

Viscous parameters 3

5

5

d0 (mm)

f0 (N)

k (N/mm)

k3c (N/mm )

k3t (N/mm3)

k5c (N/mm )

k5t (N/mm )

d0v (mm)

τ

μc (N s/mm)

μt (N s/mm)

I

2 3 5 7 9 11 16 18

0.26 0.17 0.1 0.11 0.19 0.17 0.08 0.17

1 9 22 14 20 11 14 11

9.67 35.52 210.1 126.83 97.79 44.22 165.51 68.24

232.9 586.95 792.36 1422.2 667.17 277.09 1757.9 601.31

41.35 345.18 341.32 516.59 214.32 120.84 904.7 12.94

215.05 289.43 267.29 0 237.58 63.38 0 250.72

0 0 0.01 0 81.71 0 0 0.96

 0.3  0.21  0.14  0.4  0.35  0.3 0.01  0.17

2.92 2.99 2.66 7.97 3.68 3.34 1.78 2.56

429.69 284.68 327.48 470.15 525.29 468.76 504.25 551.41

78.93 52.99 40.32 139.89 116 105.21 108.08 128.88

F

1 4 6 8 10 15 22

0.15 0.01 0.05 0.05 0.03 0.31 0.26

11 13 34 23 30 83 45

24.38 458.62 601.82 244.83 595.24 147 92.99

2201.08 1961.68 3961.17 1707.46 1112.95 1067.82 826.3

28.88 472.96 3055.95 368.17 1878.99 186.38 175.72

0 0 0 0 0 0 104.54

121.73 1511.37 147.87 2716.95 52.11 0.09 0

 0.15  0.2  0.15  0.1  0.16  0.2  0.26

4.95 3.81 8.46 2.11 2.72 1.52 2.88

378.67 347.15 377.88 351.85 343.19 651.07 556.89

64.14 110.13 247.76 32.86 63.91 67.55 117.57

FE

12 13 14 17 19 20 21

0.17 0.07 0.07 0.08 0.05 0.24 0.05

152 31 41 21 51 12 36

629.4 241.43 392.19 181.53 349.71 36.98 445.65

171.18 351.3 1181.49 1797.1 3095.29 883.7 2004.42

0 0 651.19 37.03 2956.43 195.36 958.17

152.69 32.9 0 0 0 171.49 0

0.24 838.5 2161.96 0.01 172.22 0 311.29

 0.5  0.2  0.1  0.08  0.18  0.35  0.2

21.05 2.81 1.87 1.79 1.75 10.3 3.76

391.85 580.88 343.59 440.51 598.89 497.21 384.7

291.35 294.68 28.38 102.96 112.19 193.97 103.06

Table 3 Model parameters describing the elastic and viscous response of the quasi-static tension–compression test performed after the dynamic characterization. Data were grouped according to specimen status after the dynamic characterization (specimen status, I¼ intact, F¼ failed, and FE¼ failed-endplate). Status

No.

Elastic parameters

Viscous parameters

d0 (mm)

f0 (N)

k (N/mm)

k3c (N/mm3)

k3t (N/mm3)

k5c (N/mm5)

k5t (N/mm5)

d0v (mm)

τ

μc (N s/mm)

μt (N s/mm)

I

2 3 5 7 9 11 16 18

0.46 0.23 0.2 0.31 0.02 0.21 0.14 0.19

11 10 15 31 1 8 27 11

21.96 8.95 61.85 89.23 91.9 18.84 101.99 65.33

60.33 637.29 308.61 441.93 883.02 366.27 1937.85 574.33

8.89 137.63 62.69 0 5.59 5.21 19.63 0

66.59 148.71 355.53 387.95 151.25 79.35 42.62 343.16

2.21 0 18.1 32.72 61.53 17.25 0 1.88

 0.5  0.23  0.2  0.4  0.2  0.38  0.13  0.19

2.58 4.97 3.42 2.3 2.55 3.95 5 3.22

995.54 659.55 300.3 576.56 541.76 583.67 447.49 605.81

151.96 94.28 66.25 127.5 128.12 88.46 191.21 148.85

F

1 4 6 8 10 15 22

1 0.22 0.19 0.24 0.23 0.65 0.73

14 30 3 25 35 93 52

16.46 146.12 140.73 119.46 175.81 111.02 63.71

0 217.16 1804.92 233.86 592.25 8.86 30.18

0.4 110.84 57.07 0 125.27 12.67 13.68

35.69 253.51 482.26 77.13 51.43 97.05 69.68

5.48 280.13 0 116.22 127.5 85.36 2.2

 0.05  0.22  0.19  0.24  0.2 0.2 0.05

2.54 3.76 8.66 3.04 3.06 1.45 1.84

1106.53 746.82 644.88 1009.38 671.69 753.74 925.03

218.32 226.75 340.11 265.68 212.31 56.6 166.69

FE

12 13 14 17 19 20 21

0.8  0.17 0.86 0.45 0.48 1.24 1.87

13  22 6 11 9 10 15

35.7 60.76 18.72 27.93 42.28 19.54 1.08

0 139.8 0 89.43 0 0 0

47.61 44.06 0 0 0 0 2.88

3.61 36.6 6.86 38.11 32.55 3.7 2.7

0 13.55 1.04 5.13 18.36 4.12 0.29

 0.9  0.7  0.65  0.38  0.48  0.45  0.2

1.6 3.18 1.61 2.08 1.83 1.3 1.27

920 736.2 1013.92 791.55 786.34 885.94 1609.8

215.89 233.05 136.45 130.67 148.2 99.17 273.88

NLEVM for a L1L2 specimen, compared to the response observed in the experiment [23]. Depending on the part of the response investigated (IFS or DFS), the JP-frequencies predicted by the NLEVM were sometimes different. The JP-frequencies in the IFS were lower than in the experiment (average difference,  8.78 percent) whereas in the DFS the differences were smaller (average difference,  1.4 percent). No significant effect was found when comparing the model prediction vs. the specimen status for the I and F groups

200

G. Marini et al. / Journal of Sound and Vibration 387 (2017) 194–206

Fig. 3. Peaks of the mass displacement (x) vs. input frequency for a preload (mi) of 26 kg and amplitude of the forced displacement (yi) of 0.1 mm (a, d), 0.15 mm (b, e) and 0.2 mm (c, f) predicted by the NLEVM (a–c) and related experiment (d–f) for an L1L2 specimen. The displacement recorded by the testing machine was used as input for the numerical model [23]. The plots distinguish between Increasing Frequency Sweep (black) and Decreasing Frequency Sweep (grey).

Fig. 4. Peaks of the mass displacement (x) vs. input frequency for a preload of 36 kg and different amplitude of the forced displacement (yi) predicted by a LEVM (a) and a NLEVM (b). The models represent an intact specimen (No 3, Table 1). The plots distinguish between Increasing Frequency Sweep (black) and Decreasing Frequency Sweep (grey).

G. Marini et al. / Journal of Sound and Vibration 387 (2017) 194–206

201

(unpaired, p > 0.08). For the FE-specimens, the response predicted after the endplates failure showed a lower correlation with the experiments (R < 0.73) than before fracture occurred. 3.3. Dynamic response: linear model vs. nonlinear In Fig. 4, the response predicted by LEVM is compared to the dynamics predicted by the NLEVM for different magnitude of the stimuli. In general, the LEVM predicted a symmetrical behaviour in tension–compression, with larger amplitude of oscillation compared to the NLEVM. The response of the NLEVM differed significantly from the dynamic response of the LEVM, depending on the yi magnitude. However, these differences decreased for yi < 0.1 mm(R > 0.85, p < 0.001). For the LEVM, the peak of the oscillation was not affected by changes of yi. Further, the dynamic force (mi ·x¨ ) vs. the IVD compression (x − y ) at different time increments of the experimental test and predicted by the LEVM and by the NLEVM are compared in Fig. 5. The response of the NLEVM (Fig. 5e and f) is closer to the real system response (Fig. 5a and b) than the LEVM (Fig. 5c and d). However, the NLEVM is not able to completely capture the change of the mean oscillation point (average of the 0 N dynamic force within one half cycle) observed in the experiments. This characteristic was observed when the IVD was oscillating in the area of nonlinear resonance (between JPs). Out of this region, the mean point of oscillation returned to 0 mm. 3.4. Dynamic response: numerical integration approach The maximum time step to ensure the convergence of the simulation was dependent on the integration technique. The model response was not affected by changes of the time step for time increments smaller than 0.0005 s for the ode2 and ode4 (R > 0.992, p < 0.001) and 0.001 s for ode8 (R > 0.997, p < 0.001). For larger increments, the simulation did not converge or showed a dynamic response different from expected. In Fig. 6, the effect of a continuous (Fig. 6a–c) and a discrete variation of the frequency (Fig. 6a–c) of the stimuli on the dynamic response predicted by the NLEVM for Specimen No. 1 (Table 1) are shown. The numerical approach used to characterize the model response affected the predicted dynamic response. While at low frequencies the first nonlinear mode of oscillation had similar characteristics (JP frequency range, dynamic mode) with both approaches (R > 0.95, p < 0.001), at higher frequencies the system response was sometimes different (Fig. 6, yi ¼0.2 mm). In particular, the difference between the JP-frequencies in the IFS and in the DFS reduced and in some cases was almost null. This was mostly due to changes in the JP-frequency in the DFS (average error discrete vs. continuous 9.31 percent) rather than in the IFS, where the JP-frequency did not change significantly (average error discrete vs. continuous  1.91 percent). The amplitude of oscillation at the JP showed significant differences only for the IFS (paired t test, p < 0.001). In particular, the discrete approach predicted larger amplitude of oscillation especially in elongation (average difference þ 23.6 percent).

4. Discussion 4.1. Modelling parameters The response predicted by the NLEVM was in good agreement with the experimental measurements for specimens with no endplate failures. Although the modelling approach scaled the parameters with respect to the height loss during the test sequence, this approach would not be adequate to describe the change of mechanical properties when a system failure (fracture) occurs. The JP-frequencies predicted by the NLEVM were generally lower than those observed in the experiments, which also implies that some characteristics of the system response were not yet perfectly described. However, the dynamic responses predicted by the NLEVM were mostly similar to those observed in the experiments; of particular note, the models correctly predicted only a SR for specimens that showed only a SR in the experiments. Nevertheless, a model that considers the nonlinear-elastic response of the IVD requires a careful definition of the simulation. Since the numerical integration requires the discretization of the nonlinear ordinary differential equations, the accuracy and stability of the solution depends on the chosen solver and time step [26]. The simulations showed that also the method to test the numerical response (continuous or discrete sweep) influenced the model prediction (Fig. 6). Although the main response at low frequencies had similar characteristics, after the system crossed the first JP the predicted dynamics were different, especially for large excitation (yi ¼0.2 mm). 4.2. Linear vs. nonlinear model The comparison of the dynamic response predicted by the LEVM and by the NLEVM, showed that a model that uses a local linearization of the IVD response is not able to describe the complex dynamics of the physiological system. While for small amplitudes of the excitation (yi < 0.1 mm) the NLEVM predicted a symmetrical response similar to the LEVM (Fig. 4), for larger yi the dynamics of the two models were different (R < 0.65, p < 0.001). Li et al. [27], in their study focused on the limitations of the standard linear solid model of the IVD for low-frequency vibration ( < 1 Hz) had a similar observation. Although the NLEVM showed good agreement with the experiments, it presented also some shortcomings for small

202

G. Marini et al. / Journal of Sound and Vibration 387 (2017) 194–206

Fig. 5. x − y displacement vs. time response of specimen No. 1 (Table 1) for a preload (mi) of 26 kg and forcing amplitude (yi) of 0.15 mm from the related experiment (a) and predicted by the LEVM (c) and the NLEVM (e). Below each plot the dynamic force (mi ·x¨ ) vs. relative displacement (x − y ) has been plotted (b, d, f). From left to right, the plots show the disc response at different time steps identified by the black line in the time response plots (a, c, e). The mean point of oscillation (light grey) was computed by averaging the x − y at 0 N dynamic force within half cycle.

G. Marini et al. / Journal of Sound and Vibration 387 (2017) 194–206

203

Fig. 6. Peaks of the mass displacement (x) vs. input frequency for the model of specimen No. 1 (Table 1), for a preload (mi) of 16 kg and amplitudes of the forced displacement (yi) of 0.1 mm (a, d), 0.15 mm (b, e) and 0.2 mm (c, f), predicted by the NLEVM. The effect of a continuous (a–c) and a discrete (d–f) sweep of the stimuli are compared. The plots distinguish between Increasing Frequency Sweep (black) and Decreasing Frequency Sweep (dashed grey).

amplitude of oscillation. In some cases, the model predicted an IFS which differed significantly (up to 100 percent of oscillation amplitude) from the related experiment (Fig. 3a vs. d). However, the ability to describe the sudden change of amplitude of oscillation seen in the experiments makes the NLEVM more suitable, compared to the LEVM, to describe the IVD response for different dynamic scenarios. 4.3. Viscous response The response of the NLEVM predicted also lower JP-frequencies than observed in the experiment, especially in the IFS. The differences could be due to the fact that the response of the system in elongation was characterized only up to 200 N, whereas the IVD's dynamic response showed substantially larger forces. This implies also that the characterization of the damping properties of the IVD at high frequency should be more carefully assessed. While it was shown that considering the nonlinear elasticity of the IVD drastically changes the predicted dynamics, comparison with the experiments suggests also that when entering in the region of nonlinear resonance, the damping factor also changes (Fig. 5b vs. f). In previous studies in which the effect of strain rate on the energy dissipation of the IVD has been investigated [28,25], a small decrease of the energy loss with increasing loading frequency has been reported. However, the frequency range tested ( ≤ 10 Hz) was lower than the region of nonlinear resonance measured in the experiments [23]. The low phase angle (12°) reported in these studies [25] indicates that the IVD mainly behaves as an elastic structure in that frequency range ( ≤ 10 Hz). Furthermore, Jamison et al. [29] investigated the effect of impact trauma (with frequency components ≤12.5 Hz) on lumbar IVDs and found a distinct change of the system dissipation for impact durations between 80 ms and 160 ms. Although the change was relatively small (3–7 percent), it could indicate the switching to a different mechanism of energy dissipation. In the present study, the viscous parameters (μc, μt) obtained from the analysis of the quasi-static tests were scaled by a factor α. Considering the frequency range of our experiments ( > 20 Hz) and the magnitude of the scaling (10  3 order), there is some evidence of a strain rate effect on the dissipative properties of the disc. Interestingly, Iatridis et al. [30] observed that the viscous effect linearly increases with frequency ( < 31.8 Hz) in shear stress–relaxation tests of isolated nucleus pulposus specimens. This opposite trend to that reported by Costi et al. [25] and Koeller et al. [28] suggests that the AF could also have

204

G. Marini et al. / Journal of Sound and Vibration 387 (2017) 194–206

a significant function in the damping properties of the disc. A recent numerical study showed that the collagen fibre reinforcement of the IVD can contribute to a reduction of the internal pressure provoked by impact (10 ms in duration) by changing the stress propagation in the IVD [31]. Since the response of the collagen fibres is strongly nonlinear, the preload applied on the disc has a strong effect on the system response [21]. Therefore, the effect of the preload and strain rate on the hysteresis cycle of the IVD should be further investigated to better understand the mechanism of energy dissipation of the IVD. 4.4. Clinical and biomechanical aspects Even though this modelling approach showed good agreement with the experiments, some important characteristics were still not captured which may be relevant for further study of the response of the IVD to vibration. The change of the mean point of oscillation with entering the region of nonlinear resonance, which was observed in the experiments, was not fully predicted by the numerical models. Finite element studies of impact loading on the intervertebral disc predicted a mean point of oscillation of the system different from the reference position prior to impact [31]. This characteristic could be relevant, for example to elucidate why subjects who suffer from LBP have reported a reduction of the pain perception while training on vibrating platforms [32,33]. These devices operate in the frequency regime of the nonlinear resonance of the lumbar IVD (15 Hz–30 Hz). Therefore, by stimulating the system in this frequency domain, the height of the IVDs could slightly increase and thus release compression of the nerve, with a short-term relief of the pain. Alternatively, it could induce indirectly a stiffening of the surrounding muscles by stimulating their response through a cyclic elongation (i.e. stretch reflex). This would result in a stabilization of the motion units, which could reduce the occurrence of pain due to the relative movements of the segments. Previous studies have reported increased activity of the erector spinae muscle in subjects therapeutically exposed to whole body vibration [34]. However, the benefits of using vibration as clinical treatment must consider also the bone quality of the adjacent vertebrae, since it has been already observed the potential for spinal fracture during vibrational loading [35,23]. From the biomechanics perspective, the study confirmed that the nonlinear dynamic phenomena observed in the experiment are mainly due to the nonlinear quasi-static response of the disc, and not caused by the testing equipment or by the specimen constraints. The possibility to replicate similar phenomena with a relatively simple numerical model is a powerful tool to design experiments to further investigate the disc behaviour in other dynamic loading scenarios. Furthermore, analogous approaches can be used to investigate potential nonlinear dynamic phenomena also in other biological systems where the quasi-static response is well known to be nonlinear (e.g. cartilaginous tissues).

5. Conclusions The presented study demonstrated that considering the nonlinear asymmetric tension–compression response of the IVD in a numerical model is essential to accurately capture its complex quasi-static and dynamic axial response. Depending on the amplitude of the stimuli, a standard linear-elastic model predicted a different dynamic response to a nonlinear-elastic model. Thus, the linearization of the disc quasi-static response should be avoided or restricted to study cases where the amplitude of the mechanical stimuli is relatively small.

Acknowledgement Funding for this research project was provided by the European Union through a Marie Curie action (FPT7-PITN-GA2009-238690-SPINEFX).

Appendix A. Sigmoid switch function The error function erf() belongs to the class of sigmoid functions [36]. These are functions that describe the transition between two asymptotic states by an S-curve (Fig. A1a):

erf(τu) =

1 2π

∫0

τu

2

e−t dt

(A.1)

erf() is an odd function (erf( − u) = − erf(u) and presents no singularities, with the exception of infinity. These properties make the function ideal to describe the change between two states continuously, with the derivatives also being continuous. The rate of change between the two asymptotic conditions can be adjusted by the parameter τ (Fig. A1a). By shifting erf() by one towards the positive axis and multiplying by half of a target value (μt), the function is similar to a continuous switch function (0, μt) centered about u¼0:

G. Marini et al. / Journal of Sound and Vibration 387 (2017) 194–206

205

Fig. A1. Effect of the parameter τ on the error function (a); typical transition between two parameters described by Eq. (A.4) (b); example of application of Eq. (A.4) to model a system characterized by two different sets of parameters and responding to a linear relationship with respect to the independent variable, g (u)·u(c ) .

⎛μ ⎞ ⎧ μ for u → + ∞ lim⎜ t (1 + erf(τu))⎟ = ⎨ t u ⎝ 2 ⎠ ⎩ 0 for u → − ∞

(A.2)

Assuming the application requires the function to be equal to a target value (μc) only for negative u, this can be achieved by substituting u with −u:

⎛μ ⎞ ⎛μ ⎞ ⎧ 0 for u → + ∞ lim⎜ c (1 + erf( − τu))⎟ = lim⎜ c (1 − erf(τu))⎟ = ⎨ u ⎝ 2 u ⎝ 2 ⎠ ⎠ ⎩ μc for u → − ∞ Finally, to pass continuously from

g (u) =

μc 2

(1 − erf(τu)) +

μt 2

(A.3)

μc for negative u to μt for positive u, the following relationship can be used (Fig. A1b):

(1 − erf(τu))

(A.4)

g(u) can therefore be used to describe the transit of a target parameter which allows the description of a specific response distinguishing between the positive and negative branch of the independent variable (Fig. A1c).

References [1] J.L. Kelsey, An epidemiological study of the relationship between occupations and acute herniated lumbar intervertebral discs, International Journal of Epidemiology 4 (1975) 197–205. [2] B. Netterstrøm, K. Juel, Low back trouble among urban bus drivers in Denmark, Scandinavian Journal of Public Health 17 (1989) 203–206. [3] O. Thamsuwan, R.P. Blood, R.P. Ching, L. Boyle, P.W. Johnson, Whole body vibration exposures in bus drivers: a comparison between a high-floor coach and a low-floor city bus, International Journal of Industrial Ergonomics 43 (2013) 9–17. [4] D.G. Wilder, B.B. Woodworth, J.W. Frymoyer, M.H. Pope, Vibration and the human spine, Spine 7 (1982) 243–254. [5] S. Dagenais, J. Caro, S. Haldeman, A systematic review of low back pain cost of illness studies in the United States and internationally, The Spine Journal 8 (2008) 8–20. [6] P.P.F.M. Kuijer, H.F. van der Molen, A. Schop, F. Moeijes, M.H.W. Frings-Dresen, C.T.J. Hulshof, Annual incidence of non-specific low back pain as an occupational disease attributed to whole-body vibration according to the National Dutch Register 2005–2012, Ergonomics 58 (2015) 1232–1238. [7] K.T. Palmer, C.E. Harris, E.C. Harris, M.J. Griffin, J. Bennett, I. Reading, M. Sampson, D. Coggon, Case-control study of low-back pain referred for magnetic resonance imaging, with special focus on whole-body vibration, Scandinavian Journal of Work, Environment & health 34 (2008) 364–373. [8] A. Jensen, L. Kaerlev, F. Tchsen, H. Hannerz, S. Dahl, P.S. Nielsen, J. Olsen, Locomotor diseases among male long-haul truck drivers and other professional drivers, International Archives of Occupational and Environmental Health 81 (2008) 821–827. [9] S. Lings, C. Leboeuf-Yde, Whole body vibration and low back pain: a systematic, critical review of the epidemiological literature 1992–1999, International Archives of Occupational and Environmental Health 73 (2000) 290–297. [10] J.P.G. Urban, S. Roberts, Degeneration of the intervertebral disc, Arthritis Research & Therapy 5 (2003) 123–130. [11] S.D. Boden, D.O. Davis, T.S. Dina, N.J. Patronas, S.W. Wiesel, Abnormal magnetic-resonance scans of the lumbar spine in asymptomatic subjects. a prospective investigation, The Journal of Bone & Joint Surgery 72 (1990) 403–408. [12] S. Holm, A.K. Holm, L. Ekstrm, A. Karladani, T. Hansson, Experimental disc degeneration due to endplate injury, Journal of Spinal Disorders 17 (2004) 64–71. [13] L.I. Kerttula, W.S. Serlo, O.A. Tervonen, E.L. Pkk, H.V. Vanharanta, Post-traumatic findings of the spine after earlier vertebral fracture in young patients. Clinical and MRI study, Spine 25 (2000) 1104–1108. [14] J. Sandover, The fatigue approach to vibration and health: is it a practical and viable way of predicting the effects on people? Journal of Sound and Vibration 215 (1998) 699–721. [15] C.W.S. Chan, B. Gantenbein-Ritter, S.J. Ferguson, The effects of dynamic loading on the intervertebral disc, European Spine Journal 20 (2012) 1796–1812.

206

G. Marini et al. / Journal of Sound and Vibration 387 (2017) 194–206

[16] S. Kitazaki, M.J. Griffin, A modal analysis of whole-body vertical vibration, using a finite element model of the human body, Journal of Sound and Vibration 200 (1997) 83–103. [17] L. Wei, M.J. Griffin, Mathematical models for the apparent mass of the seated human body exposed to vertical vibration, Journal of Sound and Vibration 212 (1998) 855–874. [18] T. Belytschko, L. Schwer, A. Schultz, A model for analytic investigation of three-dimensional head-spine dynamics, Technical Report, DTIC Document, 1976. [19] C.V. Toen, M.M. Sran, S.N. Robinovitch, P.A. Cripton, Transmission of force in the lumbosacral spine during backward falls, Spine 37 (2012) E519. [20] Gan Z., Hillis A.J., Darling J., Biodynamic modelling of seated human subjects exposed to uncouples vertical and fore-and-aft whole-body vibration, Journal of Vibration Engineering and Technologies 3 (3) (2015) 301–314. [21] A.A.I. White, M.M. Panjabi, Clinical Biomechanics of the Spine, Lippincott Williams & Wilkins Company J.B., Philadelphia, USA, 1978. [22] G. Schmidt, A. Tondl, Non-Linear Vibration, Cambridge University Press, Cambridge, 1986. [23] G. Marini, G. Huber, K. Püschel, S.J. Ferguson, Nonlinear dynamics of the human lumbar intervertebral disc, Journal of Biomechanics 48 (2015) 479–488. [24] L.A. Miller, Structural dynamics and resonance in plants with nonlinear stiffness, Journal of Theoretical Biology 234 (2005) 511–524. [25] J.J. Costi, I.A. Stokes, M.G. Gardner-Morse, J.C. Iatridis, 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 (2008) 1731–1738. [26] R. Hamming, Numerical Methods for Scientists and Engineers, Dover Publications Inc., New York, 1986. [27] S. Li, A.G. Patwardhan, F.M.L. Amirouche, R. Havey, K.P. Meade, Limitations of the standard linear solid model of intervertebral discs subject to prolonged loading and low-frequency vibration in axial compression, Journal of Biomechanics 28 (1995) 779–790. [28] W. Koeller, S. Muehlhaus, W. Meier, F. Hartmann, Biomechanical properties of human intervertebral discs subjected to axial dynamic compression – influence of age and degeneration, Journal of Biomechanics 19 (1986) 807–816. [29] D. Jamison, M. Cannella, E.C. Pierce, M.S. Marcolongo, A comparison of the human lumbar intervertebral disc mechanical response to normal and impact loading conditions, Journal of Biomechanical Engineering 135 (2013) 091009. [30] J.C. Iatridis, L.A. Setton, M. Weidenbaum, V.C. Mow, The viscoelastic behavior of the non-degenerate human lumbar nucleus pulposus in shear, Journal of biomechanics 30 (1997) 1005–1013. [31] G. Marini, S.J. Ferguson, Nonlinear numerical analysis of the structural response of the intervertebral disc to impact loading, Computer Methods in Biomechanics and Biomedical Engineering 17 (2014) 1002–1011. [32] J. Rittweger, K. Just, K. Kautzsch, P. Reeg, D. Felsenberg, Treatment of chronic lower back pain with lumbar extension and whole-body vibration exercise. a randomized controlled trial, Spine 27 (2002) 1829–1834. [33] S. von Stengel, W. Kemmler, M. Bebenek, K. Engelke, W.A. Kalender, Effects of whole-body vibration training on different devices on bone mineral density, Medicine & Science in Sports & Exercise 43 (2011) 1071–1079. [34] J.-A. Boucher, J. Abboud, J.-D. Dubois, E. Legault, M. Descarreaux, Y. Henchoz, Trunk neuromuscular responses to a single whole-body vibration session in patients with chronic low back pain: a cross-sectional study, Journal of Manipulative and Physiological Therapeutics 36 (2013) 564–571. [35] G. Huber, D.M. Skrzypiec, A. Klein, K. Pschel, M. Morlock, High cycle fatigue behaviour of functional spinal units, Industrial Health 48 (2010) 550–556. [36] L.C. Andrews, Special Functions for Engineers and Applied Mathematicians, MacMillan, New York, 1998.