Di Zuo

Stéphane Avril

Chunjiang Ran

Haitian Yang

S Jamaleddin Mousavi

Klaus Hackl

Yiqian He
email: heyiqian@dlut.edu.cn

Sensitivity analysis of non-local damage in soft biological tissues

Keywords: Nonlocal damage, Uncertainty analysis, Soft tissues, Surrogate model, Computational efficiency

Computational modeling can provide insight into understanding the damage mechanisms of soft biological tissues. Our gradient-enhanced damage model presented in a previous publication has shown advantages in considering the internal length scales and in satisfying mesh independence for simulating damage, growth and remodeling processes. Performing sensitivity analyses for this model is an essential step towards applications in which it is necessary to handle the uncertainty caused by patient-specific data. In this paper, a numerical analysis approach is developed by integrating two existing methods, i.e. the gradient-enhanced damage model and the surrogate model-based probability analysis, to address this need. To increase the computational efficiency of the Monte Carlo method in uncertainty propagation for the nonlinear hyperelastic damage analysis, the surrogate model based on Legendre polynomial series is employed to replace the direct FEM solutions, and the sparse grid collocation method (SGCM) is adopted for setting the collocation points to further reduce the computational cost in training the surrogate model. The effectiveness of the proposed approach is illustrated by two numerical examples, including an application to the clinical problem of dilatation of the artery.

Introduction

Understanding the damage mechanism of soft biological tissues is critical to the characterization of tissue injuries [START_REF] Mousavi | Computational predictions of damage propagation preceding dissection of ascending thoracic aortic aneurysms[END_REF][START_REF] Fleischmann | Vacuum sealing as treatment of soft tissue damage in open fractures[END_REF][START_REF] Kloner | Studies of experimental coronary artery reperfusion. Effects on infarct size, myocardial function, biochemistry, ultrastructure and microvascular damage[END_REF][START_REF] Holzapfel | Modeling of damage in soft biological tissues[END_REF]. For instance, in the treatment of cardiovascular diseases, the focus is often on the potential risk of tissue rupture, which is mainly caused by soft tissue damage, such as atherosclerotic plaques or aneurysms [START_REF] Trabelsi | Patient specific stress and rupture analysis of ascending thoracic aneurysms[END_REF]. The damage analyses of the tendons, which are important load-bearing structures and are frequently injured in both sports and work, are also important in understanding the underlying mechanisms and are helpful for treatment [START_REF] Svensson | Fracture Mechanics of Collagen Fibrils: Influence of Natural Cross-Links[END_REF].

The damage modeling of soft biological tissues and its numerical computation has been a topic of intense research. Holzapfel et al. [START_REF] Holzapfel | Changes in the mechanical environment of stenotic arteries during interaction with stents: computational assessment of parametric stent designs[END_REF], Gasser [START_REF] Gasser | Damage in vascular tissues and its modeling[END_REF] and Li [START_REF] Li | Damage Models for Soft Tissues: A Survey[END_REF] comprehensively reviewed the computational models for the damage of soft biological tissues. The developed damage models can be divided into three approaches [START_REF] Li | Damage Models for Soft Tissues: A Survey[END_REF]: (1) models based on continuum damage mechanics, (2) models based on the theory of pseudo-elasticity, and (3) the softening hyperelasticity approach. Among the continuum damage models, a gradient-enhanced non-local damage model was proposed by Dimitrijevic and Hackl [START_REF] Dimitrijevic | A method for gradient enhancement of continuum damage models[END_REF][START_REF] Dimitrijevic | A regularization framework for damage-plasticity models via gradient enhancement of the free energy[END_REF], and a nonlocal continuum healing model that combines the gradient-enhanced damage model and a temporally homogenized growth and remodeling model was originally presented in our previous works [START_REF] He | Gradient-enhanced continuum models of healing in damaged soft tissues[END_REF][START_REF] Zuo | Three-dimensional numerical simulation of soft-tissue wound healing using constrained-mixture anisotropic hyperelasticity and gradient-enhanced damage mechanics[END_REF].

This model has the following advantages: (1) good mesh independence is achieved in the simulation of damage evolution with growth and remodeling; and (2) the nonlocal damage process is realized by introducing the gradient-enhanced variable, thus allowing the effect of internal length scales of soft tissues to be considered.

Uncertainty widely exists due to the variability in human tissues, such as the patientdependent structural uncertainty and parametric uncertainty [START_REF] Paul | Quantifying the uncertainty in a hyperelastic soft tissue model with stochastic parameters[END_REF]. Parameter sensitivity analysis involving uncertainty is an essential step for developing and applying damage models to soft biological tissues. For instance, parameter sensitivity is the key information for medical device approval [START_REF] Gasser | Biomechanical modeling the adaptation of soft biological tissue[END_REF]. Some probabilistic method-based models have been presented to analyze the uncertainty of biomechanical problems. For instance, Laz and Browne [START_REF] Laz | A review of probabilistic analysis in orthopaedic biomechanics[END_REF] presented an overview of probability analysis in structural reliability, kinematics, joint mechanics, musculoskeletal modeling, and patient-specific representations. Regarding the uncertain damage analysis of soft tissues, Rodríguez et al. [START_REF] Rodríguez | A stochastic-structurally based three dimensional finitestrain damage model for fibrous soft tissue[END_REF] presented a model taking into account a stochastically-distributed waviness and a subsequent rupture of the individual fibers as the main driving force for damage evolution, and Schmidt et al. [START_REF] Schmidt | Statistical approach for a continuum description of damage evolution in soft collagenous tissues[END_REF] proposed a statistical approach to describe the evolution of microscopic damage in soft collagenous tissues.

More recently, Balzani et al. [START_REF] Balzani | Method for the quantification of rupture probability in soft collagenous tissues[END_REF] introduced a computational method for the assessment of rupture probabilities in soft collagenous tissues based on the numerical minimization and maximization of the probability of failure (PoF), which arise from random input quantities. Although the above models have been proposed to deal with the uncertainty analysis in biomechanics, there is still no work on the analysis of the uncertainty in the nonlocal damage model for soft biological tissues.

The large computational cost is a main challenge in uncertainty analysis when a probabilistic method is used [START_REF] Balzani | Method for the quantification of rupture probability in soft collagenous tissues[END_REF][START_REF] Lee | Propagation of uncertainty in the mechanical and biological response of growing tissues using multi-fidelity Gaussian process regression[END_REF]. Although the Monte Carlo method is the most straightforward way to conduct probability analysis [START_REF] Fishman | Monte Carlo: concepts, algorithms, and applications[END_REF], deterministic problems must usually be solved many times to guarantee the accuracy of the Monte Carlo method. Especially for the damage analysis for soft biological tissues, the computation of one deterministic problem usually requires considerable computational cost due to the solutions of nonlinear problems. Therefore, it is desirable to reduce the computational cost of uncertain damage analysis for soft biological tissues.

The surrogate model is an effective approach to reduce the computational cost of uncer-tainty analysis [START_REF] Ong | Evolutionary optimization of computationally expensive problems via surrogate modeling[END_REF]. The key idea of the surrogate model is to use an approximate model with rather a simple form to replace the original model by fitting the relationship between the input and output variables [START_REF] Gorissen | A surrogate modeling and adaptive sampling toolbox for computer based design[END_REF]. An overview of surrogate models can be found in the literatures provided by Díaz-Manríquez et al. [START_REF] Díaz-Manríquez | A review of surrogate assisted multiobjective evolutionary algorithms[END_REF], Razavi et al. [START_REF] Razavi | Review of surrogate modeling in water resources[END_REF], Asher et al. [START_REF] Asher | A review of surrogate models and their application to groundwater modeling[END_REF],

Zhang [START_REF] Zhang | A Review of Surrogate Models: 7-17[END_REF] or Raul et al. [START_REF] Yondo | A Review of Surrogate Modeling Techniques for Aerodynamic Analysis and Optimization: Current Limitations and Future Challenges in Industry[END_REF]. The widely-used surrogate models include the following: polynomial approximation models [START_REF] Myers | Response surface methodology: process and product optimization using designed experiments[END_REF], the radial basis function [START_REF] Hardy | Multiquadric equations of topography and other irregular surfaces[END_REF], Kriging models [START_REF] Sacks | Design and analysis of computer experiments[END_REF], etc. In this paper, the Legendre polynomial series surrogate model proposed by Wang et al. [START_REF] Wang | Collocation methods for uncertain heat convection-diffusion problem with interval input parameters[END_REF] is employed in the damage analysis of soft tissues due to its excellent error properties in the approximation of a globally smooth function [START_REF] Wang | On the convergence rates of Legendre approximation[END_REF].

To investigate the sensitivity of the uncertain material properties for the damage analysis of soft tissues, a numerical analysis approach is developed by integrating the gradientenhanced damage model [START_REF] Dimitrijevic | A method for gradient enhancement of continuum damage models[END_REF][START_REF] Dimitrijevic | A regularization framework for damage-plasticity models via gradient enhancement of the free energy[END_REF] and the Legendre polynomial series surrogate model [START_REF] Wang | Collocation methods for uncertain heat convection-diffusion problem with interval input parameters[END_REF] introduced above. The gradient-enhanced damage model is adopted to describe the behavior of the stress-softening of soft tissues and avoid mesh-dependence. The Legendre polynomial is used to replace the finite element method (FEM) solutions of deterministic problems, and the sparse grid collocation method (SGCM) is employed because of its advantage of using fewer collocation points when training surrogate models.

The remainder of this paper is organized as follows: Section 2 introduces the materials and methods, including the basic kinematics in Section 2.1, the gradient-enhanced damaged model in Section 2.2, the total potential energy and variational form in Section 2.3, the constitutive model in Section 2.4, the surrogate model in Section 2.5 and probability analysis in Section 2.6. Section 3 shows two numerical examples. Finally, discussions and conclusions are given in Section 4.

Methods and materials

Basic kinematics

Let x = ϕ(X, t) describe the motion of a body from the initial reference configuration X ∈ κ(0) to its current configuration x ∈ κ(t). The deformation gradient F and the Jacobian J that maps the referential volume dV onto the current volume dv are defined as

F = ∇ X ϕ, (1)
J = dv dV = det(F ). (2)

Gradient-enhanced damage model

The local strain energy function Ψ loc per unit of reference volume at each time is written as

Ψ loc = f (d)Ψ 0 , (3)
where Ψ 0 is the original (undamaged) strain energy density and f (d) represents a function of the damage variable d that measures the degree of material stiffness loss. In Eq. (3), f (d) = 1 means there is no damage and the total damage for f (d) → 0. It is at least twice differentiable and satisfies the following conditions:

f (d) : + → (0, 1]{f (0) = 1, lim f (d) = 0 d→∞ }. (4)
Following the approaches of Dimitrijevic and Hackl [START_REF] Dimitrijevic | A method for gradient enhancement of continuum damage models[END_REF][START_REF] Dimitrijevic | A regularization framework for damage-plasticity models via gradient enhancement of the free energy[END_REF], a gradient-enhanced nonlocal free energy term is added to the strain energy given in Eq. ( 3)

