Hemodynamics of human carotid artery bifurcations: Computational studies with models reconstructed from magnetic resonance imaging of normal subjects

Hemodynamics of human carotid artery bifurcations: Computational studies with models reconstructed from magnetic resonance imaging of normal subjects Jaques S. Milner, BESc, Jennifer A. Moore, MSc, Brian K. Rutt, PhD, and David A. Steinman, PhD, London and Toronto, Canada Purpose: The precise role played by hemodynamics, particularly wall shear stress, in the development and progression of vascular disease remains unclear, in large part because of a lack of in vivo studies with humans. Although technical challenges remain for noninvasively imaging wall shear stresses in humans, vascular anatomy can be imaged with sufficiently high resolution to allow reconstruction of three-dimensional models for computational hemodynamic studies. In this paper we present an entirely noninvasive magnetic resonance imaging (MRI) protocol that provides carotid bifurcation geometry and flow rates from which the in vivo hemodynamics can be computed. Maps of average, oscillatory, and gradients of wall shear stress are presented for two normal human subjects, and their data are compared with those computed for an idealized carotid bifurcation model. Methods: An MRI protocol was developed to acquire all necessary image data in scan times suitable for patient studies. Three-dimensional models of the carotid bifurcation lumen were reconstructed from serial black blood MR images of two normal volunteers. Common and internal carotid artery flow rate waveforms were determined from MRI phase-contrast velocity imaging in the same subjects and were used to impose fully developed velocity boundary conditions for the computational model. Subject-specific timeresolved velocities and wall shear stresses were then computed with a finite element–based Navier-Stokes equation solver. Results: Models reconstructed from in vivo MRI of two subjects showed obvious differences in branch angle, bulb size and extent, and three-dimensional curvature. Maps of a variety of wall shear stress indices showed obvious qualitative differences in patterns between the in vivo models and between the in vivo models and the idealized model. Secondary, helical flow patterns, induced primarily by the asymmetric and curved in vivo geometries, were found to play a key role in determining the resulting wall shear stress patterns. The use of in vivo flow rate waveforms was found to play a minor but noticeable role in some of the wall shear stress behavior observed. Conclusions: Conventional “averaged” carotid bifurcation models mask interesting hemodynamic features observed in realistic models derived from noninvasive imaging of normal human subjects. Observation of intersubject variations in the in vivo wall shear stress patterns supports the notion that more conclusive evidence regarding the role of hemodynamics in vascular disease may be derived from such individual studies. The techniques presented here, when combined with subject-specific MRI measurements of carotid artery plaque thickness and composition, provide the tools necessary for entirely noninvasive, prospective, in vivo human studies of hemodynamics and the relationship of hemodynamics to vascular disease. (J Vasc Surg 1998;27:143-56.) From the Imaging Research Laboratories, John P. Robarts Research Institute (Milner and Drs. Rutt and Steinman); the Department of Medical Biophysics, University of Western Ontario, London (Drs. Rutt and Steinman); the Department of Diagnostic Radiology & Nuclear Medicine, University of Western Ontario, London (Dr. Rutt); and the Department of Mechanical and Industrial Engineering, University of Toronto (Moore). Support was provided by the Heart and Stroke Foundation of Ontario (grant NA3197).

Preliminary data were presented at the 1997 American Society of Mechanical Engineers Bioengineering Conference, Sun River OR, June 1997. Reprint requests: David A. Steinman, PhD, John P. Robarts Research Institute, 100 Perth Drive, P.O. Box 5015, London, Ontario N6A 5K8, Canada. Copyright © 1998 by the Society for Vascular Surgery and International Society for Cardiovascular Surgery, North American Chapter. 0741-5214/98/$5.00 + 0 24/1/87948

143

144 Milner et al.

JOURNAL OF VASCULAR SURGERY July 1998

Fig. 1. A custom bilateral phased-array radio-frequency coil, indicating the multiple translational and rotational degrees of freedom. The inset shows the position of a volunteer in the coil.

Despite many hemodynamic studies carried out with models of arterial bifurcations, especially the carotid artery bifurcation, the precise role played by wall shear stress in the development and progression of atherosclerosis remains unclear. Early studies correlated high,1 low,2 and oscillatory3 shear with the presence of arterial disease, but later studies postulated roles for the temporal4 and spatial5 gradients of wall shear stress. The inability to conclusively identify the hemodynamic quantities that influence atherosclerosis may be attributed in part to the indirect nature of the correlations made between hemodynamics and vascular disease: flow studies are typically carried out in idealized models with averaged flow parameters, and sites predisposed to atherosclerosis

