G Dusfour
email: gilles.dusfour@gmail.com

S Lefloc'h A

P Cañadas

D Ambard

Heterogeneous mechanical hyperelastic behavior in the porcine annulus fibrosus explained by fiber orientation : an experimental and numerical approach

Keywords: annulus fibrosus, mechanical behavior, intervertebral disc

Our aim is to estimate regional mechanical properties of the annulus fibrosus (AF) using a multi-relaxation tensile test and to examine the relevance of using the transverse dilatations in the identification procedure. We collected twenty traction specimens from both outer (n = 10) and inner (n = 10) sites of the anterior quadrant of the annulus fibrosus of one pig spine. A one-hour multi-relaxation tensile test in the circumferential direction allowed us to measure the force in the direction of traction and the dilatations in all three directions. We performed a specific-sample finite element inverse analysis to identify variations, along the radial position, of material and structural parameters of a hyperelastic compressible and anisotropic constitutive law.

Our experimental results reveal that the outer sites are subjected to a significantly greater stress than the inner sites and that both sites exhibit an auxetic behavior. Our numerical results suggest that the inhomogeneous behavior arises from significant variations of the fiber angle taken into account within the hyperelastic constitutive law. In addition, we found that the use

Theoretical and numerical models of the mechanical behavior of the annulus fibrosus (AF) are of primary importance to understanding the development of the pathologies, to propose early diagnostic tools, and to evaluate possible treatment strategies for the intervertebral disc (IVD). Identifying a complete set of parameters from a structure-based finite constitutive model of the AF provides direct input for meaningful IVD models. These models are in turn a crucial prerequisite for the understanding of nonpathological and pathological IVD and provide a resource for the development of innovative prostheses [START_REF] Noailly | Finite element study of a novel intervertebral disc substitute[END_REF]). Moreover, mechanical parameters can also be useful in the field of tissue engineering and to understand pathological development through the use of growth and remodeling models [START_REF] Oomen | Growth and remodeling play opposing roles during postnatal human heart valve development[END_REF]; [START_REF] Zöllner | Growth on demand : reviewing the mechanobiology of stretched skin[END_REF]) on the IVD.

The intervertebral disc is a fibrocartilaginous joint that gives mobility, damping, and stability to the whole vertebral column. The central part is composed of the soft and gelatinous nucleus pulposus (NP), which is subject to large compressive loads. The NP transmits radial forces to the surrounding AF which transforms them into circumferential tension [START_REF] Newell | Biomechanics of the human intervertebral disc : A review of testing techniques and results[END_REF]).

The AF is a complex composite material made up of incomplete concentric lamellae which exhibit aligned fiber bundles with a crisscross angle that alternates above and below the transverse plane [START_REF] Peacock | Observations on the postnatal structure of the intervertebral disc in man[END_REF]; [START_REF] Marchand | Investigation of the laminate structure of lumbar disc anulus fibrosus[END_REF]). Furthermore, collagen type I and collagen type II fibers are embedded in the ground matrix mainly composed of proteoglycans and are respectively preponderant in the outer and inner AF [START_REF] Eyre | Types i and ii collagens in intervertebral disc. interchanging radial distributions in annulus fibrosus[END_REF].

These large macromolecules confer swelling properties from the material to the organ scales and play an important role in the radial mechanical behavior of the AF. Spacial variation of the collagen and proteoglycan contents along the radial direction have been reported by [START_REF] Adams | Biochemical aspects of development and ageing of human lumbar intervertebral discs[END_REF] and [START_REF] Iatridis | Measurements of proteoglycan and water content distribution in human lumbar intervertebral discs[END_REF] : the collagen concentration increases from inner to outer sites, whereas the proteoglycan concentration decreases.

Spatial variations of the collagen fiber angle from the outer to the inner periphery as well as from the anterior to the posterior quadrant of the AF have been densely reported by [START_REF] Holzapfel | Single lamellar mechanics of the human lumbar anulus fibrosus[END_REF], [START_REF] Cassidy | Hierarchical structure of the intervertebral disc[END_REF],