Ψ = f (d)Ψ 0 + c d 2 ∇ X φ 2 + β d 2 [φ -γ d d] 2 . (5)
In Eq. ( 5), c d represents the gradient parameter that defines the degree of gradient regularization and the internal length scales. Three other variables are introduced:

- and Hackl [START_REF] Dimitrijevic | A method for gradient enhancement of continuum damage models[END_REF][START_REF] Dimitrijevic | A regularization framework for damage-plasticity models via gradient enhancement of the free energy[END_REF] and Waffenschmidt et al. [START_REF] Waffenschmidt | A gradient-enhanced large-deformation continuum damage model for fibre-reinforced materials[END_REF].

Total potential energy and variational form

The general total potential energy for nonlocal damage model is

Π = Ω Ψ dV - Ω B • ϕ dV - ∂Ω T • ϕ dV, ( 6
)
where B is the body force vector per unit of reference volume of Ω and T is the traction on the boundary ∂Ω.

Minimization of the potential energy with respect to the primal variables ϕ and φ results in a coupled non-linear system of equations, which may be written as

Ω P : ∇ X δϕ dV - Ω B • δϕ dV - ∂Ω T • δϕ dV = 0, (7)
Ω Y : ∇ X δφ dV - Ω Y δφ dV = 0, ( 8
)
where P is the first Piola-Kirchhoff stress, Y is vectorial damage quantity related to flux terms and Y is the scalar damage quantity associated to source terms. They are defined as

P = ∂ F Ψ, Y = ∂ ∇ X φ Ψ, Y = -∂ φ Ψ. (9)
The corresponding spatial quantities in Eq. ( 9) and the body force vector B are given by

σ = P • cof(F -1 ), b = J -1 B, (10)
y = Y • cof(F -1 ), y = J -1 Y, (11)
where cof(F ) = JF -T .

Substituting Eqs. ( 2), [START_REF] Dimitrijevic | A method for gradient enhancement of continuum damage models[END_REF] and [START_REF] Dimitrijevic | A regularization framework for damage-plasticity models via gradient enhancement of the free energy[END_REF] into Eqs. ( 7) and ( 8), the variational forms of the spatial description are For the undamaged part Ψ 0 in Eq. (3), we used the Holzapfel-Gasser-Ogden (HGO) hyperelastic constitutive model [START_REF] Holzapfel | A New Constitutive Framework for Arterial Wall Mechanics and a Comparative Study of Material Models[END_REF][START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF], which can be written as

Ω σ : ∇ x δϕ dv - Ω b • δϕ dv - ∂Ω t • δϕ dv = 0, (12)
Ψ 0 = Ψ iso + Ψ aniso , (14)
where the isotropic part Ψ iso and anisotropic part Ψ aniso are

Ψ iso = 1 2 µ 0 J -2/3 (I 1 -3) + 1 2 κ 0 (J -1) 2 , (15)
Ψ aniso = k 1 2k 2 i=4,6 [exp(k 2 < E i > 2 ) -1], (16)
where µ 0 and κ 0 are the shear and bulk moduli of the soft isotropic matrix, respectively.

I 1 = tr(C)
is the first invariant of C, where C is the right Cauchy-Green tensor, which is defined as C = F T F . In Eq. ( 16), E i = κI 1 + (1 -3κ)I 4i -1 is a strain-like quantity and the

I 4i = C : A i = C : (a 0i ⊗ a 0i
) is the fourth invariant (which equals the square of the stretch of collagen fibers), the κ ∈ [0, 1/3] is a dispersion parameter and a 0i = (cosθ, sinθ, 0) is the unit direction vector of the fiber in the reference configuration. k 1 and k 2 are material parameters. The term < E i >, where

< E i >= (|E i | + E i )/2 is the Macaulay bracket,
reflects the basic assumption that fibers can only support tension.

Damage evolution

The detailed evolution of the damage variable d can be found in the works of Dimitrijevic and Hackl [START_REF] Dimitrijevic | A method for gradient enhancement of continuum damage models[END_REF][START_REF] Dimitrijevic | A regularization framework for damage-plasticity models via gradient enhancement of the free energy[END_REF] and Waffenschmidt et al. [START_REF] Waffenschmidt | A gradient-enhanced large-deformation continuum damage model for fibre-reinforced materials[END_REF]. In the following, we give a brief description of some basic and essential formulations.

Following standard thermodynamic considerations [START_REF] Dimitrijevic | A method for gradient enhancement of continuum damage models[END_REF][START_REF] Dimitrijevic | A regularization framework for damage-plasticity models via gradient enhancement of the free energy[END_REF], the damage conjugate driving force q is defined as

q = - ∂Ψ ∂d . ( 17
)
The damage condition at any time of the loading process is based on an energy release rate threshold condition [START_REF] Simo | Strain-and stress-based continuum damage models-I. Formulation[END_REF]

Φ d =q -r 1 ≤ 0, (18)
where r 1 represents the threshold value that triggers the evolution of the damage.

The differential equation of the evolution of damage is subjected to the Kuhn-Tucker optimality conditions [START_REF] Dimitrijevic | A method for gradient enhancement of continuum damage models[END_REF][START_REF] Dimitrijevic | A regularization framework for damage-plasticity models via gradient enhancement of the free energy[END_REF]

ḋ = κ∂Φ d ∂q , κ ≥ 0, Φ d ≤ 0, κΦ d = 0, ( 19
)
where κ is the Lagrange multiplier.

Verification of the damage model

An exemplary qualitative simulation is applied here to show how this non-local damage model can capture the stiffness loss of soft tissues when subjected to an overloading beyond physiological loads, and the results are shown in Fig. 1. In this simulation, an exponential damage function f (d) = e -d was adopted and was also used in following simulations. Other damage functions satisfying Eq. ( 4) could also be used [START_REF] Waffenschmidt | A gradient-enhanced large-deformation continuum damage model for fibre-reinforced materials[END_REF]. The chosen model was compared with experimental stress-strain responses of Raghavan et al. [START_REF] Raghavan | Ex vivo biomechanical behavior of abdominal aortic aneurysm: assessment using a new mathematical model[END_REF]. In these experimental data, thin strips of abdominal aortic aneurysm excised from the anterior surface during surgical repair were tested. Both the elastic and damage phases were well described by the non-local damage model as shown in Fig. 1.

Surrogate model 2.5.1. Legendre polynomial surrogate model

In this paper, the Legendre polynomial surrogate model is used to approximate the results of the FEM [START_REF] Wang | Collocation methods for uncertain heat convection-diffusion problem with interval input parameters[END_REF], and the Rodrigues formula of the N -order Legendre polynomial Φ(α) can be written as

Φ(α) = 1 2 N N ! d N d(α) ((α) 2 -1) N , ( 20
)
where α is the uncertain variable and N is the order of the polynomial series.

Based on the polynomial chaos framework [START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF], the output variable L(x, α) can be expressed as

L(x, α) = i w i (x)Φ i (α), ( 21
)
where x is the spatial coordinate, α is the uncertain variables vector, w i (x) is the corresponding expansion coefficient and i = (i 1 , i 2 , ..., i n ) stands for the multigauge. Note that L is the output variable, and it can be any variables concerned (stress and damage in this paper). In practice, a maximum order of polynomials N is selected to truncate Eq. ( 21) up to finite terms, by which Eq. ( 21) can be concisely rewritten as

L(x, α) ≈ L N (x, α) = i≤N w i (x)Φ i (α). ( 22
)

Sparse grid collocation method (SGCM)

The basic idea of a collocation method is using the response values at preselected points in uncertainty space to calculate the expansion coefficients of the surrogate model [START_REF] Wang | Collocation methods for uncertain heat convection-diffusion problem with interval input parameters[END_REF]. In this paper, the SGCM based on the Clenshaw-Curtis point-based collocation method and the Smolyak algorithm is adopted with the advantage that it can minimize the collocation error with the same number of nodes under the maximum norm [START_REF] Wang | Collocation methods for uncertain heat convection-diffusion problem with interval input parameters[END_REF].

For the univariate case, the Clenshaw-Curtis nodes distributed in the interval range [α i , α i ] are defined as

β i j = α c i α c i -cos π(j-1) m i -1 • ∆α i if m i = 1 if m i > 1 , (23)
where

α c i = (α i + α i )/2, ∆α i = (α i -α i )/2
, in which α and α represents the upper bound and lower bound of the uncertain interval, respectively, m i is the predefined number of collocation points.

For the multivariate case, the n input parameters are assumed to be independent with each other for simplification. The general approach is to extend the collocation points defined in one dimension in Eq. ( 23) to the entire interval uncertain space using a tensor product.

For each dimension, a nodal set

Θ m i i = {β i 1 , β i 2 , ..., β i m i } is selected.
On this basis, the entire nodal set named the full tensor product grids can be denoted as

Θ = Θ m 1 1 × Θ m 2 2 × . . . × Θ mn n , (24)
where the total number of collocation points M is

M = n i=1 m i = m 1 × m 2 × . . . × m i . (25)
To reduce the huge computational cost for high-dimensional problems, Smolyak-type sparse grids [START_REF] Smolyak | Quadrature and Interpolation Formulas for Tensor Products of Certain Classes of Functions[END_REF] are introduced, which are still based on the tensor product construction, but only a small subset of the full tensor grids Eq. ( 24) are selected. By using the Smolyak algorithm, the total nodal set of collocation points in Eq. ( 24) can be reconstructed in the following form:

Θ = l-n+1≤|i|≤l (Θ i 1 1 × Θ i 2 2 × . . . × Θ in n ), (26)
where l ≥ n is an integer representing the overall level of the construction and i j represents the separate level along the j th direction.

By iterating through all the nodes on full grids, only the ones whose sum order |i| = i 1 + i 2 + ... + i n across all dimensions is from ln + 1 to l are retained. To define the collocation points, an index k is introduced to denote the level of collocation points in onedimensional uncertainty space. The number m k i and positions β i j of the Clenshaw-Curtis nodes in the interval range [α i , α i ] are defined as

m k i = 1 2 k-1 + 1 if k = 1 if k > 1 β i j = α c i α c i -cos π(j-1) m k i -1 • ∆α i if k = 1 j = 1, 2..., m k i if k > 1 . (27)
After selecting the SGCM nodal set, a set Θ = {β 1 , β 2 , ..., β M } is given to specifically denote all the collocation points in the interval uncertain space, and L(x, β j ) is the solution of the objective variable in the deterministic problem at the preselected points β j that could be obtained by the FEM. Based on the surrogate function in Eq. ( 22), a group of linear equations with respect to the unknown expansion coefficients w i (x) can be derived as

      Φ 1 (β 1 ) Φ 2 (β 1 ) • • • Φ C n n+N (β 1 ) Φ 1 (β 2 ) Φ 2 (β 2 ) • • • Φ C n n+N (β 2 ) . . . . . . . . . . . . Φ 1 (β M ) Φ 2 (β M ) • • • Φ C n n+N (β M )            w 1 (x) w 2 (x) . . . w C n n+N (x)      =      L(x, β 1 ) L(x, β 2 ) . . . L(x, β M )      . ( 28
)
Once the polynomial bases Φ i and the SGCM nodal set Θ are given, the expansion coefficients w i (x) in Eq. ( 22) can be solved by Eq. ( 28). Note that the number of collocation points is required not to be smaller than the number of polynomial expansion terms, such as M ≥ C n n+N , to prevent the problem from being undetermined.

