^{1}Mechanical, Industrial and Manufacturing Engineering Department, University of Toledo, Toledo, OH, USA^{2}Engineering Center for Orthopaedic Research Excellence, University of Toledo, Toledo, OH, USA^{3}Bioengineering Department, University of Toledo, Toledo, OH, USA

Besides the biology, stresses and strains within the tissue greatly influence the location of damage initiation and mode of failure in an intervertebral disk. Finite element models of a functional spinal unit (FSU) that incorporate reasonably accurate geometry and appropriate material properties are suitable to investigate such issues. Different material models and techniques have been used to model the anisotropic annulus fibrosus, but the abilities of these models to predict damage initiation in the annulus and to explain clinically observed phenomena are unclear. In this study, a hyperelastic anisotropic material model for the annulus with two different sets of material constants, experimentally determined using uniaxial and biaxial loading conditions, were incorporated in a 3D finite element model of a ligamentous FSU. The purpose of the study was to highlight the biomechanical differences (e.g., intradiscal pressure, motion, forces, stresses, strains, etc.) due to the dissimilarity between the two sets of material properties (uniaxial and biaxial). Based on the analyses, the biaxial constants simulations resulted in better agreements with the *in vitro* and *in vivo* data, and thus are more suitable for future damage analysis and failure prediction of the annulus under complex multiaxial loading conditions.

## Introduction

Intervertebral disks are interposed between the two adjoining vertebral bodies along the spine. They impart stability and flexibility to the human spine. A disk comprises three different components: annulus fibrosis (AF), nucleus pulposus (NP), and cartilaginous endplate (EP). NP is the central part of disk enclosed in the annulus and bonded to superior and inferior cartilaginous EPs. The water content of the hydrated NP is 90% at birth. It decreases with age to 80% at 20 years and 70% by 60 years and beyond as a part of the aging process (Iatridis et al., 1996). The change in the water concentration with age may lead to the disk degeneration.

Intervertebral disks are one of the primary sources of acute low back pain because of disk degeneration and disk herniation (Kuslich et al., 1991). However, the underlying mechanisms for the disk degeneration are still unclear, especially the relationship between the external loads and damage to disk annulus. Magnetic resonance images (MRI) show annular tears, mostly in the lumbar intervertebral disks and especially at the L4–L5 level (Rajasekaran et al., 2013).

Three types of tears are prevalent: radial, circumferential, and the rim lesions (Osti et al., 1992; Vernon-Roberts et al., 2007). These tears or cracks can be produced by combinations of several loads experimentally. For example, an *in vitro* model with physiological stress produced under a cyclic combination of compression, rotation, and flexion supports that annular protrusion occurs in combined loading condition, and AF is the site of primary pathologic changes (Gordon et al., 1991). In another experiment, the effect of combination of flexion and torsion on the ovine lumbar disk herniation has been studied. The results of this study indicate that this combination effectively increased the appearance of radial tear in the AF (Veres et al., 2010).

To investigate the initiation location and then the propagation of cracks in the intervertebral disk, we need to have a good understanding of the stresses and strains in the intervertebral disk. Finite element analysis (FEA) is one of the most effective tools to evaluate the internal response of the disk due to external loads, since it is difficult, if not impossible, to experimentally measure stresses under several loading conditions. Stress and strain distributions obtained from such analysis depend on the geometry, material model, and material properties used in the simulations. For example, modeling geometry varies from disk alone to full ligamentous functional spinal unit (FSU) simulations (Natali, 1991; Kumaresan et al., 1999; Rohlmann et al., 2006; Schroeder et al., 2010). Different material models and techniques are also used to simulate each component of the intervertebral disk. Nucleus varies from linear and non-linear elastic to fluid in the simulations (Wang et al., 1997; Yao et al., 2006). Different techniques are also used for AF modeling.

Layer by layer investigation of the human AF structure using microscopic techniques has shown that AF consists of 15–25 distinct layers, depending on spine level and age (Marchand and Ahmed, 1990). The collagen fibers have two defined orientations that change from one layer to the next. The typical angle of fibers has been reported to have a wide range (Marchand and Ahmed, 1990), with an average of ±30° with respect to the transverse plane implemented in some previous studies. Tensile properties of the human AF from single and multiple lamellae samples suggest non-linear and anisotropic behavior of this tissue (Green et al., 1993; Skaggs et al., 1994; Ebara et al., 1996; Fujita et al., 1997; Holzapfel et al., 2005).

Founded on the AF complex structure, different material models have been proposed to represent this structure. Mixture of fibers and matrix in a composite structure with linear and non-linear properties of fibers and matrix has been a commonly used modeling technique (Shirazi-Adl et al., 1984; Lu et al., 1996; Wang et al., 1997; Ferguson and Steffen, 2003; Schmidt et al., 2006; Little et al., 2007). In addition, the roles of porosity, permeability (Natarajan et al., 2007; Galbusera et al., 2011; Qasim et al., 2012), or viscoelasticity (Schroeder et al., 2006) have also been investigated. In these models, fibers are modeled using rebar or spring elements and anisotropy is introduced by specifying their orientation. Modeling of such composite structures is, however, highly dependent on model parameters, which are often difficult to determine (for instance, fiber and matrix properties, or fiber volume fraction). These models have extensively been used for spine kinematics under various loading conditions (Rohlmann et al., 2009; Dreischarf et al., 2014). Failure analysis based on these classical models has been conducted using ultimate tensile strain in the fibers (Shirazi-Adl, 1989) and also stresses in the ground substance (Qasim et al., 2012).

