Haemodynamical stress in mouse aortic arch with atherosclerotic plaques: Preliminary study of plaque progression.

Atherosclerotic plaques develop at particular sites in the arterial tree, and this regional localisation depends largely on haemodynamic parameters (such as wall shear stress; WSS) as described in the literature. Plaque rupture can result in heart attack or stroke and hence understanding the development and vulnerability of atherosclerotic plaques is critically important. The purpose of this study is to characterise the haemodynamics of blood flow in the mouse aortic arch using numerical modelling. The geometries are digitalised from synchrotron imaging and realistic pulsatile blood flow is considered under rigid wall assumptions. Two cases are considered; arteries with and without plaque. Mice that are fed under fat diet present plaques in the aortic arch whose size is dependent on the number of weeks under the diet. The plaque distribution in the region is however relatively constant through the different samples. This result underlines the influence of the geometry and consequently of the wall shear stresses for plaque formation with plaques growing in region of relative low shear stresses. A discussion of the flow field in real geometry in the presence and absence of plaques is conducted. The presence of plaques was shown to alter the blood flow and hence WSS distribution, with regions of localised high WSS, mainly on the wall of the brachiocephalic artery where luminal narrowing is most pronounced. In addition, arch plaques are shown to induce recirculation in the blood flow, a phenomenon with potential influence on the progression of the plaques. The oscillatory shear index and the relative residence time have been calculated on the geometry with plaques to show the presence of this recirculation in the arch, an approach that may be useful for future studies on plaque progression.