Probability analysis

The prediction of the rupture of soft tissues with uncertain parameters can be captured by the rupture probability, and PoF is defined as the probability that a mechanical rupture criterion holds. Such a criterion can be formulated as the inequality f > f max , where f represents a mechanical quantity of interest, for example, the stress, strain, or damage, which is considered relevant for initializing a rupture [START_REF] Balzani | Method for the quantification of rupture probability in soft collagenous tissues[END_REF].

Let R be a set of input parameters and z(R) is the mechanical quantity of interest. A suitable failure criterion can be defined as

z(R) ≥ z max , (29)
where z max is a maximally admissible value. Generally, the PoF can be defined as

PoF := P[z(R) ≥ z max ], (30)
where P(•) is the probability of (•).

In this paper, the beta distribution is used for uncertain parameters to avoid some possible unreasonable values of the distribution, that could be produced in a Gaussian distribution [START_REF] Schmidt | Statistical approach for a continuum description of damage evolution in soft collagenous tissues[END_REF]. The probability density function (PDF) of the beta distribution is defined as [START_REF] Johnson | Continuous univariate distributions[END_REF]:

PDF(α; a, b) = 1 B(a, b) α a-1 (1 -α) b-1 , ( 31
)
where α is the uncertain parameter, a and b denote two shape parameters and B(•) is the beta function. Note that some alternative distributions, e.g. the lognormal distribution [START_REF] Biehler | The impact of personalized probabilistic wall thickness models on peak wall stress in abdominal aortic aneurysms[END_REF],

could also be used for uncertain parameters.

In this paper, the damage function f (d) is used to describe the risk of the rupture, and a lower value of f (d) means a higher risk of rupture. Here, the PoF is redefined as

PoF := P[f (d) ≤ f (d) max ]. ( 32
)
The PoF can be calculated as the integral of the PDF of the damage function f (d) from

f (d) = 0 up to f (d) = f (d) max as PoF := P[f (d) ≤ f (d) max ] = f (d) max 0 PDF(f (d)) df (d). ( 33
)
An example is shown for the computation of the PoF in Fig. 2. The damage function is assumed to obey a beta distribution as f (d) ∼ B(4, 4) within the interval (0, 1] and the maximum damage is set to f (d) max = 0.5. With these assumptions, the value of the PoF can be calculated as 0.5 according Eq. ( 33), which is equaling to the area of gray part in Fig. 2.

Numerical examples

The gradient-enhanced damage model is implemented within the commercial finite ele-

Deterministic analysis

In this part, a series tests on 286 elements are performed, which the mesh-dependence had been already examined in our previous works [START_REF] He | Gradient-enhanced continuum models of healing in damaged soft tissues[END_REF][START_REF] Zuo | Three-dimensional numerical simulation of soft-tissue wound healing using constrained-mixture anisotropic hyperelasticity and gradient-enhanced damage mechanics[END_REF], to investigate the influence of the shear modulus and the internal length scales on the damage analysis in deterministic problem. The geometry, hyperelastic and damage parameters are reported in Table 1, in which the bulk modulus k 0 is set to 100 times of the shear modulus to make the material nearly incompressible. for Case 1 when the stretch is λ x = 0.1 are listed in Table 3. A good agreement can be seen in Table 3 such that the surrogate model with the SGCM can approximate the results of the FEM with the maximum difference being less than 1%. The good performance of the surrogate model with the SGCM is addressed again by comparing the results of the surrogate model with the SGCM and FEM at Node 1 when λ x = 0.1 for Case 2 and Case 3, and the results are shown in Fig. 8.

Probability analysis

In this example, shear modulus µ e and internal length scales c d are assumed to be uncertain parameters that obey the beta distributions in the given intervals (µ e ∈ [0.07, 0.13] MPa and c d ∈ [0.5, 1.5] MPa -1 • mm 2 ), and seven cases for different levels of uncertainty for µ e and c d listed in Table 2 are analyzed. The bulk modulus k 0 is set to 100 times of the mean value of the shear modulus to make the material nearly incompressible. Apart from µ e , c d and k 0 , other geometric and material parameters are shown in Table 1.

First, we investigate the computational cost for using the Monte Carlo method based on the direct FEM and the Monte Carlo method based on the surrogate model. When employing the Monte Carlo method, 10 5 samples are computed to ensure the accuracy of the probability analysis, as suggested by Wu et al. [START_REF] Wu | Robust fuzzy structural safety assessment using mathematical programming approach[END_REF]. The comparison of the CPU times is shown in Table 4. Note that the computational time for the Monte Carlo method with the direct FEM is estimated by multiplying the time of one FEM by 10 5 . A significant reduction of the computational cost can be observed in Table 4 such that the computing time of the Monte Carlo method with the surrogate model is about 0.06% of that using the Monte Carlo method with the direct FEM but the relative error is less than 1%. that a larger mean value of the input parameters leads to more serious damage and the distribution of the damage function will be more concentrated.

Radial dilatation of an idealized artery

The second example is damage analysis for radial dilatation of an idealized artery. The two-dimensional geometry and FEM mesh shown in Fig. 10 were previously established by Badel et al. [START_REF] Badel | Numerical simulation of arterial dissection during balloon angioplasty of atherosclerotic coronary arteries[END_REF] according to histological pictures of epicardial coronary arteries from Viles-Gonzalez et al. [START_REF] Viles-Gonzalez | Acute and chronic effects of epicardial radiofrequency applications delivered on epicardial coronary arteries[END_REF]. The artery is assumed to consist of a single medial layer containing an atherosclerotic plaque, and the balloon used for the angioplasty is modeled as a thin circular structure whose diameter increases during the angioplasty process. In this example, the damage is assumed to only occur in the artery. Based on the experiment data for human thoracic and abdominal aortas in Weisbecker et al. [START_REF] Weisbecker | Layer-specific damage experiments and modeling of human thoracic and abdominal aortas with non-atherosclerotic intimal thickening[END_REF], the shear modulus of the artery is assumed obeys a beta distribution µ e ∼ B(2.98, 8.72) within the interval [START_REF] Gasser | Biomechanical modeling the adaptation of soft biological tissue[END_REF]55] In this simulation, some simplifications are assumed, i.e., the influence of the residual stresses, smooth muscle, lipid pool, and the image-based layer-specific structure of the arteries are not considered, and an idealized geometric shape is used. The only boundary conditions to be assigned are the nodal displacements of the balloon. A radial displacement is imposed on each node from its initial position, D i = 1 mm, to give a final deformed diameter, D f = 2 mm, where the center of the balloon is fixed. In the following, we use an

inflation progress parameter λ = (D -D i )/(D f -D i ) × 100%
, where D is the current diameter of the balloon, to denote the inflation progress. Note that all the degrees of freedom of the balloon are prescribed as displacement boundary conditions, and the artery is inflated by the contact between the balloon and artery. The contact problems (balloon and plaque, balloon and artery, plaque and artery) are assumed to be frictionless following Badel et al.

[44] and are modeled as 'Hard contact' in Abaqus. The material parameters used in this example shown in Table 5 are according to Badel et al. [START_REF] Badel | Numerical simulation of arterial dissection during balloon angioplasty of atherosclerotic coronary arteries[END_REF] and Weisbecker et al. [START_REF] Weisbecker | Layer-specific damage experiments and modeling of human thoracic and abdominal aortas with non-atherosclerotic intimal thickening[END_REF].

Verification of the surrogate model

Probability analysis

First, the computational cost of probability analysis for using the Monte Carlo method with the direct FEM and the Monte Carlo method with the surrogate model are investigated in Table 7. The results again demonstrate that the surrogate model can effectively reduce the computational cost with good accuracy, as shown in Fig. 11. Obviously, the PoF of the artery increases as the balloon inflation size increases.

Discussion and conclusions

Due to the widely-existing uncertainty in damage analysis for soft biological tissues, the sensitivity analysis of nonlocal damage for soft biological tissues is an essential step in model development and application. To investigate the influence of the uncertain input parameters on the mechanical environment of the soft biological tissues, a numerical analysis is conducted by combining the gradient-enhanced damage model [START_REF] Dimitrijevic | A method for gradient enhancement of continuum damage models[END_REF][START_REF] Dimitrijevic | A regularization framework for damage-plasticity models via gradient enhancement of the free energy[END_REF] and the surrogate model-based probability analysis method [START_REF] Wang | On the convergence rates of Legendre approximation[END_REF].

The first advantage of the presented model is that it can assess nonlocal damage with uncertain material parameters, including internal length scales. In our previous works [START_REF] He | Gradient-enhanced continuum models of healing in damaged soft tissues[END_REF][START_REF] Zuo | Three-dimensional numerical simulation of soft-tissue wound healing using constrained-mixture anisotropic hyperelasticity and gradient-enhanced damage mechanics[END_REF], it was proved that the internal length scales have obvious effects on the localization of damage, e.g., a larger internal length scales lead to a larger lower level of damage and activated zone. In this paper, the influence of uncertain internal length scales can be analyzed by the proposed model. Therefore, it is possible to provide probabilistic results of the levels of damage, stress and deformation for soft tissue if the internal length scales are uncertain variables obeying a probability distribution. Moreover, the proposed model can consider the combined effects of different uncertain material parameters. For instance, both the shear modulus and an internal length scales parameter can be set as uncertain variables at the same time. This will be useful for dealing with the cases in which more uncertain parameters exist together.

The other advantage is that this paper's method provides a new approach for the uncertain damage analysis of soft biological tissues based on the surrogate model. Although the Monte Carlo method is an accurate and robust method for uncertainty analysis, the huge computational cost is a key challenge for its application to damage analyses since the nonlinear problem has to be repeatedly solved by the FEM [START_REF] Neumaier | Interval methods for systems of equations[END_REF][START_REF] Ran | A gradient based algorithm to solve inverse plane bimodular problems of identification[END_REF]. Due to the huge computational cost, the Monte Carlo method is commonly introduced as a referenced approach, but it is rarely used in the practical engineering. In this paper, the surrogate model based on the Legendre polynomial series was adopted to approximate the results of the FEM. When training of the surrogate model, the SGCM is employed to accelerate the procedure. First, it is demonstrated that the surrogate model has good accuracy compared with direct FEM solutions. For example, the maximum error is less than 1% compared with the results (Cauchy stress and the damage function) obtained from the surrogate model and the FEM at some feature nodes in the open-hole plate. Second, the surrogate model significantly improved the computational efficiency of the Monte Carlo method. For example, in the open-hole plate problem, the computational cost for the Monte Carlo method with the direct FEM was about 40 days, but the time for the surrogate model with the SGCM including training the surrogate model and probability analysis was only 0.63 h with a 286 elements mesh.

