Mohamed Amine Chetoui
email: mohamed-amine.chetoui@centrale-marseille.fr

Dominique Ambard

Patrick Cañadas

Pascal Kouyoumdjian

Pascale Royer

Simon Le Floc'

Patrick Canãdas

Simon Le Floc'h

Impact of extracellular matrix and collagen network properties on the cervical intervertebral disc response to physiological loads: A parametric study

Keywords: Intervertebral disc, FE analysis, Fibrous soft tissue, Porous media, Cervical spine kinematics

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Tissue engineering, replacement and regeneration techniques are increasingly used in the treatment of intervertebral disc diseases. While synthetic material replacements are still limited by the biocompatibility hurdle, biomaterials designed for annulus fibrosus (AF) and nucleus pulposus (NP) repair and replacement have managed to mimic IVD biology and have shown promising results [START_REF] Zhuang | Construction of tissue-engineered composite intervertebral disc and preliminary morphological and biochemical evaluation[END_REF][START_REF] D'este | Lessons to be learned and future directions for intervertebral disc biomaterials[END_REF]. However, without focusing on mechanical compatibility, these biomaterials-based solutions have not yet shown a long-term performance [START_REF] D'este | Lessons to be learned and future directions for intervertebral disc biomaterials[END_REF][START_REF] Buckley | Critical aspects and challenges for intervertebral disc repair and regeneration-harnessing advances in tissue engineering[END_REF]. Mechanical compatibility is not less critical than biocompatibility since it permits to restore biomechanical behavior of the motion segment, and thus reduces re-herniation and promotes longevity [START_REF] D'este | Lessons to be learned and future directions for intervertebral disc biomaterials[END_REF][START_REF] Sharabi | Towards intervertebral disc engineering: Bio-mimetics of form and function of the annulus fibrosus lamellae[END_REF]. Therefore, a sophisticated comprehension of the biomechanical behavior of the intervertebral disc is invaluable to enhancing the performance of biomaterial-based techniques.

The intervertebral disc (IVD) biomechanical behavior is complex and it can not be explored using only experiments or analytical methods [START_REF] Guo | Impact of material properties of intervertebral disc on dynamic response of the human lumbar spine to vertical vibration: A finite element sensitivity study[END_REF]. Spine and disc biomechanics are commonly studied using finite element (FE) analysis. FE models developed for the IVD are increasingly improved. However, they become more difficult to implement and validate.

Recent models take into account complex structural phenomena taking part in this soft tissue such as the osmotic role of proteoglycans and the mechanical contribution of collagen fibers including fiber cross-links, fiber/matrix interaction, and interlamellar behavior (e.g. [START_REF] Hollingsworth | Modeling shear behavior of the annulus fibrosus[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][START_REF] Mengoni | Derivation of inter-lamellar behaviour of the intervertebral disc annulus[END_REF][START_REF] Masni-Azian | Biomechanical investigation on the influence of the regional material degeneration of an intervertebral disc in a lower lumbar spinal unit: A finite element study[END_REF]). These models generally managed to provide a good agreement with experimentation. However, the uniqueness of the identified mechanical properties is not guaranteed given their relatively large number and their possible interdependency [START_REF] Oreskes | Evaluation (not validation) of quantitative models[END_REF][START_REF] Chetoui | Assessment of intervertebral disc degeneration-related properties using finite element models based on ρ h -weighted MRI data[END_REF]. Furthermore, the impact of the experimentation-related uncertainties on the identification studies remains significant and it is reflected in the large scatter of the IVD mechanical properties values. Most of the biomechanical experimental tests used in the parameter identification are performed ex-vivo. Cadaveric tissue, generally used to perform biomechanical ex-vivo experimental tests, undergoes several changes in properties related to the conservation conditions which leads to variability in the experimental data [START_REF] Olson | Designing a biomechanics investigation: choosing the right model[END_REF]. In addition, the relevance of each parameter depends on the load type and the nature of the experimental data used for the identification.

For example, the micro-channels configuration of a porous medium may increase the apparent stiffness of the tissue in the case of confined consolidation compared to the unconfined one [START_REF] Urcun | Cortex tissue relaxation and slow to medium load rates dependency can be captured by a two-phase flow poroelastic model[END_REF].

Computational validation issues could be treated by several approaches. Multi-objective optimization and the use of huge experimental databases could be suggested. However, these solutions tend to complicate the studies and raise their time-memory cost. Different approaches, based upon the exploitation of quantitative MRI for the construction and the validation of FE models were introduced [START_REF] Stadelmann | Integrating MRI-based geometry, composition and fiber architecture in a finite element model of the human intervertebral disc[END_REF][START_REF] Chetoui | Assessment of intervertebral disc degeneration-related properties using finite element models based on ρ h -weighted MRI data[END_REF]. This technique seems efficient but is limited by the MRI resolution and needs to be clinically validated. Recent studies introduced probabilistic and deterministic deep learning approaches to reduce the cost of complex behavior tissue models including porous media, hyperelastic anisotropic tissue, and multiscale models.

