A Mathematical Model on the Feedback between Wall Shear Stress and Intimal Hyperplasia

In this paper, we present a mathematical model linking blood flow, shear-dependent endothelium permeability and intimal thickening (hyperplasia) of the arterial wall, which is an initial stage in the development of atherosclerosis. The key concepts are that the intimal layer swells in response to the presence of excess oxidised LDL (OxLDL) in foam cells. The hyperplasia disturbs blood flow, affecting endothelial permeability via the wall shear stress (WSS). These changes produce a feedback mechanism. LDL is transported through the arterial wall by advection and diffusion, and the concentration of LDL at each time step is assumed to be quasi-steady since it equilibrates on a fast time scale. The process is controlled by the slow timescale of the increase in concentration of OxLDL. We consider a section of uniform axisymmetric artery, and impose an initial local injury or ‘hotspot’ of relatively high permeability that enhances the influx of LDL, triggering the development of a bump-shaped lesion. In the absence of further inflammatory processes, the lesion eventually decays back to the homeostatic state. The model is used to explore how the shape of the lesion changes over time, its effect on WSS, influx rates of LDL and the sensitivity of these processes to oxidation parameters. The lesion is shown to propagate downstream driven by regions of high and low WSS on either side of the bump, and it persists for some time after the hotspot has vanished, leaving ample time for further pro-atherogenic processes to develop.


Introduction
Atherosclerosis is a disease of the large and medium sized systemic arteries.It is a major cause of death and morbidity in Western societies.The disease is characterized by the development of lesions containing lipids, cell debris and immune cells focally within the arterial wall.Risk factors, such as elevated cholesterol, that predispose an individual to develop atherosclerosis are well known.Advances in physiology, imaging and experimental techniques have produced vast amounts of data.Despite these efforts, the detailed pathogenesis of the disease remains unclear.Endothelial injury or dysfunction [Gimbrone Jr. et al., 2000] is a widely accepted hypothesis for the initiation of atherosclerotic lesions.Pro-atherogenic alterations in endothelial function include the increase in adhesion of monocytes which, once attached to the wall, pass through the endothelium developing into macrophages.Responding to infection and foreign material, macrophages normally act as scavenging cells releasing cytokines and growth factors to promote subsequent healing.In the pathogenesis of atherosclerosis, this inflammatory response is abnormal and exaggerated.
The arterial wall comprises of three main layers, namely the intima, media and adventitia.The innermost layer is the intima and includes the endothelium, a monolayer of flattened cells lining the lumen.The endothelium offers selective permeability to water, electrolytes and other substances in the blood.It senses the wall shear stress (WSS) due to the flow of blood and mediates a complex immune response to the presence of altered LDL [Libby, 2006;Steinberg, 2005;Ross, 1999]; LDL oxidation has received particular interest as being a key reaction [Steinberg, 1997;Siennicka and Zapolska-Downar, 2003].In experiments, LDL oxidation is quantified using lag times and oxidation rates [Hendrickson et al., 2005].A more complex description of the oxidation process in vivo was given by Cobbold et al. [2002], who modeled intermediate steps in the oxidation process.The haemodynamics, especially localized regions of disturbed or oscillating flow, have been implicated as being a key factor in the disease influencing mass transport into the vascular wall [Caro et al., 1971;Cunningham and Gotlieb, 2005;Berger and Jou, 2000;Nerem, 1995;Resnick et al., 2003;Heise et al., 2003].Specifically regions of low WSS acting on the endothelium have been shown to correlate with sites of enhanced macromolecule permeability [Ogunrinade et al., 2002].The focal nature of the of LDL accumulation and the role of mass transfer is reviewed by Tarbell [2003].Experimental measurements have been made to measure these changes [Fry, 2002;Jo et al., 1991] using macromolecules such as albumin in attempts to correlate lesion location with flow features.Timmins et al. [2011] also confirmed that the pathobiologic differences (e.g., intimal thickening) are a direct result of the biomechanical environment (e.g., change of wall shear stress) based on quantitative micro-CT data on the pathobiologic response of vascular stents.
Several numerical and mathematical models of LDL transport through the arterial wall have been developed, and Prosi et al. [2005] identifies three main classifications for these models, each increasing in complexity and the number of parameters that need to be determined.Model parameters are often very difficult to obtain experimentally and are usually derived using pore or fibre-matrix models which themselves have limitations [Yang and Vafai, 2006].Imaging techniques such as magnetic resonance imaging (MRI) and ultrasound do provide some information about flow characteristics and geometry.However, it is not yet practical to map blood flow in all but the simplest cases [Oyre et al., 1997].Image-based computational fluid dynamics (CFD) [Steinman, 2002] methods are continually developing and provide details of blood flow in realistic arterial geometries.Most existing models assume transport parameters are fixed and independent of each other and their environment.This is clearly not the case.Parameters such as permeability or hydraulic conductivity are very sensitive to stimuli either under the action of mechanical forces or chemical agents.
In this study, we present a novel mathematical model describing the response of the wall in the form of a locally enhanced region of permeability or "hotspot", due to the temporary loss of endothelial cells caused by e.g., disease, a temporary change in wall shear stress, or an invasive procedure such as angioplasty or stent insertion.Our aim is to relate the transport of LDL through the arterial wall to intimal hyperplasia.In particular, we allow the wall shear stress to directly control the endothelial permeability.Wherever possible we assume the simplest functional form for each process depending on the information available.The advantage of our approach is that it is flexible enough to consider a range of interacting processes both in diseased and healthy states.Each mechanism can have its own functional form which can interact with other processes.As new mechanisms and interactions are identified the model can be easily refined making the best use of experimental observation and measurement.Since oxidation occurs on a much shorter time scale than the growth of intimal hyperplasia, a quasi-steady approximation is used to derive the concentrations of LDL and oxidised LDL.Olgac et al. [2008] developed a three-pore model for the transport of LDL through the endothelium, which was extended to a patient specific model [Olgac et al., 2009] and used by Kim and Giddens [2015] on carotid bifurcations.However, intimal hyperplasia was not explicitly modeled.Budu-Grajdeanu et al. [2008] have developed an axisymmetric, disc model for neointimal hyperplasia at a venous anastomosis, that describes the interaction of growth factors, smooth muscle cells and extracellular matrix, and predicts the decrease in the lumen radius over time.The production of the growth factors is simply assumed to be proportional to the radius of the lumen.Our tube model differs in that the swelling is local and depends on the axial coordinate, and the WSS-dependent permeability of the endothelium to LDL is explicitly modeled.
In Secs.2-5, we present the mathematical model including the link between permeability and WSS.As the intima layer of the wall is thin compared to the lumen radius, an asymptotic method is used to describe LDL concentrations within the wall.To simplify the analysis a perturbation solution [Smith, 1976] is used to calculate the wall shear stress due to the flow of blood over a small bump in the wall.In Sec. 6, an iterative numerical scheme is developed.The selection of parameter values is described in Sec. 7, followed by results in Sec. 8. Discussion and conclusions are given in the final two sections.Further details can be found in Goodman [2008].