The proposed approach was applied to the simulation of the radial dilatation of an idealized artery by considering the uncertain shear modulus and internal length scales. Although some simplifications are assumed, some interesting results are found: (1) the spatial positions and level of damage could have significant impact on the uncertainty propagation of damage, for instance, the mean value and standard deviation are obviously different at different locations (Node 1, Node 2 and Node 3) in Fig. 13; (2) different parameters could also have obviously different impact on the uncertainty propagation of damage. For instance, if only the uncertain shear modulus is considered (Case 1), a smaller standard deviation of damage is produced at Node 1 with larger damage compared with Node 3, while the internal length scales has an opposite trend with a larger standard deviation at Node 1 compared with Node 3 in Case 2; (3) although the PoF is increased with the increase of inflation size as expected, the relation between uncertain parameter and resulted PoF could still provide a potential helpful approach for both directly predicting the PoF in clinical operation and inversely choosing the inflation size under a prescribed PoF.

The main limitation in this paper is that only 2D models and an idealized geometric model was considered. The development of a 3D patient-specific model is currently under progress in order to address more realistic applications. Meanwhile, some important aspects, such as the influence of the residual stresses, the smooth muscle activation, the lipid pool should also be considered. Besides, some uncertainty analyses about the geometry and the material parameters of the plaque, which is important in predicting the rupture risk of the artery, should be conducted, and the spatial uncertainty should also be considered as arteries are usually heterogeneous materials [START_REF] Bersi | Novel Methodology for Characterizing Regional Variations in the Material Properties of Murine Aortas[END_REF][START_REF] Bersi | Multimodality Imaging-Based Characterization of Regional Material Properties in a Murine Model of Aortic Dissection[END_REF]. Moreover, there is still a pressing need to develop further experiments to identify the distribution of shear moduli, the internal length scales and other unknown material parameters. Furthermore, the use of UEL presents some limitations such as, for instance, the definition of slave surfaces in contact analyses.

In summary, we have coupled a nonlocal damage model and a surrogate model to investigate the influence of uncertainty on damage in soft tissues. The Legendre polynomial surrogate model with the SGCM was adopted in this paper to speed up the computation efficiency. The evolution of the PDF of the damage function f (d) was illustrated through two examples where the input parameters were assumed to be uncertain, and the method's clinical application to balloon angioplasty was also addressed. Future work will consist of applying the proposed model to practical problems by collecting clinical data.

  the field variable φ, which transfers the values of the damage parameter across the element boundaries to make it non-local by nature; -the energy-related penalty parameter β d , which approximately forces the local damage field and the nonlocal field to coincide; and -parameter γ d , which is set as a switch between the local and enhanced models. When γ d = 0 and c d = 0, it leads to a local model, while setting γ d = 1 and c d = 0 leads to the enhanced non-local model. A typical value of this parameter is set to 1 following Dimitrijevic

Ωy: 1 .

 1 ∇ x δφ dv -Hyperelastic part of the free energy

  ment software Abaqus/Standard by means of a user subroutine UEL and the simulation of the probability analyses based on the surrogate model are carried out by Matlab R2018a. Both computations are conducted on a PC with an Intel Xeon E5-2650 CPU @ 2.40 GHz and 32 GB of RAM. The first example is used to illustrate the nonlocal damage model and the advantages of the surrogate model with the SGCM through an open-hole plate case. The second example is a radial dilatation of an idealized artery, which is used to treat atherosclerotic plaque. Note that only homogeneous uncertain material parameters without spatial correlation are assumed for simplification, and plane strain elements (Bilinear 8 node element CPE8 in Abaqus) are used in all simulations. 3.1. Open-hole plate The first numerical example is an open-hole plate under displacement loading. The geometry and the FEM model are shown in Fig. 3. Due to the symmetry, only one-fourth of the plate is analyzed. In this example, only the isotropic contribution of the soft tissues is considered. The result of the deterministic analysis are shown in Section 3.1.1, the verification of the surrogate model is shown in Section 3.1.2 and the results of the probability analysis are shown in Section 3.1.3.

3 . 1 . 2 .

 312 The contours of the damage function f (d) for different shear modulus are shown in Fig.4, while the internal length scales is c d = 1.0 MPa -1 • mm 2 . It can be obviously seen that a larger value of the shear modulus resulting a larger damage, and the same results can be also observed in Fig. 5. The value of the parameter c d represents the degree of regularization and it depends on the microstructure of the tissues. The influence of the internal length scales c d is investigated by three different values and the results are shown in Figs. 6 and 7, while the shear modulus is a constant that µ e = 0.1 MPa. As shown in Figs. 6 and 7, the levels of damage are strongly depending on the value of c d . It can be explained by the width of the activated zone [10, 11, 34]. With a smaller c d , a smaller region is influenced by the damage variable, resulting in more concentrated distribution of the damage. Verification of the surrogate model In the training of the surrogate model, the 5-order Legendre polynomial series and the SGCM with k = 4 are used. The comparison of the values of f (d) and σ x computed by the surrogate model and the FEM at six feature nodes (the locations are shown in Fig. 3(b))

  Fig. 9(b), Node 1 is the location where the damage is largest and the distribution of f (d) varies in a smaller range (the standard deviation is σ = 0.030). Comparing Nodes 1, 2 and 3, it can be found that the location with less damage is more sensitive to the uncertain parameters, and the sensitivity increases as the damage level increases such that the standard deviation σ of the damage function f (d) increases from 0.03 to 0.087 when the mean value μ of f (d) increases from 0.066 to 0.41. Finally, the levels of uncertainty for parameters are also investigated and the results are shown in Figs. 9(c) and 9(d). Cases 4-6 are tested to illustrate the influence of the standard deviation σ when the mean value μ is the same. It can be seen in Fig. 9(c) that when the mean value of the uncertain parameters is constant, the mean value of the damage function f (d) changes slightly, and the standard deviation of f (d) increases as the standard deviation of the uncertain input parameters increases. The influence of the mean value of the input parameters is also investigated in Fig. 9(d). It can be seen that the change of the mean value of the input parameters may affect the distribution of the damage function f (d) such

  kPa. The internal length scales c d is assumed to obey a beta distribution c d ∼ B(4, 4) within the interval [0.005, 0.015] kPa -1 • mm 2 .

  Similar with the previous example, three different cases, including (1) µ e is uncertain with a constant c d , (2) c d is uncertain with a constant µ e and (3) both µ e and c d are uncertain, are tested and the distributions of the uncertain input parameters are shown in Table 6. The 5-order Legendre polynomial surrogate model is used and the SGCM is set by setting k = 5. To check the accuracy of surrogate model, the damage function f (d) is calculated by the surrogate model and the direct FEM for 3 cases, and the results shown in Fig. 11 illustrate the good performance of the surrogate model. The contours of the damage function f (d) calculated by the surrogate model with SGCM and FEM for different shear modulus and the internal length scales when λ = 100% in Fig. 12 demonstrate the good performance of the surrogate model once again.

  Fig. 14(b). It can be seen that the PoF of f (d) is directly influenced by the inflation progress parameter λ. For instance, the PoF are 0.03 and 0.91 for λ = 90% and λ = 95%, respectively.

28 2 29 3 30 4 31 5 32 6 33 7 34 8 35 9Figure 1 :Figure 2 :

 282930313233343512 Figure 1: Qualitative simulation based on the damage model used in this paper and the experimental stressstrain data of the abdominal aortic aneurysm in Raghavan et al. [38]. Material parameters are chosen as µ e = 0.3 MPa, κ 0 = 100 MPa, k 1 = 0.22 MPa, k 2 = 1.20, θ = 50 deg, κ = 0.18, r 1 = 0.05 MPa.

Figure 3 :Figure 4 : 15 Figure 5 :Figure 6 :c d = 5 Figure 7 :Figure 8 :Figure 9 :

 3415565789 Figure 3: Reference model for an open-hole plate. (a) Geometric, dimension and boundary conditions; and (b) FEM mesh.

Figure 10 :Figure 11 :

 1011 Figure 10: Reference model in the radial dilatation of an idealized artery. (a) Geometry, dimensions and boundary conditions; and (b) FEM mesh.

Figure 12 :

 12 Figure 12: The contours of the damage function f (d) obtained by the surrogate model with the SGCM and FEM for different shear modulus and the internal length scales when λ = 100%. (a) µ e = 15 kPa, c d = 0.01 kPa -1 • mm 2 , (b) µ e = 55 kPa, c d = 0.01 kPa -1 • mm 2 , (c) µ e = 35 kPa, c d = 0.005 kPa -1 • mm 2 , and (d) µ e = 35 kPa, c d = 0.015 kPa -1 • mm 2

Table 1 :

 1 Geometry, hyperelastic and damage parameters used in the open-hole plate.

	Type	Description	Symbol	Values	Units
		Height	H	200	mm
	Geometric	Width	W	200	mm
		Radius	R	50	mm
	Hyperelastic	Shear modulus Bulk modulus	µ 0 κ 0	[0.5,1.0,1.5] [5,10,15]	MPa MPa
		Damage threshold	r 1	0.001	MPa
	Damage	Internal length scales Penalty parameter	c d β d	[0.5,1,5] 0.005	MPa -1 • mm 2 MPa -1
		(Non)local switch	γ d	1	-

Table 3 :

 3 Comparison of the results obtained by the FEM and surrogate model with the SGCM when λ x = 0.1 for Case 1.

			FEM		SGCM
	Node	µ e	f (d)	σ x	f (d)	error	σ x	error
	1	0.08 0.12	0.142469 0.028766	0.028497 0.026306	0.141972 0.028787	0.35% 0.07%	0.028479 0.026207	0.06% 0.38%
	2	0.08 0.12	0.38338 0.161063	0.031512 0.031508	0.382846 0.161157	0.14% 0.06%	0.031431 0.03152	0.26% 0.04%
	3	0.08 0.12	0.607423 0.293042	0.014506 0.015549	0.605692 0.294017	0.28% 0.33%	0.014729 0.015552	0.05% 0.02%
	4	0.08 0.12	1 1	0.001851 0.001331	1 1	0.00% 0.00%	0.001855 0.001342	0.23% 0.80%
	5	0.08 0.12	1 1	-0.00159 -0.00086	1 1	0.00% 0.00%	-0.00159 -0.00086	0.10% 0.20%
	6	0.08 0.12	1 1	0.022613 0.021469	1 1	0.00% 0.00%	0.022606 0.02149	0.03% 0.10%

Table 4 :

 4 The comparison of the computational cost between using the Monte Carlo method with the direct FEM and the Monte Carlo method with the surrogate model.

		Monte Carlo method with the direct FEM (Estimated, t 1 )	Monte Carlo method with the surrogate model Training the Probability Total time surrogate model analysis (t 2 )	t 2 /t 1
	79 elements	14.70 days	825 sec	11.82 sec	836.82 sec	0.06%
	286 elements	40.05 days	2249 sec	36.28 sec	2285.28 sec	0.07%
	793 elements	115.27 days	6474 sec	98.73 sec	6572.73 sec	0.07%

