The distribution of fluid forces on model arterial endothelium using computational fluid dynamics,

Numerical calculations are used in conjunction with linear perturbation theory to analyze the problem of laminar flow of an incompressible fluid over a wavy surface which approximates a monolayer of vascular endothelial cells. These calculations model flow conditions in an artery very near the vessel wall at any instant in time, providing a description of the velocity field with detail that would be difficult to identify experimentally. The surface pressure and shear stress distributions are qualitatively similar for linear theory and numerical computations. However, the results diverge as the amplitude of surface undulation is increased. The shear stress gradient along the cell model surface is reduced for geometries which correspond to aligned endothelial cells (versus nonaligned geometries). I Introduction The vascular endothelium is the simple epithelium that lines the cardiovascular system. It consists of a cellular monolayer which rests on a complicated matrix of cells and intercellular material. Intact endothelium provides a selectively permeable barrier to the passage of macromolecules from the bloodspace to the extravascular space. Moreover, the vascular endothelium bears the shear stress imparted by blood flow. The structure and function of the monolayer is affected by these mechanical factors (see A detailed description of the flow is needed. An endothelial cell in vivo witnesses flow which passes in rhythmic waves. The cells are only 1-2 fim thick and 20-50 ]x,m in the circumferential and axial dimensions. They are primarily affected by detailed flow behavior very near the wall. In this region, the fluid velocity profile is nearly linear, although the magnitude of the average velocity may vary by a factor of four or more from systole to diastole, and from point to point on the microscopically rough surface. Because of the small size of endothelial cells, flow at any instant may be considered quasisteady near the wall, as described by the local linear shear flow. One of the most difficult problems in fluid mechanics is studying flow details near a rough surface. The disturbed wall region is typically buried inside a boundary layer or is difficult to access in experiments. Additional complications include: [Davics et al. PNAS X3:2114 tiny dimensions (typically microns), wall probe interference effects, and the intrinsic difficulty of accurately measuring wall shear stress. Previous studies demonstrate that numerical solutions of the Navier-Stokes equations yield accurate predictions of flow characteristics in such circumstances • there is increased ATP utilization in cells exposed to shear stress, presumably in part due to contracting stress fibers 9 the permeability of a monolayer with no previous exposure to shear stress is transiently and acutely increased in response to flow DNA synthesis is altered in cells by exposure to flow affecting gene expression, and apparently producing a different phenotype. Few (if any) other cells in the body experience shearing forces of similar magnitude on only one side. Thus, it is difficult to identify analogous cellular models for comparison. An unrelated phenotypic modulation has been observed in microvessel endothelial cells in response to chemical factors In order to understand how shear stress produces such profound effects in endothelial cells, the detailed distribution of forces on a cell and monolayer of cells must first be known. In early studies, investigators did not consider cell shape. The average shearing force imparted by bulk fluid flow was considered determinant: high shear stress caused direct desquamation, and low shear stress caused concentration polarization effects at the wall The wavy wall problem has been studied in two dimensions by others Recent studies include numerical treatments of low Re flows over objects in shear flows. In one study, finite elements are used to estimate forces acting on a thrombus II Computational Methods A model surface was chosen which represents the cell monolayer (see Re Vw = 0 (2) The unknowns are the vector velocity u and the pressure/?. The single parameter appearing is the Reynold's number, Re = vr\ 2 /v, where a is the shear rate in the linear shear flow far from the surface, and rj is the surface undulation amplitude (characteristic length scale of the problem). The kinematic viscosity of the fluid is v = /x/p, where p is the fluid density and ix is the fluid viscosity. The equations are parabolic so velocity boundary conditions must be provided on all sides of the computational domain. The complete theoretical solution is included in the Appendix Our calculations simulate conditions in large arteries (such as the aorta) very near the vessel wall at any instant in time. The upper boundary can be represented as a shear flow at infinite distance from the surface. However, the computational code did not explicitly provide for boundary conditions at infinity. Instead we specified the upper boundary to be a rigid surface which is moved far enough away (at least 4 times the cell surface modulation amplitude) so that wall effects are no longer important. The velocity at the upper surface is fixed at the value corresponding to a linear shear flow field. The cell surface is the lower boundary of the domain. It is rigid and extends infinitely in x and z. Thus, the solution is determined by solving for one full period of the cell model surface in the relevant directions {x and z). Boundary conditions are expressed below: \u\ y^a> = ay (shear flow at large distance) \u\y= y =0 (zero velocity at wall) where a is the undisturbed shear rate far away from the wall. A shear rate of a = 800 s~' was specified for all calculations. Unsteady motion dynamics for physiologic frequencies are such that a quasi-steady approximation can be made (a = hsfcJv = 0.1 to 0.001). The computational code NEKTON was used for numerical solution of the problem 310/Vol. 114, AUGUST 1992 Transactions of the ASME surface as the lower boundary of the computational domain. NEKTON has powerful pre-and post-processing packages for mesh generation and visualization of results. The code runs on a wide variety of computers (from workstations to supercomputers). Thus, computational experiments can be performed on smaller machines, while production runs can be directed to the most efficient computers available The spectral element method for partial differential equations is the basis for spatial discretization. The method is summarized briefly in what follows. For an extensive description, one should see references 29 and 31. Spectral elements combine high-order (spectral) accuracy with the geometrical flexibility of low order finite-element methods. The computational domain is divided into K nondegenerate macro-quadrangles (spectral elements). In our problem, three-dimensional domains were broken up into "bricks," in which the two horizontal parallel faces are nondegenerate quadrangles The data, geometry, and solution, are approximated by high order polynomial expansions within each macro-element. A local Cartesian mesh is constructed within each element which corresponds to N x N x N tensor-product Gauss-Lobatto Legendre collocation points. The Gauss-Lobatto points are clustered near elemental boundaries; an arrangement which gives accurate approximation, and favorable interpolation and quadrature properties. Dependent variables are expanded in terms of (N -l) th order polynomial Lagrangian interpolants (through the Gauss-Lobatto Legendre collocation points) Spatially discrete equations are generated by inserting assumed forms of dependent variables into the governing equations, and requiring that the residual vanish in some integral and weighted sense. The computed numerical variables correspond to values occurring at the collocation points of the mesh. Convergence is obtained by increasing the number of macro-elements (K) or the order of the interpolants (TV) in the elements. The error decreases algebraically (like K~N) as K is increased; and exponentially for smooth solutions (like e~a N ) as N is increased Ill Results An analytical solution for linearized flow over a wavy wall is given in the Appendix Ty X is the normalized surface shear stress in the x-direction; and r yz is the normalized surface shear stress in the z-direction. The term jxa is the mean wall shear stress imposed by flow far (i.e., many times the cell height) from the endothelial surface. The solution predicts: 1. The surface shear stress in the x-direction consists of the sum of the average shear stress imposed by flow and a spatially varying stress perturbation due to cell shape. The magnitude of the shear stress perturbation depends on q and TJ/X X . AS T)/\ X (dimensionless surface amplitude) increases, the perturbation increases linearly. For q » 1 it is proportional to q. T yx is in phase with surface variations in x and z-it is maximum at the highest point on the cell surface, and minimum at the lowest point on the surface. 2. The presence of surface waviness introduces a lateral shear stress perturbation which is linear with rj/\ x . It is caused by the transverse flow away from surface peaks and toward surface valleys. As q becomes large (» 1), there is no dependence on q. T yz is 7r/2 out of phase with the surface waviness in the streamwise and transverse directions. It is maximum or minimum at points of maximum surface slope. 3. The pressure perturbation is linear with ^A*, but does not depend on q. It is asymmetric along the cell longitudinal axis, tending to increase the pressure on the proximal side and reduce it on the distal side. The pressure is 7r/2 out of phase with the surface variations in the direction of flow. Pressure is maximum or minimum at points of maximum slope in the cos(ax)cos(/3z) (5) fi ow direction. Numerical and analytical computations were compared for For a limited range, numerical results and linear theory predictions agree (not illustrated). Both numerical and theoretical methods predict that the wall shear stress r yx is maximum at the highest points of the coordinate surface, and minimum at the lowest points. The pressure distribution is 7r/2 out of phase in the direction of flow, and the wall shear stress and pressure distributions are periodic. Numerical magnitudes no longer agree with linear theory after the onset of separated flow. Maximum (r yXimax ) and minimum (r yx , m i n ) shear stress magnitudes for both numerical and analytical solutions are plotted in Figs. 5(a) and 6(a) for a range of parameter values. Groups of points corresponding to a particular geometry (fixed length/ width value) fall along the same line when shear stress and pressure are plotted vs. surface amplitude We did not resolve the exact amplitude where recirculation begins; however, the range which contains the critical amplitude is recorded in the table in •Transverse ribs **Vortices do not appear: streamwise ribs The analytical solution for surface pressure predicts a linear dependence on ?)/X x , and no dependence on q. The numerical result exhibits dependence on q\ namely, the slope increases with q IV Discussion The flow fields predicted by the numerical and analytical solutions are qualitatively similar. The wall shear stress and pressure distributions vary periodically at the wavy wall surface. However, the results from the two methods diverge as the amplitude of the surface waviness increases. ear theory predictions can be observed by comparing surface pressure distributions • a departure from linear growth of peak-to-peak pressure. • variation in the phase of pressure distribution. • contributions from higher harmonics of the pressure distribution. 9 dependence on the length/width ratio (parameter q). Others have obtained similar predictions At the highest surface points, the wall shear stress grows almost linearly with increasing surface amplitude as predicted by linear theory. Flow acceleration occurs along streamlines toward the peaks due to the constraint provided by the continuity equation. These processes are different than those producing flow separation in the lower surface regions. Nonuniform shear stress gradients of significant magnitude across a cell surface could be of potential importance for explaining flow induced morphological changes. The forces which result from a cell-to-cell variation on the order of the perturbation shear stress are sufficient to disturb protein-protein interactions. Bell [1] has determined that a noncovalent interaction is disrupted by a critical force of 10~5 dyne. The difference in shear force on 2 adjacent cells in laminar flow can be of order 10~4dyne, which corresponds to 10 protein-protein interactions. Experimental studies of others indicate that a force of ~ 1 dyne (10 5 protein-protein interactions) can detach a cell from a monolayer The predicted forces acting on the aligned monolayer are reduced in comparison to nonaligned endothelium. For a surface approximating a nonaligned monolayer, the perturbation shear stress can be as large as 34 percent of the average shear stress imposed by the primary flow. This decreases to 20 percent for aligned monolayers since the height/length ratio is reduced (essentially, the surface is less "bumpy"). Perhaps the monolayer is able to achieve stability by reconfiguring the actin filament system so that stress fibers attach to the apical membrane. Nonaligned cells do not have stress fibers in the proper arrangement to experience this stabilizing effect. Modeling the distribution of forces on cells also aids in understanding the role of shear stress in the pathophysiology of atherosclerosis. Endothelium exposed to large shear stress gradients displays dramatic changes in cell shape, density, and rate of division Acknowledgments We thank Prof. A. T. Patera of M.I.T. and Dr. Einar Ronquist of Nektonics for assistance with the computational program