Seungik Baek

Ataollah Ghavamian
email: a.ghavamian2@gmail.coms.jamaleddin

Stéphane Avril
email: avril@emse.fr

S Jamaleddin Mousavi

Computational Study of Growth and Remodeling in Ascending Thoracic Aortic Aneurysms Considering Variations of Smooth Muscle Cell Basal Tone

Keywords: finite element, constrained mixture theory, growth and remodeling, smooth muscle cells, active stresses

In this paper, we investigate the progression of Ascending Thoracic Aortic Aneurysms (ATAA) using a computational model of Growth and Remodeling (G&R) taking into account the composite (elastin, four collagen fiber families and Smooth Muscle Cells-SMCs) and multi-layered (media and adventitia) nature of the aorta. The G&R model, which is based on the homogenized Constrained Mixture theory, is implemented as a UMAT in the Abaqus finite-element package. Each component of the mixture is assigned a strain energy density function: nearly-incompressible neo-Hookean for elastin and Fung-type for collagen and SMCs. Active SMCs tension is additionally considered, through a length-tension relationship having a classic inverted parabola shape, in order to investigate its effects on the progression of ATAA in a patient-specific model. A sensitivity analysis is performed to evaluate the potential impact of variations in the parameters of the length-tension relationships. These variations reflect in variations of SMCs normal tone during ATAA progression, with active stress contributions ranging between 30% (best case scenario) and 0% (worst case scenario) of the total wall circumferential stress. Low SMCs active stress in the worst case scenarios, in fact, affect the rates of collagen deposition by which the elastin loss is gradually compensated by collagen deposition in the simulated ATAA progression, resulting eventually in larger aneurysm diameters. The types of length-tension relationships leading to a drop of SMCs active stress in our simulations reveal a critical condition which could also result in SMCs apoptosis.

INTRODUCTION

Nowadays, the paramount importance of mechanobiology is widely acknowledged in morphogenesis and pathogenesis [START_REF] Humphrey | Mechanisms of arterial remodeling in hypertension: coupled roles of wall shear and intramural stress[END_REF]. The adaptation ability of soft tissues relies on the existence of a preferred load-bearing mechanical state, the so-called homeostasis, across multiple length/time scales. At the tissue scale, this is exhibited through continuous mass changes of the components of the Extracellular Matrix (ECM), such as elastin, collagen, and proteoglycans in blood vessels [START_REF] Humphrey | A constrained mixture model for growth and remodeling of soft tissues[END_REF]Cyron et al., 2016). Mechanobiology implies that Growth and Remodeling (G&R) of the biological tissues is mediated by mechanical stresses [START_REF] Grossman | Cardiac hypertrophy: useful adaptation or pathologic process?[END_REF][START_REF] Humphrey | Mechanisms of arterial remodeling in hypertension: coupled roles of wall shear and intramural stress[END_REF].

In order to numerically simulate G&R in arteries, the constrained mixture model (CMM) was developed by [START_REF] Humphrey | A constrained mixture model for growth and remodeling of soft tissues[END_REF] and, thereafter, has been increasingly utilized [START_REF] Watton | A mathematical model for the growth of the abdominal aortic aneurysm[END_REF][START_REF] Baek | A theoretical model of enlarging intracranial fusiform aneurysms[END_REF][START_REF] Zeinali-Davarani | Medical image-based simulation of abdominal aortic aneurysm growth[END_REF][START_REF] Valentín | A finite element-based constrained mixture implementation for arterial growth, remodeling, and adaptation: theory and numerical verification[END_REF]Cyron et al., 2016;[START_REF] Braeu | Homogenized constrained mixture models for anisotropic volumetric growth and remodeling[END_REF][START_REF] Lin | Growth and remodeling with application to abdominal aortic aneurysms[END_REF][START_REF] Famaey | Numerical simulation of arterial remodeling in pulmonary autografts[END_REF]Latorre and Humphrey, 2018b). [START_REF] Baek | A theoretical model of enlarging intracranial fusiform aneurysms[END_REF] proposed a two-dimensional CMM for arterial G&R. A three-dimensional representative straight cylindrical artery was also introduced by [START_REF] Karsaj | A 3-d framework for arterial growth and remodeling in response to altered hemodynamics[END_REF] with evolving geometry, structure, and mechanical properties under changes in hemodynamics (i.e., mean blood flow/pressure). [START_REF] Valentín | A finite element-based constrained mixture implementation for arterial growth, remodeling, and adaptation: theory and numerical verification[END_REF] established a non-linear model using the Finite Element Method (FEM) on the basis of the Constrained Mixture Theory (CMT) aimed at facilitating numerical analyses of arterial adaptation and maladaptation. They could predict changes in fiber orientations and quantities, degradation of elastin and loss of Smooth Muscle Cells (SMCs), including disease progression and evolving geometries. [START_REF] Famaey | Numerical simulation of arterial remodeling in pulmonary autografts[END_REF] implemented the same model in Abaqus [START_REF] Hibbit | Abaqus-Theory Manual[END_REF] and predicted adaptation of a pulmonary autograft. Cyron et al. (2016), and then [START_REF] Braeu | Homogenized constrained mixture models for anisotropic volumetric growth and remodeling[END_REF], introduced the homogenized CMM framework for G&R using an informal temporal averaging approach, bringing simplicity and computational efficiency. [START_REF] Lin | Growth and remodeling with application to abdominal aortic aneurysms[END_REF] combined homogenization and the CMT to simulate the dilatation of abdominal aortic aneurysms. Their methodology could capture important aspects, such as mass turnover in arterial walls with a low computational cost. Recently, Latorre and Humphrey (2018a) introduced a new rate-based CMM formulation suitable for studying mechanobiological equilibrium and stability of soft tissues exposed to transient or sustained changes, permitting direct resolution of G&R problems with a quasi-static approach.

Despite the significant insight of the CMM methodology in arterial G&R, it has generally been employed in the case of canonical problems [START_REF] Baek | A theoretical model of enlarging intracranial fusiform aneurysms[END_REF]Cyron et al., 2016;[START_REF] Laubrie | A new finite-element shell model for arterial growth and remodeling after stent implantation[END_REF] or single-layer thick-wall axisymmetric [START_REF] Braeu | Homogenized constrained mixture models for anisotropic volumetric growth and remodeling[END_REF][START_REF] Lin | Growth and remodeling with application to abdominal aortic aneurysms[END_REF] approximations. The extension to ATAAs was challenging due to simultaneous and regionspecific evolutions of geometry, material properties [START_REF] Farzaneh | Identifying local arterial stiffness to assess the risk of rupture of ascending thoracic aortic aneurysms[END_REF], and hemodynamic loads [START_REF] Humphrey | Mechanics, mechanobiology, and modeling of human abdominal aorta and aneurysms[END_REF][START_REF] Condemi | Fluid-and biomechanical analysis of ascending thoracic aorta aneurysm with concomitant aortic insufficiency[END_REF]. Only recently, [START_REF] Mousavi | Patient-specific predictions of aneurysm growth and remodeling in the ascending thoracic aorta using the homogenized constrained mixture model[END_REF] proposed a non-linear FEM solution based on the homogenized CMT to simulate G&R in a layer-specific and patient-specific ATAA. They showed that elastin loss leads to a transfer of stress to the adventitia and continuous adaptation of the stress distribution through changes in ATAA shape.