Despite their promising results in model reduction and simulation acceleration, these methods present some limitations related to the handle with irregular mesh and some particular boundary conditions [START_REF] Thakolkaran | NN-EUCLID: Deep-learning hyperelasticity without stress data[END_REF][START_REF] Krokos | A Bayesian multiscale CNN framework to predict local stress fields in structures with microscale features[END_REF][START_REF] Deshpande | Probabilistic deep learning for real-time large deformation simulations[END_REF]. Therefore, sensitivity studies remain an essential tool to improve and validate IVD FE models [START_REF] Anderson | Verification, validation and sensitivity studies in computational biomechanics[END_REF][START_REF] Guo | Impact of material properties of intervertebral disc on dynamic response of the human lumbar spine to vertical vibration: A finite element sensitivity study[END_REF]. Studying the impact of varying the parameter values permits to better understand the IVD model biomechanics and manage the validation studies. It allows sorting parameters by relevance and fixing the values of those with no significant effect, simplifying the study by focusing on the most important parameters and therefore helping to obtain a more accurate result.

Several sensitivity analyses have been conducted on lumbar spine FE models. The impact of structural, mechanical, and morphological properties of vertebrae, IVD, and ligaments in the lumbar spine was studied under different loads [START_REF] Malandrino | Statistical factorial analysis on the poroelastic material properties sensitivity of the lumbar intervertebral disc under compression, flexion and axial rotation[END_REF][START_REF] Jebaseelan | Sensitivity studies of pediatric material properties on juvenile lumbar spine responses using finite element analysis[END_REF][START_REF] Yang | Effect of collagen fiber orientation on intervertebral disc torsion mechanics[END_REF][START_REF] Zander | Impact of material and morphological parameters on the mechanical response of the lumbar spine -a finite element sensitivity study[END_REF][START_REF] Guo | Impact of material properties of intervertebral disc on dynamic response of the human lumbar spine to vertical vibration: A finite element sensitivity study[END_REF]. The obtained results have contributed to the understanding of lumbar spine biomechanics and the studying of some FE models accuracy. However, fewer sensitivity studies have been carried out for the cervical spine and they are limited on the sensitivity of the stiffness of the spinal components, assumed to be linear elastic, on spine biomechanics [START_REF] Kumaresan | Finite element analysis of the cervical spine: A material property sensitivity study[END_REF][START_REF] Ng | Statistical factorial analysis on the material property sensitivity of the mechanical responses of the C4-C6 under compression, anterior and posterior shear[END_REF]. The current study represents a sensitivity analysis of C6-C7 IVD mechanical properties using a nonlinear poroelastic FE model with fiber reinforcement. The objective is to determine and analyze the sensitivity of the model to both the matrix mechanical properties and the collagen network configuration. For this purpose, physiological loads are used in the simulations. By prioritizing the contribution of the model properties in the IVD biomechanics for different load types, this work aims to create a benchmark for numerical studies.

Material and methods

Constitutive formulation

The definition of quantities of interest permits to perform the mesh sensibility test and to choose the constitutive model of the study [START_REF] Duprez | Quantifying discretization errors for soft tissue simulation in computer assisted surgery: A preliminary study[END_REF]. In our case, we were interested in the maximum stress linked to the solid phase, the fluid flow linked to the fluid phase and the apparent stiffness which is a global quantity linked to the entire domain.

The choice of the constitutive model has been based on the quantities of interest and on the mechanical aspects studied in the literature. In vivo measurements have shown that the cervical IVD undergoes large strain during physiological motion [START_REF] Anderst | Cervical Spine Disc Deformation During In Vivo Three-Dimensional Head Movements[END_REF]. It was also shown that the collagen fibers do not contribute only to the tissue response to mechanical loads but also to its swelling behavior [START_REF] Sharabi | The Mechanical Role of Collagen Fibers in the Intervertebral Disc[END_REF][START_REF] Yang | Swelling of fiber-reinforced soft tissues is affected by fiber ori-entation, fiber stiffness, and lamella structure[END_REF]. On the other hand, The compressibility of the extracellular matrix (ECM) has been largely investigated and discussed. Although several finite element studies have adopted a nearly-incompressible model for the extracellular matrix, recent experimental studies have highlighted the auxetic compressible behavior of the annulus fibrosus [START_REF] Dusfour | Mesenchymal stem cells-derived cartilage micropellets: A relevant in vitro model for biomechanical and mechanobiological studies of cartilage growth[END_REF][START_REF] Derrouiche | The two poisson's ratios in annulus fibrosus: Relation with the osmo-inelastic features[END_REF]. These results revealed the importance of studying the IVD behavior in large deformation, taking into account the annulus anisotropy, the compressibility of the extracellular matrix (ECM), and the fluid flow in the porous network.

The constitutive formulation implemented for this study is based on the biphasic swelling model [START_REF] Mow | Biphasic creep and stress relaxation of articular cartilage in compression: Theory and experiment[END_REF][START_REF] Lanir | Biorheology and fluid flux in swelling tissues. I. Bicomponent theory for small deformations, including concentration effects[END_REF][START_REF] Wilson | A comparison between mechanoelectrochemical and biphasic swelling theories for soft hydrated tissues[END_REF]. The IVD is assimilated to a superposition of two immiscible and isothermal phases: a porous solid skeleton describing the fiber-reinforced extracellular matrix saturated by an intrinsically incompressible fluid. The total Cauchy stress σ is the summation of the solid effective stress tensor σ e and the interstitial fluid stress derived from the hydrostatic pressure p and the osmotic pressure ∆π:

σ = σ e -(p + ∆π)1 (1)
The hydrostatic pressure p is determined by Darcy and mass conservation laws while respecting the intrinsic incompressibility assumption. The permeability k of the medium is taken isotropic and strain-dependent as proposed in [START_REF] Argoubi | Poroelastic creep response analysis of a lumbar motion segment in compression[END_REF]. The osmotic pressure ∆π of the IVD is expressed in terms of the fixed charge density of the proteoglycans C f c which is also taken strain-dependent [START_REF] Wilson | A comparison between mechanoelectrochemical and biphasic swelling theories for soft hydrated tissues[END_REF].

