# Evaluation of the Use of the Yeoh and Mooney-Rivlin Functions as Strain Energy Density Functions for the Ground Substance Material of the Annulus Fibrosus

Authors: Héctor E. Jaramillo
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Due to the importance of the intervertebral disc in the mechanical behavior of the human spine, special attention has been paid to it during the development of finite element models of the human spine. The mechanical behavior of the intervertebral disc is nonlinear, heterogeneous, and anisotropic and, due to the low permeability, is usually represented as a hyperelastic model. The intervertebral disc is composed of the nucleus pulposus, the endplates, and the annulus fibrosus. The annulus fibrosus is modeled as a hyperelastic matrix reinforced with several fiber families, and researchers have used different strain energy density functions to represent it. This paper presents a comparative study between the strain energy density functions most frequently used to represent the mechanical behavior of the annulus fibrosus: the Yeoh and Mooney-Rivlin functions. A finite element model of the annulus fibrosus of the L4-L5 segment under the action of three independent and orthogonal moments of 8 N-m was used, employing Abaqus software. A structured mesh with eight divisions along the height and the radial direction of annulus fibrosus and tetrahedron elements for the endplates were used, and an exponential energy function was employed to represent the mechanical behavior of the fibers. A total of 16 families were used; the fiber orientation varied with the radial coordinate from 25° on the outer boundary to 46° on the inner boundary, measuring it with respect to the transverse plane. The mechanical constants were taken from the reported literature. The range of motion was obtained by finite element analysis using different values of the mechanical constants and these results were compared with the reported experimental data. It was found that the Yeoh function showed a better fit to the experimental range of motion than the Mooney-Rivlin function, especially in the nonlinear region.


---

## Body

## 1. Introduction

Generally, biomechanics research uses the finite element (FE) method [1–6]. The first FE model of the spine was reported by Belytschko and Kulak [7] and included two vertebral bodies and the intervertebral disc (IVD). This model was axisymmetric and the tissues were represented as isotropic and linear elastic material. The nucleus pulposus (NP) was modeled as incompressible liquid and the annulus fibrous (AF) as orthotropic linear material. However, the most important contributions to the modeling of the IVD were made by Shirazi-Adl et al. [8]. They reported a model of the L2-L3 segment that used geometry and nonlinear materials. The AF was represented as a matrix reinforced with collagen fibers, and this was possibly the first model of the AF in which it was modeled as a composite material. Also, the model included the differences in the three kinds of bones of the vertebral bodies.A few years later, Goel et al. [9] developed an FE model of the AF which takes into account the anisotropy of the fibers and their amount and orientation. However the fibers were considered as elastic material. More recently, Schmidt et al. [10] reported a new method of calibration of the properties of the AF (matrix and collagen fibers) using the experimental curves of the L4-L5 segment without NP. This method allows the individual contributions of the fibers and the matrix of the AF to be determined. On the other hand, a model developed by Ezquerro et al. [11] represented the NP as incompressible and hyperelastic using the Mooney-Rivlin energy function with two constants, and the AF was modeled as a fiber setup embedded into a matrix. They used a polynomial hyperelastic function for the matrix and stress-strain curves for the fibers. Guan et al. [12] reported an FE model in which the vertebrae were modeled as a rigid bodies, the NP as a hyperelastic Neo-Hookean material, the AF as a hyperelastic matrix fiber reinforced using a UMAT subroutine from Abaqus (www.3ds.com, Dassault Systèmes, Francia), and the ligaments as Truss elements under tension, while Gap elements were used to simulate the contact between the joint facets.In conclusion, many researchers have used different strain energy functions and techniques to model the mechanical behavior of the AF. Actually, the AF can be modeled as a reinforced matrix fiber. Particularly, the matrix or ground substance has been modeled using the energy functions of Mooney-Rivlin [13–18] and Yeoh [19–26]. However, we did not find any study that supported the use of these energy functions and how they affect the mechanical behavior of the AF. In this direction, this work compared the mechanical behavior of the AF determined using the two energy functions with the experimental data in order to find which strain energy function fit the experimental results better.


## 2. Methods and Materials

### 2.1. Geometry

The geometry of the AF of the L4-L5 segment (Figure1) was taken from the FE model of the L4-L5-S1 segment reported by Jaramillo et al. [27]. The percentage of the transverse section of the AF with respect to the total area of the IVD was 50% according to reported works [28–30].Figure 1
Shape and dimensions in mm for the annulus fibrosus of the L4-L5 segment.The AF model has a mesh with eight divisions along the height and eight in the radial direction (Figure2). The AF was divided into five sections in order to take into account the regional variation of the material properties and to obtain a better fit to the experimental data [10]. Each radial division contains a total of 16 families of fibers [31–33], where the families were crisscrossed with one another. The orientation of the fibers varied radially from 25° on the outer boundary to 46° on the inner boundary using an increase of 3°[10, 23, 34, 35], and the angle was measured with respect to the horizontal axis.Figure 2
Final mesh and sections of the annulus fibrosus used to define their different properties.The endplates were considered as a rigid body in order to apply the moments. A structured mesh of 7448 hexahedral elements with eight nodes was used for every model analyzed (Table1). Due to the differences between the types of elements used to represent the AF and the endplates, the Tie option was used to join the two parts (Figure 3).Table 1
Type and number of elements used in the finite element model.
Part Element type Element number Annulus fibrosus Hexahedral, 8 nodes, C3D8H 7448 End plates Tetrahedral, 10 nodes, C3D4 74134Figure 3
Section view of the finite element model.