A model, which adequately describes disk mechanics and stress distribution under multiaxial loading, is essential for annulus damage and failure analysis. In the present study, because of paucity of a clear picture of clinically observed annulus damage vis-à-vis the mechanical loads, we have used a hyperelastic anisotropic material model for the AF to investigate the biomechanical behavior of the ligamentous L4–L5 segment. This material model was initially proposed for arterial layers (Holzapfel et al., 2000) and is based on the strain energy function described in Spencer (1984). Later, this material model has been used by various investigators for different purposes such as modeling intervertebral disk (Eberlein et al., 2001). In this material model, hereafter referred to as Holzapfel–Gasser–Ogden (HGO) model, AF is presented as a continuum material and the effects of fibers and matrix are considered in the strain energy. Strain distributions from hyperelastic anisotropic material model of AF were compared to a model based on classical method. Unexpected strain concentrations in the posterior part of the matrix were observed in the classical material model based on a composite structure with fibers and matrix, but not in the HGO material model (Eberlein et al., 2001). These concentrations were explained to be owing to two deficiencies of the classical material model: matrix considered as linear elastic, and fibers considered in short-fiber arrangement.

The hyperelastic anisotropic continuum model has been used to study the kinematics of the lumbar spine (Del Palomar et al., 2008; Moramarco et al., 2010). However, uniaxial properties of the AF were used in those studies. By contrast, *in situ* conditions that AF experiences is a biaxial stress state. Therefore, based on experimental data for this material model obtained from uniaxial and biaxial tests, the aim of this paper is to highlight the biomechanical differences [e.g., intradiscal pressure (IDP), motion, and stresses] resulting from dissimilarities between the two sets of material properties (uniaxial and biaxial). It is shown that the material properties from biaxial tests simulate the annulus fibrosus with hyperelastic anisotropic material model more accurately than properties from uniaxial tests and are more relevant to the *in situ* condition. A discussion and analysis of model properties relevant to understanding clinically relevant disk mechanics under different loading conditions is also presented. A model, which adequately describes disk mechanics and stress distribution under such complex loading, is essential for AF damage and failure analysis.

## Materials and Methods

### Finite Element Model

A previously valid three-dimensional (3D) FE model of ligamentous L4–L5 FSU (Goel et al., 1995) was used. This model is one of the well-established models, which has been used for modeling lumbar spine (Dreischarf et al., 2014). The geometry was developed from computer tomography (CT) scans. The slice thickness was approximately 1.5 mm and the thickness of the cortical layer was 0.5 mm across the model. The model was validated using *in vitro* kinematic data, facet loads, ligaments strains, and disk bulge under several loading conditions (Goel et al., 1995, 2007; Dooris et al., 2001). The model consists of cortical and cancellous vertebral bones, posterior bony elements, annulus, nucleus, facet joints, and seven major ligaments, Figure 1A. This model contains 15,136 elements and 19,148 nodes and is symmetric with respect to the mid-sagittal plane (Monore, 1993; Kiapour et al., 2012; Dreischarf et al., 2014). For bony structure and intervertebral disk components, C3D8 brick elements were used. All of the ligaments were modeled using two nodes tension-only truss elements (T3D2) with non-linear material properties.

**Figure 1. (A)** Ligaments in the FE model of the functional spinal unit, **(B)** FE model of intervertebral 539 disc consisting of AF and NP.

Intervertebral disk was divided in two separate anatomical parts: AF and NP, the same as the previous model, Figure 1B. Linear isotropic material constitutive relationship was used in NP region and a hydrostatic pressure was simulated. Several approaches have been used to simulate the incompressibility of the nucleus, including the one used in the present study. All of the approaches have yielded clinically relevant results. This approach for modeling has been well established. We have also used a hyperelastic neo-Hookean model for the annulus and incompressible fluid simulation for the nucleus. The results were similar. The range of elasticity modulus that has been considered for NP in the literature is within 1.5–10 MPa (Wang et al., 2000; Chen et al., 2001, 2008; Ayturk and Puttlitz, 2011) and here it was assumed to be 2 MPa. AF complex structure was modeled by the hyperelastic and anisotropic material model of HGO, which is available in ABAQUS/Standard™ version 6.10 [Simulia Inc, Rhode Island-USA]. The facet joints had a gap of 0.5 mm and could only transmit compressive forces. Facet contact stiffness was modeled with a non-linear exponential function with GAPUNI elements. The elastic properties of the different parts are listed in Table 1.

**Table 1**. **Material properties and element type of bony structures, ligaments, intervertebral disk, and facet joint (Agarwal et al., 2013a,b)**.

### Annulus Material Model

The annulus was modeled in a different manner, as opposed to fibers embedded in a ground substance in the classical model, as explained in the Section “Introduction.” The HGO material model is based on the strain energy potential for distributed collagen fibers in the ground substance (Holzapfel et al., 2000; Gasser et al., 2006). The primary advantage of using this material model rather than phenomenological models is that the fibers and matrix material properties are associated with the material constituents (histological structure). In addition, the material parameters for the model can be experimentally determined, as described earlier.

In this material model, exponential material behavior for fibers and non-linear hyperelastic neo-Hookean isotropic material model for the ground substance were used. Unlike previous modeling that considered discontinuous arrangement for annulus fibers, in HGO material model collagen are arranged as continuous fibers. These models for the matrix and fiber are represented by Eqs 1 and 2, respectively:

