## Abstract

The work herein represents a novel approach for the modeling of low-density lipoprotein (LDL) transport from the artery lumen into the arterial wall, taking into account the effects of local wall shear stress (WSS) on the endothelial cell layer and its pathways of volume and solute flux. We have simulated LDL transport in an axisymmetric representation of a stenosed coronary artery, where the endothelium is represented by a three-pore model that takes into account the contributions of the vesicular pathway, normal junctions, and leaky junctions also employing the local WSS to yield the overall volume and solute flux. The fraction of leaky junctions is calculated as a function of the local WSS based on published experimental data and is used in conjunction with the pore theory to determine the transport properties of this pathway. We have found elevated levels of solute flux at low shear stress regions because of the presence of a larger number of leaky junctions compared with high shear stress regions. Accordingly, we were able to observe high LDL concentrations in the arterial wall in these low shear stress regions despite increased filtration velocity, indicating that the increase in filtration velocity is not sufficient for the convective removal of LDL.

- low-density lipoprotein transport
- lipid accumulation
- atherosclerosis
- leaky junction

mass transport of large molecules such as low-density lipoprotein (LDL) has received considerable interest in recent years because of its significance in the early stages of atherosclerosis. Locally elevated concentrations of LDL in the arterial wall are considered to be the initiator of atherosclerotic plaque formation (20, 26). Various experimental (2, 19, 21, 37, 38) and computational (1, 13–15, 25, 31, 33, 40) studies have been conducted to comprehend the pathways and mechanisms governing LDL transport from the artery lumen into the arterial wall and within the arterial wall. There are two pathways of LDL transport through the endothelium: by vesicular transcytosis, which is regulated by receptors on the endothelial cells (30), and through leaky junctions, most of which are located at the sites of dying or replicating cells (18, 19).

Early experimental studies on LDL transport were performed in animals. The main method was systemic injection of LDL solution with subsequent harvesting and examination of the animal's arteries. Lin et al. (19) observed thoracic aortas of rats with this method and counted the number of leakage sites of LDL, as well as the portion of leakage sites associated with mitotic cells. Their results support the hypothesis that transiently opened leaky junctions provide pathways for LDL transport. Truskey et al. (38) used the same method on rabbits to quantify the apparent permeability of the endothelium to LDL. Tompkins et al. (37) measured LDL concentration at various sites in the circulatory system of monkeys. Meyer et al. (21) used a different method and excised aortic segments of rabbits to incubate them in vitro at various pressure levels with LDL solutions. With this method, they were able to assess the effect of intraluminal pressure on LDL transport. Recently, Cancel et al. (2) used endothelial cell cultures to study LDL transport under pressurized conditions and found that LDL transport through leaky junctions constitutes the major part of total LDL transport (90.9%), whereas only a small portion occurs via the vesicular pathway (9.1%). Because the occurrence of leaky junctions is governed by the adjacent endothelial cells, and because the behavior of the endothelial cells is in turn influenced by the shear stress acting on their surface, LDL transport in the arterial wall is in part governed by local blood flow characteristics that determine the level of wall shear stress. The effects of shear stress on endothelial cells have been investigated in various studies. Although there is no study that directly investigates the relationship between shear stress and the fraction of leaky junctions, there have been studies on the effect of shear stress on cell shape (4, 17, 27), proliferation rate (12, 41), and hydraulic conductivity and permeability (27, 28). Levesque et al. (17) observed a stenosed dog aorta and correlated the cell shape with levels of shear stress. Suciu (32), Chiu et al. (4), and Sakamoto et al. (27) used cell cultures to observe the cell alignment of endothelial cells in response to shear stress. These experimental findings are used in this study to obtain a single correlation between shear stress and cell shape. Chien (3) observed rabbit thoracic aortas and counted the number of mitotic cells in regions of specific cell shapes. Their findings are used in this study to correlate cell shape to the number of mitotic cells. Finally, the percentage of leaky cells associated with mitotic cells reported by Lin et al. (19) is used as a last step to obtain a correlation between shear stress and fraction of leaky junctions.

On the computational side, there have been several approaches to model LDL transport. In simple wall-free models, transport is solely calculated in the artery lumen, and flow characteristics are correlated with the solute distribution. Wada et al. (40) used a wall-free model to study the concentration polarization phenomenon, and Kaazempur-Mofrad and Ethier (13) calculated the mass transport in a realistic human right coronary artery. The main drawback of wall-free models is that they do not provide any information on the transmural flow and solute dynamics in the arterial wall. Lumen-wall models, on the other hand, account for both the arterial lumen and wall transport: Stangeby and Ethier (31) treated the wall as a homogeneous single-layer porous medium and solved the coupled luminal blood flow and transmural fluid flow using Brinkman's equations. Multilayer models have also been used to represent intima and media separately (1, 14, 25), which better represent the arterial wall, but are rather complex compared with the single-layer models.