are identified from averaged postmortem measurements. More direct studies, in which the presence or absence of disease can be compared with wall shear stress patterns throughout individual carotid bifurcations, probably are required to remove these remaining ambiguities. Magnetic resonance imaging (MRI) can produce high-resolution images of the human carotid bifurcation noninvasively and in scan times suitable for patient studies.6,7 Although MRI is also capable of directly measuring blood flow velocity, acquiring velocity data at spatial and temporal resolutions necessary for computing three-dimensional maps of wall shear stress requires scan times sufficiently long to preclude human studies. The studies that have

JOURNAL OF VASCULAR SURGERY Volume 28, Number 1

Milner et al.

145

extracted wall shear stress information from in vivo MRI velocity measurements have been limited to relatively thick slices in straight sections of the abdominal aorta.8,9 Wall shear stress patterns in complex, three-dimensional geometries such as the carotid bifurcation can be computed indirectly in models reconstructed from the MR-imaged geometry. Studies that have computed hemodynamics in models reconstructed from in vivo imaging have relied on invasive imaging techniques such as intravascular ultrasound10 and contrast-enhanced xray computed tomography.11 Although such techniques typically provide superior image quality compared with noninvasive imaging, which simplifies the task of model reconstruction, they have associated risks, making them inappropriate for repeated studies of otherwise healthy individuals. In this study we present an entirely noninvasive MRI protocol for deriving individual carotid bifurcation geometry and inlet-outlet flow rate waveforms in humans. Models reconstructed from this image-based information and subjected to computational hemodynamic analysis are used to produce subject-specific maps of a variety of important wall shear stress indices, which are contrasted with those from an idealized model. This combination of noninvasive imaging and modeling represents the critical first step toward the long-term, prospective studies we believe are required to elucidate the key role hemodynamics plays in vascular disease. MATERIALS AND METHODS Magnetic Resonance Imaging. MRI of two presumed-healthy male volunteers (28 and 32 years old) was performed using a clinical 1.5 T MRI scanner (Signa Horizon EchoSpeed v5.5; General Electric Medical Systems; Milwaukee, Wis.). The experimental protocol was approved by the University Review Board for Health Sciences Research Involving Human Subjects, and both subjects gave informed consent. Because the imaging protocol was designed ultimately for repeated studies on patients, emphasis was placed on minimizing individual and overall scan time while maintaining a target in-plane spatial resolution of 0.3 mm. A custom-mounted, bilateral, phased-array coil was developed to receive signals simultaneously from both carotid arteries (Fig. 1). Each two-element phased array of 5 cm, overlapping rectangular surface coils12 was mounted on a rigid “paddle” designed to conform to the jawline beneath which the carotid bifurcation is typically found. These paddles were mounted on a custom-

Fig. 2. A, Black blood MR image of a transverse slice through the neck of a normal human volunteer at a location just above the right carotid bifurcation apex. B, Cropped black blood images show the right carotid artery at selected locations inferior (a through e) and superior (f through i) to the bifurcation apex. JV identifies the jugular vein.

designed rig with multiple rotational and translational degrees of freedom, allowing the coil arrays to be positioned immediately adjacent to the neck and around the jaw to maximize the received signal. Compared with a conventional 8 cm surface coil, the phased-array coil showed a 50% improvement in signal-to-noise ratio in tests on human volunteers, equivalent to a more than twofold reduction in scan time, assuming the same spatial resolution and signal-to-noise ratio. After locating the carotid bifurcations with conventional phase contrast angiography, the protocol consisted of two pairs of scans: one to image vascu-

146 Milner et al.

JOURNAL OF VASCULAR SURGERY July 1998

B

A

C

Fig. 3. Magnitude (A) and phase (B) images at a location approximately 2 cm inferior to the bifurcation apex. Shown also are cropped phase and velocity images of the right common carotid artery (CCA) at 16 times that are equally spaced throughout the cardiac cycle (C). Notice the asymmetric velocity patterns, particularly evident in the third and fourth images, indicating the presence of undeveloped, helical flow in the CCA.