∆π = RT (φ i C 2 f c + 4C 2 e -2φ e C e ) (2)
C f c = C f c0 ϕ 0 ϕ 0 -1 + J (3)
where ϕ 0 is the initial porosity, J is the deformation gradient determinant, C f c and C f c0 are the current and the initial fixed charge densities respectively, C e the external salt concentration, R the universal gas constant, T the absolute temperature, and φ i and φ e respectively the internal and the external osmotic coefficients. The extracellular matrix of the medium is defined by its total strain energy density W . W was composed of a compressible Neo-Hookean isotropic part W N H for the non-fibrillar matrix and an anisotropic part W f i describing the contribution of each family of fibers i:

W = W N H (C) + i W f i (C, ⃗ e i ) (4)
W N H = µ 2 (I 1 -3) -µ ln (J) + λ 2 (ln (J)) 2 (5)
where C is the right Cauchy-Green strain tensor, I 1 is its first invariant and µ and λ are the parameters of Lamé. The vector ⃗ e i represents the director vector of the fiber family i orientation. The anisotropic part was taken null in the NP and the CEP. Two families of fibers oriented by ±α were considered in the AF. The fiber strain energy density of the i th fibers population (i=1,2) was defined according to [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientation[END_REF].

W f i = a i 2b i e b i[ (1-3K)I f i +KI 1 -1] 2 -1 (6)
where a i (MPa) and b i (unitless) are the fiber rigidity and non linearity coefficient, K is their dispersion and I f i is a fiber direction invariant written as

I f i = C : (⃗ e i ⊗ ⃗ e i ) (7)
The same mechanical parameters were taken for both families of fibers: a i = a AF and

b i = b AF .

Finite element model

The unknowns of the equation system were the hydrostatic pressure p and the threecomponent displacement field. We have used Lagrange quadratic (P2) and linear (P1) polynomial functions for the interpolation of the displacement and the pressure, respectively.

For this study, we have developed a three-dimensional finite element model of human C6-C7 IVD. We have constructed a parametric geometry based on six experimental studies conducted on the cervical spine morphology [START_REF] Gilad | A study of vertebra and disc geometric relations of the human cervical and lumbar spine[END_REF][START_REF] Frobin | Vertebral height, disc height, posteroanterior displacement and dens-atlas gap in the cervical spine: Precision measurement protocol and normal data[END_REF][START_REF] Wiegand | Cervical spine geometry correlated to cervical degenerative disease in a symptomatic group[END_REF][START_REF] Pitzen | Variation of endplate thickness in the cervical spine[END_REF][START_REF] Lou | Geometry of inferior endplates of the cervical spine[END_REF][START_REF] Yu | Ranges of cervical intervertebral disc deformation during an in vivo dynamic flexion-extension of the neck[END_REF]. The geometry is composed of a central elliptical nucleus pulposus (NP) representing 26% of the total volume, a peripheral annulus fibrosus (AF) surrounding the NP and two cartilaginous endplates (CEP) (Fig. 1).

The geometry was meshed using the software Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in preand post-processing facilities[END_REF]. A mesh refinement study was performed using 4 20-node hexahedral element meshes with different element number to check the convergence of our quantities of interest against element size. Local refinement was not considered given the absence of geometrical irregularities, high-stress gradient, or common stress concentration areas for the studied loads. The selected mesh has been composed of 5940 hexahedral elements (26106 nodes). The model was implemented in the open-source software LMGC90 [START_REF] Dubois | Lmgc90[END_REF].

Definition of the fiber orientation

The fiber orientation angle α decreases from the inner to the outer lamella of the AF [START_REF] Holzapfel | Single lamellar mechanics of the human lumbar anulus fibrosus[END_REF].

To respect this property, we have defined the fiber orientation using local spatial distributions progressing continuously from the inner to the outer AF. Two parameters were defined, α i and α o , the fiber angle to the transverse plane, at the inner and at the outer lamella, respectively (Fig. 2.B). Given the shape of the IVD geometry, it was necessary to define a local coordinate system centered in the NP center O N P and satisfying two criteria: i) For each element of the AF, fibers are defined in lamellae plane, ii) Elements of the same lamella should have the same fiber orientation. At this stage, the shape of the annular mesh has been useful. To define this coordinate system we started with the resolution of a thermal problem in which we have imposed two different temperatures in the lateral boundary of the AF (T AF ) and in the NP (T N P ) with T AF < T N P (Fig. 2.A). The heat flux vector in each node defines its local radial vector ⃗ e r AF . Then, we defined the axial vector ⃗ e a AF which corresponds to the Z-axis of the global coordinate system. The tangential vector ⃗ e t AF was determined by the cross product (⃗ e a AF ∧ ⃗ e r AF ). Finally, the vectors ⃗ e 1 and ⃗ e 2 were defined in each node according to this relation:

⃗ e 1 = cos(α).⃗ e t AF + sin(α).⃗ e a AF ( 8
) ⃗ e 2 = -cos(α).⃗ e t AF + sin(α).⃗ e a AF (9)

Model validation