In this study, we hypothesize that the dependence of fluid and LDL transport through the endothelium on local blood flow characteristics can be modeled more accurately using a three-pore model than with a single pathway approach (15, 25, 33, 34). To test this hypothesis, we use an axisymmetric representation of a stenosed coronary artery, where the arterial wall is treated as a single-layer porous medium. The volume flux and the solute flux across the interface between the fluid and the porous domains are governed by a three-pore model representing the vesicular pathway, normal junctions, and leaky junctions. The transport properties of the leaky junctions, i.e., hydraulic conductivity, diffusive permeability, and reflection coefficient, are calculated locally in dependence of the fraction of leaky junctions. The fraction of leaky junctions is correlated to wall shear stress using the available experimental data as outlined above.

## Glossary

WSS

*a*- Radius of LDL
*A*_{p}- Total pore area
- c
- LDL concentration
- C
_{o} - LDL concentration at the inlet of the artery
*D*- Diffusion coefficient
*d*- Artery diameter
*F*- Function for the restricted diffusivity in a pore
*J*_{v}- Volume flux
*J*_{s}- Solute flux
*k*- Local mass transfer coefficient
*K*_{lag}- LDL lag coefficient
*K*_{p}- Darcy's permeability in the arterial wall
*l*_{lj}- Length of a leaky junction
*L*_{p}- Hydraulic conductivity
*n*- Local mass flux
- p
- Pressure
- p
_{o} - Pressure at the outlet of the artery
*P*- Diffusive permeability
*P*_{app}- Apparent permeability
*Pe*- Péclet number
*r*- Radial location
*r*_{w}- LDL degradation rate
*R*- Resistance
*R*_{a}- Radius of the artery
*R*_{a,t}- Radius of the artery at the throat section of the stenosis
*R*_{cell}- Radius of single endothelial cell
*S*- Surface area
*Sc*- Schmidt number
*Sh*- Sherwood number
- SI
- Shape index
*t*_{a}- Thickness of the arterial wall
*u*- Velocity
*U*_{o}- Mean velocity at the inlet of the artery
*w*- Half-width of a leaky junction
- WSS
- Wall shear stress
*z*- Axial location
*z*_{b}- Base point
- Z
- Fractional reduction factor in solute concentration gradient at the pore entrance
- μ
- Blood viscosity
- μ
_{p} - Plasma viscosity
- ξ
- Half-distance between leaky junctions
- ρ
- Blood density
- ρ
_{p} - Plasma density
- σ
_{d} - Osmotic reflection coefficient
- σ
_{f} - Solvent-drag reflection coefficient
- ø
- Fraction of leaky junctions
- Φ
- Partition coefficient
*#MC*- Number of mitotic cells
*#LC*- Number of leaky cells

## Superscripts and Subscripts

WSS.

- end
- Endothelium
- adv
- Adventitia
- l
- Artery lumen
- w
- Arterial wall
- v
- Vesicular pathway
- nj
- Normal junctions
- lj
- Leaky junctions
- slj
- Single leaky junction
*i*- Single pathway

## METHODS

#### Mathematical model.

A schematic of the mathematical model of LDL transport from the artery lumen into the arterial wall is shown in Fig. 1. A description of the model integration is given in the section *Model flowchart*. Although blood flow has a pulsatile nature, we assumed steady-state conditions in this study for simplicity. Although these will yield a different wall shear stress distribution than transient conditions, the fundamental properties of our model are independent of the type of flow. The current model could be modified to take into account transient conditions in cases where quantification of the results for transient conditions is desired. Under the assumption of steady-state conditions and incompressible blood with Newtonian rheology, the fluid dynamics in the artery lumen are governed by the Navier-Stokes and continuity equations: (1) (2) where ▿ and ▵ denote the gradient and the Laplacian operator, respectively, **u**_{l} is the velocity vector field in the lumen, *p*_{l} is the lumen pressure, and μ is the blood viscosity. The transmural velocity vector field **u**_{w} in the arterial wall is calculated with Darcy's Law, (3) (4) where κ_{p}, μ_{p}, and p_{w} are the permeability of the arterial wall, viscosity of the blood plasma, and pressure within the arterial wall, respectively. The viscosity and density are taken as μ = 0.0035 Pa·s and ρ = 1,050 kg/m^{3} for blood and as μ_{p} = 0.001 Pa·s and ρ_{p} = 1,000 kg/m^{3} for blood plasma (22).

The LDL transport in the fluid domain is governed by the convection-diffusion equation: (5) where c_{l} and *D*_{1} are the LDL concentration and diffusion coefficient in the lumen, respectively. The transport of LDL in the arterial wall is modeled similarly with an additional term that represents LDL degradation, i.e., (6) where c_{w} and *D*_{w} are the LDL concentration and diffusion coefficient in the arterial wall, respectively, *K*_{lag} is the solute lag coefficient, and *r*_{w} is the LDL degradation rate in the arterial wall. The diffusivity of LDL in the blood and in the wall are taken as *D*_{l} = 5.0 × 10^{−12} m/s and *D*_{w} = 8.0 × 10^{−13} m/s. *D*_{l} is calculated with the Stokes-Einstein Equation as described previously (6), and *D*_{w} is the value reported by Prosi et al. (25) and the references therein. The solute lag coefficient *K*_{lag} is taken as 0.1486, the value proposed by Sun et al. (33). Published LDL degradation rates are 1.4 − 6.05 × 10^{−14} s^{−1} (1, 25, 33, 34). We have chosen here an average value of *r*_{w} = 3.0 × 10^{−4} s^{−1}.

