- Email: [email protected]
MRI and CFD studies of pulsatile flow in healthy and stenosed carotid bifurcation models
ARTICLE IN PRESS
Journal of Biomechanics 37 (2004) 679–687
MRI and CFD studies of pulsatile flow in healthy and stenosed carotid bifurcation models Ian Marshalla,*, Shunzhi Zhaob, Panorea Papathanasopouloua, Peter Hoskinsa, X Yun Xub a
Department of Medical and Radiological Sciences, Medical Physics, University of Edinburgh, Western General Hospital, Edinburgh EH4 2XU, UK b Department of Chemical Engineering and Chemical Technology, Imperial College, London, UK Accepted 18 September 2003
Abstract Pulsatile flow was studied in physiologically realistic models of a normal and a moderately stenosed (30% diameter reduction) human carotid bifurcation. Time-resolved velocity measurements were made using magnetic resonance imaging, from which wall shear stress (WSS) vectors were calculated. Velocity measurements in the inflow and outflow regions were also used as boundary conditions for a computational fluid dynamics (CFD) model. Experimental flow patterns and derived WSS vectors were compared qualitatively with the corresponding CFD predictions. In the stenosed phantom, flow in the bulb region of the ‘‘internal carotid artery’’ was concentrated along the outer wall, with a region of low and recirculating flow near the inner wall. In the normal phantom, the converse was found, with a low flow region near the outer wall of the bulb. Time-averaged WSS and oscillatory shear index were also markedly different for the two phantoms. r 2003 Elsevier Ltd. All rights reserved. Keywords: Carotid blood flow; Velocity measurement; Wall shear stress; CFD; Stenosis
1. Introduction Over the last few decades, research has demonstrated the influence of vessel geometry and haemodynamic forces on the development of vascular pathology (Ku et al., 1985; Glagov et al., 1988). In regions of disturbed blood flow, which are found around vessel bifurcations, there is an increased chance of atheroma deposition. The viscous drag of the blood on the vessel wall (wall shear stress; WSS) is believed to play a part in regulating arterial structure (Caro et al., 1971; Giddens et al., 1993; Gnasso et al., 1997; Malek et al., 1999), and is implicated in plaque development. The exact relationship between atheroma deposition and WSS remains uncertain, with low WSS, oscillating WSS, and WSS gradient all having been considered (Ku et al., 1985; Caro et al., 1971; Giddens et al., 1993). There is considerable interest in the development of in vivo techniques that would allow accurate blood flow quantification and estimation of WSS, but they are not *Corresponding author. Tel.: +44-131-537-1661; fax: +44-131-5371026. E-mail address: [email protected] (I. Marshall). 0021-9290/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2003.09.032
yet available for clinical use. The two most promising techniques are Doppler ultrasound and magnetic resonance imaging (MRI). Doppler ultrasound is a real time technique and does not suffer from data loss during turbulent flow. Current commercial systems are based on acquisition of two-dimensional (2D) image data in which only a single component of velocity is acquired. Acquisition of 2 or 3 velocity components requires multi-beam systems which are the subject of current research (Dunmire et al., 2000). Acquisition of real time three-dimensional (3D) information is also possible using 2D array systems (Light et al., 1998), but these are at an early stage of development. Ultrasound cannot be used on arteries that lie behind bony or air-filled structures. MRI is intrinsically 3D, but suffers from long acquisition times and signal loss in turbulent flow, which limits its use in diseased arteries. Both techniques have spatial resolutions that are only just high enough (0.5 mm) for in vivo studies of carotid arteries. WSS estimation methods based on MRI have to date been mainly restricted to measurement of the axial component at a single plane (Oshinski et al., 1995; Oyre et al., 1998; Stokholm et al., 2000). Similarly, ultrasound methods have been based on estimation of axial WSS at
ARTICLE IN PRESS 680
I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
a single location (Brands et al., 1995; Hoeks et al., 1995). Vector MRI information was obtained at the intersection of selected planes in the aorta by Suzuki et al. (1998). Time-resolved measurement of in vivo flow and WSS vectors remains a challenging goal. Therefore, the study of fluid dynamics is often carried out in physical or computational models (Zhao et al., 2000). For useful reviews, see Ku (1997) and Berger and Jou (2000). There have been very few experimental studies of pulsatile flow in carotid models, especially in stenosed models. Palmen et al. (1994) used hydrogen bubble visualisation, Gijsen et al. (1996) used laser Doppler anemometry, and Botnar et al. (2000) used MRI. We have previously reported on the estimation of WSS vectors throughout a flow region from MRI velocity measurements of steady . (Kohler et al., 2001) and pulsatile (Papathanasopoulou et al., 2003) flow in phantoms. An alternative approach has been to use MRI structural data and velocity measurements as the boundary conditions to solve the Navier–Stokes equations in a computational fluid dynamics (CFD) package, from which WSS can be calculated (Long et al., 2000; Steinman et al., 2002). In the present work, we used a time-resolved 3D phase-contrast (PC) MRI sequence to acquire all three velocity components of pulsatile flow in models of a human carotid bifurcation with and without a 30% (diameter) axisymmetrical stenosis. We show velocity vectors at selected time frames, and WSS vectors calculated from the MRI measurements. These are compared qualitatively with the predictions of CFD calculations. The comparison of normal and stenosed phantoms under physiologically realistic flow conditions is novel, as is the qualitative comparison of MRI and CFD results for the stenosed phantom.
(vx ; vy ; vz ) to obtain a temporal resolution of 50 ms. We . implemented a retrospectively gated method (Wigstrom et al., 1996; Papathanasopoulou et al., 2003). The sequence was run with 18 repetitions to ensure that each period (‘‘RR-interval’’) of the pulsatile waveform was adequately sampled for each line of k-space. The pixel size was 0.5 mm, and 48 slices of 1.4 mm thickness were used to cover the bifurcation region. The encoding velocity (venc ) was 120 cm/s. Raw k-space data files were saved for offline reconstruction. Trigger pulses from the pump were logged by a timing system (IFIS; Psychology Software Tools, Pittsburgh, PA) which was synchronised to the scanning. Acquisition time was approximately 90 min for each velocity component. 2.2. Image reconstruction and processing
2. Methods
The raw k-space data from the scanner and the data from the timing system were transferred to a UNIX workstation for reconstruction. Data were sorted into the 18 desired time frames by linear interpolation between the two nearest available samples. k-space data were zero-filled to 64 slices, and a 3D FFT applied to generate images with 1.05(=1.4 48/64) mm effective slice thickness. Finally, velocity images were generated by phase subtraction of the reference images from the corresponding velocity-encoded images. Reconstructed magnitude and velocity images for each time frame were processed as described previously, using software devel. oped in-house (Kohler et al., 2001). Flow regions were automatically segmented, and the velocity field fitted by a fifth order polynomial model. The downstream branches (‘‘internal carotid artery’’ (ICA) and ‘‘external carotid artery’’ (ECA)) were fitted separately from the ‘‘common carotid artery’’ (CCA). Vessel walls were determined from the segmentation results, and WSS vectors were calculated from the spatial derivatives of the fitted velocity model.
2.1. Data acquisition
2.3. CFD
MR scanning was performed on a 1.5 T GE Signa scanner (GE Medical, Milwaukee), using a 75 mm surface coil. The flow phantoms were anthropomorphically realistic models of a human carotid bifurcation with an axisymmetric 30% diameter (50% area) stenosis and a healthy carotid bifurcation (Shelley Medical Imaging Technologies, London, Ont., Canada). A computer-controlled pump (UHDC Flow System, Shelley) was used to generate a realistic carotid waveform, and pumped blood-mimicking fluid (Shelley) through the phantoms. The fluid viscosity was measured with a capillary viscometer and determined to be 3.4 103 N s/m2. Velocity measurements were made using a 3D-PC sequence, run separately for each velocity direction
The ‘‘magnitude’’ images from the PC sequence were used to construct the computational flow model for the stenosed phantom, whereas separately acquired T1weighted images were used for the normal phantom. The flow domain was subdivided into approximately 42,000 hexahedral elements for the normal carotid and 48,000 for the stenosed carotid bifurcation model, using an in-house mesh generation programme. The grid distribution was non-uniform to allow a finer grid in the bifurcation and near-wall regions. The length of the computational domain was 102 mm for the stenosed phantom, which was more than 12 times the common carotid diameter. The Navier–Stokes equations for 3D time-dependent flow with rigid walls were solved using the CFX4 finite-volume-based CFD solver. In order to
ARTICLE IN PRESS I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
match the experimental flow conditions, the 3D-PC velocity images acquired at the phantom inlet and outlet planes were applied as boundary conditions. Only the axial velocity components were used, since the secondary velocity components had amplitudes in the velocity noise region, which we estimated at 5% of venc : For the upstream boundary conditions, measured velocities at the model inlet were processed and smoothed before being mapped onto the computational grid inlet. Linear interpolation was performed in space and time between
681
the measurement points. For the downstream boundary conditions, zero normal gradients were assumed at both outlets, together with a time-varying flow division ratio derived from the MR velocity data. The computations were performed using the third-order accurate quadratic upwind (QUICK) differencing scheme in space and fully implicit backward Euler differencing in time. A fixed computational time step (8.29 ms) was set. Three cycles were computed to reach a cyclically repeated solution.
3. Results
22 20 18
flow (ml/s)
16 14 12 10 8 6 4 2 0 2
t1
t2
4
6
8
10
12
14
16
18
frame Fig. 1. Flow waveform, as measured by MRI in the CCA of the stenosed phantom. The arrows labelled t1 and t2 indicate the time frames studied in Figs. 3–6.
The flow waveform is shown in Fig. 1. It was measured from the axial velocity component in the common carotid vessel of the stenosed phantom, and has a mean of 7.4 ml/s, compared with a pump setting of 7.2 ml/s. Time points t1 (peak flow) and t2 (decelerating flow) are studied in more detail in Figs. 3–7. Fig. 2 shows maximum intensity projection (MIP) images generated from the phase contrast ‘‘reference’’ images, i.e. images without velocity sensitisation. An inverted grey scale is used for clarity, with fast flow showing as dark grey and slow flow showing as pale grey. Four of the 18 time frames are shown, together with the true phantom geometry constructed from data files supplied by the manufacturer. MR measured velocity vectors in the plane of symmetry are shown in Fig. 3, and the corresponding
Fig. 2. Time-resolved MIP images for the stenosed (upper row) and normal (lower row) carotid bifurcation models. Left to right: during acceleration; at peak flow; during deceleration; at minimum flow; true geometry. An inverted grey scale is used for the MIPs, with fast flow showing as dark grey and slow flow showing as pale grey. Note the regions of slow flow adjacent to the inner wall of the ICA in the stenosed model (arrow) and adjacent to the outer wall of the ICA in the normal model (arrowhead).
ARTICLE IN PRESS 682
I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
Fig. 3. Velocity vectors measured by MRI in the plane of symmetry of the stenosed phantom (upper row) and the normal phantom (lower row), at peak flow (t1 ) (left); during deceleration (t2 ) (right). Note the regions of slow flow adjacent to the inner wall of the ICA in the stenosed model (arrow) and adjacent to the outer wall of the ICA in the normal model (arrowhead). The true wall geometry has been superimposed for clarity.
CFD predictions in Fig. 4. Regions of reversed flow, predicted by CFD, are shown separately in Fig. 5. MR-derived WSS vectors are shown in Fig. 6, and the corresponding CFD predictions in Fig. 7. CFD-derived, time-averaged WSS values along the outer and inner walls of the ICA are shown in Fig. 8(a) and (b), respectively. Calculated oscillatory shear index (OSI; defined as Abkwd =ðAfwd þ Abkwd Þ; where Afwd and Abkwd are the areas under the WSS-time curve for forward and backward WSS, respectively (Suzuki et al., 1998)) is displayed in Fig. 9.
4. Discussion Flow in the symmetrically stenosed phantom is shifted towards the outer (non-divider) wall of the ICA. There is a small flow separation zone immediately downstream of the stenosis on the outer wall of the ICA, and a larger one further downstream and on the inner wall (arrow, Figs. 2 and 3). This is quite different from flow in the normal phantom, in which most of the ICA flow follows the inner wall, and there is a single flow separation zone in
ARTICLE IN PRESS I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
683
Fig. 4. Velocity vectors predicted by CFD. Stenosed phantom (left) and normal phantom (right), at t1 (upper row) and t2 (lower row).
the outer (bulb) region (arrowhead, Figs. 2 and 3). These effects are clearly seen in both the experimental (Figs. 2 and 3) and calculated (Figs. 4 and 5) results. ECA flow is similar for the two phantoms. For both phantoms, reversed flow is greater and more extensive during deceleration than at peak flow. The predicted WSS results (Fig. 7) show high values in the throat of the stenosis and in the ECA. Low values are seen downstream of the stenosis, especially on the inner wall of the ICA. Estimation of WSS vectors from experimental data (Fig. 6) requires wall identification and modelling of the velocity field from which the spatial derivatives can be calculated. Image segmentation, and hence wall determination, has broken down in the flow reversal region along the divider wall downstream of the stenosis, and hence WSS cannot be calculated there. This occurs because of the low image intensities caused by magnetic saturation of the slowly
moving fluid. Spurious flaps of ‘‘wall’’ occur along the outer ICA. The sharp flow divider has not been completely resolved, and the high WSS values predicted by CFD around there are missing from the MR results. There is reasonable agreement between the MRI and CFD results for the CCA and ECA. The plots of time-averaged WSS (Fig. 8) clearly bring out the differences between the two phantoms. The normal phantom has very low WSS magnitude along the outer wall of the ICA (Fig. 8(a)), extending from just before the bifurcation for more than 2 cm. In contrast, the stenosed phantom has a high maximum value (4 N/ m2) almost at the bifurcation, corresponding to the throat of the stenosis. Along the inner wall of the ICA (Fig. 8(b)), mean WSS in the stenosed phantom is very low (o0.5 N/m2) from just beyond the bifurcation for nearly 3 cm, corresponding to the flow reversal region. In this region, the WSS values for the normal phantom
ARTICLE IN PRESS 684
I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
Fig. 5. Regions of reversed flow predicted by CFD. Left to right; stenosed phantom at t1; normal phantom at t1; stenosed phantom at t2; normal phantom at t2.
Fig. 6. WSS vectors estimated from MR imaging of the stenosed phantom, shown superimposed on the segmented wall: (a) at time t1 ; and (b) at time t2 : WSS values above 3 N/m2 are coloured red.
ARTICLE IN PRESS I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
685
Fig. 7. WSS vectors predicted by CFD. Details as for Fig. 6.
are much higher. The flow reversal regions also show up well in the plots of OSI (Fig. 9). For the normal phantom, OSI values in the range 0.2–0.56 are present near the outer wall of the carotid bulb. On the other hand, higher values of OSI (up to 0.85) are found along the inner wall of the ICA of the stenosed phantom. Regions of low WSS and high OSI are thought to be susceptible to intimal thickening and plaque formation. Results for the normal phantom suggest that wall thickening might occur on the outer wall of the ICA. If this happened in a very mild form, it would reduce the size of the bulb and consequently the size of the flow separation zone would be reduced, WSS magnitude would be increased and OSI reduced (haemodynamically beneficial). However, this prediction is not relevant to the stenosed case, as the phantom studied has a 30% symmetrical stenosis, which causes narrowing on both the outer and inner walls of the ICA. This geometrical alteration can give rise to a series of haemodynamic changes, such as a flow separation zone along the inner wall (hence low WSS) as well as high followed by slow and reversed flow along the outer wall (spatially oscillating WSS). Our results confirm those in the literature, and extend the literature by providing qualitative comparisons between MRI and CFD for both normal and stenosed cases. Palmen et al. (1994) used a hydrogen bubble technique to visualise velocity profiles in large-scale models of healthy and mildly stenosed (25% area reduction) bifurcations. The main features observed in the carotid bulb were a shear layer that moved toward the divider wall during the deceleration phase and as the Womersley number was increased. Along the outer wall of the carotid bulb was a region of small negative velocities. These findings were confirmed by the same group (Gijsen et al., 1996) in a similar study using laser
Doppler anemometry, and are the accepted features of flow at the carotid bifurcation. These features are also found in our own experimental and numerical results. The stenosis used by Palmen et al. (1994) and Gijsen et al. (1996) was placed in a low shear region, and so had very little influence on the flow patterns. This is quite different from the effect caused by our more realistic stenosis, which shifted the predominant flow towards the outer wall in the carotid bulb region. Botnar et al. (2000) made MRI measurements of velocity in a model of a carotid bifurcation, and compared the results with the predictions of CFD. High flow occurred along the divider wall, with lower values along the outer wall, in agreement with previous work and the results reported here. They found that secondary flow developed during mid-systole, and remained until early diastole. There was good agreement of axial velocities between MRI and CFD, but not so good for the secondary flow. The maximum secondary velocities were about one-third of the axial velocities, and difficult to measure with MRI. We also had difficulty with measuring low velocities, and hence our WSS estimates are unreliable in some regions. Direct estimation of WSS from imaging remains technically challenging, with errors up to 40% expected from pixellimited data (Moore et al., 1998). Sub-pixel estimation has been used successfully in measurements of axial velocity components (Oyre et al., 1998). The combination of imaging and CFD is a powerful method that holds promise for research studies. Flow patterns are known to be highly dependent on the exact geometry. The common and branch vessels of our phantoms lay in a single plane, and so the model geometries were simpler than typical human arterial geometry. Ideally, an in vivo imaging method would be able to cope with complex geometry and extreme
ARTICLE IN PRESS 686
I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
Fig. 9. CFD predictions of OSI for the normal (left) and stenosed (right) phantoms.
Acknowledgements
Fig. 8. CFD predictions of the time-averaged WSS magnitude along the (a) outer and (b) inner wall of the ICA of the normal and stenosed phantoms. The apex of the bifurcation is at an axial distance of 0.0.
stenoses such as studied in a 2D computational model by Stroud et al. (2002). The presence of turbulent flow will, however, make MR measurement difficult. The most significant drawback of the MR acquisition method used in this work is the extended imaging time. Faster imaging sequences are clearly needed before vector WSS determination can become part of routine clinical practice.
This work was carried out at the SHEFC Brain Imaging Research Centre for Scotland, Edinburgh, and was funded in part by the EPSRC. The viscosity measurements were kindly made by Donald Easton of the Department of Haematology, Royal Infirmary of . Edinburgh. Dr. Uwe Kohler developed the segmentation and WSS calculation software. Quan Long developed software for computational grid generation. Martin Connell provided support for computing and visualisation. Bob Gravett (Shelley Medical Imaging Technologies) provided true geometry data for the phantoms.
References 5. Conclusions We have demonstrated MRI measurement of fluid flow and subsequent estimation of WSS vectors during pulsatile flow in normal and stenosed carotid phantoms. The estimates were in qualitative agreement with CFD predictions except immediately downstream of the stenosis, where image segmentation was unsatisfactory.
Berger, S.A., Jou, L.-D., 2000. Flows in stenotic vessels. Annual Review of Fluid Mechanics 32, 347–382. Botnar, R., Rappitsch, G., Scheidegger, M.B., Liepsch, D., Perktold, K., Boesiger, P., 2000. Hemodynamics in the carotid artery bifurcation: a comparison between numerical simulations and in vitro MRI measurements. Journal of Biomechanics 33, 137–144. Brands, P.J., Hoeks, A.P.G., Hofstra, L., Reneman, R.S., 1995. A noninvasive method to estimate wall shear rate using ultrasound. Ultrasound in Medicine and Biology 21, 171–185.
ARTICLE IN PRESS I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687 Caro, C.G., Fitzgerald, J.M., Schroter, R.C., 1971. Atheroma and arterial wall shear: observation, correlation and proposal of a shear-dependent mass transfer mechanism for atherogenesis. Proceedings of the Royal Society of London 117, 109–159. Dunmire, B., Beach, K.W., Labs, K.H., Plett, M., Strandness, D.E., 2000. Cross-beam vector Doppler ultrasound for angle-independent velocity measurements. Ultrasound in Medicine and Biology 26, 1213–1235. Giddens, D.P., Zarins, C.K., Glagov, S., 1993. The role of fluid mechanics in the localisation and detection of atherosclerosis. Journal of Biomechanical Engineering 115, 588–594. Gijsen, F.J.H., Palmen, D.E.M., van der Beek, M.H.E., van de Vosse, F.N., van Dongen, M.E.H., Janssen, J.D., 1996. Analysis of the axial flow field in stenosed carotid artery bifurcation models— LDA experiments. Journal of Biomechanics 29, 1483–1489. Glagov, S., Zarins, C.K., Giddens, D.P., Ku, D.N., 1988. Haemodynamics and atherosclerosis—insights and perspectives gained from studies of human arteries. Archives of Pathology and Laboratory Medicine 112, 1018–1031. Gnasso, A., Irace, C., Carallo, C., De Franceschi, M.S., Motti, C., Mattioli, P.L., Pujia, A., 1997. In vivo association between low wall shear stress and plaque in subjects with asymmetrical carotid atherosclerosis. Stroke 28, 993–998. Hoeks, A.P.G., Samijo, S.K., Brands, P.J., Reneman, R.S., 1995. Noninvasive determination of shear-rate distribution across the arterial lumen. Hypertension 26, 26–33. . Kohler, U., Marshall, I., Robertson, M.B., Long, Q., Xu, Y.X., Hoskins, P.R., 2001. MRI measurement of wall shear stress vectors in bifurcation models and comparison with CFD predictions. Journal of Magnetic Resonance Imaging 14, 563–573. Ku, D.N., 1997. Blood flow in arteries. Annual Review of Fluid Mechanics 29, 399–434. Ku, D.N., Giddens, D.P., Zarins, C.K., Glagov, S., 1985. Pulsatile flow and atherosclerosis in the human carotid bifurcation—positive correlation between plaque location and low and oscillating shear stress. Atherosclerosis 5, 293–302. Light, E.D., Davidsen, R.E., Fiering, J.O., Hruschka, T.A., Smith, S.W., 1998. Progress in two-dimensional arrays for real-time volumetric imaging. Ultrasonic Imaging 20, 1–15. Long, Q., Xu, Y.X., Ariff, B., Thom, S.A., Hughes, A.D., Stanton, A.V., 2000. Reconstruction of blood flow patterns in a human carotid bifurcation: a combined CFD and MRI study. Journal of Magnetic Resonance Imaging 11, 299–311. Malek, A.M., Alper, S.L., Izumo, S., 1999. Hemodynamic shear stress and its role in atherosclerosis. Journal of the American Medical Association 282, 2035–2042.
687
Moore, J.A., Steinman, D.A., Ethier, C.R., 1998. Computational blood flow modelling: errors associated with reconstructing finite element models from magnetic resonance images. Journal of Biomechanics 31, 179–184. Oshinski, J.N., Ku, D.N., Mukundan, S., Loth, F., Pettigrew, R.I., 1995. Determination of wall shear stress in the aorta with the use of MR phase velocity mapping. Journal of Magnetic Resonance Imaging 5, 640–647. Oyre, S., Ringgaard, S., Kozerke, S., Paaske, W.P., Scheidegger, M.B., Boesiger, P., Pedersen, E.M., 1998. Quantitation of circumferential subpixel vessel wall position and wall shear stress by multiple sectored three-dimensional paraboloid modeling of velocity encoded cine MR. Magnetic Resonance in Medicine 40, 645–655. Palmen, D.E.M., van de Vosse, F.N., Janssen, J.D., van Dongen, M.E.H., 1994. Analysis of the flow in stenosed carotid artery bifurcation models—hydrogen bubble visualisation. Journal of Biomechanics 27, 581–590. . Papathanasopoulou, P., Zhao, S., Kohler, U., Robertson, M.B., Long, Q., Hoskins, P., Xu, X.Y., Marshall, I., 2003. MRI measurement of time-resolved wall shear stress vectors in a carotid bifurcation model, and comparison with CFD predictions. Journal of Magnetic Resonance Imaging 17, 153–162. Steinman, D.A., Thomas, J.B., Ladak, H.M., Milner, J.S., Rutt, B.K., Spence, J.D., 2002. Reconstruction of carotid bifurcation hemodynamics and wall thickness using computational fluid dynamics and MRI. Magnetic Resonance in Medicine 47, 149–159. Stokholm, R., Oyre, S., Ringgaard, S., Flaagoy, H., Paaske, W.P., Pedersen, E.M., 2000. Determination of wall shear rate in the human carotid artery by magnetic resonance techniques. European Journal of Vascular and Endovascular Surgery 20, 427–433. Stroud, J.S., Berger, S.A., Saloner, D., 2002. Numerical analysis of flow through a severely stenotic carotid artery bifurcation. ASME Journal of Biomechanical Engineering 124, 9–20. Suzuki, J., Shimamoto, R., Nishikawa, J., Tomaru, T., Nakajima, T., Nakamura, F., Shin, W.S., Toyo-oka, T., 1998. Vector analysis of the hemodynamics of atherogenesis in the human thoracic aorta using MR velocity mapping. American Journal of Roentgenology 171, 1285–1290. . . Wigstrom, L., Sjoqqvist, L., Wranne, B., 1996. Temporally resolved 3D phase-contrast imaging. Magnetic Resonance in Medicine 36, 800–803. Zhao, S.Z., Xu, Y.X., Hughes, A.D., Thom, S.A., Stanton, A.V., Ariff, B., Long, Q., 2000. Blood flow and vessel mechanics in a physiologically realistic model of a human carotid arterial bifurcation. Journal of Biomechanics 33, 975–984.
Journal of Biomechanics 37 (2004) 679–687
MRI and CFD studies of pulsatile flow in healthy and stenosed carotid bifurcation models Ian Marshalla,*, Shunzhi Zhaob, Panorea Papathanasopouloua, Peter Hoskinsa, X Yun Xub a
Department of Medical and Radiological Sciences, Medical Physics, University of Edinburgh, Western General Hospital, Edinburgh EH4 2XU, UK b Department of Chemical Engineering and Chemical Technology, Imperial College, London, UK Accepted 18 September 2003
Abstract Pulsatile flow was studied in physiologically realistic models of a normal and a moderately stenosed (30% diameter reduction) human carotid bifurcation. Time-resolved velocity measurements were made using magnetic resonance imaging, from which wall shear stress (WSS) vectors were calculated. Velocity measurements in the inflow and outflow regions were also used as boundary conditions for a computational fluid dynamics (CFD) model. Experimental flow patterns and derived WSS vectors were compared qualitatively with the corresponding CFD predictions. In the stenosed phantom, flow in the bulb region of the ‘‘internal carotid artery’’ was concentrated along the outer wall, with a region of low and recirculating flow near the inner wall. In the normal phantom, the converse was found, with a low flow region near the outer wall of the bulb. Time-averaged WSS and oscillatory shear index were also markedly different for the two phantoms. r 2003 Elsevier Ltd. All rights reserved. Keywords: Carotid blood flow; Velocity measurement; Wall shear stress; CFD; Stenosis
1. Introduction Over the last few decades, research has demonstrated the influence of vessel geometry and haemodynamic forces on the development of vascular pathology (Ku et al., 1985; Glagov et al., 1988). In regions of disturbed blood flow, which are found around vessel bifurcations, there is an increased chance of atheroma deposition. The viscous drag of the blood on the vessel wall (wall shear stress; WSS) is believed to play a part in regulating arterial structure (Caro et al., 1971; Giddens et al., 1993; Gnasso et al., 1997; Malek et al., 1999), and is implicated in plaque development. The exact relationship between atheroma deposition and WSS remains uncertain, with low WSS, oscillating WSS, and WSS gradient all having been considered (Ku et al., 1985; Caro et al., 1971; Giddens et al., 1993). There is considerable interest in the development of in vivo techniques that would allow accurate blood flow quantification and estimation of WSS, but they are not *Corresponding author. Tel.: +44-131-537-1661; fax: +44-131-5371026. E-mail address: [email protected] (I. Marshall). 0021-9290/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2003.09.032
yet available for clinical use. The two most promising techniques are Doppler ultrasound and magnetic resonance imaging (MRI). Doppler ultrasound is a real time technique and does not suffer from data loss during turbulent flow. Current commercial systems are based on acquisition of two-dimensional (2D) image data in which only a single component of velocity is acquired. Acquisition of 2 or 3 velocity components requires multi-beam systems which are the subject of current research (Dunmire et al., 2000). Acquisition of real time three-dimensional (3D) information is also possible using 2D array systems (Light et al., 1998), but these are at an early stage of development. Ultrasound cannot be used on arteries that lie behind bony or air-filled structures. MRI is intrinsically 3D, but suffers from long acquisition times and signal loss in turbulent flow, which limits its use in diseased arteries. Both techniques have spatial resolutions that are only just high enough (0.5 mm) for in vivo studies of carotid arteries. WSS estimation methods based on MRI have to date been mainly restricted to measurement of the axial component at a single plane (Oshinski et al., 1995; Oyre et al., 1998; Stokholm et al., 2000). Similarly, ultrasound methods have been based on estimation of axial WSS at
ARTICLE IN PRESS 680
I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
a single location (Brands et al., 1995; Hoeks et al., 1995). Vector MRI information was obtained at the intersection of selected planes in the aorta by Suzuki et al. (1998). Time-resolved measurement of in vivo flow and WSS vectors remains a challenging goal. Therefore, the study of fluid dynamics is often carried out in physical or computational models (Zhao et al., 2000). For useful reviews, see Ku (1997) and Berger and Jou (2000). There have been very few experimental studies of pulsatile flow in carotid models, especially in stenosed models. Palmen et al. (1994) used hydrogen bubble visualisation, Gijsen et al. (1996) used laser Doppler anemometry, and Botnar et al. (2000) used MRI. We have previously reported on the estimation of WSS vectors throughout a flow region from MRI velocity measurements of steady . (Kohler et al., 2001) and pulsatile (Papathanasopoulou et al., 2003) flow in phantoms. An alternative approach has been to use MRI structural data and velocity measurements as the boundary conditions to solve the Navier–Stokes equations in a computational fluid dynamics (CFD) package, from which WSS can be calculated (Long et al., 2000; Steinman et al., 2002). In the present work, we used a time-resolved 3D phase-contrast (PC) MRI sequence to acquire all three velocity components of pulsatile flow in models of a human carotid bifurcation with and without a 30% (diameter) axisymmetrical stenosis. We show velocity vectors at selected time frames, and WSS vectors calculated from the MRI measurements. These are compared qualitatively with the predictions of CFD calculations. The comparison of normal and stenosed phantoms under physiologically realistic flow conditions is novel, as is the qualitative comparison of MRI and CFD results for the stenosed phantom.
(vx ; vy ; vz ) to obtain a temporal resolution of 50 ms. We . implemented a retrospectively gated method (Wigstrom et al., 1996; Papathanasopoulou et al., 2003). The sequence was run with 18 repetitions to ensure that each period (‘‘RR-interval’’) of the pulsatile waveform was adequately sampled for each line of k-space. The pixel size was 0.5 mm, and 48 slices of 1.4 mm thickness were used to cover the bifurcation region. The encoding velocity (venc ) was 120 cm/s. Raw k-space data files were saved for offline reconstruction. Trigger pulses from the pump were logged by a timing system (IFIS; Psychology Software Tools, Pittsburgh, PA) which was synchronised to the scanning. Acquisition time was approximately 90 min for each velocity component. 2.2. Image reconstruction and processing
2. Methods
The raw k-space data from the scanner and the data from the timing system were transferred to a UNIX workstation for reconstruction. Data were sorted into the 18 desired time frames by linear interpolation between the two nearest available samples. k-space data were zero-filled to 64 slices, and a 3D FFT applied to generate images with 1.05(=1.4 48/64) mm effective slice thickness. Finally, velocity images were generated by phase subtraction of the reference images from the corresponding velocity-encoded images. Reconstructed magnitude and velocity images for each time frame were processed as described previously, using software devel. oped in-house (Kohler et al., 2001). Flow regions were automatically segmented, and the velocity field fitted by a fifth order polynomial model. The downstream branches (‘‘internal carotid artery’’ (ICA) and ‘‘external carotid artery’’ (ECA)) were fitted separately from the ‘‘common carotid artery’’ (CCA). Vessel walls were determined from the segmentation results, and WSS vectors were calculated from the spatial derivatives of the fitted velocity model.
2.1. Data acquisition
2.3. CFD
MR scanning was performed on a 1.5 T GE Signa scanner (GE Medical, Milwaukee), using a 75 mm surface coil. The flow phantoms were anthropomorphically realistic models of a human carotid bifurcation with an axisymmetric 30% diameter (50% area) stenosis and a healthy carotid bifurcation (Shelley Medical Imaging Technologies, London, Ont., Canada). A computer-controlled pump (UHDC Flow System, Shelley) was used to generate a realistic carotid waveform, and pumped blood-mimicking fluid (Shelley) through the phantoms. The fluid viscosity was measured with a capillary viscometer and determined to be 3.4 103 N s/m2. Velocity measurements were made using a 3D-PC sequence, run separately for each velocity direction
The ‘‘magnitude’’ images from the PC sequence were used to construct the computational flow model for the stenosed phantom, whereas separately acquired T1weighted images were used for the normal phantom. The flow domain was subdivided into approximately 42,000 hexahedral elements for the normal carotid and 48,000 for the stenosed carotid bifurcation model, using an in-house mesh generation programme. The grid distribution was non-uniform to allow a finer grid in the bifurcation and near-wall regions. The length of the computational domain was 102 mm for the stenosed phantom, which was more than 12 times the common carotid diameter. The Navier–Stokes equations for 3D time-dependent flow with rigid walls were solved using the CFX4 finite-volume-based CFD solver. In order to
ARTICLE IN PRESS I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
match the experimental flow conditions, the 3D-PC velocity images acquired at the phantom inlet and outlet planes were applied as boundary conditions. Only the axial velocity components were used, since the secondary velocity components had amplitudes in the velocity noise region, which we estimated at 5% of venc : For the upstream boundary conditions, measured velocities at the model inlet were processed and smoothed before being mapped onto the computational grid inlet. Linear interpolation was performed in space and time between
681
the measurement points. For the downstream boundary conditions, zero normal gradients were assumed at both outlets, together with a time-varying flow division ratio derived from the MR velocity data. The computations were performed using the third-order accurate quadratic upwind (QUICK) differencing scheme in space and fully implicit backward Euler differencing in time. A fixed computational time step (8.29 ms) was set. Three cycles were computed to reach a cyclically repeated solution.
3. Results
22 20 18
flow (ml/s)
16 14 12 10 8 6 4 2 0 2
t1
t2
4
6
8
10
12
14
16
18
frame Fig. 1. Flow waveform, as measured by MRI in the CCA of the stenosed phantom. The arrows labelled t1 and t2 indicate the time frames studied in Figs. 3–6.
The flow waveform is shown in Fig. 1. It was measured from the axial velocity component in the common carotid vessel of the stenosed phantom, and has a mean of 7.4 ml/s, compared with a pump setting of 7.2 ml/s. Time points t1 (peak flow) and t2 (decelerating flow) are studied in more detail in Figs. 3–7. Fig. 2 shows maximum intensity projection (MIP) images generated from the phase contrast ‘‘reference’’ images, i.e. images without velocity sensitisation. An inverted grey scale is used for clarity, with fast flow showing as dark grey and slow flow showing as pale grey. Four of the 18 time frames are shown, together with the true phantom geometry constructed from data files supplied by the manufacturer. MR measured velocity vectors in the plane of symmetry are shown in Fig. 3, and the corresponding
Fig. 2. Time-resolved MIP images for the stenosed (upper row) and normal (lower row) carotid bifurcation models. Left to right: during acceleration; at peak flow; during deceleration; at minimum flow; true geometry. An inverted grey scale is used for the MIPs, with fast flow showing as dark grey and slow flow showing as pale grey. Note the regions of slow flow adjacent to the inner wall of the ICA in the stenosed model (arrow) and adjacent to the outer wall of the ICA in the normal model (arrowhead).
ARTICLE IN PRESS 682
I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
Fig. 3. Velocity vectors measured by MRI in the plane of symmetry of the stenosed phantom (upper row) and the normal phantom (lower row), at peak flow (t1 ) (left); during deceleration (t2 ) (right). Note the regions of slow flow adjacent to the inner wall of the ICA in the stenosed model (arrow) and adjacent to the outer wall of the ICA in the normal model (arrowhead). The true wall geometry has been superimposed for clarity.
CFD predictions in Fig. 4. Regions of reversed flow, predicted by CFD, are shown separately in Fig. 5. MR-derived WSS vectors are shown in Fig. 6, and the corresponding CFD predictions in Fig. 7. CFD-derived, time-averaged WSS values along the outer and inner walls of the ICA are shown in Fig. 8(a) and (b), respectively. Calculated oscillatory shear index (OSI; defined as Abkwd =ðAfwd þ Abkwd Þ; where Afwd and Abkwd are the areas under the WSS-time curve for forward and backward WSS, respectively (Suzuki et al., 1998)) is displayed in Fig. 9.
4. Discussion Flow in the symmetrically stenosed phantom is shifted towards the outer (non-divider) wall of the ICA. There is a small flow separation zone immediately downstream of the stenosis on the outer wall of the ICA, and a larger one further downstream and on the inner wall (arrow, Figs. 2 and 3). This is quite different from flow in the normal phantom, in which most of the ICA flow follows the inner wall, and there is a single flow separation zone in
ARTICLE IN PRESS I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
683
Fig. 4. Velocity vectors predicted by CFD. Stenosed phantom (left) and normal phantom (right), at t1 (upper row) and t2 (lower row).
the outer (bulb) region (arrowhead, Figs. 2 and 3). These effects are clearly seen in both the experimental (Figs. 2 and 3) and calculated (Figs. 4 and 5) results. ECA flow is similar for the two phantoms. For both phantoms, reversed flow is greater and more extensive during deceleration than at peak flow. The predicted WSS results (Fig. 7) show high values in the throat of the stenosis and in the ECA. Low values are seen downstream of the stenosis, especially on the inner wall of the ICA. Estimation of WSS vectors from experimental data (Fig. 6) requires wall identification and modelling of the velocity field from which the spatial derivatives can be calculated. Image segmentation, and hence wall determination, has broken down in the flow reversal region along the divider wall downstream of the stenosis, and hence WSS cannot be calculated there. This occurs because of the low image intensities caused by magnetic saturation of the slowly
moving fluid. Spurious flaps of ‘‘wall’’ occur along the outer ICA. The sharp flow divider has not been completely resolved, and the high WSS values predicted by CFD around there are missing from the MR results. There is reasonable agreement between the MRI and CFD results for the CCA and ECA. The plots of time-averaged WSS (Fig. 8) clearly bring out the differences between the two phantoms. The normal phantom has very low WSS magnitude along the outer wall of the ICA (Fig. 8(a)), extending from just before the bifurcation for more than 2 cm. In contrast, the stenosed phantom has a high maximum value (4 N/ m2) almost at the bifurcation, corresponding to the throat of the stenosis. Along the inner wall of the ICA (Fig. 8(b)), mean WSS in the stenosed phantom is very low (o0.5 N/m2) from just beyond the bifurcation for nearly 3 cm, corresponding to the flow reversal region. In this region, the WSS values for the normal phantom
ARTICLE IN PRESS 684
I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
Fig. 5. Regions of reversed flow predicted by CFD. Left to right; stenosed phantom at t1; normal phantom at t1; stenosed phantom at t2; normal phantom at t2.
Fig. 6. WSS vectors estimated from MR imaging of the stenosed phantom, shown superimposed on the segmented wall: (a) at time t1 ; and (b) at time t2 : WSS values above 3 N/m2 are coloured red.
ARTICLE IN PRESS I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
685
Fig. 7. WSS vectors predicted by CFD. Details as for Fig. 6.
are much higher. The flow reversal regions also show up well in the plots of OSI (Fig. 9). For the normal phantom, OSI values in the range 0.2–0.56 are present near the outer wall of the carotid bulb. On the other hand, higher values of OSI (up to 0.85) are found along the inner wall of the ICA of the stenosed phantom. Regions of low WSS and high OSI are thought to be susceptible to intimal thickening and plaque formation. Results for the normal phantom suggest that wall thickening might occur on the outer wall of the ICA. If this happened in a very mild form, it would reduce the size of the bulb and consequently the size of the flow separation zone would be reduced, WSS magnitude would be increased and OSI reduced (haemodynamically beneficial). However, this prediction is not relevant to the stenosed case, as the phantom studied has a 30% symmetrical stenosis, which causes narrowing on both the outer and inner walls of the ICA. This geometrical alteration can give rise to a series of haemodynamic changes, such as a flow separation zone along the inner wall (hence low WSS) as well as high followed by slow and reversed flow along the outer wall (spatially oscillating WSS). Our results confirm those in the literature, and extend the literature by providing qualitative comparisons between MRI and CFD for both normal and stenosed cases. Palmen et al. (1994) used a hydrogen bubble technique to visualise velocity profiles in large-scale models of healthy and mildly stenosed (25% area reduction) bifurcations. The main features observed in the carotid bulb were a shear layer that moved toward the divider wall during the deceleration phase and as the Womersley number was increased. Along the outer wall of the carotid bulb was a region of small negative velocities. These findings were confirmed by the same group (Gijsen et al., 1996) in a similar study using laser
Doppler anemometry, and are the accepted features of flow at the carotid bifurcation. These features are also found in our own experimental and numerical results. The stenosis used by Palmen et al. (1994) and Gijsen et al. (1996) was placed in a low shear region, and so had very little influence on the flow patterns. This is quite different from the effect caused by our more realistic stenosis, which shifted the predominant flow towards the outer wall in the carotid bulb region. Botnar et al. (2000) made MRI measurements of velocity in a model of a carotid bifurcation, and compared the results with the predictions of CFD. High flow occurred along the divider wall, with lower values along the outer wall, in agreement with previous work and the results reported here. They found that secondary flow developed during mid-systole, and remained until early diastole. There was good agreement of axial velocities between MRI and CFD, but not so good for the secondary flow. The maximum secondary velocities were about one-third of the axial velocities, and difficult to measure with MRI. We also had difficulty with measuring low velocities, and hence our WSS estimates are unreliable in some regions. Direct estimation of WSS from imaging remains technically challenging, with errors up to 40% expected from pixellimited data (Moore et al., 1998). Sub-pixel estimation has been used successfully in measurements of axial velocity components (Oyre et al., 1998). The combination of imaging and CFD is a powerful method that holds promise for research studies. Flow patterns are known to be highly dependent on the exact geometry. The common and branch vessels of our phantoms lay in a single plane, and so the model geometries were simpler than typical human arterial geometry. Ideally, an in vivo imaging method would be able to cope with complex geometry and extreme
ARTICLE IN PRESS 686
I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687
Fig. 9. CFD predictions of OSI for the normal (left) and stenosed (right) phantoms.
Acknowledgements
Fig. 8. CFD predictions of the time-averaged WSS magnitude along the (a) outer and (b) inner wall of the ICA of the normal and stenosed phantoms. The apex of the bifurcation is at an axial distance of 0.0.
stenoses such as studied in a 2D computational model by Stroud et al. (2002). The presence of turbulent flow will, however, make MR measurement difficult. The most significant drawback of the MR acquisition method used in this work is the extended imaging time. Faster imaging sequences are clearly needed before vector WSS determination can become part of routine clinical practice.
This work was carried out at the SHEFC Brain Imaging Research Centre for Scotland, Edinburgh, and was funded in part by the EPSRC. The viscosity measurements were kindly made by Donald Easton of the Department of Haematology, Royal Infirmary of . Edinburgh. Dr. Uwe Kohler developed the segmentation and WSS calculation software. Quan Long developed software for computational grid generation. Martin Connell provided support for computing and visualisation. Bob Gravett (Shelley Medical Imaging Technologies) provided true geometry data for the phantoms.
References 5. Conclusions We have demonstrated MRI measurement of fluid flow and subsequent estimation of WSS vectors during pulsatile flow in normal and stenosed carotid phantoms. The estimates were in qualitative agreement with CFD predictions except immediately downstream of the stenosis, where image segmentation was unsatisfactory.
Berger, S.A., Jou, L.-D., 2000. Flows in stenotic vessels. Annual Review of Fluid Mechanics 32, 347–382. Botnar, R., Rappitsch, G., Scheidegger, M.B., Liepsch, D., Perktold, K., Boesiger, P., 2000. Hemodynamics in the carotid artery bifurcation: a comparison between numerical simulations and in vitro MRI measurements. Journal of Biomechanics 33, 137–144. Brands, P.J., Hoeks, A.P.G., Hofstra, L., Reneman, R.S., 1995. A noninvasive method to estimate wall shear rate using ultrasound. Ultrasound in Medicine and Biology 21, 171–185.
ARTICLE IN PRESS I. Marshall et al. / Journal of Biomechanics 37 (2004) 679–687 Caro, C.G., Fitzgerald, J.M., Schroter, R.C., 1971. Atheroma and arterial wall shear: observation, correlation and proposal of a shear-dependent mass transfer mechanism for atherogenesis. Proceedings of the Royal Society of London 117, 109–159. Dunmire, B., Beach, K.W., Labs, K.H., Plett, M., Strandness, D.E., 2000. Cross-beam vector Doppler ultrasound for angle-independent velocity measurements. Ultrasound in Medicine and Biology 26, 1213–1235. Giddens, D.P., Zarins, C.K., Glagov, S., 1993. The role of fluid mechanics in the localisation and detection of atherosclerosis. Journal of Biomechanical Engineering 115, 588–594. Gijsen, F.J.H., Palmen, D.E.M., van der Beek, M.H.E., van de Vosse, F.N., van Dongen, M.E.H., Janssen, J.D., 1996. Analysis of the axial flow field in stenosed carotid artery bifurcation models— LDA experiments. Journal of Biomechanics 29, 1483–1489. Glagov, S., Zarins, C.K., Giddens, D.P., Ku, D.N., 1988. Haemodynamics and atherosclerosis—insights and perspectives gained from studies of human arteries. Archives of Pathology and Laboratory Medicine 112, 1018–1031. Gnasso, A., Irace, C., Carallo, C., De Franceschi, M.S., Motti, C., Mattioli, P.L., Pujia, A., 1997. In vivo association between low wall shear stress and plaque in subjects with asymmetrical carotid atherosclerosis. Stroke 28, 993–998. Hoeks, A.P.G., Samijo, S.K., Brands, P.J., Reneman, R.S., 1995. Noninvasive determination of shear-rate distribution across the arterial lumen. Hypertension 26, 26–33. . Kohler, U., Marshall, I., Robertson, M.B., Long, Q., Xu, Y.X., Hoskins, P.R., 2001. MRI measurement of wall shear stress vectors in bifurcation models and comparison with CFD predictions. Journal of Magnetic Resonance Imaging 14, 563–573. Ku, D.N., 1997. Blood flow in arteries. Annual Review of Fluid Mechanics 29, 399–434. Ku, D.N., Giddens, D.P., Zarins, C.K., Glagov, S., 1985. Pulsatile flow and atherosclerosis in the human carotid bifurcation—positive correlation between plaque location and low and oscillating shear stress. Atherosclerosis 5, 293–302. Light, E.D., Davidsen, R.E., Fiering, J.O., Hruschka, T.A., Smith, S.W., 1998. Progress in two-dimensional arrays for real-time volumetric imaging. Ultrasonic Imaging 20, 1–15. Long, Q., Xu, Y.X., Ariff, B., Thom, S.A., Hughes, A.D., Stanton, A.V., 2000. Reconstruction of blood flow patterns in a human carotid bifurcation: a combined CFD and MRI study. Journal of Magnetic Resonance Imaging 11, 299–311. Malek, A.M., Alper, S.L., Izumo, S., 1999. Hemodynamic shear stress and its role in atherosclerosis. Journal of the American Medical Association 282, 2035–2042.
687
Moore, J.A., Steinman, D.A., Ethier, C.R., 1998. Computational blood flow modelling: errors associated with reconstructing finite element models from magnetic resonance images. Journal of Biomechanics 31, 179–184. Oshinski, J.N., Ku, D.N., Mukundan, S., Loth, F., Pettigrew, R.I., 1995. Determination of wall shear stress in the aorta with the use of MR phase velocity mapping. Journal of Magnetic Resonance Imaging 5, 640–647. Oyre, S., Ringgaard, S., Kozerke, S., Paaske, W.P., Scheidegger, M.B., Boesiger, P., Pedersen, E.M., 1998. Quantitation of circumferential subpixel vessel wall position and wall shear stress by multiple sectored three-dimensional paraboloid modeling of velocity encoded cine MR. Magnetic Resonance in Medicine 40, 645–655. Palmen, D.E.M., van de Vosse, F.N., Janssen, J.D., van Dongen, M.E.H., 1994. Analysis of the flow in stenosed carotid artery bifurcation models—hydrogen bubble visualisation. Journal of Biomechanics 27, 581–590. . Papathanasopoulou, P., Zhao, S., Kohler, U., Robertson, M.B., Long, Q., Hoskins, P., Xu, X.Y., Marshall, I., 2003. MRI measurement of time-resolved wall shear stress vectors in a carotid bifurcation model, and comparison with CFD predictions. Journal of Magnetic Resonance Imaging 17, 153–162. Steinman, D.A., Thomas, J.B., Ladak, H.M., Milner, J.S., Rutt, B.K., Spence, J.D., 2002. Reconstruction of carotid bifurcation hemodynamics and wall thickness using computational fluid dynamics and MRI. Magnetic Resonance in Medicine 47, 149–159. Stokholm, R., Oyre, S., Ringgaard, S., Flaagoy, H., Paaske, W.P., Pedersen, E.M., 2000. Determination of wall shear rate in the human carotid artery by magnetic resonance techniques. European Journal of Vascular and Endovascular Surgery 20, 427–433. Stroud, J.S., Berger, S.A., Saloner, D., 2002. Numerical analysis of flow through a severely stenotic carotid artery bifurcation. ASME Journal of Biomechanical Engineering 124, 9–20. Suzuki, J., Shimamoto, R., Nishikawa, J., Tomaru, T., Nakajima, T., Nakamura, F., Shin, W.S., Toyo-oka, T., 1998. Vector analysis of the hemodynamics of atherogenesis in the human thoracic aorta using MR velocity mapping. American Journal of Roentgenology 171, 1285–1290. . . Wigstrom, L., Sjoqqvist, L., Wranne, B., 1996. Temporally resolved 3D phase-contrast imaging. Magnetic Resonance in Medicine 36, 800–803. Zhao, S.Z., Xu, Y.X., Hughes, A.D., Thom, S.A., Stanton, A.V., Ariff, B., Long, Q., 2000. Blood flow and vessel mechanics in a physiologically realistic model of a human carotid arterial bifurcation. Journal of Biomechanics 33, 975–984.