Pauline Assemat

Gwenael Hannema

Josie Carberry

Danielle Mitchell

Karen Andrews

Anthony Dart

Jaye Chin-Dusting

Kerry Hourigan

D Michell

Numerical and in vitro experimental study of arterial deformation and buckling under hypertension and atherosclerotic conditions

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

Introduction

Atherosclerosis occurs due to the build-up and infiltration of lipid streaks in artery walls, leading to plaques. It is an inflammatory disease characterized by a succession of biochemical and biomechanical complex processes [START_REF] Chatzizisis | Role of Endothelial Shear Stress in the Natural History of Coronary Atherosclerosis and Vascular Remodelling: Molecular, Cellular and Vascular Behaviour[END_REF][START_REF] Pello | A Glimpse of the Phenomenon of Macrophage Polarisation During Atherosclerosis[END_REF]. In plaques, the lipid core is separated from the circulating blood by a fibrous cap of variable thickness. As plaques develop, they induce a remodelling of the arteries [START_REF] Stone | Regions of Low Endothelial Shear Stress are the Sites where Coronary Plaque Progresses and Vascular Remodelling Occurs in Humans: an in Vivo Serial Study[END_REF]. Plaques can be divided into two distinct types: those that rupture (vulnerable) and those that are less likely to rupture (stable), a characteristic that is associated in some papers to the respective aneurysm-like and stenotic shapes [START_REF] Feldman | Vulnerable Anatomy; The Role of Coronary Anatomy and Endothelial Shear Stress in the Progression and Vulnerability of Coronary Artery Lesions: Is Anatomy Destiny?[END_REF]. While it is accepted that plaque vulnerability is influenced by fibrous cap thickness and the size of the lipid core [START_REF] Avgerinos | Biomarkers for Diagnosis of the Vulnerable Atherosclerotic Plaque[END_REF], it also depends on artery remodelling, the loads applied on the solid and the forces induced by the fluid for which no general agreement has been found.

The complexity of the present problem comes from the coupling of biological, chemical and mechanical effects. In the following, the literature review will be restricted to the influence of mechanical effects on plaque rupture. Most of the numerical studies in real arteries performed until now focus on the calculation of the flow around plaques, considering the arteries as rigid walls [START_REF] Chaichana | Computational Fluid Dynamics Analysis of the Effect of Plaques on the Left Coronary Artery[END_REF][START_REF] Groen | Plaque Rupture in the Carotid Artery is Localized at the High Shear Stress Region[END_REF]. These papers indicate that the high wall shear stress (WSS) found on plaques could be responsible for plaque rupture; however, Sadat et al. [START_REF] Sadat | Biomechanical Structural Stresses of Atherosclerotic plaques[END_REF] showed recently that structural stresses are thought to play a more important role, being several orders of magnitude higher than the WSS induced by blood flow. To calculate the structural stresses, a complex FSI approach taking into account large deformations is required. A three-dimensional (3D) fluid structure interaction model has re-vealed an association between regions of high plaque critical stress and fibrous cap disruption in human carotid [START_REF] Tang | Sites of Rupture in Human Atherosclerotic Carotid Plaques are Associated with High Structural Stresses: An In Vivo MRI-Based 3D Fluid-Structure Interaction Study[END_REF]. Furthermore, maximum principal stresses were identified on shoulder regions of the fibrous caps of carotid plaques according to 2D FSI simulations [START_REF] Kock | Mechanical Stresses in Carotid Plaques Using MRIbased Fluid-Structure Interaction Models[END_REF]. An histology based finite element analysis was used to investigate the peak circumferential stresses in aortic and brachiocephalic plaques of ApoE -/-mice [START_REF] Vengrenyuk | Computational Stress Analysis of Atherosclerotic Plaques in ApoE Knockout Mice[END_REF]. Their results showed the aortic plaque stresses only slightly depend on the cap thickness whereas for the brachiocephalic artery lesions, it is shown that there is an exponential growth of peak cap stress with decreasing cap thickness. Altogether, these recent papers underline the multiplicity of the origins of plaque rupture and the complexity of finding a general agreement of the main parameters that engender this event. In addition, very few computational studies have considered the elastic properties of arteries or have investigated the hemodynamic parameters and mechanical forces induced in realistic vessel geometry by a blood flow past plaque. Therefore, developing simplified FSI models under those requirements is a starting point towards further understanding of atherosclerotic development.