The Mathematical Model
We consider a section of healthy artery of length L = 4 cm as shown in Fig. 1 and described in cylindrical polar coordinates (x, r, θ) .The endothelium is considered to be a membrane.The intima and media are assumed to be homogeneous continuum layers of thickness H I and H M respectively.Blood flow in the lumen is assumed to be lamina and a perturbation solution Smith [1976] is used to calculate the distribution of WSS along the endothelium surface at any instance of time.The flow is assumed to be fully developed with no entrance or exit effects.Based on existing transport models [Tarbell, 2003;Fry, 1985;Kedem and Katchalsky, 1958], the net flux J of LDL through the intima and media is where the subscripts I and M denote the initma and media layers.For each layer, D is the diffusion tensor, c the LDL concentration, u the radial component of the velocity of water moving through the wall, χ the slip coefficient for the transport of LDL, and r the outward radial unit vector.In the absence of data, we assume that the diffusion tensors are constant isotropic tensors and do not vary with position within each layer, i.e., The slip coefficient 0 < χ < 1 accounts for the resistance offered to LDL advection by the internal structure of the arterial layers.Osmotic effects can be neglected in the thicker walls of the large and medium sized arteries.We do not explicitly model the adventitia because LDL in this layer is assumed to be drained away by the vasa vasorum, so we apply an appropriate boundary condition at the media-adventitia interface.
LDL transport and oxidation within the intima is described by where c o is the concentration of oxidized LDL (OxLDL) produced and k is the oxidation rate.LDL is oxidized in situ at a rate kc; once OxLDL is produced it is assumed to be fixed in position and to decay at a background rate γ, which is included to prevent the accumulation of OxLDL at homeostasis.No oxidation occurs in the media and the transport of LDL satisfies In our model, hyperplasia is assumed to be the result of LDL oxidation and its accumulation within foam cells.The timescale of this process (days) is much longer than that associated LDL transport through the wall (seconds).Consequently, the concentration profile of LDL throughout the intima and media is taken to be quasisteady so that (3) and ( 5) become respectively.At the intima-media boundary, continuity of LDL concentration and flux require that c(x, r, θ, t) is continuous and At the media-adventitia boundary, we assume a uniform and fixed LDL concentration c(x, r, θ, t) = c adv at r = H adv . (9) In the steady state of healthy homeostasis, from (4), so that the distribution of LDL and OxLDL is defined by the geometry of the wall and the transport parameters k, γ, u, D and χ.In the absence of a hotspot and, since there is no LDL oxidation within the media, the concentration of LDL in the media is uniform, c = c adv .The value of c adv is found by solving for the concentrations of LDL and OxLDL in the intima, and then used in the boundary condition (9) for the concentration of LDL at the media-adventitia interface when the system is perturbed by the introduction of the hotspot.We also note that, in the mathematical model used by Prosi et al. [2005], the upper wall LDL concentrations profile was found to be insensitive to the type of boundary condition imposed at the adventitia, which supports our use of (9).The final boundary condition determines the flux of LDL into the arterial wall from the lumen, and we introduce a phenomenological relationship between the wall shear stress τ and the endothelial permeability κ, as follows.

