Evaluation of Endothelial Shear Stress and 3D Geometry as Factors Determining the Development of Atherosclerosis and Remodeling in Human Coronary Arteries in Vivo
Combining 3D Reconstruction from Angiography and IVUS (ANGUS) with Computational Fluid Dynamics
Abstract The predilection sites of atherosclerotic plaques implicate rheologic factors like shear stress underlying the genesis of atherosclerosis. Presently no technique is available that enables one to provide 3D shear stress data in human coronary arteries in vivo. In this study, we describe a novel technique that uses a recently developed 3D reconstruction technique to calculate shear stress on the endothelium with computational fluid dynamics. In addition, we calculated local wall thickness, the principal plane of curvature, and the location of plaque with reference to this plane, relating these results to shear stress in a human right coronary artery in vivo. Wall thickness and shear stress values for the entire vessel for three inflow-velocity values (10 cm/second, 20 cm/second, and 30 cm/second equivalents with the Reynolds numbers 114, 229, and 457) were as follows: 0.65±0.37 mm (n=1600) and 19.6±1.7 dyne/cm2; 46.1±8.1 dyne/cm2 and 80.1±16.8 dyne/cm2 (n=1600). Curvature was 25±9 (m−1), resulting in Dean numbers 20±8; 46±16, and 93±33. Selection of data at the inner curvature of the right coronary artery provided wall thickness values of 0.90±0.41 mm (n=100), and shear stress was 17±17, 38±44, and 77±54 dyne/cm2 (n=100), whereas wall thickness values at the outer curve were 0.37±0.17 mm (n=100) and shear stress values were 22±17, 60±44, and 107±79 dyne/cm2 (n=100). These findings could be reconciled by an inverse relationship between wall thickness and shear stress for each velocity level under study. For the first time for human vessels in vivo, evidence is presented that low shear stress promotes atherosclerosis. As the method is nondestructive, it allows repeated measurements in the same patient and will provide new insights in the progress of atherosclerosis.
- Received August 26, 1996.
- Accepted January 27, 1997.
Atherosclerotic plaques are located more frequently near bifurcations, near trifurcations, and in curvatures of vessels compared with other sites of the arterial tree,1 2 3 4 5 implying that local factors are necessary to explain this condition fully. In search of these factors, several animal and human postmortem studies have provided evidence for a relationship between local shear stress, imposed by the velocity of the blood on the endothelium, and the localization of these predilection sites of atherosclerosis.6 7 8 9 10 11 On a cellular level, other studies have revealed a relationship between shear stress and the production of nitric oxide (NO or endothelium-derived relaxing factor), endotheline-1, and multiple growth factors.12 13 14
Despite the abundance of evidence relating shear stress to the normal adaptation of blood vessels and the pathophysiology of atherosclerosis, no 3D techniques are presently available that enable one to evaluate echographic-derived parameters as plaque size and location in conjunction with shear stress in human coronary arteries in vivo. In an earlier study, we (R.K., J.J.W., J.A.F.O., J.C.H.S., P.J.F., P.W.S., C.J.S.) reported on the possibility of reconstructing human coronary arteries in three dimensions applying a combination of ANGiography and intravascular UltraSound (ANGUS abstracts15 16 17 ). The present study focuses on the possibility of calculating regional vessel wall thickness and curve in combination with shear stress. Because these data were obtained in human vessels in vivo, this new method enables us to evaluate rheologic factors and their relationship to the progression and regression of atherosclerosis and vessel remodeling after interventions.
A detailed description of ANGUS, the 3D geometric reconstruction method, has been presented elsewhere.15 16 17 Briefly, the IVUS catheter (CVIS 2.9F, Sunnyvale, Calif) was visualized with a biplane angiographic system (Siemens, Bicor, Erlangen, Germany). The two views obtained from the x-ray system were digitized, and a custom written algorithm, developed in MATLAB (The Mathworks, Natick, Mass), was used to reconstruct the catheter path in 3D space. The echocardiographic vessel cross-sections, obtained from a motorized IVUS pullback at a speed of 0.75 mm/second, were selected offline at end diastole to avoid motion artifacts. These selected frames were digitized and analyzed with a semiautomatic contour-detecting algorithm as described before.18 Output of the program consisted of lumen contours, signifying the blood-vessel interface, and the media contours, representing the media-adventitia interface.15 16 17 18 These contours were filtered and subsequently positioned perpendicularly onto the 3D catheter path while the rotational position of the ultrasound transducer was selected on the basis of a comparison of the reconstruction with the coronary angiogram.15 16 17
The 3D lumen geometry of the vessel obtained with ANGUS was applied to define a mesh that was subsequently used for the computational fluid dynamics (Fig 1⇓). Generation of the mesh occurred in two steps. First, a straight cylindrical unit tube was created, ie, a straight tube with a circular cross-section with diameter and length equal to 1 m. In this tube, a mesh was defined consisting of 3D cells (“bricks”). A single cross-section of the tube consisted of 32 cells; and in the axial direction, the total tube included 100 cross-sections (Fig 1⇓). Second, based on the 3D-vessel geometry obtained with ANGUS,15 16 17 a transformation matrix was generated that was applied to each point of the unit tube to modify the lumen mesh until it resembled the 3D-lumen geometry of the original coronary artery (Fig 1⇓). The resolution for each cell in the mesh was 0.136 mm2 in the cross-sectional plane (range, 0.216 mm2 to 0.051 mm2) and 0.75 mm in the axial direction.
Wall Thickness Algorithm
The lumen and media contours for each cross-section were used to determine local vessel wall thickness at locations corresponding to the lumen mesh points. With a mesh cross-section consisting of 32 cells, we had 16 cells on the lumen-wall interface. With two nodes per cell, this resulted in 32 points on the border of each lumen cross-section of which 16 were used for the present analysis. Hence, with 100 cross-sections, the present analysis resulted in 1600 wall thickness values. Local vessel wall thickness was defined in the cross-sectional plane as the vector, directed from lumen to media, whose intersection with the local lumen and outer media cross-section was the most perpendicular as possible. In other words, the summed cosine of the intersection angles should be minimal.
We derived from the 3D spatially reconstructed vessel geometry a 3D idealized vessel centerline that passed through the center of the area of each cross-section. Polynomial curve fitting along the vessel centerline was applied to each of the coordinates separately. In general, a fifth order polynomial was necessary to describe the centerline accurately. At the point of intersection of the centerline with each cross-sectional mesh plane, the local curvature of the centerline was determined. The radius vector of a circle fitting through three subsequent centerline points and pointing from the origin of the circle to the midst of these points was selected for this purpose (Fig 2⇓). The so-defined local 3D radius vector was used (1) to calculate curvature (κ=radius−1) and (2) to select wall thickness at the location of the inner and the outer curvature of the vessel. Wall thickness located in the inner curvature were defined as that wall thickness values located closest to the radius vector (Fig 2⇓). The wall thickness located in the outer curvature was located at 180° from the location of the wall thickness of the inner curvature.
A local right-handed coordinate system was set up on the 3D vessel centerline according to Frenet.19 The three-unit vectors of this right-handed system are the tangent, the normal, and the binormal vectors.19 It can be shown that if the direction of the binormal vector19 remains unchanged, this 3D space curve remains in plane. Hence, changes in the direction of this unit vector imply out of plane curvatures (“twist”) of the coronary artery.
Computational Fluid Dynamics
Shear stress on the endothelium of a human vessel was calculated using the 3D mesh described above and consisting of 3D cells (“bricks”). In each cell, 27 nodes are positioned; and in each node, a differential equation describing the full incompressible Navier-Stokes equations are implemented. The nonlinear convective term in the Navier-Stokes equations is linearized with a Newton-Raphson approach. To obtain the pressure unknowns from the discrete Navier-Stokes equations, a penalty function approach is used. We have defined the following boundary conditions: uniform stationary inflow velocity of 10 cm/second, 20 cm/second, and 30 cm/second at the entrance of the vessel; constant zero stress at the outlet; and no-slip conditions at the vessel wall. The values of 10 cm/second, 20 cm/second, and 30 cm/second were chosen to represent data from patients under baseline conditions. We further assumed that blood was a newtonian fluid with a density of 1050 kg/m3 and a viscosity of 3 cP.20 21 22 Elimination of the pressure unknowns with the penalty method and linearization of the convective term results in a set of linear equations with velocity unknowns. These equations were solved by a numerical solver, applying a direct profile method of a well-validated,20 21 22 commercially available finite element package (Sepran, Sepra, Leiden, The Netherlands), which was implemented on a work station (HP 715/80, 256 Mbyte RAM). Details and elaborate validation of the solver have been presented before by different laboratories under different conditions.20 21 22 23
The present approach resulted in the calculation of the full 3D velocity, pressure, and velocity gradient fields. Because a quadratic interpolation scheme was used, the calculations resulted in approximately 25 000 nodal points in which the full 3D solution of the vector fields were evaluated. From these points, we selected for each cross-section the 16-velocity gradient values located near the vessel wall. These numbers were subsequently multiplied with viscosity to obtain shear stress at 16 points for each cross-section. For the entire vessel, with 100 cross-sections, this resulted in 1600 values of shear stress.
At corresponding locations on the vessel wall, we calculated shear stress, wall thickness, and curvature with the algorithms described above. The spatial resolution of the shear stress values was approximately 0.3 mm in the circumferential direction and 0.75 mm in the axial direction.
Analysis and Statistics
The data are presented as means±standard deviations. Differences in mean values are tested by the Student’s t test, and significance levels are set to P<.05. To relate wall thickness quantitatively to shear stress, we applied simple linear regression analysis using a standard software package (SPSS, Md). Quality of the fit was evaluated with analysis of variance and r2 values.
The 3D reconstruction of the lumen contours of a human right coronary artery by the ANGUS method is shown in the left panel of Fig 1⇑, and the right panel displays the surface of the mesh necessary for the 3D-velocity calculations. The gradient of the velocity profiles, ie, the shear rate, was used for calculations of shear stress.
Average diameter and wall (intima+media) thickness for the entire vessel was 2.3±0.4 mm and 0.65±0.37 mm (n=1600). The average maximal and minimal wall thicknesses per cross-section were 1.2±0.36 mm (n=100) and 0.26±0.08 mm (n=100, P<.05), respectively. The radius (R) of curvature of the reconstructed vessel was 40±21 mm, resulting in a curvature (1/R) of 25.09±8.96 (m−1) and a curvature ratio of approximately 1 to 17. Local twist of the catheter path was calculated from the rotation of the binormal, as defined by Frenet.18 The rotation per step of 0.75 mm along the vessel centerline was 0.1734±0.147°. Cumulative rotation was approximately 13°.
Based on the curvature calculations, we calculated average vessel wall thickness located at the inner and outer location of the curve as 0.90±0.41 mm (n=100) and 0.37±0.17 mm (n=100, P<.05), respectively. The location of the largest wall thickness value was close to the location of the inner curve, and the smallest wall thickness was close to the outer curve (Fig 3a⇓).
Shear stress values were calculated for three velocities (10 cm/second, 20 cm/second, and 30 cm/second), which are reported in this order. Averaged values of Reynolds number, per cross-section, for these inflow velocities were 114±22, 229±44, and 457±87 (n=100); the values of Dean number were 20±8, 46±16, and 93±33 (n=100). Average shear stress values for the entire right coronary artery were 20±2, 46±8, and 80±17 dyne/cm2 (n=1600). When the shear stress values were selected on the basis of the location of the inner and outer curvature, these values were 17±17, 38±44, and 77±54 dyne/cm2 (n=100) and 22±17, 60±44 and 107±79 dyne/cm2, respectively (n=100, P<.05). The location of maximal shear stress was closer to the outer wall, and the minimal shear stress was closer to the inner curve (Fig 3B⇑).
To evaluate the contribution of shear stress to the asymmetrical distribution of plaque in this right coronary artery, we averaged the wall thickness and shear stress values over the axial direction of the vessel (n=100). Shear stress and wall thickness appeared out of phase when plotted against the location of the circumference, ie, a low shear stress corresponded to a high wall thickness at the same location (Fig 4⇓). The resulting relationship between shear stress and wall thickness was inverse and could be described as linear for all three relationships (Table 1⇓ and Fig 4⇓).
Atherosclerosis is located nonuniformly over the arterial vascular system.1 2 3 4 5 The predilection sites of atherosclerosis as well as the eccentric localization of most advanced forms of atherosclerosis may be explained on the basis of hemodynamic factors.1 2 3 4 5 6 7 8 9 10 11 12 The most elaborate, studied factor associated with atherosclerosis is the shear stress, imposed by streaming of the blood along the endothelium.6 7 8 9 10 11 12 The exact implications of changes in shear stress are not yet known as different disturbances of the laminar velocity pattern have been associated with the degree and localization of atherosclerosis: low shear stress,9 10 11 low and oscillating shear stress,20 and high shear stress.21
The discrepancies between shear stress and degree and/or localization of atherosclerosis might be related to collection of the data.9 10 11 For instance, some studies related an average geometric model of the carotid artery to average histological parameters of the vessel wall.4 25 This averaging process might underestimate related phenomena as has been pointed out recently.9 10 11 To overcome this problem, Friedman et al and Gibson et al devised flowthrough casts in which they related locally derived shear rates with locally derived vessel wall thickness from postmortem material of humans.9 10 11 Although this approach improved the relationship between shear stress and intimal thickness of each specimen, this method is very time consuming, allows the collection of only a relatively small amount of data per specimen, and applies to postmortem material only. The latter constraint is especially important as it does not allow one to evaluate the effect of shear stress on different time points in a single human vessel. With the present technique, it is possible to evaluate the hemodynamic effect induced by interventions on restenosis and atherosclerosis without the complicating factor of a measurement device located inside the vessel disturbing the flow pattern.
Despite the fact that numerical techniques do not include the disadvantages of the hydraulic methods, so far only 2D geometries and recently idealized 3D geometries of vessels have been studied with this technique.20 21 As we studied these idealized geometries, it became clear that for the local shear stress calculations, the local geometry is very important. This parameter becomes even more important when local irregularities of a sick vessel induce local vertices, stagnation points, and separation of stream lines (“secondary velocity fields”). Furthermore, since the coronary vessel anatomy and geometry are not confined to a single 2D plane, it is necessary for the secondary velocity fields to be analyzed with 3D techniques as symmetry conditions are no longer applicable under these conditions.20 21 In addition, although idealized geometries are useful in understanding the relation of hemodynamics to atherosclerosis, coupling of a numerical technique to realistic 3D vessel geometries appears to be a necessary condition for understanding induction, progression, and regression of atherosclerosis in humans in vivo.
In this study, a numerical approach has been introduced to calculate shear stress on the vessel wall. Such methods do not have the constraints of the hydraulic in vitro setups described above; however, they rely on the grid spacing and accuracy of the vessel geometry obtained. Accuracy of the reproduction of the 3D geometry has been reported on recently15 16 17 and is in the order of 0.1 mm, which is accurate enough for the present approach.20 21 22 23 Accuracy of the numerical technique, related to the interpolation scheme, iterative solver used, and dimensions of the individual cells, have been compared with laser Doppler techniques; and it was shown that the method produced accurate results in 90° curved tubes, carotid bifurcations, and models of the basilar artery.20 21 22 23 Reduction of the mesh size by a factor of two increased the average shear stress by 7% and the slope of the inverse relationship of wall thickness and shear stress by only 2%. Hence, because in the present analysis we used a mesh resolution finer than in the previous studies, and a reduction of the mesh size by a factor of two in the present study did not affect relevant parameters, we conclude that our present mesh size was accurate.
We analyzed only stationary 3D-velocity fields in the present study. It has been shown that oscillating, low shear stress also may contribute to the development of atherosclerosis.25 With the present analysis, we cannot exclude that shear stress oscillation could also be involved.
We assumed that the coronary vessel wall is rigid. However, the coronary vessel wall is elastic and its radius changes by 5% to 10% over the cardiac cycle. When the effect of elasticity was included in a numerical analysis of the human carotid artery, it appeared that the effect of elasticity on shear stress is mainly periodic and does not affect the time-averaged shear stress.24 25 26 Hence, if this situation is applicable for the coronary artery, it could imply that time-averaged shear stress values are less affected by movement of the vessel wall. This may explain why a high correlation between spatially averaged shear stress and atherosclerosis was still found in the present study, despite the simplifications used.
We used uniform inflow velocity as one of the boundary conditions. As a result, the first few values of shear stress might be artificially high. However, a reanalysis of the relationship between shear stress and wall thickness omitting the first 10 points revealed a similar relationship, making this effect of minor importance.
As a first application, we studied a human right coronary artery with known atherosclerosis. This patient was sent to our hospital for anginal complaints and consequently, a percutaneous transluminal coronary angioplasty (PTCA) was performed to relieve these symptoms. The 3D reconstruction presented in this study was performed after the PTCA. To evaluate the effect of the PTCA on our findings, we also evaluated the relationship between wall thickness and shear stress for the segment of vessel upstream of the known location of the PTCA. The results still revealed a significant inverse relationship between wall thickness and shear stress.
Because of its nondestructive nature, the present method allows multiple measurements in the same patient. However, because coronary arteries are difficult to image noninvasively, we still need an invasive method when atherosclerosis and shear stress are evaluated in coronary arteries in humans in vivo.
This is the first study presenting a technique that makes it possible to evaluate the relationship of 3D coronary vessel geometry to shear stress in humans in vivo. With this technique, we were able to show that (1) wall thickness is thicker on the inner wall than on the outer wall in a right coronary artery in vivo, (2) the shear stress is low on the inner wall and high on the outer wall and (3) wall thickness is inversely related to shear stress. The present method enables one to evaluate the role of shear stress during the development of restenosis after coronary interventions or plaque progression/regression as it allows multiple measurements in the same patient.
We gratefully acknowledge funding (to J.J.W.) by the Inter-University cardiology Institute of the Netherlands (project ICIN 18) and the gift of the transformation matrix by Dr F. van de Vosse.
Tjotta E. The distribution of atheromatosis in the coronary arteries. J Atheroscler Res.. 1963;3:253-261.
Willis CG. Localizing factors in atherosclerosis. Can Med Assoc J.. 1954;70:1-9.
Smedby O, Johanson J, Molgaard J, Osen JA, Waldius G, Erikson U. Predilection of atherosclerosis for the inner curvature in the femoral artery. Arterioscler Tromb Vasc Biol.. 1995;15:912-917.
Sabbah HN, Khaja F, Hawkins ET, Stein P. Relation of atherosclerosis to arterial wall shear in the left anterior descending artery in man. Am Heart J. 86;112:453-458.
Friedman MH, Deters OJ, Bargeron CB, Hutchins GJM, Mark FF. Shear-dependent thickening of the human arterial intima. Arteriosclerosis.. 1986;60:161-171.
Gibson CM, Dlaz L, Kandarpa K, Sacks FM, Pasternak RC, Sandor T, Feldman CH, Stone PH. Relation of vessel wall shear to atherosclerosis progression in human coronary arteries. Arterioscler Tromb.. 1993;13:310-315.
Davies PF, Tripathi SC. Mechanical stress mechanism and the cell. Circ Res.. 1993;72:239-245.
Zho S, Sucin A, Ziegler T, Mone J, Burki E, Meister J, Brunner H. Synergistic effect of fluid shear stress and cyclic circumferential stretch on vascular endothelial cell morphology and cytoskeleton. Arterioscler Tromb Vasc Biol.. 1995;15:1781-1786.
Slager CJ, Laban M, von Birgelen C, Krams R, Oomen JAF, den Boer A, Wenguang L, de Feyter PJ, Serruys PW, Roelandt JRTC. ANGUS: a new approach to three-dimensional reconstruction of geometry and orientation of coronary lumen and plaque by combined use of coronary angiography and IVUS. J Am Coll Cardiol.. 1995;25:144A.
Laban M, Oomen JAF, Slager CJ, Wentzel JJ, Krams R, Schuurbiers JCH, den Boer A, Birgelen von C, Serruys PW, Feijter de PJ. ANGUS: a new approach to three-dimensional reconstruction of coronary vessels by combined use of angiography and intravascular ultrasound. Computers in Cardiology 1995, 95CH35874. Piscataway, NJ: IEEE Comp. Soc., 1995:325-328.
Slager CJ, Laban M, Oomen JAF, von Birgelen C, Li W, Krams R, Schuurbiers JCH, den Boer A, Serruys PW, Roelandt JRTC, de Feijter PJ. Three-dimensional geometry and orientation of coronary lumen and plaque: reconstruction from angiography and ICUS (ANGUS). Thoraxcenter J.. 1995;7:36-37.
Li W, Bom N, Van Egmond FC. Three-dimensional quantification of intravascular ultrasound images. J Vasc Invest.. 1995;1:2:57-61.
Davies HF, Snider AD, eds. Introduction to Vector Analysis. 5th ed. C. Brown, Dubuque, Iowa; 1988.
Rindt CCM, Steenhoven AA. Unsteady flow in a rigid 3-D model of the carotid artery bifurcation. Transaction of the ASME.. 1996;118:90-96.
Moore JF, Ku DN, Zarins CU, Glagov S. Pulsatile flow visualization in the abdominal aorta under differing physiological conditions: implications for increased susceptibility to atherosclerosis. J Biomed Eng.. 1992;114:391-397.