Table 5 :

 5 Hyperelastic and damage parameters used in the radial dilatation of an idealized artery[START_REF] Badel | Numerical simulation of arterial dissection during balloon angioplasty of atherosclerotic coronary arteries[END_REF][START_REF] Weisbecker | Layer-specific damage experiments and modeling of human thoracic and abdominal aortas with non-atherosclerotic intimal thickening[END_REF].

	Type	Part	Description	Symbol Values	Units
			Shear modulus	µ e	uncertain kPa
			Bulk modulus	κ e	5	MPa
		Medial layer	Material parameter Material parameter	k 1 k 2	5.15 8.64	MPa -
	Hyperelastic	Plaque	Dispersion parameter Fiber orientation angle θ κ 0 Shear modulus µ p Bulk modulus κ p	0.24 ±38.8 20 34	-deg kPa kPa
		Balloon	Shear modulus Bulk modulus	µ b κ b	10 10	kPa kPa
			Damage threshold	r 1	1.0	kPa
	Damage	Medial layer	Internal length scales Penalty parameter	c d β d	uncertain kPa -1 • mm 2 10 kPa -1
			(Non-)local switch	γ d	1.0	-

Table 7 :

 7 The comparison of the computational cost between using the Monte Carlo method with the direct FEM and the Monte Carlo method with the surrogate model.

	Monte Carlo method with	Monte Carlo method with the surrogate model
	the direct FEM (Estimated t 1 )	Training the surrogate model	Probability analysis	Total time (t 2 )	t 2 /t 1
	230.09 days	19422 sec	35.43 sec	19457.43 sec	0.09%

Acknowledgments

The research leading to this paper was funded by the NSFC Grant [12072063], ERC-2014-CoG-BIOLOCHANICS [647067], grants from the State Key Laboratory of Structural Analysis for Industrial Equipment [S18402, GZ19105], the Liaoning Provincial Natural Science Foundation [2020-MS-110] and the Visiting Scholar Foundation of Key Laboratory of Biorheological Science and Technology (Chongqing University) from Ministry of Education [CQKLBST-2020-002].


Baptiste Pierrat
email: pierrat@emse.fr

Vít Novácěk

Frédéric Stéphane Avril

F Turquier

Mechanical characterization and modeling of knitted textile implants with permanent set

Keywords: abdominal wall hernia, knitted textile-based implant, constitutive modeling

come

Mechanical characterization and modeling of knitted textile implants with permanent set

Introduction

Ventral incisional hernia (VIH) is the protrusion of abdomen contents through a weakness of the abdominal wall at an incision closure site. VIH repair is one of the most common surgical procedures worldwide. Despite the general use of textile based implants over the past two decades, complications have remained very high. Recurrence rates have been reported close to 11% and over 44% respectively for general and high risk patient populations [START_REF] Cassar | Surgical treatment of incisional hernia[END_REF][START_REF] Bosanquet | Systematic review and meta-regression of factors affecting midline Incisional hernia rates: Analysis of 14 618 Patients[END_REF]. The average hospital costs associated with ventral and incisional hernia repair for a patient have been estimated to 10 700$ in the US [START_REF] Plymale | Ventral and incisional hernia: the cost of comorbidities and complications[END_REF] and 6 451€ in France [START_REF] Jf Gillion | The economic burden of incisional ventral hernia repair: a multicentric cost analysis[END_REF]. This should be put in perspective with the very high number of interventions, around 350 000 annually in the US [START_REF] Smith | Ventral hernia[END_REF] and 65 000 in France [START_REF] Jf Gillion | The economic burden of incisional ventral hernia repair: a multicentric cost analysis[END_REF]. The absence of consensus has led to tremendous variation in treatment outcomes, costs and resource utilization. Moreover, VIH repair suffers drastically from a lack of treatment personalization. Numerical modeling and simulation have appeared as one of the most promising approach to investigate biomechanics related to VIH repair, to identify influential parameters, to optimize and customize textile implants in combination with surgical techniques, and finally, to improve clinical outcomes. Surgical textiles are highly conformable and porous structures that result from the knitting of fibers. Their mechanical behavior and strength depend upon the fiber diameter, the material composition, the effect of the manufacturing process, the knitting patterns, the shapes and dimensions, and finally upon the loading conditions. Surgical textiles exhibit an anisotropic non-linear behavior with high (nonlinear) Poisson's ratio. Hysteresis is also commonly observed when submitting these textiles to uniaxial and planar biaxial tensile loads [START_REF] Deeken | Biaxial analysis of synthetic scaffolds for hernia repair demonstrates variability in mechanical anisotropy, non-linearity and hysteresis[END_REF][START_REF] Deeken | Mechanical properties of the abdominal wall and biomaterials utilized for hernia repair[END_REF][START_REF] Cordero | Biaxial Mechanical Evaluation of Absorbable and Nonabsorbable Synthetic Surgical Meshes Used for Hernia Repair: Physiological Loads Modify Anisotropy Response[END_REF][START_REF] Est | Multi-directional mechanical analysis of synthetic scaffolds for hernia repair[END_REF].

Computational simulations of such complex behavior have been attempted at the fiber scale [START_REF] Duhovic | Simulating the deformation mechanisms of knitted fabric composites[END_REF][START_REF] Liu | On the role of material architecture in the mechanical behavior of knitted textiles[END_REF][START_REF] Mckee | Finite element simulation of ballistic impact on single jersey knit fabric[END_REF][START_REF] Dinh | Prediction of mechanical properties of knitted fabrics under tensile and shear loading: Mesoscale analysis using representative unit cells and its validation[END_REF]. These models have revealed a very high potential as they can take into account the impact of the knitting process on the fibers (geometry and material properties), deal with fiber-fiber interactions and open the way for multi-scale modeling (micro/meso/macro). However, due to prohibitive computational costs, only the standard Jersey fabric has been considered so far. This simple pattern is too far apart from the ones of regular surgical textiles used in VIH repair.

Numerous homogenized models have also been developed at the macroscopic scale. Textiles are then assumed as homogenous and continuous media and modeled using the theory of hyperelasticity [START_REF] Yeoman | A constitutive model for the warp-weft coupled non-linear behavior of knitted biomedical textiles[END_REF][START_REF] Hernández-Gascón | Mechanical response of the herniated human abdomen to the placement of different prostheses[END_REF][START_REF] Morch | Modeling of anisotropic hyperelastic heterogeneous knitted fabric reinforced composites[END_REF]. Their constitutive equations derive from strain energy functions developed initially for biological [START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF] or rubber-like soft tissues [START_REF] Arruda | A three-dimensional constitutive model for the large stretch behavior of rubber elastic materials[END_REF]. Although anisotropic, their characterization remains exclusively based on monotonous uniaxial tensile tests. Therefore, there is an urgent need to carry out multi-axial cyclic tests reproducing physiological loading conditions in order to calibrate computational models.

In this study, advanced uni-axial and multi-axial experiments were conducted to identify and validate a novel anisotropic hypo-elastoplastic constitutive model specifically established for surgical knitted textile implants, including all the major mechanical characteristics of these materials such as anisotropy and permanent set.

cm 1 mm 5 mm

Figure 1: Illustration of the different scales and complex structure of knitted hernia meshes.

Materials and Methods

Constitutive modelling

Surgical knitted textiles were modeled using the membrane theory for very thin shells, neglecting bending and twisting moments. A bi-dimensional plane stress formulation was adopted, relating the vector made of the 3 in-plane components of the logarithmic strain tensor {ε} = {ε 11 ; ε 22 ; ε 12 } to the vector made of the 3 in-plane components of the Cauchy stress tensor {σ} = {σ 11 ; σ 22 ; σ 12 }, where directions 1 and 2 are the local warp and weft directions of the knitted textile respectively. Anisotropy, non-linearity and permanent set of the surgical knitted textiles were captured with an original anisotropic hypo-elastoplastic membrane model.

The constitutive model was based on an additive split of the logarithmic strain tensor into a plastic strain {ε pl } representing the permanent set of the mesh, and an elastic strain {ε el } such as {ε} = {ε el } + {ε pl }.

The stress rate is obtained using an incremental formulation such as:

{ σ} = [C(ε el )]{ εel } = g(ε)[C(ε el )]{ ε} (1)
where { εel } is the elastic strain rate, { ε} is the total strain rate, { σ} = { σ} + {Ω.σ} -{σ.Ω} is the Green-Naghdi rate of the Cauchy stress and Ω is the spin tensor.

Despite the anisotropy of the surgical textiles, it is assumed (assumption verified in the Results section, see Fig. 2) that there exists a unique scalar

function f (ε ij ) = a(εij ) 2
(εij ) 2 +bεij +c such as ε el ij = f (ε ij ) and:

g(ε ij ) = ∂f ∂ε ij = εel ij εij (2)
[C(ε el )] is the strain-dependent tangent stiffness matrix which may be written for a membrane such as:

[C(ε el )] =      1 E 0 1 +E 1 1 ε el 11 -ν12 E 0 1 +E 1 1 ε el 11 0 -ν12 E 0 1 +E 1 1 ε el 11 1 E 0 2 +E 1 2 ε el 221 0 0 0 1 G 0 12 +G 1 12 ε el 12      -1 (3)
Eventually, there are 10 parameters to identify in the constitutive model:

• The 6 stiffness constants E 0 1 , E 1 1 , E 0 2 , E 1 2 , G 0 12 and G 1 12 ,
• the Poisson's ratio ν 12 ,

• the 3 parameters of the yield evolution function a, b and c.

The methodology used to identify these parameters is described in Section 2.3.

Mechanical characterization

Uniaxial and equibiaxial quasi-static cyclic tests were carried out in a 37°C

water bath, on specimens presenting a gauge area of 25x25 mm 2 . Uniaxial quasi-static cyclic tests were carried out with an Instron © ElectroPuls E3000 system fitted with a 250 N load cell. Tension was applied in 3 different directions of the textile: warp, weft and diagonal (45°). The machine was controlled in displacement with a jaw velocity of 0.25 mm/s. Cyclic displacement ramps were performed. Each stage consisted of 5 cycles, with a loading reaching x% maximum strain and an unloading to 0 N. The maximum strain applied at each stage x was successively increased, taking the following values: 5%, 10%, 15%, 20%, 25%.

A TA Instruments ElectroForce ® planar biaxial test machine was used to carry out the equibiaxial tensile tests and assess the mechanical coupling between warp and weft directions. The force was measured in each direction with two 225 N load cells and the machine being controlled in load. Ten cycles were applied between 0 and 40 N on each sample, inducing a tension of 16 N/cm similar to the effects of intra-abdominal pressure loading reported by Klinge et al. [START_REF] Klinge | Modified mesh for hernia repair that is adapted to the physiology of the abdominal wall[END_REF].

For both uniaxial and biaxial tensile tests, logarithmic strains were assessed from the length separating the holding jaws ε = log( L L0 ) where L was the current length and L 0 the initial length.

The Cauchy stress values were calculated assuming a thickness t of 0.5 mm.