Shear-dependent endothelial permeability
In the large arteries, the uptake of LDL into the wall is controlled by the endothelium and is not limited by fluid-phase transport in the lumen [Tarbell, 2003].The LDL flux J in from the lumen to the intima is mediated by κ and the concentration difference across the endothelium, together with the slow pressure-driven advective flux across the arterial wall (see ( 1)).The net radial LDL flux is given by where c b is the concentration of LDL in the blood.On combining (1) and ( 11), the final boundary condition specifies the gradient of LDL on the luminal surface as The endothelial permeability κ is composed of a shear-dependent permeability κ(τ ) and an imposed temporary local hotspot κ h , thus in which where a is the hotspot permeability coefficient and the parameters w, b and g in ( 14) define the shape and duration of the hotspot.Initially the hotspot increases, reaching a maximum at time t m .After t m , the hotspot reduces exponentially until its contribution to the total permeability is negligible.Endothelial permeability is known to decrease as the WSS increases from low to moderate values, and we model this by setting where κ min is the minimum value of permeability.The values of κ 0 and τ 0 are derived by requiring that a normal wall shear stress τ n produces a normal uniform permeability κ n in homoeostasis, and by specifying the sensitivity ∂κ/∂τ of the endothelial permeability with respect to WSS at the point (κ n , τ n ).The values of these parameters are given in Table 3.

Hyperplasia
In homoeostasis, in the absence of a hotspot, the initial intimal concentration c * of oxidized LDL is uniformly distributed along the length of the vessel.We assume that hyperplasia is proportional to the accumulation of OxLDL through the thickness of the intima, i.e., 1640011-6 Int.J. Appl.Mechanics Downloaded from www.worldscientific.comby UNIVERSITY OF GLASGOW on 12/21/16.For personal use only.
where h(x, t) is the increase in intima thickness and η determines how sensitive thickening is to the accumulation of excess OxLDL.Intima thinning is prevented by the Heaviside step function As the intima locally thickens, it protrudes into the lumen of the artery modifying the wall profile, as indicated in Fig. 1.In turn, the new wall profile alters the distribution of WSS, thus modifying the flux of LDL into the wall via the endothelial permeability κ(τ ).

Scaled equations
We scale the LDL transport equations ( 6) and ( 7), and define the dimensionless variables where H I and H M are the initial thickness of the intima and media, H 0 and H 1 are the radial coordinates of the lumen-endothelium and intima-media boundaries, and α is a typical value of the aspect ratio of the bump produced by intimal hyperplasia.
The flux into the intima is defined in (12); after scaling this boundary condition becomes We define Fourier transforms and inverses in the scaled axial coordinate X by 1640011-7 Int.J. Appl.Mechanics Downloaded from www.worldscientific.comby UNIVERSITY OF GLASGOW on 12/21/16.For personal use only.
etc.After scaling and applying the Fourier transforms, the intima and media equations reduce to respectively.Finally, we assume a separable solution of the form in the intima and media, respectively.After substituting ( 29) and ( 30) into ( 27) and ( 28) and performing the separation, we solve for ) and Θ M m (f ), as explained below.C(X, r I , θ) and C(X, r M , θ) are then found by applying the inverse transform (26).