and [START_REF] Guerin | Degeneration affects the fiber reorientation of human annulus fibrosus under tensile load[END_REF] : the collagen bundles are oriented at about ± 30°in the outer periphery to about ± 45°in the inner periphery to the transverse plane. This structural heterogeneity generates a strong anisotropic and an inhomogeneous stress-strain behavior, which has been well documented by [START_REF] Skaggs | Regional variation in tensile properties and biochemical composition of the human lumbar anulus fibrosus[END_REF] and [START_REF] Elliott | Anisotropic and inhomogeneous tensile behavior of the human anulus fibrosus : experimental measurement and material model predictions[END_REF]. However, data are lacking for the corresponding transverse radial and vertical strains. Several studies used experimental data to identify material and structural properties from constitutive models. [START_REF] Skaggs | Regional variation in tensile properties and biochemical composition of the human lumbar anulus fibrosus[END_REF] performed traction tests of single lamallae in the fiber direction and used a linear model to show that outer specimens are stiffer than inner specimens. [START_REF] Elliott | Anisotropic and inhomogeneous tensile behavior of the human anulus fibrosus : experimental measurement and material model predictions[END_REF] achieved tensile tests of multiple lamallae samples in the vertical, radial, and circumferential direction of the AF and used a linear anisotropic model to find a significantly greater stiffness in the circumferential direction than in the vertical or radial directions. However, the application and results of the previous studies are limited by the assumption of linearity. [START_REF] Wu | Mechanical behavior of the human annulus fibrosus[END_REF] were the first to propose an investigation of the mechanical behavior of AF samples with a theoretical analysis using a finite deformation theory taking into account two families of fibers. These latter authors made the assumption of incompressibility and found a close correlation between the experimental and theoretical results. However, the experimental conditions did not reproduce the in vivo temperature and humidity of the AF. On the other hand, recent studies mimic in vivo conditions [START_REF] O'connell | Theoretical and uniaxial experimental evaluation of human annulus fibrosus degeneration[END_REF]; [START_REF] Cortes | Extra-fibrillar matrix mechanics of annulus fibrosus in tension and compression[END_REF]). These latter authors combined a hyperelastic isotropic and compressible material with a fiber strain energy function to describe the impact of AF degeneration on its mechanical properties. Furthermore, [START_REF] Baldit | Experimental analysis of the transverse mechanical behaviour of annulus fibrosus tissue[END_REF] are the only authors to have investigated the transverse radial and vertical strain behavior of AF specimens in order to estimate the permeability from relaxation curves. Also, [START_REF] Acaroglu | Degeneration and aging affect the tensile behavior of human lumbar anulus fibrosus[END_REF] are among the few authors to combine quasi-static conditions and to present original results on the transverse behavior of AF samples. Unfortunately, no inverse analysis was carried out to identify mechanical parameters. These deficits in simultaneously acquiring the force in the direction of traction and the strains in all three directions combined to a hyperelastic modeling suggest the need for further studies to determine a complete set of material and structural parameters for the AF in quasi-static conditions. The completeness of the characterization method proposed here should enable us to identify variations of mechanical and structural parameters along the radial position of the AF.

The aims of this study were to simultaneously observe the force response and the deformations in all three directions of AF samples subjected to an original tensile test in order to better characterize the AF mechanical behavior. To this aim, we carried out mechanical parameter identification in the context of the global mechanical behavior. We hypothesized that the AF displays a compressible behavior and that it would be appropriate to model it in this way. Moreover, it is assumed that the gradient in fiber angle from the outer to the inner AF is responsible for the differences in mechanical behavior.

Materials and methods

Specimen preparation

A total of twenty annulus sheet specimens were harvested from the inner (n = 10) and outer (n = 10) periphery of the annulus fibrosus of the spine of one large white pig of three month old (≈ 30 kg). Specimens were taken from the anterior region of the spine. Cervical (n = 4), thoracic (n = 4), and lumbar (n = 2) IVD were extracted. The spine was obtained from the abattoir of the University of Nimes. Each functional spinal unit was extracted and immediately immerged in a 50x50x30 mm 3 bath containing a 0.15 M NaCl solution to be frozen at -20°C (maximal frozen time : 4 months). Prior to extracting the disc, the spinous and transverse processes were cut using a saw. The discs were separated from the adjacent vertebral bodies with a custom double blade cutter to obtain a planar one millimeter thick annulus fibrosus frozen sheet. Subsequently, all specimens were cut from the outer and inner periphery of the anterior quadrant with a surgical blade, along the circumferential direction (Fig 1). processes are cut with a saw (1), sheet of annulus fibrosus is obtained with a homemade double blade cutter (2), outer and inner specimens are extracted with a surgical knife (3).

Mechanical testing

Prior to testing, each specimen was equilibrated in a 0.

Image analysis

Image analysis was performed to estimate experimental displacement and to compute the deformation gradient tensor

F = I + ∂ U ∂ X F = I + ∂ U ∂ X F = I + ∂ U ∂ X
, where I I I is the identity tensor, X X X the initial coordinates, and U U U the displacement from the material to the spatial configuration. The right Cauchy-Green tensor image contrast, and therefore we deposited a graphite powder on the upper surface as proposed by [START_REF] Boyce | Fullfield deformation of bovine cornea under constrained inflation conditions[END_REF]). Digital image correlation software was used with a measurement grid of about 1000 evenly spaced points and using template and analysis subset sizes of 81x81 pixels used to cross-correlate two numerical arrays and the open cv2 library used for image segmentation.

