## Abstract

A high concentration of low-density lipoprotein (LDL) is recognized as one of the principal risk factors for development of atherosclerosis. This paper reports on modelling and simulations of the coupled mass (LDL concentration) and momentum transport through the arterial lumen and the multi-layered arterial wall of an anatomically realistic carotid bifurcation. The mathematical model includes equations for conservation of mass, momentum and concentration, which take into account a porous layer structure, the biological membranes and reactive source/sink terms in different layers of the arterial wall, as proposed in Yang & Vafai (2006). A four-layer wall model of an arterial wall with constant thickness is introduced and initially tested on a simple cylinder geometry where realistic layer properties are specified. Comparative assessment with previously published results demonstrated proper implementation of the mathematical model. Excellent agreement for the velocity and LDL concentration distributions in the arterial lumen and in the artery wall are obtained. Then, an anatomically realistic carotid artery bifurcation is studied. This is the main novelty of the presented research. We find a strong dependence between underlying blood flow pattern (and consequently the wall shear stress distributions) and the uptake of the LDL concentration in the artery wall. The radial dependency of interactions between the diffusion, convection and chemical reactions within the multi-layered artery wall is crucial for accurate predictions of the LDL concentration in the media. It is shown that a four-layer wall model produced qualitatively good agreement with the experimental results of Meyer *et al*. (1996) in predicting levels of LDL within the media of a rabbit aorta under identical transmural pressure conditions. Finally, it is demonstrated that the adopted model represents a good initial platform for future numerical investigations of the initial stage of atherosclerosis for patient-specific geometries.

## 1. Introduction

The process of forming atheromous plaques in the inner lining of arteries (atherosclerosis) is an extremely complex biochemical process and its exact origin still remains unknown. Atherosclerosis primarily occurs in large- and medium-sized arteries, such as the coronary and carotid arteries, especially at bifurcations and curvatures [1,2]. It is now widely accepted that atherosclerosis can be seen as an inflammatory disease resulting from endothelial dysfunction caused by injury of the endothelial layer. Possible causes of injury of the endothelial layer include an elevated low-density lipoprotein (LDL) concentration, hypertension, diabetes, obesity, smoking and infection [3]. Although a series of cellular and molecular responses are involved in the process of developing atherosclerosis, the uptake of LDL in the artery wall plays a key role, because the uptake of LDL in the artery wall and its partial oxidation will enhance the inflammatory response [2]. It is found that local haemodynamics also play a significant role in atherosclerosis [4]. Additionally, it is experimentally demonstrated that low wall shear stress (WSS), oscillating shear stress and long residence-times of the blood influence the formation of atherosclerosis [5].

In the modelling of the transport of macromolecules in arteries, such as LDL, a few categories of models can be identified: the wall-free, the lumen-wall and the multi-layer models [6]. In the wall-free models, only the lumen region is directly simulated, whereas the influence of the arterial wall is included through imposed boundary conditions. It is obvious that such models cannot provide any information regarding distributions of requested variables (filtration velocity or concentration) within the arterial wall (e.g. [7]). The lumen-wall models also consider the arterial wall structure. Here, the complex multi-layer heterogeneous structure of the arterial wall is approximated by a single homogeneous layer (e.g. [8]). This mono-layer approximation significantly reduces the complexity of the model and it can give reasonably accurate concentration distributions within the artery wall, especially when applied with the multi-pore approach across the lumen–endothelium interface (e.g. [9]). The most complex models are the multi-layer models, which take into account several wall layers with specific physiological properties (e.g. [6]). Obviously, these multi-layer models are geometrically the most realistic because they mimic the internal arterial wall structure in detail. However, the multi-layer models also require the largest number of transport parameters and quite challenging numerical mesh generation.

In this study, we focus on the multiple-layer wall model with constant thickness as proposed in earlier studies [10–12]. We performed an integral (by using a single computational code) steady-state analysis of combined flow and mass transfer. In contrast to the earlier works [10,11,13], in which simplified two-dimensional geometries were considered, we extend our analysis also to an anatomically realistic carotid artery bifurcation.

The patient-specific simulation of LDL accumulation in a human left coronary artery is studied in [14]. Similarly, computational modelling of LDL transport an *in vivo* computed tomography (CT) image-based human right coronary artery is analysed in [15]. However, a homogeneous (single-layer) wall model was applied in both studies.