### 2.2. Boundary Conditions and Loads

The inferior endplate was fixed and three independent moments of 8 N-m were applied to the superior endplate. The moments were applied through the orthogonal axes in order to simulate the flexion, extension, lateral flexion, and axial rotation movements.The range of motion (ROM) and the maximum stress in the families of fibers were selected as the output variables. The ROMs of the FE models were compared with the experimental data of the L4-L5 segment reported by Jaramillo et al. [36]. To compare the ROMs of the FE models with the experimental ROM, Equation (1) was employed, and the average ROM was calculated for 1, 2, 3, 4, 5, 6, 7, and 8 N-m in order to find how similar their behaviors were.(1)%Aprox.=ROMExp.-ROMFEAROMExp.x100%Also, the influence of the mechanical constants on the maximum stress for the fiber families of the AF was obtained.


### 2.3. Properties and Constitutive Equations

A strain energy function was defined for the fibers (Equation (2)) and their mechanical constants (Table 2) were used for all the models. (2)Wf=a1a2ea2I4-12+ea2I6-12-2,where a1 and a2 are the material constants andI4 andI6 are the deviatoric invariants associated with the two families of fibers, which are defined as(3)I4=N1.C.N1(4)I6=N2.C.N2where N(1) and N(2) are the unit vectors along the two fiber directions in the nondeformed configuration and C is the deviatoric right Green deformation tensor. The reinforcing fibers only work under a positive strain. This strain energy function was selected because it has been used by many researchers [22, 23, 37], and the experimental values for the mechanical constants have been reported [38–41].Table 2
Mechanical constants used to the fibers of the annulus fibrous.
Section a1 [MPa] a2 [Unit less] Anterior 2.10 25.0 Antero-lateral 0.16 4.0 Lateral 0.16 4.0 Postero-lateral 0.16 4.0 Posterior 0.24 5.64Mooney-Rivlin and Yeoh strain energy functions were used to represent the ground substance. The Mooney-Rivlin function was defined as(5)Wm=c10I1-3+c01I2-3wherec10 andc01 are the material constants andI1 andI2 are the first and second deviatoric invariants of the deviatoric right Green deformation tensor. The mechanical constants (Table 3) were taken from values reported by Heuer et al. [42–44] and Rohlman et al. [45]. For the FE analysis, the constants were multiplied by 0, 25, 50, 75, 100, 200, 300, 400, and 500%. In the first step,c10 was varied andc01 was constant, and in the second step,c10 was constant andc01 was varied. For instance, for model 3 (M3),c10 = 0.56 × 50% = 0.28 MPa andc01 = 0.14 × 100% = 0.14 MPa. In this way, in models M1 to M9,c10 was varied andc01 was constant; in models M10 to M18,c10 was constant andc01 was varied.Table 3
Mechanical constants for the strain energy function Mooney-Rivlin.
Section c 10 [MPa] c 01 [MPa] Anterior 0.56 0.14 Antero-lateral 0.56 0.14 Lateral 0.56 0.14 Postero-lateral 0.56 0.14 Posterior 0.56 0.14The Yeoh function was defined as(6)Wm=c1I1-3+c2I1-32+c3I1-33wherec1,c2, andc3 are the mechanical constants of the material andI1 is the first deviatoric invariant of the deviatoric right Green deformation tensor. The experimental constants were taken from Ayturk et al. [23] (Table 4) and multiplied by 0, 25, 50, 75, 100, 200, 300, 400, and 500%. In the first step,c1 was varied andc2 andc3 were constant; in the second step,c2 was varied andc1 andc3 were constant; finally, in the third step,c3 was varied andc1 andc2 were constant. For instance, in model 22 (M22), c1 = 0.0146 × 100% = 0.0146 MPa,c2 = –0.0189 × 100% = –0.0189 MPa, andc3 = 0.041 × 75% = 0.03075 MPa. In this way, in models M1 to M9, c1 was varied and c2 and c3 were constant; in models M10 to M18, c2 was varied and c1 and c3 were constant; and in models M19 to M27, c3 was varied and c1 and c2 were constant.Table 4
Mechanical constants used for the Yeoh strain energy function.
Section c 1 [MPa] c 2 [MPa] c 3 [MPa] Anterior 0.0146 -0.0189 0.041 Antero-lateral 0.0146 -0.0189 0.041 Lateral 0.0146 -0.0189 0.041 Postero-lateral 0.0146 -0.0189 0.041 Posterior 0.0146 -0.0189 0.041The Abaqus subroutine Uanisohyper was developed to implement the aforementioned energy functions. Also, due to the amount of data used, a subroutine in Python [46] was developed to execute every analysis automatically.


## 2.1. Geometry