C = F T F C = F T F C = F T F

Mechanical modeling

Constitutive behavior

The anisotropic nonlinear behavior of the annulus fibrosus was modeled by a hyperelastic constitutive model split into a purely isotropic part, which describes the ground matrix and an additional anisotropic part, which describes the contribution of the collagen fibers. The constitutive behavior of the isotropic matrix part was described by a compressible Neo-Hookean energy function :

Ψ m = µ 2 (I 1 -3 -2ln(J)) + λ ln(J) 2 2 , (1)
where

I 1 = trC C C is the first invariant of the right Cauchy-Green tensor C C C, J
=dv/dV is the volume change, µ the shear modulus, and λ the bulk modulus.

The contribution of the collagen fibers was modeled with the strain energy density proposed by [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF]; [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF]) :

Ψ f = i=4-6 k 1 2k 2 (e k 2 [κI 1 +(1-3κ)Iα i -1] 2 -1), (2)
where k 1 > 0 is a stress-like parameter and k 2 > 0 is a dimensionless parameter. Moreover α i α i α i is the unit vector representing the direction of the family of fibers i, and κ represents the fiber dispersion. An isotropic fiber distribution corresponds to κ = 1/3 while the ideal aligned fiber distribution corresponds to κ = 0. Here,

I α i α i α i = C C C : ( α i α i α i ⊗ α i α i α i
) is the squared elastic fiber stretch and is taken into account in the constitutive model only when

I α i α i α i > 1.
Finally, the total strain energy stored during elastic dilatation in the annulus fibrosus is made up of two terms : the extrafibrillar matrix which is supposed to be isotropic and compressible and the contribution of the anisotropic collagen fiber network. This type of energy density for soft tissues has already been proposed by [START_REF] Menzel | Anisotropic remodelling of biological tissues[END_REF]

ψ = ψ m (J, I 1 ) + ψ f (I α 4 , I α 6 ).
(3)

Finite element model

The mechanical test we designed was clearly not a simple tensile test in that we applied the displacement only on the lower surface of the sample, inducing a flexion together with the traction of the sample. In order to take advantage of such a rich mechanical test, we measured the dilatations in all three directions and developed a sample-specific FE-based identification procedure. A numerical traction specimen was created using the 3D finite element mesh generator GMSH ( and dilatations in all three directions were measured on the same physical surfaces and in the same region of interest as for the experiment.

Material parameter identification procedure

We used three experimental datasets to perform the material and struc-

f = n i=1 [(F exp θ -F num θ ) 2 i + (λ exp R -λ num R ) 2 i + (λ exp Z -λ num Z ) 2 i ]. (4)
To investigate the influence of the account of the dilatations in the objective function 4, a second objective function was used to obtain a second set of material properties :

f = n i=1 [(F exp θ -F num θ ) 2 i , ( 5
)
where n is the number of data points for each experiment. The optimi-zation was terminated when i) a change of the mechanical property values had an effect smaller than 10 -6 on the objective function, ii) a change of the mechanical property values had an effect smaller than 10 -5 on the gradient of the objective function, or iii) the differences between mechanical property values through consecutive optimization iterations were less than 10 -5 . During optimization, the permissible range of values for the bulk modulus λ, the shear modulus µ, the fiber dispersion κ, the fiber angle α 1 , the initial fiber rigidity k 1 , and the rigidification parameter k 2 were restricted to 0.0 -0.5 MPa, 0.0 -0.8 MPa, 0.0 -1./3., 20°-50°, 0.0 -0.5 MPa, and 0 -5000 respectively, thus ensuring that the resulting values would be physiologically acceptable. For each parameter, the permissible range of values was chosen according to the work of [START_REF] Eberlein | An anisotropic model for annulus tissue and enhanced finite element analyses of intact lumbar disc bodies[END_REF], [START_REF] Malandrino | Regional annulus fibre orientations used as a tool for the calibration of lumbar intervertebral disc finite element models[END_REF], and [START_REF] Cassidy | Hierarchical structure of the intervertebral disc[END_REF]. The material and structural parameters were assumed to be constant throughout the entire specimen (λ, µ, κ, k 1 , k 2 and α 6 α 6 α 6 = -

α 4 α 4 α 4 ).
Finally the normalized root mean square error (NRMSE) was used to assess the fitting quality between the predicted and measured forces and dilatations :

N RM SE = n i=1 (σ num i -σ exp i ) 2 n σ exp max -σ exp min × 100% (6) N RM SE = n i=1 (λ num i -λ exp i ) 2 n λ exp max -λ exp min × 100%. (7)

Statistical analysis

The numerical data are presented as median ± interquartile range. All statistical analyses were performed at 95% confidence level using R statistical