Introduction
Atherosclerosis is a chronic and progressive inflammatory disease of the arterial wall involving various processes such as lipid deposition and oxidation, leukocyte infiltration, smooth muscle cell migration and extracellular matrix production [1][2][3]. The lipids form the core of the atherosclerotic plaque, while the smooth muscle cells and extracellular matrix produce an overlying fibrous cap. Rupture of the fibrous cap can induce thrombus formation on the plaque surface and can result in myocardial infarction or stroke. Plaque vulnerability to rupture depends on the fibrous cap thickness [4], in addition to luminal remodelling and blood haemodynamics. Depending on various factors, not yet all well identified, the plaques can remain stable for the whole life of the subject or can become vulnerable [5], leading eventually to a fatal acute syndrome. Therefore, understanding the plaque development processes and the vulnerable phenotype selection mechanisms is crucial.
Atherosclerosis is a complex and multifaceted disease. Aside from the complexity of the biological and biochemical processes involved in plaque development, it has been shown in the literature that plaque evolution is highly related to mechanical effects (see Assemat and Hourigan [6] for a review). Two main types of mechanical factors contributing to plaque development are reported in the literature: the haemodynamical stresses and the tissue (structural) stresses. Their relative contribution to plaque development is not well understood [7] and depends on the plaque age so that these contributions vary with the formation, progression or rupture stages.
Concerning the rupture, Maehara et al. [8] have shown the existence of two classes of rupture sites for vulnerable plaques: 63% of ruptures occur in the shoulder region, where the plaque joins the healthy intima, while the remainder (37%) occurs in the centre of the fibrous cap. Two types of vulnerable plaques have been identified in the literature: the rupture-prone and the erosion-prone [5]. For the first type, two theories have emerged from the mechanical point of view. Plaque rupture is due to a chronic fatigue structural damage (structural cyclic load) according to some [9,10], or is due to an acute structural fracture failure (structural stress greater than the ultimate tensile stress limit [11,12] or delamination [13]) according to others. Furthermore, besides the pure tissue stress aspect, the presence of microcalcifications [14][15][16], or intraplaque haemorrhages [17,18] are thought to induce additional stresses that play a role in the rupture process. For the erosion-prone type, little is known, but in general, eroded plaques are scarcely calcified, rarely associated with expansive remodelling, and only sparsely inflamed [5,19].
Concerning the formation aspect, it has been shown in the literature that plaques develop at particular sites in the arterial tree [20,21] and it is widely accepted that this localization is influenced predominantly by the geometry [22,23] and haemodynamic forces. Those forces are illustrated in the literature mainly by the calculation of haemodynamic parameters such as wall shear stress (WSS) [24,25], oscillatory shear index (OSI) [20,26] or relative residence time (RRT) [27,28]. Whereas some pioneer studies [21] indicated that low WSS for plaque formation was inferior to 0.4 Pa and high WSS preventing the plaque formation superior to 1.5 Pa in human arterial tree, this hypothesis has been intensively discussed and there is currently no general agreement on the optimal haemodynamical factor or set of factors that can be correlated with precision with the plaque distribution in the arterial tree [27,29]. Generally, plaques are found near bifurcations of large and medium size arteries and high curvature areas [20,30]. In particular, plaques are mainly found in the carotid bifurcation, in the coronary bed, in the femoral, abdominal and iliac arteries, and in the aorta [31]. In the present paper, the region of interest is the aortic arch, which presents both a high curvature zone and bifurcations.
As plaques develop, they may cause luminal narrowing, leading to a reduction in vessel volume, or undergo expansive remodelling to maintain lumen diameter [2,24,32,33]. This remodelling has been first underlined by Glagov et al. [34] for coronary arteries in which Chatzizisis et al. [2] suggest that low WSS causes intense inflammation, excessive wall and lumen expansion inducing outward plaque formation, and that arterial sections with slightly low WSS induce limited inflammation and lipid accumulation, leading to stenotic shape remodelling. In that study, the authors indicate that the vulnerable plaques are associated with a biomechanical self-feeding cycle causing an excessive expansive remodelling (outward remodelling) whereas stenotic remodelling induces stable plaques. Similar conclusions related to the link between vulnerability and remodelling are presented by Phinikaridou et al. [33] in a study conducted on white rabbit abdominal aortas. However, WSS amplitude does not seem to be the only criterion for plaque stability and mouse cuff models have been developed to study the oscillatory characteristics of the wall shear stress [35]. Cheng et al. [36], in a study involving an ApoE−/− cuff model (right common carotid artery), have suggested that vulnerable plaques are related to areas of lowered WSS (upstream of the cast), whereas stable ones are found in zones of high oscillatory shear stress (downstream of the cast). Their study indicates that the instantaneous dynamics also seems to play a role in the phenotype selection. This instantaneous dynamics can be calculated using numerical methods as described in the present paper.
Finally, the details of the processes of the growth of plaques after their initiation and, in particular, the influence of their presence on the blood flow, remain poorly understood. In their position paper, the European Society of Cardiology Working Group on Atherosclerosis and Vascular Biology [5] underlines that "the natural history of vulnerable plaques such as speed of development, lifetime (persistency) and fate is presently mostly unknown." Few papers give some insight into the question of plaque progression, and limitations are found in these papers. These limitations are mainly related to the number of time points used to define the plaque progression, to the number of samples, to the image resolution or to the excessive simplification of the numerical model or to the diversity of locations in the arterial tree. In the aortic arch region, Harloff et al. [37] have conducted an in vivo study at one time point and revealed a good correlation between cross section planes with critical WSS (defined as WSS below an individual threshold) and plaque positions but limited correlation between critical WSS and plaque angular positions. Hoi et al. [38] and Tomita et al. [39] have, for their part, commented on the importance of the geometry specificities for the location of plaques, in addition to haemodynamical factors. Although it is well known that haemodynamical factors are coupled with the geometrical features, in complex geometry this relation is non-trivial and investigations of the specific geometry could also give a clue to understanding plaque growing processes. The three previously cited articles do not consider the time evolution of plaques, and two of them underline the necessity of developing numerical methods that calculate blood flow in aortic arch with plaques.
To summarize this review introduction, the calculation of wall shear stress seems to have an influence on the formation, on the development and on the rupture of plaques so its knowledge is a first step towards the prediction of acute events. In addition, according to the literature, the phenotype of the plaques seems to be dependent on haemodynamical stresses, so that it is important to develop methods to calculate their distribution in realistic geometries in the presence of plaques. To the authors' knowledge, except for the related work published by our team [40], no numerical study has been yet conducted for this purpose in aortic arches with plaques, a restriction perhaps due to the challenge of generating high resolution 3D geometries with plaques (small size of plaques compare to the resolution). Another difficulty is to generate meshes that enable the numerical methods to converge in highly complex geometries (bifurcations, high curvature zones, plaques). The present paper will give some insight into these questions and address some of the challenges described in Tang et al. [41], who underline the lack of the data necessary to develop clinical efficient predictive models of plaque vulnerability.
The purpose of the present study is, in addition to develop a methodology, to investigate the haemodynamics of blood flow through the mouse aortic arch with plaques. The focus is to improve understanding of plaque development and vulnerability in a realistic geometry and provide insight into the influence of plaque on haemodynamics. In addition, plaques in aortic arches have been shown to increase stroke frequency in patients suffering from cardiovascular diseases [42,43]. In particular, it multiplies by six the risk of stroke after cardiopulmonary bypass operations [44], therefore understanding plaque development in this region is highly important. Furthermore, while the results presented here are based on the rigid wall assumption, this paper seeks to validate a numerical approach that is adaptable to the implementation of fluid structure interaction methods. The results show the plaque changes observed between mice fed after 10, 12, 14 and 16 weeks of fat diet, and the blood flow dynamics in a geometry with and without plaque obtained with numerical calculations. The emphasis of the present work is to develop tools and undertake a preliminary study that will allow the detailed investigation of plaque progression in future studies with a larger number of subjects.