The geometry of the AF of the L4-L5 segment (Figure1) was taken from the FE model of the L4-L5-S1 segment reported by Jaramillo et al. [27]. The percentage of the transverse section of the AF with respect to the total area of the IVD was 50% according to reported works [28–30].Figure 1
Shape and dimensions in mm for the annulus fibrosus of the L4-L5 segment.The AF model has a mesh with eight divisions along the height and eight in the radial direction (Figure2). The AF was divided into five sections in order to take into account the regional variation of the material properties and to obtain a better fit to the experimental data [10]. Each radial division contains a total of 16 families of fibers [31–33], where the families were crisscrossed with one another. The orientation of the fibers varied radially from 25° on the outer boundary to 46° on the inner boundary using an increase of 3°[10, 23, 34, 35], and the angle was measured with respect to the horizontal axis.Figure 2
Final mesh and sections of the annulus fibrosus used to define their different properties.The endplates were considered as a rigid body in order to apply the moments. A structured mesh of 7448 hexahedral elements with eight nodes was used for every model analyzed (Table1). Due to the differences between the types of elements used to represent the AF and the endplates, the Tie option was used to join the two parts (Figure 3).Table 1
Type and number of elements used in the finite element model.
Part Element type Element number Annulus fibrosus Hexahedral, 8 nodes, C3D8H 7448 End plates Tetrahedral, 10 nodes, C3D4 74134Figure 3
Section view of the finite element model.


## 2.2. Boundary Conditions and Loads

The inferior endplate was fixed and three independent moments of 8 N-m were applied to the superior endplate. The moments were applied through the orthogonal axes in order to simulate the flexion, extension, lateral flexion, and axial rotation movements.The range of motion (ROM) and the maximum stress in the families of fibers were selected as the output variables. The ROMs of the FE models were compared with the experimental data of the L4-L5 segment reported by Jaramillo et al. [36]. To compare the ROMs of the FE models with the experimental ROM, Equation (1) was employed, and the average ROM was calculated for 1, 2, 3, 4, 5, 6, 7, and 8 N-m in order to find how similar their behaviors were.(1)%Aprox.=ROMExp.-ROMFEAROMExp.x100%Also, the influence of the mechanical constants on the maximum stress for the fiber families of the AF was obtained.


## 2.3. Properties and Constitutive Equations

A strain energy function was defined for the fibers (Equation (2)) and their mechanical constants (Table 2) were used for all the models. (2)Wf=a1a2ea2I4-12+ea2I6-12-2,where a1 and a2 are the material constants andI4 andI6 are the deviatoric invariants associated with the two families of fibers, which are defined as(3)I4=N1.C.N1(4)I6=N2.C.N2where N(1) and N(2) are the unit vectors along the two fiber directions in the nondeformed configuration and C is the deviatoric right Green deformation tensor. The reinforcing fibers only work under a positive strain. This strain energy function was selected because it has been used by many researchers [22, 23, 37], and the experimental values for the mechanical constants have been reported [38–41].Table 2
Mechanical constants used to the fibers of the annulus fibrous.
Section a1 [MPa] a2 [Unit less] Anterior 2.10 25.0 Antero-lateral 0.16 4.0 Lateral 0.16 4.0 Postero-lateral 0.16 4.0 Posterior 0.24 5.64Mooney-Rivlin and Yeoh strain energy functions were used to represent the ground substance. The Mooney-Rivlin function was defined as(5)Wm=c10I1-3+c01I2-3wherec10 andc01 are the material constants andI1 andI2 are the first and second deviatoric invariants of the deviatoric right Green deformation tensor. The mechanical constants (Table 3) were taken from values reported by Heuer et al. [42–44] and Rohlman et al. [45]. For the FE analysis, the constants were multiplied by 0, 25, 50, 75, 100, 200, 300, 400, and 500%. In the first step,c10 was varied andc01 was constant, and in the second step,c10 was constant andc01 was varied. For instance, for model 3 (M3),c10 = 0.56 × 50% = 0.28 MPa andc01 = 0.14 × 100% = 0.14 MPa. In this way, in models M1 to M9,c10 was varied andc01 was constant; in models M10 to M18,c10 was constant andc01 was varied.Table 3
Mechanical constants for the strain energy function Mooney-Rivlin.
Section c 10 [MPa] c 01 [MPa] Anterior 0.56 0.14 Antero-lateral 0.56 0.14 Lateral 0.56 0.14 Postero-lateral 0.56 0.14 Posterior 0.56 0.14The Yeoh function was defined as(6)Wm=c1I1-3+c2I1-32+c3I1-33wherec1,c2, andc3 are the mechanical constants of the material andI1 is the first deviatoric invariant of the deviatoric right Green deformation tensor. The experimental constants were taken from Ayturk et al. [23] (Table 4) and multiplied by 0, 25, 50, 75, 100, 200, 300, 400, and 500%. In the first step,c1 was varied andc2 andc3 were constant; in the second step,c2 was varied andc1 andc3 were constant; finally, in the third step,c3 was varied andc1 andc2 were constant. For instance, in model 22 (M22), c1 = 0.0146 × 100% = 0.0146 MPa,c2 = –0.0189 × 100% = –0.0189 MPa, andc3 = 0.041 × 75% = 0.03075 MPa. In this way, in models M1 to M9, c1 was varied and c2 and c3 were constant; in models M10 to M18, c2 was varied and c1 and c3 were constant; and in models M19 to M27, c3 was varied and c1 and c2 were constant.Table 4
Mechanical constants used for the Yeoh strain energy function.
Section c 1 [MPa] c 2 [MPa] c 3 [MPa] Anterior 0.0146 -0.0189 0.041 Antero-lateral 0.0146 -0.0189 0.041 Lateral 0.0146 -0.0189 0.041 Postero-lateral 0.0146 -0.0189 0.041 Posterior 0.0146 -0.0189 0.041The Abaqus subroutine Uanisohyper was developed to implement the aforementioned energy functions. Also, due to the amount of data used, a subroutine in Python [46] was developed to execute every analysis automatically.


