Parameter Dependencies of a Biomechanical Cervical Spine FSU - The Process of Finding Optimal Model Parameters by Sensitivity Analysis

Abstract

The knowledge about the impact of structure-specific parameters on the biomechanical behavior of a computer model has an essential meaning for the realistic modeling and system improving. Especially the biomechanical parameters of the intervertebral discs, the ligamentous structures and the facet joints are seen in the literature as significant components of a spine model, which define the quality of the model. Therefore, it is important to understand how the variations of input parameters for these components affect the entire model and its individual structures. Sensitivity analysis can be used to gain the required knowledge about the correlation of the input and output variables in a complex spinal model. The present study analyses the influence of the biomechanical parameters of the intervertebral disc using different sensitivity analysis methods to optimize the spine model parameters. The analysis is performed with a multi-body simulation model of the cervical functional spinal unit C6-C7.

1. Introduction

Biomechanical modeling offers a non-invasive possibility to analyze and answer kinematic and kinetic questions. A distinction is made between finite element (FE) simulation and multi-body simulation (MBS). The difference between FE and MBS modeling lies in the basic model structure and thus in the field of application. Further information on FE and MBS are described in [1]. Due to their complexity, FE models make an important contribution to understand the biomechanical function of the spine and the behavior of spinal structures in the state of health, illness or damage [2, 3, 4] as well as the influence of the material parameters of various implants and fusion techniques [5, 6, 7, 8]. However, the complexity of the FE models usually requires high computing times for each simulation case. If the aspect of predicting kinematic and dynamic reactions of the whole or a larger part of the spinal column during complex movement sequences is the focus of interest, MBS is a suitable simulation method due to the highly efficient short computing times [1]. The existing FE models are mostly based only on a specific or an idealized average model with unique mechanical and geometrical characteristics. According to [9], a better insight into the influence of the biomaterial and the geometrical diversity on the biomechanical behavior of the spine is essential for a better understanding of the spine mechanics and the patient care. Because a model contains numerous parameters that are often only vaguely known and too complex to implement, their effect on the responses is a priori unknown and full validation is largely impossible. Therefore there is a need for sensitivity analysis [10]. Sensitivity analysis can be used to gain the required knowledge about the correlation of the input and output variables in a complex spinal model, which has an essential meaning for the realistic modeling and system optimization. Especially the biomechanical parameters of the intervertebral discs, the ligamentous and muscular structures and the facet joints are significant components of a spine model, which define the quality of the model. Hence, it is important to understand, how the variations of input parameters for these components affect the entire model and its individual structures. The present study analyses the influence of the biomechanical parameters of the intervertebral disc using different sensitivity analysis methods, which enables the direct optimization of the spine model parameters. The analysis is performed with a multi-body simulation model of the cervical functional spinal unit C6-C7.

2. Model configuration

The MBS model of a functional spinal unit (FSU) consists of the vertebrae C7 and C6 represented by rigid bodies. Furthermore, an intervertebral disc, ligamentous structures and facet joints are implemented with specific biomechanical characteristics. When configuring the model, the focus is on the creation of the simplest possible model so that all biomechanical parameters could be adequately defined. In the case of models with many parameters, there is a risk that the parameters cannot be determined sufficiently. Therefore, a model with a high parameter dependency does not necessarily lead to better results.

2.1 FSU setup

The 3D surface of the vertebrae based on artificial vertebrae (Sawbones) and were implemented as triangular meshes. The 3D geometry data of the C7 vertebra can be taken out of Figure 1 and of C6 out of Figure 2. The abbreviations for different vertebral parts are made up of three capitalized letters and adapted from [11]. The first two describe the corresponding vertebral part and the third represents the dimension to be measured. These three letter combinations can be supplemented by a lower case letter that indicates a direction, such as right (r), upper (u) and depending on the content, lower or left (l).

Figure 1.

Geometry data of the vertebra C7.

Figure 2.

Geometry data of the vertebra C6.