Our goal is to provide a feasibility study of simulating transport of the LDL into the multi-layered arterial wall of an anatomically realistic carotid artery bifurcation. This represents the first step towards *a priori* identification of regions in the vascular system where initial development of atherosclerosis could take place.

## 2. Mathematical model

We will present separately equations for the lumen (blood flow) and for the artery wall (porous layers and membranes) regions.

### 2.1. Lumen

Blood is considered to be an incompressible Newtonian fluid in the lumen. This is valid for the typical values of the Reynolds number for a blood flow in the carotid artery (i.e. *Re* = *𝒪*(10^{2})). Its behaviour is described by conservation of mass and momentum as
2.1
2.2where the LDL concentration is modelled with the standard convective–diffusive mass transport as
2.3Here, * ρ* is the density,

*the dynamical viscosity and*

*μ**𝒟*the diffusion coefficient.

### 2.2. Arterial wall

Additional equations for transport of the blood plasma and of the LDL concentration through a multi-layer artery wall are introduced next. They are based on a combination of equations for the porous layer and for biological membranes, with a chemical-reaction source/sink term [10]. The flow and mass transport through biological membranes are modelled by using the Staverman–Kedem–Katchalsky (SKK) expressions (as given in [16,17]) and are written as
2.4
2.5where *J _{v}* is the total flow per unit area,

*J*

_{c}the concentration flux per unit area,

*K*′ is the permeability (per unit length), the effective diffusivity (per unit length),

*Δ*

*p*is the pressure difference across the membrane,

*Δ*

*=*

*π**RT*

*Δ*

*c*the osmotic pressure difference,

*R*the universal gas constant,

*T*the absolute temperature,

*σ*_{d}and

*σ*_{f}the Staverman osmotic reflection and filtration reflection coefficients, respectively. In this study, it is assumed that

*σ*_{d}=

*σ*_{f}, similar to [10,11]. The blood plasma in the arterial wall is also considered to be an incompressible Newtonian fluid, however with a different viscosity than blood in the lumen. The volume-averaged (within a porous layer) continuity equation for the plasma can be written as 2.6Note that ‘〈〉’ stands for the volume-averaged variables. To account for the effects of the biological membranes, the volume-averaged momentum equations need to be extended with the above-introduced SKK membrane expressions. The final form of the volume-averaged momentum equation can be written as 2.7where

*is the porosity (defined as a fraction of the void volume over the total volume;*

*ε**= 1 for entirely free fluid and*

*ε**= 0 for the total flow blockage). The third term on the right-hand side is the flow resistance by the solid matrix of a porous medium given by Darcy's law and the fourth term is the osmotic pressure term. Note that the presence of the osmotic pressure term makes a fully two-way coupling between the momentum and concentration equations within the artery wall.*

*ε*Similarly, the volume-averaged convective–diffusive–reactive mass transport equation can be written as
2.8where *𝒟*_{e} is the effective diffusivity and *k*_{react} is the effective volumetric reaction rate coefficient. In order to assure continuity of the convective–diffusive transport across the interfaces of different layers of the arterial wall, the following boundary condition is used
2.9where *v _{n}* is the velocity perpendicular to the interface. Note that for the lumen–endothelium interface,

*σ*_{f}= 0 on the lumen side and

*𝒟*

_{e}is the diffusivity coefficient of the LDL in blood.

A schematic of the multi-layered arterial wall structure is shown in figure 1.

In this study, we adopted the four-layer model, consisting of the endothelium, intima, internal elastic layer (IEL) and media. The endothelial luminal glycocalyx layer (‘sugar coat’) that covers the plasma membrane of a single layer of endothelial cells (with a typical thickness of 60 nm) is not considered in this study owing to its negligible effects [10].

An overview of the characteristic thickness of the four selected artery wall layers is shown in table 1. The endothelium and IEL are modelled as a combination of the porous media and biological membranes, i.e. the effects of the osmotic pressure are taken into account for these layers. In contrast to that, the intima and media are treated as (macroscopically) homogeneous porous media and the osmotic pressure effects are negligible, as shown in [18]. The LDL uptake in the media is modelled using a first-order chemical reaction (the last term in equation (2.8)). Instead of calculating mass transport and flow of plasma in the adventitia, a simple zero-gradient boundary condition is imposed at the media–adventitia interface, i.e. . The characteristic physiological parameters of each individual wall layer used in the numerical simulations are listed in table 2.

### 2.3. Numerical method

