Abstract
Atherosclerosis is a vascular disease caused by inflammation of the arterial wall, which results in the accumulation of lowdensity lipoprotein (LDL) cholesterol, monocytes, macrophages and fatladen foam cells at the place of the inflammation. This process is commonly referred to as plaque formation. The evolution of the atherosclerosis disease, and in particular the influence of wall shear stress on the growth of atherosclerotic plaques, is still a poorly understood phenomenon. This work presents a mathematical model to reproduce atheroma plaque growth in coronary arteries. This model uses the Navier–Stokes equations and Darcy's law for fluid dynamics, convection–diffusion–reaction equations for modelling the mass balance in the lumen and intima, and the Kedem–Katchalsky equations for the interfacial coupling at membranes, i.e. endothelium. The volume flux and the solute flux across the interface between the fluid and the porous domains are governed by a threepore model. The main species and substances which play a role in early atherosclerosis development have been considered in the model, i.e. LDL, oxidized LDL, monocytes, macrophages, foam cells, smooth muscle cells, cytokines and collagen. Furthermore, experimental data taken from the literature have been used in order to physiologically determine model parameters. The mathematical model has been implemented in a representative axisymmetric geometrical coronary artery model. The results show that the mathematical model is able to qualitatively capture the atheroma plaque development observed in the intima layer.
1. Introduction
Cardiovascular diseases related to atherosclerosis continue to be the leading cause of death worldwide [1]. Atherosclerosis is the process in which plaques—consisting of deposits of cholesterol and other lipids, calcium and large inflammatory cells called macrophages—are built up in the walls of the arteries causing narrowing (stenosis of the lumen), hardening of the arteries and loss of elasticity. This leads to a reduction in the blood flow through the vessels. Nevertheless, the most serious damage occurs when the plaque becomes fragile and ruptures (vulnerable plaque). Plaque rupture causes the formation of blood clots that can block blood flow or break off and travel to another part of the circulatory system, thus producing heart attacks, strokes, difficulty in walking and eventually gangrene [2–4].
The mechanical factors that could initiate of atherosclerosis lesions have been widely explored by many authors [5–7]. Cyclic stretch, laminar and oscillatory shear stress, effects of vessel compliance, curvature, pulsatile blood flow or cardiac motion are considered the main triggers of atherosclerosis initiation. Atherosclerosis is directly related to nanoscale fluid dynamics and macromolecular transport at the arterial endothelium [8]. The behaviour of atherosclerosis is timedependent, and the dominant phenomena may occur in the range of 0–100 Å from the endothelium [9]. The location of atherosclerosis is associated with flow separation and turbulence. Therefore, many studies have identified haemodynamic shear stresses as an important determinant of endothelial function and phenotype in atherosclerosis disease (see [10,11] and references therein). Wall shear stress (WSS), which is defined as the force per unit area exerted on the vessel wall by the blood flow, is related to the plaque rupture, because this phenomenon has been more frequently observed at the proximal upstream side of the minimal lumen diameter, which is supposedly exposed to higher WSS. Then, arteriallevel shear stress approximately greater than 1.5 Pa induces endothelial quiescence and an atheroprotective gene expression profile [10]. On the other hand, atherosclerosis disease preferentially affects predisposed areas of the blood vessel where the haemodynamic shear stresses are low such as bifurcations or zones with flow recirculation. Shear stress lower than approximately 0.4 Pa stimulates an atherogenic phenotype and eases atherosclerotic lesion initiation [10].
The onset and progression of atherosclerosis involves many processes. However, the fundamental cause of plaque development is believed to be the abnormal enlargement of the intima by the infiltration and accumulation of macromolecules such as lipoproteins and the associated cellular and synthetic reactions. It is well accepted that early atherogenesis tends to be hallmarked by an abnormally high accumulation of macromolecules, i.e. LDL, within the arterial wall. This suggests that the macromolecular transport in the arterial wall must have some impact on the initiation and development of atherosclerosis [12]. Therefore, the general picture is that zones of elevated LDL tend to localize with areas of atherosclerotic lesion development and intimal thickening, and such zones also tend to have ‘abnormal’ WSS patterns. The endothelium, which is the internal lining of the entire cardiovascular system, is uniquely positioned to be a sensor, responding and transducing these haemodynamic signals. Taken together, these results suggest, but do not prove, a role for mass transport in atherogenesis [13]. Furthermore, fatty streaks, the earliest detectable lesions in atherosclerosis, contain macrophagederived foam cells that differentiate from recruited blood monocytes. Monocytes are recruited to tissues via constitutive signals and in response to inflammatory mediators [14]. More advanced atherosclerotic lesions, called fibro fatty plaques, are the result of continued monocyte recruitment, together with smooth muscle cell (SMC) migration and proliferation [15]. In addition, chemokines or chemoattractant cytokines constitute a family of over 40 different cell signalling molecules important for constitutive trafficking and recruitment of leucocytes in response to inflammatory mediators. Some chemokineses that can act as potent mediators of monocyte and SMC migration and macrophage differentiation are expressed in human atherosclerotic lesions [16,17].
The Staverman–Kedem–Katchalsky membrane equations [18] are usually used to model the transport processes in the endothelium and internal elastic lamina (IEL) in the majority of the previous multilayer models. However, these equations have at least two substantial disadvantages despite the fact that they are widely used to model the mass transport through the endothelium: (i) they are derived based on the existence of a steadystate condition, which is in conflict with the real physiological conditions existing within both the endothelium and IEL, and (ii) they ignore the boundary effects on the flow across the membrane, which is not valid when the boundaries of the porous membrane have to be accounted for [12,19]. Furthermore, there are three transport pathways through the endothelium: vesicular transcytosis, which is regulated by receptors on the endothelial cells [20], and through leaky and normal junctions, most of which are located at the sites of dying or replicating cells [21]. Apart from the LDL molecules, many studies have focused on the governing mechanics interaction of the different biological species which play a role in atheroma plaque development from both experimental (see [22–24], among others) and computational (see [25–28], among others) points of view. Furthermore, there are greatly varying degrees of complexity in these computational studies depending on the number of species considered and the development of the equations proposed. Zohdi et al. [29] modelled the adhesion of monocytes to the endothelial surface, which is controlled by the intensity of the blood flow and the adhesion molecules stimulated by the excess of LDL, the penetration of the monocytes into the intima and subsequent inflammation of the tissue, and the rupture of the plaque accompanied with some degree of thrombus formation or even subsequent occlusive thrombosis. Their modelling approach predicts a priori the time to rupture as a function of arterial geometry, diameter of the monocyte, adhesion stress, bulk modulus of the ruptured wall material, blood viscosity, flow rate and mass density of the monocytes. Cobbold et al. [30] and Gessaghi et al. [31] studied the oxidation process of LDL cholesterol within the context of an in vitro framework. In addition, Cobbold et al. [30] considered the action of different vitamins such as vitamin E or C. Di Tomaso et al. [26] considered the interaction between just two species, LDL and monocytes, but monocyte behaviour was modelled in a very simple way. Fok [32] proposed a mathematical model of intimal thickening, posed as a free boundary problem. Intimal thickening was driven by damage to the endothelium, resulting in the release of cytokines and migration of SMCs. More complex studies were presented by Siogkas et al. [33], who included in their model oxidized LDL, macrophages and cytokines, considering that all the LDL molecules and the monocytes were oxidized and differentiated, respectively, at the instant in which these agents pass through the endothelium. A similar study was presented by Calvez et al. [34] from a mathematical point of view, but their study also included the foam cell. Ougrinovskaia et al. [25] explored the uptake of cholesterol by different scavenger receptors of macrophages during earlystage atherosclerosis using an ordinary differential equation (ODE) model. It was found that macrophage proliferation rather than an increased influx of LDL particles drives lesion instability. Finally, Bulelzai & Dubbeldam [35] presented a qualitative mathematical model consisting of a number of ODEs for the concentrations of the most relevant constituents of the atherosclerotic plaque: macrophages, monocytes, foam cell and oxidized LDL. Chung & Vafai [27] obtained the atheroma plaque properties based on microstructure information using a pore theorem and fibre matrix model, and taking into account the LDL transport. Additionally, another study by Chung & Vafai [28] described the effects of fluid–structure interaction and pulsation on the LDL transport.
With all these considerations at hand, this work presents a computational modelling based on reaction–convention–diffusion equations coupled with blood–wall mass transport of the main biological species which lead to atheroma plaque development. These biological agents and substances are LDL, oxidized LDL, monocytes, macrophages, foam cells, SMCs, cytokines and collagen. There are two important challenges for mass transport modelling within the arteries: providing a proper set of governing equations and the appropriate choice of boundary conditions [36–38]. Therefore, the dependence of fluid and mass transport through the endothelium on local blood flow characteristics has been modelled using a threepore model, because this is more accurate than the use of a single pathway approach [39–41]. This model is represented by the vesicular pathway, normal junctions and leaky junctions. The transport properties of the leaky junctions, i.e. hydraulic conductivity, diffusive permeability and the reflection coefficient, are calculated locally in the dependence of the fraction of leaky junctions. The fraction of leaky junctions is correlated with WSS, using the available experimental data as outlined earlier [41]. Furthermore, the constants which define the biological interaction between the different species have been obtained from experimental studies found in the literature.
2. Hypothesis of the model
The endothelium is a single layer of endothelial cells, which are elongated in the direction of the blood flow. Endothelial cells are interconnected through intercellular junctions. IEL is an impermeable elastic tissue with fenestral pores that lies between the intima and media layers. The processes that lead to atheroma plaque evolution take place mainly in the intima and media layers. Therefore, it should be noted that in the context of arterial mass transport and atheroma plaque formation, the arterial wall refers to the intima and the media layers, with the adventitia as the outer boundary layer (not included in the model).
This model focuses on the purely mechanical factors which promote endothelial dysfunction such as low shear stress [42]. Therefore, endothelium permeability has been linked to WSS distribution along the internal artery wall in this model. The arterial mass transport is coupled with both the bulk blood flow in the lumen and the transmural flow in the wall. Therefore, fluid dynamic models and solute dynamic models have been included. A standard assumption made in arterial flow modelling is that blood is a homogeneous, isothermal and incompressible fluid [43]. Despite the fact that blood is formed by a suspension of particles, it is also valid to treat it as a homogeneous fluid. This is because the diameter of the suspended particles, such as red blood cells, is much smaller than the internal diameter of the artery [43,44]. Moreover, the assumption that blood acts as a Newtonian fluid is also an established approach, owing to the predominantly high shear rates in medium to large arteries. Under the earliermentioned assumptions, the laminar flow of blood is governed by the Navier–Stokes and continuity equations.
Additionally, the arterial wall can be treated as a porous medium. As such, it is of great importance to characterize the porous media transport models used in describing biological phenomena. Most human tissues can be treated as porous media as they are composed of dispersed cells separated by connective voids where blood flows [45].
Concerning the solute dynamics part of the model, mass transfer across the arterial wall occurs via two mechanisms: convection associated with the gradient of pressures which promotes transmural flow and mass diffusion caused by concentration gradients. Molecular diffusion is driven by solute concentration gradients within the arterial walls. These gradients are caused by metabolic uptake and the production of proteins and cells within tissue. The substances transported from blood through the endothelial, intimal and media layers usually encounter some mass transfer resistance. The resistance can be quite high, depending on the size and charge of the protein or cell [21].
The main processes taken into account in the model are the migration, proliferation, differentiation and apoptosis of the cells, secretion and degradation of the substances and the interconversion of species.
3. Governing equations
3.1. Fluid dynamics
3.1.1. Blood flow along the artery lumen
The pulsatile nature of blood flow has not been taken into account. Thus, the blood flow was assumed to be steady, incompressible, laminar, Newtonian and, hence, governed by the Navier–Stokes equation 3.1and the continuity equation 3.2
The subscript ‘l’ is used to represent the lumen, then, u_{l} is the velocity vector field in the lumen, p_{l} is the lumen pressure and μ_{b} and ρ_{b} the blood density and dynamic viscosity, respectively. The internal forces have been disregarded in this model. A given parabolic velocity profile is imposed at the inlet of the artery and a given pressure at the outlet 3.3where u_{0} is the mean inlet velocity as the average of the mean diastolic and systolic velocities in the left anterior descending coronary artery, in accordance with the experimental data given by Iliceto et al. [46]. R is the internal radius of the artery and r is the radius of the artery at the axial location. Finally, a noslip boundary condition was applied between the lumen and the endothelium interface.
3.1.2. Transmural plasma flow throughout the artery wall
The transmural velocity vector field u_{w} in the arterial wall is calculated with Darcy's law 3.4and the corresponding continuity equation 3.5where the subscript ‘w’ is used to represent the artery wall, κ_{w} is the Darcian permeability of the arterial wall, ρ_{p} and μ_{p} are the density and dynamic viscosity of the blood plasma, respectively, p_{w} is the pressure within the arterial wall which is determined by the blood flow simulation along the artery lumen, and ε_{p} is the porosity of the artery wall. In order to specify a pressure drop across the arterial wall, a constant pressure of 2333.022 Pa (17.5 mmHg) was assumed at the intima–media interface [41,47].
The coupling of the flow dynamics in the artery lumen with that in the arterial wall is achieved by a threepore model, and a transmural volumetric flux (J_{v} (m s^{−1})) through the endothelium is calculated by means of this theory, which takes into account the vesicular pathway, normal endothelial junctions and leaky junctions [41] 3.6where J_{v,lj}, J_{v,nj} and J_{v,v} are the volume flux through leaky and normal junctions and the vesicular pathway, respectively. Volume flux occurs mainly through normal and leaky junctions [48,49], therefore, J_{v,v} = 0. Moreover, the volume flux through the leaky and the normal junctions can be expressed as 3.7where L_{p,nj} is the hydraulic conductivity for the normal junction and L_{p,lj}, the hydraulic conductivity for the leaky junction defined as 3.8where represents the fraction of the surface area S occupied by the leaky junctions A_{p} with the endothelial cell radius, R_{cell}, w_{l}, the halfwidth of a leaky junction and ϕ_{lj} is defined as the ratio of the area of leaky cells to the area of all the cells [50,51]. L_{p,slj} is the hydraulic conductivity of a single leaky junction. Once ϕ_{lj} is known, the transport properties of the leaky junctions can be calculated accordingly. In the presented model, ϕ_{lj} is assigned locally at the endothelium as a function of WSS. The correlation between WSS and ϕ_{lj} is obtained using published experimental data 3.9where the unit area is considered to be 0.64 mm^{2} and LC is the correlation between the number of leaky cells per unit area of 0.64 mm^{2}, and the number of mitotic cells that is directly related to the WSS. For more details, see Olgac et al. [41].
Furthermore, the transport properties of the leaky junctions, i.e. hydraulic conductivity, diffusive permeability and the reflection coefficient, are estimated using the pore theory [41,52]. Accordingly, the hydraulic conductivity of a single leaky junction is given by 3.10where l_{lj} is the length of the leaky junction [50,51]. Finally, given all these considerations, the volume flux J_{v} through a single pathway can be estimated.
3.2. Solute dynamics
Two different steps have been taken into account. The first step corresponds to a steadystate diffusion and convection process along the lumen of the LDL particles and the monocyte cells. The second step is related to the transmural timedependent process of mass transport, where different substances and species are transferred between the lumen and the vessel wall and react between them. The main phenomena modelled using diffusion–convection–reaction equations are: monocyte transmigration and later differentiation into macrophages; LDL accumulation and further oxidation; ingestion of oxidized LDL by the macrophages and their apoptosis, leading to foam cells; differentiation of contractile into synthetic SMCs, the migration of synthetic SMCs and their secretion of collagen.
3.2.1. Step 1: steadystate diffusion and convection along the lumen
3.2.1.1. Lowdensity lipoprotein particles
The LDL transport in the blood fluid domain along the artery lumen is governed by the convection–diffusion equation: 3.11where D_{LDL,l} is the LDL diffusion coefficient in the lumen blood flow and C_{LDL,l} is the LDL concentration. For the mass transport calculations at the lumen, a constant concentration boundary condition is used at the lumen inlet. A convective flux condition is used to ensure that LDL is convected out of the domain with flow . Finally, an outflow and no flux boundary conditions were applied at the lumen outlet and lumen artery wall interface, respectively.
3.2.1.2. Monocytes
Analogously, the monocyte diffusion and convection processes along the lumen are regulated by 3.12where D_{m,l} is the diffusion coefficient of the monocytes in the lumen, and C_{m,l} the monocyte concentration in the lumen. Because the normal range for the monocyte count is from 200 to 950 cells mm^{−3} [53], an intermediate value of 550 cells mm^{−3} has been chosen as the initial concentration of monocytes at the inlet of the artery. Then, as for the LDL substance, the monocyte concentration distribution at the lumen inlet was assumed to have a uniform profile. Again, an outflow and noflux boundary conditions were applied at the lumen outlet and lumen–artery wall interface, respectively; with the flux definition .
3.2.2. Step 2: timedependent diffusion and convection along the artery wall
3.2.2.1. Lowdensity lipoprotein particles
The LDL diffusion and convection through the artery wall can be modelled as 3.13where C_{LDL,w} and D_{LDL,w} are the LDL concentration in the artery wall and the LDL diffusive coefficient in the plasma, respectively. Apart from the diffusive and convective terms of equation (3.13), the term relating to the reaction part, d_{LDL}C_{LDL,w}, corresponds to the LDL oxidation formation, indicating the concentration of LDL per second which is oxidized. d_{LDL} is the degradation rate of LDL. Besides, insulation boundary conditions were applied at both axial ends of the artery wall, where the flux of LDL substance is defined as .
Moreover, at the media–adventitia interface, a constant LDL concentration value of 0.005 times the LDL concentration retained at the lumen intima interface was assigned as Meyer et al. [54] reported. The lumen–artery wall interface coupling of solute dynamics is achieved by a threepore model that takes into account the vesicular pathway, normal endothelial junctions and leaky junctions. Thus, the solute flux of LDL at the endothelium was defined as a boundary condition by the Kedem & Katchalsky [18] equations.
The total LDL solute flux, J_{s,LDL} (mol m^{−2} s^{−1}), through a single pathway is given by Olgac et al. [41] 3.14with c_{l} being the luminal LDL concentration at the endothelium. The LDL deposition coefficient, which indicates the percentage of LDL particles from the blood flow deposited at the endothelium, LDL_{dep}, was chosen as 1 × 10^{−}^{3} c_{l} [41,55]. Furthermore, because the entrance of LDL into the arterial wall cannot grow indefinitely, there is a natural saturation incorporated by the in the denominator of equation (3.14). This means that the LDL retention increases as the cholesterol level of the blood flow increases. However, it has been considered that this natural saturation occurs when the LDL levels in the blood are pathological, i.e. there are high levels of total cholesterol (270 mg m^{−3}) [56]. Finally, P_{app} is the apparent permeability coefficient given as [41,57] 3.15where P_{lj} is the diffusive permeability of the leaky junction pathway, σ_{f,lj} the solvent drag coefficient for the leaky junctions, P_{lj}, the diffusive permeability of the leaky junction pathway and Z_{lj} a fractional reduction factor in the solute concentration gradient at the pore entrance. These were estimated by Olgac et al. [41]. P_{v} is the permeability of the vesicular pathway. Applying the proportion between the leaky junction pathway and vesicular transport reported by Cancel et al. [58], P_{v} is calculated to be 1.92 × 10^{11} m s^{−1}. Following all these considerations, the total LDL solute flux, J_{s,LDL} (mol m^{−2} s^{−1}), can be estimated.
3.2.2.2. Monocytes
The diffusion and convection governed equations for the monocyte cells are 3.16and , where C_{m,w} and D_{m,w} are the concentration and diffusive coefficients of the monocytes, respectively. Owing to the size of the monocytes, the convection term has been disregarded [34,35]. Equation (3.16) has several terms. The first corresponds to the monocytes which differentiate into macrophages, d_{m} being the constant rate that modulates this differentiation. The second term represents the monocytes natural death. The monocytes undergo random motion, are chemotactically attracted to the chemicals in the presence of oxidized LDL, and grow up to a maximal value . Therefore, the last term on the righthand side of equation (3.16) corresponds to the chemotaxis of the monocytes due to the presence of oxidized LDL. The natural saturation at this maximum value has been included by means of a Gaussian function. This Gaussian function is set to 1, whereas there are no monocytes in the artery wall, and it decreases to 0. The Gaussian function tends to 0 when the monocyte concentration is bigger than , a threshold of monocyte concentration which is set to the same value as the monocyte blood concentration.
Moreover, regarding the boundary conditions, insulation conditions were applied at both axial ends of the wall and the media–adventitia interface, and a perpendicular monocyte inflow to the endothelium has been defined. The solute flow dependence on the WSS has been modelled based on the study by Bulelzai & Dubbeldam [35]. The monocytes solute flux per unit of surface, J_{s,m} (cells m^{−2} s^{−1}), is modelled as 3.17which is consistent with the fact that larger shear rates imply a smaller coefficient for the monocytes [10,59]. The value of WSS_{0} designates the WSS at which the growth rate of the monocyte concentration owing to the signalling response by the endothelium is reduced by a factor of two compared with the zero wall shear rate response by the endothelium, and it is taken to be 1 Pa [35]. Because the monocyte captation for the artery wall is due to the presence of oxidized LDL in the intima, the second term on the righthand side of equation (3.17) contemplates in the numerator. The entrance of monocytes into the arterial wall will be regulated by the oxidized LDL concentration, the monocyte flux being unnecessary when the oxidized LDL is completely eliminated. Moreover, the monocyte flux depends on the lumen concentration of this species, indicating that the flux increases as the monocyte blood concentration increases. Finally, the factor m_{r} is a constant which determines the rate at which monocytes enter the intima for small WSS and per mol m^{−3} of oxidized LDL. According to the experimental test performed by Steinberg et al. [22], 90% of the monocytes which were injected into a mutant rabbit disappeared from the blood flow at 72 h. Thus, each day 30% of the monocytes of the blood flow are recruited by the artery wall. Then, considering a unit surface of the model, i.e. where dh is the differential of the height (h), the m_{r} constant should satisfy that divided per unit of surface is equal to 0.3 m^{3} mol^{−1} day^{−1}. Thus, m_{r} has been taken to be 5.5 × 10^{−4} m^{3} mol^{−1} day^{−1}.
3.2.2.3. Oxidized lowdensity lipoprotein
For the oxidized LDL process, we came to the following equations 3.18and , where and are the concentration and diffusive coefficient of oxidized LDL, respectively. The value has been considered as similar to that used for the LDL substance by Prosi et al. [39]. d_{LDL}C_{LDL,w} modelled the oxidized LDL process from the LDL. Finally, the quantity of LDL uptake by the macrophages per second (mol m^{−3} s^{−1}) is defined by the second term on the righthand side of equation (3.18), where C_{M,w} is the concentration of macrophages and LDL_{oxr} is the rate of LDL uptake by a macrophage. Finally, insulation boundary conditions were applied at all the walls, i.e. at the lumen–intima and media–adventitia interface and at the axial ends of the wall.
3.2.2.4. Macrophages
The evolution of the macrophages is modelled by 3.19and where D_{M,w} is the diffusive coefficient of macrophages in the arterial wall taken equal to the monocyte diffusion coefficient. The first term on the righthand side of equation (3.19) corresponds to the monocytes that have differentiated into macrophages, and the second term models foam cell formation. Foam cell formation depends on the quantity of oxidized LDL ingested by the macrophages (), which has been multiplied and divided by the constants M_{r1} and M_{r2}. The former corresponds to the oxidized LDL concentration per second that a single macrophage should ingest to convert into a foam cell. The latter refers to the rate of foam cell formation per second depending on the macrophage concentration. Insulation boundary conditions were again applied at all the walls of the intima–media layers. As in the case of the monocyte behaviour, the convective term for the macrophages has been disregarded.
3.2.2.5. Cytokines
Cytokines are secreted by numerous cells and are a category of signalling molecules used extensively in intercellular communication. By means of the cytokines, macrophages express the inflammation process occurring in atheroma plaque formation. The behaviour of the cytokines can be modelled as 3.20where C_{c}_{,w} is the cytokine concentration and D_{c,w} the diffusion coefficient for the cytokines in the intima. However, diffusion and convection effects have been disregarded, because it was considered that cytokines are retained in the macrophage membrane [60]. The cytokine concentration is the result of a dynamic balance between degradation and production. Thus, d_{c} is the degradation rate of cytokines, and C_{r} the production ratio of cytokines owing to the presence of oxidized LDL and macrophages [33]. Moreover, because macrophages are responsible for cytokine secretion owing to the oxidized LDL, the concentrations of both species are included in the cytokine production term.
3.2.2.6. Foam cells
The production of foam cells is governed by 3.21and where C_{F,w} and D_{F,w} are the concentration and diffusive coefficients of foam cells in the arterial wall. The convection and diffusion of the foam cells is disregarded, because their dimensions are too large [34]. Consequently, the convection transmural velocity flow cannot drag these cells. The reaction term corresponds to macrophage apoptosis, whose constant has already been defined. Noflux boundary conditions have been defined at the contours of the artery wall volume.
3.2.2.7. Contractile smooth muscle cells
Initially, SMCs were assumed to be in a contractile phenotype. The presence of cytokines modulates their differentiation towards a synthetic phenotype by means of 3.22where is the contractile SMC concentration and S_{r} the differentiation rate of these cells owing to the presence of cytokines [61]. The contractile smooth cell differentiation is modelled as a reverse sigmoid function where at the maximum concentration of cytokines, , the differentiation is maximum. Because the mitosis and death of these cells are regulated by the body of the organism and both rates tend to be in equilibrium, the proliferation and apoptosis of these cells are not included in the equations. Again, the convention of these cells has been disregarded owing to their size, as with other cells included in this model. Moreover, these cells are considered quiescent, and therefore there is no diffusion process. An initial concentration of contractile SMCs of 29.28 × 10^{12} mol m^{−3} is defined as the initial condition [62], and insulation boundary conditions were applied.
3.2.2.8. Synthetic smooth muscle cells
The behaviour of differentiated SMCs can be modelled as 3.23where is the synthetic SMC concentration. As in the case of contractile SMCs, there is neither convection nor diffusion for this species. The convection is disregarded due to the large size of the SMCs, and the diffusion is disregarded because these cells migrate instead of being spread by diffusion. Furthermore, synthetic SMCs are induced into the intima by the chemoattractant (cytokines), being the migration rate of these cells [63].
3.2.2.9. Collagen
Finally, the last process corresponds to the collagen formation which is synthesized by the SMCs. The intima and media layers are formed by collagen, among many other substances. In addition, the contractile SMCs produce collagen which is degraded by the normal mechanism of the body. However, this model considers just the collagen secreted by the synthetic SMCs, because this is the only collagen that contributes to the development of the atheroma plaque. It can be described by 3.24where C_{G,w} is the concentration of the collagen. The collagen is secreted at the rate G_{r}, and degraded at the rate d_{G}. The natural degradation of the collagen as the age increases has not been taken into account. Moreover, insulation boundary conditions have been applied and diffusion and convection processes neglected.
3.3. Plaque formation
The mass balance for open systems can be written as 3.25where ρ_{o} is the total density of the tissue in the reference configuration [64], are the source/sinks and M_{i} the mass fluxes of the i arbitrary species. are related to migration, proliferation, differentiation and apoptosis of the cells and secretion and degradation of the substances. The concentrations of each species have the property being the total material density of the tissue as the sum over all i. The densities, , change as a result of mass transport and interconversion of species, implying that the total density in the reference configuration, ρ_{o}, changes with time.
As mass transport alters the reference density, assuming that these volume changes are isotropic, it leads to the following growth kinematics where means the original concentration of a species in the undeformed configuration [64] and I is the second order unit tensor. For small strain hypothesis and isotropic growth, we can write 3.26where v is the velocity of the material points.
Finally, as the volume of the atheroma plaque is mainly due to the foam cells, SMCs and collagen, the volume contribution of the rest of the species has been neglected. It has been assumed that monocyte and macrophage sizes (before ingesting oxidized LDL) as well as cytokine, LDL and oxidized LDL sizes are small in comparison with the sizes of foam cells and SMCs [34]. Thus, the isotropic growth of the atheroma plaque can be estimated as 3.27where ΔC_{S,w} is the variation of both contractile and synthetic SMCs with respect to the initial concentration of these species before atheroma plaque initiation, and are the volume of one foam cell and one SMC, respectively, which can be estimated knowing the radius of a foam cell R_{F} and a SMC R_{SMC}, and ρ_{G} the density of the collagen, taken as 1 g ml^{−1} [65]. The foam cell has been considered as spherical and the SMC shape as ellipsoidal [66].
For the sake of clarity, the values and physiological meaning of the parameters in this and the previous equations are given in table 1.
4. Numerical examples
In order to demonstrate the performance of the scheme of equations presented earlier, a computational model of a left descending coronary artery is presented. A commercial FE code, Comsol Multiphysics, v. 5.2 (COMSOL), was used to implement the model.
4.1. Geometrical model and mesh
The computational geometry corresponds to an idealized axisymmetric representation of a coronary artery. The general dimensions of the model are 1.85 and 0.1702 mm for the inner radius and arterial wall thickness, respectively, in accordance with the dimensions of a standard left descending coronary artery [79]. The length of the model is determined for the distance needed to obtain a completely developed flow, and this has been assessed as 60 mm. Furthermore, a stenosis rate of 55% is included in the geometry in order to reproduce a recirculation zone comparable with those that develop at arterial bifurcations. The artery wall is composed of the three porous media of the endothelium and one layer which corresponds to the intima and media layer together. The adventitia layer has not been modelled, because this layer does not play a role in the atherosclerosis process. Figure 1a shows the threedimensional coronary artery model used, and the axisymmetric section which has been revolved in order to obtain the threedimensional model. Moreover, as figure 1b shows, a fine mesh was created in the lumen and artery wall. The lumen and the arterial wall were meshed with 32 848 and 48 610 triangular elements, respectively. Furthermore, grid independence tests were performed successfully for both fluid flow and LDL transport calculations. Previous sensitivity analyses were performed to choose the definitive mesh. These mesh sensitivity analyses were mainly tested on the velocity field, by varying the number of grid cells. The computational domain is considered to be sufficient for this study when further mesh refinement can only result in less than 1% change in velocity, and an adequate convergence of the maximum WSS is obtained.
Finally, the main defined boundary conditions are outlined in figure 2.
4.2. Results
The equations were computed on the presented geometrical model simulating the plaque growth for three different scenarios; (i) normal blood cholesterol level (1.2 mg ml^{−1}) for a period of 10 years, (ii) high blood cholesterol level (2.7 mg ml^{−1}) for 10 years, and (iii) 5 years of high cholesterol followed by 5 years of a normal cholesterol level, simulating that the patient is on a diet and therefore that the blood cholesterol concentration decreases. The results obtained for these three cases are presented in the following.
4.2.1. Step 1: steadystate diffusion and convection along the lumen
The first step consists of a steadystate flow along the lumen which transports the LDL and monocyte substances. The velocity obtained along the artery lumen showing the flow streamlines is presented in figure 3a, where it can be observed that the flow is fully developed. The maximum velocity reached is 1.18 m s^{−1} and, as expected this occurs in the central part of the blood vessel after the stenosis. Furthermore, the streamlines of the velocity are shown in order to identify the recirculation zones downstream of the stenosis. These zones of reverse flow indicate the place where the initiation of the atheroma plaque would probably occur, because the WSS in these zones is low. In addition, figure 3b shows the pressure along the coronary artery where a pressure drop is observed owing to the stenosis and the applied boundary conditions.
However, the most important result of the step 1 computation is the WSS distribution along the coronary artery model, because this distribution triggers the initiation of the atherosclerosis lesion. Figure 4a shows the contour map of this variable, showing a detail of the area where the WSS is lower. The WSS along the lumen wall varies from 1 to 2 Pa, except in the zone of the stenosis, where a WSS of 8 Pa is reached. For the sake of clarity, the WSS along the zcoordinate has been plotted in figure 4b, showing that the minimum WSS occurs approximately at the model length of 22 mm.
4.2.2. Step 2: timedependent diffusion and convection along the artery wall lowdensity lipoprotein and oxidized lowdensity lipoprotein
Figures 5⇓⇓⇓–9 show the development of the concentration of each of the species over time for the third computed case, i.e. 5 years of a high cholesterol level followed by 5 years of a normal cholesterol level, as well as the contour map of each species at the end of the process. The curves of each of these figures are depicted for two points, at 25 and 50%, respectively, of the thickness of the wall and placed at the length where the WSS is minimum. These points have been chosen in order to show the different concentration of each species along the artery wall. Figure 5 shows the results for the LDL and oxidized LDL. Viewing figure 5a where the LDL results are shown, it can be seen that at the beginning the LDL concentration reaches 1.4 × 10^{−}^{3} and 8.6 × 10^{−4} mol m^{−3} for points 1 and 2, respectively. These values remain constant during the 5 years of high cholesterol level. Then, when the blood LDL level decreases, the LDL concentration in the artery wall also decreases until 6 × 10^{−}^{4} and 4 × 10^{−4} mol m^{−3}, respectively, again remaining constant over the other 5 years. The LDL solute flux is independent of time, varying only with the WSS fluctuations and the blood LDL concentration.
Moreover, the concentration contour map at the end of the process (10 years) shows that the maximum LDL concentration is located at the minimum WSS values, being higher near to the endothelium, because the LDL molecules penetrate through it. The LDL becomes oxidized as it is diffused and convected. Figure 5b shows the results for the oxidized LDL, which increases until 1.6 mol m^{−3} during the first half year, because there are few macrophages at the beginning of the process. However, as the macrophage concentration increases, the phagocyted oxidized LDL concentration also increases. Despite this increase in the macrophage concentration and the blood LDL concentration decreasing at the fifth year of computation, not all the LDL becomes oxidized. An equilibrium state of approximately 0.18 mol m^{−3} is reached. The oxidized LDL evolution over time is similar for both the studied points, showing that the concentration of this substance is constant in the artery wall cross section owing to diffusion and convection processes. This is also corroborated by the concentration contour map of this substance.
4.2.2.1. Monocytes and macrophages
Figure 6a,b shows the results for the monocytes and macrophages, respectively. The monocyte flux depends on the oxidized LDL concentration, because it is regulated by the concentration of LDL. Thus, the tendency of the monocyte evolution over time is similar to that obtained for the oxidized LDL, reaching an equilibrium of 2.5 × 10^{12} and 2 × 10^{12} cell m^{−3} for points 1 and 2, respectively, for the same computation time as that stabilized for the oxidized LDL. Furthermore, the monocyte flux and removal, either owing to natural death or macrophage differentiation, modulates the obtained concentration. The concentration contour map at 10 years of atheroma plaque growth shows that the monocyte distribution across the artery wall is not uniform, being higher close to the monocyte source, i.e. the endothelium, as occurs for the LDL substance. Macrophage evolution undergoes a rapid increase over the first days of the process and subsequently decreases slightly owing to the initiation of foam cell formation. Moreover, the slope of the curve decreases as the blood cholesterol level decreases. The macrophage concentration contour map shows that these cells are concentrated around the endothelium lesion, being uniform along the artery cross section.
4.2.2.2. Cytokines and foam cells
The results for the cytokines and foam cells are shown in figure 7a,b, respectively. The cytokines increase at the beginning of the process, because their concentration depends on the macrophage and oxidized LDL concentrations, which also increase at this time. Once the macrophage production (monocyte differentiation) and removal (foam cell formation) reach an equilibrium, the cytokines remain constant at 3 × 10^{9} mol m^{−3}, decreasing until 1.2 × 10^{9} mol m^{−3} when the blood LDL decreases. The foam cell concentration, figure 7b, increases over time. There is a variation in the slope at 5 years when the LDL blood concentration decreases. Furthermore, the concentration contour maps of both substances at the end of the process are uniform along the artery wall cross section, showing higher concentrations in those areas where the WSS is at a minimum.
4.2.2.3. Contractile and synthetic smooth muscle cells
Figure 8a,b shows the results for the contractile and synthetic SMCs, respectively. Both types of SMCs are inversely related by means of the cytokine concentration. Therefore, the contractile SMCs whose initial concentration was defined as the initial condition at 29.28 × 10^{12} mol m^{−3}, decrease over time increases in order to obtain synthetic SMCs which could migrate to the endothelium, reaching an equilibrium at around the fourth year. By contrast, synthetic SMCs increase with time until approximately 3 × 10^{13} cell m^{−3}. Moreover, the concentration contour maps of these species show the differentiation of contractile SMCs into synthetic SMCs in those areas with low WSS.
4.2.2.4. Collagen
Finally, the collagen variation over time is shown in figure 9. Because the synthetic SMCs which migrate to the endothelium are responsible for the collagen segregation, the collagen concentration follows the same tendency as that of the synthetic SMCs, increasing over time until the equilibrium value of 190 g m^{−3} is reached. The concentration contour map of this substance is uniform along the artery wall section, reaching its maximum values where the WSS is at a minimum.
4.2.2.5. Plaque formation
Figure 10 shows the volumetric growth for the three studied scenarios: (i) normal blood cholesterol level (1.2 mg ml^{−1}) over a period of 10 years, (ii) high blood cholesterol level (2.7 mg ml^{−1}) over a period of 10 years, and (iii) 5 years of a high cholesterol level followed by 5 years of a normal cholesterol level. The colour legend indicates the displacements of each point in m. Defining the degree of stenosis as the ratio between the area of atheroma plaque which has growth plus the original healthy arterial wall divided by the area of a normal artery wall, then the stenosis ratio after 10 years is approximately 45, 85 and 70%, respectively, for the three simulated cases, showing clearly the influence of the blood LDL concentration on atheroma plaque evolution. The total volume of the plaque at the final stage is 18.15, 38.42 and 25.95 mm^{3}, respectively, and the corresponding change in flow resistance increases respect to the initial stage by 22, 271 and 48%, respectively. It should be noted that high LDL cholesterol in the blood is a risk factor for atherosclerosis disease. However, even if the blood cholesterol level is acceptable, atheroma plaque formation occurs.
5. Discussion
A mathematical and computational model of atheromatous plaque emergence has been presented that allows the lesion to grow in particular areas in relation to the haemodynamics of the blood flow. Atheroma plaque growth is a complex process involving a high number of biological processes and species. In connection with the effects of blood flow on plaque rupture, the complex blood flow dynamics in the arterial system has been considered to play a key role in the initiation and development of plaques [80,81]. When lumen narrowing occurs, the upstream part of the plaque is exposed to a high WSS, whereas the downstream part is still exposed to low WSS and inverse flow. Thus, the change in lumen diameter owing to plaque growth will influence the pattern of blood flow, and therefore the WSS [82]. Even in an inflammatory environment, endothelial cells are capable of responding to signals induced by the WSS [83,84]. Therefore, the low shear stress region of the plaque will remain atheroprone, whereas the high shear stress region is thought to activate alternative signalling pathways leading to matrix degradation and eventually resulting in destabilization of the cap covering the plaque [82]. Weakening of the cap will eventually result in plaque rupture.
Although the model developed is relatively simple, it captures some of the main features of atherosclerosis lesion formation. The model is in good agreement with the clinical hypothesis that correlates atherosclerosis occurrence with low WSS, as the numerical example presented here illustrates. Moreover, the results obtained are reasonable in terms of the biological process of atheroma plaque growth. In particular, the LDL concentration profile throughout the arterial wall at the areas with low WSS and the experimental results of Meyer et al. [54] are in good agreement. However, any experimental studies reproduce experimentally the changes in concentration of the rest of the biological agents and species in order to validate the obtained results. In addition, the model realistically predicts the stenosis evolution owing to the appearance of low WSS, because a stenosis of approximately 45, 85 and 70% of the lumen area has been obtained after 10 years, for the three studied cases, respectively, which is clinically viable [85]. However, the results cannot be quantitatively validated, because the appropriate experiments that directly measured the concentration of each species at the human artery wall over time are not currently available.
This model can be considered as a step forward towards understanding atheroma plaque development. However, it requires further improvement and several limitations should be mentioned. Concerning the mathematical model, only the most important biological species have been taken account. Other species such as different types of cytokines (IL4, IL10, IL13 or TGF beta) or Tcells have been omitted. Substances such as free radicals, which oxidize the LDL, and the processes in which they are involved have also been disregarded. There are so many cells and substances which play a role in the development of atherosclerotic lesions that including all of them is practically impossible. Only the main processes have been included in the model, whereas other important processes in the development of the atheroma plaque such as the degradation of collagen with age or the apoptosis of the SMCs have been ignored. The dynamics of macrophage and oxidized LDL infiltration from the arterial wall into the lumen during lesion progression remains poorly understood, because no correlation between oxidized LDL or macrophages and low WSS zones as well as other washout mechanics has been found in the literature. Therefore, these processes were not included in the model. The complexity of the model should be increased by including more biological agents and substances as well as the processes that interconnect all the species in order to provide more precise information about the atherosclerotic disease process. The WSS has been considered as the main trigger of atherosclerosis initiation. Thus, the cyclic stretch effects of vessel compliance, curvature, pulsatile blood flow or cardiac motion have been disregarded. Furthermore, other variables such as the OSI, the WSSG or the WSSAG should also be included. The model parameters were obtained from a wide variety of experiments on many different human or animal models. However, few parameters have been estimated according to the order of magnitude of these parameters and making choices that gave biologically reasonable results. A sensitivity study of the parameters used in the growth model should be carried out in order to assess their influence on the results obtained. The transport properties are set as constant, but in fact they are very likely to change during plaque formation [27]. This model assumes that the substances can move from the lumen to the arterial wall, but not in the reverse direction. Atheroma plaque growth according to the volume of species variation from the healthy to the pathological artery is measured in the paper, and it has therefore been considered that a small quantity of these species at the initial moment is not likely to dramatically modify the results. Finally, concerning the axisymmetric computational example, the model should be applied to more complex geometries such as coronary or carotid arteries with bifurcations and to real patientspecific geometries. Furthermore, an analysis should be carried out in order to assess the influence of different factors such as the geometry, the curvature or the movement of the artery wall.
Despite these limitations, given the results obtained it can be concluded that the functional regulation of the endothelium by local haemodynamic shear stress provides a model for understanding the focal propensity of atherosclerosis in the setting of systemic factors and this may help to guide future therapeutic strategies.
Funding statement
Support from the Spanish Ministry of Research and Innovation through the research project DPI201020746C0301 and financial support to M. Cilla from the Diputación General de Aragón through grant B137/09 are gratefully acknowledged.
Footnotes

↵† Present address: Julius Wolff Institute, CharitéUniversitätsmedizin Berlin, Berlin, Germany.
 Received September 23, 2013.
 Accepted October 14, 2013.
 © 2013 The Author(s) Published by the Royal Society. All rights reserved.