The body reference frame of two vertebrae C7 and C6 are located at the same position. The location of the center of gravity (CG) of both vertebrae is visualized in Figure 3 and the data can be taken out of Table 1. The coordinates of the CG are given relative to the reference frame of the corresponding vertebra. The center of gravity is defined as a point at which the entire mass of the vertebra is united and where the earth gravitational acceleration with $g=-9.81m/s^2$ along the z-axis is applied. The mass properties of the rigid vertebrae are automatically calculate from their 3D geometry. For each single vertebra, the tessellated volume of its 3D data is multiplied with the specific density. The density is specified by [12] with 473 $\mathit{mg}/{\mathit{cm}}^3$ for C6 and 414 $\mathit{mg}/{\mathit{cm}}^3$ for C7.

Figure 3.

Center of gravity (CG), reference frame of the vertebra C7 and C6, location of center of rotation (CR) and load application point.

Vertebra	Mass [kg]	CGx [m]	CGy [m]	CGz [m]
C7	0.0070	−8.23 × 10^−9	−9.13 × 10^−9	−8.37 × 10^−9
C6	0.0057	6.2 × 10^−4	4.1 × 10^−3	2.0 × 10^−2

Table 1.

Mass and Center of Gravity of vertebrae C7 and C6.

Position of gravity center CG is given with respect to local body coordinates system.

The information of the inertia moment (Table 2) relates to the body’s center of gravity. The moment of inertia (I) is defined as a symmetric matrix whose entries are mirror-symmetric with respect to the main diagonal and relative to center of mass of the corresponding vertebra.

Vertebra	Ixx [kg m^2]	Iyy [kg m^2]	Izz [kg m^2]	Ixy [kg m^2]	Ixz [kg m^2]	Iyz [kg m^2]
C7	1.41 × 10^−6	1.48 × 10^−6	2.50 × 10^−6	−3.87 × 10^−8	6.89 × 10^−8	−2.44 × 10^−8
C6	6.77 × 10^−7	1.24 × 10^−6	1.65 × 10^−6 1	−1.47 × 10^−8	6.25 × 10^−8	1.36 × 10^−9

Table 2.

Moments of inertia with respect to local body center of gravity for the vertebrae of FSU C7-C6.

Moments of inertia I are given with respect to the local body center of gravity.

2.2 Intervertebral disc modeling

The biomechanical characteristic of the intervertebral disc between C7 and C6 is represented by a simple stiffness-deformation relation and a velocity-dependent damping term. If a load is applied to the model, the disc is deformed and develops reaction forces that depend on the deformation value and the deformation velocity. The forces $F_x$, $F_y$, $F_z$ are interacting between two defined markers, one refers to C7 and another to C6 in three translation directions x, y, z. The corresponding force equation is determined by four main components: stiffness constant $c$, damping constant $d$, disc deformation and deformation velocity. The stiffness term $c$ as well as the damping constant $d$ are represented separately for each translation direction $c_x$, $c_y$, $c_z$ and $d_x$, $d_y$, $d_z$ respectively. The disc deformation value is calculated as a distance between two points and is represented by the variables $x_F$, $y_F$, and $z_F$, where the axis-wise velocities of the markers are $\dot{x}$, $\dot{y}$ and $\dot{z}$. The disc force is defined in Eq. (1).

$$\left(\begin{array}{c}{F}_x\\ F_y\\ F_z\end{array}\right)=\left(\begin{array}{l}{c}_x{x}_F+d_x\dot{x}\\ c_y{y}_F+d_y\dot{y}\\ c_z{z}_F+d_z\dot{z}\end{array}\right)$$E1

The disc force is implemented in such a way, that its responds depend on specific movement scenarios: if the disc is deformed by an external load and the deformation velocity vector is negative, then the disc force is determined by both the stiffness and the damping terms. If the intervertebral disc is still deformed but begins to relax, then the deformation velocity vector changes into a positive direction. In this state, the disc force is only determined by the stiffness term. If the intervertebral disc is stretched, both terms are set to zero (Figure 4).