The system of equations (2.1)–(2.8) is numerically discretized by employing a finite-volume approach for general non-structured geometries, Fluent/Ansys Inc. 12.0. The diffusive terms in transport equations are discretized by the second-order central differencing scheme (CDS). The bounded second-order CDS is used for convective terms of all equations. The coupling between the velocity and pressure field is achieved by using a standard iterative SIMPLE algorithm. The flow and mass equations (equations (2.6)–(2.8)) within the arterial wall are programmed and inserted via the user-defined functions into the main lumen flow solver. The simulations are executed in the parallel mode by employing the domain-decomposition MPI directives. For typical runs, 4–8 CPUs are used on the local Linux cluster.

## 3. Results and discussion

### 3.1. Simplified geometry

The numerical implementation of the mathematical model given by equations (2.1)–(2.8) is firstly validated against simplified straight two-dimensional (polar-cylindrical) and three-dimensional (full cylinder) geometries and results are compared to available data in literature [10]. The characteristic radius of the simplified artery is *R* = 3.1 × 10^{−3} m and its length is *L* = 0.125 m. The multi-layers within the artery wall are generated to match exactly the characteristic thickness per layer given in table 1. A schematic overview of the geometry under consideration, with boundary conditions, is shown in figure 2. At the inlet, a fully developed laminar velocity profile is imposed, , with *U*_{0} = 0.338 m s^{−1}, giving *Re* = 100. The non-dimensional LDL concentration at the inlet is *c* = 1. On the axis of symmetry (lumen side), the zero-gradient boundary conditions are applied for the concentration and tangential velocity component, while the normal velocity component is zero (the symmetry boundary conditions). At the lumen outlet, the zero gradient is imposed for the mass transport (∂*c*/∂*n* = 0) and the outlet pressure is specified to be *p* = 100 mmHg (13332.2 Pa). At the media–adventitia interface, the pressure boundary condition with *p* = 30 mmHg is imposed for the momentum, and a zero-gradient boundary condition for concentration. This gives the transmural pressure gradient of 70 mmHg, which is close to the typical physiological conditions within a realistic artery. At the inlet and outlet section of the arterial wall layers, the (endothelium–media section) symmetry boundary conditions are applied for both momentum and concentration.

To validate implementation of the momentum transfer within the porous layer, the radial velocity profiles along the lumen–endothelium interface in the streamwise direction are compared to results presented in [10]. For example, for the monitoring point at *x*/*L* = 0.5 after the artery inlet, the obtained value of the radial velocity component of *v _{r}* = 2.29 × 10

^{−8}m s

^{−1}agrees very well with

*v*= 2.31 × 10

_{r}^{−8}m s

^{−1}(at the same location) as reported in [10]. The quality of the agreement is also very good for other locations (the difference is always less than 2%), proving the validity of the velocity calculations in the porous layer. The profiles of the non-dimensional LDL concentration (

*c*/

*C*

_{0}) along the same lumen–endothelium interface are shown in figure 3

*a*. To check the grid dependency, three different numerical grids are used, including approximately 10

^{4}, 2 × 10

^{4}and 4 × 10

^{4}control volumes, with 200, 400 and 800 control volumes applied in the radial direction, respectively. The endothelium and IEL are represented by 4, 12 and 30 control volumes for the coarse, intermediate and fine numerical mesh, respectively. It can be seen that a good agreement between the present results and results of Yang & Vafai [10] are obtained. It is important to note that simulation revealed the concentration polarization (or Staverman filtration mechanism) (an increase in the concentration levels compared with its inlet values) at the lumen–endothelium interface (figure 3

*a*). This is caused by the higher mass transfer resistance imposed by the endothelium, as demonstrated in [25–28]. A radial profile of the LDL concentration (started at

*x*/

*L*= 0.5 from the artery inlet) is shown in figure 3

*b*. The characteristic radial profile shows that uptake of LDL in the multi-layered arterial wall exhibits a complex behaviour. Both endothelium and IEL act as a buffering layer which generates the concentration polarization at the lumen–endothelium and intima–IEL interface, respectively. The largest decrease in the radial transfer of the LDL concentration occurs across these layers and it is important to apply a sufficiently fine numerical mesh in order to properly capture concentration gradients. It can be seen that values of the LDL concentration in the media dropped to values that are two orders of magnitude smaller, compared with the initial level at the lumen–endothelium interface. Finally, uptake of the LDL in the media is taken into account through the simple linear chemical reaction (with a negative chemical-reaction rate coefficient,

*k*

_{react}), reducing these values further. By comparing the radial profiles of the LDL concentration obtained and of results presented in [10], it can be concluded that a good agreement is obtained. Note that this good agreement is obtained despite the fact that different discretization techniques (our finite volume versus the finite element of Yang & Vafai [10]), different numerical mesh distribution and discretization schemes and solvers are used.