## 3. Results

### 3.1. Range of Motion (ROM) Behavior

Compared with the experimental ROM, the mechanical behavior of the Mooney-Rivlin function (Figures4–7) fitted less well than the behavior obtained with the Yeoh function (Figures 8–11). In this direction, the Mooney-Rivlin approximations to the experimental data were 26.1% for the flexion movement, 35.7% for the extension movement, 19.2% for the lateral flexion movement, and 39.4% for the axial rotation movement, while the Yeoh function approximations were 87.0, 95.1, 78.6, and 91.1% for the flexion, extension, lateral flexion, and axial rotation, respectively.Figure 4
ROM versus Moment for the flexion movement using Mooney-Rivlin function for different values ofc10 andc01.Figure 5
ROM versus Moment for the extension movement using Mooney-Rivlin function for different values ofc10 andc01.Figure 6
ROM versus Moment for the lateral flexion movement using Mooney-Rivlin function for different values ofc10 andc01.Figure 7
ROM versus Moment for the axial rotation movement using Mooney-Rivlin function for different values ofc10 andc01.Figure 8
ROM versus Moment for the flexion movement using Yeoh function for different values ofc1, c2, andc3.Figure 9
ROM versus Moment for the extension movement using Yeoh function for different values ofc1, c2, andc3.Figure 10
ROM versus Moment for the lateral flexion movement using Yeoh function for different values ofc1, c2 andc3.Figure 11
ROM versus Moment for the axial rotation movement using Yeoh function for different values ofc1, c2, andc3.


#### 3.1.1. Using the Mooney-Rivlin Function as the Strain Energy Function for the Ground Substance

The ROM behavior was similar for the four movements: asc10 andc01 increased, the ROM decreased, moving away from the experimental data. That is to say, when the mechanical constants increased, the ground substance contributed a greater part of the stiffness of the AF and as a consequence a lower ROM. Particularly, the set of mechanical constants with the better approximation to the experimental ROM wasc10 = 0 MPa andc01 = 0.14 MPa, with experimental approximations of 74.2, 82.4, 59.0, and 77.9% for the flexion, extension, lateral flexion, and axial rotation, respectively (Figures 4–7). Then, the left sides of Figures 4–7 show the ROM behavior whenc01 (0.14 MPa) is constant butc10 increases and the right sides show the ROM behavior whenc10 (0.56 MPa) is constant butc01 increases, for all movements. In general, thec10 variation has a major impact thanc01 over the ROM behavior.


#### 3.1.2. Using the Yeoh Function as the Strain Energy Function for the Ground Substance

When the Yeoh function was used, the ROM showed a better fit to the experimental behavior. In this direction, the fits obtained were 94.1% for flexion movement usingc1 = 0.0146 MPa,c2 = –0.0189 MPa, andc3 = 0.0205 MPa; 98.8% for extension movement usingc1 = 0.0292 MPa,c2 = –0.0189 MPa, andc3 = 0.041 MPa; 88.4% for lateral flexion movement usingc1 = 0.0146 MPa,c2 = –0.0189 MPa, andc3 = 0.01025 MPa; and 96.0% for axial rotation movement usingc1 = 0.0438 MPa,c2 = –0.0189 MPa, andc3 = 0.041 MPa (Figures 8–11). Just four curves were outside the experimental range for flexion (Figure 8), extension (Figure 9), and lateral flexion (Figure 10); these curves used the combination of c1 = 0.0146 MPa, c2 = –0.0189 MPa, and c3 = 0.082, 0.123, 0.164, and 0.205 MPa. For the axial rotation movement, all curves were within the experimental range (Figure 11). Then, Figures 8–11 are composed by three groups of figures, the first group (left side) show the ROM behavior whenc2 (-0.0189 MPa) andc3 (0.041) are constant butc1 increases, for the second group (middle)c1 (0.0146 MPa) andc3 (0.043) are constant butc2 increases, finally the third group (right side) show the ROM behavior whenc1 (0.0146 MPa) andc2 (-0.0189 MPa) are constant butc3 increases, for all movements.In general, the intervertebral disc stiffness increased due to the increase of the mechanical constants; however, a high impact ofc3 on the ROM was found in all movements. This can be explained by the high dispersion area found in the ROM behavior whenc3 varied (curves on the right side, Figures 8–11). In this case, the fits to the obtained experimental data were 84.6% for flexion, 93.4% for extension, 77.8% for lateral flexion, and 90.5% for axial rotation; these values of fit were under the average values obtained for all curves.


### 3.2. Maximum Stress in the Families of Fibers

In general, it was found that the ROM decreased due to increases of the mechanical constants of the ground substance (c1,c2, andc3). This can be explained by the fact that when the moment is applied, part of it is taken by the fibers and another by the matrix or ground substance, so when the mechanical constants of the ground substance increase, the matrix takes a greater part of the moment.


#### 3.2.1. Using the Mooney-Rivlin Function as the Strain Energy Function for the Ground Substance