with:

where, ${\overline{I}}_{1}$ is the first deviatoric strain invariant, *J* is the elastic volume ratio, and ${\overline{I}}_{4\left(\text{\alpha \alpha}\right)}$ is pseudo-invariants of the distortion part of the right Cauchy–Green strain and unit vectors in the direction of fiber families. Matrix compressibility and matrix stiffness are defined by *D* and *C*_{10} in Equation 1, respectively, *K*_{1} > 0 in Eq. 2 is a material parameter with a dimension of stress and relates to the stiffness of fibers. *K*_{2} > 0 in this equation is a dimensionless material property that is related to fiber non-linear behavior. Fibers in this structural material model can support only tensile stresses, therefore, Eq. 2 applies in fiber extension mode (Gasser et al., 2002). *N* is the number of fiber families in Eq. 2. In the AF, there are two fiber families with ±30° orientation with respect to the horizontal plane. Value of κ in Eq. 3 is between ${\scriptscriptstyle \frac{1}{3}}$ for the randomly oriented fibers and 0 for aligned fibers. Here, we considered the two fiber families of AF to be aligned fibers for simplicity.

For incompressible materials, *D* is assumed to approach infinity. Because we can assume the matrix as an incompressible material (Natali, 1991), the second terms of the matrix formulation are eliminated. Only three material constants (*C*_{10}, *K*_{1}, and *K*_{2}) are needed. The equation for strain energy potential will then simplify as Eq. 4. The values of these three constants for both uniaxial tension (O’Connell et al., 2009) and biaxial tension (O’Connell et al., 2012) stress states are given in Table 2.

It should be noted that while hyperelastic neo-Hookean isotropic material with the associated material constant *C*_{10} is used as the material model for the ground substance, the material model used by O’Connell et al. (2012) from which the value of *C*_{10} was taken is the Mooney–Rivlin model. However, with assuming an incompressible material as mentioned previously, the value of constant *C*_{10} is identical in both the neo-Hookean and Mooney–Rivlin models.

As can be observed from Table 2, there is a significant difference between the biaxial parameters and uniaxial parameters. This is mentioned in the literature as well (Gregory and Callaghan, 2011). In this work, we used both sets of parameters in two separate 3D finite element models of FSU that were previously described. It was desired to evaluate and compare FSU kinematics and stress behavior in response to this change to determine which set of material properties produce more accurate results and are more suitable for damage analysis in future.

### Boundary and Loading Conditions

The inferior surface of the L5 vertebra was rigidly fixed. Compressive follower loads were applied up to 3400 N. Using follower loads in a motion segment, the compressive load is tangent to the spinal curve and, therefore, intervertebral disk would be loaded in almost pure compression. The predicted average load vs. axial displacements in the anterior and posterior regions of the disk were compared to the *in vitro* experimental data in which pure axial compression loads were applied only to the intervertebral disk (Markolf and Morris, 1974). Therefore, both FE analysis and experiments would produce the same loading scenario in the intervertebral disk (pure compression) for comparison.

The annulus axial stresses along antero-posterior direction at the middle disk plane for 2000 N were compared to the experimental data (Adams et al., 1996). They also indicated that extension caused the apophyseal joints to become load-bearing, and damage could occur at compressive loads as low as 500 N (Adams et al., 1994). The models were also subjected to compressive forces of 300, 700, 850, and 3400 N simulating equivalent compressive loads on the supine, standing, walking, and lifting 44 lb with back bent and knees straight positions, respectively (Nachemson, 1966, 1975, 1991).

The IDP was defined as the hydrostatic pressure [−1/3*tr(σ), where tr(σ) is the first invariant of the stress tensor σ] in the center of the nucleus, and was compared with available *in vivo* results (Wilke et al., 1999). However, the results of these studies are presented for different spinal levels and may not be directly comparable. Nachemson (1966) introduced a correction factor for estimation of the compressive force. The problem of the estimation of compressive loads has been discussed in Dreischarf et al. (2013). Pure moments in different directions on the L4–L5 level were also applied. The axial displacement, flexion, extension, lateral bending, and axial rotation values for the aforementioned loads/moments were computed for comparison with data reported in the literature and evaluation of the hypothesis.

The range of motion for different loading conditions for the two FE models used in this study are compared to experimental data from Agarwal et al. (2013c). For these experiments, spine segments with ligaments attached were used. Each specimen was subjected to pure moments using rods and pulleys to simulate flexion, extension, lateral bending, and axial rotation. The spatial coordinate data obtained were used to compute the 3D segmental rotations.

## Results

The mean and range of load vs. axial displacement for an *in vitro* experiment (Markolf and Morris, 1974) on healthy intervertebral disk only, and the results of current FE study of both models with uniaxial and biaxial material properties are presented in Figure 2. The predicted load-axial displacement of the model with biaxial material properties is between the experimental ranges. However, load-axial displacement for the uniaxial model is below the lower limit of the experimental data. Range of motion for different loading conditions for the two FE models are compared to data from Agarwal et al. (2013c) in Figure 3 and indicate the results from the FE model based on biaxial properties are closer to the experimental data.

**Figure 2. Comparison of load-displacement curves between uniaxial and biaxial models and as compared to experimental results from Markolf and Morris (1974) for compressive force applied as follower load to the L4–L5 level**. Dash lines indicate the experimental range.