As the next step in the validation of the numerical method, a three-dimensional cylindrical geometry is constructed (here a segment containing one-eighth of the full cylindrical geometry is considered) with an interchanged coordinate system (the *z*-direction is now the streamwise direction) and results compared with the two-dimensional axi-symmetrical case already discussed. Although this step seems quite simple, it is an important test to validate proper model implementation for a fully three-dimensional non-orthogonal geometry. The total number of the control volumes used is 6.8 × 10^{5}, where 200 control volumes were employed in the radial direction. The calculated radial velocity at the same location as mentioned previously for the two-dimensional geometry at *x*/*L* = 0.5 after artery inlet showed a negligible difference of 1%. Similarly, the LDL concentration levels at the centre of the intima showed a difference of less than 1.5%. This confirms a proper implementation of the model for the momentum and LDL concentration transfer in the multi-layered arterial wall for three-dimensional geometries.

### 3.2. An anatomically realistic carotid artery bifurcation

#### 3.2.1. Geometry and numerical mesh

Next, we consider an anatomically realistic carotid artery bifurcation. The smoothed segmentation of a carotid bifurcation of a healthy person was provided by Erasmus Medical Center, Rotterdam, and has been constructed from the CT scans (F. J. Gijsen 2010, private communications; figure 4*a*). The characteristic circumferentially averaged diameters of the common carotid artery (CCA), internal carotid artery (ICA) and external carotid artery (ECA) are approximately 6.5, 4 and 3.2 mm, respectively. This geometry still needs to be modified by clipping the end caps and by adding flow extensions to the inlet and outlets (figure 4*b*). These flow extensions make it easier to apply appropriate boundary conditions to the inlet and outlets. This is done by using the Vascular Modelling Toolkit developed by Antiga *et al*. [29]. Despite the constant thickness assumption, owing to a high aspect ratio between different layers, it is quite challenging to generate the patient-specific numerical mesh of a good quality for the problem under consideration. Here, we will provide some practical guidelines for an efficient mesh generation containing both the lumen and the multi-layered arterial wall. The geometry of the carotid bifurcation is provided in the STL (STereoLithography) file format, widely used in stereolithography and rapid prototyping. The STL file format identifies only the surface of a three-dimensional geometry. The surface is triangulated and described by the vertices and unit normal of the triangles within a predefined Cartesian coordinate system. In our first attempt at a proof-of-concept demonstration of the multi-layered patient-specific arterial wall mesh generation, we assumed constant thickness of the layers. This implies that relatively simple offsetting of the lumen-wall surface in the wall-normal direction can serve as a basis for generating additional layers. However, simple offsetting of all triangles in their outward normal direction will result in unconnected and overlapping surfaces. To solve this problem, instead of the cell-faces, the vertices are translated in the wall-normal direction. The translation vector is calculated by averaging normals of the faces connected to each vertex. A large offset can cause an intersection of triangles at the bifurcation. To localize and to remove these intersections, a remeshing is performed to fill and smooth gaps created by removing the intersecting triangles.

To properly capture the expected large gradients of the LDL concentration in the endothelium and IEL layers, a refined numerical mesh must be applied in these regions. The computation mesh that includes the lumen and artery wall layers is shown in figure 5. At least eight control volumes are required to properly resolve the endothelium and IEL layers. The final numerical mesh consists of approximately 2 × 10^{6} hexagonal control volumes that was proven to be sufficient to obtain grid-independent results (where the second-order quadratic upwind and the second-order linear upwind discretization schemes are used for convective terms of momentum and concentration equations, respectively).

#### 3.2.2. Boundary conditions

The pressure boundary conditions need to be applied at outlets to ensure medically recorded 0.6 : 0.4 (ICA/ECA) mass-flow ratios for the present carotid artery ([30], figure 6). To obtain correct outlet pressures, the mass-flow ratio boundary conditions are imposed and simulation is performed. Then, the pressure difference between outlets of the internal and external arteries can be easily obtained. By specifying the pressure outlet of the ICA as the reference pressure of 100 mmHg, the value of the pressure correction *Δ**p* at the outlet of the ECA can be easily adjusted. The adjusted pressure difference of *Δ**p* = 6.35 Pa at the outlet of the ECA is used in this study. Also, the pressure values along the media–adventitia interface are specified to be 30 mmHg. The fully developed laminar flow profile is imposed at the inlet with the typical Reynolds number of *Re* = 105. The reference LDL concentration at the inlet is specified to be *c*_{0} = 2.5 mol m^{−3} (a recommended normal value for an adult person). The zero-gradient boundary conditions are applied for the LDL concentration at the internal and ECA outlets, and along the media–adventitia interface. An isothermal system is assumed with a reference temperature of 37°C. A summary of all boundary conditions for the present geometry is shown in figure 6. A steady-state simulation is performed.

#### 3.2.3. Results

To portray the most salient features of the blood flow in the lumen, the streamtraces coloured by the velocity magnitude are shown in figure 7*a*. It can be seen that, in contrast to a simplified geometry, a complex three-dimensional helical flow pattern is generated at the bifurcation, followed by a rapid recovery farther downstream. The non-dimensional (relative) LDL concentration contours at the lumen–endothelium interface are shown in figure 7*b*. Non-uniform distributions of the LDL concentration are obtained (as expected) owing to the local dependency of the mass transfer and the underlying flow pattern inside the lumen. It can be seen that a maximum increase of 3.5% of the LDL concentration is obtained at the carotid sinus. As discussed in the previous section, values of *c*/*C*_{0} larger than one indicate the polarization. The values observed here are in a good agreement with values reported in literature. For instance, values of the concentration polarization of 1.026 are reported in [6,11,31], albeit for the relatively simple (straight) arteries, but for identical transmural pressure gradient (i.e. the radial pressure difference between the lumen side and the end of the intima) of 70 mmHg. In the study of Yang & Vafai [10], it is demonstrated that for higher values of the lumen outlet pressure of 160 mmHg (that mimic effects of hypertension), values of the concentration polarization increased up to 7.5%. Similarly, in their computational study of LDL transport in a human right coronary artery, Sun *et al*. [15] observed that although the LDL concentration varied along the artery, the concentration polarization was not more than 6% of the inlet concentration. A higher transmural pressure of 120 mmHg and shear-dependent endothelium permeability is analysed in [15]. This non-constant permeability of the endothelium can be one of the reasons for slightly enhanced values of the LDL polarization. In their simulations of a patient-specific left coronary artery, Olgac *et al*. [14] introduced a multi-pore-based model for the mass transfer of the LDL. Here, the endothelium is represented by a three-pore model, which includes local WSS effects on the LDL transfer. Their results showed polarization levels up to 20%. These increased levels were the consequence of the shear-dependent permeability of the endothelium.