The coupling of flow and solute dynamics in the artery lumen with that in the arterial wall is achieved by a three-pore model that takes into account the vesicular pathway, normal endothelial junctions, and leaky junctions. Following Curry (6), the volume flux *J*_{v,i} through a single pathway is given by (7) where *L*_{p,i} and σ_{d,i} are the hydraulic conductivity and the osmotic reflection coefficient of the single pathway, respectively, and Δp and ΔΠ are the hydraulic and osmotic pressure differences across the endothelial layer, respectively. In the present study, the osmotic pressure difference is neglected to decouple the fluid dynamics from the solute dynamics. Therefore, *Eq. 7* can be written as (8) Again following Ref. 6, the solute flux *J*_{s,s} through a single pathway is given by (9) where c and c are the luminal and arterial wall side LDL concentrations at the endothelium, respectively, *P*_{i} and σ_{f,i} are the diffusive permeability and the solvent-drag reflection coefficient of the single pathway, respectively, and *Pe*_{i} is a modified Péclet number that is defined as the ratio of convective to diffusive fluxes through the pathway: (10)

The endothelium has a high resistance to LDL transport (23, 31, 35) and, as a result, c is much smaller compared with c. Neglecting c, *Eq. 9* can be arranged as (11) where *P*_{app,i} is the apparent permeability coefficient given as (12) and Z_{i} is a fractional reduction factor in solute concentration gradient at the pore entrance given as (13) In the absence of flow through the pore, i.e., when convection vanishes, *Pe*_{i} → 0 and Z_{i} → 1, which represents the pure diffusion case. As *Pe*_{i} → ∞, i.e., convection dominates over diffusion, Z_{i} → 0.

In our three-pore model, the total volume flux *J*_{v} and the total solute flux *J*_{s} can be obtained by summing the fluxes over each single pathway, i.e., (14) (15) where *J*_{v,v}, *J*_{v,nj}, and *J*_{v,lj} are the volume flux through vesicular pathway, normal junctions, and leaky junctions, respectively. *P*_{app} is the total apparent permeability of the endothelial layer given as the sum of apparent permeabilities of each single pathway: (16) where *P*_{app,nj} = *P*_{nj}Z_{nj} + *J*_{v,nj}(1 − σ_{f,nj}) and *P*_{app,lj} = *P*_{lj}Z_{lj} + *J*_{v,lj}(1 − σ_{f,lj}) are the apparent permeabilities of the normal junctions and the leaky junctions, respectively, and *P*_{v} is the permeability of the vesicular pathway.

Volume flux occurs through normal and leaky junctions (23, 35). Therefore, with *J*_{v,v} = 0, *Eq. 14* becomes (17) LDL flux, herein referred to as solute flux, occurs by vesicular transport and through leaky junctions, but not through normal junctions since they block the passage of solutes with a radius larger than 2 nm (LDL has a radius of 11 nm) (23, 35). Consequently, the diffusive permeability and the solvent-drag reflection coefficient become *P*_{nj} = 0 and σ_{f,nj} = 1, respectively, resulting in *P*_{app,nj} = 0. Accordingly, *Eq. 16* can be rewritten as (18)

