### Mathematical model

Practically, there is generally too much noise in the raw load-deflection data to numerically determine compliance, which is the first derivative of the load-deflection curve. Also the fact that data points are not necessarily equally spaced in time makes data processing harder. To overcome these problems, we propose to mathematically fit the experimental load-deflection data. The mathematical fit filters the noise and allows for an analytical calculation of the non-linear segment compliance. To create a mathematical model, attention has to be paid to several points. First, we look for the region of minimal stiffness. This implies that the borders of the NZ must be objectively defined in terms of stiffness or its inverse, compliance. Thus, the NZ should not be dependent on the arbitrarily chosen maximum load. Secondly, the asymmetry of the load-deflection curve should be taken into account, because of the asymmetrical structure of the intervertebral disc. Finally, the NZ should be independent from the starting point of the load-deflection curve or the position of the spinal segment after embedding. As a full load-deflection curve shows hysteresis due to poro-visco-elastic behaviour (Figure 1), the fit should be made once over the curve along one loading direction (e.g. full flexion to full extension) and once over the curve along the other loading direction (e.g. full extension to full flexion).. The analysis thus results in two values for both, the neutral zone magnitude and its stiffness, which are hypothesised to be equal.

There exists no theory that prescribes the type of function, which describes the mechanical behaviour of the segment (*e.g*. under flexion-extension). However, it is important to note that the NZ exists and actually only can be defined by virtue of the characteristic S-shape of the load-deflection curve. The double sigmoid function has such a shape and was found to provide an excellent fit over the entire curve (Figure 2):

D=\frac{1}{1+{e}^{-(a1+b1L)}}*c1+\frac{1}{1+{e}^{-(a2+b2L)}}*c2+d

(1)

Deflection D of the motion segment is expressed as a function of external load L. The sigmoid function is chosen because of its S-like shape and it is summed to allow for the natural asymmetry in the segment's flexion-extension curve. The parameters a and d in the equation determine the location of the load-deflection curve by shifting it horizontally (along the x- or loading axis) and vertically (along the y- or deflection axis), respectively. Parameters b and c reflect the biomechanical properties of the spinal motion segment. Parameter b determines the slope of the curve and thus the NZ stiffness. Parameter c determines the range of motion (ROM) of the segment but by changing the slope of the curve also the NZ stiffness. Eq.1 was fitted to experimental data, using unconstrained nonlinear minimization of the root mean squared error in a custom made matlab (the Mathworks Natick MA, USA) routine, employing the "fminsearch" function from the matlab optimization toolbox.

The NZ is that part of the load-deflection curve where compliance is maximal. The summed sigmoid function (Eq.1) can be used to determine compliance analytically as the first derivative (dD/dL; Figure 2). To obtain the range of maximal compliance (or minimal stiffness, *i.e*. the NZ), the *change* in compliance should be determined, which can be done with the second derivative (d^{2}D/dL^{2}). This curve provides unique inflection points that indicate the beginning and the end of the NZ, defined by the maximum and minimum in the second derivative (Figure 2).

### Neutral Zone Stiffness

The neutral zone stiffness itself is defined by the inverse of the compliance, which in turn is defined by the slope of the load-deflection curve over the full range of the NZ. To quantify this, a straight line was fitted over the load-deflection curve in the range of the NZ.

### Experimental data

To test the new mathematical definition of the NZ, we used experimental load-deflection curves from other (unpublished) studies. These load-deflection curves were determined from eight porcine lumbar spinal motion segments, which were loaded in flexion and extension to a maximum of 2 Nm in each direction. To this end, segments were embedded in woods metal and mounted on a four-point bending device controlled by a hydraulic material testing machine (Instron 8872, High Wycombe, UK)[12, 13] (Figure 3). Bending moments were applied by a constant rate displacement of the cross-head of the materials testing machine, providing a bending rate of 1.5°/s. During the experiments, the segments were kept wet by frequently spraying saline. The position of the vertebrae was determined with a Sonometrics Digital Ultrasonic Measurement System (Sonometrics Corporation, London, Canada), consisting of piezo-electric crystals sending and receiving pulses. Three crystals were mounted on a jib attached to each vertebra and submerged in a water bath (Figure 3). Three more crystals were connected to earth for reference. The distance between the crystals was calculated from the time difference between sending and receiving the pulses. Custom-made software was used to calculate vertebral motions from the distance between the separate crystals. In line with the recommendations by Wilke et al[7]., the third loading cycle in each experiments was used to quantify the NZ magnitude and stiffness.

### Compressive loading

To determine the effect of prolonged axial loading, 250N of axial compression was applied on the same eight segments for a period of seven hours (Figure 3). During this period, the specimens were wrapped in a saline soaked cloth to prevent dissication. After removal of the axial load the experiments were repeated. NZ magnitude and stiffness were calculated by both, the double sigmoid function and the angulation difference. Hystersis of the segment was calculated from the area between the flexion-extension and extension-flexion curves.

### Statistics

Goodness of fit of the double sigmoid function to experimental load-deflection data of procine SMS tested before and after 7 hours of axial loading was quantified by means of the coefficient of determination (r^{2}, the square of Pearson's coefficient of correlation) between predicted and measured defection data. Linearity of the load-deflection data in the NZ was evaluated by the coefficient of determination (r^{2}) of the linear least-squares fit.

To test to what extent the same results were obtained for NZ magnitude and NZ stiffness from flexion-extension and extension-flexion curves, the data were compared by means of paired t-tests and in addition the intra-class correlation (ICC) between these observations was determined. An ICC < 0.60 was qualified as poor, 0.60≤ ICC < 0.80 as moderate and ICC≥0.80 as good.

Next, the results were averaged for both movement directions and the NZ magnitudes and NZ stiffness obtained with the proposed method were compared to those obtained with the method of Wilke[7] using paired t-tests. In addition, the extent to which the two methods yielded common information was expressed by means of Pearson's coefficient of correlation.

The effect of axial loading on NZ magnitude and NZ stiffness as determined with both methods was determined using paired t-tests and a similar comparison was made between the values for hysteresis before and after loading.