Animal preparation
An apolipoprotein E knockout (ApoE −/−) mouse model sourced from Dr Grant Drummond colony (Monash University) was used as a model of spontaneous atherosclerosis development. Mice were fed with a high fat diet at 5 weeks of age (HFD SF00-219, 21% fat; 0.15% cholesterol) for 10, 12, 14 or 16 weeks. Mice were perfusion fixed at 100 mm Hg via the left ventricle of the heart. The first perfusion solution was heparinised saline (50 mL) to clear blood from the vessel and this was followed perfusion fixation with 50 mL of a low osmolarity modified Karnovsky's fluid, (2% formalin, 2.5% glutaraldehyde, 0.15 M cacodylate buffer). Tissues were post fixed for 24 h then embedded in low melting point agarose and shored in 1% formalin to maintain hydration until scanning.
Following dissection and dehydration through graded butanol, the tissues were embedded in paraffin (12 week HFD) or agarose (10, 14 and 16 weeks HFD). Preserving vessel morphology as near to native status is of paramount importance and the tissue collection protocol has been designed to maintain vessel geometry in the best possible way. Other protocols however can be found in the literature where MicroFil is used to fix the tissues [35]. Because the samples may suffer shrinkage artefact during the fixation procedure, numerical calculations have been conducted in a geometry scaled up by 30% (see Computational fluid dynamics (CFD) simulations section). The experiments were conducted following the Guidelines from the 2004 Australian Code of Practice for the Care and Use of Animals for Scientific Purposes and relevant Victorian legislation and were approved by the local Animal Ethics Commission.