Homeostasis -The Steady State Solution
For homeostasis, we seek steady state, axisymmetric solutions for the LDL concentration R I (f, r I ) in the intima, and R M (f, r M ) in the media.The ratio of the media thickness to lumen radius is Because the intimal and medial layers are thin compared to the lumen radius, ≈ 0.2, which justifies using an asymptotic expansion in powers of .We expand to order O( ) to resolve the effect of curvature of the vessel, so that Eqs. ( 27) and ( 28) become in which L j [•] is the linear operator, 1640011-8 Int.J. Appl.Mechanics Downloaded from www.worldscientific.comby UNIVERSITY OF GLASGOW on 12/21/16.For personal use only.
Table 1.Expansions for the LDL concentration terms R I and R M and for the coefficients A and F in each layer in Fourier space.

LDL concentration
The expansions of R I and R M , together with the terms A and F , for each layer are listed in Table 1.Note that the oxidation rate k enters the intima equation through an additional term in F , but the contribution from the separation constant does not enter until O( 2).The coefficient A in the intima and media is equal to the appropriate Péclet number Pe I or Pe M which indicate the relative importance of advection to diffusion in each layer, These are based on the fluid flow rates u 0 = u(H 0 ) and u 1 = u(H 1 ).
Using the expansions for R I and R M gives at O( 1) and at O( ) We now determine the general solution in each layer up to O( ).Within the intima, R 0 and R 1 take the form The coefficients A 10 , B 10 , A 11 and B 11 are given in the Appendix and determined by applying the boundary conditions, while λ 1 and λ 2 are the roots of the characteristic equation L[R 0 ] = 0. Similarly for the media, R 2 and R 3 take the form 1640011-9 Int.J. Appl.Mechanics Downloaded from www.worldscientific.comby UNIVERSITY OF GLASGOW on 12/21/16.For personal use only.
where the coefficients A 20 , B 20 , A 22 and B 22 are given in the Appendix and determined by the boundary conditions, while λ 3 and λ 4 are the roots of the characteristic equation L[R 2 ] = 0.

Boundary Conditions
As noted previously, in homeostasis, in the absence of a hotspot, the endothelial permeability is uniform, and the concentration of LDL in the media is uniform.The boundary conditions correct to O( ) are: from ( 12), continuity of the flux of LDL at r I = 0 from ( 8) 1 continuity of concentration of LDL at the intima-media boundary and from (8) 2 continuity of the flux of LDL at the intima-media boundary We note that in homeostasis, within the media, since the LDL concentration is uniform and R M simplifies to Applying these boundary conditions gives the steady state solution for homeostasis and determines the value of the concentration of LDL at the media-adventitia boundary, C adv , so that, when a hotspot is present, boundary condition (9) becomes The boundary conditions are used in deriving the coefficients listed in the Appendix.

Wall Shear Stress
To close the model, we calculate the wall shear stress (WSS) over the indentation using analytical results from Smith [1976].As the time scale for the evolution of the 1640011-10 Int.J. Appl.Mechanics Downloaded from www.worldscientific.comby UNIVERSITY OF GLASGOW on 12/21/16.For personal use only.
growth of the intima is much greater than the period of one heartbeat, we approximate the time-averaged WSS by assuming that the time-averaged flow is steady Poiseuille flow upstream of (i.e., proximal to) the region of intimal hyperplasia.We scale lengths on the undisturbed radius a of the blood vessels so that the axial and radial coordinates are x * = ax and r * = ar * , respectively (where asterisks denote unscaled quantities), pressure on ag where −g is the axial pressure gradient so that p * = agp, and fluid velocity on ga 2 /ρν where ρ is the density of blood and ν is its kinematic viscosity and thus u * = ga 2 u/ρν.The dimensionless Navier-Stokes equations for steady flow in a rigid tube are where Re = ga 3 /ρν 2 is the Reynolds number.The oncoming flow is assumed to be fully-developed steady Poiseuille flow with axial velocity which satisfies the no-slip boundary condition When Re 1 and for a small indentation of the form G m (x) cos(mθ), using a Fourier series in the circumferential (θ) coordinate, the flow over the bump consists of a boundary layer disturbance of height Re −1/3 and an unperturbed core flow.Writing the axial component of velocity as where Z = (1 − r)Re −1/3 is the boundary layer coordinate, and taking a Fourier transform the unscaled WSS is given by evaluated at the boundary.When 0 < h 1, Smith [1976] shows that at Z = 0 where 1640011-11 Int.J. Appl.Mechanics Downloaded from www.worldscientific.comby UNIVERSITY OF GLASGOW on 12/21/16.For personal use only.
Ai(ξ) is the Airy function, and G m (ω) is the Fourier transform of G m (x).Note that when m = 0 and ω = 0, T m is zero and so makes no contribution to WSS, and ( 55) is the form used in the numerical scheme described below.