Nevertheless, the significant role of SMCs in regularizing the wall tension of hollow organs should not be neglected. This occurs through contraction and relaxation of SMCs which is regulated by phosphorylation of the myosin motors associated with the smooth muscle contractile units. SMCs also help to control the arterial wall stiffness (e.g., during a cardiac cycle) and, then, regularizing the blood pressure (Murtada et al., 2010b). It is reported by [START_REF] Murtada | Investigating the role of smooth muscle cells in large elastic arteries: a finite element analysis[END_REF] that the large elastic deformations of arteries undergoing a physiological loading cannot rely only on a passive response. The active contribution (contractility) of SMCs can generate significant tension. [START_REF] Rachev | Theoretical study of the effects of vascular smooth muscle contraction on strain and stress distributions in arteries[END_REF] proposed a phenomenological model for SMCs in the arterial wall with a length-tension relationship depending on the circumferential stretch. Similarly, [START_REF] Schmitz | On a phenomenological model for active smooth muscle contraction[END_REF] proposed a steadystate model of SMCs activation containing a phenomenological explanation of the active length-tension behavior. [START_REF] Zulliger | A constitutive formulation of arterial mechanics including vascular smooth muscle tone[END_REF] extended this method, but this time for a vascular SMCs. They investigated the pressure-radius curve for three different SMC states (i.e., fully relaxed, maximally contracted, normal tone) under physiological and pathological conditions with varying levels of SMCs tone. [START_REF] Standley | Cyclic stretch induces vascular smooth muscle cell alignment via NO signaling[END_REF] proposed a thermodynamically consistent constitutive model as a function of stretches, a variable accounting for the SMC state (free unphosphorylated myosin, phosphorylated crossbridges, phosphorylated and dephosphorylated cross-bridges attached to actin), a variable related to the Ca 2+ concentration, and temperature. They illustrated how SMCs contraction and related stretch are co-dependent. This work inspired Murtada et al. (2010a) to propose a structural description of the SMCs contractile unit including elastic elongation of the cross-bridges [START_REF] Hai | Cross-bridge phosphorylation and regulation of latch state in smooth muscle[END_REF][START_REF] Yang | The myogenic response in isolated rat cerebrovascular arteries: smooth muscle cell model[END_REF] and myofilament sliding. They showed that the driving stress, in general, depends on the state of the muscle and, during contraction, depends only on the attached cycling cross-bridges. The model was extended (Murtada et al., 2010a[START_REF] Murtada | Experiments and mechanochemical modeling of smooth muscle contraction: significance of filament overlap[END_REF] to capture the active length-tension relationship, the evolution law of filament sliding and the isotonic force-shortening velocity relationship. [START_REF] Murtada | Investigating the role of smooth muscle cells in large elastic arteries: a finite element analysis[END_REF] studied the role of SMCs in large elastic arteries using a FEM simulation and demonstrated that, yet, changes in the intracellular calcium amplitudes barely affect the circumferential stress. An increase in the mean intracellular calcium value and, then, in the medial wall thickness, clearly results in changes in stress distribution and the overall deformation. [START_REF] Murtada | Adaptation of active tone in the mouse descending thoracic aorta under acute changes in loading[END_REF] investigated the active tone of SMCs in the specific case of murine descending thoracic aortas undergoing acute loading changes and showed the variation in active SMCs tone depending on different induced stretches. They realized that the active tone was reduced when the artery adapted below the optimal stretch with no significant change in passive behavior.

Although the aforementioned accomplished contributions have brought a significant insight on the active role of SMCs into soft tissues undergoing physiological loading, there is still a pressing need to have computational models investigating the possible variations of SMCs contractility in realistic multi-layer patient-specific geometries with irregular boundary conditions and complex deformations. In this work, the G&R model developed by [START_REF] Mousavi | Patient-specific predictions of aneurysm growth and remodeling in the ascending thoracic aorta using the homogenized constrained mixture model[END_REF] is extended to include the variations of SMCs contractility and a sensitivity analysis is performed to evaluate their impact on ATAA progression in a patient-specific model.

MATERIALS AND METHODS

Mathematical Model

Kinematics

We consider the three dimensional deformation of a continuum (e.g., in vivo healthy configuration of a blood vessel) moving from its initial undeformed configuration, say before any G&R happens, occupying a volume 0 , of boundary ∂ 0 , to a time dependent deformed configuration occupying a volume (t), of boundary ∂ (t) at time t. The motion of the continuum is defined through a deformation mapping, such as x = φ(X, t) by which a solid particle at the reference state is transformed to the spatial configuration. The total deformation gradient F of a mixture, consisting of different constituents (e.g., elastin, collagen fiber families, and SMCs), can be expressed as

F(X, t) = ∂φ(X, t) ∂X . ( 1
)
The deformation gradient essentially relates a fiber at material configuration to its spatial counterpart (dx = FdX).

Additionally, the volume map J is defined as

J = det(F), (2)
linking a reference volume element into the deformed states (d = Jd 0 , J > 0). Two main assumptions are made for the deformation gradient on the basis of CMT, namely, (1) all constituents deform together resulting in a unified deformation gradient F and (2) the deformation gradient of each constituent i of the mixture (i ∈ {e, c j , m}, where e stands for elastin, c j stands for collagen family j where j ∈ {1, 2, 3, 4} and m stands for SMCs) is decomposed into an elastic F i el and an inelastic F i gr part as

F = F i el F i gr = F i el F i r F i g , (3)
where the elastic deformation gradient F i el is responsible for the generation of the stress field and the inelastic deformation gradient F i gr considers the differential mass turnover in F i g and the changes in the microstructure in F i r . Initially F(t 0 ) = I where I is the identity tensor, each constituent i of the mixture has an individual elastic deformation gradient F i el (t 0 ) corresponding to its deposition stretch [START_REF] Mousavi | Patient-specific stress analyses in the ascending thoracic aorta using a finite-element implementation of the constrained mixture theory[END_REF] in which G i h are second order tensors defined, such as

F i el (t 0 ) = G i h = F i r (t 0 ) -1
G e h = diag 1 λ e 0θ λ e 0z
, λ e 0θ , λ e 0z , (4a)

G i h = λ i 0 a i 0 ⊗ a i 0 + 1 λ i 0 (I -a i 0 ⊗ a i 0 ); i ∈ {c j , m}, (4b)
where λ e 0θ and λ e 0z represent the initial (in the undeformed configuration before any G&R) deposition stretches of elastin in the circumferential and longitudinal directions, respectively, and λ i 0 stands for the initial deposition stretch of constituent i in the fiber direction whose unit vector is denoted by a i 0 .

Balance of Linear Momentum

The G&R deformation process is typically governed by the balance of the linear momentum equation

DIV(P) + ̺ 0 b 0 = 0, (5)
in which DIV is the material divergence operator, P denotes the first Piola-Kirchhoff stress tensor, b 0 is the body forces per unit mass and ̺ 0 = i ̺ i 0 represents the mass density (per unit undeformed volume) of the mixture, which is equivalent to the sum of mass densities (per unit undeformed volume) of all constituents. Note here that the dynamic effects (i.e., inertia or viscoelasticity) are neglected since the G&R progression occurs at very slow time scales (days to months) [START_REF] Braeu | Homogenized constrained mixture models for anisotropic volumetric growth and remodeling[END_REF].

Mechanobiological Constitutive Model

A strain energy density function W (per unit undeformed volume) is introduced [START_REF] Braeu | Homogenized constrained mixture models for anisotropic volumetric growth and remodeling[END_REF][START_REF] Mousavi | Patient-specific stress analyses in the ascending thoracic aorta using a finite-element implementation of the constrained mixture theory[END_REF])

W = ̺ e 0 (t)W e elastin + n j=1 ̺ c j 0 (t)W c j
Collagen fiber families

+ ̺ m 0 (t)W m SMCs . ( 6
)
In the above equation, ̺ i 0 (t) and W i (i ∈ e, c j , m ) represent, respectively, the mass densities and strain energy densities (per unit mass) of each individual constituent.

A nearly incompressible Neo-Hookean material model is introduced for elastin as [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF][START_REF] Mousavi | Patient-specific stress analyses in the ascending thoracic aorta using a finite-element implementation of the constrained mixture theory[END_REF])

W e = W e (I e 1 , J) = µ e 2 (I e 1 -3) W e Deviatoric + κ 2 (J e el -1) 2 W e Volumetric , (7)
where µ e and κ stand for the shear modulus and the bulk modulus (stress per unit mass dimensions) of the elastin, respectively. J e el denotes the elastic contribution of the Jacobian for elastin and I i 1 is the isochoric first invariant of the elastic Cauchy-Green tensor, defined, such as

I e 1 = tr(C e el ), (8)
where

C e el = F e el
T F e el is the elastic part of the isochoric right Cauchy-Green tensor for elastin. The isochoric deformation gradient is defined as F e el = J e el -1/3 F e el . The strain energy density of collagen families is described with a Fung-type exponential expression [START_REF] Mousavi | Patient-specific stress analyses in the ascending thoracic aorta using a finite-element implementation of the constrained mixture theory[END_REF])

W c j (I c j 4 ) = k c j 1 2k c j 2 exp k c j 2 (I c j 4 -1) 2 -1 , (9)
where k c j 1 and k c j 2 are stress per unit mass and dimensionless material parameters, respectively, taking different values depending whether the fibers are under compression or tension [START_REF] Bersi | Novel methodology for characterizing regional variations in the material properties of murine aortas[END_REF][START_REF] Mousavi | Patient-specific stress analyses in the ascending thoracic aorta using a finite-element implementation of the constrained mixture theory[END_REF], and I c j 4 is the fourth invariant which may be written, such as

I c j 4 = C c j el :(a c j ⊗ a c j ) = (λ c j el ) 2 ; a c j = F c j gr a c j 0 F c j gr a c j 0 , ( 10
)
where C c j el is the elastic part of the right Cauchy-Green stretch tensor, a c j is the unit vector along the fiber direction in the inelastically deformed intermediate configuration and λ c j el is the elastic stretch of collagen fibers.

The passive and active strain energy density of SMCs is also introduced as [START_REF] Wilson | Parametric study of effects of collagen turnover on the natural history of abdominal aortic aneurysms[END_REF][START_REF] Murtada | Investigating the role of smooth muscle cells in large elastic arteries: a finite element analysis[END_REF][START_REF] Murtada | Adaptation of active tone in the mouse descending thoracic aorta under acute changes in loading[END_REF])