Results

Experimental results

All the measured data for each specimen are presented in figures 5, 6, and 7, as well as the median exterior and interior experimental responses. The AF quasi-static traction response exhibited nonlinearity on both stress-strain (fig 5) and strain-strain curves (Figs. 6 and7 ). The exterior specimens show a significantly greater Cauchy stress for λ θθ > 1.06 compared to the inner specimens (p < 0.05, fig 5). Furthermore, the magnitude of circumferential stress σ θθ was found to be 4.4 times larger in the outer periphery than in the inner periphery for λ θθ = 1.08. More remarkably, we noticed that most of the samples undergo a transition from an initial swelling to shrinkage for both radial and vertical dilatations (Figs. 6 and7). After an initial swelling up to 1.1 in the e ZZ e ZZ e ZZ direction and 1.06 in the e RR e RR e RR direction, most samples showed a remarkable contraction in the e ZZ e ZZ e ZZ direction up to λ ZZ = 0.75 and a minor contraction in the e RR e RR e RR direction up to 0.9. The exterior specimens exhibited a greater shrinkage, by up to a factor of 1.08, than the interior specimens in the e ZZ e ZZ e ZZ direction (no significant differences after application of the Bonferroni correction ; Fig 6). Similarly, no significant differences were observed for the dilatation in the e RR e RR e RR direction (Fig. 7). It should be noted that a few samples experience only swelling or contraction (Figs. 6, 7 and 8). Such conflicting behavior is mainly visible in interior samples. The measurement of dilatations in all three directions allowed us to estimate the variation of volume as a function of dilatation in the e θθ e θθ e θθ direction (Fig. 8). The volume of both exterior and interior specimens increased in median by 3.3%±2.9 and 4.6%±5.1, respectively, for low strain values and decreased in median by 6.5% in the inner periphery and by 13.8% ±9.2 in the outer periphery for large strain values.

No significant differences were observed for the volume variation between the inner and outer specimens.

Material parameter identification

The estimated average parameters from the exterior and interior part of the annulus fibrosus are summarized in Table 1 for both the force minimization (function 5) and the force and dilatation minimization (function 4). All parameters show no significant differences between the exterior and interior specimens (p > 0.05) except the fiber angle which is significantly lower in the exterior part (28.04°± 7.35 °) than in the interior part (34.11°± 12.05°, p < 0.05) for the force and dilatation minimization (function 4), corresponding to an increase by a factor of 1.22. As for the minimization on both force and dilatation, the minimization on force (function 5) shows a significant In addition, Table 2 gives a proper response on the stress-strain behavior and better results than the force minimization on the strain-strain behavior.

Discussion

Our study aims to observe regional mechanical variations from the outer to the inner sites of the anterior annulus fibrosus and to explain these variations through mechanical and structural parameter identification. Taking inspiration from [START_REF] Elliott | Anisotropic and inhomogeneous tensile behavior of the human anulus fibrosus : experimental measurement and material model predictions[END_REF], who used a linear anisotropic materiel model, we perform, to our knowledge, the first specific sample inverse finite element analysis on the annulus fibrosus to identify both mechanical and structural parameters of a hyperelastic, compressible, and anisotropic constitutive model on a large number of porcine IVD specimens (n = 10).

In order to characterize the true quasi-static behavior of annulus fibrosus, and to avoid temporal effects, we used a multi-relaxation test, as proposed authors for IVD [START_REF] Elliott | Anisotropic and inhomogeneous tensile behavior of the human anulus fibrosus : experimental measurement and material model predictions[END_REF]) and other tissues [START_REF] Yoshida | Material properties of mouse cervical tissue in normal gestation[END_REF]). Such a multi-relaxation test allows us to examine the two mechanical responses of soft tissues : the transient response, which deals with viscoelastic and/or poroelastic behavior, and the quasi-static response (independent of time), which only deals with the anisotropic hyperleastic behavior of the sample. Our experimental results, in accordance with [START_REF] Elliott | Anisotropic and inhomogeneous tensile behavior of the human anulus fibrosus : experimental measurement and material model predictions[END_REF], show that a tensile test of up to 8% of strain should last at least 1hr to obtain a time-independent mechanical response. This corresponds to a mean strain rate of 2.×10 -5 s -1 , much lower than commonly used strain rates [START_REF] Newell | Biomechanics of the human intervertebral disc : A review of testing techniques and results[END_REF] We observed the usual strain-stiffening of the samples in the circumferential direction e θθ e θθ e θθ , as already reported in a series of studies on porcine IVD [START_REF] Gregory | A comparison of uniaxial and biaxial mechanical properties of the annulus fibrosus : a porcine model[END_REF]; [START_REF] Monaco | A comparison between porcine, ovine, and bovine intervertebral disc anatomy and single lamella annulus fibrosus tensile properties[END_REF]) and human IVD [START_REF] Newell | Biomechanics of the human intervertebral disc : A review of testing techniques and results[END_REF]). In addition, our experimental results show that pig annulus fibrosus exterior specimens are 4.4 times stiffer than interior specimens in the circumferential direction, as has been observed previously by [START_REF] Ebara | Tensile properties of nondegenerate human lumbar anulus fibrosus[END_REF], Elliott &[START_REF] Elliott | Anisotropic and inhomogeneous tensile behavior of the human anulus fibrosus : experimental measurement and material model predictions[END_REF][START_REF] Skaggs | Regional variation in tensile properties and biochemical composition of the human lumbar anulus fibrosus[END_REF] on human samples.

