Abstract

The study of the mechanical properties of the annulus fibrosus of the intervertebral discs is significant to the study on the diseases of lumbar intervertebral discs in terms of both theoretical modelling and clinical application value. The annulus fibrosus tissue of the human intervertebral disc (IVD) has a very distinctive structure and behaviour. It consists of a solid porous matrix, saturated with water, which mainly contains proteoglycan and collagen fibres network. In this work a mathematical model for a fibred reinforced material including the osmotic pressure contribution was developed. This behaviour was implemented in a finite element (FE) model and numerical characterization and validation, based on experimental results, were carried out for the normal annulus tissue. The characterization of the model for a degenerated annulus was performed, and this was capable of reproducing the increase of stiffness and the reduction of its nonlinear material response and of its hydrophilic nature. Finally, this model was used to reproduce the degeneration of the L4L5 disc in a complete finite element lumbar spine model proving that a single level degeneration modifies the motion patterns and the loading of the segments above and below the degenerated disc.

1. Introduction

The annulus fibrosus (AF) of the human intervertebral disc (IVD) presents a complex structure. Located on the radial periphery of the IVD, the AF is believed to experience a combination of compressive, tensile, and shear stresses during weight-bearing and intervertebral joint motions [1, 2]. For these reasons, its structure combines a highly orientated fibres network with a multiphasic behaviour  [35]. Furthermore, the presence of osmotic forces, due to the hydrophilic nature of the proteoglycans aggregates (PGs) [1, 2], ensures that the fibre network of the tissue works under tensile prestress and that the fluid inside the tissue is under pressure. This is essential for the correct functioning of the tissue since the fluid typically resists compression and fibres resist tension [6].

In the last years several analytical and numerical models have been developed to simulate the AF tissue behaviour. Some of them, like those proposed by Simon and coworkers [7], implemented its poroelastic behaviour and swelling. Others authors, like Lai et al.  [8] and Sun et al.  [9], or Frijns et al.  [10] and van Loon et al.  [11], developed multiphasic mixture models that permit a good modelling of poromechanics of the tissue. All these models, nevertheless, neglected the mechanical behaviour of the extracellular matrix related to collagen network. In this regard, others authors  [12, 13], who considered that the fibres network implies a linear viscoelastic response, did not consider the high nonlinearity behaviour of the fibres, though.

Furthermore, in spite of the fact that AF tissue is subjected to several changes during life  [1, 1416], there are not many models to simulate the effects of degeneration. Some of them [17] simply removed some elements from the annulus. Others, like Rohlmann et al. [18], Polikeit et al. [19], and Natarajan et al. [20], showed the influence of geometrical, mechanical, or poromechanical parameters on disc behaviour but did not show clearly how the mechanical properties change with degeneration. More recently, Schmidt et al. [21] developed a poroelastic finite element model of the lumbar spine to investigate spinal response during physiological activities. In any case, nobody of the aforementioned authors has directly validated the influence of these changes on AF tissue behaviour.

On the other hand, the degenerative disc disease seems to be related to the modification of the biomechanical functioning of the spine. Abnormal mechanical loads and/or motion patterns have been related to risk of injury to the spine [22, 23]. Despite numerous studies on the mechanics of the degenerated disc, there is limited data on how this condition affects the adjacent caudal and cephalic segments or the lumbar spine stability contributing to the progression of disc degeneration  [1, 4, 14, 15, 18, 20, 2428]. The results of these studies vary considerably. Similarly to Kirkaldy-Willis and Farfan [29], some authors report instability during the early stages of degeneration [30, 31] while others rather show the opposite  [32, 33]. These partial contradictive results are probably enforced by the relative small number of specimens. Recently, Kettler et al. [34] analyze the influence of the degree of disc degeneration on the flexibility of lumbar spine segments based on the data from a large in vitro database. There are several studies that numerically analyze the degeneration of the adjacent levels after lumbar fusion [35, 36] but not much has been done to computationally reproduce the influence of a single-level lumbar degenerative disc disease on the behaviour of the adjacent segments except from the work of Ruberte et al. [37].