A validation process was conducted to find the reference parameter vectors to perform the sensitivity study. We have compared the behavior of our model in lateral bending, flexion, and axial rotation to experimental previous works by reproducing the applied load described in [START_REF] Yoganandan | Leveldependent coronal and axial moment-rotation corridors of degeneration-free cervical spines in lateral flexion[END_REF][START_REF] Wheeldon | Experimental flexion/extension data corridors for validation of finite element models of the young, normal cervical spine[END_REF][START_REF] Yoganandan | Normative segment-specific axial and coronal angulation corridors of subaxial cervical column in axial rotation[END_REF]. These experimental studies provide functions describing the C6-C7 IVD rotation when it is exposed to a 0 to 2Nm pure moment. To compare the behavior of the current model to experimental curves, we have performed 4 tests for each motion by applying 0.5, 1, 1.5, and 2Nm moments and determined the IVD orientation angle for each test. Then we have interpolated these points to find the numerical evolution function of each motion (lateral bending, flexion, and axial rotation). The interpolation function used here has been defined in the experimental studies and has the following form:

θ = γ v + α v ln(β v M + 1)
. where θ and M are the rotation angle and the applied moment, respectively, and γ v , α v and β v are the function parameters. In our case, we have taken γ v = 0 for all the loads to start from 0 • of rotation when a null moment is applied. We have manually varied the model parameters to finally find a parameter vector that permits to closely reproduce the experimental behavior of the C6-C7 IVD (Fig. 3). This parameter vector will constitute the reference vector of the parametric study called for the rest of this paper "the basic model". Here, the basic model is not necessarily the best fit combination with a minimized error compared to the experimental curve. However, we are sure that, with the retained combination, the model reproduces the experimental mechanical response used in this validation step.

Boundary conditions for the parametric study

The simulations were performed under uniaxial compression and three physiological independent loads of the cervical spine: lateral bending, axial rotation, and flexion. Each load has been applied in a 3-step sequence (Fig. 4). The first step was a 2-hour preconditioning step. The bottom surface has been constrained in all directions, the remaining boundaries being constraint-free. This step permits to obtain the unloaded IVD equilibrium in which the initial osmotic pressure distribution and the unloaded IVD equilibrium swelling were established from the fixed charge density and the porosity initial fields (C f c0 and ϕ 0 table 2).

The second step permits to simulate the head weight. A 100N creep-compressive force was applied on the top surface of the upper CEP in 8 minutes then remains constant until the equilibrium. Finally, starting from the result of the compression-creep step, we have applied a displacement-driven load to simulate the C6-C7 translations and rotations during a physiological cervical spine lateral bending, axial rotation, and flexion (Fig. 4). Data for C6-C7 kinematics in the cervical segment are obtained from in vivo measurements taken from the literature (table 1). To reproduce the experimental measurements, we have used a y-z-x sequence (corresponding to flexion/extension-axial rotation-lateral bending) to apply the C6-C7

rotations during spine movement. The center of rotation is the center of the most posterio-inferior point of the subjacent vertebrae (C7) for lateral bending and axial rotation [START_REF] Ishii | Kinematics of the cervical spine in lateral bending: In vivo three-dimensional analysis[END_REF][START_REF] Lin | In vivo three-dimensional intervertebral kinematics of the subaxial cervical spine during seated axial rotation and lateral bending via a fluoroscopy-to-CT registration approach[END_REF].

For the flexion, the rotation center is the center of the subjacent vertebral body [START_REF] Crawford | Methods for determining spinal flexion/extension, lateral bending, and axial rotation from marker coordinate data: Analysis and refinement[END_REF][START_REF] Anderst | Three-dimensional intervertebral kinematics in the healthy young adult cervical spine during dynamic functional loading[END_REF].

Studied parameter range

The parametric study has been carried out by varying 10 parameters (table 2). For the NP and the AF ground substance tissues, the varied parameters were the shear modulus µ defining the ECM stiffness, the first Lamé parameter λ defining its compressibility, and the initial permeability k 0 . Collagen fiber network properties have also been studied by varying their rigidity a, their nonlinearity b, and their directions (α o and α i ). We have chosen the interval of variation of each parameter to be always in the range of values in previous numerical studies found in the literature (table 2). This choice permits testing the sensitivity of the IVD mechanical behavior to model parameters in extended ranges used in previous studies. It is important to mention that even if some cited studies do not use the same constitutive law as the actual work, they were used to define variation range for parameters equivalent to ours in their formulation. For example, we have used the value of a AF in [START_REF] Jacobs | Validation and application of an intervertebral disc finite element model utilizing independently constructed tissue-level constitutive formulations that are nonlinear, anisotropic, and time-dependent[END_REF] even if this study does not take into consideration the fiber dispersion defined by K AF in our model. We also calculated the equivalent λ N P and µ N P from Young modulus and Poisson's ratio provided in [START_REF] Ferguson | Fluid flow and convective transport of solutes within the intervertebral disc[END_REF] using elastic parameter conversion formulae.

For each parameter, the reference value was managed to be at the center of the studied range. The higher and the lower values were chosen in a way to be equally distant from the reference value and to form an interval that covers as many previously used values as possible.

The first simulation has been performed using the reference parameter values (basic model).

Then, for each test, we have varied only one parameter. Only the two fiber orientation angles have been simultaneously varied. To make this simultaneous variation simple we have fixed the mean fiber direction at 35 o and defined the parameter ∆α which represents the difference of the mean angle from inner and outer angles: α i = ±(35 + ∆α) and α o = ±(35 -∆α).

Finally, 37 tests have been performed for each load type. A total of 1320h of computing time was needed in this study (about 220h with parallel computation) using an HPC cluster with 28 cores and 128GB of RAM.

