**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–4] as well as the influence of the material parameters of various implants and fusion techniques [5–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).

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

*Parameter Dependencies of a Biomechanical Cervical Spine FSU - The Process of Finding…*

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

mass of the corresponding vertebra.

*Geometry data of the vertebra C6.*

**Figure 1.**

**Figure 2.**

**79**

*Geometry data of the vertebra C7.*

*DOI: http://dx.doi.org/10.5772/intechopen.98211*

**2.2 Intervertebral disc modeling**

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*:*81 *m=s* <sup>2</sup> 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=cm*<sup>3</sup> for C6 and 414 *mg=cm*<sup>3</sup> for C7.

*Parameter Dependencies of a Biomechanical Cervical Spine FSU - The Process of Finding… DOI: http://dx.doi.org/10.5772/intechopen.98211*

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.

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

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),

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

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

<sup>2</sup> along the z-axis is

**2. Model configuration**

*Recent Advances in Numerical Simulations*

**2.1 FSU setup**

**78**

does not necessarily lead to better results.

*mg=cm*<sup>3</sup> for C6 and 414 *mg=cm*<sup>3</sup> for C7.

upper (u) and depending on the content, lower or left (l).

where the earth gravitational acceleration with *g* ¼ �9*:*81 *m=s*

**Figure 2.** *Geometry data of the vertebra C6.*

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.

### **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 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

*Parameter Dependencies of a Biomechanical Cervical Spine FSU - The Process of Finding…*

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 coeffi-

intervertebral disc is stretched, both terms are set to zero (**Figure 4**).

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

cients have been reported in the literature [14].

*DOI: http://dx.doi.org/10.5772/intechopen.98211*

**Figure 4.**

**81**

#### **Figure 3.**

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


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

#### **Table 1.**

*Mass and Center of Gravity of vertebrae C7 and C6.*


#### **Table 2.**

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

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 axiswise velocities of the markers are *x*\_, *y*\_ and *z*\_. The disc force is defined in Eq. (1).

$$
\begin{pmatrix} F\_x \\ F\_y \\ F\_x \end{pmatrix} = \begin{pmatrix} c\_x \mathbf{x}\_F + d\_x \dot{\mathbf{x}} \\ c\_\mathcal{I} \mathbf{y}\_F + d\_\mathcal{I} \dot{\mathbf{y}} \\ c\_x \mathbf{z}\_F + d\_x \dot{\mathbf{z}} \end{pmatrix} \tag{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 *Parameter Dependencies of a Biomechanical Cervical Spine FSU - The Process of Finding… DOI: http://dx.doi.org/10.5772/intechopen.98211*

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**).

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

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 axiswise velocities of the markers are *x*\_, *y*\_ and *z*\_. The disc force is defined in Eq. (1).

C7 1.41 � <sup>10</sup>�<sup>6</sup> 1.48 � <sup>10</sup>�<sup>6</sup> 2.50 � <sup>10</sup>�<sup>6</sup> �3.87 � <sup>10</sup>�<sup>8</sup> 6.89 � <sup>10</sup>�<sup>8</sup> �2.44 � <sup>10</sup>�<sup>8</sup> C6 6.77 � <sup>10</sup>�<sup>7</sup> 1.24 � <sup>10</sup>�<sup>6</sup> 1.65 � <sup>10</sup>�<sup>6</sup> <sup>1</sup> �1.47 � <sup>10</sup>�<sup>8</sup> 6.25 � <sup>10</sup>�<sup>8</sup> 1.36 � <sup>10</sup>�<sup>9</sup>

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

**Vertebra Mass [kg] CG***<sup>x</sup>* **[m] CG***<sup>y</sup>* **[m] CG***<sup>z</sup>* **[m]** C7 0.0070 �8.23 � <sup>10</sup>�<sup>9</sup> �9.13 � <sup>10</sup>�<sup>9</sup> �8.37 � <sup>10</sup>�<sup>9</sup> C6 0.0057 6.2 � <sup>10</sup>�<sup>4</sup> 4.1 � <sup>10</sup>�<sup>3</sup> 2.0 � <sup>10</sup>�<sup>2</sup>

**] I***zz* **[kg m<sup>2</sup>**

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

0

B@

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

*cxxF* þ *dxx*\_ *cyyF* þ *dyy*\_ *czzF* þ *dzz*\_

**] I***xy* **[kg m<sup>2</sup>**

**] I***xz* **[kg m<sup>2</sup>**

**] I***yz* **[kg m2**

**]**

1

CA (1)

*Fx Fy Fz* 1

CA <sup>¼</sup>

0

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

**] I***yy* **[kg m2**

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

*Mass and Center of Gravity of vertebrae C7 and C6.*

*Recent Advances in Numerical Simulations*

**Vertebra I***xx* **[kg m2**

**Figure 3.**

**Table 1.**

**Table 2.**

**80**

*application point.*

B@

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.4 Ligament modeling**

*DOI: http://dx.doi.org/10.5772/intechopen.98211*

**Figure 6.**

**83**

*Representation of the ligament attachment points.*

The spinal ligaments provide stability to the motion segments allowing motion within physiological limits. Ligaments are uniaxial structures that resist only tensile

*Parameter Dependencies of a Biomechanical Cervical Spine FSU - The Process of Finding…*

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

or distractive forces becoming slack in compression [14, 20].

#### **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.

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

*Parameter Dependencies of a Biomechanical Cervical Spine FSU - The Process of Finding… DOI: http://dx.doi.org/10.5772/intechopen.98211*