In this work a constitutive model to describe the mechanical behaviour of the normal and degenerated AF tissue of human IVD has been developed. The complete mathematical formulation is presented with particular focusing on the PGs behaviour modelling. Then, complete characterization and validation of the osmotic and elastic response for a normal and degenerated AF have been carried out  [4, 3840]. Afterwards, a complete finite element model of the lumbar spine was analyzed inferring the mechanical consequences on the overall behaviour of the spine taking into account some degree of degeneration at a single-level.

2. Methods

The annulus fibrosus (AF) of human IVD presents a very particular structure. It consists of a solid porous matrix, saturated with water, which mainly contains proteoglycans and collagen fibres. To simulate its behaviour a 3D osmo-hyperelastic model reinforced with two families of fibres was constructed and implemented. Several validation tests were performed in order to validate the accuracy of this model and its ability of reproducing some of the consequences of degeneration.

2.1. Constitutive Modeling of the Annulus Fibrosus

The strain energy density function (1), initially presented by Eberlein et al. [17] to characterize fiber reinforced materials extensively used for biological tissues, was modified to introduce the contribution of the osmotic pressure. This term takes into account the effect of the electric charged proteoglycans and it can be coupled with biphasic formulation.

As widely known, the Helmholtz free energy function can be divided into different components corresponding to the ground substance, the fibres network, and the material compressibility (see Nomenclature for notation). In this case it can be written as Accordingly, the 2nd Piola-Kirchhoff stress tensor is where the contribution of each part can be expressed as It can be seen that in (5), the term relative to the material compressibility has been modified to incorporate the osmotic pressure contribution. Thus, in agreement with Wilson et al. [13, 40], the water chemical potential can be defined as while the osmotic pressure gradient is calculated by means of the external and internal osmotic pressure: Both external and internal osmotic pressures depend on the external and internal osmotic coefficients, respectively, and on the external salt concentration, while the internal osmotic pressure depends also on the fixed charged density that is related to the proteoglycans content, Furthermore in (9), the fixed charged density associated with the affinity to water is calculated as: Finally, it can be considered that the hydraulic permeability of the ground substance matrix also decreases as the tissue is compressed. It can be assumed to be dependent on porosity changes in the following way [41]; therefore, where    is a positive coefficient equal to 15.

2.2. Material and Testing Simulations for Healthy Annulus Fibrosus Tissue

To characterize and validate the model presented in the previous section, we reproduced in silico some in vitro experimental results available in the literature for healthy samples of annulus.

2.2.1. Osmotic Swelling Behaviour

To point out the principal effects of the osmotic swelling on the behaviour of soft hydrated tissues, the most common validation test to assess the contribution of electrochemical effects to the response of the tissue is a pure swelling test.

A 1D simplified numerical model was run in ABAQUS 6.11 (Figure 1(a)). The constitutive equations were implemented in a user defined UMAT subroutine. A specimen of tissue ( mm), equilibrated in a NaCl solution with a salt concentration of , was placed inside an impermeable confined compression chamber (Figure 1(b)). After an equilibrium phase, at the external salt concentration was changed to (Figure 1(c)) and, at , the external salt concentration was returned to . To simulate the swelling effect, only the axial displacements were allowed for all nodes while free flux condition ( bar) was assumed only at the lower nodes. See Table 1 for material properties.

2.2.2. Characterization of the Normal IVD Annulus Tissue

To characterize the solid matrix AF behaviour, the polynomial strain energy function shown before was used. The values of the elastic constants (, and ) were determined using the stress-strain response under a traction axial load for a specimen with two families of fibres, placed at [42]. The theoretical model behaviour was fitted to experimental data presented by Ebara et al. [38] using MatLab v.7.1. It can be observed (Figure 2) how the resultant curve can be considered a good average fit of experimental data of the whole annulus [17, 43, 44]. These constants are summarized in Table 2. Here, we have to point out that each part of the disc could have been considered with different mechanical properties but we were interested in the overall behaviour of the annulus.