**Figure 3. Comparison of range of motion between models with uniaxial and biaxial material properties with experimental results from Agarwal et al. (2013c) at 10 Nm pure moment in several directions**.

The outcomes shown in Figure 3 indicate that the model with uniaxial material properties resulted in higher values in extension by 57%, flexion by 30%, axial rotation by 82%, and lateral bending by 125%, compared to the model with biaxial material properties. In flexion, extension, and lateral bending, the results of the predicted biaxial model based data were within the SD of the experimental data. Axial rotation data for the biaxial model were close to the lower limits of the experimental data. However, the corresponding data for the uniaxial material properties for all loading conditions exceeded the upper limit of the experimental data, Figure 3.

Corresponding IDPs for the biaxial and uniaxial models vs. the *in vivo* experimental data for various compressive forces are plotted in Figure 4. The difference between these two models (uniaxial vs. biaxial) was <10%. Both of the models show close results to the experimental data (Wilke et al., 1999). Figure 5 shows a more detailed picture of IDP for biaxial-based FE model vs. applied loads and as compared to experimental data with the same loading conditions (Adams et al., 1994).

**Figure 4. Comparison of intradiscal pressure (IDP) between simulations with uniaxial and biaxial material property models as compared with mean values of experimental results (Wilke et al., 1999) under the same loading conditions**.

**Figure 5. Comparison of intradiscal pressure (IDP) with FE model with biaxial material properties and experimental results from Adams et al. (1994)**.

Axial stress distribution along the anterior–posterior direction (Figure 6A) in the middle disk height plane for 2000 N compression in the AF is shown in Figure 6B. The difference between the model results with biaxial properties and experimental data (McNally and Adams, 1992) is 25%. However, this difference in the uniaxial model is 60%. Thus, the model based on the biaxial material properties in both cases exhibited closer trend to the experimental results, as compared with the model based on uniaxial properties.

**Figure 6. (A)** Different anatomical regions in the FE model of AF. **(B)** Axial stress distribution obtained from biaxial and uniaxial FE models in the mid-height plane of intervertebral disk from posterior to anterior midline and comparison with the experimental data from McNally and Adams (1992) at 2000 N compressive force.

To compare the results in this work with an example of a study using a continuum damage mechanics methodology based on a composite model with fibers and matrix, the results from Qasim et al. (2012) were considered. Our prediction based on biaxial model with no refinement for porosity and other parameters are in agreement with the results from this study, Figure 7.

**Figure 7. Comparison of range of motion between the model with biaxial material properties and the detailed FE model from Qasim et al. (2012) as well as with experimental results from Fujiwara et al. (2000) for 6.6 Nm moment**.

## Discussion

The aim of the current study was to use an anisotropic material model for the AF region and compare the predicted FSU biomechanics based on using uniaxial or biaxial material properties assigned to the AF. The predicted data were also compared to relevant experimental data from the literature. These experiments were among the most accepted and cited experiments in the field.

The load vs. axial displacement results obtained from biaxial material properties showed good agreement with experimental data (Figure 2). Uniaxial properties led to more axial deformation, in comparison to the biaxial material properties. This lower stiffness is also observed in the angular rotation under applied moments in several directions (Figure 3). The higher stiffness of the motion segment with biaxial material properties of AF led to higher axial stresses in the AF, much closer to the experimental data (Figure 6B). This indicates that biaxial material properties represent a stress state closer to the *in situ* condition. In addition, *in situ* AF fibers are attached to the upper and lower vertebral bodies, which provide a constrained condition. Therefore, the stress state of AF *in situ* condition is in agreement with biaxial tests. This condition has also been suggested in similar studies on other soft tissues such as skin, when subjected to *in vivo* multi axial stress state (Tong and Fung, 1976), and arteries (Debes and Fung, 1995).

Non-linear behavior of the lumbar segment under axial compressive force, which has been observed in experiments (Kulak et al., 1976) was predicted by non-linear material model used for AF region for both the uniaxial and biaxial material property models. Experiments have shown that IDP has a linear relationship in pure compression force (Adams et al., 1994), and our predictions are in agreement for both uniaxial- and biaxial-based models, Figure 4.

Instead of considering separate stresses in the matrix and fibers, a combined nominal stress as obtained in this study for AF gives the opportunity to compare the model results directly with experimental data. Axial stress distribution resulting from compressive force leads to a peak stress on the posterior region of the disk according to the model results, Figure 6B. Posterior region is the area that has the most likelihood of damage, according to clinical data (Vernon-Roberts et al., 2007). Tensile stresses in outer annulus regions observed in Figure 6B may be regarded as contribution of the annulus fibers to the total stress tensor, which is not measurable experimentally.

It was shown that incorporating biaxial material properties in a hyperelastic anisotropic material model result in internal stresses in the intervertebral disk that are in agreement with the *in vitro* data, as compared with the properties obtained from uniaxial tests. The magnitude of stresses has a dominant effect on damage evaluation and failure prediction of AF. Further work is being pursued to document damage prediction using the biaxial material property model.