Figures12–15 show the stress behavior versus the mechanical constants variation for a moment maximum of 8 N-m applied. Then, for a maximum moment of 8 N-m, the stress in the fibers varied from 0.17 to 47.1 MPa in flexion (Figure 12), from 0.57 to 125.5 MPa in extension (Figure 13), from 0.17 to 17.0 MPa in lateral flexion (Figure 14), and from 0.015 to 12.0 MPa in axial rotation (Figure 15). The percentage differences in the stresses between the fiber families were 16.5% for flexion, 12.8% for extension, 5.4% for lateral flexion, and 97.7% for axial rotation whenc10 was varied andc01 was equal to 0.14 MPa. Whenc10 was constant (0.56 MPa) andc01 was varied, the differences between the stresses for the two families of fibers were 9.4% for flexion, 16.8% for extension, 7.8% for lateral flexion, and 99.8% for axial rotation. In general, for all values ofc10 andc01, the stress differences between the two families of fibers were 12.9% for flexion, 14.8% for extension, 6.6% for lateral flexion, and 98.7% for axial rotation.Figure 12
Maximum stress in the fibers family of AF versus mechanical constants values for the flexion movement, using Mooney-Rivlin function.Figure 13
Maximum stress in the fibers family of AF versus mechanical constants values for the extension movement, using Mooney-Rivlin function.Figure 14
Maximum stress in the fibers family of AF versus mechanical constants values for the lateral flexion movement, using Mooney-Rivlin function.Figure 15
Maximum stress in the fibers family of AF versus mechanical constants values for the axial rotation movement, using Mooney-Rivlin function.


#### 3.2.2. Using the Yeoh Function as the Strain Energy Function for the Ground Substance

Figures16–19 show the stress behavior versus the mechanical constants variation for a moment maximum of 8 N-m applied. Then, for a maximum moment of 8 N-m, the fiber stress varied from 26.6 to 67.7 MPa for flexion (Figure 16), from 48.7 MPa to 250.0 MPa for extension (Figure 17), from 14.6 to 33.9 MPa for lateral flexion (Figure 18), and from 10.3 to 158.9 MPa for axial rotation (Figure 19). In the case whenc1 was varied andc2 (–0.0189 MPa) andc3 (0.041 MPa) were constant, the average differences in stresses between the two fiber families were 2.2% for flexion, 36.0% for extension, 14.1% for lateral flexion, and 69.9% for axial rotation. In the second case, whenc1 (0.0146 MPa) andc3 (0.041 MPa) were constant andc2 was varied, the stress differences between the two fiber families were 2.3% for flexion, 16.4% for extension, 15.1% for lateral flexion, and 68.7% for axial rotation. In the last case, whenc1 (0.0146 MPa) andc2 (–0.0189 MPa) were constant andc3 was varied, the stress differences between the two families of fibers were 3.3% for flexion, 26.2% for extension, 8.8% for lateral flexion, and 73.6% for axial rotation. In general, for all values ofc1, c2, andc3, the stress differences between the two families of fibers were 2.6% for flexion, 32.9% for extension, 12.7% for lateral flexion, and 70.5% for axial rotation.Figure 16
Maximum stress in the fibers family of AF versus mechanical constants values for the flexion movement, using Yeoh function.Figure 17
Maximum stress in the fibers family of AF versus mechanical constants values for the extension movement, using Yeoh function.Figure 18
Maximum stress in the fibers family of AF versus mechanical constants values for the lateral flexion movement, using Yeoh function.Figure 19
Maximum stress in the fibers family of AF versus mechanical constants values for the axial rotation movement, using Yeoh function.


## 3.1. Range of Motion (ROM) Behavior

Compared with the experimental ROM, the mechanical behavior of the Mooney-Rivlin function (Figures4–7) fitted less well than the behavior obtained with the Yeoh function (Figures 8–11). In this direction, the Mooney-Rivlin approximations to the experimental data were 26.1% for the flexion movement, 35.7% for the extension movement, 19.2% for the lateral flexion movement, and 39.4% for the axial rotation movement, while the Yeoh function approximations were 87.0, 95.1, 78.6, and 91.1% for the flexion, extension, lateral flexion, and axial rotation, respectively.Figure 4
ROM versus Moment for the flexion movement using Mooney-Rivlin function for different values ofc10 andc01.Figure 5
ROM versus Moment for the extension movement using Mooney-Rivlin function for different values ofc10 andc01.Figure 6
ROM versus Moment for the lateral flexion movement using Mooney-Rivlin function for different values ofc10 andc01.Figure 7
ROM versus Moment for the axial rotation movement using Mooney-Rivlin function for different values ofc10 andc01.Figure 8
ROM versus Moment for the flexion movement using Yeoh function for different values ofc1, c2, andc3.Figure 9
ROM versus Moment for the extension movement using Yeoh function for different values ofc1, c2, andc3.Figure 10
ROM versus Moment for the lateral flexion movement using Yeoh function for different values ofc1, c2 andc3.Figure 11
ROM versus Moment for the axial rotation movement using Yeoh function for different values ofc1, c2, andc3.


### 3.1.1. Using the Mooney-Rivlin Function as the Strain Energy Function for the Ground Substance