To validate the global behaviour of the proposed model two different compression tests were performed in order to compare our results with previous data in the literature  [4, 39]. A 3D finite element (FE) mesh was created (Figure 3(a)). The chosen geometry was in agreement with that proposed by Wilson et al. [13]. Each family of fibres was placed radially with an angle orientation of with respect to the -plane [39, 42] (Figure 3(b)). In both tests the displacements for the lower nodes were restricted in the -direction while for lateral perimetric nodes the radial displacements were confined. In the first test, a zero pore pressure condition was imposed at the top while the remaining surfaces were assumed to be impermeable. In this case, the loading protocol implied that the mesh was axially compressed by 5% strain followed by a relaxation period (Figure 3(c)). The second test was conducted following the guidelines proposed by Iatridis et al. [4]. In this case, free flux condition ( bar) was assumed on the bottom surface. Initially a compression strain of 10% (velocity strain ratio 0.001 ) was applied on top followed by 2500 s of relaxation at fixed displacement. Then, this loading cycle was repeated three times adding a 5% of deformation strain at each increment (Figure 3(d)).

2.3. Material and Testing Simulations for Degenerated Annulus Fibrosus Tissue

Several experimental studies have shown how degeneration process produces serious changes on AF tissue behaviour  [14, 15, 24, 28]. Particularly, Iatridis et al.  [4] showed that the residual stress at the equilibrium point under a compression strain of for the degenerated (grade IV) AF tissue was nearly twice in comparison with the normal AF tissue. Therefore, the same test procedure made for normal tissue was done for the degenerated one.

Some authors  [4, 15, 27] agree with relating the reduction of proteoglycans content to disc degeneration. In our simulations, this phenomenon was taken into account reducing the tissue initial fixed charge density . Furthermore, AF tissue permeability has also shown a decrease with degeneration process [4] that also was considered into this model.

For sure, also the elastic behaviour of the AF would be subjected to changes [4, 14, 24]. In particular, an increase of the stiffness has been addressed [4, 24]. However, it has been seen that annulus fibrosus stiffening is not due to the reduction of water content which slightly changes with degeneration [15]. Furthermore, other experimental works [45] have shown that this stiffening is not related either to fibres behaviour. Therefore, one reasonable assumption would be to suppose that this increase could be related to the ground substance material behaviour. To simulate this phenomenon the value of the elastic constant was modified in order to reproduce a stress increase of 100% for a compression strain of at equilibrium [4] (see Table 2).

To validate these assumptions the same test procedure as for healthy AF was done for the degenerated one following Iatridis et al. [4] simulating a Grade IV degeneration.

2.4. FE Simulation of the Complete Lumbar Spine

A complete finite element of the lumbar spine (see [43] for model details) was used. This model is composed of the lumbosacral segment (L1-S1) including the five intervertebral discs and the most important ligaments. The validation of this model with experimental data of the literature is presented [46]. The degeneration of the AF of the L4L5 disc as in the previous section was simulated and the movement and tensional response of the overall spine were compared with those of the healthy one during flexion-extension loading (see Figure 4).

The mechanical properties of the different tissues involved are summarized in Tables 3 and 4. It can be seen that degeneration of L4L5 disc has been considered by modifying the permeability, the void ratio, and the stiffness of the annulus fibrosus as in previous section.

3. Results

3.1. Healthy Annulus Fibrosus

Regarding the free swelling test, an excellent agreement between these calculations and those presented by Wilson et al. [40] for human cartilage was achieved (Figure 5). It can be seen that our simulations are almost equal to the mechanoelectrochemical model of Wilson et al.  [40]. Only a slight difference can be appreciated in the swelling part of the test, since the shrinking part perfectly fits.

For the complete validation set described in Section 2.2.2, two different simulations were performed. In the first, a relaxation test, the stress-time behaviour, was compared with experimental results presented by Schroeder et al. [39]. The trends of the numerical and experimental curves showed a very good agreement (Figure 6(a)). It can be noted how the numerical curve is almost completely falling within the experimental interval found by Schroeder et al. [39]. The differences between numerical and experimental results could be related to the permeability value not measured in the experimental set-up and the neglecting of the intrinsic viscoelasticity of the matrix and collagen fibres in our numerical model.