lar anatomy, the other to image entrance and exit velocity fields. In the first pair of scans, cardiac-gated (by peripheral probe) black blood MRI was used to acquire a total of 32 × 2 mm thick, transverse, contiguous slices (0.313 mm in-plane spatial resolution) centered approximately at the bifurcation apex. The first scan acquired 16 slices through the common carotid arteries (CCA); the second acquired images of the internal and external carotid arteries (ICA and ECA). This particular arrangement of slices proved optimal for avoiding plaque-mimicking artifacts we have identified in black blood MRI.7 Scan parameters included two-dimensional fast spin echo (FSE), 8 cm superior and inferior saturation bands, 16 × 12 cm field of view, 512 × 384 acquisition, 2 R-R (approximately 2 seconds) repetition time, 15 millisecond (ms) effective echo time, 4 echo train length, and two scans that were averaged. Each 16slice scan was 4 to 5 minutes long, depending on the subject’s heart rate. This particular combination of imaging parameters provided the highest-quality images in the least amount of scan time. Gating parameters were selected such that all slices were acquired at least 100 ms after peak systole to ensure

a consistent reconstruction of the artery in its diastolic state. Images were then corrected for nonuniform signal intensity (because of use of the surface radio-frequency coils) with the technique of Murakami et al.13 Representative images are shown in Figure 2. In the second pair of scans, retrospectively gated cine phase-contrast MRI was used to image the timevarying velocity field through 4 mm thick transverse slices of 0.469 mm in-plane resolution ± 2 cm from the bifurcation apex. Scan parameters included a 12 × 12 cm field of view, 256 × 256 acquisition, 18 ms repetition time, 6.8 ms echo time, first-order flow compensation, and 100 cm/sec maximum encoded velocity. Images were reconstructed at 16 equally spaced phases of the cardiac cycle. Each scan was approximately 4 minutes long. To avoid low-pass filtering of the imaged velocities while maintaining a minimum repetition time,14 only the superiorinferior (axial) component of velocity was encoded; only the axial component of velocity was required to compute flow rates through the arteries, because the image planes were perpendicular to the axial direction. Representative images are shown in Fig. 3.

JOURNAL OF VASCULAR SURGERY Volume 28, Number 1

Milner et al.

147

Fig. 4. Stages in the reconstruction of model II. A, The lumen contours were derived from the MR images; the contour identified by the arrow corresponds to the terminal common carotid artery (CCA) contour. B, Interpolated and extension contours are added, highlighting the intersecting internal carotid artery (ICA) and external carotid artery (ECA) terminal contours. C and D, Compare the matched surface-spline representation and surface of the quadratic, tetrahedral-volume finite element mesh.

Reconstruction of Lumen Geometry. The lumen boundary from each black blood image was segmented in a two-step process. First, the boundary was estimated with conventional region-growing in the vessel of interest to produce a set of points ordered counterclockwise around the circumference, as required for the subsequent solid modeling stage of the reconstruction. These data points were then fitted to a B-spline surface, the smoothness of which was controlled by varying the number of control points. The points were manually adjusted to provide the (operator-perceived) best fit to the lumen. Typically, 8 to 12 control points provided an optimal balance between smoothing and user control. When the external carotid artery itself branched beyond the bifurcation, only those lumen contours

from the ECA inferior to the branch point were included in the model. The lumen contours were imported into a computer-aided design program (ICEM CFD/CAE; ICEM Engineering, Berkeley, Calif.), and “lofted” together to produce a continuous surface representation of the lumen (Fig. 4). Because it is beyond the capability of a single B-spline surface to conform to a bifurcating geometry, we fit by approximating single-segment B-spline surfaces (typically of order 7 to 10, depending on the branch curvature and tortuosity) to each of the external, internal, and common branches. Because single-segment splines are required only to follow the contour data approximately, they provide an additional level of threedimensional smoothing, depending on their order.15

148 Milner et al.

JOURNAL OF VASCULAR SURGERY July 1998

Fig. 5. Surface-shaded views of the two in vivo carotid bifurcation models reconstructed from black blood MRI and an idealized model. The common carotid artery (CCA), internal carotid artery (ICA), and external carotid artery (ECA) flow rate waveforms were reconstructed from the phase-contrast images (models I and II) or assumed from the literature (idealized model).