The ROM behavior was similar for the four movements: asc10 andc01 increased, the ROM decreased, moving away from the experimental data. That is to say, when the mechanical constants increased, the ground substance contributed a greater part of the stiffness of the AF and as a consequence a lower ROM. Particularly, the set of mechanical constants with the better approximation to the experimental ROM wasc10 = 0 MPa andc01 = 0.14 MPa, with experimental approximations of 74.2, 82.4, 59.0, and 77.9% for the flexion, extension, lateral flexion, and axial rotation, respectively (Figures 4–7). Then, the left sides of Figures 4–7 show the ROM behavior whenc01 (0.14 MPa) is constant butc10 increases and the right sides show the ROM behavior whenc10 (0.56 MPa) is constant butc01 increases, for all movements. In general, thec10 variation has a major impact thanc01 over the ROM behavior.


### 3.1.2. Using the Yeoh Function as the Strain Energy Function for the Ground Substance

When the Yeoh function was used, the ROM showed a better fit to the experimental behavior. In this direction, the fits obtained were 94.1% for flexion movement usingc1 = 0.0146 MPa,c2 = –0.0189 MPa, andc3 = 0.0205 MPa; 98.8% for extension movement usingc1 = 0.0292 MPa,c2 = –0.0189 MPa, andc3 = 0.041 MPa; 88.4% for lateral flexion movement usingc1 = 0.0146 MPa,c2 = –0.0189 MPa, andc3 = 0.01025 MPa; and 96.0% for axial rotation movement usingc1 = 0.0438 MPa,c2 = –0.0189 MPa, andc3 = 0.041 MPa (Figures 8–11). Just four curves were outside the experimental range for flexion (Figure 8), extension (Figure 9), and lateral flexion (Figure 10); these curves used the combination of c1 = 0.0146 MPa, c2 = –0.0189 MPa, and c3 = 0.082, 0.123, 0.164, and 0.205 MPa. For the axial rotation movement, all curves were within the experimental range (Figure 11). Then, Figures 8–11 are composed by three groups of figures, the first group (left side) show the ROM behavior whenc2 (-0.0189 MPa) andc3 (0.041) are constant butc1 increases, for the second group (middle)c1 (0.0146 MPa) andc3 (0.043) are constant butc2 increases, finally the third group (right side) show the ROM behavior whenc1 (0.0146 MPa) andc2 (-0.0189 MPa) are constant butc3 increases, for all movements.In general, the intervertebral disc stiffness increased due to the increase of the mechanical constants; however, a high impact ofc3 on the ROM was found in all movements. This can be explained by the high dispersion area found in the ROM behavior whenc3 varied (curves on the right side, Figures 8–11). In this case, the fits to the obtained experimental data were 84.6% for flexion, 93.4% for extension, 77.8% for lateral flexion, and 90.5% for axial rotation; these values of fit were under the average values obtained for all curves.


## 3.1.1. Using the Mooney-Rivlin Function as the Strain Energy Function for the Ground Substance

The ROM behavior was similar for the four movements: asc10 andc01 increased, the ROM decreased, moving away from the experimental data. That is to say, when the mechanical constants increased, the ground substance contributed a greater part of the stiffness of the AF and as a consequence a lower ROM. Particularly, the set of mechanical constants with the better approximation to the experimental ROM wasc10 = 0 MPa andc01 = 0.14 MPa, with experimental approximations of 74.2, 82.4, 59.0, and 77.9% for the flexion, extension, lateral flexion, and axial rotation, respectively (Figures 4–7). Then, the left sides of Figures 4–7 show the ROM behavior whenc01 (0.14 MPa) is constant butc10 increases and the right sides show the ROM behavior whenc10 (0.56 MPa) is constant butc01 increases, for all movements. In general, thec10 variation has a major impact thanc01 over the ROM behavior.


## 3.1.2. Using the Yeoh Function as the Strain Energy Function for the Ground Substance

When the Yeoh function was used, the ROM showed a better fit to the experimental behavior. In this direction, the fits obtained were 94.1% for flexion movement usingc1 = 0.0146 MPa,c2 = –0.0189 MPa, andc3 = 0.0205 MPa; 98.8% for extension movement usingc1 = 0.0292 MPa,c2 = –0.0189 MPa, andc3 = 0.041 MPa; 88.4% for lateral flexion movement usingc1 = 0.0146 MPa,c2 = –0.0189 MPa, andc3 = 0.01025 MPa; and 96.0% for axial rotation movement usingc1 = 0.0438 MPa,c2 = –0.0189 MPa, andc3 = 0.041 MPa (Figures 8–11). Just four curves were outside the experimental range for flexion (Figure 8), extension (Figure 9), and lateral flexion (Figure 10); these curves used the combination of c1 = 0.0146 MPa, c2 = –0.0189 MPa, and c3 = 0.082, 0.123, 0.164, and 0.205 MPa. For the axial rotation movement, all curves were within the experimental range (Figure 11). Then, Figures 8–11 are composed by three groups of figures, the first group (left side) show the ROM behavior whenc2 (-0.0189 MPa) andc3 (0.041) are constant butc1 increases, for the second group (middle)c1 (0.0146 MPa) andc3 (0.043) are constant butc2 increases, finally the third group (right side) show the ROM behavior whenc1 (0.0146 MPa) andc2 (-0.0189 MPa) are constant butc3 increases, for all movements.In general, the intervertebral disc stiffness increased due to the increase of the mechanical constants; however, a high impact ofc3 on the ROM was found in all movements. This can be explained by the high dispersion area found in the ROM behavior whenc3 varied (curves on the right side, Figures 8–11). In this case, the fits to the obtained experimental data were 84.6% for flexion, 93.4% for extension, 77.8% for lateral flexion, and 90.5% for axial rotation; these values of fit were under the average values obtained for all curves.