W m (I m 4 , λ act ) = k m 1 2k m 2 exp k m 2 (I m 4 -1) 2 -1 W m pass + σ actmax ̺ 0 λ act + 1 3 (λ max -λ act ) 3 (λ max -λ 0 ) 2 W m act , (11)
where k m 1 and k m 2 are similar to the k c j 1 and k c j 2 parameters, I m 4 is also identical to I c j 4 defined in Equation ( 10) (replacing c j by m) and λ act = λ m λ m 0 is the active stretch in the fiber direction, where λ m = Fa m 0 = C : a m 0 ⊗ a m 0 . Moreover, σ actmax denotes the maximal active Cauchy stress and λ max and λ 0 are parameters denoting the active stretches at maximum and zero active stresses, respectively.

It is worthwhile noting that the same strain energy density functions are assumed all across the aorta. However, different material properties and mass densities of the individual constituents are assigned to each layer (media and adventitia). Note that the intima layer is disregarded in this work, due to its relatively thin thickness.

Stress Evaluation

Having Equation (6), it is now possible to evaluate the second Piola-Kirchhoff stress tensor using the following relations

S = 2 ∂W ∂C = 2 i ̺ i 0 ∂W i ∂C = i ϕ i S i (12)
where S i is the Piola-Kirchhoff stress of constituent i and ϕ i its volume fraction in the current configuration, which may be written, such as ϕ e S e = ̺ e 0 (t) (J e el ) -2/3 µ e (C e gr ) -1 -

1 3 tr C e el C -1 + κJ e el (J e el -1)C -1 ; (13a) ϕ c j S c j = ̺ c j 0 (t) 2 k c j 1 (I c j 4 -1) exp k c j 2 (I c j 4 -1) 2 F c j gr a c j 0 2 a c j 0 ⊗ a c j 0 ; (13b)
ϕ m S m = ̺ m 0 (t) 2 k m 1 (I m 4 -1) exp k m 2 (I m 4 -1) 2 F m gr a m 0 2 a m 0 ⊗ a m 0 (13c) + ̺ m 0 (t) ̺ 0 (0) σ actmax C : a m 0 ⊗ a m 0 1 - (λ max -λ act ) 2 (λ max -λ 0 ) 2 a m 0 ⊗ a m 0 .
The obtained second Piola-Kirchhoff stress can simply be pushed forward in order to compute the Cauchy stress tensor

σ i = 1 J FS i F. (14)
Therefore, the Cauchy stress for each individual constituent can be expressed as

ϕ e σ e = ̺ e (t)(J e el ) -2/3 µ e B e el - 1 3 tr C e el I
+ ̺ e 0 (t)(J e gr ) -1 κ(J e el -1)I;

(15a)

ϕ cj σ cj = 2̺ cj (t) k cj 1 (I cj 4 -1) exp k cj 2 (I cj 4 -1) 2 (Fa cj 0 ) ⊗ (Fa cj 0 ) F cj gr a cj 0 2 ; ( 15b
)
ϕ m σ m = 2̺ m (t) k m 1 (I m 4 -1) exp k m 2 (I m 4 -1) 2 (Fa m 0 ) ⊗ (Fa m 0 ) F m gr a m 0 2 (15c) + ̺ m (t) ̺ 0 (0) σ actmax 1 - λ max -λ m λ m 0 2 (λ max -λ 0 ) 2 (Fa m 0 ) ⊗ (Fa m 0 ) Fa m 0 2
.

where B e el = F e el F e el T is the elastic part of the left Cauchy-Green stretch tensor for elastin.

In Equation (15c), λ act was replaced by λ m λ m 0 , which means that we assume there is a continuous increase of the active stretch of SMCs during ATAA progression ( λm > 0). This is a different assumption of the one made by [START_REF] Braeu | Homogenized constrained mixture models for anisotropic volumetric growth and remodeling[END_REF] who assumed fast remodeling of SMCs, which implied that SMCs were able to keep the same active stretch throughout the G&R process (λ act = 1 meaning that λact = λ m ). In the current work, we want to investigate computationally the opposite situation, where SMC would proliferate at an extremely low rate [START_REF] Owens | Molecular regulation of vascular smooth muscle cell differentiation in development and disease[END_REF], which implies that the active stretch has to accommodate ATAA expansion.

Growth and Remodeling

In this work, a two-layer arterial model was considered where the rate of mass degradation or deposition of collagen in both layers was computed, such as

̺c j 0 (t) = ̺ c j 0 (t)k c j σ σ m -σ m h σ m h + Ḋc j g in the media, ( 16a
)
̺c j 0 (t) = ̺ c j 0 (t)k c j σ σ c j -σ c j h σ c j h + Ḋc j g in the adventitia, ( 16b
)
where k c j σ stands for (typically constant) collagen growth parameter, σ m h and σ c j h denote the average SMCs and collagen fiber stresses at homeostasis, σ m and σ c j are the current stress of the extant SMCs and collagen fibers, respectively. Moreover, it is assumed that elastin can be only subjected to degradation if its mass loss cannot be compensated by new elastin deposition ( ̺e 0 (t) = Ḋe g ). In this case, Ḋe g is the so-called generic function of the local elastin degradation rate, defined as

Ḋe g (X, t) = - ̺ e 0 (X, t) T e - D max t dam ̺ e 0 (X, 0)e -t t dam , ( 17
)
whose objective is to describe additional mass deposition or mass degradation (e.g., as a result of any damage) in elastin. D max is the maximum damage length, X is the material position, L dam and t dam are the spatial and the temporal damage spread parameters, respectively, and T e is the average turnover time for elastin constituent.