Moreover, the sample-specific inverse finite element analysis shows that none of the mechanical parameters are significantly different between the exterior and interior specimens (p > 0.05). However, the structural parameter of fiber angle is significantly lower in the exterior specimens than in the interior specimens (p < 0.05), revealing the existence of a fiber angle gradient from the outer to the inner periphery such as that already directly observed in many vertebral species [START_REF] Elliott | Anisotropic and inhomogeneous tensile behavior of the human anulus fibrosus : experimental measurement and material model predictions[END_REF]; [START_REF] Holzapfel | Single lamellar mechanics of the human lumbar anulus fibrosus[END_REF]; Matcher et al. ( 2004)) but only supposed for pig through the study of [START_REF] Hsu | Diffusion tensor microscopy of the intervertebral disc anulus fibrosus[END_REF]. Furthermore, our results in terms of fibre angle range are in agreement with the measures made by [START_REF] Holzapfel | Single lamellar mechanics of the human lumbar anulus fibrosus[END_REF] and [START_REF] Guerin | Degeneration affects the fiber reorientation of human annulus fibrosus under tensile load[END_REF] on human samples, by [START_REF] Vergari | Lamellar and fibre bundle mechanics of the annulus fibrosus in bovine intervertebral disc[END_REF] on bovine samples, and by [START_REF] Hsu | Diffusion tensor microscopy of the intervertebral disc anulus fibrosus[END_REF] on porcine samples. Referring to the strain energy function used in our analysis, the rigidity gradient between the outer and inner periphery of the annulus fibrosus seems to be mainly due to the fiber angle gradient. However, the fibers of the AF do not only follow a gradient of angulations but also a concentration gradient of type I and II collagen from the outer to the inner part [START_REF] Eyre | Types i and ii collagens in intervertebral disc. interchanging radial distributions in annulus fibrosus[END_REF]. Such a concentration gradient is also visible for proteoglycans along the radial direction [START_REF] Iatridis | Measurements of proteoglycan and water content distribution in human lumbar intervertebral discs[END_REF]. Moreover, the number of complete and incomplete lamellae has not been taken into account in the inverse analysis method whereas it could greatly influence the global mechanical behavior. One possible further development of this study would be to experimentally determine complete and incomplete lamallae as well as inter-lamallae positions and to recreate them on the numerical sample to better predict the local mechanical behavior as experimentally identified by [START_REF] Vergari | Lamellar and fibre bundle mechanics of the annulus fibrosus in bovine intervertebral disc[END_REF]. Finally, while [START_REF] Guerin | Degeneration affects the fiber reorientation of human annulus fibrosus under tensile load[END_REF] showed that AF fiber reorientation can be calculated using an affine prediction for pure tensile loading -i.e. the fiber re-orientation behavior is exclusively driven by the macroscopic strain state -it is not certain that this prediction is true in our case, as already observed for different types of loading by [START_REF] Krasny | Kinematics of collagen fibers in carotid arteries under tension-inflation loading[END_REF] on arteries. It would be interesting to perform an identical experiment under a multiphoton microscope to identify whether fiber reorientation is affine or not. Moreover, models which can deal with a nonaffine assumption in hyperelasticity, such as models using homogenization techniques [START_REF] Michel | Microscopic and macroscopic instabilities in finitely strained fiberreinforced elastomers[END_REF][START_REF] Greco | A study of stability and bifurcation in micro-cracked periodic elastic composites including self-contact[END_REF][START_REF] El Nady | Computation of the homogenized nonlinear elastic response of 2d and 3d auxetic structures based on micropolar continuum models[END_REF], or more generally multi-scale modeling [START_REF] Morin | Non-affine fiber kinematics in arterial mechanics : a continuum micromechanical investigation[END_REF], could overcome the limitations of our study.