The studied outputs for the creep-compression step were the maximal NP pressure, the maximal normal stress, the axial and radial outflow rates, the final volume change, and the final axial displacement of the upper surface. For the physiological loads, we have studied the maximal normal and shear stress, the axial and radial outflow rates, the final volume change, and the final corresponding moment applied on the upper IVD surface. We have also studied the impact of model parameters on a global output which is the apparent stiffness for the different loads. The apparent stiffness for the rotational loads was calculated about the center of rotation (Appendix).

Results

Compression step

During the compression step, we can differentiate two stages. The first one is the establishment of the applied force. At this stage, the axial displacement pressurizes the fluid content due to the low permeability of the tissue. The IVD swells following this pressurization until the total establishment of the force. At this instant, the pressure and the radial swelling reach their maximal values. The second stage begins when the applied force is stabilized. In this stage, the water expulsion from the IVD becomes more important, the swelling decreases and the displacement of the loaded CEP continues until the equilibrium. We measured the influence of modifying the model parameters on the NP pressure, the normal stress, and the axial and radial flow rates at the end of the first stage. We studied also two equilibrium results which are the volume loss caused by the water expulsion and the final displacement of the loaded face which will determine the apparent stiffness of the IVD. With the basic model, the relative change in IVD height was 4.8% and the relative change in its volume was 4.3%.

Figure 5 shows the effect of modifying the model parameters on the IVD response to compression. Results show that the shear modulus of the AF ground substance µ AF is the most influential parameter on IVD response to compressive load. The AF compressibility λ AF , affects significantly the radial flow rate. This parameter has a lower influence on the axial flow rate, the final displacement, and the volume loss but it contributes clearly to their determination. The NP and the AF permeabilities k0 N P and k0 AF play also an important role in the pressurization stage. The pore pressure is very sensitive to the NP permeability and the radial flow rate is highly influenced by the AF permeability. However, the permeabilities do not contribute to the equilibrium state. Apart from the fiber rigidity, no noticeable influence was found on fiber parameters. The reference apparent stiffness obtained with the basic model was about 337N mm -1 . This value depends on model parameters in a nearly linear way with high sensitivity to the AF stiffness µ AF and less noteworthy dependency on fiber rigidity a AF and AF compressibility λ AF . No significant effect on the IVD apparent stiffness was found for the other parameters. Contrary to compression, we have noticed a significant dependency of the IVD behavior on fiber parameters in lateral bending. The fiber non-linearity parameter b AF has a crucial role in the determination of the shear stress, the lateral bending moment, and the volume change.

Lateral bending

Modifying the fiber rigidity a AF has a noticeable effect on the same outputs. Although fiber direction parameter ∆α has not the same importance in the IVD response to lateral bending, we have found that varying this parameter may affect the shear stress and the radial outflow rate. On the other hand, apart from µ N P , which slightly affects the radial flow rate, the NP parameters varying in the study range has not a real influence on the IVD response to lateral bending. The apparent stiffness in lateral bending found with the basic model was about 16N mrad -1 . This value is affected by modifying b AF in the first row then AF stiffness µ AF and fiber rigidity a AF . The dependency of the LB apparent stiffness to these parameters is nonlinear.

Flexion

Axial rotation

The results of the effect of modifying model parameters on IVD response to axial rotation are shown in figure 8. As for the other types of load, varying µ AF has a high effect on the maximal stress and the flow rates in axial rotation. In addition, the permeability of the AF slightly contributes to the determination of the radial flow rate where the role of NP permeability remains negligible. The fiber parameters b AF and a AF role remains crucial in the determination of shear stress, the axial rotation moment, and the volume change. We note that fiber orientation contributes more significantly to the IVD response to axial rotation than its contribution to flexion and lateral bending. Similarly to the two latter motion results, the apparent stiffness value is dominated by b AF , µ AF and a AF with a more noticeable effect of the NP shear modulus µ N P . The axial rotation apparent stiffness value found with the basic parameters was 2.83N mrad -1 .

The apparent stiffness values found with the three physiological rotations are significantly smaller than the apparent stiffness computed with 2N m pure moments applied in the validation study with (25.8N mrad -1 in lateral bending, 15.3N mrad -1 in flexion, and 24.9N mrad -1 in axial rotation). This comparison reveals that coupling the different rotations and translations during a physiological movement, together with the specific anisotropic hyperelastic properties of the IVD, help reduce the resulting moment that is necessary to turn the head.

Discussion

This study was aimed to answer a key question that precedes any study treating the IVD mechanical modeling: given the data set and objectives, to which mechanical parameters specific attention should be paid? The answer to this question permits reducing the model numerical cost by considering the appropriate assumptions or by reducing the number of parameters in the identification studies. To prioritize model parameters we have conducted a parametric study using an anisotropic poro-hyperelastic model of the IVD and we have examined the effect of varying 10 parameters, in the literature range, on the model response to physiological day-to-day loads.

Our methodology in the definition of parameters range has been oriented to cover the most used values in literature. Therefore, it does not consist of a sensibility study where we try to find the effect of each parameter on a limited interval but of a parametric study where we have examined the parameters contribution to the IVD behavior based on the variety of these parameter values used by previous works. Although this choice imposes a very varied width of parameters range, it allows examining the utility of using very high or very low parameter values compared to the most commonly used order of magnitude. For example, the common value of b AF is under 200 in the most studies [START_REF] Ayturk | Modeling degenerative disk disease in the lumbar spine: a combined experimental, constitutive, and computational approach[END_REF][START_REF] Jacobs | Validation and application of an intervertebral disc finite element model utilizing independently constructed tissue-level constitutive formulations that are nonlinear, anisotropic, and time-dependent[END_REF][START_REF] Chetoui | Assessment of intervertebral disc degeneration-related properties using finite element models based on ρ h -weighted MRI data[END_REF][START_REF] Honegger | Development of a multiscale model of the human lumbar spine for investigation of tissue loads in people with and without a transtibial amputation during sit-to-stand[END_REF] but it takes a value of 300 in [START_REF] Castro | Intervertebral disc creep behavior assessment through an open source finite element solver[END_REF] and 1045 in [START_REF] Mengoni | Derivation of inter-lamellar behaviour of the intervertebral disc annulus[END_REF]. The high dependency of our model on the value of b AF can be explained by the use of those two extreme values.