The second test, which implies a cycled loading to the sample, was conducted following the guidelines proposed by Iatridis et al. [4]. In this test, the normalized stress-time behaviour of the material was recorded. The comparison between the experimental results, reported by Iatridis et al. [4], and our calculations are shown in Figure 7(b). In general a good agreement can be found. Furthermore, if we consider the ratio between the pick and the valley of each cycle, it can be noted that it is quite perfect for the second cycle (exp. 4.07, FE 4.11) while the maximal difference was found for the last cycle (exp. 2.55, FE 3.13).

3.2. Degenerated Annulus Fibrosus

Using the same constitutive equations for healthy AF but taking into account some modifications of different constants as mentioned in Section 2.3, the same experimental test performed by Iatridis et al. [4] was numerically reproduced. The equilibrium elastic stress-stretch responses of a normal and degenerated model are shown in Figure 8. The numerical results showed a good fit with experimental data, and these prove the matrix stiffening with degeneration.

Subsequently, the stress-strain curves in traction for normal and degenerated cases are plotted in Figure 9. The results showed an increase of the elastic modulus in the first part of the curve but no important changes are observed in the second part . These results were in agreement with the experimental ones found by Guerin and Elliott [24], who obtained an increase of and , respectively.

3.3. Finite Element Simulation of Healthy and Degenerated Lumbar Spine

First of all the finite element model of the healthy lumbar spine [43] was validated using data from the literature. The degree of rotation of the different segments was compared with Guan et al.  [46]. It can be seen in Figure 10 how this model mimics the behaviour of the healthy spine both under flexion and extension moments. The relative rotation between each pair of vertebrae is plotted and compared with the experimental results obtained from the literature. It can be seen that the response of the numerical simulation fits within the dispersion limits of the experimental protocol.

Now the influence of the degeneration of the disc at the L4L5 level is analyzed. First of all, the movement of the spine when the D45 disc is degenerated is compared with the healthy spine shown before. The comparison between both models is presented in Figure 11. It can be seen that the effect of degeneration is not only located at the damaged level L45 since the adjacent levels are also modified with respect to the healthy scenario. It can be appreciated that with degeneration the pair L4L5 gets stiffer and the degree of rotation is smaller in both flexion and extension. This effect also extends to the adjacent levels (L3L4 and L5S1) but the loss of rotation is less pronounced.

Moreover, the stresses in the discs can also be analyzed to prove the influence of the degeneration of one disc on the overall behaviour of the spine. In Figure 12 the maximum and minimum principal stresses are plotted at the levels affected by the degeneration of D45. The remaining discs, D12 and D23, are not shown because these are not affected by the degenerated D45. The stresses are shown for the posterior and anterior part of each disc taking into account both the flexion and extension moment. It can be seen that the stresses are not very different in any case, but it is interesting to highlight that the stresses are more similar in D45 disc between healthy and degenerated situation than in D34 and D51 which are only influenced by degeneration in L4L5 level.

4. Discussion

The goal of this work was to construct a suitable constitutive model to characterize the mechanical behaviour of normal and degenerated annulus fibrosus of the IVD. This model takes into account the biphasic nature and also the clear preferential orientation of the collagen fibres in the annulus. The phenomenon of swelling, indeed, is really important for a correct simulation of the stress-strain behaviour of cartilage as it has evidently shown in previous numerical [7, 40] and experimental works [39, 47, 51]. Furthermore, it plays an important role in the stress relaxation, in fluid regulation and therefore in tissue permeability. Here, the model for fibred reinforced materials, originally proposed by Holzapfel [52], was modified to incorporate the swelling contribution due to the PGs being coupled with the biphasic formulation.

Due to the hydrophilic behaviour of cartilage, an increase of PG’s concentration (e.g., due to a compression of the tissue) produces an influx of fluid into the cartilage. This implies an increase of the internal pressure giving the tissue more capacity to absorb loads. To simulate this effect, the component of the 2nd Piola-Kirchhoff stress tensor related to the material compressibility was modified introducing the osmotic gradient pressure. This term considers the tissue fixed charge density as a function of the volume variation [26, 40].