## 3.2. Maximum Stress in the Families of Fibers

In general, it was found that the ROM decreased due to increases of the mechanical constants of the ground substance (c1,c2, andc3). This can be explained by the fact that when the moment is applied, part of it is taken by the fibers and another by the matrix or ground substance, so when the mechanical constants of the ground substance increase, the matrix takes a greater part of the moment.


### 3.2.1. Using the Mooney-Rivlin Function as the Strain Energy Function for the Ground Substance

Figures12–15 show the stress behavior versus the mechanical constants variation for a moment maximum of 8 N-m applied. Then, for a maximum moment of 8 N-m, the stress in the fibers varied from 0.17 to 47.1 MPa in flexion (Figure 12), from 0.57 to 125.5 MPa in extension (Figure 13), from 0.17 to 17.0 MPa in lateral flexion (Figure 14), and from 0.015 to 12.0 MPa in axial rotation (Figure 15). The percentage differences in the stresses between the fiber families were 16.5% for flexion, 12.8% for extension, 5.4% for lateral flexion, and 97.7% for axial rotation whenc10 was varied andc01 was equal to 0.14 MPa. Whenc10 was constant (0.56 MPa) andc01 was varied, the differences between the stresses for the two families of fibers were 9.4% for flexion, 16.8% for extension, 7.8% for lateral flexion, and 99.8% for axial rotation. In general, for all values ofc10 andc01, the stress differences between the two families of fibers were 12.9% for flexion, 14.8% for extension, 6.6% for lateral flexion, and 98.7% for axial rotation.Figure 12
Maximum stress in the fibers family of AF versus mechanical constants values for the flexion movement, using Mooney-Rivlin function.Figure 13
Maximum stress in the fibers family of AF versus mechanical constants values for the extension movement, using Mooney-Rivlin function.Figure 14
Maximum stress in the fibers family of AF versus mechanical constants values for the lateral flexion movement, using Mooney-Rivlin function.Figure 15
Maximum stress in the fibers family of AF versus mechanical constants values for the axial rotation movement, using Mooney-Rivlin function.


### 3.2.2. Using the Yeoh Function as the Strain Energy Function for the Ground Substance

Figures16–19 show the stress behavior versus the mechanical constants variation for a moment maximum of 8 N-m applied. Then, for a maximum moment of 8 N-m, the fiber stress varied from 26.6 to 67.7 MPa for flexion (Figure 16), from 48.7 MPa to 250.0 MPa for extension (Figure 17), from 14.6 to 33.9 MPa for lateral flexion (Figure 18), and from 10.3 to 158.9 MPa for axial rotation (Figure 19). In the case whenc1 was varied andc2 (–0.0189 MPa) andc3 (0.041 MPa) were constant, the average differences in stresses between the two fiber families were 2.2% for flexion, 36.0% for extension, 14.1% for lateral flexion, and 69.9% for axial rotation. In the second case, whenc1 (0.0146 MPa) andc3 (0.041 MPa) were constant andc2 was varied, the stress differences between the two fiber families were 2.3% for flexion, 16.4% for extension, 15.1% for lateral flexion, and 68.7% for axial rotation. In the last case, whenc1 (0.0146 MPa) andc2 (–0.0189 MPa) were constant andc3 was varied, the stress differences between the two families of fibers were 3.3% for flexion, 26.2% for extension, 8.8% for lateral flexion, and 73.6% for axial rotation. In general, for all values ofc1, c2, andc3, the stress differences between the two families of fibers were 2.6% for flexion, 32.9% for extension, 12.7% for lateral flexion, and 70.5% for axial rotation.Figure 16
Maximum stress in the fibers family of AF versus mechanical constants values for the flexion movement, using Yeoh function.Figure 17
Maximum stress in the fibers family of AF versus mechanical constants values for the extension movement, using Yeoh function.Figure 18
Maximum stress in the fibers family of AF versus mechanical constants values for the lateral flexion movement, using Yeoh function.Figure 19
Maximum stress in the fibers family of AF versus mechanical constants values for the axial rotation movement, using Yeoh function.


## 3.2.1. Using the Mooney-Rivlin Function as the Strain Energy Function for the Ground Substance

Figures12–15 show the stress behavior versus the mechanical constants variation for a moment maximum of 8 N-m applied. Then, for a maximum moment of 8 N-m, the stress in the fibers varied from 0.17 to 47.1 MPa in flexion (Figure 12), from 0.57 to 125.5 MPa in extension (Figure 13), from 0.17 to 17.0 MPa in lateral flexion (Figure 14), and from 0.015 to 12.0 MPa in axial rotation (Figure 15). The percentage differences in the stresses between the fiber families were 16.5% for flexion, 12.8% for extension, 5.4% for lateral flexion, and 97.7% for axial rotation whenc10 was varied andc01 was equal to 0.14 MPa. Whenc10 was constant (0.56 MPa) andc01 was varied, the differences between the stresses for the two families of fibers were 9.4% for flexion, 16.8% for extension, 7.8% for lateral flexion, and 99.8% for axial rotation. In general, for all values ofc10 andc01, the stress differences between the two families of fibers were 12.9% for flexion, 14.8% for extension, 6.6% for lateral flexion, and 98.7% for axial rotation.Figure 12
Maximum stress in the fibers family of AF versus mechanical constants values for the flexion movement, using Mooney-Rivlin function.Figure 13
Maximum stress in the fibers family of AF versus mechanical constants values for the extension movement, using Mooney-Rivlin function.Figure 14
Maximum stress in the fibers family of AF versus mechanical constants values for the lateral flexion movement, using Mooney-Rivlin function.Figure 15
Maximum stress in the fibers family of AF versus mechanical constants values for the axial rotation movement, using Mooney-Rivlin function.