Our study outlines the out most influence of changing the stiffness of the extracellular matrix of AF on the poro-hyperelastic behavior of the cervical disc, by taking into account creep compression, physiological relaxation movements during flexion, bending and axial rotation.

It is important to cite that we have faced numerical processing issues when using a low value of µ AF and we suggest using µ AF ≥ 0.8M P a for the same type of model. In the current study, we have found that even the AF compressibility parameters λ AF had a relevant role in the determination of many output results. This influence highlights the necessity of using a compressible extracellular matrix in FE models. These results can explain why many studies using a nearly incompressible model for the IVD are validated with a low bulk modulus which does not satisfy the near incompressibility assumption for the AF (1.45 MPa in [START_REF] Schmidt | Is the ovine intervertebral disc a small human one? A finite element model study[END_REF], 0.8 to 1.4

MPa in [START_REF] Chetoui | Assessment of intervertebral disc degeneration-related properties using finite element models based on ρ h -weighted MRI data[END_REF] and 0.67 MPa in [START_REF] Galbusera | Comparison of four methods to simulate swelling in poroelastic finite element models of intervertebral discs[END_REF]). Moreover, very few studies take into account transverse deformation to the traction direction, which are crucial to estimate the compressible behavior of AF [START_REF] Acaroglu | Degeneration and aging affect the tensile behavior of human lumbar annulus fibrosus[END_REF][START_REF] Elliott | Anisotropic and inhomogeneous tensile behavior of the human annulus fibrosus: Experimental measurement and material model predictions[END_REF][START_REF] Baldit | Experimental analysis of the transverse mechanical behaviour of annulus fibrosus tissue[END_REF][START_REF] Dusfour | Mesenchymal stem cells-derived cartilage micropellets: A relevant in vitro model for biomechanical and mechanobiological studies of cartilage growth[END_REF][START_REF] Derrouiche | The two poisson's ratios in annulus fibrosus: Relation with the osmo-inelastic features[END_REF]. The last two experimental studies have clearly outlined, on bovine and porcine disks respectively, the auxetic behavior of AF, and both report the dependency of the apparent Poisson's ratios to the deformation amplitude, with negative Poisson's ratios in the circumferential radial plane. Therefore, on the compressible behavior of the AF, we conclude that more experimental, numerical, and theoretical studies are needed to understand the complex, non-linear auxetic behavior of the AF. On the other hand, the relevance of the λ N P has not been shown in the current study. We think that the role of this parameter was hidden by its small range of study (0.1 to 0.5M P a). However, the compressible behavior of the human NP has been clearly outlined in several experimental studies [START_REF] Cloyd | Material properties in unconfined compression of human nucleus pulposus, injectable hyaluronic acid-based hydrogels and tissue engineering scaffolds[END_REF][START_REF] Johannessen | Effects of degeneration on the biphasic material properties of human nucleus pulposus in confined compression[END_REF]. This incites us to suggest using a compressible swelling material for the NP to catch its mechanical behavior except at very high strain rates or under dynamic conditions where a viscoelastic incompressible behavior could be sufficient [START_REF] Iatridis | Is the nucleus pulposus a solid or a fluid? mechanical behaviors of the nucleus pulposus of the human intervertebral disc[END_REF][START_REF] Iatridis | The viscoelastic behavior of the nondegenerate human lumbar nucleus pulposus in shear[END_REF][START_REF] Leahy | Viscoelastic properties of the nucleus pulposus of the intervertebral disk in compression[END_REF].

As we have shown, the permeabilities of the AF and the NP do not contribute to the determination of the equilibrium state. However, the radial flow rates are sensitive to AF permeability and less sensitive to NP permeability. The latter contributes more to the determination of the pore pressure. It is known that the water exchange within the IVD is a key factor for nutrient transport and cell activity. Therefore, inaccurate estimation of the permeabilities, especially of the AF, may lead to a loss of sufficiency in the biomechanical models of treating degeneration and growth. The current strain-dependent permeability laws

proposed by [START_REF] Argoubi | Poroelastic creep response analysis of a lumbar motion segment in compression[END_REF] and [START_REF] Gu | New insight into deformation-dependent hydraulic permeability of gels and cartilage, and dynamic behavior of agarose gels in confined compression[END_REF] are sufficient to investigate the IVD solid-like behavior as reported by many studies. However, given their isotropic formulation, they are not able to describe precisely the water exchange rate in the fiber-reinforced tissues like AF [START_REF] Chetoui | Assessment of intervertebral disc degeneration-related properties using finite element models based on ρ h -weighted MRI data[END_REF]. We anticipate that a more accurate permeability expression that depends on strain but also on collagen fiber network organization will have a higher impact on water flow rates.

Another important finding is that collagen fiber properties including rigidity, nonlinearity, and orientation have no significant effect on the response of the IVD to compression. However, they are a determinant factor in the cases of flexion, lateral bending, and axial rotation.

