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

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 elementmodels of the human spine.Themechanical 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-mwas used, employing Abaqus software. A structuredmesh 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.Themechanical 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.


Introduction
Generally, biomechanics research uses the finite element (FE) method [1][2][3][4][5][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 2 Mathematical Problems in Engineering 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][14][15][16][17][18] and Yeoh [19][20][21][22][23][24][25][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.

Methods and Materials
2.1.Geometry.The geometry of the AF of the L4-L5 segment (Figure 1) 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][29][30].
The AF model has a mesh with eight divisions along the height and eight in the radial direction (Figure 2).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][32][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.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 (Table 1).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).

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.
Also, the influence of the mechanical constants on the maximum stress for the fiber families of the AF was obtained.2)) and their mechanical constants (Table 2) were used for all the models.
where a 1 and a 2 are the material constants and I 4 and I 6 are the deviatoric invariants associated with the two families of fibers, which are defined as I 4 = N (1) .C.N (1)  (3) where  (1) and  (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][39][40][41].
Mooney-Rivlin and Yeoh strain energy functions were used to represent the ground substance.The Mooney-Rivlin function was defined as where c 10 and c 01 are the material constants and I 1 and I 2 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][43][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, c 10 was varied and c 01 was constant, and in the second step, c 10 was constant and c 01 was varied.For instance, for model 3 (M3), c 10 = 0.56 × 50% = 0.28 MPa and c 01 = 0.14 × 100% = 0.14 MPa.In this way, in models M1 to M9, c 10 was varied and c 01 was constant; in models M10 to M18, c 10 was constant and c 01 was varied.
The Yeoh function was defined as where c 1 , c 2 , and c 3 are the mechanical constants of the material and I 1 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, c 1 was varied and c 2 and c 3 were constant; in the second step,  The 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.

Range of Motion (ROM) Behavior.
Compared with the experimental ROM, the mechanical behavior of the Mooney-Rivlin function (Figures 4-7) fitted less well than the behavior obtained with the Yeoh function .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.

Using the Mooney-Rivlin Function as the Strain Energy
Function for the Ground Substance.The ROM behavior was similar for the four movements: as c 10 and c 01 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 was c 10 = 0 MPa and c 01 = 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 when c 01 (0.14 MPa) is constant but c 10 increases and the right sides show the ROM behavior when c 10 (0.56 MPa) is constant but c 01 increases, for all movements.In general, the c 10 variation has a major impact than c 01 over the ROM behavior.

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 using c 1 = 0.0146 MPa, c 2 = -0.0189MPa, and c 3 = 0.0205 MPa; 98.8% for extension movement using c 1 = 0.0292 MPa, c 2 = -0.0189MPa, and c 3 = 0.041 MPa; 88.4% for lateral flexion movement using c 1 = 0.0146 MPa, c 2 = -0.0189MPa, and c 3 = 0.01025 MPa; and 96.0% for axial rotation movement using c 1 = 0.0438 MPa, c 2 = -0.0189MPa, and c 3 = 0.041 MPa .Just four curves were outside the experimental range for flexion (Figure 8), extension (Figure 9), and lateral flexion  In general, the intervertebral disc stiffness increased due to the increase of the mechanical constants; however, a high impact of c 3 on the ROM was found in all movements.This can be explained by the high dispersion area found in the ROM behavior when c 3 varied (curves on the right side, Figures [8][9][10][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.

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 (c 1 , c 2 , and c 3 ).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.

Using the Mooney-Rivlin Function as the Strain Energy
Function for the Ground Substance.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 when c 10 was varied and c 01 was equal to 0.14 MPa.When c 10 was constant (0.56 MPa) and c 01 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 of c 10 and c 01 , 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.

Using the Yeoh Function as the Strain Energy Function
for the Ground Substance.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 when c 1 was varied and c 2 (-0.0189MPa) and c 3 (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, when c 1 (0.0146 MPa) and c 3 (0.041 MPa) were constant and c 2 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, when c 1 (0.0146 MPa) and c 2 (-0.0189MPa) were constant and c 3 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  of c 1 , c 2 , and c 3 , 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.

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 c 10 = 0 and c 01 = 0.14 MPa; here c 01 is equal to the value reported by Heuer et al. [43,47] and Rohlmann et al. [45].However, these researchers reported values of c 10 that differed from zero.On the other hand, the c 01 value is very close to that reported (c 01 = 0.25) by Campbell and Petrella [48].With regard to the stress in the fiber families (flexion: 35.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 c 1 = 0.0292 ± 0.0146 MPa, c 2 = -0.0189±0 MPa, and c 3 = 0.025625 ± 0.015375 MPa.The last values of c 1 , c 2 , and c 3 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 c 1 = 0.0438 MPa, c 2 = -0.0189MPa, and c 3 = 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].

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.

Figure 1 :
Figure 1: Shape and dimensions in mm for the annulus fibrosus of the L4-L5 segment.

Figure 2 :Figure 3 :
Figure 2: Final mesh and sections of the annulus fibrosus used to define their different properties.

c 2
was varied and c 1 and c 3 were constant; finally, in the third step, c 3 was varied and c 1 and c 2 were constant.For instance, in model 22 (M22), c 1 = 0.0146 × 100% = 0.0146 MPa, c 2 = -0.0189× 100% = -0.0189MPa, and c 3 = 0.041 × 75% = 0.03075 MPa.In this way, in models M1 to M9, c1 was varied and c 2 and c 3 were constant; in models M10 to M18, c 2 was varied and c 1 and c 3 were constant; and in models M19 to M27, c 3 was varied and c 1 and c 2 were constant.

Figure 4 :Figure 5 :Figure 6 :
Figure 4: ROM versus Moment for the flexion movement using Mooney-Rivlin function for different values of c 10 and c 01 .

Figure 7 :c1
Figure 7: ROM versus Moment for the axial rotation movement using Mooney-Rivlin function for different values of c 10 and c 01 .

Figure 8 :
Figure 8: ROM versus Moment for the flexion movement using Yeoh function for different values of c 1 , c 2 , and c 3 .

Figure 9 :Figure 10 :Figure 11 :
Figure 9: ROM versus Moment for the extension movement using Yeoh function for different values of c 1 , c 2, and c 3 .

Figure 12 :
Figure 12: Maximum stress in the fibers family of AF versus mechanical constants values for the flexion movement, using Mooney-Rivlin function.

Figure 13 :
Figure 13: Maximum stress in the fibers family of AF versus mechanical constants values for the extension movement, using Mooney-Rivlin function.

Figure 14 :
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 :
Figure 15: Maximum stress in the fibers family of AF versus mechanical constants values for the axial rotation movement, using Mooney-Rivlin function.

Figure 16 :
Figure 16: Maximum stress in the fibers family of AF versus mechanical constants values for the flexion movement, using Yeoh function.

Figure 17 :Figure 18 :
Figure 17: Maximum stress in the fibers family of AF versus mechanical constants values for the extension movement, using Yeoh function.

Figure 19 :
Figure 19: Maximum stress in the fibers family of AF versus mechanical constants values for the axial rotation movement, using Yeoh function.

Table 1 :
Type and number of elements used in the finite element model.

Table 2 :
Mechanical constants used to the fibers of the annulus fibrous.

Table 3 :
Mechanical constants for the strain energy function Mooney-Rivlin.

Table 4 :
Mechanical constants used for the Yeoh strain energy function.