Abstract
We describe a threedimensional numerical simulation of interstitial flow through the medial layer of an artery accounting for the complex entrance condition associated with fenestral pores in the internal elastic lamina (IEL) to investigate the fluid mechanical environment around the smooth muscle cells (SMCs) right beneath the IEL. The IEL was modeled as an impermeable barrier to water flow except for the fenestral pores, which were assumed to be uniformly distributed over the IEL. The medial layer was modeled as a heterogeneous medium composed of a periodic array of cylindrical SMCs embedded in a continuous porous medium representing the interstitial proteoglycan and collagen matrix. Depending on the distance between the IEL bottom surface and the upstream end of the proximal layer of SMCs, the local shear stress on SMCs right beneath the fenestral pore could be more than 10 times higher than that on the cells far removed from the IEL under the conditions that the fenestral pore diameter and area fraction of pores were kept constant at 1.4 μm and 0.05, respectively. Thus these proximal SMCs may experience shear stress levels that are even higher than endothelial cells exposed to normal blood flow (order of 10 dyn/cm^{2}). Furthermore, entrance flow through fenestral pores alters considerably the interstitial flow field in the medial layer over a spatial length scale of the order of the fenestral pore diameter. Thus the spatial gradient of shear stress on the most superficial SMC is noticeably higher than computed for endothelial cell surfaces.
 computer simulation
 fenestral pore
it has been widely accepted that a confluent layer of endothelial cells forms a mechanosensitive barrier separating blood and the underlying tissue space and regulates many complex processes in the vasculature (6,19). In contrast to extensive studies describing the response of endothelial cells to fluid mechanical stimulation (shear stress), much less attention has been focused on the biochemical responses of smooth muscle cells (SMCs) to fluid flow forces. The level of shear stress on SMC, in the presence of an intact endothelium, was thought to be too low to influence SMC function due to very low interstitial flow velocities (order of 10^{−6} cm/s) in the medial layer. Recent studies, however, indicate that SMCs can be subjected to significant levels of shear stress associated with interstitial flow even though the superficial velocity of interstitial flow is very low (27).
Wang and Tarbell (27) developed a twodimensional (2D) model describing the medial layer of an artery as a periodic array of cylindrical SMCs residing in a uniform matrix composed of proteoglycan and collagen fibers and simulated interstitial flow in the media. Their results showed that under typical physiological conditions, the level of shear stress on SMCs could be of the order of 1 dyn/cm^{2}, a level that has been shown to affect SMC biology (1, 15,28). Tada and Tarbell (22) subsequently extended the 2D modeling of interstitial flow by incorporating the complex entrance condition induced by the funneling of flow through leaky fenestral pores in the internal elastic lamina (IEL). Their 2D results indicated that the average shear stress on the most superficial SMCs, immediately below the IEL, could be 10–60 times higher than on cells located far from the IEL.
Because of the growing awareness of the importance of fluid shear stress in modulating SMC function and the particular relevance of the most superficial layer in atherosclerosis and intimal hyperplasia (12, 25), the present study describes a full threedimensional (3D) simulation of the flow field around the SMC right beneath a fenestral pore. We find that previous 2D calculations provide useful upper bounds on shear stress levels, but full 3D calculations are required to determine the distribution of shear stress along the cell axis.
Glossary
 a
 Distance between IEL and upstream end of SMC
 D
 SMC diameter
 d
 Fenestral pore diameter
 F
 Volume fraction of SMC
 f
 Area fraction of fenestral pore
 K_{p}
 Permeability of media
 L
 Distance between the center of neighboring SMC
 l
 Fenestral pore spacing
 P
 Pressure
 Re
 Reynolds number
 r
 Radial coordinate
 U
 Superficial flow velocity
 U_{0}
 Flow velocity at fenestral pores
 u
 Velocity vector
 u_{x}
 xcomponent of velocity
 u_{r}
 rcomponent of velocity
 u_{θ}
 θcomponent of velocity
 x
 Positional vector
 x
 xaxis of Cartesian coordinate system
 y
 yaxis of Cartesian coordinate system
 θ
 Angular coordinate
 μ
 Viscosity
 ρ
 Fluid density
 τ
 Shear stress
 δ
 Boundary layer thickness
 κ
 Dimensionless permeability
 wall
 Value at the SMC surface
 *
 Dimensionless value