Finally, a strain dependent permeability function was implemented in the model. A power-low function of the porosity rate change was chosen relating the actual and the initial permeability [53]. The equivalence between this function and an exponential one [54] was demonstrated by Riches et al. [41].

To validate the behaviour of the poromechanical part of the model a 1D free swelling test was reproduced in silico. Because of the lack of data regarding the AF tissue free swelling behaviour, human cartilage tissue behaviour was simulated and compared with previous results found in the literature. Since this tissue did not present a fibred structure, no fibres were considered in the model in agreement with Wilson et al. [40]. The response of the model showed an excellent agreement with these authors  [40].

Once the model was validated for free swelling, different tests were performed to guarantee the accuracy of the results. The material constants of the annulus fibrosus related to the fluid and ionic contribution () were obtained from the literature [4, 26, 40, 47]. On the other hand, we were not interested in detecting differences across the zones of the annulus (anterior, posterior, inner, or outer). Therefore, only one set of elastic constant values (, and ) was fitted to the experimental data by Ebara et al. [38] representing a good average behaviour for the whole annulus. This assumption was fully demonstrated in a previous work [43].

First, the experimental confined compression test proposed by Schroeder et al. [39] for the healthy AF tissue was reproduced in silico. The FE model showed a good agreement with the experimental results. The numerical curve, in fact, fills within the experimental curves range. Furthermore this first analysis showed that our model was capable of predicting values of stresses equivalent to those measured experimentally. However, the found value of stress at the end of the relaxation period was rather different in comparison with the average value reported in the literature. Many explanations could be made for this behaviour. For example the PGs content, or equivalently the fixed charge density, strongly influences the material response [9, 25, 55]. Also the distribution and orientation of the fibres direction could influence the behaviour of the tissue  [39, 52].

Then, a transient response behaviour analysis as proposed by Iatridis and coworkers [4] was carried out. Since the experimental curve presented by Iatridis et al. [4] is only related to a specific specimen and not to an average curve; the value of force was normalised and compared. This procedure, furthermore, is in agreement with that proposed by DiSilvestro and Suh [56] and Wilson et al. [13]. From the other hand, the results obtained with the experimental result proposed by Schroeder and coworkers [39] are adequate to guarantee the goodness of found internal stress values. The results showed a good accordance with experimental data. In particular the model showed the ability to recover the value of the stress with the same velocity of experimental data.

In Section 4 a criterion to define a degenerated AF tissue behaviour was defined. In particular the experimental evidences of Iatridis et al. [4] were used to modify the constant values of the model. In Figure 8 the normalised value of the stress was compared between experimental and numerical results for normal and degenerated AF tissue. It can be noted how the proposed procedure was able to reproduce the effects of degeneration: increase of stiffness and reduction of nonlinearity of the behaviour. Iatridis et al. [4] mentioned that the increase in the elastic modulus with degeneration is likely related to an increase in tissue density resulting from the loss of water content. However, the reduction in nonlinearity behaviour with degeneration could suggest a diminished compaction effect of the degenerate tissues at large deformations which could be related to structural changes in matrix. Besides, the effect of the increase of the ground substance elastic modulus was analysed under a pure axial tensile load. The obtained results are also in agreement with the experimental behaviour showed by Guerrin et al. [24]. Experimental results, in fact, showed an increase of the elastic modulus in the toe part, before fibres act, but when the fibres act there is no additional stiffening.