Numerical Scheme
The structure of the numerical simulation is shown in Fig. 2. A fourth order Runge-Kutta method is used to propagate the solution forward in time using the oxidation process described by (4).At each timestep, a discrete fast Fourier transform (FFT) Press et al. [1988] is used to calculate a steady state distribution of LDL in response to a permeability boundary condition ∂ C/∂r I at r I = 0.The adventitia boundary condition ( 52) is discretely transformed as follows: and N is the number of points in the x-direction.The hotspot is kept constant throughout each time step.Both the permeability and boundary conditions are transformed into Fourier space before being solved.Afterwards an inverse FFT restores the solution to real space.Discrete Fourier series were used in the circumferential (θ) direction.
The permeability boundary condition ( 12) is made explicit in time by using the value of C at the previous time step when evaluating dC/dr I .As κ * and C are  generally functions of X and therefore f , dC/dr I is calculated in real space and then transformed to d C/dr I , avoiding the need for convolution.Table 2 lists typical values of the computational parameters.The numerical scheme was tested for time step convergence and dependence on mesh size N .We also tested that our numerical scheme was stable in all parameter regions.

Parameter Selection
Table 3 shows the range of parameters used in our initial simulation.The vessel diameter and wall structure are based on the dimensions of the common carotid Table 3. Parameter values.Where a range is given, the value in parentheses gives the nominal value used in the basic simulation.
In our model hyperplasia is assumed to be an instantaneous response to the presence of oxidized LDL in the intima.In vivo hyperplasia is in response to the uptake of oxidized LDL by macrophages and the formation of foam cells.This uptake by macrophages is a dynamic process and not included in the current model.We chose the oxidation and removal constants k and γ to give LDL profiles across the intima consistent with measurements for rabbits [Meyer et al., 1996].
Based on Solokis [1995], the retardation or slip coefficients of the intima and media are 1.0 and 0.6, respectively.The media is known to provide a much greater hydraulic resistance and this is reflected in the lower slip coefficient of the media.
The effective diffusivity for solutes in the fiber matrix of the intima D I = 5.0 × 10 −10 cm 2 s −1 is based on experimental data as used in Tada and Tarbell [2004] who investigated the effect of the internal elastic laminar and its effect on the distribution of macromolecules in the arterial wall.For the media D M = 1.0 × 10 −8 cm 2 s −1 .In the absence of more detailed data, we assume the diffusion tensors for the intima and media are diagonal and that the diffusion parameters are the same in the (x, r, θ) coordinates for each layer.
The thickening constant η used in ( 16) regulates the extent of hyperplasia in response to excess OxLDL in the intima.The analytical solution for wall shear stress [Smith, 1976] is valid only for small indentations of height O(aRe −1/3 ), where a is the lumen radius and Re is the Reynolds number, so η is chosen to produce thickening of this order.
For our purposes it is sufficient to know what Reynolds number would be physiologically representative.Techniques to measure the temporal and spatial distribution of shear stress in blood vessels [Cheng et al., 2002] are continually improving.Among the most comprehensive are measurements of the carotid and abdominal arteries using MRI or ultrasound.Using values calculated from Zhao et al. [2000] we know that wall shear stresses in the common carotid artery are approximately, peak ≈ 29.5±8.2dyn cm −2 , mean ≈ 12.1±3.1 dyn cm −2 .Using these values together with a typical lumen diameter gives 491 ≤ Re ≤ 1115 in the common carotid artery (CCA).