Note that this thickness is an average value at the macroscopic scale where the textile is modelled as a continuum medium. The Cauchy stress was derived as σ = F t w where F is the measured load and w is the current width measured using a camera and the VIC-3D TM (Correlated Solutions) digital image correlation system.

Identification technique

After performing the experiments, the obtained experimental stress-strain curves (see Figure 4) were separated into two regions for post-processing:

• the envelope, i.e. the part when the plastic strain increases with the load,

• the cycles, i.e. the unloading and reloading part obtained after the first loading of each stage. During the cycles, the plastic strain remains constant.

Tangent moduli in warp and weft directions

The unloading and reloading stage of the last cycle was extracted and averaged (neglecting hysteresis effects) for the warp and weft directions in the uniaxial tensile tests. It is assumed that during this stage the response is purely elastic (Figure 5a). A 2 nd order polynomial function was fitted to these elastic stress-strain curves with a least-squares algorithm implemented in Matlab ® (Figure 5a), yielding the 4 stiffness constants E 0 1 , E 1 1 , E 0 2 and E 1 2 .

Yield evolution function

Under uniaxial tension, the accumulated permanent strain after unloading depended on the maximum strain value which was previously reached. The evolution law was similar regardless of the loading direction (Figure 2). This evolution was determined by plotting the elastic strain as a function of total strain for the envelope, using the previously determined elastic stress-strain response: at a given σ, (i) the elastic strain was determined from Figure 5a and (ii) the total strain was determined from the envelope. Then, the elastic strain was plotted against the total strain for warp and weft directions, again showing an isotropic behavior in this response. A 2 nd order rational polynomial was used to fit these data points with a least-square algorithm in Matlab ® (Figure 5b), The permanent set was derived as the residual strain obtained after unloading. The curves show that the evolution of the permanent set is similar in the warp, weft and diagonal directions.

Poisson's ratio

Equibiaxial stress-strain curves were used to identify the Poisson's ratio ν 12 , after identifying the 4 stiffness constants E 0 1 , E 1 1 , E 0 2 and E 1 2 . A linear leastsquare fitting was used to determine ν 12 from the stress-strain response curves of equibiaxial experiments (Figure 4d). Note that uniaxial experiments were not used for this identification; the free boundary conditions on the edges made the textile curl slightly, over-estimating the effective lateral contraction.

Shear modulus

The diagonal tensile experiments could not be used directly for identifying G 0 12 and G 1 12 , as this experiment involves a mix of tension and shear. Another issue with this experiment is that it is not statically determined for finite deformations (there is no analytical expression of the local stress as a function of the applied load). Consequently, an inverse problem must be formulated and solved. A single finite-element was coded in Matlab ® with the previously identified constitutive parameters, and submitted to the diagonal tensile test. At every strain increment, the tangent modulus G 12 (ε el 12 ) was identified to fit the diagonal stress-strain response. Then, this function was approximated with a linear function to derive G 0 12 and G 1 12 .

Abaqus implementation

The models were implemented and solved using the commercial finite-element (FE) software Abaqus ® /Explicit.

The constitutive behaviour relating the strain rates to the Cauchy stress rates is the core of the explicit solver and was implemented as a VUMAT subroutine in Abaqus ® /Explicit. The following pseudo-code was used:

• from the current strain increment, update the total strain ε and the maximum strain ε max ,

• using f (ε), calculate the elastic strain:

-if ε i ≥ ε max i : ε el i = f (ε i ), -else ε el i = ε i -ε pl i
where ε pl i is kept the same as in the previous increment,

• build the stiffness matrix [C] using E 0 1 , E 1 1 , E 0 2 , E 1 2 , G 0 12 , G 1 12 and ν 12 ,
• compute the Cauchy stress increment and update the total Cauchy stress.

For the sake of verification, the implemented subroutine was assigned to a single element of dimensions 1 × 1 × 1 mm 3 and submitted to the loading history of the uniaxial and biaxial experiments. The stress-strain responses of the model were successfully compared to those of the experiments.

Validation

Previously acquired functional experimental data [START_REF] Belzacq | Assessment of shear forces distribution at fixation points of textile based implants[END_REF] were used to validate the constitutive models. The experimental set-up consisted of a 100 mm diameter hemispherical plunger loading a 250 mm diameter circular mesh through a 110 mm diameter circular patch, with fixation points represented by metal rods fitted with markers (more details available in [START_REF] Belzacq | Assessment of shear forces distribution at fixation points of textile based implants[END_REF]). The distribution of reaction forces at the mesh fixation points was determined via rod continuum mechanics and the positions of markers were recorded by marker tracking sys-

tem. An FE quarter model of this experiment was built in Abaqus ® . This

Results

Tensile test results

The response of the mesh subjected to tensile loading exhibited an anisotropic, non-linear behaviour (Figure 4). The warp direction is the stiffest, followed by the diagonal and weft directions. Under biaxial tension, the mesh appeared much stiffer, due to the strong coupling between the directions. The cyclic tests showed a permanent set after unloading, which increased when the maximum applied strain increased (although it also slightly increased with the number of cycles, but to a minor extent). Note that further testing showed that this permanent set was not due to viscoelastic effects, as the residual strain was permanent and did not relax with time. A strong hysteresis was also observed during the unloading-reloading cycles.

Model identification

The polynomial functions fitted to the response curves of the uniaxial tensile tests are shown in Figure 5a. They are in very good agreement with the 2 nd order polynomial functions.

Regarding the yield evolution function, the ratio between the elastic strain and the total strain was the same in the weft and in the warp direction:

εel 11 ε11 = εel 22 ε22
, confirming the isotropic nature of this phenomenon (Figure 5b).

All the identified material parameters are reported in

Material parameter Value Unit

E 0 1 13.2 [MPa] E 1 1 1205 [MPa] E 0 2 13.3 [MPa] E 1 2 187 [MPa] G 0 12 0.87 [MPa] G 1 12 148 [MPa] ν 12 0.86 [-] a 0.12 [-] b 0.089 [-] c 0.

[-]

Table 1: Identified material parameters for the hypo-elastoplastic constitutive model.

Model evaluation

The response curves predicted with a single finite-element to reproduce the uniaxial and biaxial cyclic tests are shown in Figure 6. As expected, the model captures the non-linearity, anisotropy and permanent set of the mesh behavior. Both uniaxial and biaxial responses are well reproduced. The unloadingreloading hysteresis is neglected, as a purely elastic behavior was assumed without viscous effects. The prediction of the model with respect to the experiments was quantified using the coefficient of determination (R 2 ). A satisfactory value of 0.94 was obtained for the overall data, despite hysteresis effects not reproduced by the model.

Validation results

Static equilibrium (defined as ratio of kinetic energy to external work lower than 1%) was reached with the FE model up to a 231 N force. Note that this is equivalent to a 182 mmHg intra-abdominal pressure, that is just above the mean pressure reported by Cobb et al. [START_REF] Cobb | Normal intraabdominal pressure in healthy adults[END_REF] for jumping in place (the most severe activity of his study). The strain field at this maximum force is shown in Figure 7a whereas the strain field after unloading is shown in Figure 7b.

Strain concentrations are visible near the fixations, reaching values as high as 37%, which is out of the strain range covered in the experimental tests. The region under the plunger exhibits strains in the 10-20% range. After unloading, a residual strain field is present, which leads to a residual bulge as seen in Figure 7b.

Fixation forces were measured at an intermediate value of 60 N force, and were compared against their counterpart predicted by the FE model (Figure 8).

A very good agreement can be noted, with a mean difference of 0.55 N. The numerical model successfully captured the heterogeneous force distribution, with higher forces along the warp direction.

Finally, the global force vs. displacement response curve in the plunger experiment was compared to the FE predictions (Figure 9). As previously mentioned, the FE model could satisfy mechanical equilibrium up to 231 N,

Discussion

In this work, we have developed and validated an original constitutive model for knitted textiles which takes into account the major features needed to simulate numerically VIH interventions using these meshes: large deformations, non-linearity, anisotropy, and permanent set.

To our best knowledge, this is the first constitutive model of a knitted textile including all these features. Several homogeneous constitutive models of knitted textile implants were previously proposed in simulations of abdominal hernia repair [START_REF] Hernández-Gascón | Mechanical response of the herniated human abdomen to the placement of different prostheses[END_REF][START_REF] Morch | Modeling of anisotropic hyperelastic heterogeneous knitted fabric reinforced composites[END_REF][START_REF] Tomaszewska | Computer Methods, Imaging and Visualization in Biomechanics and Biomedi-cal Engineering[END_REF][START_REF] Yeoman | A constitutive model for the warp-weft coupled non-linear behavior of knitted biomedical textiles[END_REF]. They were based on anisotropic hyperelasticity and the material parameters were identified using uniaxial tensile tests. Effective porosity of the textile was previously taken into account [START_REF] Horbach | Modelling of compressible and orthotropic surgical mesh implants based on optical deformation measurement[END_REF]. Although the authors mention residual strains and corresponding experimental data in [START_REF] Morch | Modeling of anisotropic hyperelastic heterogeneous knitted fabric reinforced composites[END_REF],

their constitutive model did not account for permanent set.

Structural effects at the scale of pores were neglected in our continuum mechanics approach. The experiments used to fit material coefficients were done on small specimens with dimensions of about 10 pores (25 mm). It was shown that the continuum assumption is validated for predicting reaction forces at the fixation points. Although all our sample dimensions were designed to be at least 10 times larger than the unit cell size of the mesh (pores), scale effects may still exist and future work will also investigate the mechanical response of larger specimens. Another limitation of the identification procedure is the assumption of pure uni-and bi-axial stress and strain state in the center of the specimen, whereas full-field measurements could be employed in the future.

However, as the plunger validation yielded satisfactory results, this limitation proved to have a limited impact.

Our model relies on exhaustive experimental data and may be particularly sensitive to the range of boundary conditions applied. In the current study, the range of stress ratios was limited to the two extremes (40:0 and 0:40 N ratios for uniaxial, 40:40 N for biaxial tests). When the material experiences stresses out of this range in FE simulations, the constitutive model extrapolates the response, which may deviate from the actual response. Accordingly, strain concentrations have to be considered with caution. However, the validation of the FE model through very good predictions of reaction forces at the fixation points in the plunger test shows that the domain of validity of our model remains wide enough to consider various applications. Note that the plunger force at which fixation forces have been extracted (60 N) is equivalent to a 47 mmHg intra-abdominal pressure, which is slightly above the range of pressures that have been reported for a Valsalva manoeuvre [START_REF] Cobb | Normal intraabdominal pressure in healthy adults[END_REF]. These moderate conditions are consistent with what is expected for patients in hospital the first days or weeks after an abdominal wall repair surgery.

Note that fixation reaction forces were previously estimated using computational simulations for Dualmesh Gore ® and DynaMesh ® implants [START_REF] Lubowiecka | Behaviour of orthotropic surgical implant in hernia repair due to the material orientation and abdomen surface deformation[END_REF]. They were of the same order of magnitude compared to the present study. Present data are also comparable to data from porcine tissues despite strong dependence on mesh implant and fixation device [START_REF] Rieder | Mesh fixation in laparoscopic incisional hernia repair: Glue fixation provides attachment strength similar to absorbable tacks but differs substantially in different meshes[END_REF].