To ensure continuity of the three separate surfaces, each was required to terminate at the same contour level. To this end, the last lumen contour inferior to the bifurcation apex was first manually segmented from the MR image as two individual, overlapping, closed contours, representing the inferior of the ICA and ECA spline surfaces (Fig. 4B). A single, closed contour was then derived from these two intersecting contours, representing the end of the CCA carotid spline surface (Fig. 4A). The relatively coarse 2 mm spacing between the lumen contours was initially found to produce poor matching of the surface splines, often resulting in a sharp seam dividing the CCA and the ICA or ECA surfaces. This was ultimately avoided by providing additional constraints for approximating spline surfaces in the form of contours linearly interpolated between each pair of MR-derived lumen contours (Fig. 4B).15 Combined with cylindrical flow extensions to facilitate the prescription of fully developed flow boundary conditions, a continuous, smoothed surface model was produced (Fig. 4C). An octree-

based, unstructured, finite element mesh generator (TETRA; ICEM Engineering, Berkeley, Calif.) was used to produce a volume mesh of approximately 35,000 to 40,00010-node, quadratic, tetrahedralvolume finite elements with between and the 55,000 and 65,000 nodes (Fig. 4D). This spatial resolution was deemed sufficient, because studies with a symmetric idealized carotid bifurcation model showed that doubling the number of nodes produced no qualitative change in the wall shear stress patterns. Reconstruction of Flow Waveforms. The retrospectively gated phase-contrast scans produced 16 magnitude (anatomic) and phase (velocity) images, equally spaced throughout the cardiac cycle, each at the level of the CCA and ICA away from the apex. At each time, the lumen was identified in the magnitude image by masking pixels with a signal above a user-defined threshold and then applying the resulting mask to the corresponding phase image. Mean velocity was computed as the average pixel value within this mask. Background phase roll, determined by averaging the phase in several small

JOURNAL OF VASCULAR SURGERY Volume 28, Number 1

Milner et al.

149

Table I. Common carotid artery parameters for the study models. Model Idealized In vivo I In vivo II

Diameter (mm) 6.3 6.2 6.8

Mean flow (ml/s) 6.0 7.6 6.7

regions of known static tissue near the masked region, was subtracted from the computed mean velocity from each image. The flow rate was computed as the mean velocity multiplied by the area of the masked region, and a cubic spline fit was assumed among flow rate data points. To ensure conservation of mass in the finite element model, flow rates in the ECA were calculated by taking the difference between the measured CCA and ICA flow rates. Direct measurement of flow in the ECA was avoided for two reasons: away from the bifurcation, the ECA branches multiple times, making it extremely difficult to image flow in all of the small adjoining vessels, and close to the bifurcation, artifacts arising from complex flow can introduce errors in the velocity measurements.16 Computational Fluid Dynamic Modeling. Modeling of the in vivo hemodynamics was carried out using a well-validated, in-house finite element–based Navier-Stokes solver.17 Fully developed flow was applied at the model CCA and ECA by imposing the time-varying velocities computed by the solution for pulsatile flow in a straight, rigid tube on the basis of the respective Fourier-decomposed flow rate waveforms.18 Traction-free boundary conditions were applied at the ICA, and no-slip (i.e., zero velocity) boundary conditions were applied at the walls. Rigid walls and a constant blood viscosity of 3.5 centipoise were assumed for all studies, because vessel wall distensibility and non-Newtonian viscosity have previously been shown to have only minor effects on the distribution of wall shear stress patterns.19-22 For each model, the time-dependent NavierStokes equations were solved at 480 time-steps per cardiac cycle, and two cardiac cycles were required to fully damp initial transients. Computations were carried out on a 180 MHz R10000 Indigo/2 workstation (Silicon Graphics, Mountain View, Calif.), with each study requiring approximately 60 CPU hours. Wall shear stresses and derived indices were determined directly from the computed velocity field and visualized using Tecplot software (Amtec Engineering, Bellevue, Wash.). All dimensional wall shear stress quantities were normalized to the

Peak flow (ml/s) 23.5 26.0 21.3

Reynolds number 346 444 360

Womersley number 4.2 4.3 4.4