Figure 4.

Schematic representation of the intervertebral disc characteristics under different stress scenarios.

The initial value for the stiffness bases on [13] and the damping value is set to 10% of the stiffness value, because no actual cervical spine disc damping coefficients have been reported in the literature [14].

In reality the intervertebral disc is not only deformed by loads, but also bent by external torques. Depending on the action direction of the external torque the intervertebral disc performs a flexion and extension movement, an axial rotation or a lateral flexion. To counteract this rotations, the intervertebral disc develops a counter-torque. This non-linear disc torque is defined by two-dimensional functions that describe the relationship between the disc torque and the relative angle. A specific input function is assigned to the torques acting around three axes of rotation x, y and z. The applied input function bases on [15].

2.3 Facet joint modeling

Through the facets, adjacent vertebrae are connected via a thin layer of cartilage. In the model the facet cartilage layers are approximated by an unilateral contact spring-damper element, whose contact area is determined by the facet geometry. The contact area is a rectangular region, which represents the facet width and height. With an additional dimension the cartilage layer of the facet joint is simulated. The cartilage layer thickness of the lower cervical spine bases on [16] and is determined to be 0.00045 m for the superior layer and 0.00049 m for the inferior. The parameterization of the geometry, positioning and orientation of the 3D facet contact area is determined with respect to the C7 upper facet surface. The modeled facet contact surface is assumed to be an average facet width and facet height of the superior facet surface C7 and the inferior facet surface C6. The average model geometry results in a facet width (FW) of 0.0094 m and a facet height (FH) of 0.009 m. Comparison of the approximate facet area (FCA) of the current model with FCA = 0.000085 m² with the average facet area superior C7 and inferior C6 reported in [17], a discrepancy of FCA = 0.000089 m² (Figure 5) can be observed. This information is given at this point in order to show the extent to which the model assumption differs from the experimental measurements with regard to the geometry.

Figure 5.

Representation of the facet width and height, which builds the basis area for the facet contact simulation.

The stiffness coefficients are taken from [18] and the damping values is defined as 10% of the stiffness term. The damping coefficient is used to obtain a better attenuation of the maximum linear and angular accelerations of the head [19].

2.4 Ligament modeling

The spinal ligaments provide stability to the motion segments allowing motion within physiological limits. Ligaments are uniaxial structures that resist only tensile or distractive forces becoming slack in compression [14, 20].

In the FSU model the following ligaments are incorporated: anterior and posterior longitudinal ligament (ALL and PLL), flava ligament (FL), interspinous ligament (ISL), nuchal ligament (NL) and the left and right capsular ligaments (CL) (Figure 6). Ligaments, which have a broad structure, are represented by several fiber bundles. For instant, ALL and PLL are composed of a right, left and middle ligament structure. CL is approximated by four individual ligament structures that attach to the top, bottom, left and right surfaces of the articular processes. The ISL extends over the entire edge of the spinous process and is therefore modeled using three bundles of ligaments. The LF attaches to the proximal edge of the lamina and is represented by six ligament fibers. The NL is an extension of the SSL which extends from the external occipital protuberance to the spinous process of C7 and attaches all the posterior tips of the spinous processes in between [21].

Figure 6.

Representation of the ligament attachment points.

The determination of the ligament attachment points is carried out on the basis of the vertebral geometry and is checked by an expert.

The ligament’s characteristic is modeled by the load displacement curves [13, 22, 23]. When a ligament is stretched, it develops a force that is specific to the ligament in question. It acts against the direction of the stretch with no resistance in compression.

2.5 Load case configuration

In order to analyze the reaction of the spinal structures to a load, an external force of 80 N is applied to the endplate of the vertebra C6. This loading case is chosen because the cervical spine is permanently loaded by the weight of the head [24]. To prevent additional torques, the y-coordinates of the external load markers have the same position as the y-coordinates of the disc joint, so that there is no initial lever arm that could lead to unintentional torques.