The endothelium and the arterial wall exhibit a combined resistance against inflow from the artery lumen. We use an electrical analogy for the fluid dynamics as shown in Fig. 2. The flow, which is analogous to current in an electric circuit, is driven by the pressure difference, which is analogous to potential difference. The flow resistances for the normal and leaky junction pathways are derived from *Eq. 8* as *R*_{nj} = 1/*L*_{p,nj} and *R*_{ij} = 1/*L*_{p,lj}, respectively. Similarly, the flow resistance for the arterial wall is derived from Darcy's Law, *Eq. 3*, as *R*_{wall} = μ_{p}*t*_{a}/κ_{p}, where *t*_{a} is the thickness of the artery. Therefore, the filtration velocity can be expressed as (19) where *R*_{T} = *R*_{end} + *R*_{wall} is the total flow resistance formed by the endothelial resistance *R*_{end} and wall resistance *R*_{wall} connected in series, and, p and p_{adv} are the lumen side pressure over the endothelium and the pressure at the media-adventitia interface, respectively. The normal and leaky junctions are two parallel flow pathways with a combined flow resistance across the endothelium, *R*_{end}, given as (20) Accordingly, the volume fluxes through normal junctions and leaky junctions are (21) and (22)

#### Calculation of transport properties.

The endothelium is modeled as a boundary between the artery lumen and the arterial wall with endothelial clefts and leaky junctions in between the endothelial cells. Adopting the approach in Refs. 10, 42, and 43, the normal endothelial junctions are considered to be cylindrical pores on the junction strands at the endothelial clefts, whereas the leaky junctions are modeled as ring-shaped pores surrounding the leaky cells. The leaky junctions are assumed to be randomly distributed and spaced 2ξ apart, as shown in Fig. 3. A periodic unit can be constructed by drawing a circle of radius ξ centered at each leaky junction. The fraction of leaky junctions, ø, is defined as the ratio of the area of leaky cells to the area of all cells (10). Using this definition, (23) where *R*_{cell} is the endothelial cell radius taken to be 15 μm (10, 42, 43) (Fig. 3).

The transport properties of the leaky junctions, i.e., hydraulic conductivity, diffusive permeability, and reflection coefficient, are estimated using the pore theory (5, 6). Accordingly, the hydraulic conductivity of a single leaky junction is given by (24) where *w* and *l*_{lj} are the half-width and the length of the leaky junction taken to be 20 nm and 2 μm, respectively (10, 42, 43). The diffusive permeability of a single leaky junction is given by (25) where *D*_{lj} is the diffusion coefficient in a leaky junction, which is related to the free diffusion coefficient in the lumen, *D*_{l}, by (26) where α_{lj} = *a*/*w* is the ratio of the radius *a* of the transported molecule to the half-width of the leaky junction, and *F*(α_{lj}) is a function describing the restricted diffusivity in a pore. The solvent drag coefficient σ_{f,lj} for the leaky junctions is given by (27) The hydraulic conductivity and diffusive permeability of the leaky junction pathway are given by (28) (29) respectively, where Φ = 1 − α_{lj} is the partition coefficient, i.e., the ratio of the pore area available for solute transport to the total pore area. *A*_{p}/*S* represents the fraction of the surface area *S* occupied by the leaky junctions. Referring to Fig. 3, *A*_{p}/*S* is *A*_{lj}/πξ^{2} in the periodic unit circle of radius ξ, where *A*_{lj} = (2π*R*_{cell})(2*w*) is the area of a single leaky junction. Using *Eq. 23* for the fraction of leaky junctions ø, *A*_{p}/*S* is related to ø as (30) Therefore, once ø is known, the transport properties of the leaky junctions are calculated accordingly. In our model, ø is assigned locally at the endothelium as a function of WSS. The correlation between WSS and ø is obtained using published experimental data in several steps as described in the following paragraphs. To the best of our knowledge, no direct correlations between WSS and ø have been published to date.

Levesque et al. (17) studied stenosed dog aortas and correlated the cell shape with levels of shear stress. Cell shape is generally expressed in terms of a shape index (SI) defined as (31) where *area* and *perimeter* are the area and perimeter of a single cell. SI is one for a circle and approaches zero for a line. Suciu (32), Chiu et al. (4), and Sakamoto et al. (27) used cell cultures to observe the cell alignment of endothelial cells as a response to shear stress. As shown in Fig. 4, we have used the data from these studies to establish an empirical correlation between wall shear stress and endothelial cell shape index, SI: (32)

Chien (3) observed rabbit thoracic aortas and counted the number of mitotic cells in regions of specific shape indexes. As shown in Fig. 5, we have used Chien's data to correlate the shape index to the number of mitotic cells, *#MC*, per unit area of 0.64 mm^{2}: (33)

Lin et al. (19) observed thoracic aortas of 10 rats to determine the total number of mitotic cells and the total number of LDL leakage sites. He found that 45.3% of total LDL leakage sites are formed by mitotic cells and that 80.5% of total mitotic cells are leaky. In our model, we take the normalized axial location z*_{b} = z/*d* = 25 as our base point, where the WSS reaches a constant value after the disturbances introduced by the stenosis, and, at this base point, we calculate the number of mitotic cells using *Eqs. 32* and *33* with the local WSS value. Adopting the proportion by Lin et al. (19) that 80.5% of total mitotic cells are leaky, we calculate the number of leaky cells associated with mitotic cells by multiplying the number of mitotic cells with 0.805. Again, using the ratio by Lin et al. (19) that 45.3% of total LDL leakage sites are formed by mitotic cells, we calculate the total number of leaky junctions by dividing the number of leaky cells associated with mitotic cells by 0.453. Because we know the total number of leaky cells and the number of leaky cells associated with mitotic cells, we determine the number of leaky junctions associated with nonmitotic cells to be 0.307. Because there are no available experimental data on the change of the number of leaky junctions associated with nonmitotic cells with respect to WSS, we keep this portion constant in our correlation between WSS and the fraction of leaky junctions. With these assumptions, we obtain a correlation between the number of leaky cells, *#LC* per unit area of 0.64 mm^{2}, and the number of mitotic cells: (34)

Because the fraction of leaky junctions ø is defined as the ratio of the area of leaky cells to the area of all cells (10), next to the form given in *Eq. 23*, ø can also be written as (35) *Equations 32–35* constitute the correlation between wall shear stress and fraction of leaky junctions used in our simulations.

Among the three transport properties of normal junctions, only the hydraulic conductivity has to be determined, since the diffusive permeability *P*_{nj} = 0 and the solvent reflection coefficient σ_{f,nj} = 1 as described in the previous section. To set the hydraulic conductivity of the normal junctions, we use the electrical analogy as introduced in the previous section at the base point z*_{b} and set the filtration velocity to 1.78 × 10^{−8} m/s as reported by Meyer et al. (21). There are three unknowns in *Eq. 19* for the filtration velocity (*L*_{p,nj}, *L*_{p,lj}, and κ_{p}). *L*_{p,lj} is calculated with the pore theory as described earlier based on the fraction of leaky junctions at z*_{b}. Darcy's permeability of the arterial wall is estimated using the electrical analogy by employing the equation for filtration velocity between the wall side of the endothelium and adventitia, i.e., (36) where p is the wall side pressure at the endothelium, which is set using the pressure drop of p− p= 18 mmHg given by Tedgui and Lever (36), since p is known from fluid dynamics calculations. The remaining unknown in *Eq. 36* is Darcy's permeability of the arterial wall, and it is estimated to be *K*_{p} = 1.2 × 10^{−18} m^{2}, which agrees well with the literature values of 1.0 − 2.0 × 10^{−18} m^{2} (10, 25, 31, 34, 39). Once the permeability of the wall is fixed, the last remaining unknown, *L*_{p,nj}, is estimated with *Eq. 19* as *L*_{p,nj} = 1.16 × 10^{−9} m/(s·mmHg), which is within the range of estimated values of *L*_{p,nj} by Yuan et al. (43) and close to the value of *L*_{p,nj} = 1.58 × 10^{−9} m/(s·mmHg) by Tedgui and Lever (36) and *L*_{p,nj} = 1.19 × 10^{−9} m/(s·mmHg) by Vargas et al. (39).

The only transport property to be determined for the vesicular pathway is the vesicular permeability *P*_{v}. To find out *P*_{v}, *P*_{app,lj} in *Eq. 18* is calculated at z*_{b}. *P*_{lj} and σ_{f,lj} are calculated from the pore theory as described earlier, and *J*_{v,lj} is calculated from the electrical analogy for the fluid dynamics with *Eq. 22*. Accordingly, *P*_{app,lj} assumes a value of 1.9 × 10^{−10} m/s. Applying the proportion between the leaky junction pathway and vesicular transport reported by Cancel et al. (2), *P*_{v} is calculated to be 1.92 × 10^{−11} m/s.

#### Model flowchart.

A general flowchart of the model and its modules is shown in Fig. 6 and described in this section. The two modeled domains (artery lumen and arterial wall) are separated by an interface, the endothelium. The gray-colored modules in Fig. 6 labeled with lowercase letters and connected with solid lines are used over the entire artery, whereas the white-colored modules labeled with Roman numerals and connected with dashed lines are used at a single point z*_{b}. The governing Navier-Stokes, continuity, and convection-diffusion equations (Fig. 6*b*) in artery lumen are solved with proper boundary conditions and material properties (Fig. 6*a*) to obtain the blood flow velocity, pressure, and LDL concentration fields (Fig. 6*c*). WSS (Fig. 6*d*) is calculated from the velocity field, and a base point z*_{b} (I in Fig. 6) is chosen far downstream of the stenosis, where local WSS reaches a constant value after the disturbances introduced by the stenosis. The WSS value at this point is used with the WSS-SI and SI-*#MC* correlations and *#MC*-*#LC* proportion in the module (II in Fig. 6) to calculate the fraction of leaky junctions associated with mitotic and nonmitotic cells. The fraction associated with nonmitotic cells is kept constant with respect to WSS in the overall correlation between WSS and fraction of leaky junctions (Fig. 6*e*). Values of the fraction of leaky junctions (Fig. 6*f*) are assigned at the endothelium depending on the local WSS using the overall correlation between WSS and the fraction of leaky junctions (Fig. 6*e*). These local ø values are fed into the pore theory (Fig. 6*g*) to calculate the local transport properties of the leaky junction pathway. The calculated transport properties are fed into the three-pore model (Fig. 6*h*), and their values at the base point are used to determine Darcy's permeability of the wall, hydraulic conductivity of normal junctions, and vesicular permeability. The experimental values of Meyer et al. (21) for the filtration velocity and of Tedgui and Lever (36) for the pressure drop across the endothelium (III in Fig. 6) are used at the base point, and an electrical analogy for the fluid dynamics in the wall is adapted to calculate Darcy's permeability of the wall and the hydraulic conductivity of the normal junctions. In addition, the ratio between the leaky junction pathway and vesicular transcytosis for LDL transport reported by Cancel et al. (2) (IV in Fig. 6) is used to calculate the vesicular permeability. These parameters are fed into the three-pore model (Fig. 6*h*), which uses the luminal side pressure and concentration at the endothelium from the module (Fig. 6*c*) to calculate the solute flux *J*_{s} and, with the help of an electrical analogy, the volume flux *J*_{v} and its distribution among the normal and leaky junction pathways. Volume and solute flux are used as boundary conditions for the solution of the governing equations in the artery lumen (Fig. 6*b*) and the arterial wall (Fig. 6*j*). Darcy's law and the convection-diffusion equations (Fig. 6*j*) governing in the arterial wall are solved with proper boundary conditions and material properties (Fig. 6*i*) to obtain the transmural velocity, pressure, and LDL concentration fields in the artery wall (Fig. 6*k*).

#### Computational implementation.

The computational geometry is an axisymmetric representation of a stenosed coronary artery as shown in Fig. 7. The inner radius and wall thickness of the artery are chosen to be *R*_{a} = 1.85 mm and *t*_{a} = 0.34 mm, respectively, in accordance with the dimensions of the human left anterior descending (LAD) coronary artery (7, 9). A stenosis with 40 and 64% reduction of the diameter and area, respectively, is included in the geometry to produce a recirculation zone comparable to those that develop at arterial bifurcations. The stenosed region is 2*R*_{a} long and is placed 8*R*_{a} downstream of the inlet. The entire domain is 60*R*_{a} long, ensuring that fully developed flow is reached before the end of the domain.

The interface between the lumen and the arterial wall at the stenosed part (8*R*_{a} < z < 10*R*_{a}) is expressed as (37) where *r*(z) is the radius of the artery at the axial location z, *R*_{a,t} is the radius of the artery at the throat of the stenosis, and z_{1} = 8*R*_{a} and z_{2} = 10*R*_{a} are the start and end points of the stenosis, respectively.

For the flow calculations in the lumen, a parabolic velocity profile is used at the inlet: (38) where *U*_{0} is the mean inlet velocity taken to be 24 cm/s as the average of the mean diastolic and systolic velocities in the LAD (11). To match the experimental conditions of Meyer et al. (21), a constant pressure boundary condition of 70 mmHg is used at the outlet. The volume flux *J*_{v}, calculated with *Eq. 19*, is specified in the normal direction to the endothelium for both the artery lumen and the arterial wall.

For the calculation of transmural flow, a constant pressure boundary condition is used at the media-adventitia interface with the value of p_{adv} = 17.5 mmHg [Ai and Vafai (1) assumed p_{adv} = 30 mmHg for an intraluminal pressure of 100 mmHg]. At both axial ends of the wall, zero flux boundary conditions are applied in the axial direction.

For the mass transport calculations at the lumen, a constant concentration boundary condition of C_{o} = 3.12 mol/m^{3} is used at the inlet in accordance with the common human blood LDL concentration level (8). To ensure that LDL is convected out of the domain with flow, at the lumen outlet, a convective flux condition is used. Because the solute dynamics in the artery lumen are convection dominated, the diffusive flux is set to zero: (39)

The solute flux, *J*_{s}, calculated with *Eq. 15* is specified in the normal direction to the endothelium for both the artery lumen and the arterial wall. At the media-adventitia interface, a constant LDL concentration value of c_{adv}/C_{o} = 0.005 as reported by Meyer et al. (21) is used. At both axial ends of the wall, insulation boundary conditions are applied: (40)

A commercial finite element code, Comsol Multiphysics, Version 3.3 (COMSOL), was used to implement the model. The artery lumen was meshed with 29,100 quadrilateral elements and the arterial wall was meshed with 25,220 quadrilateral elements for the solution of fluid dynamics. Second-order elements for velocity and linear elements for pressure were used in the lumen. Linear elements were used for flow in the arterial wall. For the LDL transport calculations, the artery lumen was meshed with 86,400 quadrilateral elements and the arterial wall was meshed with 144,000 quadrilateral elements. Linear elements were used to solve the convection-diffusion equations in both domains. The flow calculations were performed first. Next, the obtained velocity field was used for the mass transport calculations. Grid independence tests were performed successfully for both fluid flow and LDL transport calculations. The maximum relative change in velocity, pressure, and LDL concentration values was under 0.1% compared with a coarser grid with one-half the number of elements.

## RESULTS

To ensure that the computational mesh used in the simulations is adequate, the wall shear stress distribution and Sherwood number (*Sh*) along the endothelium in the axial direction are compared with the corresponding analytical solution to the Graetz-Nusselt problem for a straight tube in Figs. 8 and 9, respectively. Leveque's approximate solution (16) to the Graetz-Nusselt problem is used, since the Schmidt number is large (*Sc* = μ/ρ*D* = 666,666). The Sherwood number is a dimensionless mass transfer coefficient defined as *Sh* = *kd*/*D*, where *k* = *n*/(C_{o} − c) is a local mass transfer coefficient and *n* is the local mass flux. The simulation results for the wall shear stress and Sherwood number agree very well with the analytical solution at the inlet section. Downstream of the stenosis, the simulation results again converge to the analytical solution after having deviated from them due to the disturbances induced by the stenosis. This is an indication that the computational mesh is sufficiently fine.