Vessel geometry
In order to study the geometry of the plaques with accuracy, an ex vivo high resolution synchrotron imaging technique has been used. Indeed, as the characteristic thickness of plaques and of aorta walls in mice is of the order of 100 μm, high resolution images are necessary (see Assemat and Hourigan [6] for a discussion of different image technique resolutions). The mouse aortas were imaged by micro-computed tomography (μ-CT). This employs 15 keV up to 20 keV monochromatic synchrotron X-rays (PSI, Swiss Light Synchrotron, Switzerland) that are detected by a CCD detector. The μ-CT data were reconstructed into two-dimensional (2D) vessel slices using the programme X-TRACT developed at the CSIRO [45] (http://www.ts-imaging. net/Services/AppInfo/X-TRACT.aspx), which employs the Paganin algorithm [46]. Image reconstruction and segmentation of the vessel wall, plaque and lumen were conducted in the software AVIZO (Fig. 1a, b). Due to the high resolution of the images (1.6 to 5.92 μm depending on the samples), it is possible to differentiate not only the plaque from the wall component but also the fibrous cap from the lipid core if required (Fig. 1c). This differentiation enables us to study the disease state (with plaques) and the healthy state (in which plaques have been removed numerically [40]). This approach can be justified in the present study as Glagov remodelling [34] does not seem to occur in young ApoE −/− mice [47] as well as intimal thickening seems to be minor or absent in ApoE−/− mice [48]. The geometry was compared with histological data (obtained by staining 5 μm cross sections of the aorta with Masson's Trichrome [40]) to verify that plaques were well defined. A smoothing factor of 10 has been used in the software NX 7.5 Siemens on the geometry obtained in AVIZO software. The outlets were then cropped in ICEM to obtained straight boundary conditions and extensions were added to the end of the vessel following the normal direction of the cropped section.
Two approaches have been applied in this study. First, the characteristics of the geometry have been calculated using the software AVIZO for mice after 10, 12, 14 and 16 weeks under HFD. Then, the haemodynamical forces applied on the vessel wall and plaques have been calculated using the CFD software ANSYS CFX for one mouse (12 weeks under HFD). To this aim, the different materials (plaque, lumen) were imported into ICEM CFD in order to generate the mesh.

Mesh generation
The ends of the vessels are extended for the flow to be fully developed in the aortic arch zone. This approach is discussed in more detail in [40]. B-blocking and O-grid methods are used to generate the mesh with hexahedral elements. These methods allow a tight control of the mesh density and topology; in particular, the size and number of the boundary layers along the geometry can be controlled. An example of a hexahedral mesh is given in Fig. 2.

Computational fluid dynamics (CFD) simulations
In this study, boundary conditions at the inlets and outlets will be defined based on redistributed oscillating volumetric flow rates (see [40] for a review), as described by Trachet et al. [49], who noticed that, under physiological conditions, there is a time lag between the waves measured in the ascending and descending aortas due to the elasticity of the walls. When, rigid boundary conditions are considered, the mass balance needs to be corrected at all times. To achieve this, a b Trachet et al. [49] determined the transit time and shifted all the waveforms to make sure that they all concurred through the aortic arch. In addition, they showed that mice specific geometries had a little impact on the volumetric flow rate and they concluded that a percentage of flow rates equal to 22.5%, 11.2%, and 11.3% in the secondary branches can be suitable for computations in mice aortic arches. In the present work, a Fourier series was used to obtain the equations of the flow rates from the data in Trachet et al. [49] ( Fig. 2c) obtained also for male ApoE −/− mice. It is of importance to notice although atherosclerotic plaque formation was not the purpose of this recent study, the mice were quite old (34 ± 5 weeks) and plaques could have be expected in the arterial tree [50] even if the spatial resolution of the imaging techniques might not be able to show them. As a consequence, we think that these boundary conditions are suitable for an ApoE −/− mouse aortic arch with plaques. In addition, in [40], a comparison in a healthy configuration has been performed between two sets of boundary conditions and the spatial distribution of the wall shear stress has been shown to be preserved so that similar boundary conditions will be used in the healthy configuration. Finally, a heart rate of 400 beats per minute was assumed [49]. Walls were modelled as rigid with a no slip boundary condition.
Numerical modelling was performed in ANSYS CFX. The Navier-Stokes equations were solved using the element based finite volume methods. Blood was considered Newtonian, with ρ a density and μ a dynamic viscosity of 1060 kg/m 3 and 0.0035 Pa·s, respectively [49]. The calculation of the time evolution of the solution was conducted for three cycles (0.45 s), with a time step of 0.001 s. Optimisation of the mesh was performed to ensure the convergence of the solution and its independence from the numerical parameters [40]; calculations were performed on a larger number of cycles to ensure the time periodicity of the results.
As the present study considers a geometry with rigid walls, only the haemodynamical factors are of interest. These parameters are the wall shear stress vector τ w , with τ w = μ∂U/∂n| wall , where U is the velocity field, n the outwards normal to the wall, and μ the dynamic viscosity, and the time average wall shear stress (TAWSS = ∫ 0 T |τ w |dt/T).
The current geometry average inlet diameter is 1.3 mm, in agreement with the literature; however, to check the influence of potential uniform shrinkage, some numerical calculations have been conducted in geometries with different scaling (Fig. 3). For example, when the geometry is scaled up by 30%, the spatial distribution of the TAWSS is preserved whereas the local amplitude is reduced.

Plaque progression with age
To our knowledge, the progression of plaques in time, under in vivo conditions, is not well documented in the literature, particularly for the aortic arch. In this article, data obtained from a preliminary study were used to indicate the progression of plaques in mice fed for 10, 12, 14 and 16 weeks under a fat diet. Fig. 4 shows the aortic arch plaque distribution at these different stages. Even if a larger number of samples were necessary to extract statistic quantitative measurements, it is already noticeable that 5 main plaques (indicated by arrows in Fig. 4b) are observed in those images corresponding to plaques in the main curvature (one in the AA and one in the DA) and in the opposite surface to branch bifurcation (in BCA, LCC, LS). Because lumen restriction can be important for plaque stability determination [42], Table 1 indicates the maximum blockage in % on the different branches for all the samples. This maximum blockage corresponds to the ratio between plaque area and healthy lumen area in a section perpendicular to the skeleton of the geometry. Whereas plaque size globally increases in time, no real tendency can be extracted from Table 1 as the plaque progression rate is animal dependent, and because some of the plaques have been cut during the dissection so that more samples would be necessary to assess the evolution of the maximum blockage in time. Although some information is missing, some general features already appeared and can be compared to the literature.
One of the few existing 3D time studies on this region of the arterial tree was conducted on the BCA [51]. In this article, Mc Ateer et al. indicate that the evolution of the BCA plaques in ApoE −/− mice occurs continuously (3 time points used) from the proximal to the distal region of the BCA and that the plaques are observed preferentially on the right Table 1 Maximum blockage in % due to the presence of plaques for the ascending aorta (AA), brachiocephalic artery (BCA), left common carotid artery (LCC), left subclavian artery (LS) and descending aorta (DA). The blockage was calculated at the point where maximum plaque protrusion occurs in each branch. It is based on measurements of cross sectional area in the diseased and healthy vessel. The cross section is perpendicular to the skeleton branches. The initials NA mean not available as this plaque does not appear on the geometry (Fig. 4e). side of the artery, a fact also observed in the present work ( Fig. 4b-e). In addition, it is noticeable that the plaque along the aortic arch is mainly present in the ascending aorta and is not centred along the arch, a phenomenon already observed by Maeda et al. [52] in a 2D study comparing plaque development in ApoE with different genetic backgrounds. According to Zhu et al. [23], the variability in the plaque distribution a b  between C57BL/6 and 129/SvEv comes from the difference in the geometrical features and consequently the wall shear stress distribution.

AA
In the C57BL/6 background mice, a local maximum of wall shear stress is observed close to the highest curvature point of the arch, which could explain the spatial limitation of the plaque in this zone. A correlation between the TAWSS and the plaque size has been shown in Tomita et al. [39], but as indicated by the authors, the TAWSS is estimated using velocity average values so that the knowledge of the local distribution of the WSS is necessary. This local distribution can be calculated using numerical methods as detailed in the following. The geometry used in the next section corresponds to Fig. (4b).
The geometry used in the next section corresponds to Fig. 4b. Because the dissection practice and the aortic arch geometry can be animal-dependent, some discrepancy can be observed in the size and the regularity of the samples.

Computational fluid dynamics calculations
In the following, a particular attention will be paid to the flow structure in the aortic arch, specifically in the brachiocephalic region where mouse plaques are the most vulnerable and whose biological characteristics are the most similar to human plaques [53]. In addition, the change in the blood flow due to the presence of the plaque in the inner curvature of the arch is examined, as the presence of this plaque is known to increase considerably the risk of emboli in cardiac surgery patients [44].
The complexity of blood flow in the aortic arch comes not only from the complexity of the geometry (bifurcations and curvature) but also from the potential high Reynolds number (Re = ρUD / μ, where D is the average diameter) that can be found in this region of the arterial tree. In humans, the mean Reynolds number in the aorta is 1200 and the instantaneous Reynolds number can reach values up to 6500 [54]. In this range of parameters, the blood flow can become turbulent. In mice, the mean Reynolds number in the different branches of the aorta remains under 100 so laminar (although unsteady) blood flow is expected [55]. While being laminar, this flow can nonetheless be complex. Fig. 5 shows the swirling instantaneous streamlines obtained at the systole (maximum of flow rate).
The presence of plaques (Fig. 5b), compared to the healthy configuration (Fig. 5a), induces more complexity in the flow (distortion of the instantaneous streamlines) and accelerate locally the flow (indicated in red on the streamlines Fig. 5b). In addition, the presence of plaques may generate recirculation zones as indicated in Fig. 6b behind the plaque at the inner curvature of the aorta. These recirculation zones change the blood flow locally and may have an impact on the preferential direction for the plaque growth. Indeed, correlating this result with Fig. 4 and the results of Maeda et al. [52], this region seems to be prevented from plaque formation; however, images for mice after longer times under fat diet would be necessary to assess the relation between the recirculation and the plaque expansion in time. The presence of these recirculation zones is related to the presence of separation and stagnation regions in the flow. Whereas those regions can be well determined by locating the null WSS zones when the flow is steady [56], they are not well defined when the flow is unsteady [57]. To locate the presence of a recirculation, it is however possible to calculate the instantaneous WSS (Fig. 7a) and to locate the region where the WSS is null. It is noticeable that the TAWSS (Fig. 7b) shows a region of very low amplitude around a similar position. Other parameters such as the oscillatory shear index OSI ¼ 0:5 1− ∫ T 0 τ w dt T =TAWSS or the relative residence time (RRT = 1/(1 − 2OSI) * TAWSS) can be used to locate the recirculation zone. These parameters are represented Fig. 7c and d for the geometry with plaque. In particular, the OSI, which represents a measure of the fraction of one period during which the sign of the WSS is reversed, gives a good indication of the recirculating zone. It takes values in the range of 0 to 0.5, with 0 corresponding to unidirectional flow, and 0.5 to purely oscillatory flow. The RRT has been found particularly high in two points around the sharp edge of the plaques. These regions can be described as a sort of separation points for a time averaged flow. In the present geometry, only the plaque generated at the arch generates a circulation zone. In addition, it is noticeable that contrary to bigger animals (rabbits or dogs), no separation region is observed in the healthy configuration (Fig. 6a), indicating that the presence of separation and stagnation regions is not initially necessary for the formation of plaques as suggested in some studies [56].
In addition to the previous results, the presence of plaques was found to alter the TAWSS distribution on the vessel wall, a phenomenon particularly visible on the BCA (Fig. 8). The magnitude of the highest WSS is increased significantly, and the position of local high WSS is altered. The high WSS is located on the upstream face of the brachiocephalic artery, corresponding to the site of maximal luminal narrowing. The WSS distribution on the brachiocephalic branch at different stages throughout the cardiac cycle is shown in Fig. 9. This highlights the variation in WSS that occurs throughout the cardiac cycle, in addition to the variation induced by the presence of the plaque. Several studies suggest that plaque rupture may be associated with high WSS [58,59], and hence the local maximum of WSS on the surface of the brachiocephalic plaque could contribute to plaque rupture. It is, however, likely that structural stresses are a more important indicator for predicting plaque rupture [60]. Such structural stresses can be investigated using fluid structure interaction simulations, a method to which the present code can be adapted.
The modelling of such a complex physical system has involved some simplifications (rigid walls, no animal specific flow rate boundary conditions, limited number of samples) that may limit the reach of our conclusions but which we aim to address in the future.

Conclusions
Given the limited numerical data on blood haemodynamics in the presence of atherosclerosis, the main goal of this study was to develop methods for modelling pulsatile blood flow in a real artery with plaque. An X-ray synchrotron imaging technique has been successfully applied to follow up the plaque progression in ApoE−/− mice. Whereas general features have been observed, such as regions prone to or prevented from plaque development, a larger number of samples would be necessary to assess the preferential direction of plaque progression and to relate it to haemodynamical parameters.
Concerning the flow distribution, the presence of plaques was found to increase WSS significantly, particularly in regions with large luminal narrowing, such as the brachiocephalic artery. It is also shown that the plaque development alters the blood flow and can induce recirculation zones. In particular, these changes in the flow could lead to phenomena such as delamination [13] or erosion [19], both identified as related to plaque rupture; however, further study would be necessary to relate those phenomena to WSS distribution.