2.6 Model validation

An important step in the simulation process is the model validation, with which the simulation results are checked for correctness. The correctness of the FSU is proven by comparing the intervertebral disc pressure and disc deformation to existing published data. After researching the literature, it turned out that there is only a limited possibility of validation data that exactly depicts the simulation scenario we have modeled at the moment. In general, there is the difficulty that the own model configuration does not necessarily exactly match to that of other researchers, since different specific research questions have to be answered. In order to get the response of the FSU model to different loads, the FSU is exposed to small and large external loads. The disc pressure and deformation are compared (Table 3).

Model	C7 EPWu [mm]	C7 EPDu [mm]	C6 EPWl [mm]	C6 EPDl [mm]	C7 EPAu [mm^2]	C6 EPAl [mm^2]	C6-C7 DW [mm]	C6-C7 DD [mm]	C6-C7 DH [mm]	C6-C7 DA [mm^2]
Current FSU model	21.7	16.9	19.2	15.4	288.0	232.2	20.45	16.15	0.0068	260.1
Hueston et al. [23]	19.0	15.1	19.5	15.7	220.8	316.3				268.6
Tan et al. [25]	19.0	15.1	19.5	15.7	220.8	316.3
Yoganandan et al. [20]									0.005–0.0075	168–502
Pooni et al. [26]										200–502

Table 3.

Comparison of the vertebra C7 and C6 anthropometry.

The width (W), depth (D) and area (A) of upper (u) and lower (l) endplates (EP) are presented of different studies. Further, the disc width (DW), the disc depth (DD), the disc area (DA) and the disc height (DH) is presented. The idea is to present the various possible measures to be able to assess the model parameters of the current model.

2.7 Motion segment response to small loads

A validated intact FE model of the C4-C5-C6 cervical spine to simulate progressive disc degeneration at the C5-C6 level is presented by [24]. The intact and three degenerated cervical spine models are exercised under the compression load of 80 N. The results of the intact spine model are used to compare the intervertebral disc pressure between vertebrae C5-C6 in the current FSU model. The motion segments were subjected to a small static compression load of 80 N in z-direction. While in the current model the resulting displacement of the intervertebral disc is measured, in the FSU model the overall force displacement response of C4 with respect to C6 is determined. Therefore, the comparison can only be taken as a rough evaluation of the models deformation.

2.8 Motion segment response to large loads

In the second stage of validation, the FSU model is subjected to larger loads of 200 N, 500 N and 673 N to determine its intervertebral disc pressure and disc deformation. The load of 200 N is chosen to represent the combined effects of head weight and muscle tension [27]. The human cervical disc pressure using a pressure transducer, side-mounted in a 0.9 mm diameter needle is investigated by [27]. Forty-six cadaverous cervical motion segments aged 48–90 years are subjected to a compressing load of 200 N for 2 s. Due to the lack of data available for high load cases, these data are used to analyze the characteristics of the intervertebral discs. The deformation value under a certain load is only provided for the specific healthy disc segment C7-T1. These results are used to compare the characteristics of the intervertebral discs in the current model.

A MBS model of human head and neck C7-T1 is presented by [14]. The MBS model comprise soft tissues, i.e. muscles, ligaments, intervertebral discs and supported through facet joints. Also eighteen muscle groups and 69 individual muscle segments of the head and neck are included in the model. For load–displacement testing, each motion segment is mounted so that the inferior vertebra is rigidly fixed whereas the superior vertebra is free to move in response to the applied loads. The response of model motion segments C5–C6 to the applied translational load of 500 N is shown.

A review with the focus on soft tissue structural responses with an emphasis on finite element mathematical models is done by [20]. Biomechanical data of intervertebral disc under compression test are provided for the FSU C6-C7. Under a load of 673 N the intervertebral disc between vertebrae C6 and C7 is deformed with 1.7 mm.

The comparison of the intervertebral disc pressure under the loads of 80 N and 200 N are shown in Figure 7.

Figure 7.

