Open access peer-reviewed chapter

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

Written By

Sabine Bauer and Ivanna Kramer

Submitted: 27 November 2020 Reviewed: 30 April 2021 Published: 08 June 2021

DOI: 10.5772/intechopen.98211

From the Edited Volume

Recent Advances in Numerical Simulations

Edited by Francisco Bulnes and Jan Peter Hessling

Chapter metrics overview

336 Chapter Downloads

View Full Metrics


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.


  • multi-body simulation
  • sensitivity analysis
  • cervical spine FSU model
  • intervertebral disc pressure
  • stiffness and damping coefficients

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/s2 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 mg/cm3 for C6 and 414 mg/cm3 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.

VertebraMass [kg]CGx [m]CGy [m]CGz [m]
C70.0070−8.23 × 10−9−9.13 × 10−9−8.37 × 10−9
C60.00576.2 × 10−44.1 × 10−32.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.

VertebraIxx [kg m2]Iyy [kg m2]Izz [kg m2]Ixy [kg m2]Ixz [kg m2]Iyz [kg m2]
C71.41 × 10−61.48 × 10−62.50 × 10−6−3.87 × 10−86.89 × 10−8−2.44 × 10−8
C66.77 × 10−71.24 × 10−61.65 × 10−6 1−1.47 × 10−86.25 × 10−81.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 Fx, Fy, Fz 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 cx, cy, cz and dx, dy, dz respectively. The disc deformation value is calculated as a distance between two points and is represented by the variables xF, yF, and zF, where the axis-wise velocities of the markers are ẋ, ẏ and ż. The disc force is defined in Eq. (1).


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 m2 with the average facet area superior C7 and inferior C6 reported in [17], a discrepancy of FCA = 0.000089 m2 (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).

ModelC7 EPWu [mm]C7 EPDu [mm]C6 EPWl [mm]C6 EPDl [mm]C7 EPAu [mm2]C6 EPAl [mm2]C6-C7 DW [mm]C6-C7 DD [mm]C6-C7 DH [mm]C6-C7 DA [mm2]
Current FSU model21.716.919.215.4288.0232.220.4516.150.0068260.1
Hueston et al. [23]
Tan et al. [25]
Yoganandan et al. [20]0.005–0.0075168–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 δPi=PiPinit, where the initial pressure value is 0.301315MPa.

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 50000Ns/m. For the sake of simplicity, the same value ci is assigned to cx,cy,cz 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 f0.1,0.2,0.3,0.4,0.5.

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×105MPa 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 di is assigned to dx,dy,dz in each i-th simulation run. The stiffness term is set to be constant, i.e. ci=500000N/m for all N trials. The damping parameter di is increased and decreased by a fraction f0.1,0.2,0.3,0.4,0.5 of its initial value dinit=50000Ns/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 di are reflected in the linear changes of the disc pressure pi. 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.00593MPa for d=25000Ns/m, which is approximately 1.991% of the initial pressure value. In comparison, for d=75000Ns/m the change is 0.00226MPa 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.301315MPa, 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 sv is defined to be an average quotient of the disc pressure change pi to the i-th change in the parameter value vi:

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%.


where N is a total number of trials and δvi is the i-th change in the observing parameter, vicidi.

Determined coefficients for the stiffness sc and damping term sd 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.

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 lL, where L=100N200N800N 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: c2.5×106,3.0×106 as well as for c4.0×106,4.5×106 and 6.5×106,7.0×106.

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=50000Nm/s.

Figure 13 illustrates the results of the simulations where the applied loads lL 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×104Ns/m for the applied force of 200N, 3.5×104Ns/m for 100N load, another peaks are observed for the exerting force of 500N at damping value of 5.0×104Ns/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 mm2 and maximum of 502 mm2 are published in [20]. Estimated disc areas of 180 mm2, 230 mm2 and 295 mm2 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 mm2 and 269 mm2 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 R2 is calculated. R2 is defined to be the square of Pearson correlation coefficient rx,y, i.e. R2=rx,y2. Pearson correlation coefficient is determined for data pairs x1y1(xnyn) 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].


where N is sample size, xi,yi are the individual sample points, x¯=1Ni=1Nxi the sample mean for x, which is calculated for y¯ analogically and iN.

The calculated R2 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:


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×106 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.



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.