The complexity of the task comes from the large number of influential parameters involved in the mechanical model (geometrical parameters, transmural pressure (external internal pressure), oscillating flow rate, variable elasticity and density properties of the tissues, non-Newtonian effects). In the present work, the geometries considered are restricted to straight vessels and straight vessels with stenosis (stable shape). In addition to the physical parameters, a large number of numerical parameters must be checked in order to insure the convergence of the numerical methods as large deformations are involved. Indeed, in the physiological range of parameters, it has been shown [START_REF] Han | A biomechanical model of artery buckling[END_REF] that straight vessels buckle when a critical transmural pressure is reached. In addition, whilst the general dynamics of fluidconveying systems have been already studied for "rigid" materials [START_REF] Paidoussis | Fluid-Structure Interaction[END_REF], in the present study, materials with weak Young's modulus are considered, as well as the coupling of flow and high transmural pressure, effects which together have implications on the mechanical stability of the system. Hence, within physiological range, vessels are readily prone to large deformation.

In the following, after the presentation of the model, the numerical approach will be described and compared to linear theory. In addition, a description of the method used to extract the elastic properties from the experiments will be presented. Some key results will be presented at the end of the article.

Model

Geometry

The model used in the present study follows the work of Young and Tsai [START_REF] Young | Effect of a time dependent stenosis on flow through a tube[END_REF] who developed an axisymmetric model to represent the stenosis:

r(z) = D 2 - δ s 2 1 + cos πz Z s , (1)
where r(z) represents the radius of the vessel, D s the stenosis diameter at the maximum constriction location (z = 0), Z s half of the stenosis length and δ s the stenosis height

δ s = (D -D s )/2.
This geometry is represented figure 1. Most of the paper will concern straight vessel for which δ s = 0.

Model equations

The fluid and solid behaviour are described by the following coupled equations:

ρ f ∂U ∂t + (U • ∇)U = ∇ • σ f , (2)
∇ • U = 0, (3) ρ s (D tt η -f s ) = ∇ • σs , (4)
where ρ f is the fluid density, U the 3D velocity field and σ f the fluid stress tensor:

σ f = -P Ī + µ ∇U + ∇U T , (5)
with Ī the unit tensor, P the pressure and µ the dynamic viscosity of the fluid. In the solid equation, ρ s represents the material density, D t the material derivative, η the 3D structure displacement, f s the body force and σs the structural stress tensor. Concerning the boundary conditions: a Poiseuille velocity profile is imposed at the inlet when flow is considered whereas ∂P/∂z = 0 at the outlet. At the fluid solid interface, the boundary conditions are:

σ f • n f = -σs • n s , (6) U = ∂ t η. (7)
In addition, a contant atmospheric pressure is applied along the external surface of the solid. The remaining term to be determined is σs which depends on the material properties of the solid and on the strain (ε). The calculation of this term will be detailed in the next section. For the numerical calculation, the dynamic viscosity, fluid and solid densities remain constant with µ = 4.10 -3 Pa.s [START_REF] Trachet | The Impact of Simplified Boundary Conditions and Aortic Arch Inclusion on CFD Simulations in the Mouse Aorta: A Comparison with Mouse-Specific Reference Data[END_REF], ρ f = ρ s = 1060kg/m 3 [START_REF] Feintuch | Hemodynamics in the Mouse Aortic Arch as Assessed by MRI, Ultrasound, and Numerical Modeling[END_REF]. As a consequence, the fluid is considered Newtonian, a reasonable assumption provided the shear rate is greater than 150s -1 [START_REF] Wells | Shear rate dependence of the viscosity of whole blood and plasma[END_REF], which is the case in the range of parameters considered for the present study. The density of the solid has been considered equal to the fluid in order to prevent the influence of the gravity effects. This approximation is acceptable considering that in a recent study [START_REF] Shinohara | Atherosclerotic plaque imaging using phasecontrast X-ray computed tomography[END_REF], Shinohara et al. showed that the density parameter range for the plaque and the artery walls is 1010 -1080kg/m 3 .

Stuctural stress tensor and theoretical divergence model