respective values at the CCA inlet to facilitate intermodel comparisons. Wall shear stress spatial gradients, computed using the definition of Kleinstreuer et al.,23 were further normalized to the inlet radius. RESULTS Pulsatile hemodynamics were computed for three different carotid bifurcation geometries: in vivo models I and II, representing the right carotid bifurcations of the 28- and 32-year-old volunteers, respectively, and the idealized normal carotid bifurcation model of Smith et al.24 (Fig. 5). For the two in vivo models, the measured flow rates were imposed as described in the Methods section. For the idealized model, two different flow rate waveforms were applied: a representative CCA flow rate waveform averaged from Doppler ultrasound measurements in normal subjects,25 with a constant ICA:ECA flow split of 56:44, and the measured flow rates from model II. The latter case was included to identify the relative contribution of geometry and flow rate waveform to differences observed in the computed hemodynamics. Measured and assumed quantities at the common carotid for each of these models are summarized in Table I. Compared with model I, model II had a more pronounced curvature of the internal branch in and out of the plane of the CCA, a wider branch angle, a larger bulb, and more gradual tapering of the CCA toward the bifurcation. In contrast to the reconstructed models, branches in the idealized model were co-planar, circular, and symmetric about the plane of the bifurcation. The CCA flow rate waveforms all showed qualitatively similar behavior. The ICA:ECA flow rate division, however, was highly time dependent, in contrast to the constant flow split typically assumed for idealized model studies. During systole, flow in the ECA exceeded that in the ICA by a factor of two or more, an observation attributed to differences in the impedance of the vascular beds fed by the ICA and ECA.26 During diastole this trend was reversed, with most of the flow entering the ICA. Averaged over the cardiac cycle, the ICA:ECA flow split in the in vivo models was approximately 75:25.

150 Milner et al.

The time-averaged, normalized wall shear stress magnitude shown in Fig. 6 revealed in both in vivo models peak values at the bifurcation apex and moderately elevated values spiraling around the ICA from the inner (flow divider) wall along the superior direction. The lowest average shear was seen primarily at the base of the bulb. In model I, the region of lowest average shear was isolated to a thin strip along the posterior wall and around to the outer wall at the base of the bulb, but it extended around the circumference at the corresponding axial location in model II. A region of elevated average shear was present on the anterior wall of the CCA in model I, a feature entirely absent in model II. In marked contrast to the in vivo models, the idealized model showed much more elevated values of average shear in the ECA and ICA, regardless of the flow waveform used. As with the in vivo models, minimum values of average shear were present along the outer wall of the ICA bulb, the distribution of which was affected by the choice of flow rate waveform; in contrast to the in vivo models, values were generally higher on the posterior and anterior walls at this axial location. Perhaps the most striking difference between the in vivo and idealized models was the “striping” or helical appearance of the contours in the ICA in the former, particularly near the bulb. Inspection of the average velocity patterns in the in vivo models (not shown) revealed path lines twisting up through the ICA, following the same direction and overall pitch as these contours. The “pitch” of these helical contours was steeper in model I than in II. Maps of the dimensionless oscillatory shear index (OSI)27 shown in Fig. 7 revealed elevated values on the outer walls of the ICA and ECA proximal to the bifurcation apex in both in vivo models. However, model I showed low oscillatory shear on the anterior wall proximal to the bifurcation apex, whereas this was a region of relatively high OSI in model II. The region of elevated oscillatory shear in the CCA extended further inferiorly along the posterior wall of model I than model II. In general, the posterior-anterior asymmetry in OSI in both models was similar. In the idealized model, oscillatory shear in the ICA was similar to that on the anterior wall of the in vivo models. Although patterns of elevated oscillatory shear on the inner wall of the ECA were also similar to those of the in vivo models, the elevated values seen throughout the rest of the ECA of the in vivo models, as well as on the posterior walls of the CCA, were entirely absent in the idealized model. Overall, the use of the in vivo flow rate waveform in the idealized model had

JOURNAL OF VASCULAR SURGERY July 1998

little effect on OSI patterns, apart from slightly reducing the extent of the high OSI region in the bulb, and introduced modest oscillatory shear in the ECA. The peak wall shear stress temporal (WSST) gradient, mapped in Fig. 8, revealed much higher values in the ICA of model I than in model II or the idealized model. All models, however, predicted similar WSST gradient levels in the CCA proximal to the bifurcation apex. Contrary to what was observed for average and oscillatory shear, introduction of the model II flow waveform in the idealized model produced larger WSST gradient values than observed for the idealized waveform case. Wall shear stress spatial gradients, shown in Fig. 9, were consistently high around the bifurcation apex of all models and at the sharp transitions of the idealized model, regardless of the flow waveform used. Moderately high values were also seen on the outer walls of the ICA and ECA of the in vivo models; this trend was not evident in the idealized geometry. DISCUSSION With the use of a small sample, we demonstrated that the hemodynamics of models derived from in vivo imaging can vary from subject to subject. These models produced wall shear stress patterns markedly different from those of a conventional idealized carotid bifurcation geometry. Because realistic geometry and flow rates were employed, an obvious question arises about whether these variations were attributable to geometric or flow rate differences (or both equally). Our observations indicate that geometry plays the key role in determining the nature of the wall shear stress patterns. For example, the curvature of the in vivo models induced secondary, helical flow patterns that were apparent in the spiraling contours seen for the average wall shear stresses and that contributed to the strong anterior-posterior asymmetry of the hemodynamics patterns. Such features were absent in the idealized models. Further evidence is provided by the data from the idealized model with the model II flow rate waveform; these were much more similar to those from the same geometry with the idealized waveform than those from the model II geometry with the same waveform. Several differences among the models were directly attributable to flow rate waveform effects. These included the clear relationship between oscillatory shear magnitude and the amount of negative flow present in the ECA and the increased magnitude of temporal gradients in the idealized model using the model II waveform. Because of this limit-