To develop uniaxial- and biaxial-based models, experimental data from other studies were used and their limitations, highlighted in those studies, will be reflected in our model as well. In these experimental studies, the properties were calculated for outer anterior site of AF (O’Connell et al., 2009, 2012). Therefore, we used the same data for the model used in this current study as a homogenous model. In the experimental studies, degeneration effects have been considered as well. It has been shown that for biaxial stress state matrix parameters (*C*_{10} and *D*) and fiber stiffness (*K*_{1}) in Eq. 4 do not correlate with degeneration, whereas fiber non-linearity (*K*_{2}) decreases with increasing grade of degeneration (O’Connell et al., 2012). However, for uniaxial stress state matrix parameters contribute to represent the effect of degeneration (O’Connell et al., 2009). In the current study, we focused on the healthy AF properties, but the models can be modified to account for changes in parameters as a function of degeneration.

Qasim et al. (2012) refined a basic model by adding properties like porosity, strain dependent permeability, and osmotic pressure into their model. Agreement of the results in this study without refinement for porosity and other parameters with those from Qasim et al. (2012), Figure 7, suggests that effects of porosity and other parameters such as water content may not significantly change the findings due to the static and short-time loading conditions that were considered in this study. For more robust model predictions and for long-time loading conditions, however, these refinements should be included to evaluate their effects, as they contribute to viscous behavior. The loading rate has also been shown to greatly influence mechanical response of spinal structures to external loads (Ochia et al., 2003). Therefore, viscoelasticity, which is a characteristic of soft tissues, can affect stress and strain distributions and is important for damage analyses and failure prediction. Although incorporation of this effect into material modeling is very complex, it may need to be considered for higher levels of accuracy.

The annulus region was considered as a homogeneous material in term of fiber orientation and density. Even though the local variations of the fiber orientation were not considered, the model based on biaxial material properties predicts stresses close to those measured experimentally. For more accurate analyses, however, it is necessary to consider this heterogeneous property of annulus.

## Conclusion

This study highlighted the biomechanical differences resulting from dissimilarities between two sets of material properties, uniaxial and biaxial. It was shown that the material properties from biaxial tests simulate the annulus fibrosus behavior more closely, as compared with the experimental data, than properties from uniaxial tests. The predicted load-axial displacement curve of the model with biaxial material properties was shown to be between the experimental ranges, while it is below the lower limit of the experimental data for the uniaxial model. Range of motion for different loading conditions (compression, extension, flexion) from the FE model based on biaxial properties was also shown to be closer to the experimental data. Under a compressive force, axial stress distribution based on the biaxial material properties exhibited closer trend to the experimental results. It is, therefore, concluded that simulation of AF mechanical behavior with biaxial material properties may yield better predictions than with uniaxial material properties under different loading conditions.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## References

Adams, M. A., McNally, D. S., Chinn, H., and Dolan, P. (1994). The clinical biomechanics award paper 1993 posture and the compressive strength of the lumbar spine. *Clin. Biomech.* 9, 5–14. doi: 10.1016/0268-0033(94)90052-3

Adams, M. A., McNally, D. S., and Dolan, P. (1996). “Stress” distributions inside intervertebral discs the effects of age and degeneration. *J. Bone Joint Surg. Br.* 78, 965–972. doi:10.1302/0301-620X78B6.1287

Agarwal, A., Agarwal, A. K., and Goel, V. K. (2013a). The endplate morphology changes with change in biomechanical environment following discectomy. *Int. Clin. Med.* 4, 8–17. doi:10.4236/ijcm.2013.47A1002

Agarwal, A., Palepu, V., Agarwal, A. K., Goel, V. K., and Yildirim, E. D. (2013b). Biomechanical evaluation of an endplate-conformed polycaprolactone-hydroxyapatite intervertebral fusion graft and its comparison with a typical nonconformed cortical graft. *J. Biomech. Eng.* 135, 61005–61009. doi:10.1115/1.4023988

Agarwal, A. K., Kodigudla, M., Desai, D., Agarwal, A., Palepu, V., and Goel, V. K. (2013c). *ISASS13 – Regular Posters – Lumbar Therapies and Outcomes – 396 – PEEK Rod and Ti Screw Fixation Provides Comparable Construct Stability*. Vancouver: International Society for the Advancement of Spine Surgery-New York.

Ayturk, U. M., and Puttlitz, C. M. (2011). Parametric convergence sensitivity and validation of a finite element model of the human lumbar spine. *Comput. Methods Biomech. Biomed. Engin.* 14, 695–705. doi:10.1080/10255842.2010.493517

Chen, C. S., Cheng, C. K., Liu, C. L., and Lo, W. H. (2001). Stress analysis of the disc adjacent to interbody fusion in lumbar spine. *Med. Eng. Phys.* 23, 485–493. doi:10.1016/S1350-4533(01)00076-5

Chen, S. H., Tai, C. L., Lin, C. Y., Hsieh, P. H., and Chen, W. P. (2008). Biomechanical comparison of a new stand-alone anterior lumbar interbody fusion cage with established fixation techniques: a three-dimensional finite element analysis. *BMC Musculoskelet. Disord.* 9:88. doi:10.1186/1471-2474-9-88

Debes, J. C., and Fung, Y. C. (1995). Biaxial mechanics of excised canine pulmonary arteries. *Am. J. Physiol.* 38, H433–H442.

Del Palomar, A. P., Calvo, B., and Doblaré, M. (2008). An accurate finite element model of the cervical spine under quasi-static loading. *J. Biomech.* 41, 523–531. doi:10.1016/j.jbiomech.2007.10.012