## 3.2.2. Using the Yeoh Function as the Strain Energy Function for the Ground Substance

Figures16–19 show the stress behavior versus the mechanical constants variation for a moment maximum of 8 N-m applied. Then, for a maximum moment of 8 N-m, the fiber stress varied from 26.6 to 67.7 MPa for flexion (Figure 16), from 48.7 MPa to 250.0 MPa for extension (Figure 17), from 14.6 to 33.9 MPa for lateral flexion (Figure 18), and from 10.3 to 158.9 MPa for axial rotation (Figure 19). In the case whenc1 was varied andc2 (–0.0189 MPa) andc3 (0.041 MPa) were constant, the average differences in stresses between the two fiber families were 2.2% for flexion, 36.0% for extension, 14.1% for lateral flexion, and 69.9% for axial rotation. In the second case, whenc1 (0.0146 MPa) andc3 (0.041 MPa) were constant andc2 was varied, the stress differences between the two fiber families were 2.3% for flexion, 16.4% for extension, 15.1% for lateral flexion, and 68.7% for axial rotation. In the last case, whenc1 (0.0146 MPa) andc2 (–0.0189 MPa) were constant andc3 was varied, the stress differences between the two families of fibers were 3.3% for flexion, 26.2% for extension, 8.8% for lateral flexion, and 73.6% for axial rotation. In general, for all values ofc1, c2, andc3, the stress differences between the two families of fibers were 2.6% for flexion, 32.9% for extension, 12.7% for lateral flexion, and 70.5% for axial rotation.Figure 16
Maximum stress in the fibers family of AF versus mechanical constants values for the flexion movement, using Yeoh function.Figure 17
Maximum stress in the fibers family of AF versus mechanical constants values for the extension movement, using Yeoh function.Figure 18
Maximum stress in the fibers family of AF versus mechanical constants values for the lateral flexion movement, using Yeoh function.Figure 19
Maximum stress in the fibers family of AF versus mechanical constants values for the axial rotation movement, using Yeoh function.


## 4. Discussion

In general, increases in the values of the mechanical constants produce an increase in the stiffness of the disc and decreases in the ROM and stress of the fibers. The decrease in stress in the fibers as a result of the increase of the mechanical constants of the ground substance can be explained by the fact that the matrix takes a greater part of the applied moment, decreasing the stress in the fiber.Using the Mooney-Rivlin function to represent the ground substance, the values of the mechanical constants that obtained the best fit to the experimental data were c10 = 0 and c01 = 0.14 MPa; here c01 is equal to the value reported by Heuer et al. [43, 47] and Rohlmann et al. [45]. However, these researchers reported values of c10 that differed from zero. On the other hand, the c01 value is very close to that reported (c01 = 0.25) by Campbell and Petrella [48]. With regard to the stress in the fiber families (flexion: 35.58 and 41.09 MPa; extension: 90.15 and 125.5 MPa; lateral flexion: 16.97 and 16.93 MPa; axial rotation: 11.95 and 82.5 MPa) for previous values of the mechanical constants, the values of the stress were within the range of the experimental ultimate strength reported by Iatridis et al. [49] of 88 ±61 MPa, but out of the range of the data reported by Fujita et al. [50] (0.37 ±0.20 MPa), Skaggs et al. [51] (10.3 ±8.4), and Green et al. [52] (8.6 ±4.3 MPa).Using the Yeoh function to represent the mechanical behavior of the ground substance produced several set of values of the mechanical constants with a good fit to the experimental data. Calculating the median and standard deviation for the set of values with the best fit to the experimental ROM in all movements, we obtain c1 = 0.0292 ± 0.0146 MPa, c2 = –0.0189 ±0 MPa, and c3 = 0.025625 ± 0.015375 MPa. The last values of c1, c2, and c3 are within the experimental range reported by O’Connell et al. [53] but are out of the experimental range reported by Cortes et al. [38]. The differences found may be due to the regional variation of the mechanical properties of the matrix [39] and the significant differences between the values reported in the literature [50–52, 54, 55] for these mechanical properties which are produced by the different protocols used to obtain them. Now, if we wish to obtain a good fit to the experimental ROM and the fiber stress below the ultimate strength, it is necessary to use the obtained values of superior limits of the mechanical constants, so c1 = 0.0438 MPa, c2 = –0.0189 MPa, and c3 = 0.041 MPa. With these constants, the fiber stresses are 47.44 and 48.56 MPa, which are within the experimental range reported by Iatridis et al. [49].


## 5. Conclusion

It is necessary to underline that the Yeoh function produced, as a result, an important set of values that fit the experimental ROM better than the Mooney-Rivlin function. Besides, the Yeoh function has a better fit in the nonlinear region of the experimental ROM than the Mooney-Rivlin function, due to the fact that the Yeoh function has the invariants of second and third grade, whereas the Mooney-Rivlin function has invariants of first grade.



---

*Source: 1570142-2018-11-07.xml*