The use of transverse dilatations in the objective function revealed a si-gnificant improvement in the ability of our model to predict the transverse swelling behavior of the samples. The differences between the mechanical parameters identified with and without taking into account the transverse dilatations were significant for three parameters (k 1 , k 2 and κ). This highlights the importance of capturing the experimental transverse behaviors in addition to the circumferential one. Although some significant differences were noted between these parameters, taking into account these changes is unlikely to affect our conclusions, namely the spatial variation of the mechanical behavior could be explained by the significant variation of identified angles from inner parts (34.11°± 12.05°) to outer parts (28.04°± 7.35°).

All other parameters showed no significant differences from inner to outer sites when using the simple minimization procedure (function 5). However, it is necessary to focus on the small sample size of our pilot study. In order to ensure that all the parameters are not statistically different with a margin of error of less than 10% (type II error), a power study shows that it would be necessary to increase the sample size to at least 200. An in-depth human study with a larger number of samples and at least six different spines should be considered.

Measurement of transverse dilatations also highlighted an unusual increase in specimen width and thickness at low strain values for 16 samples in the e RR e RR e RR direction and for 18 samples in the e ZZ e ZZ e ZZ direction over the 20 samples. Such an unexpected behavior has only been reported once on the AF by [START_REF] Acaroglu | Degeneration and aging affect the tensile behavior of human lumbar anulus fibrosus[END_REF], who is one of the rare authors to have combined quasi-static conditions with the observation of strain in one transverse direction to traction. Interestingly, such an auxetic behavior has also been observed on skin by [START_REF] Veronda | Mechanical characterization of skin-finite deformations[END_REF] and more generally on composite materials by [START_REF] Clarke | Negative poisson's ratios in angle-ply laminates : theory and experiment[END_REF] and [START_REF] Alderson | Auxetic materials[END_REF], Indeed, the low values of the initial fiber rigidity parameter k 1 compared to the shear modulus µ permits the mechanical response of the isotropic part to be preponderant for low strain values. Therefore, the compressible constitutive model proposed in this study was able to qualitatively reproduce the initial experimental swelling and the following shrinkage of the AF samples.

Quantitatively, as our results show, the model used within this study still has some limitations in its ability to fully reproduce the amplitude of the initial experimental swelling and the shrinkage that follows. The modeling part could be greatly enhanced using multi-scale and nonlinear homogenization approaches that take into account the reorientation and instabilities that may occur in the micro structure, as proposed by [START_REF] Ayyalasomayajula | A computational model for understanding the micro-mechanics of collagen fiber network in the tunica adventitia[END_REF]; [START_REF] Greco | A study of stability and bifurcation in micro-cracked periodic elastic composites including self-contact[END_REF][START_REF] El Nady | Computation of the homogenized nonlinear elastic response of 2d and 3d auxetic structures based on micropolar continuum models[END_REF].

During our experimental procedures, we made sure that we applied the same NaCl concentration and bath temperature and imposed the same preload of 0,01 N on each sample. Nevertheless, one limitation is the variability of sample dimensions (L θθ × L RR × L ZZ = 5.29 ± 0.5 × 3.49 ± 0.87 × 1.36 ± 0.4 mm 3 ), which is difficult to avoid, and led to differences in the initial pre-stress. This inconsistency can probably explain the large dispersion in both radial and vertical initial swelling. Another limitation is the assumption of the vertical plane symmetry of the two families of fibers during the inverse analysis process. This assumption is usually made when modeling disc mechanics [START_REF] Malandrino | Regional annulus fibre orientations used as a tool for the calibration of lumbar intervertebral disc finite element models[END_REF]; [START_REF] Reutlinger | Specimen specific parameter identification of ovine lumbar intervertebral discs : On the influence of fibre-matrix and fibre-fibre shear interactions[END_REF]) but does not seem to be valid in all cases [START_REF] Disney | Visualising the 3d microstructure of stained and native intervertebral discs using x-ray microtomography[END_REF]). A further analysis could be carried out, comparing our inverse analysis to image analysis in order to evaluate fiber angle measurement [START_REF] Matcher | The collagen structure of bovine intervertebral disc studied using polarization-sensitive optical coherence tomography[END_REF]; [START_REF] Disney | Visualising the 3d microstructure of stained and native intervertebral discs using x-ray microtomography[END_REF]) but also to increase the sample quantity.

Conclusion

The results of this study are particularly relevant to models of both the heterogeneous and the compressible behavior of the annulus fibrosus. Further studies are required, using multi-relaxation tensile tests and combining measurements of transverse dilatations with an adequate inverse analysis in order to capture the compressible heterogeneous behavior of the AF of human IVD. Our present study suggests that if the mechanical behavior of the human IVD proves to be equivalent to that of the porcine IVD, structural parameters such as the fiber angle could be sufficient to explain the