Dooris, A. P., Goel, V. K., Grosland, N. M., Gilbertson, L. G., and Wilder, D. G. (2001). Load-sharing between anterior and posterior elements in a lumbar motion segment implanted with an artificial disc. *Spine* 26, E122–E129. doi:10.1097/00007632-200103150-00004

Dreischarf, M., Rohlmann, A., Zhu, R., Schmidt, H., and Zander, T. (2013). Is it possible to estimate the compressive force in the lumbar spine from intradiscal pressure measurements? A finite element evaluation. *Med. Eng. Phys.* 35, 1385–1390. doi:10.1016/j.medengphy.2013.03.007

Dreischarf, M., Zander, T., Shirazi-Adl, A., Puttlitz, C. M., Adam, C. J., Chen, C. S., et al. (2014). Comparison of eight published static finite element models of the intact lumbar spine: predictive power of models improves when combined together. *J. Biomech.* 47, 1757–1766. doi:10.1016/j.jbiomech.2014.04.002

Ebara, S., Iatridis, J. C., Setton, L. A., Foster, R. J., Mow, V. C., and Weidenbaum, M. (1996). Tensile properties of nondegenerate human lumbar anulus fibrosus. *Spine* 21, 452–461. doi:10.1097/00007632-199602150-00009

Eberlein, R., Holzapfel, G. A., and Schulze-Bauer, C. A. J. (2001). An anisotropic model for annulus tissue and enhanced finite element analyses of intact lumbar disc bodies. *Comput. Methods Biomech. Biomed. Engin.* 4, 209–229. doi:10.1080/10255840108908005

Ferguson, S. J., and Steffen, T. (2003). Biomechanics of the aging spine. *Eur. Spine J.* 12, S97–S103. doi:10.1007/s00586-003-0621-0

Fujita, Y., Duncan, N. A., and Lotz, J. C. (1997). Radial tensile properties of the lumbar annulus fibrosus are site and degeneration dependent. *J. Orthop. Res.* 15, 814–819. doi:10.1002/jor.1100150605

Fujiwara, A., Lim, T. H., An, H. S., Tanaka, N., Jeon, C. H., Andersson, G. B., et al. (2000). The effect of disc degeneration and facet joint osteoarthritis on the segmental flexibility of the lumbar spine. *Spine* 25, 3036–3044.

Galbusera, F., Schmidt, H., Noailly, J., Malandrino, A., Lacroix, D., Wilke, H. J., et al. (2011). Comparison of four methods to simulate swelling in poroelastic finite element models of intervertebral discs. *J. Mech. Behav. Biomed. Mater.* 4, 1234–1241. doi:10.1016/j.jmbbm.2011.04.008

Gasser, T. C., Ogden, R. W., and Holzapfel, G. A. (2006). Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. *J. R. Soc. Interface* 3, 15–35. doi:10.1098/rsif.2005.0073

Gasser, T. C., Schulze-Bauer, C. A. J., and Holzapfel, G. A. (2002). A three-dimensional finite element model for arterial clamping. *J. Biomech. Eng.* 124, 355–363. doi:10.1115/1.1485284

Goel, V. K., Mehta, A., Jangra, J., Faizan, A., Kiapour, A., Hoy, R. W., et al. (2007). Anatomic facet replacement system (AFRS) restoration of lumbar segment mechanics to intact: a finite element study and in vitro cadaver investigation. *SAS J.* 1, 46–54. doi:10.1016/SASJ-2006-0010-RR

Goel, V. K., Monroe, B. T., Gilbertson, L. G., and Brinckmann, P. (1995). Interlaminar shear stresses and laminae separation in a disc. Finite element analysis of the L3-L4 motion segment subjected to axial compressive loads. *Spine* 20, 689–698. doi:10.1097/00007632-199503150-00010

Gordon, S. J., Yang, K. H., Mayer, P. J., Mace, A. H. Jr., Kish, V. L., and Radin, E. L. (1991). Mechanism of disc rupture: a preliminary report. *Spine* 16, 450–456. doi:10.1097/00007632-199104000-00011

Green, T. P., Adams, M. A., and Dolan, P. (1993). Tensile properties of the annulus fibrosus II. Ultimate tensile strength and fatigue life. *Eur. Spine J.* 2, 209–214. doi:10.1007/BF00299448

Gregory, D. E., and Callaghan, J. P. (2011). A comparison of uniaxial and biaxial mechanical properties of the annulus fibrosus: a porcine model. *J. Biomech. Eng.* 133, 024503. doi:10.1115/1.4003327

Holzapfel, G. A., Gasser, T. C., and Ogden, R. W. (2000). A new constitutive framework for arterial wall mechanics and a comparative study of material models. *J. Elasticity Phys. Sci. Solids* 61, 1–48. doi:10.1023/A:1010835316564

Holzapfel, G. A., Schulze-Bauer, C. A. J., Feigl, G., and Regitnig, P. (2005). Single lamellar mechanics of the human lumbar anulus fibrosus. *Biomech. Model. Mechanobiol.* 3, 125–140. doi:10.1007/s10237-004-0053-8

Iatridis, J. C., Weidenbaum, M., Setton, L. A., and Mow, V. C. (1996). Is the nucleus pulposus a solid or a fluid? Mechanical behaviors of the nucleus pulposus of the human intervertebral disc. *Spine* 21, 1174–1184. doi:10.1097/00007632-199605150-00009

