We describe a three-dimensional 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/cm2). 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 two-dimensional (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/cm2, 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 three-dimensional (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.
- Distance between IEL and upstream end of SMC
- SMC diameter
- Fenestral pore diameter
- Volume fraction of SMC
- Area fraction of fenestral pore
- Permeability of media
- Distance between the center of neighboring SMC
- Fenestral pore spacing
- Reynolds number
- Radial coordinate
- Superficial flow velocity
- Flow velocity at fenestral pores
- Velocity vector
- x-component of velocity
- r-component of velocity
- θ-component of velocity
- Positional vector
- x-axis of Cartesian coordinate system
- y-axis of Cartesian coordinate system
- Angular coordinate
- Fluid density
- Shear stress
- Boundary layer thickness
- Dimensionless permeability
- Value at the SMC surface
- Dimensionless value
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 elastin-rich 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,Kp .
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 no-slip boundary conditions on the surfaces of SMCs and the IEL Equation 1and the equation of continuity Equation 2In the above, ∇ is the Nabla operator, Δ is the Laplacian, u is the local flow velocity vector, P is the pressure, μ is the viscosity of the fluid, andKp is the Darcy permeability of the fiber matrix. The term on the left-hand side of Eq.1 is the pressure gradient that drives the flow, the first term on the right-hand side of Eq.1 represents the viscous force that allows satisfaction of the no-slip condition, and the last term on the right-hand side, the Darcy-Forchheimer term, characterizes flow in the porous medium away from the solid boundaries. After defining dimensionless variables Equation 3 where D is the diameter of the SMC, U is the superficial velocity of the fully developed interstitial flow, ρ is the fluid density, and x is the position vector in Cartesian coordinates, Eqs. 1 and 2 are rewritten in dimensionless form as Equation 4 Equation 5Re is Reynolds number defined as Equation 6
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 Kp was set at 1.432 × 10−14 cm2, 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.
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, boundary-fitted coordinate system. The boundary-fitted coordinate was generated by using FI-GEN, 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 piece-wise constant approximation was applied. Numerical simulations were carried out on the parallel supercomputing machine, Silicon Graphics Origin-2000, 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 half-spacing of the fenestral pores. At the downstream end, a gradient-free (∂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 no-slip boundary condition was applied. Because thin Brinkman boundary layers arise in the vicinity of no-slip 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 , or in dimensionless form, δ* = O( ). The order of magnitude of δ* is estimated to be 1 × 10−3; hence, a minimum mesh size of 1 × 10−5 was taken in the proximity of both the SMC and IEL surfaces. After it was established that the numerical results were independent of mesh density, a computational mesh consisting of ∼48,000 finite elements was applied. As the convergence criteria for implicit Gaussian elimination iterations, a relative error of velocity Equation 7was applied to the value of the nth iteration.
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.
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 Equation 8where U is the superficial velocity of the interstitial flow far removed from the IEL, l is the distance between neighboring fenestral pores, and d is the fenestral pore diameter. Equation 8 ensures that the volumetric flow rate distal to the IEL is equal to that supplied from the fenestral pores (see Fig. 3). The velocity profile at the fenestral pore was assumed to be a uniform flow distribution. This assumption is reasonable because the flow velocity at the pore exit has a flat distribution except in the immediate vicinity of the pore edge within the Brinkman boundary layer (22). The distance between neighboring SMCs (L) was determined from the volume fraction of the SMC (F) and the SMC diameter (D) as Equation 9 L = 5.6 μm was determined by using the physiological values F = 0.4 and D = 4.0 μm. The fenestral pore spacing (l) was taken to be the same asL. The use of this assumption and a fenestral pore diameter of 1.4 μm leads to a calculated fenestral pore area fraction off = 0.049, which is physiological (8).
Shear stress distributions on SMC.
Shear stress values are nondimensionalized with respect to the reference value introduced by Wang and Tarbell (26) Equation 10where μ, U, and Kp are water viscosity, the superficial flow velocity, and hydraulic permeability of the media, respectively. The nondimensional local shear stress on the surface of the SMC is defined as Equation 11where τθ and τx are defined as Equation 12where the asterisk indicates that the value is nondimensional, and the subscript “Wall” represents the value taken on the SMC surface (r* = 1). uθ is the angular component of the velocity (see Fig. 3), andur is the radial component of the velocity when the origin of the polar coordinate system is taken at the center of the SMC. τ is calculated to first-order discretization accuracy using the numerically obtained velocity variables uy and uz . Moreover, the wall shear stress averaged over the entire SMC surface is defined as Equation 13where the integration is taken with respect to the unit length of cylindrical SMC surface (0 ≤ θ ≤ 2π, 0 ≤ x≤ 2.8 μm, r = 2 μm).
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 (τ ) over the cell surface. For each value of a, a focal elevation of shear stress occurs on the upstream portion of the first SMC, because this surface is exposed to high-velocity flow issuing from the fenestral pore (fenestral pores are indicated with thick solid lines drawn on the IEL). Because the point on the surface of the SMC right beneath the fenestral pore center is a stagnation point of the interstitial flow, the wall shear stress is zero at this point. Moving away from the stagnation point, wall shear stress increases steeply along the direction of flow but then decreases as fluid spreads into the media. The maximum value of the shear stress on the surface of the first SMC (above τ* = 18 for a = 0.2 μm) appears within the projection of the fenestral pore periphery onto the SMC surface, and this value decreases markedly with an increase in the distancea. Hence, the influence of the IEL on the shear stress distribution over the SMC extinguishes within a spatial length scale of the order of the fenestral pore diameter (≈1 μm in this case). The shear stress on the second SMC is much lower than on the first and exhibits the same magnitude and distribution as that on SMCs exposed to fully developed interstitial flow. Thus both the shear stress level and spatial distribution of shear stress on the first SMC are greatly altered due to the presence of the IEL.
Figure 5 shows the distribution of τ around the periphery of the first and second SMCs at x* = 0 (the fenestral pore center axis) for three different values of the distance between the bottom of IEL and upstream end of the first SMC. The circumferential coordinate θ is taken in the counterclockwise direction from the downstream end of the SMC as indicated in Fig. 5. The flow direction is also indicated with an arrow in the figure. To obtain a sense of the dimensional magnitude of τ , the normalization factor in Eq. 10 , μU/ takes on a value of 0.52 dyn/cm2 based on data for the rabbit thoracic aorta (24). For each value of a, the maximum shear stress occurs at θ ≅ 7π/8, on the upstream side of the first SMC. The maximum value predicted in the present study is ∼9.0 dyn/cm2 for a = 0.2 μm. This order of magnitude corresponds to shear stress on endothelial cell surfaces exposed to blood flow. On the downstream side, for θ ≤ 5π/8, the shear stress profiles are independent of a and identical to the distribution on the second SMC, indicating that the entrance effect associated with the fenestral pore has been dissipated
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/cm2) 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/cm2, 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 τ around the periphery of the first SMC at x* = 0 for the case a= 0.36 μm. The 2D and 3D shear stress profiles are of similar shape; both profiles display a maximum at θ ≅ 7π/8 and another local maximum at θ ≅ π/2 before monotonically decreasing to zero at θ ≅ 0. The level of the shear stress predicted in the 3D model, however, is lower than that of the 2D model. Because fluid is not allowed to flow out of the SMC cross-sectional plane in the 2D model, a higher level of fluid velocity and shear stress are induced over the SMC surface relative to the 3D model, which allows for the spreading of fluid along the SMC axis.
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 τ becomes larger with increasing a. For instance, whena = 0.2 μm, the 2D model predicts τ ≅ 55, which is about 1.5 times greater than the value predicted by the 3D model. Whena = 0.8 μm, the 2D model predicts τ ≅ 17, which is about fourfold greater than the 3D prediction. With these comparisons in mind, it should be possible to estimate the influence of fenestral pore diameter (d) and area fraction (f) on maximum shear stress values by referring to the 2D simulation results in (22).
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.
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 100-fold 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/cm2 (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/cm2) (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/cm2. 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 cell-scale 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/cm2 and increase their synthesis of transforming growth factor-β, tissue plasminogen activator (25), heme oxygenase-1 (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 intimal-medial 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 balloon-injured rabbit carotid arteries (3) have demonstrated that increasing blood flow on the injured artery inhibits matrix metalloproteinase-2 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.
This work was supported by National Heart, Lung, and Blood Institute Grant HL-35549.
Address for reprint requests and other correspondence: J. M. Tarbell, 155 Fenske Laboratory, The Pennsylvania State Univ., Univ. Park, PA 16802-4400 (E-mail:).
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 © 2002 the American Physiological Society