Response of model motion segments to applied compressive loads of 80 N and 200 N. the orange bar shows the intervertebral disc pressure for the corresponding motion segment as reported by [24] and the yellow one as reported by [27]. The results of the current FSU model is highlighted in blue.

Figure 8 compares the effect on the intervertebral disc deformation of the loads 80 N, 500 N and 673 N. The intradiscal pressure almost agrees with the published pressure at the loads of 80 N and 673 N. The current model has about one-third lower disc displacement than the comparison model [14]. One reason for this can be, that under certain circumstances the muscles, that are not taken into account in the current FSU model, are accompanied by a lack of muscle tension, which leads to the less compression of the intervertebral discs.

Figure 8.

Comparison of the intervertebral disc displacement with experimental results [20], FE model [24] and MBS model [14]. The results of the currant FSU model are presented in blue bars, die orange bar represent the results of [24], the yellow one of [14] and the green one of [20].

3. Effects of stiffness and damping variations

3.1 Method

In the biomechanical modeling the quality of the simulation can be considered valid only when the model and the input parameters are accurate and robust [28]. To examine the robustness of the modeled system a method called sensitivity analysis is the first choice. Generally speaking, sensitivity analysis is collection of approaches, that determine, quantify and analyze the impact of the input parameters on the model outcome [29]. The sensitivity analysis can also identify those components of the model that might need additional studies to be performed. Further, in the model optimization the sensitivity analysis can be used to refine the values of the critical parameters as well as to simplify or ignore those factors, which do not show any impact on the model response [30].

One of the simplest and effective techniques used to determine the level of the sensitivity or insensitivity of the model outputs to the plausible variation of one particular parameter is one-way sensitivity analysis [31, 32].

In order to identify the effect of both stiffness and damping variations on the current FSU model, one-way sensitivity analysis is performed. After every simulation run the corresponding changes in the intervertebral disc pressure in the C6-C7 segment are reported as a difference between initial and current disc pressure $\delta P_i=P_i-P_{\mathit{\text{init}}}$, where the initial pressure value is $0.301315\mathit{MPa}$.

3.2 Simulation results of stiffness term alternation

The first series of simulations consider the variation of the stiffness term $c$, which is expressed in Eq. (1), however the damping value $d$ is held constant by $50000\mathit{Ns}/m$. For the sake of simplicity, the same value $c_i$ is assigned to $c_x,c_y,c_z$ in every $i$-th simulation. Starting from the initial value of $500000N/m$ in each experiment repetition the stiffness parameter is increased and decreased by a fraction $f\in \left\{\mathrm{0.1,0.2,0.3,0.4,0.5}\right\}$.

It can be seen in Figure 9, that the variation in the stiffness term results in the linear change of the disc pressure, which points out to the linear relationship between these two variables. However, the change plot indicates the opposite linear relationship between the stiffness and the disc pressure, where increasing of the stiffness causes decreasing of the pressure. Note, that the course of the disc pressure changes is symmetrical, i.e. the minimum and maximum changes in the pressure value are of the same magnitude. According to the revealed results the maximal absolute pressure change is reported to be $7.0935249\times {10}^{-5}\mathit{MPa}$ or 0.02% of the initial value.

Figure 9.

Decreasing (left side) and increasing (right side) the stiffness term $c$ by factor $f$ impact the intradiscal pressure. The pressure change at the initial point is $0$ and is marked green.

3.3 Simulation results of damping term alternation

The second part of the experiments aims the analysis of the system sensitivity with respect to alternations of the damping constant $d$ (see Eq. (1)). In order to simplify the experiment execution, the same value $d_i$ is assigned to $d_x,d_y,d_z$ in each $i$-th simulation run. The stiffness term is set to be constant, i.e. $c_i=500000N/m$ for all $N$ trials. The damping parameter $d_i$ is increased and decreased by a fraction $f\in \left\{\mathrm{0.1,0.2,0.3,0.4,0.5}\right\}$ of its initial value $d_{\mathit{\text{init}}}=50000\mathit{Ns}/m$.