METHODS
Numerical analysis.
In the present study, under realistic physiological flow conditions, 3D numerical simulations of the interstitial flow in the tunica media were performed to estimate the magnitude and spatial distribution of shear stress on SMCs taking into account the influence of the IEL in distributing flow to the media.
As described in detail elsewhere (23), but not shown here, transmural volume (water) flow driven by the transvascular pressure gradient enters the intima after passing through interendothelial junctions and then spreads laterally in the subendothelial intima that is separated from the media by the IEL. The elastinrich IEL is relatively impermeable to water flow except for its fenestral pores, which are filled with extracellular matrix material that is permeable to water (9).
A schematic illustration of the medial layer right beneath the endothelium is shown in Fig. 1. The IEL separates the subendothelial intimal from the medial layer and provides a complex entrance flow condition through the fenestral pores. In the present 3D computational model, the medial layer consists of a periodic array of cylindrical SMCs embedded in a continuous, porous interstitial phase. Transmural flow is distributed into the medial layer from the intimal layer through the fenestral pores in the IEL. The fenestral pores are modeled as circular openings distributed in a square array over the IEL. The IEL is assumed to be an impermeable rigid wall except for its pore openings. Justification for this assumption has been discussed extensively in previous studies (9, 23). SMCs are treated as circular cylinders impermeable to fluid because of the low hydraulic conductance of the cell membrane relative to that of interstitium (27). Under the assumptions described above, flow in the interstitial phase filled with fiber matrix is treated as Newtonian fluid flow (viscosity μ) through a homogeneous porous medium having a Darcy permeability coefficient,K_{p} .
Mathematical formulation.
The governing equations for fluid flow are Brinkman's equation (4), a generalization of Darcy's law for flow in saturated porous medium, which allows for satisfaction of noslip boundary conditions on the surfaces of SMCs and the IEL
Physiological parameters.
All of the physiological parameters used in the present study were defined and justified in the earlier 2D study of Tada and Tarbell (22). The SMC diameter (D) was 4.0 μm based on data for the rabbit thoracic aorta. The volume fraction of SMC (F) was always set at 0.4. For the fiber matrix in the interstitial phase, the K_{p} was set at 1.432 × 10^{−14} cm^{2}, which is consistent with the value used by Wang and Tarbell (27) based on the data of Tedgui and Lever (24) for the rabbit aorta. The superficial velocity of the fully developed interstitial flow, U= 5.8 × 10^{−6} cm/s, was obtained by applying Darcy's law to physiological data (10). The Reynolds number (Eq. 6) was 3.4 × 10^{−7}. In addition to the constants mentioned above, the area fraction of fenestral pores (f), the diameter of fenestral pores (d), and the distance between the IEL and the upstream end of the most proximal SMC (a) are important parameters defining the flow field. According to recent experimental data (rat thoracic aorta at 0–150 mmHg pressure) presented by Huang et al. (8), f and d are set at 0.05 and 1.40 μm, respectively. A range of values for each of these parameters was considered in our earlier 2D study (22), but we are restricted here by the computational demands of 3D simulations. The parameter a was estimated to be 0.36 μm (rabbit thoracic aorta) (22). However, because the distance a is thought to be the most critical parameter affecting the distribution of flow between the IEL and the first layer of SMC, the value ofa was varied in the range from 0.2 to 0.8 μm to examine its influence on the flow field and associated shear stress on the SMCs in the immediate vicinity of the IEL. All of the constants and parameters described above are listed in Table1.
Numerical analysis.
The FIDAP software package (8.52, parallel computing version; FLUENT Lebanon, NH, http://www.fluent.com) was used for the numerical simulations. The set of governing equations was integrated using a direct Gaussian elimination method with a structured, boundaryfitted coordinate system. The boundaryfitted coordinate was generated by using FIGEN, a submodule of the FIDAP solver. The global matrix arising from the finite element method discretization was decomposed into smaller submatrices to save memory space. Velocity grid points were arranged on each nodal point of the coordinate system, whereas the pressure was solved at the centroid of every finite element hexahedral element. In every finite element, biquadratic interpolation functions were used to approximate the velocities. For the pressure variable, a piecewise constant approximation was applied. Numerical simulations were carried out on the parallel supercomputing machine, Silicon Graphics Origin2000, at the National Center for Supercomputing Applications (NCSA, Champaign, IL).
The geometric configuration of the present 3D flow model is shown in Fig. 2. The fenestral pores were assumed to be aligned with SMC centers thus generating the highest possible shear stress on SMC surfaces. SMCs were arranged in a square array configuration leading to a relatively simple computational domain. Therefore, symmetry boundary conditions could be applied on the lateral surfaces of the rectangular computational domain the width of which was equal to the halfspacing of the fenestral pores. At the downstream end, a gradientfree (∂u*/∂z* = 0) outlet boundary condition was applied. The effect of the position of the outlet boundary on the resulting flow field was examined by increasing the number of columns of SMC parallel to the IEL surface (e.g., moving the boundary A to boundary B as shown in Fig. 2). The distribution of the shear stress over the last column of SMCs was used as an indicator of the fully developed outlet flow condition. Results were not affected by the exit length when two or more layers of SMCs were present. The value of the shear stress distribution over the second SMC layer, calculated by applying the outlet boundary described as boundary A in Fig. 2, showed less than a 2% maximum difference from that of the fully developed interstitial flow field. Therefore, the boundary A in Fig. 2 (two columns of SMC present in the computational domain) was taken as the outlet boundary throughout the remaining numerical analysis.
The complete computational domain used in the present study is shown in Fig. 3. A uniform velocity profile was imposed on the interstitial flow at the fenestral pore inlet. On all surfaces, i.e., the bottom surface of IEL and the periphery of SMCs, a noslip boundary condition was applied. Because thin Brinkman boundary layers arise in the vicinity of noslip surfaces, a finer mesh of grid points taken perpendicular to surfaces was employed to ensure sufficiently high resolution within the boundary layer of the interstitial flow. In the present problem, the Brinkman boundary layer thickness (δ) is of the order
Numerical validation.
To validate the present interstitial flow simulation method, the wall shear stress on a SMC (τ) was computed under the physiological condition that the SMC volume fraction (F) for an “infinite” (without IEL) periodic square array of SMC was 0.4. The computed results were compared with the analytic solution of Wang and Tarbell (26), and the 2D numerical simulation result of Tada and Tarbell (22). In the present 3D calculations, perfect 2D profiles of the shear stress distribution (no variation of shear stress profile along the axes of the SMC) were obtained on the second SMC for each value of the distance a. The maximum value of the dimensionless shear stress on the second SMC obtained in the present 3D calculation was 2.59, whereas the analytically obtained value for fully developed interstitial flow (26) was 2.57, about a 0.7% difference. On the other hand, the average value of shear stress around the second SMC was 1.333, which is identical to that of the 2D simulation result (22) and only 0.3% greater than the analytic value (1.329) obtained by Wang and Tarbell (26). Both comparisons show favorable agreement and confirm the validity of the present numerical simulation procedure.
RESULTS
The mathematical model and computational methods described in the previous sections were applied in analyzing dynamics of flow through the medial layer near the IEL for different values of a, the distance between the IEL and the first layer of SMCs. The incoming flow rate at the fenestral pore (U
_{0}) used in the study was derived from the relation
Shear stress distributions on SMC.
Shear stress values are nondimensionalized with respect to the reference value introduced by Wang and Tarbell (26)
Figure 4 shows a perspective rendering of the shear stress distributions over a SMC for three different values of the distance between the IEL and the upstream end of the first SMC,a. The color scale represents the magnitude of the local shear stress (τ
Figure 5 shows the distribution of τ
The shear stress distribution transverse to the primary flow direction is of particular interest because this could not be estimated in our previous 2D calculations (22). For exactly the same interstitial flow conditions, the shear stress distribution along the leading edge (meridian at θ = π) of the first SMC is shown in Fig. 6. As expected, the shear stress gradient is steeper for smaller values of a, and the maximum shear stress is higher for smaller a. When Figs. 5 and6 are compared, it is clear that the fenestral pore has a similar effect on shear stress distribution over the first SMC surface in both the axial and circumferential direction, but the maximum shear stress in the direction of the SMC axis is slightly lower than in the circumferential direction
Figure 7 shows the variation of the average (Eq. 13 ) and maximum shear stress on the first SMC as the distance a changes. Elevation of the local shear stress on the first SMC due to concentrated flow through the fenestral pore does not alter the average shear stress very much because there are also regions of low shear stress (e.g., near stagnation points). The average shear stress, therefore, attains a low value (order of 1 dyn/cm^{2}) regardless of the distance a. The results in Fig. 7 demonstrate that the influence of the IEL on the shear stress distribution over the first SMC is extinguished at distances beyond about 1 μm from the IEL. For spacing greater thana ≅ 1 μm, the magnitude of the maximum and average shear stress are on the order of 1 dyn/cm^{2}, corresponding to fully developed interstitial flow.
Comparison with 2D simulations.
A 3D calculation requires far greater memory capacity and computation time than a 2D calculation, even for the idealized 3D geometry demonstrated in the present study. Therefore, comparisons between results obtained using the 2D model previously proposed by Tada and Tarbell (22) and those of the present 3D model are shown in Figs. 8 and 9 to assess the accuracy of 2D predictions. The physiological parameters U, d, f, andF applied in the 2D calculation by Tada and Tarbell (22) are 5.8 × 10^{−6} cm/s, 0.8 μm, 0.016, and 0.4, respectively. These parameters again correspond to a case where the fenestral pores are centered on the SMC axes. The same value of U _{0} (Eq. 8 ) was used in both the 2D and 3D calculations.
Figure 8 compares the 2D and 3D distribution of τ
The difference between the 2D and 3D predictions depends strongly on the distance between the IEL and the upstream end of the first SMC (a). Figure 9 shows a comparison of the maximum local shear stress as a function ofa for the 2D and 3D models. For both the 2D and the 3D calculations, the maximum shear stress decreases with increasinga, showing a similar trend. However, the discrepancy between the 2D and 3D predictions of τ
It should be pointed out that in another 2D study of flow distribution to a network of fenestral pores (23), we modeled the pores as rectangular slits having the same area as the circular pores. That model is equivalent to assuming the flow has been spread over the length of the SMC upon entrance to the media. This model results in lower shear stresses than the 3D predictions but will not be considered here.
DISCUSSION
The presence of IEL with fenestral pores can greatly alter the flow field around the SMCs bordering the subendothelial intima. Because the area fraction of fenestral pores is in the range 0.002–0.01, and the mean diameter of the fenestral pores ranges from 0.4 to 2.1 μm (8, 17), the velocity of water issuing from an individual pore could be more than 100fold greater than the superficial flow velocity in the medial layer, resulting in a significant change in the flow and shear stress distribution on the surface of SMCs lying in the immediate vicinity of the IEL. The elevated shear stresses on the cells in the most proximal layer may affect SMC biology and could play a role in SMC proliferation and migration that occur in atherosclerosis and intimal hyperplasia (20).
In this paper, particular emphasis was placed on determination of the shear stress distribution on the 3D SMC, which differs significantly from that of the 2D model reported previously by Tada and Tarbell (22). The detailed 3D simulation provides more definitive data on the distribution and magnitude of shear stress on the surface of SMCs residing right beneath the fenestral pore entrance.
Depending on the distance between the IEL bottom surface and the upstream end of this proximal layer of SMCs (a), local shear stress on the most proximal layer of cells could be 3–14 times higher than on the cells far removed from the IEL, which have been estimated to experience an average shear stress on the order of 1 dyn/cm^{2} (23). The 2D results in Ref.22 indicate even higher levels of shear stress elevation on the first SMC for fenestral pore diameters, and area fractions smaller than were considered in the present study (d = 1.4 μm; f = 0.05). Thus the first layer of SMCs may experience shear stress levels that are even higher than endothelial cells exposed to normal blood flow (order of 10 dyn/cm^{2}) (13). On the other hand, the simulations show that the effect of the IEL on shear stress levels is reduced sharply with increasing values of a. For a greater than ≈1 μm, both the shear stress level and the shear stress distribution on the first SMC approach the characteristics of SMCs far removed from the IEL. In contrast to the significantly elevated SMC shear stress in the focal area right beneath the fenestral pore, the shear stress averaged with respect to the whole cell surface remains on the order of 1 dyn/cm^{2}. This indicates the presence of steep spatial gradients of shear stress on the SMC surface due to the jet issuing from the fenestral pore. The calculated maximum shear stress gradient on the first layer of SMCs is ≈200 dyn · cm^{−2} · μm^{−1}, whereas for the distal layers the maximum gradient is only ≈1 dyn · cm^{−2} · μm^{−1}. It is interesting to note that the maximum shear stress gradient over the surface of endothelial cells in blood flow has been computed to be on the order of 5 dyn · cm^{−2} · μm^{−1}(2). Although there is no direct evidence available to indicate that shear stress gradients affect SMC biology, analogies to endothelial cells in which the cellscale variations in shear stress may activate mechanotransduction signaling path ways (2,6) suggest a potential role for shear stress gradients on SMC.
Recent in vitro studies have shown that vascular smooth muscle cells are responsive to shear stress in the range 1–25 dyn/cm^{2} and increase their synthesis of transforming growth factorβ, tissue plasminogen activator (25), heme oxygenase1 (26), nitric oxide (7, 15), prostaglandins (1, 28), and basic fibroblast growth factor (16). These results and our simulations suggest that in a normal artery with intact IEL, the innermost layer of SMC in the intimalmedial border may be the most active biochemically because of elevated shear stresses. In addition, several in vitro studies have shown that shear stress is an inhibitor of vascular SMC proliferation (21, 25) and migration (18). Other in vivo studies in ballooninjured rabbit carotid arteries (3) have demonstrated that increasing blood flow on the injured artery inhibits matrix metalloproteinase2 mRNA and intimal hyperplasia. Taken together, these studies suggest that high levels of shear stress on SMC are beneficial because they suppress SMC proliferation and migration, which might otherwise contribute to intimal hyperplasia.
Acknowledgments
This work was supported by National Heart, Lung, and Blood Institute Grant HL35549.
Footnotes

Address for reprint requests and other correspondence: J. M. Tarbell, 155 Fenske Laboratory, The Pennsylvania State Univ., Univ. Park, PA 168024400 (Email: jmt{at}psu.edu).

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.

10.1152/ajpheart.00751.2001
 Copyright © 2002 the American Physiological Society