The computations were performed for both shear-dependent transport properties and constant transport properties. WSS is the same in both cases, since the changes in transmural flow are negligible compared with the flow in the artery lumen. The absolute value of WSS is 1.79 Pa at the base point z*_{b}. With the correlations given in the Section *Calculation of transport properties,* this WSS value gives us a fraction of leaky junctions of ø = 6.2 × 10^{−4}, which is acceptably close to the published value of 5.0 × 10^{−4} in Refs. 10 and 43. For the constant transport properties case, the fraction of leaky junctions at the base point, ø = 6.2 × 10^{−4}, is used to calculate the transport properties, which are then applied at the entire endothelium. In the shear-dependent transport properties case, all of the transport properties are calculated locally with respect to WSS. A jump in the absolute values of WSS and Sherwood number occur at the stenosis. WSS values go back down after a peak is reached at the normalized axial location z* = 4.434 and reach zero at the separation point, z*_{sep} = 7.510. A recirculation zone exists between the separation point and the reattachment point, z*_{reatt}, where WSS becomes zero again. The Sherwood number also approaches zero at the separation and reattachment points, since the convection of molecules decreases significantly at these points. The recirculation zone and the separation and reattachment points, as well as velocity streamlines, are shown in Fig. 10.

The volume flux over the endothelium in the vicinity of the stenosis is shown in Fig. 11. The increase in thickness of the wall at the stenosis causes an increase of flow resistance through the wall, *R*_{wall}, which leads to a decrease in the volume flux. At the reattachment point, there is an increase of the volume flux for the shear-dependent case, which is a result of higher hydraulic conductivity of the leaky junction pathway, *L*_{p,lj}, therefore of lower endothelial resistance to flow, *R*_{end}. The reason for this is the increased number of leaky junctions in this low WSS region. For the shear-dependent case, the volume flux can be decomposed into two parts: volume flux through normal junctions, *J*_{v,nj}, and volume flux through leaky junctions, *J*_{v,lj}. It is observed that the flux through leaky junctions is increased in low WSS areas, where the flux through normal junctions is decreased. The total filtration velocity assumes the value of 1.78 × 10^{−8} m/s for both cases at z*_{b}.