A possible limitation of our membrane model is that the bending stiffness of the textile was neglected. Future work should characterize it using the Kawabata evaluation system [START_REF] Kawabata | A guide line for manufacturing "ideal fabrics[END_REF] and consider it in the FE model as it may affect the prediction of creases, as shown in Figure 7. These creases may also occur during VIH surgical interventions and affect the biological integration process after surgery.

The main goal of our FE model is to be used in a simulation tool for assisting Stiffness change of an implant as a consequence of cyclic loading experienced every day by a patient was quantitatively described with an attempt to simulate a preconditioned state of the mesh [START_REF] Tomaszewska | Mechanics of mesh implanted into abdominal wall under repetitive load. Experimental and numerical study[END_REF]. The authors observed stiffening induced by uniaxial cyclic load in DynaMesh ® -IPOM implant. They remarked that the stiffening effect due to reload could possibly lead to higher fixation forces, fixation failure and hernia recurrence. In our work, stress softening was observed for a different textile implant under uniaxial cyclic load. While this might possibly lead to lower fixation forces, the permanent set data suggests possible functional issues related to the bulging. Once the implant is loaded to a certain level, e.g. during coughing or exercise, it would remain permanently deformed creating a bulge in the abdominal wall. Further loading would even increase the bulge, with dramatic consequences.

Stiffening and permanent set effects after cyclic loading were also observed on mesh implants in air and at 37 • C in physiological saline in four commonly used commercial hernia meshes [START_REF] Velayudhan | Evaluation of dynamic creep properties of surgical mesh prostheses -uniaxial fatigue[END_REF]. Unlike frequent bench-top testing conditions (air, room temperature), body temperature and environment increase polymer chain mobility and facilitate the relative movement of the yarns. This leads to higher deformations for the same load and, by consequence, to a decrease in the overall stiffness. This behavior is modified as the abdominal wall tissues integrate biologically the implant.

As a macroscopic continuum, the proposed model does not aim at simulating and understanding yarn-scale phenomena. The mechanical behavior of the yarns (stretching and bending response), the frictional effects, as well as the knitting pattern and manufacturing process all contribute to the global behavior of the implant. This highlights the need for comprehensive yarn-scale models, which would be extremely valuable to understanding the relationship between yarn mechanics, mesh structure and global behavior. However, due to their high computational cost, they have been limited to simple patterns [START_REF] Duhovic | Simulating the deformation mechanisms of knitted fabric composites[END_REF][START_REF] Liu | On the role of material architecture in the mechanical behavior of knitted textiles[END_REF][START_REF] Mckee | Finite element simulation of ballistic impact on single jersey knit fabric[END_REF][START_REF] Dinh | Prediction of mechanical properties of knitted fabrics under tensile and shear loading: Mesoscale analysis using representative unit cells and its validation[END_REF].

Linking yarn-scale approaches to continuous models would be a major contribution to the development of meshes and simulation of their in-situ mechanical performance and possible failure mechanisms.

As the implant is loaded, the pores tighten until the yarns become locked.

"Structural plasticity" may occur in this phase when the implant recovers its form only partially after unloading. Once the yarns are completely blocked, the material of the yarns takes over the structure and starts to deform under In addition to permanent set accumulation, mechanical fatigue of knitted textile implants for hernia repair is often neglected. It has been addressed

in bench-top tests applying several hundreds or thousands of loading cycles [START_REF] Rynkevic | In vitro study of the mechanical performance of hernia mesh under cyclic loading[END_REF][START_REF] Sahoo | Clinically relevant mechanical testing of hernia graft constructs[END_REF][START_REF] Patterson | The effect of cyclic loading on the mechanical performance of surgical mesh[END_REF]. However, hernia mesh fatigue testing raises many experimental challenges. Millions of loading cycles would be needed to simulate real life conditions. As polymers used to knit hernia mesh exhibit frequency-dependent mechanical behavior, the mechanical fatigue test should not be accelerated; it should be performed at some physiological frequency, e.g. that of the breathing or coughing. Such long-lasting characterization would ignore the process of tissue interaction and integration and provide information limited to the standalone implant. The material parameters of our model could possibly be modified to evolve with mechanical "ageing". However, the capacity of such constitutive model to take into account the accumulation of plastic deformation possibly leading to fatigue failure would strongly depend on the availability of relevant experimental data.

Although the constitutive behavior presented in this study has been developed for a specific mesh, it is representative of the general mechanical behavior of all meshes. By adjusting the model parameters and updating the evolution laws, it can reproduce mechanical features shared by all knitted meshes: nonlinearity, anisotropy, strong coupling and permanent [START_REF] Deeken | Biaxial analysis of synthetic scaffolds for hernia repair demonstrates variability in mechanical anisotropy, non-linearity and hysteresis[END_REF]. Therefore, adjusting the model to other meshes should not raise major difficulties (just increasing polynomial order for instance) and should not require to develop additional mechanical tests.

Conclusion

Knitted textile implants are critical in many surgical procedures, as for instance hernia repair. For the first time ever, we have reported an extensive mechanical characterization of one of these biomaterials, including cyclic uniaxial tension, planar equibi-axial tension and plunger type testing. These tests showed the complex mechanical behavior with strong non linearity, anisotropy and permanent set. In addition, an anisotropic hypo-elasto-plastic constitutive modeling has been developed based on these experiments. Satisfactory identification and promising preliminary validation have been found.

In the future, combined characterization at fiber (micro), pore (meso) and implant (macro) scales should be performed to better understand the knitted textile mechanical behavior and strength. In particular, such experiments should enable to study separately the respective roles of the material and its structure. They should also lead to more advanced modeling of rupture.

Besides, implementation of this characterization and modeling approach to additional surgical knitted textiles, coupled with abdominal wall models, should enable to simulate abdominal wall repair biomechanics. Immediate post-op performance of repairs, varying surgical techniques and knitted textile implants, could then be assessed in terms of force at primary fixations (failure mode criterion) and residual bulging (patient satisfaction). The final objective would be to provide recommendations to clinicians.

yielding constitutive parameters a, b and c.Figure 2 :

 2 Figure 2: Permanent set evolution for uniaxial tensile tests as a function of maximum applied strain. The permanent set was derived as the residual strain obtained after unloading. The

  model featured 12553 M3D3 and M3D4 membrane elements and 340 B31 beam elements. The plunger and circular patch surfaces were modelled as frictionless analytical contact surfaces. Symmetry boundary conditions were imposed to the edges of the mesh, while the lower end of the rods were fixed. A vertical displacement was imposed to the plunger. The simulation was performed with the Abaqus ® /Explicit solver. The numerical reaction forces at the fixation points were compared to three repeated experiments of a single cycle up to 310 N plunger force.

  (a) Experimental setup (b) Corresponding FE model

Figure 3 :

 3 Figure 3: Experimental setup of the plunger experiment allowing to measure fixation forces (a) and corresponding quarter FE model (b).

Figure 4 :

 4 Figure 4: Experimental results obtained from cyclic tensile tests: tension vs nominal strain for uniaxial loading in warp (a), weft (b) and 45°diagonal (c) direction, and equibiaxial loading (d).

Figure 5 :

 5 Figure 5: Polynomial fits used to build the constitutive model: (a) polynomial functions used to derive 6 stiffness constants E 0 1 , E 1 1 , E 0 2 , E 1 2 , G 0 12 and G 1 12 , in warp, weft and shear directions respectively; (b) polynomial function relating the part of elastic strain as a function of total strain: ε el = f (ε).

Figure 6 :

 6 Figure 6: Identification of the constitutive model and comparison to experimental results in uniaxial tension (warp, weft and diagonal direction), and equibiaxial tension. The individual R 2 (coefficient of determination) are presented for each loading direction. The global R 2 of the model is 0.94. .

  (a) At maximum (231 N) force (b) After unloading

Figure 7 :

 7 Figure 7: Logarithmic membrane principal strain in the mesh at maximum loading (a) and after unloading (b), showing the residual strain and remaining bulge. Strain magnitudes are within the experimental domain of uni and biaxial tests, except in the close vicinity of fixation points.

