Madame Le Professeur Faiza

Keywords: Multi-layered soft tissues, Microstructure, Volumetric strain, Constitutive model, Time-dependent response., Transpose  Simple contraction Lamellar ILM: Inter-lamellar ECM: Extracellular matrix OCF: Oriented collagen fibers PG: Proteoglycan I.1 Full disc model, Multiscale structure, Regional variation, Nucleus-annulus interaction

Elle a été soutenue financièrement par le programme de bourses PROFAS B+ et le programme de bourses de mobilité internationale de l'Université de Lille (MOBILEX). Je profite de l'occasion pour remercier les deux programmes pour le support financier. Je souhaite, en premier lieu, exprimer toute ma gratitude et mes remerciements à mes directeurs, le Professeur Fahmi ZAÏRI et le Professeur Amar MESBAH pour m'avoir encadré, guidé et encouragé tout au long de mes travaux. J'ai beaucoup appris à leurs côtés, que ce soit au niveau scientifique qu'au niveau humain. * Merci Fahmi, d'avoir toujours montré un grand intérêt pour mes travaux. Merci de m'avoir transmis votre passion pour la recherche et votre rigueur de travail qui m'ont permis de mener à bien ce projet. * Merci Amar pour votre disponibilité et votre grande réactivité face à tous mes problèmes. Un remerciement tout particulier à l'ensemble des membres de mon jury pour avoir accepté de contribuer à l'évaluation de mes travaux de recherche et de participer à ma soutenance de thèse. Je remercie Monsieur le Professeur Jean-François GANGHOFFER pour l'honneur qu'il me fait d'accepter la présidence de ce jury. Mes vifs remerciements sont adressés à Monsieur le Professeur Grégory CHAGNON et Monsieur le Professeur Saïd RECHAK pour avoir accepté de juger ce travail en tant que rapporteurs et pour leurs rapports pertinents

General introduction

Back pain is a common health problem affecting about 50 to 80% of humans at least one time during their life. The pain is sometimes related to muscles and ligaments, but when it is intense and chronic it is often linked to the intervertebral disc injuries and dysfunctions. This soft tissue located in the vertebral column gives to the body the mobility needed to achieve the different movements while resisting in parallel to the loads coming from the body weight and the daily activities. Any small issue in the intervertebral discs may cause many troubles in the human back varying from a simple back pain to disc hernia and sometimes could lead to paralysis. For this reason, the importance of studying the intervertebral disc conduct and predicting its response under the different physiological movements is mandatory in order to understand the chemo-mechano-biological mechanisms taking place inside the disc soft tissues which could help avoiding damage and producing better treatment procedures if damage is already present. Great progresses in this field were achieved over the recent years. Although important differences could be witnessed between the discs of different individuals, many experimental contributions succeeded to identify the mechanical response of the disc at different scales: complete spine scale, functional spine unit scale and material volume element scale.

However, accessing the core of the disc in-vivo while maintaining unchanged its natural response is very hard to achieve which represents an obstructing point in the intervertebral disc research. That is why multi-physics simulations relating the microstructure of the chemo-biological tissue to its multiaxial response become mandatory for understanding the disc behavior under normal activity.

However, the multi-subject problematic translated by the complex interactions between the different microstructural elements of the disc makes the task very complex. Also, due to the high heterogeneity of the biological tissues intervening in the intervertebral disc and the varying response between different individuals, a huge amount of experimental works of hardly available human cadavers should be established in order to construct and validate a reliable model of the disc.

Objective and thesis plan

The damage in annulus fibrosus soft tissues is a complex multiscale phenomenon due to a complex structural arrangement of collagen network at different scales of hierarchical organization. A fully three-dimensional constitutive representation that considers the regional variation of the structural complexity to estimate annulus multiaxial mechanics till failure has not yet been developed. In the present PhD dissertation, a model, formulated within the framework of nonlinear continuum mechanics, is developed to predict deformation-induced damage and failure of annulus under multiaxial loading histories considering as time-dependent physical process both chemical-induced volumetric effects and damage accumulation. In order to realize our objective, the actual thesis was planned following several well-constructed steps that will be described in details through the following chapters:

In Chapter I, a microstructure-based model is proposed to connect structural features, intrinsic mechanics and electro-chemical properties of annulus soft tissues. The multi-layered lamellar/interlamellar annulus model is constructed by considering the effective interactions between adjacent layers and the chemical-induced volumetric strain. The model/experiments comparison demonstrates that the evaluation of the overall time-dependent response involves considering stress, volumetric change and auxetic feature simultaneously in relation to structural features.

In Chapter II, the model is enriched by considering the hierarchical structure of the soft tissue from the nano-sized collagen fibrils to the micro-sized oriented collagen fibers. The stochastic process of progressive damage events operating at different scales of the solid phase is introduced for the extracellular matrix and the network of nano-sized fibrils/micro-sized fibers. The directional effects on annulus mechanics and failure are highlighted in relation to external loading mode, structure features, damage events and hydration.

In Chapter III, the model is further developed by considering the regional variation of the complex structural organization of collagen network at different scales to predict the regional anisotropic General introduction 3 multiaxial damage of the intervertebral disc. After model identification using single lamellae extracted from different disc regions, the model predictability is verified for various multiaxial elementary loading modes representative of the spine movement. The stretching along the circumferential and radial directions till failure serves to check the predictive capacities of the annulus model for the different regions. Model results under simple shear, biaxial stretching and plane-strain compression are further presented and discussed.

In Chapter IV, a full human disc model is constructed using the regional annulus model to examine the heterogeneous mechanics in the disc core. Damage fields in the disc are analyzed under axial compression, axial twist and combined loadings to assess the areas where the risk of failure is the The biological soft tissues are electro-chemically active media and many of them are fiberreinforced laminates. Their constitutive modeling including structure-property relationship is of prime importance for a better understanding of the tissue biomechanics. It may also help in a broader way to better understand injury mechanisms, to detect and monitor pathologies and to design new bio-substitutes.

The formulation of a constitutive model has to consider three main key points in order to physically describe the underlying chemo-mechanical mechanisms and the macro-response: (i) firstly, a microstructure representation has to be proposed, (ii) secondly, free energy functions have to be defined for each constituent of the microstructure, (iii) thirdly, the electro-chemical activity resulting in the fluid flow inside the tissue has to be taken into account. In existing constitutive models, the structural features of the soft tissue are basically related to the intrinsic properties of extracellular matrix (ECM) and collagen fibers (Holzapfel et al., 2000;[START_REF] Peng | An anisotropic hyperelastic constitutive model with fiber-matrix shear interaction for the human annulus fibrosus[END_REF][START_REF] Balzani | A polyconvex framework for soft biological tissues: adjustment to experimental data[END_REF] Chapter I: A microstructure-based model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus 7 [START_REF] Guo | A composites-based hyperelastic constitutive model for soft tissue with application to the human annulus fibrosus[END_REF][START_REF] Kroon | A model for saccular cerebral aneurysm growth by collagen fibre remodelling[END_REF][START_REF] Li | A structural multi-mechanism constitutive equation for cerebral arterial tissue[END_REF][START_REF] Chen | A micromechanics finite-strain constitutive model of fibrous tissue[END_REF]. The presence of the fluid phase is implicitly modeled by considering a volumetric contribution into the free energy function (Rodriguez et al., 2006;[START_REF] Pierce | A hyperelastic biphasic fibre-reinforced model of articular cartilage considering distributed collagen fibre orientations: continuum basis, computational aspects and applications[END_REF][START_REF] Swedberg | Continuum description of the Poisson's ratio of ligament and tendon under finite deformation[END_REF][START_REF] Mcevoy | Compressibility and anisotropy of the ventricular myocardium: experimental analysis and microstructural modeling[END_REF]Derrouiche et al., 2019b;Kandil et al., 2019). In other constitutive models based on bi-or tri-phasic theory, the osmotic swelling effects due to the electro-chemical activity in the soft tissue are more explicitly taken into account by considering the interactions between charged ECM and fluid phase [START_REF] Mow | Biphasic creep and stress relaxation of articular cartilage in compression: theory and experiments[END_REF][START_REF] Holmes | The nonlinear characteristics of soft gels and hydrated connective tissues in ultrafiltration[END_REF][START_REF] Lai | A triphasic theory for the swelling and deformation behaviors of articular cartilage[END_REF][START_REF] Iatridis | Influence of fixed charge density magnitude and distribution on the intervertebral disc: applications of a poroelastic and chemical electric (PEACE) model[END_REF][START_REF] Yao | Three-dimensional inhomogeneous triphasic finite-element analysis of physical signals and solute transport in human intervertebral disc under axial compression[END_REF][START_REF] Ehlers | An extended biphasic model for charged hydrated tissues with application to the intervertebral disc[END_REF][START_REF] Stalhand | A mechanochemical 3D continuum model for smooth muscle contraction under finite strains[END_REF][START_REF] Stracuzzi | Chemomechanical models for soft tissues based on the reconciliation of porous media and swelling polymer theories[END_REF].

This work is focused on a highly complex multi-layered soft tissue, that is, the annulus fibrosus of the intervertebral disc. The annulus fibrosus is a fiber-reinforced tissue in which the ECM is reinforced by oriented type-I collagen fibers (OCF) providing tensile strength and mechanical stiffness (Holzapfel et al., 2005;[START_REF] Roberts | Histology and pathology of the human intervertebral disc[END_REF]. To these "solid" components, the soft tissue contains also a fluid phase (composed of water, mobile charges and small proteins) that interacts with the fixed charges of the ECM providing the chemo-mechanical coupling (Derrouiche et al., 2019a(Derrouiche et al., , 2020b)). The successive fiber-reinforced lamellar (LM) layers, exhibiting alternate fiber angles, are inter-separated by an interlamellar (ILM) matrix, i.e. intermediate ground substance devoid of OCF (Pezowicz et al., 2006). The organization is such that the thickness of layers is graded and increases from the outer annulus towards the inner annulus (Cassidy et al., 1989;Marchand and Ahmed, 1990;[START_REF] Hsu | Diffusion tensor microscopy of the intervertebral disc anulus fibrosus[END_REF]Holzapfel et al., 2005). It is only in very recent experimental studies that the ILM zone has been appreciated from structural/mechanical viewpoint [START_REF] Michalek | Measurement of local strains in intervertebral disc anulus fibrosus tissue under dynamic shear: contributions of matrix fiber orientation and elastin content[END_REF]Vergari et al., 2016;Tavakoli et al., 2016Tavakoli et al., , 2017Tavakoli et al., , 2018)). It is considered as the main factor ensuring disc structural integrity and disc shearing resistance. To date, the modeling studies including the ILM zone are rare [START_REF] Nerurkar | Modeling interlamellar interactions in angle-ply biologic laminates for annulus fibrosus tissue engineering[END_REF]Labus et al., 2014;[START_REF] Adam | Inter-lamellar shear resistance confers compressive stiffness in the intervertebral disc: an image-based modelling study on the bovine caudal disc[END_REF][START_REF] Mengoni | Derivation of inter-lamellar behaviour of the intervertebral disc annulus[END_REF]Derrouiche et al., 2019b;Kandil et al., 2019). This ground substance is introduced solely into numerical models of the disc and its role is generally limited to a sliding zone Chapter I: A microstructure-based model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus 8 with either cohesive or non-cohesive features. Recently, Derrouiche et al. (2019b) and Kandil et al. (2019) proposed a numerical model of the annulus considering the LM/ILM interaction effects on the intrinsic properties and the fluid transfer through the multi-layered soft tissue in the aim to consider the shrinkage and swelling mechanisms. The latter may be considered as indicators of the chemo-mechanical coupling in annulus during the annulus stretching2 . This coupling in annulus is believed to be the main responsible of an unusual time-dependent transversal behavior with out-ofbounds Poisson's ratios (higher than 0.5 and smaller than 0 for isotropic media) during short-term stretching followed by equilibrium return during long-term interrupted stretching (relaxation)

towards more usual values [START_REF] Baldit | Experimental analysis of the transverse mechanical behaviour of annulus fibrosus tissue[END_REF]Derrouiche et al., 2019bDerrouiche et al., , 2020a;;Kandil et al., 2019).

To date, no analytical constitutive modeling is proposed to describe these effects in annulus.

The purpose of this chapter is to formulate a constitutive model for the annulus time-dependent response as the result of the coupling between mechanics and electro-chemical activity. We propose a two-stage strategy consisting in a decoupling stage of each layer of the multi-layered soft tissue using full-field strain data and a re-coupling stage reconnecting the layers by inter-layer ionic diffusion as well as stress/strain continuity along the interfaces. The final model is then utilized to examine the effective role of the relevant structural features (ILM zone and LM/ILM dimensions)

and the diffusion-mechanics on the time-dependent response of the multi-layered soft tissue. The response is examined under three successive steps: free swelling, quasi-static stretching and relaxation.

The outline of the present chapter is as follows. Section I.2 is devoted to the microstructure-based model formulation. The hybrid experimental/modeling decoupling strategy is described in Section I.3. The effective contribution of the LM and ILM layers on the coupling is discussed in Section I.4.

Concluding remarks are given in Section I.5.

The following notation is used throughout the text. Tensors and vectors are denoted by normal boldfaced letters and italicized boldfaced letters, respectively, while scalars and individual components of vectors and tensors are denoted by normal italicized letters. Simple contraction of two vectors or two tensors is denoted by a dot ''  '', and a direct (outer) product by the symbol  . The LM layer is regarded as a three-phase medium consisting in ECM, OCF and fluid phase whereas the ILM layer is a two-phase medium consisting in ECM and fluid phase3 . The deformation gradient F of the multi-layered soft tissue is given by: 1 __ _ 0 _ 0 1,3... 2,4... is the initial thickness of the layer i.

The condition of continuity of the deformation along the interfaces writes:

1 00

i i i i + = F .
F . nn (I.3) where 0 i n (i=1, …m-1) is the arbitrary unit vector between two adjacent layers as illustrated in Figure I.2a.

To combine chemical deformation and mechanical deformation in each layer, the conceptual sequence of configurations is used by the introduction of an intermediate chemical configuration Similarly, the invariants of the chemical part

T chem chem chem = C
F .F can be defined as follows: where 0 a is again the unit vector of the fibers direction. Note that a sixth invariant is introduced in Eq. (I.11) according to the derivations in Appendix I.A for the volumetric change in a fiberreinforced material system.

( ) ( ) ( ) ( ) ( )

I.2.2. Free energy: microstructure and chemo-mechanics decomposition

According to the chemo-mechanics decomposition, the free energy W is decomposed into two parts:

mech chem W W W =+ (I.12)
The interaction between solid and fluid phases is considered via the concept of volume fraction. For each layer, the respective (time-dependent) volume fraction of the main constituents is expressed as a function of its chemically-induced volumetric change i chem J :

Chapter I: A microstructure-based model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus

14 _0 _0 _0 1 1 i ECM i ECM i chem i OCF i OCF i chem i fluid i fluid i chem J J J       = = - =- (I.
    + + = = + = = - (I.14)
in which m is the total number of layers.

The mechanical part i mech W is written as follows for the LM and ILM layers, respectively: ( )

( ) ( ) ( ) _ _ _ __ 1,3,... 2, 4,... 1 LM i i i LM i i i LM i mech ECM ECM mech OCF OCF mech ILM i i i ILM i mech ECM ECM mech W W W i m W W i m   = + = = =
1 1 3 2 ii ECM ECM W G I =- (I.16) in which ECM G
is the ECM shear modulus.

In addition to fibers direction, the transition from initially undulated to fully stretched states of the OCF (Pezowicz et al., 2006) affects their intrinsic nonlinear response. Thus, the OCF free energy i OCF W takes the following exponential form which depends only on the fourth-invariant (Holzapfel et al., 2000): (  )

( ) ( ) ( ) 1 exp 1 1 2 ii OCF C W C I C = - - (
2 1 2 2 1 exp 1 1 2 ii fluid chem K W K J K = - - (I.19)
where 1 K and 2 K are material constants and, their ratio 12 KK represents the volumetric stiffness.

The volumetric change i chem J takes different forms for the LM and ILM layers:

6 1 5 2 4 3 _ _ 1,3,... 2, 4,... 1 L i i i i LM i chem IM i chem i i J I I I I I I im J i m = == -+ = - (I.20)
Further details are provided in Appendix I.A.

The free energy W of the multi-layered soft tissue is the sum of free energies of the LM layers

_ LM i i i mech chem W W W =+
and free energies of the ILM layers Note that, in the intervertebral disc, the proportions of individual layers vary from the outer annulus towards the inner annulus (Cassidy et al., 1989;Marchand and Ahmed, 1990;[START_REF] Hsu | Diffusion tensor microscopy of the intervertebral disc anulus fibrosus[END_REF]Holzapfel et al., 2005). In our approach, the thickness of layers varies continuously following a Chapter I: A microstructure-based model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus 16 special function in the aim to describe a gradual variation of fractions through the multi-layered soft tissue (see Appendix I.B for more details).

_ ILM i i i mech chem W W W =+: ( ) ( )

I.2.3. Constitutive stress-strain response

According to Eq. (I.21), the Cauchy stress σ in the multi-layered soft tissue is expressed as: are the Cauchy stresses in the LM and ILM layers.

( ) ( )
In order to respect the compatibility along the interfaces, the conditions of continuity of the stresses is:

1 i i i i + = σ . σ . nn (I.24)
where i n (i=1, …m-1) is the arbitrary unit vector between two adjacent layers in the current configuration  .

The mechanical and chemical contributions in the stress may be additively split in reason of the chemo-mechanical coupling:

( ) ( ) 1 _ _ _ _ 1,3,.. 2,4,.. mm i LM i LM i i ILM i ILM i layer mech chem layer mech chem ii  - == =+  σ σ σ σ σ ++ (I.25)
Each part of the Cauchy stress is obtained from the differentiation of the free energy functions with respect to the corresponding deformations:

2 T i i i i i mech mech mech mech i mech W p  =-  σ F . .F I C (I.26) i i chem chem i chem W J  =  σI (I.27)
where i p represents a hydrostatic pressure in each layer introduced in reason of the mechanical incompressibility and determined from boundary conditions. The term I represents the identity tensor.

Chapter I: A microstructure-based model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus 17

In Eq. (I.26), the derivation may be expressed as follows:

5 1 i i i mech mech i i i mech mech W W I I    =     =       CC (I.28) with ( ) 1 2 1 1 3 3 4 00 5 0 0 0 0 i i mech i ii mech i mech i ii mech i mech i i mech i ii mech mech i mech I I I I I I I -  =   =-   =   =   =  +   I C IC C C C C C . .C C aa a a a a (I.29)
After a series of straightforward derivations, the different parts of the Cauchy stress are given by:

( ) ( ) ( ) ( ) ( ) __ 2 _ _ _ _ 1 4 2 4 0 0 2 1 exp 1 T LM i i LM i i mech ECM ECM mech i LM i LM i LM i LM i OCF mech mech Gp C I C I   =- + - -  σ B I F . .F aa (I.30) __ ILM i i ILM i i mech ECM ECM mech Gp  =- σ B I (I.31) ( ) ( ) ( ) 2 12 1 exp 1 i i i i chem fluid chem chem K J K J  = - - σI (I.32)
where

T i i i mech mech mech = B F .F
is the mechanical part of the left Cauchy-Green deformation.

I.2.4. Electro-chemical diffusivity kinetics

Due to the physiological environment and the intrinsic ionic strength of the ionized ECM, the biological soft tissues are able to retain water up to several times their dehydrated volumetric state.

As the fixed charges density of the charged hydrated ECM is large, a significant electro-chemical flux occurs through the multi-layered soft tissues leading to their swelling [START_REF] Frijns | A validation of the quadriphasic mixture theory for intervertebral disc tissue[END_REF].

Because of swelling, an osmotic pressure is generated [START_REF] Urban | Diffusion of small solutes into the intervertebral disc: as in vivo study[END_REF]. According to the The phenomenon of ionic migration of the solute in the charged hydrated soft tissue is related to absorption and diffusion process as illustrated in Figure I.2b. The motion process of the ions may be classically expressed using the second Fick law depending on time and space [START_REF] Ferguson | Fluid flow and convective transport of solutes within the intervertebral disc[END_REF][START_REF] Soukane | Investigation of solute concentrations in a 3D model of intervertebral disc[END_REF]:

( ) ( ) ( ) ( ) , , , , , , c x y z t c x y z t t    =    D (I.35)
where  D is the diffusivity of the ions.

Diffusivity in biological soft tissues is anisotropic due to material heterogeneity [START_REF] Hsu | Diffusion tensor microscopy of the intervertebral disc anulus fibrosus[END_REF][START_REF] Antoniou | Apparent diffusion coefficient of intervertebral discs related to matrix composition and integrity[END_REF] and depends on the strain and the principal direction of diffusion [START_REF] Jackson | Anisotropic ion diffusivity in intervertebral disc: an electrical conductivity approach[END_REF][START_REF] Jackson | Effect of compression and anisotropy on the diffusion of glucose in annulus fibrosus[END_REF][START_REF] Jackson | A conductivity approach to measuring fixed charge density in intervertebral disc tissue[END_REF][START_REF] Travascio | Anisotropic diffusive transport in annulus fibrosus: experimental determination of the diffusion tensor by FRAP technique[END_REF]. The tensor of the diffusivity takes the following form: The diffusivity of the ions is expressed as follows [START_REF] Gu | Diffusivity of ions in agarose gels and intervertebral disc: effect of porosity[END_REF]:

( ) ( ) 12 exp free s Dr      - =- Dk (I.37)
where  and  are two coefficients, free D  is the free diffusivity of the spherical solutes in the external solution, s r  is the solute size and k is the permeability tensor given by [START_REF] Gu | Diffusivity of ions in agarose gels and intervertebral disc: effect of porosity[END_REF][START_REF] Perie | Correlating material properties with tissue composition in enzymatically digested bovine annulus fibrosus and nucleus pulposus tissue[END_REF]:

( ) _0 0 _0 1 M chem fluid chem fluid J J   -  - =   -  kk (I.38)
in which M is a coefficient and 0 k is the zero-strain permeability tensor in the x, y and z directions, and:

( ) What may seem like a reconciliation of the two theories is especially valid when the ionic equilibrium is achieved (chemical potential of the fluid is null). Indeed, the chemical disequilibrium generated by a stretching and the equilibrium return when the stretching is interrupted induce differentiations between layers swelling as illustrated in Figure I.2b.

0 _ 0 _ 0 _ 0 diag , ,

I.3. Hybrid experimental/modeling decoupling strategy

I.3.1. Database

The model is applied on annulus specimens extracted from bovine cervical spine discs using two useful macro-quantities to design the material kinetics and to identify the model parameters: the volumetric change and the macro-stress. The experimental protocol is detailed elsewhere (Derrouiche et al., 2020a) but some important information is provided succinctly here. • The first step is a free swelling: The local strains were measured until the balance of swelling.

• The second step is a quasi-static stretching: The sample was strained up to a strain of 4% at a strain-rate of 0.001 s -1 in order to study the mechanical behavior coupled with chemistry.

Before starting the second step, the sample was preconditioned by a series of ten loadingunloading cycles at low amplitude of 1% at a strain-rate of 0.001 s -1 in order to stimulate fluid exchanges inside the tissue.

• The third step is a relaxation: At the end of a stretching step, the strain is kept constant until the equilibrium of the response. The auxetic behavior of the lamellae plane is followed during the relaxation by positive values tending towards equilibrium. In the fibers plane, the Poisson's ratio presents values higher than 0.5 and reaches approximately 1.0 at the end of the stretching step. Subsequently, it decreases during the relaxation period to 0.5 at equilibrium.

I.3.2. Deformation of a layer

Under a free swelling, the deformation gradients of a layer i are:

( )

_ _ _ diag , , 1
LM i mech f f f ILM i mech i i i chem x chem x chem z x z im im           -- -- = =       = = -               F R .. F R (I.44)
where  is the applied stretch and, f  is the stretch in the fiber direction.

The quantity i z  R is the basis-change tensor in each layer, with a rotation i  around the z-axis:

cos sin 0 sin cos 0 0 0 1 i ii ii z      =-   R (I.45)
After a series of lengthy but straightforward derivations, we get the mechanical deformation gradient of the LM layers:  is the chemical expansion at the equilibrium state and i j  is the chemical expansion rate. The term t corresponds to the swelling time in Eq. (I.43) and to the stretching and relaxation times in Eq. (I.48). The identification exercise of fitting parameters of Eq. (I.48) can be formulated as an optimization problem and solved numerically using Matlab software [START_REF] Pyrz | Identification of viscoplastic parameters of phenomenological constitutive equations for polymers by deterministic and evolutionary approach[END_REF]

( ) ( ) ( ) ( ) 3 _ 2 1 2 _ _ 3
JJ F J         =        - - -  = + +                      (I.49)
where  is the number of considered data points. The identification result is presented in the top of Figure I.4a for the three successive loading steps. It is worth noticing that the free swelling is anisotropic. The transversal free swelling stretch is higher in the lamellae plane than in the fibers plane. The axial free swelling stretch is found between the two transversal free swelling stretches.

I.3.3. Electro-chemical flux through the multi-layered soft tissue

The electro-chemical flux operates as an inter-layer transfer process of the ions as illustrated in The solution of the Fickian law can be analytically obtained [START_REF] Oldham | The replacement of Fick's laws by a formulation involving semidifferentiation[END_REF][START_REF] Crank | The mathematics of diffusion[END_REF][START_REF] Tixier | Modeling of damage in cement-based materials subjected to external sulfate attack. I: Formulation[END_REF] and it may be expressed for a layer as follows:

Chapter At t =0, the concentration of the electronic charges inside each layer of the tissue is given by the Donnan equation [START_REF] Maroudas | Biophysical chemistry of cartilaginous tissues with special reference to solute and fluid transport[END_REF]:

( ) ( ) ( ) 22 * 00 2 ,0 2 FF ext i c c c cz   + + = (I.54)
where 0 F c is the fixed charge density at the initial configuration.

Table I.1 provides the diffusivity and permeability parameters taken from the literature [START_REF] Gu | Diffusivity of ions in agarose gels and intervertebral disc: effect of porosity[END_REF][START_REF] Perie | Correlating material properties with tissue composition in enzymatically digested bovine annulus fibrosus and nucleus pulposus tissue[END_REF][START_REF] Jackson | A conductivity approach to measuring fixed charge density in intervertebral disc tissue[END_REF] The stress-stretch relationships take the following forms:

( ) ( ) ( ) ( ) ( ) 2 _2 _3 _2 22 1 4 2 4 1 2 1 exp 1 cos LM i i mech xx ECM ECM i chem x i i i i OCF GS S C I C I          = + -      + - - (I.57) 2 _ _ _ _ i chem x ILM i i mech xx ECM ECM i chem x G        =-     (I.58) ( ) ( ) ( ) 2 _ 1 2 1 exp 1 i i i i chem xx fluid chem chem K J K J  = - - (I.59)
The total stress of the multi-layered annulus σ is simply given by Eq. (I.25). The re-coupling is obtained by considering the electro-chemical flux through the multi-layered soft tissue, as constitutively presented in Section I.3.3, while considering stress/strain continuity along the interfaces to account for mechanical interactions between adjacent layers, see Eqs (I.3) and (I.24).

The latter compatibility was also considered by [START_REF] Van Dommelen | Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers[END_REF] for semi-crystalline polymer systems regarded as two-phase (crystalline/amorphous) layered composite inclusions.

Several inputs related to structural features are required by the modeling. The initial volume fraction of fluid is 74% by wet weight [START_REF] Demers | Value and limitations of using the bovine tail as a model for the human lumbar spine[END_REF][START_REF] Perie | Correlating material properties with tissue composition in enzymatically digested bovine annulus fibrosus and nucleus pulposus tissue[END_REF]. The initial volume fraction of OCF is 7% [START_REF] Oshima | The use of coccygeal discs to study intervertebral disc metabolism[END_REF][START_REF] Demers | Value and limitations of using the bovine tail as a model for the human lumbar spine[END_REF]) and the angle with respect the loading axis is 26° for the inner zone (Holzapfel et al., 2005). The thickness of the ILM zone is one eighth that Chapter I: A microstructure-based model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus 31 of the LM layer (Tavakoli et al., 2017) while a gradual variation of thickness is adopted through annulus.

The parameters involved in the stress expression have been identified using the experimental stress where  is again the number of considered data points and mod  is the stress provided by the model. Although the ILM zone is devoid of OCF and represents a small fraction of the whole tissue, Figure I.5b shows that it contributes significantly to the apparent stiffness and the stabilized relaxed due to the ILM swelling inducing a significant chemically-induced volumetric stress in the ILM zone.

Physically, perpendicular fibrils bridging two adjacent LM layers are observed in the ILM zone (Tavakoli et al., 2017). These perpendicular fibrils could play the role of "elastic springs" in order to retain the significant ILM swelling. In the absence of diffusion, the apparent stiffness decreases and the stress no longer relax upon interrupted stretching. The absence of the electro-chemical activity (loss of fixed charges) simulates the dysfunction of whole PG macromolecules, their loss being the main indication of disc degeneration6 . As a consequence of the inter-layer fluid transfer mechanism, the unusual time-dependent transversal behavior is vanished without diffusion and the two Poisson's ratios become time-independent.

The effective contribution of the LM and ILM layers on the coupling can be also discussed with respect to the model. A small decrease in the LM ionic concentration leads to a strong increase in the ILM ionic concentration, resulting in significant volumetric deformation compared to that of the LM layers (Figure I.7a). This model result is in accordance with the literature [START_REF] Huyghe | 3D non-affine finite strains measured in isolated bovine annulus fibrosus tissue samples[END_REF]Vergari et al., 2016) in which a localized strain of about two to three times higher in the ILM zone was reported on bovine annulus specimens. After the imbalance of the ionic concentration, the electrolytes return, during the relaxation process (

I.4.3. Limitations

Some limitations ought to be mentioned. One limitation of this study is the use of data obtained at room temperature (22°C) which does not replicate the real physiological environment in-vivo. The latter could influence osmotic effects, since osmotic pressure is related to temperature [START_REF] Emanuel | Osmosis and viscoelasticity both contribute to time-dependent behaviour of the intervertebral disc under compressive load: a caprine in vitro study[END_REF]. The biological reactions have not been considered in the current study which could influence the predicted plots of the electro-chemical activity in the soft tissue. Also, the current theory requires a more accurate description of the collagen network at different scales (Tavakoli et al., 2017(Tavakoli et al., , 2018)). Another limitation of this study concerns the mechanical conditions used to illustrate the model. Although it is representative of the disc compression mechanics, the model application should not be limited to the uniaxial circumferential data. The application on multiaxial loading paths is indeed mandatory to further verify the model (Kandil et al., 2020). Finally, because the annulus mechanics is different between species, it would have been more clinically relevant to conduct the model application on human data by including directional and regional effects.

I.5. Conclusion

In The mathematical structure of the constitutive model is general enough to be applicable to the timedependent response of all electro-chemically active multi-layered soft tissues, whatever biological or synthetic (such as hydrogels). In particular, the model could be employed to design prostheses including regional variation [START_REF] Jiang | Biomechanical response of a novel intervertebral disc prosthesis using functionally graded polymers: a finite element study[END_REF] and interaction with electro-chemical environment to better mimic the healthy disc function. Moreover, the model could be used to better understand the tissue response changes during aging and degeneration in response to changes in the tissue composition. These aspects will be the subject of future works.

Appendix I.A. Volumetric free energy function

The Cayley-Hamilton theorem applied to the chemical part chem C of the right Cauchy-Green deformation tensor gives:

3 2 2 12 chem chem chem I I J - + = C C C I (I.A1)
The Nanson relation dv J dV = , between volumes in the actual and initial configurations, dv and dV , respectively, may be written as: ( )

32 12 chem chem chem J dv I I dV = - + Ι C C C (I.A3)
Multiplying the two sides of Eq. (I.A3) by the unit vector 0 a , we get:

( )

32 0 0 0 1 2 0 chem chem chem J dv I I dV = - + Ι C C C a a a a (I.A4)
which can be also written as:

65 4 1 32 0 0 0 0 1 0 0 2 0 0 II I chem chem chem J dv I I dV   = - +   Ι C C C a a a a a a a a (I.A5)
The invariants in Eq. (I.11) of the main chapter body are used to re-write Eq. (I.A5) as follows: An exponential form of the free energies of LM and ILM layers was proposed in Eq. (I.19) of the main chapter body using Eq. (I.A6) or Eq. (I.A7), respectively, for dv dV :

2 1 2 2 1 exp 1 1 2 fluid d W K K K v dV    = - -         (I.A8)
The behavior of the formulated free energy function (I.A8) has to be verified to get physically reasonable responses and convergence of solutions [START_REF] Balzani | A polyconvex framework for soft biological tissues: adjustment to experimental data[END_REF][START_REF] Ehret | A continuum constitutive model for the active behaviour of skeletal muscle[END_REF]. In the case of a confined stretching

Appendix I.B. Structural features

In order to account the gradual variation of the thickness of layers, which increases from the outer annulus towards the inner annulus (Cassidy et al., 1989;Marchand and Ahmed, 1990;[START_REF] Hsu | Diffusion tensor microscopy of the intervertebral disc anulus fibrosus[END_REF]Holzapfel et al., 2005), a second-order polynomial function is introduced:

_2 LM i layer z ai bi c = + + 1,3,... im = (I.B1)
in which m is the total number of layers and, a, b and c are parameters governing the graduation. According to Tavakoli et al. (2017), the thickness of the LM phase is eight times that of the ILM phase, and the ILM thickness is given by: Chapter I: A microstructure-based model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus

40 _ 1 _ 1 _ 1 82 LM i LM i layer layer ILM i layer zz z -+  + =    2, 4,... 1 im =- (I.B2)
The total thickness of the multi-layered annulus is then given by:

Abstract

This article presents a multiscale model to predict deformation-induced damage and failure of human annulus fibrosus under multiaxial loading. In the modeling approach, formulated within the framework of nonlinear continuum mechanics, the hierarchical structure of the soft tissue is considered from the nano-sized collagen fibrils to the micro-sized oriented collagen fibers. At the macroscale, the multi-layered lamellar/inter-lamellar organization of the soft tissue is introduced by considering the effective interactions between adjacent layers. The stochastic process of progressive damage events operating at different scales of the solid phase is introduced for the extracellular matrix and the network of nano-sized fibrils/micro-sized fibers. The damage is made anisotropic due to lamellar oriented collagen fibers and special orientation distribution of the inter-fibrillar and inter-lamellar network of fibrils. The chemical-induced volumetric strain is also considered in our modeling approach to take into account the osmolarity effects along with the anisotropic timedependent transversal deformations. The capacity of the model is discussed using a few available stretching datasets till failure along circumferential and radial directions. Model predictions under tilted stretching, biaxial stretching and shearing are also presented to illustrate further the efficiencies of our modeling approach. This work shows for the first time the directional effects on annulus mechanics and failure in relation to external loading mode, structure features, damage events and hydration.

Keywords: Multiscale; Osmolarity; Multiaxial; Damage; Failure.

Nomenclature F : Deformation gradient tensor L : Velocity gradient tensor C : Right Cauchy-green strain tensor E : Green-Lagrange strain tensor σ : Cauchy stress tensor P : First Piola-Kirchhoff stress tensor R : Rotational basis-change tensor Ι : Identity tensor x , y , z : Cartesian coordinate vectors a : Unit vector of collagen fiber b : Unit vector of collagen fibril n: Unit vector at the inter-layer interfaces  : Fiber angle with respect to the xy-plane  : Fibril angle with respect to the az-plane  : Load angle with respect to xz-plane q : Fibril quantity  : Fibril bundle number

II.1. Introduction

The ineluctable deterioration of the intervertebral disc annulus fibrosus is a multiscale phenomenon due to the hierarchical structure of this soft tissue (Yu et al., 2002(Yu et al., , 2015;;Pezowicz et al., 2006aPezowicz et al., , 2006b;;Melrose et al., 2008;Schollum et al., 2008;Han et al., 2012;Vergari et al., 2016;Han et al., 2016;Tavakoli at al., 2016Tavakoli at al., , 2017Tavakoli at al., , 2020aTavakoli at al., , 2020b)). The irreversible degradation in annulus fibrosus operates from the nanoscale [START_REF] Tang | A constitutive model of soft tissue: from nanoscale collagen to tissue continuum[END_REF] to the macroscale (Werbner et al.,2019) and results in the development of multiple cracks and tears (Osti et al., 1992;Vernon-Roberts et al., 2007).

There is a considerable qualitative understanding of the microstructure and chemical factors that govern the response of the soft tissue. However, the coupling between these factors and mechanics is yet far to be fully understood in annulus fibrosus.

The introduction of the tissue structure at different scales in a continuum-based modeling framework would allow to propose predictive tools for mechanical damage and failure in a physically consistent Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 47 way. Over the years, macroscopic continuum-based damage models were proposed using two main approaches. The first one considers reduction factors incorporated in the free energy function (Rodriguez et al., 2006;[START_REF] Calvo | An uncoupled directional damage model for fibred biological soft tissues. Formulation and computational aspects[END_REF][START_REF] Pena | A constitutive formulation of vascular tissue mechanics including viscoelasticity and softening behaviour[END_REF][START_REF] Blanco | A regularised continuum damage model based on the mesoscopic scale for soft tissue[END_REF]Shahraki et al., 2015;[START_REF] Li | Damage models for soft tissues: a survey[END_REF][START_REF] Holzapfel | Modeling of damage in soft biological tissues[END_REF][START_REF] Li | Multiscale modeling of fiber recruitment and damage with a discrete fiber dispersion method[END_REF]Ghezelbash et al., 2020a). The second one introduces an energy limiter in the expression of the free energy function [START_REF] Volokh | Hyperelasticity with softening for modeling materials failure[END_REF][START_REF] Volokh | Modeling failure of soft anisotropic materials with application to arteries[END_REF]. No model until now is able to predict the multiaxial behavior of annulus fibrosus till failure and such a development remains a challenging task. The main difficulty comes from the anisotropic mechanics of a tissue with a complex hierarchical solid phase starting from the nanoscale (Tavakoli et al., 2020a(Tavakoli et al., , 2020b)). The degree of knowledge of the collagen network has considerably increased thanks to the technological advances in the imaging field (Yu et al., 2002(Yu et al., , 2015;;Pezowicz et al., 2006aPezowicz et al., , 2006b;;Melrose et al., 2008;Schollum et al., 2008;Han et al., 2012;Vergari et al., 2016;Han et al., 2016;Tavakoli at al., 2016Tavakoli at al., , 2017)). For a physically consistent constitutive representation, an accurate structure of the solid phase must be taken into account (Ghezelbash et al., 2020b) in interaction with the surrounding environment. In this regard, the fluid phase transfer by osmotic effect through the tissue solid phase affects largely the biochemical volumetric behavior (Derrouiche et al., 2019a(Derrouiche et al., , 2020a;;Yang and O'Connell, 2019;Feki et al., 2020).

The coupling between deformation-induced stress and swelling due to internal fluid content variation by osmosis is mandatory to quantitatively analyze the annulus fibrosus multiaxial mechanics (Kandil et al., 2020) along with the failure (Werbner et al., 2019).

The deterioration of the annulus fibrosus is indeed a multi-physics phenomenon since it implies biochemical degradation processes interacting in a complex manner with the external loading mode and the induced mechanical damage. A loss of electrolytic properties of extracellular matrix (ECM)

proteoglycan macromolecules occurs with the tissue degeneration, which results in a decrease of osmotic pressure [START_REF] Kiani | Structure and function of aggrecan[END_REF]Urban and Robert, 2003). In the meantime, high mechanical loads cause tissue degenerative changes resulting in the degradation of the proteoglycan

Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 48 macromolecules associated with the enzymes of the ECM metalloproteinase activity [START_REF] Iatridis | Region-dependent aggrecan degradation patterns in the rat intervertebral disc are affected by mechanical loading in vivo[END_REF]. The induced metabolic phenomena in the ECM cause a degradation of the osmotic properties of the annulus fibrosus associated with a decrease in water content and structure changes (Iatridis at al., 2006). These changes may alter the resistance of the annulus fibrosus to the mechanical-induced crack propagation (Werbner et al., 2019;[START_REF] Vergari | The effects of needle damage on annulus fibrosus micromechanics[END_REF]. Establishing annulus damagestructure-osmotic relationships is a fundamental task for a detailed understanding of the degeneration mechanisms, to reduce their effects and to potentially prevent dangerous degenerative cases from happening. For an accurate annulus behavior, the structure-chemo-mechanics relationship can only be seen by considering, in addition to the annulus stretching in the circumferential direction, representative of the disc compression mechanics, the multiaxiality of the external mechanical loading representative of the body movements (Kandil et al., 2020). The latter is the key point to better understand the origin of multi-directional micro-cracks and tears commonly observed in the annulus fibrosus [START_REF] Boos | Classification of agerelated changes in lumbar intervertebral discs[END_REF].

In this chapter, a multiscale model for deformation-induced damage and failure of annulus fibrosus is proposed by considering the structure at different scales, the micromechanical deformation processes leading to the final failure, the anisotropic nonlinear mechanics and the coupling with the osmotic effects. The proposed model, fully three-dimensional, is used to examine the multiaxial mechanics of the tissue and the directional effects. Uniaxial stretching, biaxial stretching and shearing are all being considered. The multiaxial model results are correlated to the microstructure features, to their progressive deformation-induced changes, to the tissue volumetric behavior and to the final failure.

This chapter is organized as follows. In Section II.2, the annulus multiscale model is presented within the framework of nonlinear continuum mechanics. In Section II.3, the model is verified for different loading paths and is used to discuss important aspects of the annulus mechanics in relation

Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 49 to volumetric effects and to damage accumulation in the tissue due to progressive local failure events. Section II.4 closes the chapter with some concluding remarks.

II.2. Model formulation

The following notation is used throughout the text. Tensors and vectors are respectively denoted by normal boldfaced letters and italicized boldfaced letters, while scalars and individual components of vectors and tensors are denoted by normal italicized letters. The superposed dot designates the time derivative.

II.2.1. Hierarchical organization of annulus fibrosus

Let us start by the presentation of the hierarchical organization of all discrete solid components. The micro-sized type-I oriented collagen fibers (OCF) represent up to 80% of the total volume of collagen network collagen  (Yu et al., 2002;Schollum et al., 2008;[START_REF] Sharabi | The mechanical role of collagen fibers in the intervertebral disc[END_REF]. This family of fibers is described geometrically by a unit vector a in the initial configuration with respect to the global benchmark of the disc as illustrated in where  is the fibers angle with respect to the circumferential direction. The nano-sized elastic fibers (NEF) consist in a complex network of collagen fibrils, present both in the LM zone and in the ILM zone. In the LM zone, they act as inter-fibrillar connections of type-I collagen fibers (Figure II.1c). In the ILM zone, they act as inter-lamellar connections of two adjacent lamellae (Figure II.1d) (Yu et al., 2002(Yu et al., , 2015;;Melrose et al., 2008;Schollum et al., 2008;Han et al., 2012). For symmetry reasons, the distribution of NEF is described geometrically by an The volume fraction of a fibril in each layer i is given by: _ 0 ,180 45 ,90 ,135

LM ILM LM LM i LM LM qq q q q        + =     =   =      and _ 0 0 ,180 45 ,90 ,135 ILM i ILM ILM q q      =     =  =      (II.4)
where LM q and ILM q are the respective total number of NEF in each layer i:

Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 53 0 ,180 In what follows, the features of the annulus mechanics are described along with a set of constitutive equations for the solid phase behavior using a physically realistic approach accounting for microstructure organization at different scales, damage and hydration.

LM LM ILM q q q      =   =+   

II.2.2. Free energy of the fluid phase and swelling-related changes

Local ionic diffusion is responsible for dilatation effects and nonlinear time-dependent transversal stretch response of the multi-lamellar annulus [START_REF] Sarkadi | Activation of ion transport pathways by changes in cell volume[END_REF][START_REF] Guizouarn | Swelling activation of transport pathways in erythrocytes: effects of Cl-, ionic strength, and volume changes[END_REF][START_REF] De Los Heros | Role of WNK kinases in the modulation of cell volume[END_REF]Derrouiche et al., 2020b). The derivation of the constitutive model must take into account this chemical coupling in the description of the finite-strain kinematics (Derrouiche et al., 2019b). In this framework, the basic quantity is the deformation gradient tensor mapping material tangent vectors in their reference configuration to their actual position in the deformed configuration. It is designated here as i F for each layer i of the multi-layered material element. The introduction of an intermediate (virtual) configuration assumed stress-free allows to multiplicatively decompose the deformation gradient tensor i F into a chemical-induced volumetric Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus

54 1 1 1 ii chem mech i i i i i i i chem chem chem mech mech chem - - - =+ LL L F .F F .F .F .F (II.7)
The determinant (Jacobian) ( ) det 0

ii J = F of the deformation gradient tensor i F gives:

i i i chem mech J J J = , ( ) det 1 ii mech mech J == F and
( )

det ii chem chem J = F (II.8)
where the term i chem J represents the tissue chemical-induced volumetric change. In view of the mechanical incompressibility of all solid components (ECM, OCF and NEF), the term i mech J is equal to one. We may define the total Green-Lagrange strain tensor i E , the mechanical (isochoric) Green- Lagrange strain tensor i mech E and the chemical-induced volumetric strain i chem  :

( )

1 2 ii =- E C I , ( )
1 2 ii mech mech =- E C I and 1 ii chem chem J  =- (II.9)
where I is the unit tensor,

T i i i = C F .F
is the total right Cauchy-Green strain tensor and

T i i i mech mech mech = C F .F
is the mechanical (isochoric) right Cauchy-Green strain tensor.

For each layer, the respective volume fractions of the three solid components (ECM, OCF and NEF) may be expressed as:

_0 i ECM i ECM i chem J   = , _0 i OCF i OCF i chem J   = and _0 i NEF i NEF i chem J   = (II.10) where _0 i ECM  , _0 i OCF  and _0 i NEF  are the initial volume fractions in each layer i.
The volume fraction of the fluid phase i fluid  in each layer i is then given by:

_ 0 _ 0 _ 0 1 i i i ECM OCF NEF i fluid i chem J     ++ =- (II.11)
The Cauchy stress tensor i σ , defined in the deformed configuration, is additively split into a mechanical part i mech σ and a chemical-induced volumetric part i chem σ

. The latter tensor is given by:

i i chem chem i chem W J  =  σI (II.
( ) ( ) ( ) 2 1 2 2 1 exp 1 1 2 i i i chem fluid chem K W K J K  = - - (II.13)
where 1 K and 2 K are material constants, 1 K and 12 KK being respectively the initial and maximal volumetric stiffness values.

II.2.3. Free energies of the solid phase

Each with respect to the corresponding deformation:

2 i i i iT i mech mech mech mech i mech W p  =-  σ F . .F I C (II.14)
where i p is determined from boundary conditions.

The mechanical free energy i mech W is given as a function of the free energies of the three solid components (ECM, OCF and NEF) using the volume fraction concept:

_ _ _ LM i i i i LM i i LM i mech ECM ECM OCF OCF NEF NEF W W W W    = + + (II.15)
in LM layer and,

_ _ ILM i i i i ILM i mech ECM ECM NEF NEF W W W  =+ (II.16) in ILM layer.
The free energy functions of the solid components (ECM, OCF and NEF) are separately described below.

II.2.3.1. ECM free energy

Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 56

The ECM proteoglycan macromolecular network is described as an isotropic deformable solid component with a free energy i ECM W expressed by a Neo-Hookean form:

( )

1 1 3 2 ii ECM ECM W G I =- (II.17) in which ECM G
is the ECM shear modulus and

( )

1 trace ii mech I = C
is the first stretch invariant.

II.2.3.2. OCF free energy

The

( ) ( ) ( ) ( ) 2 * 1 4 2 4 2 1 1 exp 1 1 2 i i i i OCF OCF C W H I C I C      = - - -    (II.18) below * OCF  and, ( ) ( ) ( ) 2 ** 3 4 4 4 4 1 2 i i i i i OCF OCF C W H I I I I           = - - + -    (II.19) above * OCF  . The stretch * OCF 
defines the stretch transition between the 'toe' and 'linear' regions of the stressstretch curve [START_REF] Quapp | Material characterization of human medial collateral ligament[END_REF][START_REF] Werbner | A novel method for repeatable failure testing of annulus fibrosus[END_REF][START_REF] Li | Multiscale modeling of fiber recruitment and damage with a discrete fiber dispersion method[END_REF].

II.2.3.3. NEF free energy

The inter-fibrillar and inter-lamellar network of fibrils is described by the free energy

ILM i ILM i i NEF WW    = =  (II.21)
For the same reasons invoked for OCF, the free energy of each fibril i W  is given by two distinct relations below and above a stretch * NEF  : Damage induces thus a progressive evolution of the effective amount of the solid components:

( ) ( ) ( ) ( ) 2 * 4 4 5 4 5 1 1 exp 1 1 2 i i i i NEF C W H I C I C       = - - -    (II.22) below * NEF  and, ( ) ( ) ( ) 2 ** 6 4 4 4 4 1 2 i i i i i NEF C W H I I I I              = - - + -    (II.
( ) _ 1 i i i ECM d ECM ECM d  =- , ( ) _ 1 i i i OCF d OCF OCF d  =- and ( ) _ 1 i i i d d     =- (II.
( )( )( ) _ 1 1 1 i i i i i chem d chem ECM OCF NEF d d d  = - - - (II.30) in which 1 n ii NEF d n d   =
 is the damage of NEF and n is the number of fibrils.

II.2.5. Summary of the model

A new model for the microscopic description of the annulus mechanics is constructed by using the constitutive equations described above. To formulate the overall mechanical response of the multilayered annulus under general three-dimensional loading conditions, the compatibility at inter-layer interfaces (Figure II.2a) must be respected while the above descriptions and definitions are used. In this regard, the continuity conditions in deformation and stress to be satisfied are:

_ _ 1 00 LM i i ILM i i + = F . F . nn and _ _ 1 LM i i ILM i i + = σ . σ . nn (II.31)
where 0 i n and i n are the arbitrary unit vectors between two adjacent layers in the initial and current configurations, respectively.

Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 60

The macro-deformation F and the macro-stress σ of the soft tissue are expressed as: where the deformation gradients and the Cauchy stresses in LM and ILM are given by their volumetric and mechanical parts:

1 __ _ 0 _ 0 1,3... 2,
_ _ _ LM i LM i LM i chem mech = F F .F and _ _ _ ILM i ILM i ILM i chem mech = F F .F (II.34) and _ _ _ LM i LM i LM i mech chem = σ σ σ + and _ _ _ ILM i ILM i ILM i mech chem = σ σ σ + (II.

II.3. Results and discussion

Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 62

In addition to the hierarchical organization of the collagen network presented in Section II.2.1, some direct structure inputs required by the modeling (namely dimensions, collagen orientation/content and water content) are extracted from well-documented papers of the literature (Skaggs et al., 1994;Acaroglu et al., 1995;Holzapfelet al., 2005;[START_REF] Baldit | Etude des interactions hydro-chimio-mécaniques dans les tissus biologiques: application à la nutrition du disque intervertébral[END_REF]Tavakoli et al., 2017) and are listed in The in-vitro experiments used for model identification and verification do not generally precise the exact disc region whereas variations in water content and collagen content/orientation actually exist.

Due to the lack of precise information, the selected structure features are representative of the annulus anterior area. Nine LM/ILM layers are considered resulting in a parallelepipedic material element of 10 × 4 × 4 mm 3 for circumferential direction, and 2 × 2 × 4 mm 3 for radial direction, in order to respect the sample size effect on the mechanical behavior of the annulus (Zhou et al., 2019;Kandil et al., 2020) The simulations were performed by using a constant strain rate of 0.001 s -1 .

Components Parameters Values

Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus estimated via a step-by-step approach to dissociate the different phenomena: osmotic effect, directional-dependent straightening of collagen networks and damage effect. The volumetric parameters, 1 K and 2 K , were identified using the UA circumferential data of Baldit (2013) describing the time-dependent volumetric response upon stretching followed by relaxation. The Unless explicitly otherwise stated, the values in Table II.2. will be employed in all model results, especially to generate the multiaxial model predictions under other loading paths.

63 Swelling Fluid 1 K 0.14   MPa 2 K 0.075   - Toe ECM ECM G 0.01   MPa OCF 1 C 26   MPa 2 C 12   - NEF 4 C 25   MPa 5 C 0.2   - Linear OCF 3 C 70   MPa * OCF  1.09   - NEF 6 C 1.5   MPa * NEF  1.7   - Damage ECM ECM  2   - ECM  5.5   MPa OCF OCF  5   - OCF  43   MPa NEF NEF  6   - NEF  850   MPa Table II.2. Model parameters.
 =      --  =+            (II.

II.3.1. Chemical-induced swelling

The swelling response was identified using experimental data under free-swelling and UA circumferential stretching (Acaroglu et al., 1995;[START_REF] Baldit | Etude des interactions hydro-chimio-mécaniques dans les tissus biologiques: application à la nutrition du disque intervertébral[END_REF]. For the latter, only moderate levels through very recent contributions (Derrouiche et al., 2019b;Kandil et al., 2019;Tamoud et al., 2021) and was shown as determinant in the time-dependent multiaxial response (Kandil et al., 2020). The ILM anatomic function is also to confer a shearing resistance to the annulus and the inclusion of the network of inter-lamellar fibrils in the model allows a more realistic material representation. Note that the ECM reorganization may induce viscoelastic features in the annulus fibrosus (Tavakoli and Costi, 2018). Nonetheless, the quantitative estimation of the coupling between viscoelastic effects and chemical-induced volumetric effects showed that the volumetric features are first-order factors in the annulus mechanics compared to the deformation-induced viscoelastic stress in ECM (Derrouiche et al., 2019b;Kandil et al., 2019Kandil et al., , 2020)). Although the solid components are described are in agreement with other data taken form the literature (Ebara et al., 1996). A direct relation between the annulus circumferential mechanics and the microstructure at the different scales can be pointed out. The straightening of the collagen network during stretching causes a (toe) nonlinear macro-response until the disappearance of the undulations of the collagen fibers at a strain of about 13%. Then, a nearly linear response appears up to reach the limiting extensibility of the collagen fibers that is seen at a strain of about 16%. A catastrophic stiffness decrease is then observed at a strain of about 18% due to damage accumulation that arises out of a wide fibril network degradation.

The circumferential stretching creates ionic imbalance and chemical stress due to chemo-

II.3.3. Biaxial stretching path

The biaxiality effects on the annulus stress-stretch curves till failure are presented in Figure II.12 for different biaxial states. As far as we know, no experimental data are available in the literature when circumferential stretching and radial stretching are combined. These results are for predictive purpose only. The computed results are shown for different strain loading ratios (= circumferential strain : radial strain, {1:0, 0:1, 1:0.5, 1:1 and 1:2}). The stiffness and ultimate properties are strongly influenced by the biaxial stretch ratio with a larger effect on the radial behavior than on the circumferential one. The radial stretchability of the annulus is considerably decreased under a biaxial stretching state. The difference in circumferential and radial behaviors is connected to the difference in how local failure occurs due to the collagen orientation. For future works, the annulus regional microstructure will be implemented in an axisymmetric version of the present (multiscale and multiaxial) model in the main goal to predict anisotropic damage and failure in a complete intervertebral disc model including nucleus-annulus interaction.

Appendix II.A. Deformation gradients

II.A.1. Chemical-induced swelling

The hydration effect is considered in our modeling approach. Under a free-swelling, the deformation gradients of a layer i are:

( ) ( )

_ _ _ diag , , 1
diag cos , cos , i i i i i i i chem x x y y z z         =+ FI (II.A2)
Such that the chemical-induced volumetric strain i chem  is given by: ( )( )( )

1 cos 1 cos 1 1 i i i i i i i chem x x y y z z          = + + + - (II.A3)
where  is a coefficient that changes according to the loading mode.

Guided by the experimental observations of [START_REF] Baldit | Etude des interactions hydro-chimio-mécaniques dans les tissus biologiques: application à la nutrition du disque intervertébral[END_REF] and Derrouiche et al. (2020b) on circumferentially stretched annulus specimens, the chemical deformation gradient (II.A2) was initially proposed in the previous chapter. The term  is introduced here to consider the loading angle effect on tilted stretching behavior and on biaxial stretching behavior.

For tilted stretching, the term  is given by:  = (II.A4) in which  is the angle between the loading direction and the circumferential direction.

For biaxial stretching, the term  is geometrically deduced from the biaxiality ratio: ( )

1 exp i i x x t  = - - , ( )
1 exp i i y y t  = - - and
( ) In what follows, the mechanical deformation gradients of each loading case, given in Figure II.5, are separately provided.

1 exp i i z z t  = - - (II.A6)

II.A.2. UA stretching path

The UA stretching of the multi-layered material element can be carried out along the circumferential direction, along the radial direction or along a direction tilted by an angle  with respect to the Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 78 circumferential direction. The mechanical deformation gradients may be expressed as a function of the applied stretch  and the stretch in the fiber direction f  :

( )

11 22 _ _ _ diag , , 1,3,... TT LM i i i mech y z f f f z y im        --   = =   F R . R . .R .R (II.A7)
in LM layer and, After a series of lengthy but straightforward derivations, we get the mechanical deformation gradients under UA stretching path:

1 1 2 2 _ _ _ _ diag , , 2, 4,... 1 T ILM i mech y y i i i chem x chem x chem x im         --         = = -                 F R . .R (
( )

22 3 _ _ 2 2 _ _ 2 3 2 3 3 22 3 _ _ 2 2
cos sin 1 cos cos sin cos sin 1 sin cos cos sin sin

i i i i i chem x chem x i chem x LM i i i i i mech i i i i i chem x chem x S S S S S S S S S S                             + -         = +           - +         F (II.A10)
in LM layer and,

Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus

79 22 _ _ _ _ _ _ 22 _ _ _ _ cos sin 1 0 cos sin 1 0 0 1 sin cos cos sin 0 i i chem x chem x i i chem x chem x ILM i mech i chem x i i chem x chem x i i chem x chem x                                      + -               =            - +          F          (II.A11)
in ILM layer.

The terms 1 i S , 2 i S and 3 i S are given by:

13 3 4 2 1 2 _ 2 21 2 _1 3 2 2 tan 1 tan 4 3cos 2 1 3 cos 1 cos sin i i i ii chem x ii i i i chem x i i i i i S SS S SS S             = - -           =+      =-  
(II.A12)

II.A.3. Biaxial stretching path

The biaxial mechanical deformation gradients may be expressed as a function of the stretches in the fiber axes

1 x f  , 2 x f  and 3 x f  :
( )

1 2 3 _ _ _ diag , , 1,3,... T x x x LM i i i mech z f f f z im      = = F R . .R (II.A13) in LM layer and, __ _ _ _ diag , , 2, 4,... 1 ii chem x chem z ILM i circ rad mech i i chem x circ rad chem z im         = = -       F (II.A14)
in ILM layer.

Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 80

The mechanical deformation gradient in LM layer is given by:

6 _ _ 65 _ 0 0 00 i circ i chem x LM i i i mech rad i chem z S SS             =         F (II.A15)
The terms 4 i S , 5 i S and 6 i S are given by:

( ) ( ) 2 2 _ 4 __ 4 4 2 5 _ 4 6 _ 2 cos sin tan 1 tan 2 tan 2 cos sin i i i chem z i circ cir ii chem x chem x rad ii ii circ i chem x i i i circ i i i chem x S S S S S                -  = - -      -   =+     =- (II.A16)

II.A.4. Shearing path

The annulus under shearing is subjected to the following mechanical deformation gradients: __ 1 0 0 01 01

LM i ILM i mech mech     ==    FF (II.A17)
where  is the shear stretch.

Appendix II.B. Multi-layered annulus behavior

The mechanical Cauchy stress tensor i mech σ is written as follows:

_ _ _ LM i i LM i LM i mech ECM OCF NEF =+ σ σ σ σ + (II.
ILM i ILM i NEF NEF    = =  σσ (II.B3)
In Eq. (II.14), the derivation may be expressed as follows:

4 4 1 1 4 4 i i i i i i i mech mech mech mech i i i i i i i mech mech mech mech I W W W I W I I I I            = + +        C C C C (II.B4) with 1 i i mech I  =  I C , 4 i i mech I   = C aa and 4 i i mech I    = C bb (II.B5)
After a series of straightforward derivations, the different parts of the Cauchy stress tensor are given by:

_ T i i i i i ECM ECM d ECM mech mech Gp  = - σ F .F I (II.B6) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 _ * 4 _ 1 4 2 4 _ * * 4 _ 3 4 4 1 2 1 exp 1 12 T T LM i i i i i i i i OCF OCF d mech mech OCF OCF LM i i i i i i i OCF OCF d mech mech OCF OCF H I C I C I H I C I I                = - - -         = - - +     σ F . .F σ F . .F aa aa (II.B7) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 _ _ * 4 _ 4 4 5 4 _ _ * * 4 _ 6 4 4 1 2 1 exp 1 12 T T LM i i i LM i i i i i i NEF d mech mech NEF NEF LM i i i LM i i i i i NEF d mech mech NEF NEF H I C I C I H I C I I                               = - - -         = - - +     σ F . .F σ F . .F bb bb (II.B8) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 2 _ _ * 4 _ 4 4 5 4 _ _ * * 4 _ 6 4 4 1 2 1 exp 1 12 T T ILM i i i ILM i i i i i i NEF d mech mech NEF NEF ILM i i i ILM i i i i i NEF d mech mech NEF NEF H I C I C I H I C I I                               = - - -         = - - +     σ F . .F σ F . .F bb bb (II.B9) ( ) 2 1 _ 2 _ exp i i i i chem fluid chem d chem d KK    = σ I (II.B10)
where i p is written from boundary conditions in each loading case:

Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 82 ( )

2 _ _ 2 2 _ 2 3 3 i chem x LM i i i i i ECM d ECM p G S S S         = + +        and 2 _ _ _ i chem x ILM i i ECM d ECM pG     =    (II.B11)
for UA stretching path, ( )

_ 2 2 _ 5 6 LM i i i i ECM d ECM p G S S  =+ and 2 __ _ _ ii chem x chem z ILM i i ECM d ECM circ rad pG     =    (II.B12)
for biaxial stretching path and, ( )

2 _ 1 ii ECM d ECM pG  =+ (II.B13)
for shearing path.

Modeling multiaxial damage regional variation in human annulus fibrosus8

Abstract

In the present article, a fully three-dimensional human annulus fibrosus model is developed by considering the regional variation of the complex structural organization of collagen network at different scales to predict the regional anisotropic multiaxial damage of the intervertebral disc. The model parameters are identified using experimental data considering as elementary structural unit, the single annulus lamellae stretched till failure along the micro-sized collagen fibers. The multilayered lamellar/inter-lamellar annulus model is constructed by considering the effective interactions between adjacent layers and the chemical-induced volumetric strain. The regional dependent model predictions are analyzed under various loading modes and compared to experimental data when available. The stretching along the circumferential and radial directions till failure serves to check the predictive capacities of the annulus model. Model results under simple shear, biaxial stretching and plane-strain compression are further presented and discussed. Finally, a full disc model is constructed using the regional annulus fibrosus model and simulations are presented to assess the most likely failed areas under disc axial compression.

Keywords:

Annulus fibrosus; Regional dependence; Multiscale structure; Multiaxial damage and failure.

III.1. Introduction

The annulus fibrosus of the intervertebral disc is probably one of the most extraordinary composites that the nature produces. The annulus soft tissues exhibit an organized structural arrangement of collagen network at different scales of hierarchical organization (Yu et al., 2002;Pezowicz et al., 2006aPezowicz et al., , 2006b;;Melrose et al.,2008;Schollum et al., 2008). The detailed information was only appreciated very recently (Han et al., 2012;Yu et al., 2015;Vergari et al., 2016;Han et al., 2016;Tavakoli et al., 2016Tavakoli et al., , 2017Tavakoli et al., , 2020aTavakoli et al., , 2020b)). The structural arrangement of collagen network varies throughout the disc which leads to a regional variation of its mechanics. In this regards, various damage and failure mechanisms operate in the different disc regions (Osti et al., 1992;[START_REF] Boos | Classification of agerelated changes in lumbar intervertebral discs[END_REF]Vernon-Roberts et al., 2007). In order to prevent all kinds of sudden disc injuries, it is necessary to understand the failure risks considering the regional variation under various static loading conditions. Fruitful experimental studies have been undertaken over the years allowing the assessment of the annulus intrinsic constitutive response in relation to disc region and mechanical loading (Acaroglu et al., 1995;Ebara et al., 1996;Fujita et al., 1997;[START_REF] O'connell | Human annulus fibrosus material properties from biaxial testing and constitutive modeling are altered with degeneration[END_REF].

Modeling multiaxial damage regional variation in annulus can provide a better understanding of the failure onset and location in the disc, beyond experimental observations. Modeling mechanics of the healthy annulus was developed in several papers considering different approaches as physically realistic as possible, to name a few recently published (Derrouiche et al., 2019b;Kandil et al., 2019;[START_REF] Ghezelbash | Modeling of human intervertebral disc annulus fibrosus with complex multi-fiber networks[END_REF]Kandil et al., 2020;Castro and Alves, 2020;Komeili et al., 2021;Tamoud et al., 2021a;Zhou et al., 2019[START_REF] Zhou | Fiber engagement accounts for geometry-dependent annulus fibrosus mechanics: a multiscale, structure-based finite element study[END_REF]. The quantitative prediction of the annulus damage is less reported. Damage criteria were proposed to predict annulus failure under different loading conditions (Qasim et al., 2012(Qasim et al., , 2014;;Shahraki et al., 2015;Subramani et al., 2020). Nevertheless, the interest of this approach is questionable because it is intrinsically decoupled from the intrinsic properties and disconnected from the real microstructure and its evolution. The progressive loss of collagen network integrity implies a damage accumulation in annulus (Tamoud et al., 2021b).

III.2.1. Multiscale structure

The model considers the hierarchical organization of the solid phase starting from the nanoscale (Tavakoli et al., 2020a(Tavakoli et al., , 2020b) ) as well as its regional variation. As schematically represented in OCF and two adjacent LM between them (Yu et al., 2002;Melrose et al., 2008). The network of OCF and NEF will be described geometrically with respect to the global benchmark of the disc, respectively, with the angle  with respect to the circumferential direction and with the angle  with respect to the radial direction in the parallel plane of the fibers.

Figure III.2. Multiscale view of the intervertebral disc annulus. The organization, in terms of orientation and distribution, of this multiple-network medium is regional dependent.

Quantities and orientations of the constituents can be found in the Appendix III.A. The latter is based on well-documented papers for the four main disc regions (Skaggs et al., 1994;Acaroglu et al., 1995;Ebara et al., 1996;Holzapfel et al., 2005;Tavakoli et al., 2020aTavakoli et al., , 2020b)).

III.2.2. Damaged free energy functions

The randomly oriented proteoglycan macromolecules of the extracellular matrix (ECM), the OCF and the (inter-fibrillar and inter-lamellar) NEF are all being considered to model the solid phase. The damage variable d varies between 0 for the virgin state and can increase up to 1 in the case of fully damaged (cracking) state. It is taken as a monotonic increasing function of loading, i.e. the damage is considered as an irreversible phenomenon without possible recovery. Note that the anisotropy of the damage mechanisms developed in the annulus is due to the collagen network organization. The model is based on the general damage mechanics framework for the hyperelasticity theory that introduced hyperelastic laws derived from free energy functions [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF][START_REF] Miehe | Discontinuous and continuous damage evolution in Ogden-type large strain elastic materials[END_REF]. The mechanical free energy mech W is expressed as a function of the ECM free energy ECM W and the collagen networks free energies W  and W  : region is due to the collagen straightening during stretching that is followed by the "linear" region when the undulations disappear. The free energies W  and W  are expressed by two distinct relations below and above * ,   to describe the transition between 'toe' and 'linear' regions:

( ) ( ) ( ) 1 1 1 mech ECM ECM ECM OCF W d W d W d W          = - + - + -  (
( ) ( ) ( ) ( ) 2 * 1 , 4 , 2 4 , , , 2 1 1 exp 1 1 2 C W H I C I C              = - - -      (III.4)
C W H I I I I                         = - - + -      (III.5)
where 1 C , 2 C and 3

C are material constants, ( )

( ) ( ) ( ) 2 * * , 1 4 , 2 4 , 1 exp 1 C I C I        = - - (III.6) ( ) ( ) ( ) 2 * * 1 , 4 , 2 4 , 2, exp 1 1 2 C I C I C          = - - - (III.7)
The fourth stretch invariants 4 I  and 4 I  are given by: [START_REF] Karsaj | The modelling of fibre reorientation in soft tissue[END_REF][START_REF] O'connell | Human annulus fibrosus material properties from biaxial testing and constitutive modeling are altered with degeneration[END_REF][START_REF] Ghezelbash | Modeling of human intervertebral disc annulus fibrosus with complex multi-fiber networks[END_REF].

III.3. Simulation vs. experiments

In this section, after its identification using available experimental data, the predictability of the above-described model is discussed. The constitutive equations were coded in MATLAB software and reduced to different loading forms.

III.3.1. Single lamellae identification

The intrinsic constitutive response of the annulus under UA loading was experimentally determined The identification was performed numerically using the MATLAB software to minimize the following objective function: in which the letters with an overbar denote the experimental data with  the number of considered data points. Table III.1. provides the obtained values for the OCF properties. Note that the regional dependence of the single lamellae stiffness is explained by a difference in the initial undulated state and the OCF content (Skaggs et al., 1994;Holzapfel et al., 2005). III.2.

Parameters

The simulated curves of the single lamellae are in good agreement with the experiments of Holzapfel et al. (2005). The results are shown until the complete loss of OCF load-bearing capability in order to compare the ultimate points with the experimental data of Skaggs et al. (1994). Except PI, it can be observed that the simulated stress and strain thresholds fall within the experimental standard deviation intervals.

Constituents Parameters

Values ECM III.2. ECM, NEF and swelling parameters.

ECM G 0.01   MPa ECM  2   - ECM  5.5   MPa NEF 1 C 12   MPa 2 C 0.2   - 3 C 1.5   MPa *   1.7   - NEF  6   - NEF  400   MPa Swelling 1 K 0.14   MPa 2 K 0.0075   - Table

III.3.2.Multi-lamellae multiaxial predictions

Table III.3. presents the range of mechanical paths used to check the predictive capacities of the multi-lamellae model. It includes UA stretching in the circumferential direction, UA stretching in the radial direction, SS loading, BA stretching and PSC loading 10 . Except the latter loading (in which the sample is confined in the radial direction), the samples are loaded in specific directions while the others are free. All the samples are considered to be immersed in a physiological solution at a concentration of 0.15 M NaCl. When experimental data are available, a quantitative comparison with the model results is presented. The dimensions of the annulus specimens are given in Table III.3. for the different loading conditions and disc regions.

10 These different loadings are representative of the disc mechanics. UA stretching along the circumferential direction and PSC are representative of the disc axial compression (Shah et al., 1978). UA stretching along the radial direction may be related to the disc shearing whereas SS may be related to the disc axial twist. BA stretching may be caused for example by a disc bending (i.e. flexion, extension and lateral flexion) (Costi et al., 2007). Chapter III: Modeling multiaxial damage regional variation in human annulus fibrosus 101 (without statistical analysis), the model-experiments comparison can be considered as acceptable.

In accordance with experiments, the outer regions have much higher stiffness than the inner regions, due mainly to regional variation in OCF angle. The simulations are conducted until stress drop corresponding to the failure stage. The simulated ultimate properties are reported in Ebara et al. (1996), circle symbols with standard deviation intervals: experimental data of Acaroglu et al. (1995)).

The failure strain predictions fail within the standard deviation intervals. Although the predicted failure stress can be found outside the experimental intervals, the predictive trends provide a regional difference in agreement with the experimental data. This important result confirms that the model provides a quite proper description of the regional effects on the relationship between the actual microstructure, the damage events and the macroscopic ultimate properties. Recall that the model connects macroscopic response and hydration (volumetric) effects. The capacity of the model to capture the regional variation of the transverse responses is shown in Figure III.6. The LM-ILM interaction highlights the crucial role that the OCF regional orientations play in the transverse

Chapter III: Modeling multiaxial damage regional variation in human annulus fibrosus 102 response (Kandil et al., 2020). In particular, regarding the possibility of an apparent auxeticity manifested by the positive transversal strains of the lamellae plane [START_REF] Baldit | Etude des interactions hydro-chimio-mécaniques dans les tissus biologiques: application à la nutrition du disque intervertébral[END_REF]. Thanks to the incorporation of the actual nanostructure materialized by the network of fibrils, the model can predict the regional behavior of the annulus in the radial direction. (1997). By taking into account accurately the connections between layers, the model allows also to predict the shearing resistance. Figure III.7b shows such predictions for which it can be observed a strong regional dependence. In addition to amount, type and orientation of the collagen network, the regional variation of the NEF arrangement leads to different predicted stiffness and failure thresholds according to the disc region. The connection of collagen network between lamellae from nanoscale made it possible to represent precisely the damage processes in the LM and ILM connections in the different annulus regions. The model shows that the different regional thresholds

Chapter III: Modeling multiaxial damage regional variation in human annulus fibrosus 103 of radial and shear failures are caused by the regional arrangement of the NEF. These results are in agreement with the qualitative data analyzed in Fujita et al. (1997). Indeed, it was reported from image analysis that inter-layer connections such as NEF contribute to the annulus radial resistance (Pezowicz et al., 2006a;Tavakoli and Costi, 2018).

( The directional effects are quite well captured by the model, especially the much higher stretchability along the radial direction than along the circumferential direction. The BA stress-strain responses are presented in Figure III.9 for different strain ratios (= circumferential strain : axial strain, {1:1, 1:2 and 2:1}). In order to further evidence the directional effects, Figure III.10 presents the PSC stress-stretch responses for different strain ratios. The strain ratio has an important effect on the stiffness and the ultimate properties, with a stronger effect on the axial properties than the circumferential ones. The posterior zones show higher stiffness than the anterior zones. That can be justified by the fact that the OCF angles are close to 45°, and that the forces are well shared according to the two loading directions (circumferential and axial). Although the location was not

III.4. Full disc model

In this section, predictions of a full disc model are presented using the regional dependent constitutive representation of multi-lamellae and local damage mechanisms are clarified in order to get a basic insight about the starting areas of failure modes.

III.4.1. Disc construction and boundary conditions

A lot of finite element disc models were developed in the literature (Qasim et al., 2012(Qasim et al., , 2014;;Castro and Alves, 2020;Subramani et al., 2020;Komeil et al., 2021;Kandil et al., 2021). Beyond fundamental differences in the constitutive representation of the disc tissues, the time of model design and the time of calculation remain the main common weak points of all previous disc models.

In the final objective to create time and cost-efficient patient-specific models of human discs, the present disc modeling strategy is purely analytical and again coded in MATLAB. The volume element of multi-lamellae, constitutively designed in the Cartesian space (x, y, z), was transferred to the Cylindrical space (r, θ, z) in order to consider the full disc geometry as illustrated in Figure III.11a. From the theoretical point of view, the constitutive equations were re-written within a finitestrain kinematic framework in Cylindrical coordinates. The reader is referred to the work of Holzapfel et al. (2000) for the general theory on tube-like solid mechanics. In the present chapter, we restrict the loading to a uniaxial compression at the disc external boundary. The disc is thus seen as a cylindrical solid loaded in compression on the superior surface and with a pressure in the cylindrical internal wall (Figure III.11b). The latter simulates the nucleus swelling applied to the internal wall of the disc annulus; the applied vertical axial displacement was translated to an internal pressure by considering the nucleus bulk modulus. The regional variation in structure and properties of the disc annulus was taken into account using the constitutive representation of each human annulus region. In order to ensure smooth transitions of the microstructural features (ECM, NEF, OCF and fluid) and properties from a disc annulus region to another, with a radial and circumferential variation, linear interpolations were considered. An anatomically accurate three-dimensional model of a lumbar disc at the L1-L2 level was constructed using quantitative anatomical data taken from the work of Holzapfel et al. (2005). The L1-L2 disc geometric contours were approximated using special mathematical functions considering especially non-symmetric thicknesses and heights (Figure III.11a). Note that they may be adapted to consider the shape and size-changing with the spinal level and age (Amonoo-Kuofi, 1991). The annulus was subdivided into 15 different layers according to the L1-L2 disc dimensions provided by Holzapfel et al. (2005). The annulus size is considered to be approximately 50% of the disc volume (Violas et al., 2007). It is worth noticing that a notable difference exists between invivo and in-vitro responses regarding fluid exudation. In-vivo, the healthy disc completely recovers its height and intradiscal pressure after long diurnal loading. In-vitro, several studies reported losses in the disc height and intradiscal pressure (Feki et al., 2020). The chemical-induced volumetric strain was then affected by a multiplying factor. Several numerical simulations were performed with different values of the multiplying factor ranged from 0 to 100%, and a value of 50% was retained in the present work. During the simulations, a vertical axial displacement was applied on the superior surface while the inferior surface was constrained (Figure III.11b). The simulations were performed at a displacement rate of 0.01 mm/s in order to maintain quasi-static loading conditions.

III.4.2. Damage fields

The damage process in the disc core simulated by the model is recorded in Excessive damage accumulation in these zones suggests the possibility of tear creation that may lead to an anterior herniation as experimentally reported in (Osti et al., 1992;[START_REF] Schollum | The influence of concordant complex posture and loading rate on motion segment failure: a mechanical and microstructural investigation[END_REF][START_REF] Sapiee | Sagittal alignment with downward slope of the lower lumbar motion segment influences its modes of failure in direct compression: A mechanical and microstructural investigation[END_REF]. Besides, the ECM and NEF damage patterns are presented in the two layers in order to appreciate the specific role of the ILM matrix and compare it to the adjacent LM. For both ECM and NEF, two critical zones appear in the lamellae on either side of the annulus with very low intensities. With relatively higher intensities, critical damage spots in the ILM matrix (for both ECM and NEF) are firstly located in the anterior side of the annulus and are then progressively more likely located in the neighborhood of the anterior inner part of the annulus. Chapter III: Modeling multiaxial damage regional variation in human annulus fibrosus 111 Quite interestingly, a confined distribution of the ILM damage can be observed in the lateral anterior inner and posterior inner sides of the annulus. The latter is often judged as one of the most dangerous disc failure modes with the highest risk of disc delamination and failure propagation in the radial direction of the disc as experimentally reported in (Osti et al., 1992;[START_REF] Schollum | The influence of concordant complex posture and loading rate on motion segment failure: a mechanical and microstructural investigation[END_REF][START_REF] Sapiee | Sagittal alignment with downward slope of the lower lumbar motion segment influences its modes of failure in direct compression: A mechanical and microstructural investigation[END_REF]. The distribution of the fibrils makes the damage anisotropic. To illustrative this model capacity, the damage of the fibrils is plotted in

III.1. Discussion and concluding remarks

In the present chapter, a novel human annulus model was proposed to constitutively relate the regional variation in structural arrangement of collagen network at different scales and the multiaxial mechanics till failure. The model provides a useful tool to predict a variation of the annulus stiffness and ultimate properties with both disc region and loading mode. The different loading modes were simulated and chosen to better understand the regional behavior until the failure of the solid components in interaction with the chemical-induced volumetric response. The capabilities of the model were analyzed and an acceptable agreement between predicted values and experimental data was highlighted. A full disc model was constructed using the regional dependent constitutive representation of multi-lamellae and damage fields in the compressed disc were analyzed to assess the areas where the risk of failure is the highest.

The present model can be considered as the stepstone for the long-term prediction of the disc dysfunctions. During different complex loading scenarios of the spine, the disc dysfunctions may occur either by: i) mechanical damage due to a high level of loads (Acaroglu et al., 1995;Ebara et al., 1996) or relatively moderate loads repeated several times which cause damage accumulation [START_REF] Green | Tensile properties of the annulus fibrosus[END_REF][START_REF] Iatridis | Mechanical damage to the intervertebral disc annulus fibrosus subjected to tensile loading[END_REF]Qasim et al., 2012;Subramani et al., 2020), ii) by an agerelated biological degeneration [START_REF] Thompson | Preliminary evaluation of a scheme for grading the gross morphology of the human intervertebral disc[END_REF]Urban and Roberts, 2003;[START_REF] Cegoñino | A constitutive model for the annulus of human intervertebral disc: implications for developing a degeneration model and its influence on lumbar spine functioning[END_REF] or by the cellular disorders, including the cell death, the production of inflammatory mediators and a shift towards catabolism (Bruehlmann et al., 2004;[START_REF] Walter | Complex loading affects intervertebral disc mechanics and biology[END_REF]Bloom et al., 2021).

Nonetheless, it is difficult to identify or precisely model the disc failure taking into account the two causes of dysfunction at the same time. Since our study is based on the in-vitro annulus results, it cannot be directly extrapolated to the in-vivo disc response. Indeed, the understanding of the disc mechanobiology must be taken into account by the model considering especially nutrition and cellular activity along with degenerative aspects (Bruehlmann et al., 2004;[START_REF] Walter | Complex loading affects intervertebral disc mechanics and biology[END_REF]Bloom et al., 2021). The model could be then a real help to provide solutions for remodeling and

The volume fractions of fibril bundles d i   in the two layers are estimated as follows: 0 ,180 45 ,90 ,135

LM ILM LM LM LM LM qq q q q        + =     =   =      and 0 0 ,180
45 ,90 ,135

ILM ILM ILM q q      =     =  =      (III.A1)
where LM q and ILM q are the respective total number of NEF in each layer i: Quantities and orientations of the other constituents are extracted from well-documented papers for the four main disc regions (Acaroglu et al., 1995;Ebara et al., 1996;Skaggs et al., 1994;Holzapfel et al., 2005). Theses inputs required by the model are provided in To consider the effective interactions between adjacent layers, the compatibility condition on the deformation gradients is introduced at the LM/ILM interfaces: _ _ 1 00

0 ,180 LM LM ILM q q q       =   =+   
LM i i ILM i i + = F . F . nn (III.B2)
where 0 i n is the arbitrary unit vector between two adjacent layers in the initial configuration.

From the first chapter, the macroscopic Cauchy stress tensor σ is expressed as:

( ) ( ) 1 _ _ _ _ 1,
det i i chem chem J = F is the chemical-induced volumetric change.
The stress equilibrium condition at the LM/ILM interfaces is given by:

_ _ 1 LM i i ILM i i + = σ . σ . nn (III.B5)
where i n is the arbitrary unit vector between two adjacent layers in the current configuration.

The initial thickness where m is the number of layers.

For each disc region, structural parameters are extracted from well-documented papers (Tavakoli et al., 2017;Holzapfel et al., 2005) and are listed in Table III [START_REF] Boos | Classification of agerelated changes in lumbar intervertebral discs[END_REF]

IV.1. Introduction

The complex structural organization of interpenetrating collagen fibrillar networks in the intervertebral disc has been only revealed and appreciated recently (Tavakoli et al., 2020a(Tavakoli et al., , 2020b)).

The disc presents a hierarchical lamellar structure in which nano-sized elastic fibers (NEF) connect both two micro-sized type-I oriented collagen fibers (OCF) of the concentric lamellae (LM) and two adjacent LM between them by crossing the inter-lamellar (ILM) ground matrix. Besides, the regional variation in organization of collagen fibrillar networks throughout the disc leads to a heterogeneous mechanics of the disc (Tamoud et al., 2021a). In the literature, a lot of disc models were proposed (Qasim et al., 2012(Qasim et al., , 2014;;Shahraki et al., 2015;Castro et al., 2020;Subramani et al., 2020;Kandil et al., 2021;Komeili et al., 2021). A short overview shows that the disc is most naturally modeled using the finite element method. While the OCF network is commonly introduced in the current finite element models, the interpenetrating NEF network is still largely unappreciated and neglected. The consideration of the collagen fibrillar networks in disc modeling strategies is essential for an accurate prediction of the local and overall disc responses. To reflect the real disc mechanics, any model should also consider changes of the disc shape and size with the spinal level and age [START_REF] Amonoo-Kuofi | Morphometric changes in the heights and anteroposterior diameters of the lumbar intervertebral discs with age[END_REF][START_REF] Kim | Changes in spinal canal diameter and vertebral body height with age[END_REF], the disc hydration state (Urban and Maroudas, 1981;Derrouiche et al., 2019)

IV.2. Model

IV.2.1. Tube-like kinematics

This section provides the kinematics framework of the tube-like solid mechanics. The general theory is detailed by [START_REF] Ogden | Non-linear elastic deformations[END_REF] and Holzapfel et al. (2000) for arterial applications in which some adaptations to discs are here brought especially considering regional variations in morphology and in structure-property relationships. Combined axial/torsion loadings at the external boundary of the disc will be considered in the present theoretical developments as illustrated in Figure IV.1. An internal pressure simulating the nucleus swelling may be activated when the compression is considered (Shah et al., 1978;Derrouiche et al., 2020). The nucleus is referred as the central portion i = 0 and the annulus is seen as a multi-layered tube divided into n rings i = 1, … n. The half-disc topology is provided in Figure IV.2 in superior and sagittal views.

In the initial configuration, the solid is referenced in cylindrical polar coordinates (R, Θ , Z): The kinematics of an axially/torsionally loaded tube-like solid is given as follows [START_REF] Ogden | Non-linear elastic deformations[END_REF]):

, i i i inner outer R R R    ,
( ) The deformation gradient tensor i F of each ring may be thus written as (Ogden, 1997): ( )

1 , i i i i r f R r - = , i i Z H   =+
i i i i i i r r z z z i i i i i r r z z r R R Z z Z           =  +  +  +       F e E e E
_ _ _ _ 1 1 1 i i i i z mech r r z z z i i i i i i chem r chem chem chem z r r r R R H             =  +  +  +   F e E e E
0 1 chem  =+ F I (IV.8)
where  is the chemical expansion at the equilibrium state and  is a function of time ( )

1 exp t   = - - in which  is the chemical expansion rate.
After a series of straightforward derivations, the nucleus actual radius 0 r is given by: ( ) ( ) After a series of straightforward derivations, the following expression is obtained for the function i f using the continuity condition between the adjacent rings ( 1 ii inner outer rr -=

0
2 cos diag , , 3
i i i i i i i i chem r r z z         =+ FI (IV.
) and the mechanical incompressibility:

( ) ( ) ( ) ( ) 2 2 2 1 , 1 i i i i i chem inner outer inner outer i z f R R r   - + = - + (IV.11) where ( ) det 1 i i chem chem  =-F
is the chemical-induced volumetric strain in the ring i .

The annulus is assembled by considering the effective interactions between adjacent layers via the compatibility conditions on the deformations i F and on the Cauchy stresses i σ :

_ _ 1 0 0 LM i i ILM i i + = F . F . nn and _ _ 1 LM i i ILM i i + = σ . σ . nn (IV .
)( ) ( ) 4 00 1 1 1 1 1 frontal R R R f R R f f f H f H f            - = + - - + +         -        (IV.16)
in which H  is the Heaviside function:

2 1 3 H    = and 2 0 3 H    = (IV.17)
Chapter IV: A fully three-dimensional model of interpenetrating collagen fibrillar networks for intervertebral disc mechanics 130

The function f  is expressed as:

2 2 exp 3 f      = --      (IV.18)
The term R f is a shape factor of the nucleus given by the ratio between the sagittal and frontal initial radii 0 sagittal R and 0 frontal R :

0 0 1 sagittal R frontal R f R =- (IV.19)
In the limit case 00 sagittal frontal RR = , Eq. (IV.16) becomes 00 frontal RR = and the disc has a circular basis.

The initial radius i R in each ring is calculated as follows:

0 1 2 i i k outer k t RR = =+  (IV.

IV.2.3.2. Thickness regional variation

The thicknesses of the LM and ILM layers are given by the following linear equations:  = and the posterior region ()  = whereas the terms AI t and PI t are their counterparts for the inner side of the ring 1 i = . The linear interpolation is performed with an angular increment  of 1°. The same increment is fixed for all the other spatial discretizations in the circumferential direction. Note that the radial disctretization is fixed by the total number of layers n.

IV.2.3.3. Height regional variation

The regional variation of the disc height is taken into account using the following linear equation: The regional variation in quantities was taken into account using linear interpolations similar to those introduced above for the macrostructure features. The interpenetrating collagen fibrillar networks are spatially disposed in the model using two representative angles with respect to the global benchmark of the disc. The angles i  are introduced to dispose the micro-sized OCF with respect to the circumferential direction of the disc with alternate orientations between successive lamellae. The angles i  is introduced to dispose the nano-sized NEF network with respect to the radial direction in the OCF parallel plane.

The OCF in each LM ring is described geometrically by the unit vector i a expressed, in the current configuration, with respect to the global benchmark of the disc:

cos sin i i i  =+ az  (IV.29)
In order to introduce the change in the OCF angle induced by the circumferential stretch   , a measure of the reorientation may be considered geometrically using the following relation: The unit vector i b of each fibril is expressed, in the current configuration, with respect to the global benchmark of the disc as follows: sin cos cos sin cos

i i i i i      b = r + + z  (IV.31)
Again, the stretch-induced reorientation of the fibrils array is geometrically considered in each ring using the following expression: (2017,2020).

Chapter IV: A fully three-dimensional model of interpenetrating collagen fibrillar networks for intervertebral disc mechanics 135

The regional variation in network orientation was considered using linear interpolations with radial and circumferential variations as detailed in Appendix IV.A. The collagen fibrillar organization is provided in Figure IV.4 in the form of contour plots in the superior view of the disc for the OCF network and in the form of contour plots in the sagittal view of the disc for the NEF network.

IV.2.5. Constitutive equations

The local structure-mechanics relations are described using a recently developed fully threedimensional constitutive model (Tamoud et al., 2021a). The model efficiencies to reproduce regional dependent multiaxial mechanics till failure of the human lumbar disc annulus were shown.

In what follows, we present a summary of the constitutive equations of both healthy and damaged components.

IV.2.5.1. Healthy components

In virtue of the chemo-mechanical coupling, the total free energy W is additively split into a mechanical part mech W and a chemical part chem W :

mech chem W W W =+ (IV.33)
The mechanical free energy mech W is obtained using as ensemble-volume averaged homogenization procedure a volume-weighted average of the different solid components: are functions that ensure continuity between 'toe' and 'linear' regions:

mech OCF ECM ECM W W W W        = + +  (IV.
) ( ) ( ) ( ) 2 * 1 , 4 , 2 4 , , , 2 1 1 exp 1 1 2 C W H I C I C             = - - -    (IV.35) below * ,   and, ( ) ( ) ( ) 2 ** 3 , 4 , 4 , 4 , , 4 , , , , 1 2 C
W H I I I I                        = - - + -    (IV.
( ) ( ) ( ) 2 ** , 1 4 , 2 4 , 1 exp 1 C I C I        = - - (IV.37) ( ) ( ) ( ) 2 ** 1 , 4 , 2 4 , 2, exp 1 1 2 C I C I C          = - - - (IV.38)
The ECM free energy ( ) ( )

2 1 2 2 1 exp 1 2 chem fluid chem K WK K  =- (IV.40)
where 1 K and 2 K are material constants, 1 K and 12 KK being the initial and maximal volumetric stiffness values.

The overall Cauchy stress tensor annulus σ in the annulus may be calculated as: (  )

( ) ( ) 1 _ _ _ _ 1,
1 1 1 i chem i i ring n j chem j j t t    = + = +  (IV.42)
The stress in the nucleus is given by Eq. (IV.15).

IV.2.5.2. Damaged components

The above free energies may be modified to consider the internal stiffness degradation due to the development of mechanical-induced damage (see chapter II):

( ) ( ) ( ) A strong coupling exists between hydration and mechanical-induced damage (Werbner et al., 2019). In ordert to consider the latter coupling, the chemical-induced volumetric strain

1 1 1 mech OCF ECM ECM ECM W d W d W d W          = - + - + -  (

IV.3. Results and discussion

The model was entirely coded in MATLAB software. At this step of the model development, predictions on axially/torsionally loaded discs can be analyzed. The simulations were performed at a displacement rate of 0.01 mm/s and a twist rate of 0.1 °/s in order to maintain quasi-static loading conditions. The maximum ranges are changed according to the observed local fields in the disc core.

IV.3.1. Model inputs

Table IV.1 provides the anatomical data used as direct inputs for a lumbar disc model geometry taken from cadaver disc images at the L1-L2 level of the Holzapfel et al. (2005) work. The annulus size is considered to be approximately 50% of the disc volume (Violas et al., 2007) and is subdivided into 15 different lamellae that are intercalated with 14 ILM. The base radius 0 R is calibrated from the cadaver disc images of Holzapfel et al. (2005) resulting in a ratio between the sagittal and frontal radii of 00 sagittal frontal RR = 0.68 with a frontal radius of 0 frontal R = 15 mm. The values of the OCF orientation (with alternating signs between successive lamellae), the collagen volume fraction and the fluid volume fraction are extracted from previous documented experimental findings (Skaggs et al., 1994;Acaroglu et al., 1995;Ebara et al., 1996) and are listed in Table IV

IV.3.2. Overall disc response

The simulated load curve obtained on a healthy disc upon compression are compared in Figure

IV.5a with available experimental observations taken from the literature [START_REF] Virgin | Experimental investigations into the physical properties of the intervertebral disc[END_REF][START_REF] Markolf | The structural components of the intervertebral disc: a study of their contributions to the ability of the disc to withstand compressive forces[END_REF][START_REF] Adams | Abnormal stress concentrations in lumbar intervertebral discs following damage to the vertebral bodies: a cause of disc failure?[END_REF]. As far as we know, these papers are the only one that provide data on lumbar discs loaded till failure upon a displacement-controlled compression. Note that the exact level in the lumbar spine is not specified in these experimental works and important variabilities may exist regarding the disc health state of donor (age, degeneration…) and the testing protocol (loading conditions, environment…). In this regard, the predicted curve may be considered as acceptable since it falls through the experimental points. Nonetheless, when the damage process is introduced in the model prediction, the result underestimates the experimental points as observed in Figure IV.5b. When the OCF damage is excluded from the analysis, the overall disc response first follows the damaged disc response path and then converges progressively towards the healthy disc response path at higher applied displacements. Recall that the constitutive model has been identified at the scale of a unit lamella (stretched along the OCF direction) for the four main disc regions (PO, PI, AO and AI) (Tamoud et al., 2021a). It is clear that the passage from the structural unit to the whole macrostructure needs some adjustement to consider effective interactions between the different lamellae and even between the different networks by adjusting the damage severity to correct its effect on the overall disc stiffness. The latter statement holds true for the

IV.3.3. Strain fields

The predicted radial strains in the core of a compressed healthy disc are compared in For different loading cases of combined axial/torsion loadings, the local strain components can be analyzed in Figures IV.7 and IV.8. Figure IV.7 shows the low strain components that appeared under torsion are extremely amplified when the nucleus is activated by the combined effect of a compression.

Figure IV.8 highlights the very-low radial strains upon pure torsion and important values when the compression is applied. In the latter case, the ILM is found to severly swell whereas the LM shrinks. This phenomenon leads to an auxetic behavior in the radial direction (Tamoud et al., 2021b). This volumetric feature of the lamellar structure is fundamental for the disc functionality for maintaining nucleus swelling and absorbing applied loads. Interestingly, the auxetic response is also influenced by the disc morphology as evidenced in Figure IV.9 in which the contour plots of a simplified disc version are presented. Actually, the heterogeneity of the strain fields is determined by the local structure-mechanics relations but also by the morphology. Figure IV.12 shows the loading mode has a profound effect on the OCF damage. When the disc is compressed the highest OCF damage zones are mainly observed in mid-anterior/lateral anterior ring. The latter can be the sign of the starting areas of annulus circumferential tears. Combined to a compression, the twisted disc presents the highest OCF damage zones for the fibers oriented in twist direction with a progression from the AO side to the PO side. The alternate fibers are only damaged in a small area of the AO side. When the twisted disc is combined to a tension, the damage progression is reversed from the posterior side to the anterior side and, the damage of the alternate fibers is located in the PO region. The NEF network in ILM follows the same reverse compression/tension effect as shown in Figure IV.11. The latter may be due to the nucleus swelling effect but also to the regional difference in disc height. The model can be improved by acting on the constitutive representation to take into consideration the age effects by biological coupling in order to assess the interaction with damage events. Other physiological movements could be introduced in further developments of the disc kinematics in order to propose quantitive predictions under more complex loading modes.

Appendix IV.A. Structure regional variation IV.A.1. OCF regional variation

The regional variation in OCF angle is taken into account using the following linear equation:

AI

 +-are the regional OCF angles.

IV.A.2. NEF regional variation

The volume fraction of a fibril in each ring i is given by: Appendix IV.B. Material properties regional variation

The regional variation of the (regional dependent) material properties is taken into account using the following linear equation:

  -Lagrange strain tensor C : Right Cauchy-green strain tensor B : Left Cauchy-green strain tensor σ : Cauchy stress tensor R : Rotational basis-change tensor Ι : Identity tensor D : Ions diffusivity tensor k : Permeability tensor x , y , z : Cartesian coordinate vectors a : Unit vector of collagen fiber n: Unit vector at the inter-layer interfaces +/-: cationic/anionic c : Ions concentration t : time  : Fiber angle with respect to the xy-plane z : Layer thickness m : Number of layers 

  Figure I.1. Decomposition of the deformation gradient F into chemically-induced fluid transfer part chem F and purely mechanical part mech F .

Figure

  Figure I.2. Chemo-mechanical interactions: (a) Mechanical interactions between layers due to the stress/strain continuity along the interfaces (The main a0-axis of OCF in LM forms an angle θ with (a) (b)

X

  illustrated in Figure I.1. The total deformation gradient F takes the following multiplicative form: are thus considered. The first one, related to the chemical dilatation of the continuum body, is characterized by the chemically-induced volumetric deformation gradient: is the position of the material point in the chemical configuration. The second one, related to the mechanical behavior of the continuum body, is characterized by the mechanical deformation gradient: these transformations, the volumetric change J is defined as:  is the angle between the loading direction and the fibers main axis, see Figure I.2a.

  and OCF free energies, respectively. The intrinsic response of the ECM tangled PG macromolecules is supposed incompressible and fully defined by their stiffness. The ECM free energy i ECM W takes a simple Neo-Hookean form which depends only on the first-invariant:

  of a homogeneous swelling, the internal pressure int p is a function of the concentration of the cationic + and anionic -electronic charges, c + and c -, and the external pressure which R is the gas constant, T is the absolute temperature and, int  and ext  are the internal and external osmotic coefficients, respectively.

  zz D  being the diffusivity in the x, y and z directions, respectively. Chapter I: A microstructure-based model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus 19

  in the x, y and z directions, respectively.It worth noticing that the chemically-induced volumetric stress chem σ of the multi-layered annulus and the swelling pressure given by the classical Donnan theory may be equalized:

  Figure I.3 Annulus specimen extracted from intervertebral disc and stretched in the circumferential direction x with simultaneous local displacement measurements in fibers plane xy and in lamellae plane xz using a right-angle prism.

  Figure I.4. Hybrid experimental/modeling decoupling strategy: (a) from real-time strain field captured in fibers plane xy and in lamellae plane xz using digital image correlation method, the average deformation gradient F is determined in each plane and (b) the deformation of individual (homogenized) layer is extracted using the microstructure-based model. Especially, our strategy allows identifying the volumetric deformation of an individual layer from the specimen volumetric deformation upon different successive loading steps (free swelling/quasi-static stretching/relaxation). The chemo-mechanical response of an individual layer is then deduced. The reconstruction of the multi-layered tissue is finally achieved by considering the chemo-mechanical interactions described in Figure I.2.

Figure

  Figure I.2b. Only ions exchange along the lamellae z direction may be thus considered with the

  Eq. (I.52) are the concentrations of the electronic charges in the left and right border sides of the layer i, respectively, ensuring the continuity conditions of the ions diffusivity between the adjacent layers. The concentrations of the electronic charges in both ends of the multi-layered annulus, i.e. the left border side of the layer 1 and the right border side of the layer m, are equal to the concentration of the external solution *

  and involved in Eqs. (I.37) and (I.38). They are used as direct input data along with the volumetric change (Figure I.5a) to evaluate the spatio-temporal distribution of the anionic concentration through the multi-layered annulus.

  Figure I.6. Anionic concentration through the thickness of the multi-layered soft tissue upon different successive loading steps along with temporal changes: (a) free swelling, (b) quasi-static and (c) relaxation.

exp

  for the stretching and relaxation steps by means of the following objective function obj stress F

I. 3

 3 .5. Model vs. experimental data The model is compared to experimental observations of the time-dependent mechanical behavior of the multi-layered annulus in terms of volumetric change, stress and Poisson's ratio in Figure I.5. A global view at this figure shows the ability of the model to capture the macro-response upon the different successive loading steps. I.4. Model results and discussion The model/experiments comparison in Figure I.5 demonstrates that the evaluation of the overall time-dependent response involves considering volumetric change, Poisson's ratios and stress simultaneously in relation to structural features and diffusion-mechanics. The model can be now used to analyze structural features and diffusion-mechanics effects independently. Chapter I: A microstructure-based model for time-dependent mechanics of multi-

  Figure I.7. Model results: (a) Volumetric change of the LM and ILM layers, (b) Chemical stress of the layered soft tissue, (c) Chemical stress of the LM and ILM layers. It worth noticing this large difference results in a significant change in the volume fraction of the

  Figure I.7b. This feature at the specimen scale is connected to the chemical response at the scale of

  this work, we have presented a physically-based model including realistic structural features (collagen orientation/content, interlamellar thickness and dimensional graduation) in the material representation of the disc annulus along with the electro-chemical activity of the soft tissue. The model successfully reproduced a series of experimental data on bovine tissues under free swelling, quasi-static stretching and relaxation. It was demonstrated that for the evaluation of the overall timedependent response, it is fundamental to consider simultaneously the structural features and diffusion-mechanics effects on volumetric change, Poisson's ratios and stress. The final model is a powerful tool to better understand the annulus mechanics in relation to loading conditions, electrochemical environment, intrinsic properties / proportion / orientation (if any) of constituents (extracellular matrix, fibers, fluid) and structural features.

  (I.A1) into Eq. (I.A2), we obtain:

  I.A6) describes the volumetric change in the LM layers. In the ILM zones devoid of OCF, the volumetric change is simply given by:

  Figure I.A1. Surface of (a) free energy and (b) free energy derivative in the stretch-K2 space for different K1 values.

Figure I. 6

 6 Figure I.6 shows the dimensional graduation taken for this study.

  AbbreviationsLM: Lamellar ILM: Inter-lamellar ECM: Extracellular matrix OCF: Oriented collagen fibers NEF: Nano-sized elastic fibers

Figure

  Figure II.1 presents a schematic representation of the soft tissue hierarchical structure based on

  Figure II.2a. It is expressed in the Cartesian coordinates by

  Figure II.2. Coordinate systems of (a) the multi-layered material element (n0-axis is the unit vector at inter-layer interfaces and a-axis is the unit vector along OCF that forms an angle θ with the xaxis) and (b) of the five NEF bundles  = {0°, 45°, 90°, 135°, 180°} (b-axis is the unit vector along NEF that forms an angle  with the a-axis).

  Figure II.3. Quantity of NEF (a) in LM and (b) in ILM; continuous lines: Eq. (II.3), dashed lines:experimental data(Tavakoli et al., 2017).

  presence of inter-fibrillar connections of type-I collagen fibers in LM zone, a second term appears in the first equation.

F

  is the swelling expansion due to internal fluid content variation and i mech F is the mechanical (isochoric) stress-producing contribution. The sequence of configurations implies also an additive split of the gradient tensor i L of the spatial velocity into a chemical-induced volumetric part i chem L and a mechanical part i mech L as:

  12) Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 55 where i chem W is the chemical-induced volumetric free energy function, due to the internal fluid content variation, expressed in the following form:

  OCF network gives anisotropy, high stiffness and mechanical resistance to annulus. The sketch in Figure II.4 shows the different stretch stages of a stretched fiber in relation to its deformation behavior. In addition to fiber orientation and intrinsic properties, the description of the OCF nonlinear response requires two distinct relations below and above a stretch * OCF  to consider the deformation-induced morphological evolution. The OCF free energy

  Figure II.4. Behavior of a single collagen fiber stretched until failure.

C

  are material constants of the NEF, ( ) H • is again the Heaviside function, Under external mechanical loading, the nature of local failure is the scission of the ECM tangled proteoglycan macromolecules, the detachment of the internal bonds in the discrete (inter-fibrillar and inter-lamellar) fibrils and the breaking of the regularly-oriented fibers. When an individual solid component (ECM, OCF and NEF) is broken, it does no longer contribute to sustain the macro-stress.

   characterize the damage evolution occurring within the multi-layered material element microstructure. This damage representation is in line with the general damage mechanics framework for the hyperelasticity theory. The respective internal stress of each solid component is considered as the quantity controlling the progressive damage evolution from 0.0 for the virgin state to 1.0 for the fully damaged state. Two-parameter Weibull statistical distributions are used to introduce the stochastic nature of the damage process: and  are respectively the Weibull shape parameter and the Weibull scale parameter and, the organization of the collagen network at different scales (as illustrated in Figure II.1) allows to propose a proper description of the progressive degradation considering both localization and anisotropy. The coupling between damage and osmotic effects is considered by introducing the following modification of the chemical-induced volumetric change:

  35) The deformation gradients are detailed in Appendix II.A. The different parts of the Cauchy stress tensor are provided in Appendix II.B. The volume fractions in the initial and current configurations, ratio of the two zones and, a and b are constants calculated according to the maximum and minimum thicknesses of the annulus layers the number of layers. The model makes it possible to estimate the deformation-induced damage and failure of annulus under three-dimensional loading conditions. The multi-layered material element will be subjected to the different mechanical paths illustrated in Figure II.5, namely, uniaxial (UA) stretching in different directions (circumferential, radial and tilted), biaxial stretching and shearing.

Figure II. 5 .

 5 Figure II.5. Loading paths applied on the multi-layered material element.

  41) in which the letters with an overbar denote the experimental data with  the number of considered data points. The optimization of the model parameters was performed using a few well-known available datasets of the swelling response (transversal stretches and volumetric strain) and the UA Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 64 (circumferential and radial) stress responses till moderate stretching and ultimate failure. Before their final optimization using formulae (II.40) and (II.41), a first set of model parameters was

  identified using the nonlinear part (toe zone) of the UA circumferential data ofAcaroglu et al. (1995) and the radial data ofFujita et al. (1997). The stretches, * OCF  and * NEF  , as well as the remaining collagen parameters, 3 C and 6 C , were determined using the transition between the UA nonlinear response and the UA linear response of the Acaroglu et al. (1995) and Fujita et al. (1997) data. The OCF damage parameters, OCF  and OCF , were determined using the UA circumferential data ofAcaroglu et al. (1995) whereas the ECM and NEF damage parameters, determined using the UA radial data ofFujita et al. (1997).

  of stretching are available in the literature. The variation of fluid content by osmotic effect is firstly used to identify the free-swelling response. The identification result is provided in Figure II.6a until chemical equilibrium in a physiological NaCl solution (0.15 M). The volumetric tissue deformation induced by the internal fluid content variation under mechanical loading was then calibrated using Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 65 the experimental data of Baldit (2013) considering UA circumferential stretching and relaxation at constant circumferential stretch under a physiological surrounding environment. The temporal changes of transversal (radial and axial) stretches and the associated circumferential stress are given in Figure II.7. Our model considers the time-dependent chemo-mechanical coupling in the soft tissue to capture the unusual transversal strain history in relation to the actual structure and the chemical-induced internal fluid variation. The strong shrinking in the axial direction and the swelling in the radial direction upon UA circumferential stretching tend towards more usual transversal behavior upon relaxation due to chemo-mechanical equilibrium. The ILM zone is introduced as the key structural parameter governing inter-lamellae fluid exchanges under external mechanical loading and especially the auxetic behavior. This annulus feature was only appreciated

  by purely hyperelastic laws, the present constitutive representation considers the rate-dependency of the annulus mechanics implied by the deformation-induced damage mechanisms and the volumetric effects. Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 66 (a) (b) Figure II.6. Volumetric response of the multi-layered material element: (a) under free swelling (line: model, symbol: experimental data of Acaroglu et al. (1995)), (b) under UA circumferential stretching (solid line: model, dashed line model without damage). The results in Figure II.7 are presented for a moderate stretch range within which the soft tissue stays in a physiological level avoiding any damage effects. The question which arises is now how the tissue behaves under large strains till final failure by including as time-dependent physical process both volumetric effects and damage accumulation. (a) (b) Figure II.7. Temporal changes of (a) transversal stretches and (b) normalized stress of the multilayered material element under a UA circumferential stretching until 1.1 followed by a relaxation (lines: model, symbols: experimental data of Baldit (2013)).Chapter II: A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus 67II.3.2.UA stretching pathThe model efficiencies to capture the stress-stretch behavior of the annulus till failure under UA stretching is presented in Figures II.8 and II.9 for the two main directions in comparison with the experimental data extracted from the works ofAcaroglu et al. (1995) andTavakoli and Costi (2018).The numerical curves show an excellent agreement with the annulus circumferential experimental results ofAcaroglu et al. (1995) in Figure II.8a including the nonlinear stress-stretch response and the failure stage displaying stress drop. Figure II.8c shows that both moduli and ultimate properties

Figure

  Figure II.8. UA circumferential stretching of the multi-layered material element: (a) stress-stretch response (solid line: model; dashed line: model without chemical effect, symbols: experimental data of Acaroglu et al. (1995)), (b) local damage events in the different solid constituents (solid line: model; dashed line: model without chemical effect), (c) toe and linear circumferential moduli, failure strain and failure stress.

Figure

  Figure II.9a shows the UA macro-stress response in the radial direction. A fairly good agreement

  Figure II.12. Biaxial stretching of the multi-layered material element for different biaxiality ratios:(a) circumferential stress-stretch response, (b) radial stress-stretch response.

Figure II. 14 .

 14 Figure II.14. Shearing of the multi-layered material element (continuous line: model, symbols:experimental data ofTavakoli and Costi (2018)).

  stretches in the x, y and z directions, respectively. No fluid transfer occurs under shearing, i.e., i chem = FI . Under UA or biaxial stretching, the chemical deformation gradient is empirically expressed as follows:



  are the chemical expansion rates.

R

  are the basis-change tensor in each layer, with a rotation i  around the z-axis and stretch direction  around the y-axis, respectively:

Figure

  Figure III.1. Two-step strategy consisting in annulus model identification using single lamellae extracted from different disc regions followed by multi-lamellae annulus model predictions for multiaxial elementary loading modes representative of the spine movement.

Figure

  Figure III.2, the disc annulus is a multiple-network medium in which the organization in terms of

  Each individual component is assumed as isotropic, homogeneous, incompressible and hyperelastic body. During the mechanical loading, the initiation and propagation of damage occur through a continuous evolution over the entire loading history. A damage variable d is introduced to provide a mathematical description of the local damage events occurring within the collagen network. Their progressive evolution is given by an expression following a two-parameter Weibull statistical distribution governed by the internal stress: is the Frobenius norm of the effective first Piola-Kirchhoff stress P of the individual solid component (ECM, OCF and NEF).

  fractions at the chemical equilibrium state. The term chem J is the chemical-induced volumetric change.The set of free energies must be prescribed for each component of the solid phase. The simple Neo-Hookean formulation is employed for the ECM free energy exhibits 'toe' and 'linear' regions as described in Figure III.3. The 'toe'

.

  The two angles follow the following geometrical relations to consider stretch-induced re-orientation of the OCF and NEF networks: are the principal stretches, respectively, in the circumferential (x), axial (y) and radial (z) directions, and f  is the stretch in the OCF direction. The re-orientation allows to update the fourth invariant by considering the current fibers orientation in the deformed configuration

Figure

  Figure III.3. Phenomenon of reorientation of collagen fibers.

by

  Skaggs et al. (1994) andHolzapfel et al. (2005) on single annulus lamellae stretched along the fibers main direction. Figure III.4 shows the computed results for the identification of the OCF properties for the different disc regions. (a) (b) Figure III.4. Single lamellae identification for the different disc regions: (a) stress-stretch responses (lines: model simulations, symbols: experimental data of Holzapfel et al. (2005)), (b) ultimate properties (bars: model simulations, symbols with standard deviation intervals: experimental data of Skaggs et al. (1994)).

  1. OCF parameters. This structural difference is reflected by distinct values. The other model parameters come from the previous chapter and are listed in Table

  Loading modes for the multi-lamellae. The multi-lamellae predictions are presented in Figures III.5 and III.6 for UA stretching along the circumferential direction. The simulated stress-stretch curves are compared with the typical experimental curves presented by Ebara et al. (1996) in Figure III.5a. Considering the important variability in the experimental annulus response and that these experimental data are typical results

  Figure III.5b and compared to the failure points with their standard deviation intervals extracted from the works of Acaroglu et al. (1995) and Ebara et al. (1996). (a) (b) Figure III.5. Multi-lamellae UA predictions along circumferential direction for the different disc regions: (a) stress-stretch responses (lines: model simulations, symbols: experimental data of Ebara et al. (1996)), (b) ultimate properties (bars: model simulations, square symbols with standard deviation intervals: experimental data of

  6. Multi-lamellae transverse responses under UA circumferential stretching for the different disc regions (solid lines: fibers plane, dashed lines: lamellae plane): (a) model results, (b) experimental data of Baldit (2013).

  Figure III.7a shows the multi-lamellae stress-stretch responses under UA stretching along the radial direction. Again, the simulations are conducted until the appearance of the stress drop. The experimental data of Fujita et al. (1997) are reported in Figure III.8 to further verify the model predictive capacities. Both moduli and ultimate properties fall within the experimental standard deviation intervals of Fujita et al.

Figure

  Figure III.7. Multi-lamellae stress-stretch responses under (a) UA radial stretching and (b) SS loading for the different disc regions (lines: model simulations, circle symbols: experimental data of Fujita et al. (1997).

Figure

  Figure III.10. Indeed, once PO reaches the circumferential failure, the axial failure occurs. As a

Figure

  Figure III.10. Multi-lamellae stress-stretch responses under PSC loading with various strain ratios for the different disc regions.

  Figure III.11. Full disc model: (a) Passage from annulus constitutive model to anatomically accurate disc model including volumetric strain, regional and directional effects (disc dimensions at the L1-L2 level), (b) boundary conditions of the uniaxially compressed disc.

  Figure III.12 by differentiating the damage of the different solid constituents (ECM, NEF and OCF). That permits a better understanding of the origin of the damage mechanisms affecting the disc and existing at different scales, nanoscale for ECM proteoglycan macromolecules and for NEF network and microscale for OCF. A global view at the contour plots shows that damage progressively increases in intensity and in extent as the applied compression increases, which is accompanied by a disc volumetric change and an annulus-nucleus non-symmetric interaction. The different damage mechanisms may co-exist at the same time and initiate at different locations. Remind that, in the present model, damage initiation and propagation are dependent on structure-property regional variation but also on the variation in displacement fields due to local differences in disc height. The highest damage zones in OCF are mainly observed in the anterior outer/lateral anterior ring.

Figure

  Figure III.12. Damage contour plots in the superior view of the uniaxially compressed disc for different applied displacement levels (0.5, 1.0 and 1.5 mm); ECM damage ECM d

  Figure III.13. Damage of the fibrils in different directions (90° is the radial direction) at two zones of the uniaxially compressed disc: (a) posterior (in LM) and (b) anterior (in ILM) sides of the inner annulus for different applied displacement levels (from zero in the dial center to 2 mm in the dial boundary).

Figure

  Figure III.A2. Model inputs in terms of (a) quantities and (b) orientations for the different disc regions (bars: model inputs, symbols with standard deviation intervals: experimental data).

  ratio of the two zones and, a and b are constants calculated according to the maximum and minimum thicknesses of the annulus layers max z

  and the disc health state. The latter aspects are very attractive in the aim to create patient-specific models taking into account realistic structure/macrostructure features of each patient and disc changes under applied loadings at different levels including disc geometry and local degradation events. The main objective of this work is to formulate a fully three-dimensional model of a complete human intervertebral disc while taking into account the most recent discoveries about the interpenetrating collagen fibrillar networks and the structure-mechanics relationships of the annulus. The disc model is developed within a purely analytical framework in the final goal to Chapter IV: A fully three-dimensional model of interpenetrating collagen fibrillar networks for intervertebral disc mechanics 124 create time and cost-efficient patient-specific models. The predictions are used to predict the heterogeneous mechanics in the intervertebral disc under axial/torsion loading modes. The outline of the present chapter is as follows. Section IV.2 describes the theory considering general assessments of the tube-like solid mechanics along with specifities regarding disc structure and macrostructure. Predictions are presented and discussed in Section IV.3. Section IV.4 closes this chapter with concluding remarks. The following notation is used throughout the text. Tensors and vectors are respectively denoted by normal boldfaced letters and italicized boldfaced letters, while scalars and individual components of vectors and tensors are denoted by normal italicized letters. The superposed dot designates the time derivative.

  inner and outer radii of each ring, respectively, and i H is the initial height. These anatomical features vary from a region to another according to the circumferential axe Θ .

Figure

  Figure IV.1. Axially/torsionally loaded functional spine unit (i.e. an intervertebral disc and two adjacent vertebrae) and related boundary conditions on the disc.

Figure

  Figure IV.2. Disc macrostructure in the polar coordinates.

  inner and outer radii of each ring, respectively, and i h is the current height. The term  is the twist angle of the disc arising from the torsion and i z  is the axial stretch applied in the z-axis. The quantity i f is a function dependent on the initial radius i R of the ring i and the current radius 1 i rof the adjacent ring 1 i -.

F

  the unit vectors in the initial and current configurations, respectively. The chemical-induced volumetric effects are also introduced via a chemo-mechanical coupling considering the multiplicative decomposition concept and the introduction of an intermediate configuration. The total deformation i stretches along the three polar directions r,  and z, respectively. Chapter IV: A fully three-dimensional model of interpenetrating collagen fibrillar networks for intervertebral disc mechanics 127 Using Eqs. (IV.4) to (IV.5), the mechanical deformation gradient tensor i mech may be expressed as:

  The annulus volumetric change due to the internal fluid variation is expressed by the chemical deformation gradient i chem F . Considering fluid diffusivity anisotropy in annulus, the latter is expressed in chapter II as:



  are the chemical expansion rates. The term i  denotes the OCF angle with respect to the circumferential direction of the disc.

  unit vectors at the LM/ILM interfaces in the initial and current configurations, respectively.In the absence of body loads, the equilibrium equations write as the spatial divergence of the overall Cauchy stress tensor σ : , Eq. (IV.13) can be re-written using the stress components in the cylindrical coordinate system as follows: wall of the first ring of the annulus can be expressed dependent on  direction in the form(Holzapfel et al., 2000):As shown in the half-disc topology provided in Figure IV.2, we consider radial and circumferential variations of the geometry in order to propose anatomically accurate patient-specific models. For reasons of symmetry with respect to the sagittal plane of the disc, the nucleus initial radius 0 R was defined by the following special mathematical function by an angle Θ between the anterior region ( 0) Θ = and the posterior region () Θ  = following the direction Θ as illustrated in Figure IV.2:

(

  Figure IV.3. Disc structure decomposition into ECM (along with fluid) and interpenetrating collagen fibrillar networks at the nanoscale (NEF) and at the microscale (OCF). Contour plots in the superior view of the disc show the regional dependent quantities in terms of fluid phase fraction fluid  , collagen fraction

  Figure IV.3 provides the regional dependent quantities in the form of contour plots in the superior view of the disc.



  is the OCF angle with respect to the circumferential direction of the disc in the initial configuration (see Figure IV.4a) and i z  is the axial stretch.



  is the fibrils angle in the initial configuration, i r  is the radial stretch and i f  is the stretch in the OCF direction. The term ( ) H • is the Heaviside function. Organization of the interpenetrating collagen fibrillar networks: (a) OCF angle in the form of contour plots in the superior view of the disc, (b) fibrils angle in the form of contour plots in the sagittal view of the disc. The fibrils spatial distribution is obtained by the linear interpolation of approximated Gaussian-type distributions in the outer and inner boundaries of PO (Posterior-Outer), AO (Anterior-Outer), PI (Posterior-Inner) and AI (Anterior-Inner) regions. The dashed lines in the Gaussian-type distributions are extracted from the work of Tavakoli et al.

  34)in which W  is the OCF free energy, W  is the fibrils network free energy and ECM W energies in Eq. (IV.27) are expressed by considering both collagen fibrillar networks and ECM as hyperelastic bodies with isotropic, homogeneous and incompressible features. The free energies W  and W  may be expressed as a function of the fourth . To describe the 'toe'-to-'linear' transition of the collagen response, two distinct relations are retained below and above * ,

(

H

  • is the Heaviside function that ensures only stretching of the collagen and, ,

  evolution is assumed over the entire loading history without healing possibilities. The progressive damage of each individual solid component (OCF, NEF and ECM) follows a two-parameter Weibull statistical distribution governed by their respective internal stress: norm of the effective first Piola-Kirchhoff stress P.



  IV.40) is affected by the local damage events using the following expression (F of the disc is obtained from the following formula: are the axial, circumferential and radial stress components.

  al. (1996) TableIV.1. Disc macrostructure and structure features.The intrinsic material parameters, identified in the previous chapters, are listed in Tables IV.2 and IV.3. A linear interpolation is here achieved for the regional dependent fibers parameters using the formulation provided in Appendix IV.

  fibrils for which identification has been carried out at the scale of multi-layered annulus stretched along the radial direction(Tamoud et al., 2021c). Whereas the fibrils in LM (inter-fibrillar fibrils) seems to play a minor role on the disc resistance, Figure IV.5a emphasizes the influential and decisive role of the NEF network in ILM (inter-lamellar fibrils) on the disc stiffness especially at low applied displacements. As shown in Figure IV.5b, when the NEF network damage in ILM is not taken into consideration, the overall response is first that of a healthy disc and then it diverges more and more as the applied displacement increases. It can be also observed that the ECM damage has a minor effect on the overall disc resistance. (a) (b) Figure IV.5. Predicted load curves of the axially compressed disc in comparison to experiments (lines: model; symbols: experiments): (a) influence of the LM/ILM collagen networks, (b) influence of the damage events.

  Figure IV.6 with the stereo-radiographic experimental data of Costi et al. (2007). The model results were averaged in different zones and then normalized by the applied displacement. Note that the Costi et al. (2007) data are typical results averaged for the asymmetric values in the right and left sides Chapter IV: A fully three-dimensional model of interpenetrating collagen fibrillar networks for intervertebral disc mechanics 143of the sagittal plane. A weak regional variation is predicted but, except posterior and nucleus zones, the model provides acceptable predictions of the other disc regions.

Figure IV. 6 .

 6 Figure IV.6. Predicted normalized radial strains (%/mm) in a disc compressed at 1.1 mm. The comparison to Costi et al. (2007) data is performed for different regions in which the local strains are averaged (A: Anterior, P: Posterior, AL: Anterior-Lateral, L: Lateral, PL: Posterior-Lateral, NP: nucleus pulposus).

  Figure IV.7. Contour plots of the local axial, circumferential and shear strains z E , E  and z E  in axially/torsionally loaded discs.

Figure

  Figure IV.12. Contour plots of the OCF damage d  (in LM) in axially/torsionally loaded discs.

  respective total number of NEF in each ring i: ) of the outer ring in = and the anterior peak values

  regional material properties.

A microstructure-based model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus

  , 2,...

	i mech i chem FI = == i i i chem x chem y chem z im    F and _ i chem z  are the chemical stretches in the x, y and z directions, (I.43) Assuming mechanical incompressibility, the mechanical deformation gradients under an axial in which _ i chem x  , _ i chem y  respectively. stretching may be expressed as: 11 22 _ _ _ _ diag , , 2, 4,... 1 Chapter I: 25 ( ) 11 22 _ diag , , 1,3,... i iT

Chapter I: A microstructure-based model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus

	( ) ii _ chem j j t  =	1 exp   i j t  --     	i	=	1, 2,... m j x y z , , =	(I.48)
	in which i
	j
			i   0 i chem x S   
	LM i mech F	3 i  2 0 i chem x i i S S S     =+	(I.46)
			12 - 2 i 00 S 
			 
	where
	1 i S	=	13   3 42 2 _ tan tan 4 3 cos 2 ii ii chem x      --     
	21 ii SS  2 _1 1 3 cos i i i chem x S    =+  	2	(I.47)
	32  2 1 cos sin i i i i i SS S   =- 
	The chemical stretches are empirically expressed as follows:
				26

  .

	Using the experimental transversal stretches,	_ exp y 	and	exp z _ 	, and the experimental volumetric
	change,	J , the objective function obj F	to be minimized for the three loading steps is:
		exp			swelling
							2		2	2
		obj	1		_ mod y h	h exp_y		_ mod z h	h exp_z	_ mod chem h	h exp
		swelling			h			h	h
			h	1	exp_y			exp_z	exp

I: A microstructure-based model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus

	( ) ( ) 0 ,, 1 2 1 ,, i i j layer j i i j layer j c z t c i i z t c z t c z t   -= =  -     -  	( sin 2 1 ) 21 n nz ( 21 ) exp i layer i layer z z n    2 i zz Dt   +      +     =- 0 4 n   =    +      	(I.51)
	in which
				ii
	c	i	z	j	11 , j at t z z
				layer				layer
				22 jj
				ii
	c	i	z	j	,	t	at	z	z	j
				layer				layer
				11 jj
							27

Table I

 I

	.2 gives the identified values with a goodness of fit R 2 = 0.986.
				Parameters	Values
	ECM G	  KPa	4 (0.37) *
	1 C	  MPa	3.21 (0.02) *
	C	2	  -	11.40 (2) *
	1 K	  KPa	61.40 (0.24) *
	K	2	  -	0.12 (0.02) *
				Table I.2. Parameters of Eqs. (I.30)-(I.32).
			* Between brackets is the standard deviation.

  4...

			mm	-
	F	i layer ii LM i  layer i == =+  F	F	ILM i	(II.32)
	and
			mm	1 -
		σ	__ 2,4,.. LM i i ILM i layer ii 1,3,.. i layer  == =+  σ σ	(II.33)

Table II

 II

	.1.
	Parameters	Values	References
		z	max	0.76	  mm	Holzapfel et al. (2005)
		z	min	0.69	  mm	Holzapfel et al. (2005)
	ILM LM r	0.133	  -	Tavakoli et al. (2017)
				30	  deg	Holzapfel et al. (2005), Baldit 2013
	collagen 	0.07	  -	Skaggs et al. (1994)
		fluid	_0	0.7	  -	Acaroglu et al. (1995)
						Table II.1. Structural parameters.

Table II

 II

	F	obj	and obj F	to be minimized:
	swelling	stress
				obj swelling F	=	2 hh hh hh yy zz yz     2 1 1 hh chem chem h h chem   =   -2   -- + +          	(II.40)
						obj	1		circ P h	circ P h	22 rad rad P P h h
					stress F	h	1	hh circ rad PP

.2. gives the identified values of the intrinsic mechanical parameters of the solid components (ECM, OCF and NEF) and the swelling response. Their identification was solved numerically as an optimization problem using MATLAB software considering the two following objective functions

multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus

	B1) (II.B2) NEF  is given by the sum of the Cauchy stress tensor i _ _ ILM i i ILM i mech ECM NEF =+ σ σ σ  σ of all Chapter II: A 81 in LM layer and, in ILM layer. The NEF Cauchy stress tensor i NEF σ discrete fibrils as follows: _ _ LM i LM i NEF NEF    σσ and __ =  45,90,135

Table III .

 III .B1. B1. Structural parameters for the different disc regions.Acaroglu, E.R., Iatridis, J.C., Setton, L.A., Foster, R.J.,Mow, V.C., Weidenbaum, M., 1995. Degeneration and aging affect the tensile behavior of human lumbar anulus fibrosus. Spine 20, 2690-2701. Amonoo-Kuofi, H.S., 1991. Morphometric changes in the heights and anteroposterior diameters of the lumbar intervertebral discs with age. Journal of Anatomy 175, 159-168. Baldit, A., 2013. Etude des interactions hydro-chimio-mécaniques dans les tissus biologiques: application à la nutrition du disque intervertébral. Phd Thesis, University of Montpellier. Bloom, E.T., Lee, A.H., Elliott, D.M., 2021. Tendon multiscale structure, mechanics, and damage are affected by osmolarity of bath solution. Annals of Biomedical Engineering 49, 1058-1068.

		Parameters		Regions		References
			AO	PO	AI	PI
		  mm	0.76	0.42	0.86	0.46	Holzapfel et al. (2005)
		  mm	0.62	0.34	0.66	0.34	Holzapfel et al. (2005)
		  -	0.133 0.133 0.133 0.133	Tavakoli et al. (2017)
	III. References
	max min ILM LM
	r	z z

  . Classification of agerelated changes in lumbar intervertebral discs. Spine 27, 2631-2644. Borde, B., Grunert, P., Härtl, R., Bonassar, L.J., 2015. Injectable, high-density collagen gels for annulus fibrosus repair: an in vitro rat tail model. Journal of Biomedical Materials Research Part A 103, 2571-2581. Bruehlmann, S.B., Hulme, P.A., Duncan, N.A., 2004. In situ intercellular mechanics of the bovine outer annulus fibrosus subjected to biaxial strains. Journal of Biomechanics 37, 223-231.

fully three-dimensional model of interpenetrating collagen fibrillar networks for intervertebral disc mechanics

	LM i t r ILM LM i t a i b t =+ LM t = represents the thickness ratio of the ILM and LM zones and, for 1,3,5... in = ILM i t for 2, 4, 6... 1 in =-( ) ( ) ( ) 1 PO PI AI AO AO AI t t t t t t t a n   -+ -+ -= -( ) PI AI AI t t t t b a   -+ = -t , PI ILM LM r expressed as: where in which PO t , AO t and AI t are the regional ring thicknesses. The terms Chapter IV: A 131 thicknesses for the outer side of the ring in = , respectively, in the anterior region ( t a and b are functions (IV.21) (IV.22) t (IV.23) (IV.24) AO t and PO t denote the 0)

fully three-dimensional model of interpenetrating collagen fibrillar networks for intervertebral disc mechanics

  .1.

	Regions Chapter IV: A 140 AO AI Layer thickness [mm] AO t 0.69 AI t 0.76 Disc height [mm] AO H AI H 12.8 8.40	PO PO t 0.38 PO H 8.00	PI PI t 0.40	References Holzapfel et al. (2005)

fully three-dimensional model of interpenetrating collagen fibrillar networks for intervertebral disc mechanics 150

	Chapter IV: A
					a	, LM ILM _ 	=	, LM ILM _ PO  		+	, LM ILM _ AI  	-	, LM ILM _ AO  	, LM ILM _ 1 PI   n  --	, LM ILM _ AO    +	-	, LM ILM _ AI  		(IV.A10)
												( ,, ) , , , _ _ _ __ LM ILM LM ILM LM ILM PI AI AI LM ILM LM ILM ba              -+ =-	(IV.A11)
	The terms	, LM ILM _ a 	and	, LM ILM _ b 	are expressed as a function of the posterior and anterior scale factors
		, LM ILM _ PO 	,		, LM ILM _ PI 	,		, LM ILM _ AO 	and		, LM ILM _ AI 	:
					, LM ILM _  a	=	(		, LM ILM _ PO 	+		, LM ILM _ AI 	-		, LM ILM _ AO 	( ) , _ 1 LM ILM PI n    --	) (   +	, LM ILM _ AO 	-		, LM ILM _ AI 	)		(IV.A12)
													( ,, ) , , , _ _ _ __ LM ILM LM ILM LM ILM PI AI AI LM ILM LM ILM ba              -+ =-	(IV.A13)

This chapter is based on the following paper:Tamoud, A., Zaïri, F., Mesbah, A., Zaïri, F.,

A microstructurebased model for time-dependent mechanics of multi-layered soft tissues and its application to intervertebral disc annulus. Meccanica 56, 585-606.

This loading path is representative of the compression mechanics of disc. Indeed, the vertebral column is permanently subjected to an axial compression (due to muscular tension and body weight) prior to further

3D loadings that expose the annulus to a tensile stretching in the disc circumferential direction due to the swelling of the "gelatinous" nucleus pulposus retained by the adjoining concentric annulus lamellae[START_REF] Nachemson | In vivo measurements of intradiscal pressure: discometry, a method for the determination of pressure in the lower lumbar discs[END_REF].

ECM represents all non-fibrillar "solid" components and is mainly constituted by a network of randomly oriented and negatively charged proteoglycan (PG) macromolecules, conferring hydrophilic properties to the disc soft tissues(Urban and Maroudas, 1981), randomly oriented type-II collagen and other types of collagen. In the special case of the disc, the proportions of the latter constituents vary from a disc region to another[START_REF] Roberts | Biochemical and structural properties of the cartilage end-plate and its relation to the intervertebral disc[END_REF].

The experimental results have been analyzed statistically taking into account: (i) the errors related to the inter-specimen variability using fifteen samples tested under the same physiological and mechanical conditions by the analysis of variance (ANOVA) method with a variability factor of p < 0.05, and (ii) the errors related to the correlation procedure with an average error of 9.7%.

It is termed apparent Poisson's ratio since it is dependent on the fluid amount and it is not intrinsic to the tissue.

Age-related changes include decreased PG macromolecules content and a corresponding reduction in chemicallyinduced volumetric stress leading to a decrease in annulus stiffness[START_REF] Eyre | Biochemistry of the intervertebral disc[END_REF][START_REF] Buckwalter | Aging and degeneration of the human intervertebral disc[END_REF][START_REF] Roberts | Histology and pathology of the human intervertebral disc[END_REF].

This chapter is based on the following paper:Tamoud, A., Zaïri, F., Mesbah, A., Zaïri, F., 2021. A multiscale and multiaxial model for anisotropic damage and failure of human annulus fibrosus. International Journal of Mechanical Sciences, 106558.

This chapter is based on the following paper: Tamoud, A., Zaïri, F., Mesbah, A., Zaïri, F., 2021. Modeling multiaxial damage regional variation in human annulus fibrosus. ActaBiomaterialia 136,[375][376][377][378][379][380][381][382][383][384][385][386][387][388]

This chapter is based on the following paper: Tamoud, A., Zaïri, F., Mesbah, A., Zaïri, F., 2021. A fully threedimensional model of interpenetrating collagen fibrillar networks for intervertebral disc mechanics. Submitted.

Acknowledgement

Cette thèse s'inscrit dans le cadre d'une

self-healing response of the disc tissues (Likhitpanichkul et al., 2014; Borde et al., 2015;Wan et al., 2016;Peng et al., 2020). Moreover, the near-field direct local interactions within the collagen network and the ECM ought to be treated through either the incorporation of a supplementary energy term (Guerin and Elliott, 2007) or a mean-field approach (Saadedine et al., 2021).

Chapter III

86

Modeling multiaxial damage regional variation in human annulus fibrosus

Proper description of damage and failure requires to consider the latter that can be introduced in continuum-based models. [START_REF] Mengoni | Modelling the failure precursor mechanism of lamellar fibrous tissues, example of the annulus fibrosus[END_REF] developed a damage model for the non-fibrillar matrix restricted to the anterior annulus stretched in the radial direction. Later, [START_REF] Gao | An anisotropic multiphysics damage model with application to annulus fibrosus[END_REF] and Ghezelbash et al. (2020) proposed a damage model considering only the collagen fibers failure as local damage event. Tamoud et al. (2021b) developed a damage model considering the multiscale failure due to network hierarchical organization. The chemical-induced volumetric effects on the structure-damage relation were also introduced to consider the hydration state effects. Nonetheless, our previous work in (Tamoud et al., 2021b) disregards the regional variation while it is the prerequisite condition to any disc prediction. The directional effects due to the structural arrangement along with the regional dependence make damage prediction a challenging task that still remains largely an open issue. To date, no constitutive model has been proposed to predict the multiaxial mechanics of the annulus fibrosus while the regional dependence of the relation between structure and damage is taken into consideration.

The main objective of this work is to formulate a fully three-dimensional constitutive representation of the annulus fibrosus to capture the regional dependence of the structure-damage relation. In order to reproduce the multiaxial loadings as close as to those undergone in-vivo by the annulus fibrosus, the material representation considers the complex collagen network arrangement of each disc region. To consider the progressive deformation-induced changes in the tissue, the model introduces the stretch-induced network re-orientation as well as the network damage leading to final failure.

The model is identified at the scale of single annulus lamellae stretched along the fibers main direction. The regional dependent model predictions are then analyzed under various loading modes at the scale of multi-lamellae and compared to experimental data when available. The human annulus constitutive model is finally used to predict the heterogeneous mechanics of a full disc model.

III.2. Model formulation

We propose here a two-step strategy consisting in single lamellae identification followed by multilamellae multiaxial predictions as illustrated in Figure III.1. The required model parameters are identified on axial tests performed along the micro-sized collagen fibers of single annulus lamellae extracted from the different disc regions (AO: Anterior-Outer, AI: Anterior-Inner, PI: Posterior-Inner, PO: Posterior-Outer). Beyond taking into account the regional effect, this identification strategy is a valuable difference with our previous work (Tamoud et al., 2021b) in which the model identification had been performed at the scale of multi-lamellae. In this chapter, after its identification at the scale of single lamellae, the model is then verified by comparing its prediction with the experimental data of multi-lamellae performed under multiaxial loading. Uniaxial (UA), simple shear (SS), biaxial (BA) and plane-strain compression (PSC) tests are all being considered.

Appendix III.A. Quantities and orientations

The total number of NEF is reported in Figure III.A1 as a function of the orientation  . The information is extracted from the two works of Tavakoli et al. (2017Tavakoli et al. ( , 2020a)); the observed Gaussiantype shape is mathematically represented by a special function proposed in the previous chapter.

Due to the lack of information regarding the regional variation, we admitted similar NEF distributions in LM and ILM layers for outer and inner regions.

Appendix III.B. Multi-lamellae mechanics

By separating chemical-induced volumetric and mechanical quantities, the macroscopic deformation gradient tensor F is given from chapter I:

General conclusions and perspectives

In Other physiological movements could be introduced in further developments of the disc kinematics in order to propose quantitive predictions under more complex loading modes. The model could be improved by acting on the constitutive representation to take into consideration the age effects by biological coupling in order to assess the interaction with damage events. The model could be indeed used to better understand the tissue response changes during aging and degeneration in response to changes in the tissue composition. The present model can be thus considered as the stepstone for the long-term prediction of the disc dysfunctions. During different complex loading scenarios of the spine, the disc dysfunctions may occur either by: i) mechanical damage due to a high level of loads or relatively moderate loads repeated several times which cause damage accumulation, ii) by an age-related biological degeneration or by the cellular disorders, including the cell death, the production of inflammatory mediators and a shift towards catabolism. Nonetheless, it is difficult to identify or precisely model the disc failure taking into account the two causes of dysfunction at the same time. Since our study is based on the in-vitro annulus results, it cannot be directly extrapolated to the in-vivo disc response. Indeed, the understanding of the disc mechanobiology must be taken into account by the model considering especially nutrition and cellular activity along with degenerative aspects. The model could be then a real help to provide solutions for remodeling and self-healing response of the disc tissues. Moreover, the near-field direct local interactions within the collagen network and the ECM ought to be treated through either the incorporation of a supplementary energy term or a mean-field approach.

Abstract

The damage in annulus soft tissues is a complex multiscale phenomenon due to a complex structural