Finally, the constitutive model that has been presented here has been applied to a complete finite element model of the lumbar spine and degeneration of one disc has been provoked by modifying its mechanical properties. In this analysis it has been obtained that the degeneration of one level affects the biomechanics of the adjacent levels; in particular a decrease of the range of motion in flexion/extension has been seen. These results are in agreement with Kettler et al. [34] who analyzed 203 lumbar spine segments obtaining that the range of motion decreased with degeneration. This same result was radiographically shown by Mimura et al. [33] for flexion/extension moments. Regarding the stress distributions in the discs, there is no much data to compare. As mentioned in the introduction, only the work of Ruberte et al. [37] analyzed the influence of single-level lumbar degenerative disc disease on the behaviour of the adjacent levels using the finite element model. However, the conclusions of their study are opposite to what we have obtained here and other experimental results [30, 31]. They obtained that as degeneration progresses the stiffness increased but was significantly less than the healthy model; however it seems that the clinical evidence proves our results. With respect to the stresses distribution, our results showed that there is a slight modification on the adjacent levels when one single level is degenerated. However, this data cannot be transferred to the clinical evidence because it is not quantitatively known to what extent loading changes involve degenerative changes.

5. Assumptions and Limitations

The numerical model for the degeneration of the annulus fibrosus has several limitations. First, only the swelling effect has been introduced to describe its behaviour. There are more effects related to the electrochemical nature of cartilage that will be taken into account in further developments. As shown in Figure 6, the viscoelastic nature of the solid matrix can affect the overall response of the tissue and it should be included to improve this model. Numerical predictions have been compared with experimental results; however, there is a lack of data in the literature about the poromechanics of the healthy or degenerated AF and therefore more experimental data would be required to better guarantee our results. On the other hand, the stiffening of the annulus fibrosus has been reproduced by modifying the mechanical behaviour of the matrix. It is known that with degeneration, collagen fibers are reoriented and broken. This contribution should be introduced in future works.

With respect to the lumbar spine finite element model there are also limitations and underlying assumptions. The most important one is that the same finite element model has been used for healthy and damaged spine. It is known that the height of the discs is modified with degeneration and also the internal structure of the disc but here only the mechanical properties of the degenerated disc have been modified. Notwithstanding the fact that the height of the discs can modify the biomechanical response of the spine, it has been seen  [57] that when the disc height is reduced, flexibility of the motion segment decreased. Therefore, this effect would contribute in the same direction to the results that have been obtained here.

6. Conclusions

In conclusion, the developed model is capable of simulating the poromechanical behaviour of normal and degenerated AF tissue with a good approximation. The most important feature of this model is related to the coupling between a strongly fibred matrix and its intrinsic affinity to water regulated by its ionic nature. On the other hand, the degeneration of cartilage is a very complex process that has been here simplified for AF tissue with accurate predictions.

Moreover, the present study provides a qualitative analysis of the influence of single level disc degeneration on the mechanics of the adjacent segments under flexion/extension moments. It has been seen that degeneration modified the degree of motion and loading of both the degenerated disc and adjacent levels. These changes could increase the risk of progression of degeneration to the nearest segments of the spine.

Nomenclature

:Helmholtz free energy function
:Ground substance Helmholtz free energy function
:Fibres behaviour Helmholtz free energy function
:Material compressibility Helmholtz free energy function
:Right Cauchy-Green tensor
:Modified right Cauchy-Green tensor
:Relative volume change
:Structural tensor of the fibre direction
:Fiber direction
:2nd Piola-Kirchhoff stress tensor
:Ground substance 2nd Piola-Kirchhoff stress tensor
:Fibres behaviour 2nd Piola-Kirchhoff stress tensor
:Material compressibility 2nd Piola-Kirchhoff stress tensor
:Project tensor
:Fourth-order unite tensor
:Ground substance material constants
:Collagen fibres behaviour constants
:Tissue incompressibility modulus
:First invariant
:Fibres invariants
:Water chemical potential
:Fluid pressure
:Osmotic pressure gradient
:External osmotic pressure
:Internal osmotic pressure
:Universal gas constant
:Absolute temperature
:Internal osmotic coefficient
:External osmotic coefficient
:External salt concentration
:Fixed charge density
:Initial fixed charge density
:Initial fluid fraction
:Permeability
:Initial permeability
:Porosity
:Initial porosity
:Positive coefficient.

Conflict of Interests

There is not any financial and personal relationship with other people or organisations that could inappropriately influence this work.

Acknowledgment

The support of the Spanish Ministry of Economy and Competitiveness through DPI 2011-23148 project is highly appreciated.