MBSMulti-Body Simulation
FSUFunctional Spinal Unit
CoGCenter of Gravity
CRCenter of Rotation
MIMoment of Inertia
ALLAnterior Longitudinal Ligament
PLLPosterior Longitudinal Ligament
FLFlava Ligament
ISLInterspinous Ligament
NLNuchal Ligament
CLCapsular Ligament
SPLSpinous Process Length
lleft or lower (depending on the context)


  1. 1. Bauer S. Basics of Multibody Systems: Presented by Practical Simulation Examples of Spine Models. In: Lopez-Ruiz R editor. Numerical Simulation. InTech; 2016. pp. 29-49. DOI: 10.5772/62864
  2. 2. Galbusera F, Schmidt H, Neidlinger-Wilke C, Wilke HJ: The effect of degenerative morphological changes of the intervertebral disc on the lumbar spine biomechanics: a poroelastic finite element investigation. Comput Methods Biomech Biomed Engin. 2011;14(8):729-39. DOI: 10.1080/10255842.2010.493522
  3. 3. Ryang Y, Pape H, Meyer B: Degenerative Lumbale Instabilität - Definition, klinische und radiologische Zeichen, Management. Die Wirbelsäule. 2017;01(02):101-116. isbn:2509-8241
  4. 4. Bashkuev M, Reitmaier S, Schmidt H: Effect of disc degeneration on the mechanical behavior of the human lumbar spine: a probabilistic finite element study.The Spine Journal. 2018;18(10):1910-1920., isbn:1529-9430
  5. 5. Mas Y, Gracia L, Ibarz E, Gabarre S, Pena D, Herrera A: Finite element simulation and clinical follow-up of lumbar spine biomechanics with dynamic fixations. PloS one. 2017;12(11):1-19.
  6. 6. Schmidt H, Heuer F, Wilke H-J: Which axial and bending stiffnesses of posterior implants are required to design a flexible lumbar stabilization system?. Journal of Biomechanics. 2009;42(1):48-54.isbn:0021-9290
  7. 7. Wilke H-J, Heuer F, Schmidt H: Design optimization of a new posterior dynamic stabilization system. Journal of Biomechanics. 2008;41. DOI 10.1016/S0021- 9290(08)70312-9
  8. 8. Goel VK, Grauer JN, Patel TCh, Biyani A, Sairyo K, Vishnubhotla S, Matyas A, Cowgill I, Shaw M, Long R, Dick D, Panjabi MM, Serhan H: Effects of charit´e artificial disc on the implanted and adjacent spinal segments mechanics using a hybrid testing protocol. Spine (Phila Pa 1976). 2005;30(24):2755-64. DOI: 10.1097/01.brs.0000195897.17277.67.
  9. 9. Dreischarf M, Zander T, ShiraziAdl A, Puttlitz C.M, Adam C.J, Chen C.S, Goel V.K, Kiapour A, Kim, Y.H, Labus K.M, Little J.P, Park W.M, Wang Y.H, Wilke H.J, Rohlmann A, Schmidt H: Comparison of eight published static finite element models of the intact lumbar spine: predictive power of models improves when combined together.Journal of biomechanics. 2014;47:1757-1766
  10. 10. Zander T, Dreischarf M, Timm AK, Baumann WW, Schmidt H: Impact of material and morphological parameters on the mechanical response of the lumbar spine - A finite element sensitivity study. J Biomech. 2017;53:185-190. DOI: 10.1016/j.jbiomech.2016.12.014
  11. 11. Panjabi MM, Duranceau J, Goel V, Oxland T, Takata K: Cervical human vertebrae. Quantitative three-dimensional anatomy of the middle and lower regions. Spine. 1991;16.8:861-9. DOI: 10.1097/00007632-199108000-00001
  12. 12. Anderst WJ, Thorhauer ED, Lee JY, Donaldson WF, Kang JD: Cervical spine bone mineral density as a function of vertebral level and anatomic location. Spine J. 2011;11(7):659-67. DOI: 10.1016/j.spinee.2011.05.007
  13. 13. White AA, Panjabi MM: Clinical biomechanics of the spine. 2nd ed. Philadelphia: Lippincott; 1990;752 p.
  14. 14. van Lopik DW, Acar M. Development of a multi-body computational model of human head and neck. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics. 2007;221(2):175-197. DOI: 10.1243/14644193JMBD84
  15. 15. Panjabi MM, Crisco JJ, Vasavada A, Oda T, Cholewicki J, Nibu K, Shin E: Mechanical properties of the human cervical spine as shown by three-dimensional load-displacement curves. Spine (Phila Pa 1976). 2001;26(24):2692-700. DOI: 10.1097/00007632-200112150-00012
  16. 16. Yoganandan N, Knowles SA, Maiman DJ, Pintar FA: Anatomic study of the morphology of human cervical facet joint. Spine (Phila Pa 1976). 2003;28(20):2317-23. DOI: 10.1097/01.BRS.0000085356.89103.A5
  17. 17. Panjabi MM, Oxland T, Takata K, Goel V, Duranceau J, Krag M: Articular facets of the human spine. Quantitative three-dimensional anatomy. Spine (Phila Pa 1976). 1993;18(10):1298-310. DOI: 10.1097/00007632-199308000-00009
  18. 18. Yang KH, King AI: Mechanism of facet load transmission as a hypothesis for low-back pain. Spine (Phila Pa 1976). 1984;9(6):557-65. DOI: 10.1097/00007632-198409000-00005
  19. 19. Jager, de MKJ: Mathematical head-neck models for acceleration impacts. Eindhoven: Eindhoven University of Technology; 1996. 143 p.
  20. 20. Yoganandan N, Kumaresan S, Pintar FA: Biomechanics of the cervical spine Part 2. Cervical spine soft tissue responses and biomechanical modeling. Clin Biomech (Bristol, Avon). 2001;16(1):1-27. DOI: 10.1016/s0268-0033(00)00074-7
  21. 21. Kadri, PA, Al-Mefty O: Anatomy of the nuchal ligament and its surgical applications. 2007;61:301–304. DOI: 10.1227/01.neu.0000303985.65117.ea
  22. 22. Mattucci SF, Moulton JA, Chandrashekar N, Cronin DS: Strain rate dependent properties of younger human cervical spine ligaments. J Mech Behav Biomed Mater. 2012;10:216-26. DOI: 10.1016/j.jmbbm.2012.02.004
  23. 23. Hueston S, Makola M, Mabe I, Goswami T: Cervical Spine Anthropometric and Finite Element Biomechanical Analysis. 2012. IntechOpen., doi:10.5772/35524.
  24. 24. Kumaresan S, Yoganandan N, Pintar FA, Maiman DJ, Goel VK: Contribution of disc degeneration to osteophyte formation in the cervical spine: a biomechanical investigation. J Orthop Res. 2001 Sep;19(5):977-84. DOI: 10.1016/S0736-0266(01)00010-9
  25. 25. Tan SH, Teo EC, Chua HC: Quantitative three-dimensional anatomy of cervical, thoracic and lumbar vertebrae of Chinese Singaporeans. Eur Spine J. 2004;13(2):137-46. DOI: 10.1007/s00586-003-0586-z
  26. 26. Pooni J, Hukins D, Harris P, Hilton R, Davies K: Comparison of the structure of human intervertebral discs in the cervical, thoracic and lumbar regions of the spine. Surgical and Radiologic Anatomy. 2006;8:175-182. DOI: 10.5772/35524
  27. 27. Skrzypiec DM, Pollintine P, Przybyla A, Dolan P, Adams MA: The internal mechanical properties of cervical intervertebral discs as revealed by stress profilometry. Eur Spine J. 2007 Oct;16(10):1701-9. DOI: 10.1007/s00586-007-0458-z
  28. 28. Sellers WI, Crompton RH: Using sensitivity analysis to validate the predictions of a biomechanical model of bite forces. Ann Anat. 2004;186(1):89-95. DOI: 10.1016/S0940-9602(04)80132-8.
  29. 29. Wexler P, Anderson B, Gad S, Hakkinen B, Kamrin M, De Peyster A, Locey B, Pope C, Mehendale H, Shugart L: Encyclopedia of toxicology. 3rd ed. Waltham; Academic Press; 2005. 236-237p. DOI: 10.1177/1091581815586498
  30. 30. Smith E, Szidarovszky F, Karnavas W, Bahill A. Sensitivity analysis, a powerful system validation technique. The Open Cybernetics Systemics Journal. 2008; 2(1). DOI: 10.2174/1874110X00802010039
  31. 31. Taylor M: What is sensitivity analysis. Consortium YHE. University of York. 2009: 1-8.
  32. 32. Qian G, Mahdi A: Sensitivity analysis methods in the biomedical sciences. Mathematical biosciences. 2020;323:108306. DOI: 10.1016/j.mbs.2020.108306
  33. 33. Lehman S, Lawrence S: Three algorithms for interpreting models consisting of ordinary differential equations: sensitivity coefficients, sensitivity functions, global optimization. Mathematical Biosciences. 1982;62.1:107-122. DOI: 10.1016/0025-5564(82)90064-5
  34. 34. Ghosh P. The biology of the intervertebral disc. 1st ed. Boca Raton, FL: CRC press; 1988. DOI: 10.1007/978-3-7091-1535-0

Written By

Sabine Bauer and Ivanna Kramer

Submitted: 27 November 2020 Reviewed: 30 April 2021 Published: 08 June 2021