Kiapour, A., Ambati, D., Hoy, R. W., and Goel, V. K. (2012). Effect of graded facetectomy on biomechanics of Dynesys dynamic stabilization system. *Spine* 37, E581–E589. doi:10.1097/BRS.0b013e3182463775

Kulak, R. F., Belytschko, T. B., Schultz, A. B., and Galante, J. O. (1976). Nonlinear behavior of the human intervertebral disc under axial load. *J. Biomech.* 9, 377–386. doi:10.1016/0021-9290(76)90115-9

Kumaresan, S., Yoganandan, N., Pintar, F. A., and Maiman, D. J. (1999). Finite element modeling of the cervical spine: role of intervertebral disc under axial and eccentric loads. *Med. Eng. Phys.* 21, 689–700. doi:10.1016/S1350-4533(00)00002-3

Kuslich, S., Ulstrom, C., and Michael, C. (1991). The tissue origin of low back pain and sciatica: a report of pain response to tissue stimulation during operations on the lumbar spine using local anesthesia. *Orthop. Clin. North Am.* 22, 181–187.

Little, J. P., Adam, C. J., Evans, J. H., Pettet, G. J., and Pearcy, M. J. (2007). Nonlinear finite element analysis of anular lesions in the L4/5 intervertebral disc. *J. Biomech.* 40, 2744–2751. doi:10.1016/j.jbiomech.2007.01.007

Lu, Y. M., Hutton, W. C., and Gharpuray, V. M. (1996). Do bending, twisting, and diurnal fluid changes in the disc affect the propensity to prolapse? A viscoelastic finite element model. *Spine* 21, 2570–2579. doi:10.1097/00007632-199611150-00006

Marchand, F., and Ahmed, A. M. (1990). Investigation of the laminate structure of lumbar disc anulus fibrosus. *Spine* 15, 402–410. doi:10.1097/00007632-199005000-00011

Markolf, K. L., and Morris, J. M. (1974). The structural components of the intervertebral disc: a study of their contribution to the ability of the disc to withstand compressive forces. *J. Bone Joint Surg. Am.* 56, 675–687.

McNally, D. S., and Adams, M. A. (1992). Internal intervertebral disc mechanics as revealed by stress profilometry. *Spine* 17, 66–73. doi:10.1097/00007632-199201000-00011

Monore, B. T. (1993). *Interlaminar Shear Stresses as a Cause of Disc Degeneration with Age: A Non-Linear Composite Finite Element Analysis*. M.S. Thesis, Biomedical Engineering, University of Iowa, Iowa, IA.

Moramarco, V., Pérez del Palomar, A., Pappalettere, C., and Doblaré, M. (2010). An accurate validation of a computational model of a human lumbosacral segment. *J. Biomech.* 43, 334–342. doi:10.1016/j.jbiomech.2009.07.042

Nachemson, A. (1966). The load on lumbar disks in different positions of the body. *Clin. Orthop. Relat. Res.* 45, 107–122. doi:10.1097/00003086-196600450-00014

Nachemson, A. (1975). Towards a better understanding of low-back pain: a review of the mechanics of the lumbar disc. *Rheumatology* 14, 129–143. doi:10.1093/rheumatology/14.3.129

Nachemson, A. (1991). Spinal disorders: overall impact on society and the need for orthopedic resources. *Acta Orthop.* 62, 17–22. doi:10.3109/17453679109155099

Natali, A. N. (1991). A hyperelastic and almost incompressible material model as an approach to intervertebral disc analysis. *J. Biomed. Eng.* 13, 163–168. doi:10.1016/0141-5425(91)90063-D

Natarajan, R. N., Williams, J. R., Lavender, S. A., and Andersson, G. B. J. (2007). Poro-elastic finite element model to predict the failure progression in a lumbar disc due to cyclic loading. *Comput. Struct.* 85, 1142–1151. doi:10.1016/j.compstruc.2006.08.043

Ochia, R. S., Tencer, A. F., and Ching, R. P. (2003). Effect of loading rate on endplate and vertebral body strength in human lumbar vertebrae. *J. Biomech.* 36, 1875–1881. doi:10.1016/S0021-9290(03)00211-2

O’Connell, G. D., Guerin, H. L., and Elliott, D. M. (2009). Theoretical and uniaxial experimental evaluation of human annulus fibrosus degeneration. *J. Biomech. Eng.* 131, 111007. doi:10.1115/1.3212104

O’Connell, G. D., Sen, S., and Elliott, D. M. (2012). Human annulus fibrosus material properties from biaxial testing and constitutive modeling are altered with degeneration. *Biomech. Model. Mechanobiol.* 11, 493–503. doi:10.1007/s10237-011-0328-9

Osti, O. L., Vernon-Roberts, B., Moore, R., and Fraser, R. D. (1992). Annular tears and disc degeneration in the lumbar spine: a post-mortem study of 135 discs. *J. Bone Joint Surg. Br.* 74, 678–682.

Qasim, M., Natarajan, R. N., An, H. S., and Andersson, G. B. J. (2012). Initiation and progression of mechanical damage in the intervertebral disc under cyclic loading using continuum damage mechanics methodology: a finite element study. *J. Biomech.* 45, 1934–1940. doi:10.1016/j.jbiomech.2012.05.022

Rajasekaran, S., Bajaj, N., Tubaki, V., Kanna, R. M., and Shetty, A. P. (2013). ISSLS prize winner: the anatomy of failure in lumbar disc herniation: an in vivo, multimodal, prospective study of 181 subjects. *Spine* 38, 1491–1500. doi:10.1097/BRS.0b013e31829a6fa6