JOURNAL OF VASCULAR SURGERY Volume 28, Number 1

Milner et al.

151

Fig. 6. Time-averaged wall shear stress magnitude (normalized to the value at the common carotid artery) for the in vivo and idealized models. Anterior and posterior views are provided to identify asymmetries of the patterns. Notice the contour levels identified in the lower right corner.

ed data, we believe it is important to include subjectspecific geometry and flow rates in modeling in vivo flow patterns; however, as suggested by Perktold and Resch,28 geometry plays an important role in determining local hemodynamics.28

Helical flow in the ICA (induced primarily by asymmetry and branch curvature in the in vivo models) played a major role in determining the distribution of the wall shear stresses. This suggests the intriguing possibility that the conventional

152 Milner et al.

JOURNAL OF VASCULAR SURGERY July 1998

Fig. 7. Oscillatory shear index for the in vivo and idealized models. Notice the contour levels identified in the lower right corner.

assumption of fully developed flow, which imposes axisymmetric flow in the CCA entering the bifurcation, may have an effect on the computed flow patterns. In this study, fully developed velocity boundary conditions were applied to the models on the basis of the measured flow rates. However,

inspection of Fig. 3 shows that axial velocity in the CCA even two diameters inferior to the bifurcation apex was not axisymmetric; its appearance suggested helical rather than fully developed flow. Helical flow in the CCA was seen in most of the normal volunteers (n = 5) and patients with early carotid

JOURNAL OF VASCULAR SURGERY Volume 28, Number 1

Milner et al.

153

Fig. 8. Peak wall shear stress temporal gradient (normalized to the value at the common carotid artery) for the in vivo and idealized models. Notice the contour levels identified in the lower right corner.

disease (n = 4) who were imaged, and it has been observed in MRI measurements by Caro et al.29 It is therefore possible that undeveloped, helical flow in the CCA, even when imposed in an idealized model, may produce wall shear stress patterns quite

different from those computed assuming fully developed flow. Although this may appear to contradict our assertion that geometry plays the dominant role in determining wall shear stress patterns, it does support the notion that secondary, helical

154 Milner et al.

JOURNAL OF VASCULAR SURGERY July 1998

Fig. 9. Wall shear stress spatial gradient (normalized to the average inlet shear stress divided by the inlet radius) for the in vivo and idealized models. Notice the contour levels identified in the lower right corner.

flow, whether induced by geometric, flow rate, or entrance flow characteristics, may be a key factor in determining the distribution of wall shear stresses and, by extension, vascular disease. The effect of helical compared with axisymmetric entrance flow

on downstream flow patterns will be reported in another article. In both volunteer studies, ICA flow rates were larger than CCA flow rates at the time of the minimum flow, producing a short period of reverse flow

JOURNAL OF VASCULAR SURGERY Volume 28, Number 1