The arterial cross sections at three characteristic locations, shown in figures 8*a*, 9*a* and 10*a* are analysed next. The first cross section (1-1) is at the carotid sinus, i.e. the location where the maximum of the concentration polarization is observed (figure 8). The axial velocity (perpendicular to the cross section) shows a single peak distribution (figure 8*b*). The intensity of the secondary flow is depicted by contours of the axial vorticity (figures 8*c*, 9*c* and 10*c*) and by superimposed two-dimensional velocity vector projections and corresponding streamtraces (figures 8*d*, 9*d* and 10*d*). By comparing distribution of the contours of the axial vorticity (* ω_{z}*) (figure 8

*c*) and of the non-dimensional concentration within the artery wall (figure 8

*d,e*), it can be seen that the higher values of the LDL concentrations are obtained along the lumen–endothelium interface where gradients of the axial vorticity reach their extreme values. Similarly, at the cross sections 2-2 and 3-3, shown in figures 9 and 10, lower values of the LDL concentration are observed along the locations where the axial vorticity is close to zero, as shown in figures 9

*c*,

*e*and 10

*c*,

*e*.

To establish additional correlations between the local flow and concentration distributions, three arbitrarily selected profiles along the lumen–endothelium interface in the streamwise direction are extracted (figure 11*a*). Note that ‘line 1’ is extracted along the CCA, ‘line 2’ along the ICA and ‘line 3’ is extracted along the ECA. The resulting profiles of the WSS and of the non-dimensional LDL concentration are shown in figure 11*b*,*c*. Note that the horizontal axis shows a non-dimensional total distance between the beginning and the end of the line. Starting with ‘line 1’, it can be seen that a local minimum of the WSS takes place at *L/s* = 0.7 and a local maximum at *L*/*s* = 1. At the same location, profiles of LDL concentration reach local maximum and local minimum, respectively. Similarly, for ‘line 2’, the local maxima of WSS are at *L*/*s* = 0 (also for line 3) and *L*/*s* = 0.84. At these locations, the LDL concentration reaches its local minimum. It can be concluded, that many regions of low WSS correspond to regions of high LDL concentration along the lumen–endothelium interface. Similar conclusions were drawn in [7] where, generally speaking, it is stated that the points of the lowest WSS and highest surface concentration of LDL mostly coincide. In addition, they found that local surface values of the LDL concentration sharply increase as the WSS decreases from 1 to 0 Pa. We observed a similar trend for a characteristic peak of the surface concentration for the ‘line 1’ at *L*/*s* = 0.7 (figure 11*c*). The connection between WSS and surface concentration of the LDL along the same lines is shown in figure 11*d*. It can be seen that there is a significant scatter in distributions indicating multiple values of *c*/*C*_{0} for a given value of WSS. This confirms that the surface LDL concentration does not depend solely on WSS but also that additional parameters can play a significant role in these regions. Similar conclusions were drawn in [7,32]. Potential additional parameters can be previously analysed secondary flow patterns, shown in figures 8⇑–10.

The analysis of the radial profiles of the LDL concentrations is performed next. For two characteristic locations (P1 and P3), the entire radial profiles are extracted and shown in figure 12. It can be seen that the maximum resistance of the transmural transport of the LDL takes place in the membrane-like porous layers, i.e. endothelium and IEL. The local dependency of distributions is clearly visible, with values of the LDL in the intima twice as high for the P1 location in comparison with the P3 location. The radial profiles within the multi-layered arterial wall at selected points (P1–P4) are shown in figure 13. Here *r** is the non-dimensional radial distance that extends from the lumen–endothelium interface (*r** = 0) to the media/adventitia interface (*r** = 1). Also experimental data of Meyer *et al*. [31] for a rabbit aorta under identical transmural pressure of 70 mmHg are plotted too. It can be seen that a good agreement between numerical distributions and experimental data is obtained. Note that this is just a qualitative comparison demonstrating similar behaviour of the measured and simulated LDL distributions within the media.

It can be concluded that the model is able to capture differences in spatial distributions of the LDL concentration within the arterial wall in the realistic carotid bifurcation caused by local changes in WSS along the lumen–endothelium interface. The obtained results show qualitatively good agreement with experimental data of Meyer *et al*. [31] for the LDL levels within the arterial wall of a rabbit aorta.

## 4. Conclusion and outlook

It is demonstrated that a multi-layer wall model for transport of the LDL through the artery wall, proposed by Yang & Vafai [10] and tested on relatively simple two-dimensional geometries, can be successfully applied to the patient-specific carotid artery case. Despite challenging requirements in terms of the relatively complex mesh generation and relatively high number of required transport parameters for each particular layer, it is concluded that the four-layer (endothelium–IEL–intima–media) wall model represents a good basis for future numerical investigations of atherosclerosis development for patient-specific conditions. Despite this, there are still many segments of the mathematical model that can be improved. The first step will be to replace the single-pore transfer mechanism considered here with the multi-pore theory where the permeability of the endothelium is expressed as a function of the local WSS conditions [9]. This will make a unique combination of the multi-layer and multi-pore arterial wall model (we are currently developing such an integrated approach). The here analysed four-layer model can be expanded to include the glycocalyx layer effects, similar to [33]. Another limitation is that this study treated arterial walls as rigid. Additional effort is needed to have a fully coupled time-dependent fluid flow, mass transfer and flexible structure (multi-layered arterial wall) interactions, as shown in [34] for simplified geometries. The next step can be to have the variable thickness of the wall layers, which can be mimicked relatively easily with the presented meshing algorithm, similar to [35]. Future improvements should also include the inflammatory process, which requires modelling of the biochemical interactions that include transfer of the LDL, oxygen and foam cells [36]. This will result in a comprehensive integrated model of atherosclerotic plaque formation.

## Acknowledgements

Dr Frank Gijsen, Erasmus Medical Centre, Rotterdam, The Netherlands is acknowledged for providing geometry of the patient-specific carotid artery used in this study and for inspiring discussions.

- Received October 14, 2013.
- Accepted November 1, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.