Prediction of intervertebral disc mechanical response to axial load using isotropic and fiber reinforced FE models

Keywords:

• Intervertebral disc
• annulus fibrosus
• anisotropy
• collagen fiber
• finite element modeling

1. Introduction

Collagen fibers play a major role in the determination of the mechanical behavior of soft tissues. In the intervertebral disc (IVD), it was shown that they contribute in both tissue tensile and compressive stiffness (Römgens et al. 2012). IVD collagen network is progressively organized from the nucleus pulposus (NP) to the annulus fibrosus (AF). In the NP, the collagen density is low and the fibers are randomly disposed. However, they are oriented in two directions in the outer lamellae of the AF (about ±30° to the transverse plane). The organization of the collagen fiber network provides an anisotropic mechanical behavior to the AF.

Several finite element (FE) models taking into account the anisotropy of the AF have been proposed and validated for several cases of load such as compression, twist and bending of the whole IVD (Jacobs et al. 2014; Reutlinger et al. 2014) or compression, tension and shear of the AF sample (Mengoni et al. 2015; Hollingsworth and Wagner 2011).

In order to develop a new method of assessment of degeneration based on the coupling of Magnetic resonance imaging (MRI) and FE modeling, two models of IVD were developed. The first model is isotropic and the second one incorporates an anisotropic behavior of the AF. Through a comparison of these two models using an experimental compressive test, the purpose of this study is to determine the more suitable formulation to assess IVD mechanical properties using a new coupled (MRI-FE modeling) approach to investigate degeneration.

2. Methods

The two models were based on the biphasic swelling formulation (Mow et al. 1980). Two intrinsically incompressible phases were defined: a solid phase and a fluid phase to describe respectively the porous extracellular matrix (ECM) and the water filling the pores. In this formulation, the second Piola-Kirchoff (PK2) stress tensor $\mathit{S}$ is written:

$$\mathit{S}={\mathit{S}}^e-J\left(p+\mathrm{\Delta }\mathit{\pi }\right)·{\mathit{C}}^{-1}$$

where ${\mathit{S}}^e$ is the PK2 stress tensor related to the solid, J is the local volume ratio, p is the fluid pressure within the IVD, Δπ the osmotic one and $\mathit{C}$ is the right Cauchy-Green strain tensor.

The ECM was considered hyperelastic. The PK2 stress tensor related to the ECM is determined in term of W, the strain energy density:

$${\mathit{S}}^e=2·\left(\frac{\text{d}W}{\text{d}\mathit{C}}\right)$$

For the isotropic model, a nearly incompressible Mooney-Rivlin strain energy density was taken. For the fiber reinforced model, W was defined by a summation of a nearly incompressible model of Neo-Hooke (W^NH) and an anisotropic part:

$$W=W^{\text{NH}}\left(\overline{\mathit{C}},J\right)+\sum _i{W}_i^{\text{anis}}\left({\overrightarrow{e}}_i,\overline{\mathit{C}}\right)$$

where, $W_{\text{anis}}^i$ is the anisotropic strain energy density of a simple family of fibers, $\overline{C}=J^{-2/3}.\mathit{C}$ and ${\overrightarrow{e}}_i$ defines the direction of the fibers of the family i in the reference configuration. $W_{\text{anis}}^i$ is defined by a ρ-model as described in Holzapfel and Ogden 2010. This model allows to consider the anisotropy of the outer AF (OAF), the isotropy of the NP and the partial anisotropy of the intermediate zone. In the cartilaginous endplates (CEP), the tissue mechanical behavior was taken Neo-Hookean. Two families of fiber were defined. We supposed that these two families have the same mechanical properties and two different orientations ±α.

The experimental tests are detailed in Ghiss et al. 2015. Briefly, a compressive displacement was applied on porcine IVDs placed in a wet room of a compressive bench while recording the force. The first Piola-Khirchoff (PK1) stress is then deduced. The set-up was placed on a 4.7T magnet in order to acquire ρ_H-weighted MRI sequences during the tests.

The IVD geometry and the initial porosity spatial field (Figure 1) were determined using transverse (X,Y) MRI sequences (resolution: 312 × 312 μm, slice thickness: 0.75 mm). The three IVD zones (AF, NP and CEP) were delimited using the porosity field.

Figure 1. (A). IVD geometry. (B). Spatial initial porosity field in the central transverse plane.

The mechanical properties evolution from NP to OAF was progressive, continuous and linear.

Based on these previous experimental tests, an optimization process was performed in order to define the material parameters of each model.

3. Results and discussion

Comsol multiphysics was used for the resolution of the FE models. The optimized parameters were in line with the literature. For the fiber reinforced model, the fiber orientations were ±25.5° to the transverse plane.

The identification procedure was about 15 h shorter and took less memory space with the isotropic model since less parameters have to be identified compared to fiber reinforced one (on workstation with 3.6 Ghz quad core processor and 16Go of RAM).

The two models were able to reproduce the experimental PK1 stress evolution (Figure 2). The average relative error between experimental and numerical curves were 3.11% for the isotropic model and 3.61% for the fiber reinforced one.

Figure 2. Evolution of the experimental (black) and the numerically obtained (blue and red) PK1 average stress.

The measured IVD volume loss after the relaxation was 5% experimentally, 5.8% with the isotropic model and 5.2% with the fiber reinforced one.

The results show that even without taking into account the AF anisotropy, it is possible to predict with proper accuracy the IVD mechanical response to axial load.

For the isotropic model, the ECM was stiffer given that it includes the fiber rigidity but in an isotropic way.

The role of collagen fibers will be more considerable in the axial rotation or the non axial loads such as flexion or extension since the distribution of the strain in IVD is more complex and non homogenous with these loads.

The difference found in the volume loss shows that the fiber reinforced model better describes the water exchange between the IVD and the extra-discal area. This point will be enhanced for higher loads since the volume loss is related to the final porosity and the IVD bulge.

4. Conclusions

Based on MRI images and FE modeling, two IVD models were developed and compared using experimental results of a compressive test. The results show that the two models were able to predict the IVD mechanical response to axial load. For other types of load, the fiber reinforced model appears more relevant describing associated experimental results. Moreover, coupling MRI water content measurement and FE modeling seems to be an efficient tool to find a link between IVD mechanical properties and the porosity distribution. The latter link will be used in the assessment of mechanical degenerative changes.