in the derived ECA flow rate waveforms; this effect was most pronounced in model I. Although it is possible that even small errors in the extraction of the ICA and CCA flow rate data (from which ECA flow was derived) may account for this finding, reverse flow in the ECA has been reported by Aldrete et al.30,31 The findings suggest that the effects of such negative flow, most apparent in the elevated OSI and reduced average shear stresses in the ECA of the in vivo models, may be genuine and significant. Further careful studies are needed to confirm the finding of transient reverse flow in the ECA of normal human volunteers. In the absence of a reliable and noninvasive standard, it is difficult to directly assess the absolute accuracy of the in vivo geometry and flow rate waveform reconstructions in normal subjects. When we used an in vitro carotid bifurcation model imaged under conditions admittedly more ideal than those in this study (e.g., higher signal-to-noise ratio, no motion artifacts), we found that our reconstruction technique can accurately reproduce the true hemodynamics of a known complex arterial geometry.15 Qualitative features of the wall shear stress patterns remained even in the presence of modest reconstruction errors. The reconstructed flow rate waveforms compare favorably in terms of magnitude and overall shape with ultrasound measurements in normal subjects.32,33 Phenomena that might have been attributed to measurement error (e.g., reverse flow in the ECA) were found to have precedents in the literature, and results of our flow waveform sensitivity studies in the idealized model indicated only moderate differences in wall shear stress patterns induced by a large change in the applied flow rate waveform. We therefore conclude that it is extremely unlikely that errors in flow rate or geometry reconstruction could explain the pronounced qualitative differences observed in the hemodynamics of the models studied. However, we also believe that extension of this novel combined vascular imaging and modeling paradigm from the qualitative to the quantitative realms will require rigorous accuracy and precision studies. In summary, we have shown that the distributions of wall shear stress variables hypothesized to play a role in development and progression of vascular disease can vary from subject to subject. This finding supports the notion that more conclusive evidence regarding the role of hemodynamics in vascular disease may be derived from studies in individual subjects. The techniques developed and applied in this study, when combined with concomitant in

Milner et al.

155

vivo measurements of plaque distribution and composition (possibly derived from the same MR images used to reconstruct the geometry), provide the means necessary to carry out such studies. The observation that geometry plays a key role in determining local wall shear stress patterns in principle supports the notion of a “geometric risk factor” in atherosclerosis.34-36 Idealized models were found to mask interesting and potentially important in vivo hemodynamic behavior. Although such models have in the past proved invaluable for providing hemodynamic measurements with which pathologic data could be correlated, further advances in our understanding of the complex relationship between hemodynamics and vascular disease will probably come only when individual subjects are studied in a prospective manner. REFERENCES 1. Fry DL. Acute vascular endothelial changes associated with increased blood velocity gradients. Circ Res 1968;22:165-97. 2. Caro CG, Fitz-Gerald JM, Schroter RC. Atherosclerosis and arterial wall shear: observations, correlation and proposal of a shear dependent mass transfer mechanism for atherogenesis. Proc R Soc London B 1971;177:109. 3. Ku DN, Giddens DP, Zarins CK, Glagov S. Pulsatile flow and atherosclerosis in the human carotid bifurcation: positive correlation between plaque location and low oscillating shear stress. Arteriosclerosis 1985;5:293-302. 4. Ojha M. Spatial and temporal variations of wall shear stress within an end-to-side arterial anastomosis model. J Biomech 1993;26:1377-88. 5. Kleinstreuer C, Nazemi M, Archie JP. Hemodynamics analysis of a stenosed carotid bifurcation and its plaque-mitigating design. J Biomech Eng 1991;113:330-5. 6. Yuan C, Murakami JW, Hayes CE, et al. Phased-array magnetic resonance imaging of the carotid artery bifurcation: preliminary results in healthy volunteers and a patient with atherosclerotic disease. J Magn Reson Imaging 1995;5: 561-5. 7. Steinman DA, Rutt BK. On the nature and reduction of plaque mimicking flow artifacts in black blood MRI of the carotid bifurcation. Magn Reson Med. In Press. 8. Oshinski JN, Ku DN, Mukundan SJ, Loth F, Pettigrew RI. Determination of wall shear stress in the aorta with the use of MR phase velocity mapping. J Magn Reson Imaging 1995;5:640-7. 9. Oyre S, Pedersen EM, Ringgaard S, Boesiger P, Paaske WP. In vivo wall shear stress measured by magnetic resonance velocity mapping in the normal human abdominal aorta. Eur J Vasc Endovasc Surg 1997;13:263-71. 10. Chandran KB, Vonesh MJ, Roy A, et al. Computation of vascular flow dynamics from intravascular ultrasound images. Med Eng Phys 1996;18:295-304. 11. Taylor CA, Hughes TJ, Zarins CK. Computational investigations of vascular disease. Comput Physics 1996;10:224-32. 12. Hayes CE, Mathis CM, Yuan C. Surface coil phased arrays for high-resolution imaging of the carotid arteries. J Magn Reson Imaging 1996;6:109-12.

156 Milner et al.