The solute flux across the endothelium in the vicinity of the stenosis is shown in Fig. 12. The total apparent permeability assumes the value of *P*_{app} = 2.09 × 10^{−10} m/s at the base point z*_{b}, which is very close to the value of 1.9 × 10^{−10} m/s reported by Truskey et al. (38) based on findings in rabbits. In both shear-dependent and independent cases, *J*_{s} decreases slightly at the stenosis because of the decrease in filtration velocity, which causes a decrease in the advective part of the solute flux through leaky junctions. In the shear-dependent case, *J*_{s} has four local maxima: at the two ends of the stenosis and at the separation and the reattachment points. These four points constitute the centers of the low shear stress regions of the model. The Péclet number for the leaky junction pathway assumes values between 130 and 220 at which the gradient of solute concentration at the pore entrance practically vanishes, meaning that there is no diffusion through leaky junctions and the molecules are simply convected (6).

The wall side concentration of LDL at the endothelium in the vicinity of the stenosis is shown in Fig. 13. The concentration is normalized by the inlet concentration C_{o}. We observe a similar pattern as seen with the solute flux: four peak concentration locations are present at the low shear stress regions. At the reattachment point, the wall side concentration reaches up to 26% of the luminal concentration. In the lumen, accumulation of LDL at the surface of the endothelium occurs because of the high resistance of endothelium to LDL transport, which is known as the concentration polarization phenomenon (35). However, the lumen side concentration at the endothelium only reaches a maximum value of 1% larger than that in the bulk flow.