Initial steady state
Before the introduction of a hotspot, the uniform endothelium produces the concentrations of LDL and OxLDL in the intima shown in Figs.3-5  how the solutions depend on the intima thickness H I , the oxidation rate κ and the removal rate γ.In the initial steady state, c o is given by equation ( 10) and LDL in the media is uniformly distributed (not shown).
Figure 3 shows that thicker intima layers contain lower LDL concentrations that decrease form the endothelium boundary to the media boundary.It is interesting to compare the change in LDL flux distribution with depth for different intima thicknesses.For an intima less than about 10 µm thick the flux of LDL decreases almost linearly; in a deeper intima layer, the oxidation process reduces the spatial gradient 1640011-15 Int.J. Appl.Mechanics Downloaded from www.worldscientific.comby UNIVERSITY OF GLASGOW on 12/21/16.For personal use only.with increasing depth.Although concentrations are higher for thinner layers, thicker intima layers contain marginally greater quantities of OxLDL.
Lowering the oxidation rate (Fig. 4) changes the shape as well as the position of the profile of OxLDL within the intima.The differences in the associated OxLDL profiles are more pronounced.As the oxidation rate is reduced, OxLDL concentrations at the endothelium boundary decrease whilst increasing at the media boundary.At very low oxidation rates, it is possible for the LDL minimum to occur within the intima close to the media boundary.This is possible because the diffusive flux, driven by the concentration gradient, is less than the advective flux which is always radial.The net flux shown in remains positive in all cases.Furthermore, decreasing the OxLDL decay rate increases LDL concentrations more at the endothelium than media boundary.

1640011-16
Int. J. Appl.Mechanics Downloaded from www.worldscientific.comby UNIVERSITY OF GLASGOW on 12/21/16.For personal use only.Increasing the decay rate γ increases the OxLDL profile throughout the intima (Fig. 5), but has no effect on the radial flux or concentration of LDL within the layer.

Wall response to a hotspot
The response to the onset of a hotspot is shown in Fig. 6.The plots in the first column show profiles of the hotspot permeability κ h , wall thickening h, WSS τ and permeability κ(τ ) along the length of artery L at specific times after the hotspot is first applied.The flow of blood is from left to right.The second column shows the variation of the maximum and minimum values of these variables with time.Upstream, and far downstream, the luminal flow remains unperturbed, permeability and shear stress remain at their normal healthy levels.In the middle section, blood flow is perturbed by the emergence of a bump in the wall profile seen in Fig. 6(b).Initially this profile is a direct consequence of the Gaussian profile of the imposed hotspot.With time, the shape of the bump becomes increasingly asymmetrical, deformed and elongated downstream.This is because, as the hyperplasia grows, the influence of WSS controlled permeability becomes increasingly important in controlling the LDL influx, resulting oxidation and thickening process.The hotspot 1640011-18 Int.J. Appl.Mechanics Downloaded from www.worldscientific.comby UNIVERSITY OF GLASGOW on 12/21/16.For personal use only.
shown in Fig. 6(a) rises to a maximum by 3 days and then decays exponentially over 25 days.Because of the retention of OxLDL the wall profile, WSS and associated shear dependent permeability all take longer to rise and decay than the hotspot.With larger maximum hotspots, or with longer decay times, the minimum permissible permeability is reached.This further elongates the protrusion.Figure 6(c), the increase in WSS along the upstream face of the bump is larger than the decrease downstream.When t < 30 days, and the bump is pronounced, the WSS can fall below upstream values, reducing LDL influx.Later, because the bump becomes skewed, the WSS still rises but does not fall downstream as much.This explains why the minimum WSS occurs earlier than the maximum and why the maximum WSS persists for longer.The shear dependent permeability changes in response to the WSS.Within the hotspot region, permeability reduces towards a minimum value before increasing again after the maximum bump height is reached.Permeability mirrors this pattern in WSS but is prevented from falling below the minimum value of κ n /2 set in the model, as can be seen in Fig. 6(d).

LDL flux into the intima
Figure 7 describe how the total radial flux J and its three components, i.e., the fluxes due to (i) hotspot permeability J h , (ii) advection J a , and (iii) WSS controlled permeability J κ , vary with time and position.Initially the effect of hotspot dominates flux but later is replaced by the WSS component, advection makes a minimal contribution to the net flux except at much later times.After 6 days the hotspot has reduced and sufficient time elapsed for excess LDL to enter, oxidize and modify the wall shape.The asymmetric influence of shear permeability now becomes more visible, the flux downstream increases proportionately as time increases, and the hotspot flux decreases.As the bump becomes more prominent, the WSS flux can reduce below initial values as WSS rises over the bump.Similarly the lower downstream WSS region increases influx and contributes to the asymmetry.At 12 days the hotspot has reduced and the WSS flux becomes the most dominant component.By this time the overall thickening is less and as the bump shape is skewed there is no longer a high WSS region downstream.As time progresses further the wall fluxes return towards their initial values.Only in a small region, downsteam from the bump does the influx increase above its initial value.More prominent is the reduction of flux which corresponds to regions of elevated WSS.It is this shear controlled flux which is responsible for the reduced net flux visible at later times.through.LDL transport is mostly radial as the radial concentration gradient dominates.Accumulation is greatest toward the upper portion of the intima and focused around the hotspot region.Far upstream and downstream the profiles return to the uniform healthy levels show in Fig. 3.At 13.2 h the hotspot flux is dominant and the LDL profile is only slightly asymmetric.Later at 27.6 h the LDL contours are clearly asymmetric, being influenced by the WSS permeability.After 70.8 h, the skewed wall shape and reduced flux causes the LDL concentration to fall below normal levels.At later times this fall is recovered and the LDL profiles return to their original state.The OxLDL profiles are similarly asymmetric but do not fall below healthy levels.