Having mass turnover relationships, it is now possible to capture the inelastic deformation induced by G&R. The evolution of the inelastic remodeling deformation gradient of constituent i at time t can be evaluated by solving the following system of equations (Cyron et al., 2016;[START_REF] Braeu | Homogenized constrained mixture models for anisotropic volumetric growth and remodeling[END_REF] ̺i 0 (t) ̺ i 0 (t)

+ 1 T i S i -S i pre = 2 ∂S i ∂C i el : C i el L i r F mix ,F i g , ( 18
)
where S is the second Piola-Kirchhoff stress tensor,

L i r = Ḟi r F i r -1
is the remodeling velocity gradient, T i is the period within which a mass increment is degraded and replaced by a new mass increment, known as the average mass turnover time. It is worthwhile noting that S i pre denotes the deposition pre-stress, the net mass production rate is already defined in Equation ( 16). In the case of elastin, it is assumed that this constituent is not produced any longer during adulthood and it undergoes a slow degradation with a half-life time of several decades (Cyron and Humphrey, 2016;[START_REF] Braeu | Homogenized constrained mixture models for anisotropic volumetric growth and remodeling[END_REF]. Consequently, the remodeling velocity gradient and, subsequently, the remodeling deformation gradient are zero (L e r = F e r = 0). Hence, elastin growth can essentially be computed on the basis of its degradation rate Ḋe g , defined at Equation ( 17). The inelastic growth deformation gradient is derived by summing the growth-related deformation gradient rates of each individual constituent, such as

Ḟi g = Ḟg = n i=1 ̺i 0 (t) ̺ 0 (t) F i g -T : a i g ⊗ a i g a i g ⊗ a i g , (19)
where a i g is a unit vector along the growth direction per individual constituents which, for instance, can represent an anisotropic growth across the direction of the arterial wall thickness and ̺ 0 (t) = i ̺ i 0 (t).

Numerical Implementation

A three-dimensional (3D) structural mesh, made of hexahedral elements, is reconstructed across the arterial wall using the Abaqus finite-element package [START_REF] Hibbit | Abaqus-Theory Manual[END_REF]. Each element has a regionally different mass density and contains a mixture of elastin, collagen fiber families and SMCs. The structural mesh implies the edge of each element is locally aligned with the material (i.e., radial, circumferential, and axial) directions of the arterial wall. In the case of non-perfectly cylindrical geometries, the radial direction is defined as the outward normal direction to the luminal surface and the axial direction is considered parallel to the luminal centerline along the blood flow stream. Finally, the circumferential direction is defined perpendicularly to the other two introduced directions.

Parameters

Values Units α c1 , α c2 , α c3 , and α c4 are axial, circumferential, and two diagonal directions of collagen fiber families, respectively.

α cj , j = 1, 2, 3, 4 {0, π 2 , ± π 4 } [Rad] µ e 82 [J/kg] κ 100µ e [J/kg] k cj ,c 1 = k m,c 1 15 [J/kg] k cj ,c 2 = k m,c 2 1.0 [-] k cj ,t 1 105 [J/kg] k cj ,t 2 0.13 [-] k m,t 1 10 [J/kg] k m,t 2 0.1 [-] ̺ e 0 250 [kg/m 3 ] ̺ cj 0 460 [kg/m 3 ] ̺ m 0 280 [kg/m 3 ] λ e 0z 1.3 [-] λ c 0j 1.1 [-] λ m 0 1.1 [-] λ 0 0.8 [-] λ max {1, 1.1, 1.4} [-]
The proposed model is implemented in Abaqus, through a coupled User-defined MATerial subroutine (UMAT) [START_REF] Hibbit | Abaqus-Theory Manual[END_REF]. The evolution of the arterial wall deformation is obtained by approximating the solution of the non-linear system of equations using the FEM approach, incremented through the Newton-Raphson method. Each time step of the simulation represents 1 month in the real (physical) timescale. The G&R deformation gradients are obtained at each time step on the basis of the stresses assessed at the previous step. It is worthwhile noting that the initial time step is only assumed to satisfy homeostatic conditions and the next time step is triggered by arterial dilatation through the G&R progression.

Description of the Case Study

Here, the model of Mousavi et al. ( 2019) is applied on a real human ATAA geometry in order to predict G&R of a patientspecific arterial wall. To do so, the preoperative Computerized Tomography (CT) scan of the patient as obtained after informed consent from a donor undergoing elective surgery for ATAA repair at CHU-SE (Saint-Etienne, France). Although the lumen of the aneurysm was clearly visible in the Digital Imaging and COMmunications (DICOM) file, it was not trivial to detect the aneurysm surface. A non-automatic segmentation of the CT image slices was performed using MIMICS (v. 10.01, Materialize NV) to reconstruct the ATAA geometry. The obtained geometry was assigned as the reference configuration and was structurally meshed with 7,700 hexahedral elements. A wall thickness of 2.38 mm was defined evenly in the reference configuration, yielding an average thickness of 2.67 mm at zero pressure state as the measured thickness in the supplied sample [START_REF] Farzaneh | Identifying local arterial stiffness to assess the risk of rupture of ascending thoracic aortic aneurysms[END_REF]. Figure 1 shows the computational domain in which the localized elastin degradation region is illustrated along with the media and adventitia layers.

The material parameters (Table 1), such as deposition stretches of collagen and other material parameters are the ones reported in references [START_REF] Bellini | A microstructurally motivated model of arterial wall mechanics with mechanobiological implications[END_REF][START_REF] Braeu | Homogenized constrained mixture models for anisotropic volumetric growth and remodeling[END_REF] or result from previous calibration against experimental data [START_REF] Davis | Local mechanical properties of human ascending thoracic aneurysms[END_REF]. Note that 97% of total elastin, 100% of total SMC, and 15% of total axial and diagonal collagen fibers are assigned to the media layer. Conversely, 3% of total elastin, 85% of total axial and diagonal collagen, and 100% of total circumferential collagen are assigned to the adventitia layer [START_REF] Bellini | A microstructurally motivated model of arterial wall mechanics with mechanobiological implications[END_REF][START_REF] Mousavi | Patient-specific stress analyses in the ascending thoracic aorta using a finite-element implementation of the constrained mixture theory[END_REF][START_REF] Mousavi | Computational predictions of damage propagation preceding dissection of ascending thoracic aortic aneurysms[END_REF]. The geometry is subjected to a luminal pressure of 80 mmHg (diastole), an axial deposition stretch of λ e 0z = 1.3 defined for elastin while the deposition stretches of collagen families and SMCs are set to λ

c j 0 = λ m 0 = 1.1.
The spatially varying circumferential deposition stretch of elastin is determined to ensure equilibrium with the luminal pressure using the iterative algorithm presented in [START_REF] Mousavi | Patient-specific stress analyses in the ascending thoracic aorta using a finite-element implementation of the constrained mixture theory[END_REF]. The applied boundary conditions to the geometry are defined as fixed at ends of the ATAA model in axial and circumferential directions (z and θ ) while free (no-traction) boundary condition are assigned in the radial direction (r). [START_REF] Guzzardi | Valve-related hemodynamics mediate human bicuspid aortopathy: insights from wall shear stress mapping[END_REF] found that regions of the aortic wall subjected to the localized effects of the jet flow expelled from the heart underwent greater elastin degradation associated with localized vessel wall remodeling. Therefore, and on the basis of these findings, a localized elastin degradation is considered here and its G&R effects are numerically simulated on the patient-specific ATAA geometry. The localized region of elastin degradation is shown in Figure 1.

Sensitivity Analysis

A parametric study was carried out about the effects of three parameters on G&R progression:

• the time at which the artery is maximally damaged (t dam ),

• the rate of collagen deposition (k c j σ /T c j ),

• the maximum contractility of SMCs (λ max ).

The different values for t dam follow the assumption that patientspecific G&R strongly relates to the temporal and spatial distributions of elastin degradation [START_REF] Watton | A mathematical model for the growth of the abdominal aortic aneurysm[END_REF][START_REF] Bellini | A microstructurally motivated model of arterial wall mechanics with mechanobiological implications[END_REF][START_REF] Braeu | Homogenized constrained mixture models for anisotropic volumetric growth and remodeling[END_REF]. Different values of k c j σ /T c j were considered since the rate of collagen deposition determines how fast the tissue adapts to elastin loss [START_REF] Mousavi | Patient-specific predictions of aneurysm growth and remodeling in the ascending thoracic aorta using the homogenized constrained mixture model[END_REF]. As the focus of this work is on the active contribution of SMCs on G&R, different values were also considered for λ max , possibly corresponding to possibly pathologic variations of intracellular filament structure in SMCs [START_REF] Liu | Role of cyclic strain frequency in regulating the alignment of vascular smooth muscle cells in vitro[END_REF][START_REF] Mantella | Variability in vascular smooth muscle cell stretch-induced responses in 2D culture[END_REF].

RESULTS

SMCs Active Stress Contribution

In order to quantify the contribution of SMCs active tone in the stress distribution field, two ratios were considered, namely (1) the ratio of the SMCs active stress with respect to the SMCs (active and passive) stresses (R m ); and (2) the ratio between the SMCs active stress and the total stress field (R t ).

R m = σ m act σ m ; R t = σ m act σ t . ( 20
)
These ratios are shown in percentage depending on the maximum SMCs active stretches for different rates of collagen deposition (Figure 2) and damage times (Figure 3) due to elastin degradation.

Figure 2 shows that these ratios are different depending on the amount of collagen deposition. The active stress includes an average value (dotted-lines in Figure 2) between 25 and 85% of the total SMCs stresses. However, this average value has reduced to an approximate value of R m = 20% in the case of λ max = 1.0 when the SMCs lose their functionality. This contribution of active stress with respect to the total stress field of the mixture is the average value between 14 and 28%. Indeed, the SMCs contribution to the total stress field has reduced down to R t = 4% in the critical condition of λ max = 1.0. Notice here that, in general, the solution has computationally converged faster by increasing the rate of collagen deposition as the artery becomes thicker and less displacement is expected. The convergence criteria state that the average value of the artery dilation should not exceed 10 -5 m with respect to the previous time step and reads

(u Ave ) n -(u Ave ) n-1 < 0.01; u Ave = 100 × N e i=1 ||u i || N e × T , (21
) where u Ave is the average displacement of all mesh elements of the artery, ||u i || = u 2 1 + u 2 2 + u 2 3 is the norm of elemental displacement, n stands for the time step, T = 2.5 × 10 -3 denotes the thickness of the artery and N e is the number of Hexahedral elements.

A similar behavior of SMCs active contribution can be observed in Figure 3 for different values of damage time t dam . In this case, the active stress has averagely represented between 25 and 45% of SMCs stresses and between 13 and 17% of the total stress. This has reduced to <10% in the critical condition of λ max = 1.0.

SMCs Stress-Stretch Curve

Figure 4 shows SMCs stress-stretch curves. In fact, the SMCs active stress in the circumferential direction is plotted vs. the SMCs stretches for different values of maximum SMCs contractility (λ max ). The curves are comparatively depicted in Figure 4A for various speeds of collagen deposition (k c j σ /T c j ) when the maximum active stretch is larger than the initial deposition stretch (λ max > λ m 0 ). It can be observed that the active stress is initially increased while the artery is stretched due to G&R. This monotonically increasing trend remains up to the point in which SMCs reach their maximum contractility and, then, as expected from their Hill functional behavior, it begins to decrease depending on the amount of collagen deposition.

Figure 4B illustrates the stress-stretch curve for various rates of collagen deposition (k c j σ /T c j ) when the maximum active stretch is larger than the initial deposition stretch (λ max ≃ λ m 0 ). The above-mentioned explanations regarding SMCs behavior still hold by noting that the stress has almost reached its maximum value for all values of k c j σ . Moreover, the decrease in maximum contractility (λ max ) has led to a larger decrease in active stress in comparison to the previous case, once reaching the maximum active stress value.

Interestingly, the active stress even drops and reaches zero (see Figure 4III). In this case, SMCs may fail to adhere onto and pull the extracellular matrix and eventually go into apoptosis. This critical situation happens no matter what rate of collagen deposition is applied.

Figure 5 depicts the SMCs stress-stretch curves compared for different values of maximum SMCs contractility (λ max ) and plotted for various damage times t dam . As clearly shown in Figure 5I, the larger damage time has led to more active stress since the artery is allowed to be more dilated, and then develop more stress, before it is damaged due to elastin degradation. In this case in which the maximum active stretch is larger than the initial deposition stretch, the active stress has an increasing trend from the beginning and all curves reach their corresponding maximum stresses and retain it. This behavior can still be related to the inverted parabola shape in Figure 5II when the maximum SMCs contractility is equivalent to the initial deposition stretch. This is evident as the first step (homeostasis) shows no change in developed stress (horizontal line). As shown in Figure 5III, the critical condition of a zero active stress is reached again after some time of G&R.

Figure 6 shows the distribution of the first principal stress distribution in the arterial wall for different values of SMC maximum contractility and damage time. The enlargement of the stress range in the artery is evident as the maximum SMCs stretches are reduced. Increasing the damage time results in increasing the stress values as well, with maximum values reached in the region of local elastin degradation. The distribution of the normalized collagen density (with respect to the arterial wall thickness) is depicted in Figure 7 for different values of λ max and k c j σ .

Stress Evolution

Figure 8 shows the evolution of stretches, active stresses and total stresses of SMCs along with the total stress field of the mixture. The solutions are compared for different values of maximum SMCs active stretches and for different rates of collagen deposition, all for a fixed value of t dam = 40. The stretch evolution shows that the ATAA dilation evolves in the opposite direction of the rate of collagen deposition. This occurs as more collagen is deposited for compensating elastin degradation: the thicker the artery and the smaller the displacement. This can be related to the evolution of the active stresses in which the inverted parabola behavior is more evident for the bigger value of k c j σ . In fact, when deposited collagen has a larger contribution to the stress field, SMCs stresses reach a maximum and, then, decrease afterwards. In fact, the longer G&R occurs, the more stress is transferred from the media to the adventitia layer.

The evolutions shown in

Aneurysm Growth

Figure 11 depicts the evolution of the aortic diameter and thickness during G&R. The results are compared for different values of λ max where various rates of collagen deposition are investigated. As expected, faster collagen deposition results in less dilation of the aortic wall and, therefore, faster convergence of the numerical solution. This trend can be observed for both diameter and thickness whose opposite trend is inevitable. In addition, the effects of SMCs maximum contractility, especially in the case of diameter, is evident as a smaller λ max has led to an increase in the diameter. However, this increasing trend of diameter does not completely agree with the thickness of the artery. It seems that thickness decreased in the critical condition (λ max = 1.0).

The same behavior can be observed in Figure 12 for different values of damage time. However, it can be generally concluded that whatever the conditions of SMCs active stress, the change in the rate of collagen deposition or damage time does not have a major influence on changing the diameter/thickness of the arterial wall.

The contour plot of the arterial wall thickness and displacement distributions are shown in Figure 13 for the representative case of k c j σ = 0.05 and t dam = 40 in the media layer for different values of maximum SMCs contractility.

Figure 14 also illustrates the arterial wall dilation compared to its initial state through thickness and displacement distributions after 15 years of G&R progression. The maximum displacement and minimum thickness occurred at the location of the localized elastin loss.

DISCUSSION

In this paper, growth and remodeling of a patient-specific aortic geometry has been computationally modeled using the homogenized CMT in order to predict ATAA evolution while the artery has undergone proteolysis of elastin localized in regions of deranged hemodynamics. To the best of the authors' knowledge, it is the first time that the effects of possible variations in SMC contractility on ATAA progression are simulated in a multilayer patient-specific model. The results reveal a possible critical condition in which the active stress of SMCs can reduce to zero, which can be interpreted as if SMCs do not pull onto the extracellular matrix anymore, potentially leading to apoptosis.

The paramount importance of the SMCs role in short-and long-term changes of arterial responses is reported in [START_REF] Cox | Arterial wall mechanics and composition and the effects of smooth muscle activation[END_REF][START_REF] Cox | Regional variation of series elasticity in canine arterial smooth muscles[END_REF], [START_REF] Humphrey | A constrained mixture model for growth and remodeling of soft tissues[END_REF], [START_REF] Liu | Role of cyclic strain frequency in regulating the alignment of vascular smooth muscle cells in vitro[END_REF], and [START_REF] Murtada | Experiments and mechanochemical modeling of smooth muscle contraction: significance of filament overlap[END_REF][START_REF] Murtada | Adaptation of active tone in the mouse descending thoracic aorta under acute changes in loading[END_REF]. The active tension of SMCs, produced over a large range of muscle lengths, can be described through a length-tension relationship (Dorbin, 1973;[START_REF] Rachev | Theoretical study of the effects of vascular smooth muscle contraction on strain and stress distributions in arteries[END_REF][START_REF] Gunst | Cytoskeletal remodeling of the airway smooth muscle cell: a mechanism for adaptation to mechanical forces in the lung[END_REF][START_REF] Zulliger | A constitutive formulation of arterial mechanics including vascular smooth muscle tone[END_REF][START_REF] Carlson | A theoretical model for the myogenic response based on the length-tension characteristics of vascular smooth muscle[END_REF][START_REF] Herrera | Sarcomeres of smooth muscle: functional characteristics and ultrastructural evidence[END_REF][START_REF] Syyong | Adaptive response of pulmonary arterial smooth muscle to length change[END_REF][START_REF] Yamin | Deciphering actin cytoskeletal function in the contractile vascular smooth muscle cell[END_REF] and Morgan, 2012). In this work, the classic inverted parabola shape is used for the length-tension relationship as it, specifically, describes the overlap between the actin and myosin filaments [START_REF] Murtada | Experiments and mechanochemical modeling of smooth muscle contraction: significance of filament overlap[END_REF].

Compared to the well-known model proposed by [START_REF] Zulliger | A constitutive formulation of arterial mechanics including vascular smooth muscle tone[END_REF], our model, which stems from the work of [START_REF] Murtada | Experiments and mechanochemical modeling of smooth muscle contraction: significance of filament overlap[END_REF], has a number of differences which are:

1. [START_REF] Zulliger | A constitutive formulation of arterial mechanics including vascular smooth muscle tone[END_REF] have a multiplicative split between the active and passive stress of SMCs whereas we use an additive split; 2. [START_REF] Zulliger | A constitutive formulation of arterial mechanics including vascular smooth muscle tone[END_REF] considered effects of the change of calcium concentration with the stretch, which we disregarded, as these effects mainly manifest for fast deformations, whereas we modeled deformations over large timescales; 3. for the length tension relationship, [START_REF] Zulliger | A constitutive formulation of arterial mechanics including vascular smooth muscle tone[END_REF] used a Heaviside function whereas we used an inverted parabola function.

It is worth noting that it is the first time that variations of the active stress are investigated in G&R progression through computational modeling. The situations that we have considered for these variations of the active stress remain partly fictitious because there are no in vivo measurements of the active stress in vivo.

In this work, a parametric study has been performed on the basis of the time at which the artery is maximally damaged (t dam ), the rate of collagen deposition (k c j σ /T c j ) and the maximum contractility of SMCs (λ max ) in order to study the effects of active SMCs on G&R of a patient-specific ATAA geometry. Results illustrate the importance of damage time on the expansion rate of the artery such that the larger value of t dam leads the higher rate of G&R. Note that elastin degradation during ATAA growth may be related to multiple biological and mechanical parameters including abnormal distribution of Wall Shear Stress (WSS) [START_REF] Guzzardi | Valve-related hemodynamics mediate human bicuspid aortopathy: insights from wall shear stress mapping[END_REF] and circumferential stress [START_REF] Humphrey | A constrained mixture model for growth and remodeling of soft tissues[END_REF]. Our simulations highlight some dependence between the active stress and damage time, especially when the maximum active stretch was larger than the deposition stretch. However, this dependency is dramatically reduced in the other two cases when the maximum SMCs contractility is equal or less than the deposition stretch. The evolution of the SMCs total (active and passive) stresses indicates that the passive contribution is increasing all the time. In fact, increasing both the damage time and maximum SMCs contractility have resulted in an increase in SMCs stresses and in total stresses of the mixture.

It appeared that active stress decreases more when the rate of collagen deposition increases. This makes sense as collagen can take over a part of the stresses borne by SMCs. For instance, in the case of k c j σ = 0.05, the active stress increases and remains at its maximum value since the amount of collagen deposited in the media layer is not sufficient to compensate the loss of active tone in SMCs. This is an opposite trend compared to the case k c j σ = 0.30 in which the fast deposition of collagen prevents the SMCs to reach their maximum contractility and a comparative fast decrease in active stress can be observed. It must be noted that SMC contractility does not dramatically affect the arterial dilation. Considering the evolution of the SMCs total (passive and active) stresses, the monotonically increasing trend of passive stresses can be acknowledged. In some cases (say k c j σ > 0.15), though, the active stresses dominate the behavior of SMCs total stresses. This can also be observed on the evolution of the total mixture stresses which begin to increase once the G&R starts. The decrease in total mixture stresses indicates that if the rate of collagen deposition reaches a certain threshold, the artery becomes thick enough to enable a dramatic reduction of the wall stress. This certainly affects the speed of convergence of the numerical solutions as the computational time was significantly lengthened for these situations.

It must be noted that there are still several limitations and technical challenges associated with the presented model:

• The presented model relies on a number of assumptions.

One of these assumptions, stating that initial reference is at homeostasis, is typical of G&R [START_REF] Kassab | Mechanical homeostasis of cardiovascular tissue[END_REF]. Although this assumption may be actually satisfied globally at the scale of an organ, we also assume that it is satisfied point-wisely in our G&R model of the arterial wall. • Considering a uniform thickness for the arterial wall is a limitation of the current work since the distribution of material and structural parameters (e.g., thickness, fiber orientations) of a patient-specific geometry should be consistent with in vivo data. • In this work, the constitutive parameters for the patientspecific model were estimated by curve fitting from the ex vivo bulge inflation data of an ATAA segment excised after the surgical intervention of the same patient. However, the in vivo material properties of ATAAs [START_REF] Farzaneh | Identifying local arterial stiffness to assess the risk of rupture of ascending thoracic aortic aneurysms[END_REF] should be identified non-invasively in clinical applications. • We have generated fictitious situations where the active stretch of SMCs evolves with the actual stretch of the tissue and have echoed this evolution onto the active stress through the inverted parabola function. Although these situations can be justified if SMC proliferate at an extremely low rate and do not adapt to aneurysm progression, they remain fictitious and there is a pressing need to characterize the evolution of SMCs length-tension relationships in aortic aneurysms in vivo in order to improve the predictions of G&R computational models.

• Several biological pathways may be activated during aneurysm development and induce adaptation mechanisms of SMCs that were not considered in this work. More specifically, the distribution of WSS is complex in ATAAs [START_REF] Condemi | Fluid-and biomechanical analysis of ascending thoracic aorta aneurysm with concomitant aortic insufficiency[END_REF] and this may have an impact on SMCs which should be included in future developments of the model.

CONCLUSION

In this work, the contribution of SMCs contractility during aneurysmal growth and remodeling has been computationally investigated with a multi-layer patient-specific model of the human aorta subjected to localized elastin degradation. The model relies on the homogenized constrained mixture theory.

The results show that the active contribution of SMCs is of paramount importance in analyzing the total stress field of the mixture. It, in fact, affects the rates of collagen deposition by which the elastin loss is gradually compensated. Interestingly, a critical condition has been observed, under which the active stress of SMCs could drop and possibly result in apoptotic conditions. Future work will carry on investigating the major role of SMCs contractility on G&R in aneurysm development, considering supplemental biological pathways in our computational model in order to account for actual adaptation mechanisms undergone by SMCs.

FIGURE 1 |

 1 FIGURE 1 | The computational domain of the ATAA geometry composed of the localized elastin degradation region, media layer and adventitia layer shown in yellow, gray, and green colors, respectively.

FIGURE 2 |

 2 FIGURE 2 | The evolution of the stress ratios {R m , R t } with respect to physical time t. The dotted lines denote the average value for each color line correspondingly. The results obtained for a two-layer thick-wall patient-specific human ATAA responding to a regional elastin degradation for 15 years. Different rates of collagen deposition k cj σ /T cj are considered depending on the (column-wise presentation of the) maximum SMCs contractility (I) λ max = 1.4; (II) λ max = 1.1; and (III) λ max = 1.0 for a fixed value of t dam = 40.

FIGURE 3 |

 3 FIGURE 3 | The evolution of the stress ratios {R m , R t } with respect to physical time t. The dotted lines denote the average value for each color line correspondingly. The results obtained for a two-layer thick-wall patient-specific human ATAA responding to a regional elastin degradation for 15 years. Different damage times t dam are considered depending on the (column-wise presentation of the) maximum SMCs contractility (I) λ max = 1.4; (II) λ max = 1.1; and (III) λ max = 1.0 for a fixed value of k cj σ = 0.05.

FIGURE 5 |

 5 FIGURE 5 | SMCs stress-stretch (σ m 22act -λ m ) curve: the results obtained for a two-layer thick-wall patient-specific human ATAA responding to a regional elastin degradation for 15 years. Different damage times t dam [days] are considered depending on the (column-wise presentation of the) maximum SMCs contractility (I) λ max = 1.4; (II) λ max = 1.1; and (III) λ max = 1.0 for a fixed value of k cj σ = 0.05.

FIGURE 6 |

 6 FIGURE 6 | Principal stresses: the stress distributions are comparatively shown at the final physical time (180 months). The results obtained for a two-layer thick-wall patient-specific human ATAA responding to a regional elastin degradation for 15 years. Different damage times t dam [days] are considered depending on the (column-wise presentation of the) maximum SMCs contractility (I) λ max = 1.4; (II) λ max = 1.1; and (III) λ max = 1.0 for a fixed value of k cj σ = 0.05. Each column shows the stress distributions, from left to right, for the interior of the arterial wall (cut in half) and for the media layer, respectively. The range of principal stress distributions are (A) [6.4e4-3.3e5]; (B) [7.1e4-3.6e5]; (C) [7.2e4-6.3e5]; (D) [8.2e4-6.2e5]; (E) [7.3e4-7.4e5]; (F) [7.4e4-8.8e5]; (G) [1.1e5-8.1e5]; (H) [7.3e4-8.9e5]; (I) [8.0e4-1.0e6] in [Pa].

FIGURE 7 |

 7 FIGURE 7 | Normalized collagen density: the mass distributions are comparatively shown at the final physical time (180 months). The results obtained for a two-layer thick-wall patient-specific human ATAA responding to a regional elastin degradation for 15 years. Different rates of collagen deposition k cj σ /T cj are considered depending on the (column-wise presentation of the) maximum SMCs contractility (I) λ max = 1.4; (II) λ max = 1.1; and (III) λ max = 1.0 for a fixed value of t dam = 40 [days]. Each column shows the stress distributions, from left to right, for the interior of the arterial wall (cut in half) and for the adventitia layer, respectively. The range of normalized collagen density distributions are (A) [0.95-1.6]; (B) [0.99-1.8]; (C) [1.0-1.8]; (D) [0.99-1.0]; (E) [1.0-1.9]; (F) [1.1-4.5]; (G) [1.1-4.2]; (H) [1.5-4.3]; (I) [5.1-3.2].

Figure 9

 9 are similar to the ones shown in Figure 8, but this time for different values of damage time t dam . In this case, increasing the damage time increased the developed stretch due to G&R in the artery. The absolute value of the principal stress distributions are shown in Figure 10 for the representative case of k c j σ = 0.05 and t dam = 40. The results are compared for different values of the maximum SMCs stretches shown in media and adventitia layers.

FIGURE 11 |

 11 FIGURE 11 | Evolution of arterial diameter with respect to the physical time. The results obtained for a two-layer thick-wall patient-specific human ATAA responding to a regional elastin degradation for 15 years. Different rates of collagen deposition k cj σ /T cj are considered depending on the (column-wise presentation of the) maximum SMCs contractility (I) λ max = 1.4; (II) λ max = 1.1; and (III) λ max = 1.0 for a fixed value of t dam = 40 [days].

FIGURE 12 |

 12 FIGURE 12 | Evolution of artery thickness with respect to physical time t. The results obtained for a two-layer thick-wall patient-specific human ATAA responding to a regional elastin degradation for 15 years. Different damage time t dam [days] are considered depending on the (column-wise presentation of the) maximum SMCs contractility (I) λ max = 1.4; (II) λ max = 1.1; and (III) λ max = 1.0 for a fixed value of k cj σ = 0.05.

FIGURE 13 |

 13 FIGURE 13 | Arterial wall thickness and displacement evolutions: the results obtained for a two-layer thick-wall patient-specific human ATAA responding to a regional elastin degradation during 15 years. The wall thickness distributions are shown at t = {20, 80, 120} [Months] depending on the (column-wise presentation of the) maximum SMCs contractility (I) λ max = 1.4; (II) λ max = 1.1; and (III) λ max = 1.0 for a fixed value of k cj σ = 0.05 and t dam = 40 [days]. Each column shows, from left to right, the thickness [0.78-1.3] and displacement [5.5e-4-3.5e-3] in [mm], of the arterial wall in the media layers, respectively.

FIGURE 14 |

 14 FIGURE 14 | The thickness and displacement distributions of the arterial wall: the results obtained for a two-layer thick-wall patient-specific human ATAA responding to a regional elastin degradation for 15 years. In this representative case, the collagen deposition gain is k cj σ = 0.05, the maximum SMCs contractility is λ max = 1.4 and the damage time is t dam = 40 [days]. From left to right, the thickness [0.78-1.3] and displacement [5.9e-4-3.5e-3] in [mm], of the arterial wall are depicted.

TABLE 1 |

 1 Material

parameters employed for two-layer patient-specific human ATAA model adapted from

[START_REF] Mousavi | Patient-specific predictions of aneurysm growth and remodeling in the ascending thoracic aorta using the homogenized constrained mixture model[END_REF]

.

Frontiers in Bioengineering and Biotechnology | www.frontiersin.org

November 2020 | Volume 8 | Article 587376

FUNDING

This work was supported financially by the European Research Council ERC through grant ERC-2014-CoG BIOLOCHANICS.

DATA AVAILABILITY STATEMENT

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT

The studies involving human participants were reviewed and approved by CHU Saint-Etienne. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS

AG: software implementation, simulations, validation, analysis of results, and writing-original draft presentation. All authors supervision of the research, methodology, funding, analysis of results, reviewing, and editing of manuscript.

Conflict of Interest:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Ghavamian, Mousavi and Avril. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.


Fabrice Barbe
email: fabrice.barbe@insa-rouen.frquey@emse.fr

Romain Quey

Andrei Musienko
email: andrei.musienko@mat.ensmp.fr

Georges Cailletaud
email: georges.cailletaud@ensmp.fr

Three-dimensional characterization of strain localization bands in high resolution elastoplastic polycrystals

Keywords: Strain localization, Crystal plasticity, 3D polycrystal, Finite elements

In crystalline materials, the experimental observation of the localization of plastic strains in particular directions is generally restricted to the surface of a sample containing some hundreds of grains, because of the difficulties underlying microstructural analysis. In these conditions, the determination of the morphological characteristics of localization can be limited by the poor statistical representativity of the domain of observation. The purpose of this work is to extend the analysis of localization -localization bands or elseto the three-dimensional elastoplastic strain fields of a high resolution representative volume element of a polycrystal.

Introduction

Elastoplastic strains in polycrystals are known to localize in the form of bands which are oriented in the direction of maximum shear stress, i.e. at 45 • with respect to the uniaxial tension direction. Such bands can be observed at the surface of plastically deformed polycrystals, either experimentally [START_REF] Doumalin | Computational and experimental investigations of the local strain field in elastoplastic two-phase materials[END_REF][START_REF] Doumalin | Characterisation of the strain distribution in heterogeneous materials[END_REF][START_REF] Gélébart | Identification of Crystalline Behavior on Macroscopic Response and Local Strain Field Analysis: Application to Alpha Zirconium Alloys[END_REF][START_REF] Héripré | Coupling between experimental measurements and polycrystal finite element calculations for micromechanical study of metallic materials[END_REF][START_REF] St-Pierre | 3D simulations of microstructure and comparison with experimental microstructure coming from O.I.M. analysis[END_REF] or by numerical simulations [START_REF] Cordier | Plastic Deformation of mantle minerals: Multiscale numerical modelling[END_REF][START_REF] Zeghadi | Ensemble averaging stressstrain fields in polycrystalline aggregates with a constrained surface microstructure -Part 2: crystal plasticity[END_REF][START_REF] Osipov | FE modelling of bainitic steels using crystal plasticity[END_REF][START_REF] Lebensohn | Modeling viscoplastic behavior and heterogeneous intracrystalline deformation of columnar ice polycrystals[END_REF]. They appear at the early stages of plasticity and can expand over few grains with a width of the order of the grain size or less. Some features of their morphology (orientation and length) have been determined in 2D by means of point-to-point covariance analysis of strain fields measured by SEM microextensometry [START_REF] Doumalin | Microextensométrie locale par corrélation d'images numériques[END_REF][START_REF] Doumalin | Computational and experimental investigations of the local strain field in elastoplastic two-phase materials[END_REF][START_REF] Doumalin | Characterisation of the strain distribution in heterogeneous materials[END_REF]: once segmented into two classes corresponding to lower and higher strains, the strain field is considered as a 1 two-phase medium whose morphology is analyzed in terms of covariance. The point-topoint covariance provides the probability to find a given phase at two points separated by a given vector; it thus enables to detect any periodic repetition of patterns in the spatial distribution of field variables and more generally to determine the range of the interactions between grains in a polycrystal.

Such information can be a prevailing factor for the relevance of modellings based on microstructural morphology: the size of the representative volume element, the depth over which boundary conditions affect the local behavior, the distance over which fields become uncorrelated, . . . and the degree of accuracy of the morphological features which should be taken into account might depend on it. As current experimental techniques cannot, to this day, give access to this kind of information concerning a large threedimensional polycrystal, it is proposed to extract this information from FE simulations of high resolution 3D polycrystals, that provide full-field solutions for the elastoplastic intracrystalline stresses and strains in a large number of grains. This numerical approach has been proved to be successfull in reproducing the intracrystalline strain localization experimentally measured on the surface of a multicrystal, provided that the mesh is fine enough [START_REF] Musienko | Threedimensional finite element simulation of a polycrystalline copper specimen[END_REF][START_REF] Héripré | Coupling between experimental measurements and polycrystal finite element calculations for micromechanical study of metallic materials[END_REF][START_REF] Zhao | Investigation of three-dimensional aspects of grain-scale plastic surface deformation of an aluminium oligocrystal[END_REF]. For the present study of the localization in a 3D polycrystalline medium mimicking a real microstructure with its random morphological nature, we resort to Voronoi tessellations of space and a fine discretization of each Voronoi polyhedron into several tens of elements, so that local fields inside each grain be accessible (Barbe et al., 2001a;[START_REF] Barbe | Intergranular and intragranular behavior of polycrystalline aggregates. Part 2: results[END_REF].

Section 2 recalls the principles of the numerical modelling and presents contour plots of axial stress and strain. The process of computing the point-to-point covariance of a segmented variable is described briefly in section 3. Then 2D and 3D covariance fields are analyzed in section 4.

Numerical simulation

Data for the analysis of stress and strain intracrystalline localization in elastoplastic polycrystals are obtained from FE simulations. The behavior of each integration point of the FE mesh is described by a constitutive model which was first introduced in [START_REF] Méric | Single Crystal Modeling for Structural Calculations. Part 1: Model Presentation[END_REF]. It has then been extensively used in several works which have laid an important groundwork (Barbe et al., 2001a;[START_REF] Diard | Evaluation of finite element based analysis of 3D multicrystalline aggregates plasticity. Application to crystal plasticity model identification and the study of stress and strain fields near grain boundaries[END_REF][START_REF] Sai | Micro-mechanical modeling of the inelastic behavior of directionally solidified materials[END_REF][START_REF] Osipov | FE modelling of bainitic steels using crystal plasticity[END_REF] and more particularly for the direct comparison of experimental and simulated intracrystalline strain fields in [START_REF] Musienko | Threedimensional finite element simulation of a polycrystalline copper specimen[END_REF].

The constitutive model is based on the Schmid law: the plastic strain ε ∼ g results from the slip on the crystallographic slip systems γ s (eq. 1); the onset of slip and the slip rate are determined in eq. 2 from the resolved shear stress τ s and the hardening variables r s and x s (isotropic and kinematic hardenings, respectively). Whereas τ s directly results from the state of stress σ ∼ and the Schmid tensor m ∼ s (eq. 3, where n - s is the unit vector normal to the slip plane and l s is the slip direction), the hardening variables follow non-linear evolutions defined by classical phenomenological formulations of hardening (eq. 4,5). The material parameters are those of a Ni-based alloy which has been under consideration in (Barbe et al., 2001a); it has an FCC structure (slip on the {111} 110 octahedral slip systems). γs = vs sign (τ s -x s ) ; ε

∼ g = s m ∼ s γs (1) vs = | τ s -x s | -r s K n with x = M ax(x, 0) and v s (t = t 0 ) = 0 (2) τ s = σ ∼ : m ∼ s = 1 2 σ ∼ : n - s ⊗ l - s + l - s ⊗ n - s (3) x s = c α s ; r s = R 0 + Q r h rs 1 -e -b v r (4) αs = γs -d α s vs with α s (t = t 0 ) = 0 (5)
At the current state of the study, the polycrystalline domains of simulation, with cubic shape, contain 13 × 13 × 13 = 2197 grains. According to previous analyses in ( [START_REF] Barbe | Intergranular and intragranular behavior of polycrystalline aggregates. Part 2: results[END_REF], any integration point at a distance lower than 2 Mean Grain Size (MGS) from the polycrystal boundaries must be excluded from the analysis to avoid any boundary effect. The remaining region at the centre is large enough (9 × 9 × 9 = 729 grains) to provide representative properties of the polycrystalline material at the mean polycrystal scale, as checked in the previously cited studies based on the same modelling. As a major improvement regarding these previous studies, the polycrystals microstructural morphologies are obtained from Voronoi tessellations which are slightly modified prior to free meshing. This procedure, which enables a rapid and robust meshing without local overrefinements, is implemented in the NePeR software [START_REF] Neper | NePeR: a software to generate random polycrystals for the finite element method (version 1.6[END_REF][START_REF] Quey | Free meshing of microstructures based on modified Voronoi tessellations[END_REF]. The resulting meshes conform to polyhedra boundaries, as illustrated on fig. 1a. The polycrystals in this study, made of 2197 grains, are discretized into 180000 tetrahedral quadratic elements, which leads to an average of 82 elements per grain. This mesh density respects the conditions under which local fields are unaffected by the type and number of elements (Barbe et al., 2001c;[START_REF] Osipov | FE modelling of bainitic steels using crystal plasticity[END_REF]. FE simulations are performed with the software Zset (Mines Paristech, ONERA, nwnumerics.com, [START_REF] Besson | Object-Oriented Programming Applied to the Finite Element Method. Part II. Application to Material Behaviors[END_REF]) on a Linux cluster, using an implicit method for the integration of the constitutive laws and a Newton type method for the resolution of the non-linear problem.

The polycrystals are subjected to uniaxial tension up to a total imposed axial strain of the order 1 %. Their lateral boundaries are left free of mean stress; they are also imposed to remain parallel to their initial positions, such as to induce a boundary effect similar, in average, to the one which would be exerted by surrounding grains of the medium from which the considered polycrystal is extracted. The obtained axial stress and axial strain contour plots are presented on fig. 1b,c. The main following classical features can be noticed: both stress and strain are heterogeneous at the scale of grains as well as inside grains; strain localization regions do not systematically correspond to stress localization regions and their patterns do not seem to be the same; strain seems to localize in the form of bands oriented at 45 • with respect to the tension direction (vertical, on the figure).

Statistical treatments

Given a field variable, say axial strain ε(X -) where X -= (x 1 , x 2 , x 3 ) is the position vector in the polycrystal frame, and an associated threshold ε t , the strain field is first segmented into two sets of points: those where ε < ε t and the others, forming point clouds that will be denoted phase B = {X -/ ε(X -) ≥ ε t }. The 3D visualization of these points does not provide particular information on the localization of higher strains: they form more or less large connex regions in the polycrystal space which seem to be equiprobably present and to propagate in different directions. The axial stress shows the same features. So, the characterization of stress and strain localization based on direct visualization of segmented variables can only be qualitative.

We thus resort to a morphological analysis, the point-to-point covariance, which enables the characterization of two-phase materials or any binary medium. It is defined according to a vector h -and corresponds to the probability that a point and its translated according to h -belong to the same phase. Let f B (X -) be the characteristic function of the phase B.

f B (X -) = 1 if X -∈ B 0 else
Following the formulation of point-to-point covariance in a continuous space, in a discrete space as the one provided by the set of integration points of the polycrystal, the covariance S εt (h -) of B is written as a function of the translation vector h -:

S εt (h -) = 1 3 i=1 (L -h i ) 3 i=1 L-h i x i =0 f B (X -)f B (X -+ h -)
where h -= (h i ) i=1,3 , L is the length of the edge of the cubic aggregate. The sum

3 i=1 L-h i
x i =0 corresponds to the sum over all the points within the range [0, L -h 1 ] × [0, L-h 2 ]×[0, L-h 3 ] i.e. to the volume average over the intersection set {X -}∩{X -+h -}.

Analysis

The obtained covariance field exhibits a central symmetry over h -= 0 -since a translation according to h -or -h -leads to the same intersection set. Taking MGS as the length scale, the covariance is defined in the range [-13, 13] × [-13, 13] × [-13, 13] since the polycrystals under concern in this work contains 13 × 13 × 13 grains. In order to make a qualitative comparison with the experimental characterizations of [START_REF] Doumalin | Microextensométrie locale par corrélation d'images numériques[END_REF][START_REF] Doumalin | Characterisation of the strain distribution in heterogeneous materials[END_REF] on a two-phase material and on a zirconium alloy under tension, the covariance field of axial strain is first visualized in two perpendicular planes of the translation space (h 1 , h 2 , h 3 ) which include the tensile axis (x 3 ): (h 1 = 0) in fig. 2a and (h 2 = 0) in fig. 2b. For both plots, high values of covariance appear close to the centre of the translation space: at h -= 0 -, the probability of finding the same phase B is maximum and S εt (0 -) takes the value of the volume fraction of the phase B (here, 0.5). For small distances h, S εt remains high with a spherical symmetry; this short distance isotropic localization is the consequence of local interactions independent of the microstructure, as in any FE simulation with fine enough mesh. With h from 1 up to 3, the highest values of S εt are distributed according to a cross with branches oriented at 45 • with respect to the tension direction (vertical) and the smallest values appear at 0 • and 90 • angles to the tension direction. Similar characterizations can be made concerning equivalent strain or the accumulated plastic strain. This indicates that localization of highest strains will more probably happen in the form of bands, 2 or 3 MGS long and oriented at 45 • to the tension axis, which is consistent with experimental observations of Doumalin et al on a two-phase material and a zirconium polycrystal. Furthermore, axial strain is unlikely to localize in bands parallel or perpendicular to the tensile direction. For h > 3, the distributions of S εt differ from one plane to the other, which means that the isotropy of the problem, though it is accurately respected at the mean polycrystal scale, reaches its limits at the local scale, as soon as the analysis of the interactions between grains is performed above 4 MGS. This limit can evidently be extended with larger polycrystals, but it is not in the scope of this paper. As h is increased, the size of the intersection set of points decreases, thus directly affecting the statistical representativity for the estimation of the covariance. For this reason and as however local fields are influenced by the boundary conditions on a depth of the order of 2 or 3 MGS, any analysis for h ≥ 10 should be avoided with this size of polycrystal.

In fig 3a, a 3D visualization of the covariance field is performed by means of isointensity surfaces corresponding to values of covariance above the one that distinctly delimits the cross in the 2D visualizations, taken at 0.24. Decreasing the value of S εt from its maximum 0.5, the shapes of the isointensity surfaces vary from a sphere to kind of a hourglass and then takes the shape of a three-dimensional cross whose eight branches are the diagonals of a parallelepiped with a length of the order of 3 MGS, thus no longer respecting an axial symmetry. As nothing at this state enables to ensure the statistical representativity of the strain patterns for h ≥ 3, the interpretation of the geometry of this eight-branches cross would require great care. The hourglass-like shape of isointensity surfaces is however confirmed in fig 3b, where the covariance field is obtained by ensemble averaging over three covariance fields obtained from three different 2197-grains polycrystals. This shape is characteristic of an equiprobable localization of the highest strains in all directions at a 45 • angle from the tension direction.

Conclusion

The characterization of the three-dimensional distribution of highest strains in an elastoplastic polycrystal has been performed by means of point-to-point covariance of segmented strain field which was computed in a polycrystal with high intracrystalline resolution. It has shown the propensity for strains to localize equiprobably along bands 2 or 3 MGS long and oriented at around 45 • with respect to the tensile direction. As to this day experimental investigations can hardly provide measures which are statistically representative of three-dimensional localization, this work constitutes the first step ever towards any such detailed and systematic 3D characterization. Future work will extend the analysis to different loadings and materials computed on larger polycrystals and, besides, should refine the characterization as in [START_REF] Doumalin | Characterisation of the strain distribution in heterogeneous materials[END_REF] by resorting to a correlation function based on real values of strains instead of segmented values.

Figure 1 :

 1 Figure 1: Strain and stress fields in a high resolution polycrystal submitted to 0.5 % uniaxial tension. (Left: a) Polycrystal morphology and mesh (2197 grains freely meshed into 180000 tetrahedral elements). (Middle: b) axial strain; (Right: c) axial stress. Colors on the online version.

Figure 2 :

 2 Figure 2: Covariance of segmented axial strain for translations in all the directions of two planes containing the tensile axis (x 3 ): (h 1 = 0) (Left: a) and (h 2 = 0) (Right: b). The length scale of the axis are given in terms of mean grain size. Colors on the online version.

Figure 3 :

 3 Figure 3: Isointensity surfaces of covariance of segmented axial strain in the threedimensional translation space: from a single 2197-grains polycrystal (Left: a) and obtained by ensemble averaging over three 2197-grains polycrystals (Right: b). In both cases selected, covariances are superior to 0.24. Colors on the online version.

Acknowledgements

Thanks are due to Guido Lodato (CORIA, CNRS INSA Rouen) for his selfless introduction to the Paraview visualization application. CRIHAN (Centre de Ressources Informatiques de HAute Normandie) is acknowledged for providing computing ressources (Project N • 2008012).