The impact of the damping parameter changes on the disc pressure is presented in Figure 10. Similarly to the results obtained in Section 3.2, an obvious influence of the damping term on the disc pressure can be observed, where the linear changes of $d_i$ are reflected in the linear changes of the disc pressure $p_i$. However, the magnitude of the change is not symmetrical for decreasing and increasing values of the damping factor. Moreover, the smaller damping values result in higher changes of the intradiscal pressure. The maximal pressure change is found to be $0.00593\mathit{MPa}$ for $d=25000\mathit{Ns}/m$, which is approximately 1.991% of the initial pressure value. In comparison, for $d=75000\mathit{Ns}/m$ the change is $-0.00226\mathit{MPa}$ and 0.6% respectively.

Figure 10.

Decreasing (left side) and increasing (right side) the damping value by $f$ effects the disc pressure in a linear manner. The disc pressure at the initial point is $0.301315\mathit{MPa}$, the change of pressure at initial point is 0 (marked green).

In Figure 11 the disc pressure changes affected by percentage decreasing of both model factors $d$ and $c$ following the one-way sensitivity analysis approach are depicted. It can be seen, that the same alternation of the damping term causes approximately two orders higher magnitude of the disc pressure. To examine the hypothesis, that the damping parameter has much stronger influence on the system, the calculation of a further sensitivity metric called sensitivity coefficient is elaborated [33]. In our particular setting the sensitivity coefficient $s_v$ is defined to be an average quotient of the disc pressure change $p_i$ to the $i$-th change in the parameter value $v_i$:

Figure 11.

Comparison of the maximal disc pressure changes [MPa] given with respect to the variability of the input parameters $c$ and $d$. The stiffness (marked red) and damping (marked blue) terms are decreased by the factor of 10–50%.

$$s_v=\frac{1}{N}\sum \frac{\delta p_i}{\delta v_i},$$E2

where $N$ is a total number of trials and $\delta v_i$ is the $i$-th change in the observing parameter, $v_i\in \left\{c_i,d_i\right\}$.

Determined coefficients for the stiffness $s_c$ and damping term $s_d$ are shown in Table 4. The obtained results support the above statement, that the behavior of the current model is approximately 12.59 times more sensitive to the damping term than to the stiffness parameter.

	sc	sd
Sensitivity coefficient value	−1.145 × 10^−5	−9.093^−510^−7

Table 4.

Sensitivity coefficient determined for parameters $c$ and $d$ using Eq. (2).

4. Impact of different load cases and intervertebral disc areas on intervertebral disc pressure

4.1 Impact of different loads on intervertebral disc pressure

One of the main functional task of the intervertebral disc is transmitting the compressive loads through the spine [34]. Therefore, it is important to study the sensitivity of the input parameters as well as mechanical responses of the model considering multiple loading cases. For this experiment, the acting of various external loads $l\in L$, where $L=\left\{100N,200N,\dots ,800N\right\}$ on the upper endplate of the C6 vertebra (see Figure 3) is simulated. Such high forces are selected in order to investigate the model behavior under different boundary conditions.

The disc pressure responded by the current model is reported in Figure 12. It can be seen, that the stiffness alternations among the load cases do not lead to significant change in the disc pressure. An unusual pattern is observed in each particular load situation, where the stiffness variation causes the linear growth in the disc pressure followed by piece-wise non-linear regions. Please note, that this disc behavior is not detected in the simulations applying the load of 80 N (see Sections 3.2 and 3.3). In details, the non-linear pressure change is detected withing the following ranges of the stiffness term: $c\in \left[2.5\times {10}^6\mathrm{,3.0}\times {10}^6\right]$ as well as for $c\in \left[4.0\times {10}^6\mathrm{,4.5}\times {10}^6\right]$ and $\left[6.5\times {10}^6\mathrm{,7.0}\times {10}^6\right]$.

Figure 12.