Rohlmann, A., Zander, T., Rao, M., and Bergmann, G. (2009). Applying a follower load delivers realistic results for simulating standing. *J. Biomech.* 42, 1520–1526. doi:10.1016/j.jbiomech.2009.03.048

Rohlmann, A., Zander, T., Schmidt, H., Wilke, H. J., and Bergmann, G. (2006). Analysis of the influence of disc degeneration on the mechanical behaviour of a lumbar motion segment using the finite element method. *J. Biomech.* 39, 2484–2490. doi:10.1016/j.jbiomech.2005.07.026

Schmidt, H., Heuer, F., Simon, U., Kettler, A., Rohlmann, A., Claes, L., et al. (2006). Application of a new calibration method for a three-dimensional finite element model of a human lumbar annulus fibrosus. *Clin. Biomech.* 21, 337–344. doi:10.1016/j.clinbiomech.2005.12.001

Schroeder, Y., Huyghe, J. M., Donkelaar, C. C., and Ito, K. (2010). A biochemical/biophysical 3D FE intervertebral disc model. *Biomech. Model. Mechanobiol.* 9, 641–650. doi:10.1007/s10237-010-0203-0

Schroeder, Y., Wilson, W., Huyghe, J. M., and Baaijens, F. P. T. (2006). Osmoviscoelastic finite element model of the intervertebral disc. *Eur. Spine J.* 15, 361–371. doi:10.1007/s00586-006-0110-3

Shirazi-Adl, A. (1989). On the fibre composite material models of disc annulus-comparison of predicted stresses. *J. Biomech.* 22, 357–365. doi:10.1016/0021-9290(89)90050-X

Shirazi-Adl, A., Shrivastava, S. C., and Ahmed, A. M. (1984). Stress analysis of the lumbar disc-body unit in compression: a three-dimensional nonlinear finite element study. *Spine* 9, 120–134. doi:10.1097/00007632-198403000-00003

Simulia Inc. Rhode Island, USA, ABAQUS/StandardTM version 6.10. Dassault Systèmes Americas Corp., Waltham, MA.

Skaggs, D. L., Weidenbaum, M., Iatridis, J. C., Ratcliffe, A., and Mow, V. C. (1994). Regional variation in tensile properties and biochemical composition of the human lumbar anulus fibrosus. *Spine* 19, 1310–1319. doi:10.1097/00007632-199406000-00002

Spencer, A. J. M. (1984). “Constitutive theory for strongly anisotropic solids,” in *Continuum Theory of the Mechanics of Fibre-Reinforced Composites*, ed. A. J. M. Spencer (New York, NY: Springer-Verlag), 1–32.

Tong, P., and Fung, Y. C. (1976). The stress-strain relationship for the skin. *J. Biomech.* 9, 649–657. doi:10.1016/0021-9290(76)90107-X

Veres, S. P., Robertson, P. A., and Broom, N. D. (2010). The influence of torsion on disc herniation when combined with flexion. *Eur. Spine J.* 19, 1468–1478. doi:10.1007/s00586-010-1383-0

Vernon-Roberts, B., Moore, R. J., and Fraser, R. D. (2007). The natural history of age-related disc degeneration: the pathology and sequelae of tears. *Spine* 32, 2797–2804. doi:10.1097/BRS.0b013e31815b64d2

Wang, J. L., Parnianpour, M., Shirazi-Adl, A., and Engin, A. E. (2000). Viscoelastic finite-element analysis of a lumbar motion segment in combined compression and sagittal flexion. Effect of loading rate. *Spine* 25, 310–318. doi:10.1097/00007632-200002010-00009

Wang, J. L., Parnianpour, M., Shirazi-Adl, A., Engin, A. E., Li, S., and Patwardhan, A. (1997). Development and validation of a viscoelastic finite element model of an L2/L3 motion segment. *Theor. Appl. Fract. Mech.* 28, 81–93. doi:10.1016/S0167-8442(97)00032-3

Wilke, H. J., Neef, P., Caimi, M., Hoogland, T., and Claes, L. E. (1999). New in vivo measurements of pressures in the intervertebral disc in daily life. *Spine* 24, 755–762. doi:10.1097/00007632-199904150-00005

Keywords: functional spinal unit, finite element modeling of intervertebral disk, annulus material modeling, hyperelastic and anisotropic behavior, uniaxial vs. biaxial stress state

Citation: Momeni Shahraki N, Fatemi A, Goel VK and Agarwal A (2015) On the use of biaxial properties in modeling annulus as a Holzapfel–Gasser–Ogden material. *Front. Bioeng. Biotechnol.* 3:69. doi: 10.3389/fbioe.2015.00069

Received: 16 September 2014; Accepted: 30 April 2015;

Published: 03 June 2015

Edited by:

Fabio Galbusera, Ulm University, GermanyReviewed by:

Amir Abbas Zadpoor, Delft University of Technology, NetherlandsMaxim Bashkuev, Charité – Universitätsmedizin Berlin, Germany

Copyright: © 2015 Momeni Shahraki, Fatemi, Goel and Agarwal. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Ali Fatemi, Mechanical, Industrial and Manufacturing Engineering Department, University of Toledo, 2801 West Bancroft Street, Toledo, OH 43606, USA, afatemi@eng.utoledo.edu