13. Murakami JW, Hayes CE, Weinberger E. Intensity correction of phased-array surface coil images. Magn Reson Med 1996;35:585-90. 14. Frayne R, Rutt BK. Frequency response of retrospectively gated phase-contrast MR imaging: effect of interpolation. J Magn Reson Imaging 1993;3:907-17. 15. Ethier CR, Moore JA, Steinman DA, Holdsworth D. Construction of vascular models from MRI: accuracy assessment based on flow simulations in a carotid bifurcation model [Abstract]. Int J Cardiovasc Med Sci 1997; 1:73. 16. Steinman DA, Ethier CR, Rutt BK. Combined analysis of spatial and velocity displacement artifacts in phase contrast measurements of complex flows. J Magn Reson Imaging 1997;7:339-46. 17. Zhang X, Ethier CR, Karpik SR, Steinman DA. Development of a 3D unsteady Navier-Stokes solver. Proceedings of the Computational Fluid Dynamics Society of Canada; Toronto. Toronto: Computational Fluid Dynamics Society of Canada, 1994:379-86. 18. Frayne R, Steinman DA, Ethier CR, Rutt BK. Accuracy of MR phase contrast velocity measurements for unsteady flow. J Magn Reson Imaging 1995;5:428-31. 19. Ballyk PD, Steinman DA, Ethier CR. Simulation of nonNewtonian blood flow in an end-to-side anastomosis. Biorheology 1994;31:565-86. 20. Steinman DA, Ethier CR. The effect of wall distensibility on flow in a two-dimensional end-to-side anastomosis. J Biomech Eng 1994;116:294-301. 21. Perktold K, Resch M, Florian H. Pulsatile non-Newtonian flow characteristics in a three-dimensional human carotid bifurcation model. J Biomech Eng 1991;113:464-75. 22. Perktold K, Rappitsch G. Computer simulation of local blood flow and vessel mechanics in a compliant carotid artery bifurcation model. J Biomech 1995;28:845-56. 23. Kleinstreuer C, Lei M, Archie JP Jr. Flow input waveform effects on the temporal and spatial wall shear stress gradients in a femoral graft-artery connector. J Biomech Eng 1996; 118:506-10. 24. Smith RF, Rutt BK, Fox AJ, Rankin RN, Holdsworth DW. Geometric characterization of stenosed human carotid arteries. Acad Radiol 1996;3:898-911.

JOURNAL OF VASCULAR SURGERY July 1998

25. Holdsworth DW, Norley CJ, Frayne R, Steinman D, Rutt BK. Variability in common carotid blood flow waveforms [Abstract]. Radiology 1995;197P:451. 26. Spence JD. Spectral analysis of carotid vs femoral Doppler velocity patterns: a clue to genesis of flow disturbances in cerebral arteries? EEE 1981 Frontiers of Engineering in Health Care; Houston. New York City: 1981:355-9. 27. He X, Ku DN. Pulsatile flow in the human left coronary artery bifurcation: average conditions. J Biomech Eng 1996;118:74-82. 28. Perktold K, Resch M. Numerical flow studies in human carotid artery bifurcations: basic discussion of the geometric factor in atherogenesis. J Biomed Eng 1990;12:111-23. 29. Caro CG, Dumoulin CL, Graham JM, Parker KH, Souza SP. Secondary flow in the human common carotid artery imaged by MR angiography. J Biomech Eng 1992;114:147-9. 30. Aldrete JA, Narang R, Sada T, Tan LS, Miller GP. Reverse carotid blood flow—a possible explanation for some reactions to local anesthetics. J Am Dent Assoc 1977;94:1142-5. 31. Aldrete JA, Romo-Salas F, Arora S, Wilson R, Rutherford R. Reverse arterial blood flow as a pathway for central nervous system toxic responses following injection of local anesthetics. Anesth Analg 1978;57:428-33. 32. Eriksen M. Effect of pulsatile arterial diameter variations on blood flow estimated by Doppler ultrasound. Med Biol Eng Comput 1992;30:46-50. 33. Smith RF, Holdsworth DW, Rutt BK. Sophisticated test object for quantitative evaluation of MR angiographic sequences [Abstract]. Radiology 1995;197:188. 34. Caplan LR, Baker R. Extracranial occlusive vascular disease: does size matter? Stroke 1980;11:63-6. 35. Friedman MH, Deters OJ, Mark FF, Bargeron CB, Hutchins GM. Arterial geometry affects hemodynamics: a potential risk factor for atherosclerosis. Atherosclerosis 1983;46: 225-31. 36. Harrison MJ, Marshall J. Does the geometry of the carotid bifurcation affect its predisposition to atheroma. Stroke 1983;14:117-8.

Submitted Aug. 21, 1997; accepted Dec. 3, 1997.