In Fig. 14, a comparison between the wall concentration profiles obtained by our simulations at z*_{b} and those observed by Meyer et al. (21) are shown. The radial location is normalized with the thickness of the artery, which is *t*_{a} = 0.34 mm for our simulations and *t*_{a} = 0.17 mm for Meyer's results. The two profiles are in good agreement, with the exception of the measurement location of Meyer closest to the endothelium. The concentration value at the media-adventitia interface is fixed by the imposed constant concentration boundary condition. Its value is given by that measured by Meyer et al. (21).

## DISCUSSION

The transport of LDL from the artery lumen into the arterial wall has been simulated in an axisymmetric representation of a stenosed LAD coronary artery using a three-pore model to represent the endothelium. The effects of wall shear stress on the endothelial cells and the various pathways for both volume and solute flux through the endothelium are taken into account by the pore theory. Although the model geometry is quite tractable, it is sufficient to reproduce the key flow features in a coronary artery.

The electrical analogy of the endothelium and the wall serves well in determining the model parameters and explaining the velocity distribution through the normal and leaky junctions. As can be seen in Fig. 11, volume flux through normal junctions decreases at the spots where volume flux through leaky junctions increases. This is most pronounced around the separation and reattachment points, where WSS is low. The increase in volume flux through leaky junctions is the result of the increased number of leaky junctions and accordingly decreased flow resistance through the entire leaky junction pathway. Because normal junctions and leaky junctions are parallel pathways, the decreased flow resistance of the leaky junction pathway leads to a decrease in the volume flux through the normal junction pathway as can be explained with the electrical analogy. This decrease in volume flux through normal junctions was shown experimentally by Sill et al. (28). They performed experiments on endothelial cell cultures investigating the effect of shear stress on the hydraulic conductivity of the cell layer. In these experiments, cells grew under no shear conditions. When the cell layer reached confluence, it was first exposed to a pressure gradient and, after stability in filtration velocity was reached, to flow for one to three hours. It was observed that the filtration velocity through the cell layer increased up to two- to fourfold with the flow. Sill et al. (28) hypothesized that a change in the number of leaky junctions may be responsible for the shear-induced solute permeability, whereas alterations of normal junctions may be responsible for the shear-induced hydraulic conductivity. Because the experiment duration of one to three hours was not long enough for the cells to respond in a proliferative manner, we believe that this increase in filtration velocity with WSS is only due to the alteration of the structure of the normal junctions. One simple explanation could be that, with shear stress, the cells are elongated and the total length of contact between the cells is increased, which can lead to a higher available normal junction area for the volume flux. According to the findings of Sill et al., one might expect low volume flux in low WSS regions. Our simulations indicate that their findings reflect decreased volume flux through normal junctions in these low WSS regions, whereas the total volume flux is observed to be increased.