Distribution of LDL and OxLDL in the intima
1640011-20 Int.J. Appl.Mechanics Downloaded from www.worldscientific.comby UNIVERSITY OF GLASGOW on 12/21/16.For personal use only.3. LDL, and consequently OxLDL, values initially rise above homoestatic values due to the presence of a hotspot, and then decrease as intimal hyperplasia changes the wall shear stress profile on the wall of the lumen.There is "rebound" in LDL values which later fall below homeostatic values at t = 70.8h before recovering at later times.

Dependence on timescales for hotspot
The influence of hotspot magnitude and decay is shown in Fig. 9. Decreasing the rate for the decay of hotspot increases the maximum wall thickening and extends the time taken for the lesion to reach its maximum height before returning more gradually to the initial homeostatic state.

Discussion
The healthy homeostatic distribution of LDL throughout the intima is consistent with that observed in experiments [Morris et al., 1991].The initial hotspot produces localized changes in wall profile that over time increasingly influence changes downsteam of the initial injury.The change in WSS caused by localised thickening extends far downstream and to some extent upstream as well.After the initial injury has decayed, there is a period of time over which the wall takes to recover its initial state.During this later period the wall profile is skewed downstream and the WSS no longer falls below healthy levels downstream.In this way the low WSS, enhanced permeability and thus region is reduced by having a feedback mechanism between WSS and permeability.Several changes occur to the structure and function of the arterial wall with age [Nicholas and O'Rourke, 2005] as well as WSS [Bond et al., 2011].Ageing arteries generally have a thicker intima, and we find that OxLDL levels are greater in thicker walls with correspondingly lower LDL levels.Increased retention of OxLDL in older, thicker arteries predisposes the arterial wall to further pro-atherogenic changes.
It is clear that the distribution of LDL within the intima is confined to a region close to the origin of the hotspot, toward the upper surface of the intima.In the case of relatively high oxidation and removal parameters, appreciable concentrations of LDL are soon removed close to the hotspot region.The spread of LDL through the intima's thickness depends on the magnitude of the endothelium's permeability as well as the oxidation and removal parameters.A higher hotspot permeability increases the flux of LDL into the intima, supplying larger amounts of LDL to the oxidation process, transporting LDL deeper into the intima.Wall thickening can be enhanced by either increasing the sensitivity of cells (by increasing η) or by increasing LDL flux (by increasing the size of hotspot).The difference between these is that the hotspot is localized where as η affects the whole wall.In the remainder of this section, we discuss the assumptions and limitations of our study.The uptake of LDL and its oxidation within the wall is a gradual process involving many biochemical messengers.Migration of circulating monocytes into the intima adds to the number of foam cells that can form to scavenge excess OxLDL.We have neglected the influence of chemokines in directing the motion of circulating monocytes.Macrophages developing into foam cells eventually degrade after becoming engrossed on OxLDL, this later contributes to a necrotic core in the later stage of lesion formation.
In vivo arteries bend and curve introducing "secondary flows" [Caro et al., 1992].The changes in WSS caused by these features is many of magnitude greater than changes caused by the slight changes in wall profile that we model.However our interest is restricted to a localized region where the hemodynamic conditions are relatively uniform.Smith's [1976] results only hold for indentations of height order H 0 Re −1/3 .The small lumen constrictions we consider are in the order of 6%, 0.4 mm on an initial lumen of radius 6.4 mm.Our model focuses on the very initial stages of intimal hyperplasia and would not be able to capture later stages where the composition, function and profile of the wall is known to change.In this study, we have assumed steady Newtonian blood flow as a reasonable approximation to the time averaged effects of the in vivo pulsatile flow.The influence of non-Newtonian blood viscosity has previously been shown to have little effect on the hemodynamics of large arteries [Ballyk et al., 1994;Steinman and Lee, 2007].
We prescribe a relationship between WSS and permeability.Equally we could have chosen to link other parameters such as water flow rate u, hindrance χ I or decay rate of hotspot to the WSS.By varying the relationship used we can identify features typical of different WSS-permeability relations.Other investigators have considered other wall structures such as the internal and external elastic laminae on macromolecules and other vasoactive solutes (e.g., nitric oxide).Additional parameters such as hydraulic conductivity or sieving values could be included in our model.Concentration polarization [Wang et al., 2003] makes it possible for a nominally normocholesterolemic subject to have locally hypercholesterolemic conditions.[Wada and Karino, 2002] used numerical simulations to investigate this mechanism, showing it to be dependent upon the radial flux rate and vessel geometry.In our model including this phenomenon would slightly increase the concentration difference driving LDL flux through the intimal layer.However, when compared to the uncertainty in permeability values this effect would be negligible.It has been suggested that retention of LDL within the extracellular matrix of the intima and media layers (by Proteoglycans) could be important in atherogenesis [Williams and Tabas, 1995] and not just increased flux rates.Oxidation and uptake by macrophages is only one of many processes that could contribute to this retention.We assumed uniform transport properties along the length of wall, only increasing the permeability.It is plausible to suggest that pre-existing metabolic differences such as transport properties could also affect retention.
Although experiments to identify the functions and responses of the arterial layers are challenging, there is still enough information to propose mechanisms that might be at work in the development of atherosclerosis.The method outlined in this paper can be used as a framework for future work.As new mechanisms are identified and the interaction between existing ones become elucidated the model could easily be extended.In the future, it is likely that advances both in computation and experimental methods will improve the understanding of the mechanisms involved.
In this paper, we have assumed axisymmetry but this restriction is readily removed.Goodman [2008] has shown that results for the growth and decay of peak values of intimal thickening, WSS and endothelial permeability are not significantly different to the antisymmetric case, so that the results reported here are robust to changes in the geometry of the hotspot.

Conclusion
Our understanding of atherogenesis has developed from the early work of Caro and Fry.It seems likely that in the future the number/complexity of interactions, both physical and biochemical, implicated in this disease process will only increase.The modeling challenge is how to best describe and couple these interactions in the light of new discoveries.Our model is one example of how fluid dynamics (WSS), wall permeability (the properties of the endothelium) and wall profile changes (hyperplasia) can be coupled in a computationally efficient scheme.+ (E 3 e λ4 − E 4 e λ3 )(λ 1 e λ2 − λ 2 e λ1 ).

Fig. 3 .
Fig. 3. Effect of intima depth H I on LDL, OxLDL concentrations and net LDL flux across the intima layer.

Fig. 4 .
Fig. 4. Effect of oxidation rate on steady state distributions of LDL and OxLDL.

Fig. 5 .
Fig.5.Effect of decay rate on steady state distribution of OxLDL.Here 1/γ = 3.1 days, and the oxidation rate is the same for all three curves.

Fig. 6 .
Fig. 6.Response to a hotspot: (a) hotspot permeability, (b) intimal thickening, (c) wall shear stress and (d) WSS controlled permeability, left along the length of artery at 3, 6, 9 and 12 days after the initial insult, and right variation of peak values with time.1640011-17 Int.J. Appl.Mechanics Downloaded from www.worldscientific.comby UNIVERSITY OF GLASGOW on 12/21/16.For personal use only.

Figure 8
Figure8shows the concentration of LDL and OxLDL within the intima layer at three different times.LDL concentrations vary most within the intima, where LDL is being oxidized.Once in the media LDL is no longer oxidized and LDL is flushed 1640011-19 Int.J. Appl.Mechanics Downloaded from www.worldscientific.comby UNIVERSITY OF GLASGOW on 12/21/16.For personal use only.

Fig. 8 .
Fig. 8. Contour lines of the concentrations of LDL (left) and OxLDL (right) within a crosssection of the intima.Parameter values are given in Table3.LDL, and consequently OxLDL, values initially rise above homoestatic values due to the presence of a hotspot, and then decrease as intimal hyperplasia changes the wall shear stress profile on the wall of the lumen.There is "rebound" in LDL values which later fall below homeostatic values at t = 70.8h before recovering at later times.