Figure 1 :

 1 Figure 1: Schematic of the sheet specimen extraction protocol. Spinous and transverse

  15 M NaCl bath for 3 hours. The specimens (L θθ × L RR × L ZZ = 5.29 ± 0.5 × 3.49 ± 0.87 × 1.36 ± 0.4 mm 3 ) were then glued to coaxial aluminum traction trays with cyano-acrylate adhesive (Fig 2) and immersed in a 37°C controlled bath of 0.15 M NaCl. An initial tensile force of 0.01 N was applied to each sample corresponding to the traction tray weight. Ten cycles at 0.1 Hz with 1 mm displacement were imposed on the specimens as a preconditionning stage. The specimens were then tested under displacement control by incrementally applying 0.1 mm displacement eight times. The specimens were held at each strain level for 2 minutes for the first step and for 20 minutes for the last step to allow the circumferential stresses to relax to a steady-state value (Fig 3). The traction force, the displacement, and two series of images (in two perpendicular directions) were acquired simultaneously at 0,5 Hz. The force was measured with a miniature S Beam load cell (Futek Inc., model 10 N LSB210 , Irvine, USA). The cameras (IDS-UI-3360CP-M-GL Rev 2.0. 16 bits, Imaging Development Systems GmbH, Obersulm, Germany) were positioned perpendicularly to the tensile direction to observe both transverse planes and Fujifilm-Fujinon 1 :1.4/16 mm lenses were set on cameras to give a spatial resolution of 16 µm/px. The circumferential stress σ θθ was calculated by dividing the axial force by the deformed cross-sectional area of the tested specimen.

Figure 2 :

 2 Figure 2: Tensile test principle in the circumferential direction with an imposed displacement. Force, displacement, and two series of images were acquired simultaneously in the ( e θθ , e ZZ e θθ , e ZZ e θθ , e ZZ ) and in the ( e θθ , e RR e θθ , e RR e θθ , e RR ) planes. The fibers are oriented in the ( e θθ , e ZZ e θθ , e ZZ e θθ , e ZZ ) plane.

Figure 3 :

 3 Figure 3: Example experimental data set. The specimens are tested under imposed displacement by incrementally applying 0.1 mm circumferential displacement. Displacement is held 2 minutes for the first step and 20 minutes for the last step in order to let the load relax to a steady-state value. Only the red dot points are used for the elastic stress-strain curves.

  ; Cortes & Elliott (2012); Oomen et al. (2018)) :

  tural parameter identification : a) the force in the direction of traction F exp θ , b) the dilatation in the radial direction λ exp R , and c) the dilatation in the vertical direction λ exp Z . Here, F num θ , λ num R , and λ num Z represent the predicted

Figure 4 :

 4 Figure 4: Sample-specific numerical model used in the material parameter identification. For each numerical sample, only a quarter of the geometry was modeled due to the symmetry plans of the geometry.

  software. Nonparametric tests had to be applied because normality tests(Shapiro-Wilk) showed that some of the data groups do not have normal distributions. Data have to be treated as paired samples since two values, related to the radial position, come from the same AF. The studied values are the circumferential stress, the radial and vertical dilatations, and all the mechanical and structural parameters. We performed a Mann-Whitney U test for (a) comparison of the mean values at each strain data point of the circumferential stress σ θθ , the vertical dilatation λ ZZ , and the radial dilatation λ RR of the outer position with those of the inner position, and (b) comparison of the mean values of the mechanical (λ, µ, κ, k 1 , k 2 ) and structural (α i α i α i ) parameters of the outer specimens with those of the inner specimens. A Bonferroni correction on alpha to account for repeating the hypothesis at each strain data point was performed for the statistical analysis on the the circumferential stress σ θθ , the vertical dilatation λ ZZ , and the radial dilatation λ RR .

Figure 5 :

 5 Figure 5: Circumferential stress vs. circumferential dilatation for interior (blue lines) and exterior samples (geen lines) as well as median values. The filled areas around the median values represent the interquartile range.

Figure 6 :

 6 Figure 6: Vertical dilatation versus circumferential dilatation for interior and exterior samples as well as median values. The filled areas around the median values represent the interquartile range.

Figure 7 :Figure 8 :

 78 Figure 7: Radial dilatation versus circumferential dilatation for interior and exterior samples as well as median values. The filled areas around the median values represent the interquartile range.