As mentioned previously, the high dependency of IVD behavior to the value of b AF may be partially caused by the large difference between its higher and lower tested values. The fiber nonlinearity parameter is not a measurable property which makes it difficult to estimate.

Therefore, on b AF value, we suggest performing more oriented AF local identification studies to limit the range of variation of this parameter. Collagen fibers have high stiffness but they do not contribute to the mechanical behavior only when they are in tension. Contrary to compression, in the case of rotational motion, the strain principal directions are near to the fiber orientation angles due to the large torsion deformation, especially in axial rotation. This result comes back to the results of [START_REF] Argoubi | Poroelastic creep response analysis of a lumbar motion segment in compression[END_REF][START_REF] Cortes | Elastic, permeability and swelling properties of human intervertebral disc tissues: A benchmark for tissue engineering[END_REF] and [START_REF] Chetoui | Prediction of intervertebral disc mechanical response to axial load using isotropic and fiber reinforced FE models[END_REF] affirming that isotropic models are able to reproduce the IVD response to compression. However, as underlined by [START_REF] Yang | Effect of collagen fiber orientation on intervertebral disc torsion mechanics[END_REF], it is essential to consider fiber reinforcement when studying the IVD behavior in axial rotation. Furthermore, collagen network architecture should be neatly described by namely regional orientation difference between inner/outer and posterior/anterior AF.

This study also revealed that the collagen network within the AF has a crucial role to optimize the necessary resulting moment to turn the head. In combination with the relative movements of the adjacent vertebrae during axial rotation, lateral bending, and flexion, the mechanical properties of the collagen fiber network appear to be particularly well adapted to reduce the resulting moment by an order of magnitude compared to pure rotations. This key feature has to be taken into account while designing IVD implants.

To take back the main aim of this study, we suggest simplifying IVD model validation studies by reducing the number of identified parameters depending on the experimental optimization objective and load type. [START_REF] Zhuang | Construction of tissue-engineered composite intervertebral disc and preliminary morphological and biochemical evaluation[END_REF] For tests in statics, very slow loading tests, and when the experimental object is not time-dependent such as force-displacement curve, there is no need to identify the permeability value. (2) For fiber-reinforced models, fiber orientation and nonlinearity should not be included in the identified parameter except when studying IVD response to rotational motion. Fiber rigidity is a determinant parameter for all types of load.

(3) None of the mechanical parameters of the extracellular matrix of the AF could be neglected in the identification process. High attention should be paid to the determination of the stiffness of the AF ground substance.

Some sources of uncertainty could be identified from the current results. On a numerical level, the results of this study are provided without focusing on the uncertainty quantification of each parameter in the validation step. However, we can identify some parameters with infinite uncertainty such as the two permeabilities for stationary analysis or the fiber orientation for the uniaxial load. In the identification studies, the accuracy of the model needs to be investigated by uncertainty quantification. This will be useful to localize the optimal experimental fields to measure, design accurate experiments and therefore minimize model error [START_REF] Sutula | An open source pipeline for design of experiments for hyperelastic models of the skin with applications to keloids[END_REF][START_REF] Elouneg | An open-source FEniCS-based framework for hyperelastic parameter estimation from noisy full-field data: Application to heterogeneous soft tissues[END_REF][START_REF] Mazier | Inverse deformation analysis: an experimental and numerical assessment using the FEniCS Project[END_REF]. On a physiological level, only the IVD prestress related to the osmotic pressure was included. However, a second part, related to the gradient of tissue growth could be taken into account. The IVD prestress was studied for bovine and porcine tissues [START_REF] Duclos | Residual strains in the intervertebral disc annulus fibrosus suggest complex tissue remodeling in response to in-vivo loading[END_REF]76,77,[START_REF] Dusfour | Mesenchymal stem cells-derived cartilage micropellets: A relevant in vitro model for biomechanical and mechanobiological studies of cartilage growth[END_REF] but no data concerning human disc prestress are yet provided in the literature to our knowledge.

In addition, the loading conditions, taken from previous experimental studies, may affect the current results given the non linear anisotropic behavior of the model.

Our study, which focused on varying each parameter independently of the others, revealed the complex mechanical behavior in in-vivo solicitations. In future work, the variation of multiple parameters at the same time will be performed. We anticipate this scheme may reveal coupled effects and allow more efficient simplifications of the validation studies. The head-weight force in compression and the coupled displacements in the other physiological loads were applied on the upper surface of the disc which is assumed to be rigid.

The apparent stiffness is calculated using the applied load or displacement and the resulting displacement or external force, respectively. For compression, the apparent stiffness

A pp S(N m -1
) is defined as:

A pp S = F app w res
where F app (N ) is the applied force and w res (m) is the displacement of the upper surface in the z direction (the same direction as the applied force).

For the other physiological loads, the apparent stiffness A pp S(N m rad -1 ) is defined as:

A pp S = ⃗ M res .⃗ e p φ app
where ⃗ M res (N m) is the global moment of the upper surface points about the center of rotation, φ app (rad) is the global cervical segment rotation as described in the table 1 and ⃗ e p the principal direction of the rotation: ⃗ x for flexion, ⃗ y for lateral bending and ⃗ z for axial rotation. M res (N m) is obtained according to this equation: λ CEP (M P a) 10.44 - C f c0 (mol/m 3 ) 90 - [START_REF] Chetoui | Assessment of intervertebral disc degeneration-related properties using finite element models based on ρ h -weighted MRI data[END_REF] ϕ 0 0.8 -

⃗ M res = up ⃗ M n,res dS = up ⃗ OM n ∧ ⃗ F n,res dS