Structural stress tensor

The arterial walls are hyperelastic material for which the elasticity properties, the Young's modulus, depends non-linearly on the load applied on the vessel; however, in the present study, only orthotropic (linear) material will be considered in order to validate our FSI model. The method described in the following indicates how to extract the constant axial and circumferencial Young's modulus from experiments.

Following the work of Liu and Han [START_REF] Liu | Mechanical buckling of artery under pulsatile pressure[END_REF], free inflation test has been performed. The method consists of imposing an internal pressure, the vessel being plugged at one extremity (figure 2). The expansion of the vessel in the axial and longitudinal directions is recorded by a camera on the top. The pressure is increased to measure its influence on the deformation of the vessel (stretch). From the obtained stretch versus pressure curves, the structural stresses σ s can be calculated using Lame's equations.

As the Young's modulus E is a function of the strain ε, the relation between stress and strain is σ s = E(ε)ε. Figure 3 indicates the stress as a function of the strain for a rat carotid. A constant Young's modulus is then extracted by calculating the slope of the lines indicated on figure 3. Thoses lines are obtained by imposing the following equality using energy arguments.

ε p 0 σ s dε = Ē(ε p ) 1 2 ε 2 ε p 0 , ( 8
)
where ε p is chosen accordingto a pressure of interest. Further information on the procedure can be found in [START_REF] Hannema | Three-dimensional fluid structure interaction methods applied to blood flow in straight vessels with atherosclerotic plaques[END_REF]. The other parameter, the Poisson ratio, necessary to describe the elastic properties of the tissue, is fixed to ν = 0.49 (nearly incompressible material).

Theoretical buckling model

To validate the numerical results, a comparison with the theory has been performed. In [START_REF] Paidoussis | Fluid-Structure Interaction[END_REF], the author indicates that the dynamics of buckling of cylindrical structures can be modeled in first approximation by the following one dimensional equation:

EI ∂ 4 w ∂x 4 + MU 2 -T + pA(1 -2νδ) ∂ 2 w ∂x 2 +2MU ∂ 2 w ∂x∂t + (M + m) ∂ 2 w ∂t 2 = 0, ( 9
)
where E corresponds to the Young's modulus in the axial direction (E z ), I to the area moment of inertia, w to the lateral deflection, T to the externally applied tensions, p to the mean pressure, ν to the Poisson ratio and δ depends on boundary condition (in our case, closed pressure vessel δ = 0). In addition, M is the mass of the fluid flow per unit length and m the mass of the pipe per unit length.

As divergence is a static rather than a dynamic form of instability, one needs to consider only the time-independent terms of equation [START_REF] Paidoussis | Fluid-Structure Interaction[END_REF] to analyse the dynamics of the system. Hence, the critical buckling pressure can be obtained following a procedure similar to [START_REF] Timoshenko | Theory of Elastic Stability[END_REF] and is given by:

pcritical = 4 πD 2 4π 2 EI L 2 k + T -MU 2 , ( 10
)
where L k represents the buckling length, which depends on the length of the vessel, and on the boundary conditions (in the present case L k = L/2).

Numerical approach and validation

The software codes used to solve the equation are Ansys CFX for the fluid and Ansys Mechanical APDL for the structure. ANSYS CFX is an element-based finite-volume method with second-order accurate discretization schemes in space and time. Ansys Mechanical APDL uses a finite element based structural solver. In this study, hexahedral quadratic elements are used for the spatial discretisation of both fluid and solid domains. The coupling between the fluid and the structural solver is achieved using Ansys Multi-field solver capabilities for which field variables (displacement and forces) are exchanged between both solvers through the fluid-structure interface in an iterative process. For the mesh, attention has been paid to define the characteristic size of the grid cells in the fluid near the interface. This size must be compatible with the capture of the boundary layer (the first computational cell is located within the viscous sublayer) but it must prevent the folding of the mesh due to large deformations(large cells). A compromise has been found to fulfill these two requirements (figure 4) and to optimize the radial mesh distribution. The method used to define the mesh in the fluid part is an O-grid method.