Maximal intradiscal pressure for C6-C7 segment calculated for multiple compressive loads and different stiffness value. The initial stiffness term is decreased and increased by a factor up to 50% of its initial value $c=500000N/s$. The damping term is help constant $d=50000\mathit{Nm}/s$.

Figure 13 illustrates the results of the simulations where the applied loads $l\in L$ and the damping factor $d$ are varied simultaneously. The perspective view of this diagram is slightly different in order to emphasize the regions where the non-linearity in the disc pressure occurs. The dark spots in the plot indicate the jumps in the disc pressure value over the loads, where the step-wise patterns show the non-linear responses of the current model for the following cases: $3.0\times {10}^4\mathit{Ns}/m$ for the applied force of $200N$, $3.5\times {10}^4\mathit{Ns}/m$ for $100N$ load, another peaks are observed for the exerting force of $500N$ at damping value of $5.0\times {10}^4\mathit{Ns}/m$.

Figure 13.

Maximal intradiscal pressure calculated for C6-C7 segment. Different compressive loads and values of the damping term are simultaneously changed, however the stiffness factor is set constant to $c=500000N/s$. The simulation results reveal the areas (dark spot regions in the plot) with a non-linear changes of the maximal disc pressure.

4.2 Impact of different intervertebral disc areas on disc pressure

The size of the disc area is presented in the literature with different values. This leads to the question how different disc areas influence the disc pressure. In order to investigate this effect, approximated intervertebral disc areas from the literature [11, 20, 25, 26] are used as examples. Values of the FSU C6-C7 disc area with minimum of 168 mm² and maximum of 502 ${\mathit{mm}}^2$ are published in [20]. Estimated disc areas of 180 mm², 230 mm² and 295 mm² are published by [26] and represented as mean values from 3 specimens of their cervical spine. The EPAu of C7 and the EPAl of C6 is specified in [11, 25]. The mean values of EPAu of C7 and EPAl of C6 with 284 mm² and 269 mm² respectively are taken to approximate the area of the corresponding intervertebral disc.

All disc area values listed above serve as input for the analysis of the relationship between intervertebral disc pressure and disc area. The effects of different intervertebral disc area on the intradiscal pressure are considered under the load case of 80 N.

In Figure 14 it can be clearly seen that the size of the disc area has a direct effect on the disc pressure. The course can be approximated by a 3-degree polynomial. In order to assess the goodness of the polynomial fit, the coefficient of determination $R^2$ is calculated. $R^2$ is defined to be the square of Pearson correlation coefficient $r_{x,y}$, i.e. $R^2=r_{x,y}^2$. Pearson correlation coefficient is determined for data pairs $\left\{\left(x_1,y_1\right),\dots ,(x_n,y_n)\right\}$ as follows:

Figure 14.

Representation of the relationship between disc area size and intradiscal pressure. The pressure is determined for eight disc areas of different sizes and an external load of 80 N. the blue points in the plot are the data points connected by a best-fit straight line. The disc areas are based on literature data. The assignment of the data points with their specific intervertebral disc areas to the corresponding literature is as follows (starting from the top left side): Data point 1 [20], data point 2 [26], data point 3 [26], data point 4 (current FSU model), data point 5 [25], data point 6 [11], data point 7 [26], data point 8 [20].

$$r_{x,y}=\frac{{\sum }_{i=1}^N\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right)}{\sqrt{{\sum }_{i=1}^N{\left(x_i-\overline{x}\right)}^2{\sum }_{i=1}^N{\left(y_i-\overline{y}\right)}^2}},$$E3

where $N$ is sample size, $x_i,y_i$ are the individual sample points, $\overline{x}=\frac{1}{N}{\sum }_{i=1}^N{x}_i$ the sample mean for $x$, which is calculated for $\overline{y}$ analogically and $i\in N$.

The calculated R² value of 1 (see Figure 14) shows a perfect correlation of the variables. Based on this correlation, the disc pressure can be determined by means of the third degree polynomial for given disc surface areas. The resulted polynomial is:

$$r_{x,y}=-1.16\times {10}^{10}{x}^3+1.39\times {10}^7{x}^2+5.89\times {10}^3+1.05.$$E4

This method offers the possibility of comparing and checking the disc pressure calculated in the simulation model with the one determined by the polynomial.

5. Conclusions

This study should be seen as a first approach to analyze the cervical spine’s sensibility to different influencing factors. The focus is on the analysis of the effects of various stiffness and damping parameters and disc area on the intradiscal pressure of the FSU C6-C7 in order to indicate the model weaknesses and optimize the model design.

In the first part of this study an one-way sensitivity analysis is performed in order to indicate, whether one of the given input parameter, namely stiffness or damping term, has a dominant influence on the model behavior. The experimental results show, that both parameters exhibit an identical impact on the disc pressure. However the variations of the damping term indicate a slightly stronger effect on the intradiscal pressure measurements, which is reflected in relatively higher value of the calculated sensitivity coefficient. When applying compressive loads from 100 N up to 800 N on the FSU model and varying the analyzing parameters a not foreseeable response pattern in the disc pressure is explored. Simultaneous change of the load and the corresponding parameter values results in a non-linear outcome regarding the intradiscal pressure, which is not detected in the simulations that consider the exerting external force of 80 N.

Further, it could be shown that the correlation between disc area and disc pressure can be approximated by a third-degree polynomial. This allows a further possibility for model validation of the simulated intervertebral disc pressure. For this purpose, the simulation result can be compared with the intervertebral disc pressure calculated by the polynomial with a known disc surface area.

An essential point to be considered in the next step is the implementation of the musculature. This is not taken into account in this model. It is still unclear what influence other cervical parameters, e.g. the facet joints, ligaments or muscles have and how these affect the overall mechanic when changed. Therefore, following this investigation, the effect of model parameters of others spinal structures, such as facet alignment and size, on the load on the intervertebral discs will be evaluated. Further, it must be questioned critically whether these results can be transferred to a model with a larger spinal column section. In order to discuss this question, in a further step not only an FSU should be considered, but the sensitivity of model parameters in a model that contains an entire spinal column section should be analyzed.

In case when additional elements are integrated into the model and the number of input factors grows, another broadly used method called multivariate sensitivity analysis can be applied in order to investigate the model response affected by the simultaneous variations of the underlying parameters. This procedure can help to optimize the model structure by finding the variables, that primarily impact the model outcomes. Moreover, using the sensitivity analysis methods the values of the principal parameters can be determined so that realistic simulation of model behavior is possible.

The experimental design of the presented sensitivity analysis follows the recommendations found in the literature. In the future work, the boundary conditions of the experiments should be extended. For instance, the range of the stiffness value might be increased up to $8.3\times {10}^6$ as it was used in the model proposed in [14]. Then the response of the current FSU model can be compared with the outputs of the referenced model.

Acknowledgments

We like to thank Prof. Dietrich Paulus, Institute of Computer Visualistics, University Koblenz-Landau for the fruitful discussion and Dr. Francis Kilian, Head of the Clinic for Spinal Surgery, Head of the Spinal Center Catholic Clinic Koblenz-Montabaur and PD Dr. Roland Jacob, specialist in ear, nose and throat medicine for guidance on medical and anatomical questions.

Conflict of interest

The authors declare no conflict of interest.

Abbreviations

MBS	Multi-Body Simulation
FSU	Functional Spinal Unit
CoG	Center of Gravity
CR	Center of Rotation
MI	Moment of Inertia
ALL	Anterior Longitudinal Ligament
PLL	Posterior Longitudinal Ligament
FL	Flava Ligament
ISL	Interspinous Ligament
NL	Nuchal Ligament
CL	Capsular Ligament
SPL	Spinous Process Length
FC	Facet
DC	Disc
EP	Endplate
A	Area
W	Width
H	Height
D	Depth
l	left or lower (depending on the context)
r	right
u	upper