Figure 8 :Figure 9 :

 89 Figure 8: Experimental and numerical fixation forces at 60 N force: (a) force vectors at each fixation point and (b) force magnitude profile as a function of angle (0 and 180°corresponds to the warp direction and 90°to the weft direction.

  surgical interventions. A number of previous studies have highlighted the clinical significance of modeling the mechanical behavior of knitted textiles. Despite the different assumptions made in our model, it can still have clinical applications as for instance allowing the evaluation of permanent bulging and fixation forces post-surgery, before tissue integration and evolution of mesh properties. The model has been validated within a situation representing open defects, and is thus suited to evaluate techniques such as simple IPOM (Open Intraperitoneal Onlay Mesh) repair without primary closure of the defect.

  further load. Due to high complexity of the knitted patterns, both deformation mechanisms (structural and material) can coexist. Only some of the yarns can lock while the others can still move freely. Repetitive mechanical loading of the yarns can lead to alignment of polymer chains and increase the crystallinity, resulting in stiffening and possible local accumulation of damage and plastic deformations. Our constitutive model pools together material and structural plasticity and effectively reproduces the observed permanent set. This unique asset makes the model suitable to assess effects of repeated loading such as implant bulging after surgical procedures and before the beginning of tissue integration process.

Table 1 .

 1

							0.1
		5							Experimental points (warp)
			Experimental points (warp)					Experimental points (weft)
	Cauchy stress [MPa]	1 2 3 4	Polynomial fit (weft) Experimental points (shear) Polynomial fit (shear) Experimental points (weft) Polynomial fit (warp)		Elastic log. strain [-]	0.02 0.04 0.06 0.08		Polynomial fit
		0.02 0	0.04	0.06	0.08	0.1	0	0	0.05	0.1	0.15	0.2	0.25
				Elastic log. strain [-]						Total log. strain [-]
			(a) Elastic behaviour				(b) Yield evolution function

Acknowledgements

The work was supported from European Regional Development Fund-Project "Application of Modern Technologies in Medicine and Industry" (No. CZ.02.1.01 /0.0/0.0/17 048/0007280).


Valérie Declerck

Evelina Colacino

Xavier Bantreil

Jean Martinez

Frédéric Lamaty
email: frederic.lamaty@univ-montp2.fr

Vale´rie Declerck

Fre´de´ric Lamaty

Poly(ethylene

come

Grinding, dry mixing or ball-milling have proved their efficiency as unconventional techniques in the field of organic chemistry in solid state. 1 These techniques have been applied to a large number of organic transformations, but with limited use in the case of C-C bond formation, to aldol and Knoevenagel condensations, Michael additions, Baylis-Hillman and Wittig reactions, and asymmetric alkylation of Schiff bases. 2 Transition metal-catalyzed coupling reactions, especially with palladium, are very efficient and well developed C-C bond formation reactions. 3 Such reactions performed in a ball-mill have been limited to a small number of examples of transformations, mainly the Suzuki-Miyaura 4 and Mizoroki-Heck 5 cross couplings.

In the last few years we have developed a metal-based catalytic system using solid PEGs [poly(ethylene glycols)] as solvent and stabilizing/precipitating agents to separate the organic product from the metallic catalytic system. 6 We questioned whether the use of solid PEG would be adapted to the development of a useful catalytic system to perform Pd-catalyzed reactions in a ball-mill.

As a model reaction, the Mizoroki-Heck arylation 7 of tertbutyl acrylate to produce cinnamate was chosen. A mixture of PhI, tert-butyl acrylate, PEG, Pd(OAc) 2 and an inorganic base was reacted in a high energy stainless steel vibratory ball mill for 1 h at 30 Hz together with additives in some cases. Under traditional heating and microwave irradiation, PEG has already been shown to reduce Pd(OAc) 2 to Pd(0). 6d,8 However, under ball-milling conditions, it was found necessary to add sodium formate as a reducing agent for the activation of Pd(OAc) 2 (67% yield). This was the best from a practical point of view, as sodium formate is easy to handle, as well as from the results obtained with reducing agents such as H 2 (28%) or NaBH 4 (59%). At the end of the reaction, the mixture was dissolved in a small amount of CH 2 Cl 2 and precipitated in ether, filtered, evaporated and analyzed by 1 H NMR using CH 2 Br 2 as an internal standard. Integration of the signal corresponding to the H of the ester group of tertbutylcinnamate measured against the unique signal of the standard (CH 2 Br 2 ) provided the yield of the product. Results are presented in Table 1.

It was shown previously that a solid additive like sodium chloride could be valuable for the Mizoroki-Heck arylation reaction. 5 In our study, when PEG-3400-OH was used, adding NaCl was found to be detrimental to the reaction and barely any conversion was observed whilst 67% yield was obtained in the absence of NaCl (entries 1-2). Attempt to change the inorganic base by switching from K 2 CO 3 to Na 2 CO 3 (entry 3) resulted in a lower conversion most probably because the PEG-K + interaction is stronger than the corresponding PEG-Na + thus enhancing the base activity. When triethylamine was used, substrates also remained unchanged. We realized that during the course of the reaction, it appeared that PEG melted due to the slight heating of the ball-mill. 9 As the viscosity of the PEG used might have an influence on the reaction, we turned our attention to shorter polymers having a lower melting point. As a result, the yield increased by 20% by using PEG-2000-OH (entry 4). While reducing catalyst loading to 1 mol% resulted in no conversion (entry 5), increasing the quantity of acrylate from 1.2 to 5 equivalents resulted in total conversion (entry 6). The exact excess amount of acrylate needed for completion of the reaction was not studied in detail. Several PEGs were then tested under these conditions. The use of mono-and di-methylated PEG-2000 was detrimental as the yields dropped down (entries 7 and 8). Finally, shorter PEG-1100-OH gave a 74% yield (entry 9). It is worth noting that due to the stabilizing effect of oxygen atoms in PEG, no phosphine ligand was necessary for the catalysis. In the absence of PEG, no reaction occurred (entry 10). In addition, Mizoroki-Heck coupling generally requires elevated temperatures, unfriendly solvent and an inert atmosphere. In our study, the use of PEG in a vibratory ball-mill allowed the reaction to proceed in air under mild conditions.

The scope of the reaction was first explored by varying the olefin part. Unfortunately, the reaction was efficient only with tert-butyl acrylate. Either no reaction was observed or a very poor yield was obtained with acrylonitrile (0%), acrylamide (0%), methyl methacrylate (5%) or styrene (13%). The aromatic partner was then evaluated (Table 2). Different substituents could be used on the aryl moiety. [START_REF]Typical procedure for the ball-mill Mizoroki-Heck reaction: in a 10 mL stainless steel ball-mill reactor were added Pd(OAc) 2 (1.1 mg, 0.005 mmol), K 2 CO 3 (41.4 mg, 0.3 mmol), HCO 2 Na (1.4 mg, 0.02 mmol) and PEG-2000-OH (110 mg). Phenyl iodide (11.2 mL, 0.1 mmol), tert-butyl acrylate (73.2 mL, 0.5 mmol) and two stainless steel balls with a 7 mm diameter were then added. The mixture was then subjected to high speed vibration milling for 1 h at 30 Hz on a Retsch MM200. A minimum of CH 2 Cl 2 was then added and the mixture was precipitated in Et 2 O (150 mL) for 4 h at À20 1C. The solution was filtered and the filtrate concentrated in vacuo. After addition of 20 mL of CH 2 Br 2 as an internal standard, 1 H NMR analysis allowed measurement of a quantitative yield[END_REF] Electron-donating as well as withdrawing groups were well tolerated in the para position as p-MeO, p-CN, p-Br gave satisfactory yields (entries 2-4) and p-CHO was converted quantitatively (entry 5). Surprisingly, p-NO 2 substitution resulted in poor yield (entry 6). This is may be due to a physico-chemical characteristic of the starting 4-iodo nitrobenzene rather than an electronic effect of the substituent. A fluorine atom in the meta position allowed quantitative conversion (entry 7). However, having a methyl group in the ortho position reduced the yield to 29% (entry 8). Finally, the influence of the halide leaving group (I vs. Br, Cl) was explored. Only the iodine atom was reactive enough to provide the expected product, even when sodium iodide was added to the reaction mixture (entries 9-12).

It has already been demonstrated that palladium salts in the presence of polyethylene glycol are transformed into nanoparticles by classical activation. Because the order of the reduction potentials of the polyol and noble metals is not favorable at rt, this method requires high temperatures to reduce the noble metals. [START_REF] Bonet | [END_REF] This was confirmed by our previous work, using convection 6g or microwave heating. 6d Furthermore, when a long chain PEG (PEG-2000-OH or PEG-4000-OH) in the presence of Pd(OAc) 2 but no substrate was stirred at a temperature of 80-120 1C, nanoparticles of 5 nm were observed. 8 When using PEG-400-OH, phenanthroline as a reducing agent was necessary to observe the formation of 2-6 nm particles. 12 However, under vibrating ball-milling conditions, since the reaction temperature was mild, 9 an external reducing agent such as sodium formate was required for activation of Pd(OAc) 2 .

According to these results, we questioned whether formation of nanoparticles would also occur in the ball mill during the course of the reaction, knowing that physical and chemical properties would depend on the size and shape of the particles. Transition Electron Microscopy (TEM) analysis of the precipitates obtained after reactions reported in Table 1 with various PEGs was performed (Fig. 1, left column). In order to compare, experiments were run under the same reaction conditions but in the absence of phenyl iodide and tert-butyl acrylate (Fig. 1, right column). In all cases, light yellow Pd(OAc) 2 was transformed into a deep brown solid, indicating the formation of Pd/PEG nanoparticles. TEM analysis confirmed this experimental assumption (Fig. 1). To our knowledge, this represents the first example of generation of Pd-polymer nanoparticles using a ball-mill. 13 In the absence of substrates, the reaction led to significant particle aggregation, regardless of the PEG used. When substrates were present, the aggregation was less important and nice round-shaped nanoparticles could be obtained. The size of nanoparticles is directly connected to the size of the polymer. Indeed, the average size of the nanoparticles observed was 6-8 nm, 8-11 nm and 7-13 nm with PEG-1100-OH, PEG-2000-OH and PEG-3400-OH, respectively (Fig. 1a, c ande). In addition, with PEG-3400-OH, particles were found to be slightly more elliptical than with other polymers. Surprisingly, when PEG-2000-OH was mono-or di-methylated, the particle size changed to 7-13 nm and 5-7 nm, respectively (Fig. 1g andi).

Formation of palladium nanoparticles in PEG after a Mizoroki-Heck reaction activated by microwave irradiation and TEM analysis has already been reported. In PEG-400-OH, which is a liquid polymer, particles of 5-8 nm were obtained. 14 In PEG-3400-OH, following benzazepine synthesis, TEM analysis of the catalytic system revealed aggregated particles of 5-7 nm. 6d The particle aspect was quite different from those shown in Fig. 1 (left column). Thus, the activation This journal is c The Royal Society of Chemistry 2012 method seems to have a strong influence on the size and aggregation state of the nanoparticles. We reported herein a solvent-free/phosphine-free palladiumcatalyzed Mizoroki-Heck procedure in a ball mill. Under mild conditions, quantitative yields were obtained using the appropriate polymer PEG-2000-OH, with a correct tolerance toward functional groups on the aryl moiety. Nanoparticles formed during these reactions were characterized by TEM and were found to be size-dependent on the PEG used. In addition, activation using ball-milling yielded nanoparticles different from those obtained previously using convection or microwave heating.

We are grateful to the ''Service Commun de Microscopie Electronique et Analytique'' of the University Montpellier 2 for TEM analysis and fruitful discussions.

Notes and references

Table 1

 1 Optimization of ball-milling Mizoroki-Heck reaction

	1	PEG-3400-OH	1.2	NaCl	0.4
	2	PEG-3400-OH	1.2	-	67
	3 c	PEG-3400-OH	1.2	-	2
	4	PEG-2000-OH	1.2	-	88
	5 d	PEG-2000-OH	1.2	-	0
	6	PEG-2000-OH	5	-	100
	7	MeO-PEG-2000-OMe 5	-	46
	8	MeO-PEG-2000-OH 5	-	27
	9	PEG-1100-OH	5	-	74
	10	-	1.2	-	0

a Entry PEG Acrylate (equiv.) Additive Yield b (%)

a Reaction conditions: phenyl iodide (0.1 mmol), tert-butyl acrylate, Pd(OAc) 2 (5 mol%), HCO 2 Na (0.2 equiv.), K 2 CO 3 (3 equiv.), PEG (110 mg). b Yields were determined by 1 H NMR using CH 2 Br 2 as an internal standard. c Na 2 CO 3 was used instead of K 2 CO 3 . d 1 mol% Pd(OAc) 2 .

Table 2

 2 Exemplification of the ball-milling Mizoroki-Heck reaction a

	Entry	ArX	Additive	Yield b (%)
	1	PhI	-	100
	2	4-MeO-C 6 H 4 -I	-	73
	3	4-NC-C 6 H 4 -I	-	62
	4	4-Br-C 6 H 4 -I	-	68
	5	4-OHC-C 6 H 4 -I	-	100
	6	4-O 2 N-C 6 H 4 -I	-	35
	7	3-F-C 6 H 4 -I	-	100
	8	2-Me-C 6 H 4 -I	-	29
	9	PhBr	-	0
	10	PhBr	NaI	0
	11	PhCl	-	0
	12	PhCl	NaI	0

a Reaction conditions: phenyl iodide (0.1 mmol), tert-butyl acrylate (0.5 mol), Pd(OAc) 2 (5 mol%), HCO 2 Na (0.2 equiv), K 2 CO 3 (3 equiv.), PEG-2000-OH (110 mg). b Yields were determined by 1 H NMR using CH 2 Br 2 as an internal standard.

11780 Chem.

Commun., 2012, 48, 11778-11780