Figure 6

 6 Figure 6 illustrates the influence of modifying model parameters on the IVD response to

Figure 7

 7 Figure7shows the influence of modifying model parameters on the IVD response respec-

1 T 9 φ e External osmotic coefficient=0. 9 385

 199 Director vector of the i th fiber family p Fluid phase hydrostatic pressure (P a) ∆π Osmotic pressure (P a) k 0 initial permeability m 4 N -1 s -1 C f c Fixed charge density mol m -3 C f c0 Initial fixed charge density mol m -3 Fiber non linearity coefficient (unitless) K Fiber dispersion coefficient (unitlessI f iDirection invariant of the i th family fiberα i Inner fiber orientation ( • ) α o Outer fiber orientation ( • ) ∆α (α i -α o )/2 ( • ) Constants R universal gaz constant =8.314J mol -1 .K -Absolute temperature =310KC e External salt concentration =150mol m -3 φ i Internal osmotic coefficient =0.Appendix: Apparent stiffness 386

Figure 1 :

 1 Figure 1: Geometry of C6-C7 IVD. A. Sagittal section. B. Transverse section. C. Exploded view.

Figure 2 :Figure 4 :Figure 5 :

 245 Figure 2: Definition of fiber orientation. A. Thermal solution and definition of the new coordinate system. B. Fibers directions in the inner and the outer lamellae of the AF

Figure 6 :

 6 Figure 6: Impact of model parameters on IVD response to lateral bending. The bar charts: difference between the 1 and the 4 cases (the 2 and 3 for the hatched area) normalized to the reference result (ref) obtained with the basic model. The line graph corresponds to the evolution of the apparent stiffness in lateral bending with the different scenarios of 7 parameters.

Figure 7 :

 7 Figure 7: Impact of model parameters on IVD response to flexion. The bar charts: difference between the 1 and the 4 cases (the 2 and 3 for the hatched area) normalized to the reference result (ref) obtained with the basic model. The line graph corresponds to the evolution of the apparent stiffness in flexion with the different scenarios of 7 parameters.

Figure 8 :

 8 Figure 8: Impact of model parameters on IVD response to axial rotation. The bar charts: difference between the 1 and the 4 cases (the 2 and 3 for the hatched area) normalized to the reference result (ref) obtained with the basic model. The line graph corresponds to the evolution of the apparent stiffness in axial rotation with the different scenarios of 7 parameters.

Table 1 :

 1 Applied intervertebral range of motion for the relaxation step for each load type. Motion: physiological motion of the whole cervical segment. Intervertebral rotation: C6-C7 rotation with ILB:lateral bending, IFl: flexion and IAR: axial rotation. Intervertebral translation: C6-C7 rotation with P/A: posterior(-)/interior(+), L/R: left (-)/Right(+) and D/U : down (-)/ up (+). For lateral bending and axial rotation, the intervertebral motion values correspond to 10°and 6°, respectively, of the C3-C7 segment and were taken from the curves in the cited reference. For the flexion, the intervertebral motion values correspond to the half of the studied range of motion (59.5°of the C1-T1 segment) in the cited reference.MotionIntervertebral rotation ( o ) Intervertebral translation (mm) Ref ILB (X) IFl (Y) IAR(Z) P/A (X) L/R (Y) D/U (

Table 2 :

 2 Model parameters values for the basic model and the cases used for the parametric study. Variation of each parameter of this list was separately tested.

		Parameter	Basic Cases				References used in the definition
				1	2	3	4	of the intervals of study
	NP	µ N P (M P a)	0.4	0.1 0.25 0.55 0.7 [58, 53]
		λ N P (M P a)	0.3	0.1 0.2	0.4	0.5 [58, 53]
		k0 ×10 -4 N P C f c0 (mol/m 3 ) (mm 4 /(N s)) 13 300	3 -	8	18	23	[52, 57] [58]
		ϕ 0	0.8	-				[58]
	AF ground µ AF (M P a)	1.2	0.4 0.8	1.6	2	[7, 57]
	substance	λ AF (M P a)	1.3	0.2 0.75 1.85 2.4 [70, 4]
		k0 ×10 -4 AF C f c0 (mol/m 3 ) (mm 4 /(N s)) 13 150	3 -	8	18	23	[52, 57] [11]
		ϕ 0	0.7	-				[58]
	AF fibrillar a AF (M P a)	1.2	0.3 0.75 1.65 2.1 [52, 54]
	network	b AF	350	50 200 500 650 [11, 8]
		∆α	5	0	2.5	7.5	10	[52, 71]
		K AF	0.166 -
	CEP	µ CEP (M P a)	7.14	-				[70]

Source of funding

This work has been supported by MESO@LR Platform of the University of Montpellier, by Labex NUMEV (ANR-10-LABX-20 projects 2017-1-27) and by CNRS (AAP "Osez l'interdesciplinarité 2018", MoTiV Project)

where ⃗ M n,res (N m) is the resulting nodal moment, ⃗ F n,res (N ) the resulting nodal force, and ⃗ OM n (m) the final nodal position with regard to the center of rotation.

Conflict of interests

The authors declare that they have no conflict of interests to disclose.

Ethical approval

Not required

Annulus fibrosus functional extrafibrillar and fibrous mechanical behaviour: experimental and computational characterisation. R Soc Open Sci 2017;4(8):170807

[77] Michalek AJ, Gardner-Morse MG, Iatridis JC. Large residual strains are present in the intervertebral disc annulus fibrosus in the unloaded state. J Biomech 2012;45 [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]:1227-31