A numerical parameter sensitivity study, not presented here, has been performed on the vessel studied by [START_REF] Liu | Mechanical buckling of artery under pulsatile pressure[END_REF]. For L = 45mm, D = 3.38mm and h = 1mm, the optimal mesh distribution is 36 nodes in the circumferential direction and 100 in the axial one. The time step is 0.005s. The pressure is increased linearly in time to reach 100mmHg (physiological pressure). The Young's modulus calculated from [START_REF] Liu | Mechanical buckling of artery under pulsatile pressure[END_REF] and used for the present calculation are E z = 45.5KPa and E r = 308.6KPa. The code has been validated comparing the buckling critical pressure in straight vessels to the theory by setting U = T = 0 in equation 10. The numerical buckling critical pressure is defined as the transition between the axisymetric expansion and the sinusoidal shape observation. The mean variation between theoretical prediction and numerical approach is around 5% when the length and the thickness of the vessels are varied. The post buckling behaviour will be discussed in the next section.

Figure 4: Typical mesh slice representing the fluid domain (inner disk -light grey filling) and the solid domain (outer diskdark grey filling).The monitor points represented by crosses are used to track the system evolution.

Results

In this paragraph, few results will be presented and discussed.

Figure 5 shows a comparison between experimental and numerical results for an internal pressure of 95mmHg (post-buckling regime). A good qualitative agreement is obtained between the experimental and numerical approach, the deformation being of the same order as the internal pressure increased from 0 to 95 mmHg. Additional improvements are needed to obtain a quantitative matching, including taking into account the hyper-elastic properties of the vessel wall. A steady flow state with a flow rate of 160ml/min [START_REF] Liu | Mechanical buckling of artery under pulsatile pressure[END_REF] has been added to the previous configuration. Its destabilising effect induces a reduction in the critical buckling pressure. Again, the results have been compared to the linear theory and an acceptable mean error of 6% has been obtained. For the model vessel whose dimensions have been described in the previous section, the numerical critical pressure is 24.8mmHg compared to 25.303 in theory.

Finally, numerical simulations were conducted on stenotic and aneurysm-like vessels. Results are only presented here for a stenotic configuration. Figure 6 shows the Von Mises stresses for a vessel with 60% stenosis (percentage of fluid domain narrowing) for an internal pressure of 60mmHg. The maximal stress on the plaque is located along the fluid solid interface (the fibrous cap ), location at which the stenosis is the most stretched and compressed by the buckling phenomenon. A close examination of figure reveals that local maxima are observed on the inner side of the buckling at the shoulder of the stenosis (boundary between the stenosis and the healthy wall), a result compatible with the 2D simulation of [START_REF] Kock | Mechanical Stresses in Carotid Plaques Using MRIbased Fluid-Structure Interaction Models[END_REF]. In addition, it has been found that the maximum of the structural stress increases when the pressure increases (not presented here [START_REF] Hannema | Three-dimensional fluid structure interaction methods applied to blood flow in straight vessels with atherosclerotic plaques[END_REF]), a result in agree-ment with the fact that people with hypertension (higher internal pressure) suffer more cardivascular conditions.

Figure 6: Von Mises stresses in a buckled axisymmetric stenosis at 60mmHg.

Conclusion and discussion

An FSI model has been developed for simplified axisymmetric geometries. Characteristic vessel dimensions found in the literature [START_REF] Liu | Mechanical buckling of artery under pulsatile pressure[END_REF] as well as experimental samples of rat carotid aortas were used to calculate realistic Young's modulus and to validate an FSI approach. The FSI model was able to predict the buckling threshold within a range of physiological parameters under both quiescent and steady state flow. Furthermore idealised cardiovascular conditions were modelled by axisymmetric stenosis. The stability analysis permitted to highlight zones under traction or compression and the resulting stresses acting on the vessel wall. The results showed that the maximum structural stresses were located near the thin cap, a condition favorable for plaque rupture.

Figure 1 :

 1 Figure 1: Model geometry of an idealised axisymmetric stenotic vessel.

Figure 2 :Figure 3 :

 23 Figure 2: View of the experimental inflation test in rat carotid.

Figure 5 :

 5 Figure 5: Comparison between experimental (upper frame) and numerical (lower frame) buckled mode shape at a lumen pressure of 95mmHg.

Acknowledgements

The authors acknowledge financial support from the Australian Research Council, under grant no. DP110100434.