The wall side concentration profile at the endothelium in Fig. 13 reveals four high concentration locations: the two end points of the stenosis, the separation point, and the reattachment point. High concentration values are reached especially at the separation and the reattachment points, which is in agreement with the finding of Smedby (29) in an angiographic study showing that atherosclerotic plaques grow in the downstream direction of the stenoses significantly more frequently than in the upstream direction. The highest wall side concentration at the endothelium is observed at the reattachment point and is 26% of the luminal concentration. Tompkins et al. (37) measured LDL concentrations in four squirrel monkeys at a total of 280 sites. They observed elevated concentrations at 15 sites, with the highest one being 9.8% of the luminal concentration at 10% relative depth from the endothelium. At the reattachment point, we observed 14% of the luminal concentration at this depth. The discrepancy between our simulation results and the experimental observations could be explained in three ways. First, Tompkins et al. waited only for 30 min after the systemic injection of labeled LDL into the monkeys to measure the wall side concentrations. This may not be enough time for a stable concentration profile to be established at the wall. Therefore, the in vivo values of LDL may be higher than what the experiments suggest. Second, we have neglected the effects of pulsatile flow, which expands and contracts the flow recirculation zone downstream of the stenosis. This leads to a cyclic motion of the reattachment point, which would soften the local LDL peak observed at this point. Additionally, the effects of oscillatory shear stress on the formation of leaky junctions remain to be investigated. Finally, the correlation between WSS and fraction of leaky junctions in our model is based on several experimental results on various animals and cell cultures, since no direct correlation has been published to date. Especially, the correlation between the shape index and the number of mitotic cells is only based on the observations of Chien (3). The sharp increase in the number of mitotic cells with respect to shape index reported therein leads to a sharp increase of the fraction of leaky junctions with respect to WSS. This might be the reason for overestimated values of concentration at the high permeability and low WSS locations. To investigate this, we performed a sensitivity analysis by scaling the right-hand side (RHS) of *Eq. 32*, i.e., the correlation between wall shear stress and shape index and the RHS of *Eq. 33*, i.e., the correlation between shape index and number of mitotic cells by ±10%. *Equation 34* was also updated accordingly as described in *Calculation of transport properties.* As can be seen in Fig. 15, the wall side concentration at the endothelium did not change qualitatively, but quantitatively. The maximum relative change was observed at the base point z*_{b} with 66% for +10% scaling of the RHS of both *Eqs. 32* and *33*, and 36% for −10% scaling of the RHS of both *Eqs. 32* and *33*. The relative change at the reattachment point was 56 and 46% for the +10 and −10% scaling, respectively. The sensitivity analysis shows that accurate experimental data on human endothelial cells or real artery segments are necessary to quantify the fraction of leaky junctions at different shear stress levels. The analysis also shows clearly that the qualitative behavior of our model is not sensitive to the fraction of leaky junctions as a function of shear stress.

The concentration profile throughout the arterial wall at z*_{b} in Fig. 14 and the results of Meyer et al. (21) at 70 mmHg are in good agreement with the exception of the experimental measurement point that is the closest to the endothelium. The lower concentration values close to the endothelium in our simulations may be because of the lack of differentiation between the intima and the media in our model. Physiologically, the layer between the intima and the media, i.e., the internal elastic lamina, also features a high resistance to the transport of macromolecules (24). Therefore, the intimal concentrations are higher than the medial concentrations. Because there is no internal elastic lamina in our model, the LDL close to the endothelium is more easily convected to the depths of the wall, which may be the reason for the discrepancy between the simulation results and the experimental results. A multilayer model is not included in this study for the sake of simplicity, since the main focus of the work at hand is the effect of WSS on the endothelium and the volume and solute flux through it.

One previous study attempting to incorporate the effects of local WSS on LDL transport from the blood to the arterial wall was performed by Sun et al. (34), where the data from Sill et al. (28) were used to define local hydraulic conductivity with respect to shear stress. As mentioned before, this approach only reflects the effects of WSS on normal junctions, and the filtration velocity profile obtained by Sun et al. is similar to our results for volume flux through normal junctions in the sense that it decreases at low shear stress regions. However, in our model, to fully cover the effects of local WSS on volume and solute flux, different pathways for the fluxes are defined, and transport properties for each pathway are calculated using the pore theory. Therefore, the effects of WSS on leaky junctions are also included, and the total filtration velocity is found to be increased in low shear stress regions. Sun et al. (34) calculated high concentration values in the low shear stress regions where their filtration velocity is decreased and explained this as a result of weaker local convective clearance effects of the transmural flow. According to Sun's argument about convective clearance effects, one would expect low concentration values in our results at the reattachment point where the total filtration velocity *J*_{v} is increased. However, this is not the case because of the elevated values of solute flux resulting from the leaky junction pathway. The relatively low increase in total filtration velocity is not enough to clear away the molecules, and high concentration values are reached at these low shear stress regions.

In conclusion, the work at hand constitutes a novel way of modeling LDL transport from the artery lumen into the arterial wall that incorporates for the first time the effect of WSS on the endothelial cells and the pathways of volume and solute flux. Further development of the model could be achieved with the incorporation of a direct relationship between WSS and the fraction of leaky junctions that could be obtained from appropriate experiments not available currently.

## GRANTS

The financial support of the Swiss National Science Foundation through the National Center of Competence in Research in Computer-Aided and Image-Guided Medical Interventions is kindly acknowledged.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2008 by the American Physiological Society