Figure 9 :

 9 Figure 9: Comparison of the circumferential stress versus circumferential dilatation for experimental and numerical T1-T2 sample. Solid lines represent numerical predictions and dashed lines represent experimental results. The results of two minimization procedures are shown : the minimization on force only (green lines) gives outstanding results on the stress-strain behavior whereas the minimization on both force and dilatation (blue lines)

  and is mainly explained by the fibrous structure of such material. By means of experimental estimation of the transverse dilatations λ RR and λ ZZ during the quasi-static tensile test, we show that all the samples had a large compressible behavior. As proposed by several authors[START_REF] Elliott | Anisotropic and inhomogeneous tensile behavior of the human anulus fibrosus : experimental measurement and material model predictions[END_REF];[START_REF] Cortes | Extra-fibrillar matrix mechanics of annulus fibrosus in tension and compression[END_REF];[START_REF] Guerin | Quantifying the contributions of structure to annulus fibrosus mechanical function using a nonlinear, anisotropic, hyperelastic model[END_REF]), this compressible behavior of the IVD can be interpreted as the volume change of the porous solid phase, considering a biphasic media. To our knowledge, this is the first time that this behavior has been reported on the AF in the e RR e RR e RR and e ZZ e ZZ e ZZ directions simultaneously while the sample is loaded in the e θθ e θθ e θθ direction. Furthermore, the additive decomposition of the strain energy function into an isotropic compressible part and an anisotropic part allows the numerical model with optimized parameters to predict the nonlinear transverse dilatation behaviors.

  ii e ii e ii ⊗ e ii e ii e ii ) and e ii e ii e ii = e θθ e θθ e θθ , e ZZ e ZZ e ZZ , or e RR e RR e

and the circumferential (λ θθ ), radial (λ RR ), and vertical dilatations (λ ZZ ), with λ ii = C C C : ( e RR (see figure 2 for a representation of the different vectors). The circumferential and radial dilatations were estimated through incremental digital image correlation (DIC) from one relaxation step to the next. Digital image correlation requires an increased

Table 1 :

 1 Median values ± interquartile range of the parameters for both exterior and interior samples for two different objective functions. Only the fiber angle is significantly different (p < 0.05) for both objective functions.

		Force minimization	Force and dilatation
	Parameters
				minimization
		(function 4)	(function 5)
		Exterior	Interior	Exterior	Interior
	Bulk Modulus λ (MPa)	0.034 ± 0.077	0.024 ±0.036	0.016 ± 0.075 0.021 ± 0.066
	Shear Modulus µ (MPa)	0.22 ± 0.327	0.14 ± 0.18	0.14 ± 0.06	0.11 ± 0.15
	k 1 (MPa)	0.00065 ± 0.0026 0.00226±0.0074 0.043 ± 0.061	0.048±0.087
	k 2	611± 384	896±478	1216± 1237	896±2522
	κ	0.091± 0.033	0.15 ± 0.043	0.12± 0.061	0.18 ± 0.09
	α α α (°)	25.67 ± 7.93	35.61±5.35	28.04 ± 7.35	34.11±12.05
	NRMSE (%)	σ θθ	λ ZZ	λ RR	J
	Force minimization	10.35±8.84 36.01±21.41 44.68±20.61 45.20±17.70
	Force and dilatation minimization 15.67±4.71 15.05 ±18.11 33.07±24.31 39.17±20.97

Table 2 :

 2 Comparison of function minimization on the median NRMSE ± interquartile range for all measured data compared to model predictions. The lower NRMSE values demonstrate the ability of the model to predict the experimental responses.

  ). Moreover, nonconventional tensile tests were used, fixing the samples only on their lower faces. Contrary to the uniaxial tensile test, this test produces a nonconventional strain state of the sample, inducing a richness of loading. The sample is in traction in the e θθ e θθ e θθ direction, but also in shear in the e θθ e θθ e θθ -e ZZ e ZZ e ZZ and e θθ e θθ e θθ -e RR e RR e RR planes. This induces a variation of each invariant I 1 , J, I 4 , and I 6 during the test, which is mandatory to

	identify the parameters of the nonlinear anisotropic hyperelastic law. The
	use of a nonconventional tensile test, combined with the presence of the
	transverse dilatations in the cost function, are probably the key ingredients
	permitting the optimization process to robustly identify the six parameters,
	including the fiber angle.

Acknowledgements

The authors gratefully acknowledge Stephan Devic (University of Montpellier) for manufacturing the tensile test machine.

spatial variations of AF quasi-static circumferential tensile behavior. While [START_REF] Goel | Interlaminar shear stresses and laminae separation in a disc : finite element analysis of the l3-l4 motion segment subjected to axial compressive loads[END_REF] and [START_REF] Lee | Impact response of the intervertebral disc in a finite-element model[END_REF] used the same mechanical properties in the whole AF, other authors used gradients of mechanical properties [START_REF] Eberlein | Multi-segment fea of the human lumbar spine including the heterogeneity of the annulus fibrosus[END_REF]). Therefore, our study shows that numerical models of the AF of the IVD could use constant mechanical properties with only variation of the fiber angle from the inner to the outer periphery [START_REF] Schmidt | Application of a new calibration method for a threedimensional finite element model of a human lumbar annulus fibrosus[END_REF]).

Competing interest

We have no competing interests.