Assessment of Carotid Flow Using Magnetic Resonance Imaging and Computational Fluid Dynamics

Knowledge of blood flowpatterns in the human body is a critical component in cardiovascular disease research and diagnosis. Carotid atherosclerosis, for example, refers to the narrowing of the carotid arteries. One symptom of atherosclerosis is abnormal flow. The carotid arteries supply blood to the brain, so early detection of carotid stenosis may prevent thrombotic stroke. The current clinical gold standard for cardiovascular flow measurement is Doppler ultrasound. However, evaluation by ultrasound is inadequate when there is fat, air, bone, or surgical scar in the acoustic path, and flowmeasurement is inaccurate when the ultrasound beam cannot be properly aligned with the axis of flow (Hoskins, 1996; Winkler & Wu, 1995). Two alternative approaches to 3D flow assessment are currently available to the researcher and clinician: (i) direct, model-independent velocity mapping using flow-encoded magnetic resonance imaging (MRI), and (ii) model-based computational fluid dynamics (CFD) calculations. MRI is potentially the most appropriate technique for addressing all aspects of cardiovascular disease examination. MRI overcomes the acoustical window limitations of ultrasound, potentially allowing flowmeasurements to be obtained along any direction, and for any vessel in the cardiovascular system. MRI measurements are also less operator-dependent than those of Doppler ultrasound, and the true direction of flow can generally be precisely measured. The current gold standard for MRI flow quantification is phase contrast (PC) (O’Donnell, 1985). However, PC suffers from partial-volume effects when a wide distribution of velocities is contained within a single voxel (Tang et al., 1993). This is particularly problematic when flow is turbulent and/or complex (e.g., flow jets due to stenosis) or at the interface between blood and vessel wall (viscous sublayer). This issue is typically addressed by increasing the spatial resolution, which dramatically affects the signal-to-noise ratio (SNR) and increases the scan time. Therefore, PC may be inadequate for estimating the peak velocity of stenotic flow jets and for assessing wall shear rate. 23


Introduction
Knowledge of blood flow patterns in the human body is a critical component in cardiovascular disease research and diagnosis. Carotid atherosclerosis, for example, refers to the narrowing of the carotid arteries. One symptom of atherosclerosis is abnormal flow. The carotid arteries supply blood to the brain, so early detection of carotid stenosis may prevent thrombotic stroke. The current clinical gold standard for cardiovascular flow measurement is Doppler ultrasound. However, evaluation by ultrasound is inadequate when there is fat, air, bone, or surgical scar in the acoustic path, and flow measurement is inaccurate when the ultrasound beam cannot be properly aligned with the axis of flow (Hoskins, 1996;Winkler & Wu, 1995). Two alternative approaches to 3D flow assessment are currently available to the researcher and clinician: (i) direct, model-independent velocity mapping using flow-encoded magnetic resonance imaging (MRI), and (ii) model-based computational fluid dynamics (CFD) calculations. MRI is potentially the most appropriate technique for addressing all aspects of cardiovascular disease examination. MRI overcomes the acoustical window limitations of ultrasound, potentially allowing flow measurements to be obtained along any direction, and for any vessel in the cardiovascular system. MRI measurements are also less operator-dependent than those of Doppler ultrasound, and the true direction of flow can generally be precisely measured. The current gold standard for MRI flow quantification is phase contrast (PC) (O'Donnell, 1985). However, PC suffers from partial-volume effects when a wide distribution of velocities is contained within a single voxel (Tang et al., 1993). This is particularly problematic when flow is turbulent and/or complex (e.g., flow jets due to stenosis) or at the interface between blood and vessel wall (viscous sublayer). This issue is typically addressed by increasing the spatial resolution, which dramatically affects the signal-to-noise ratio (SNR) and increases the scan time. Therefore, PC may be inadequate for estimating the peak velocity of stenotic flow jets and for assessing wall shear rate.
An alternative approach for cardiovascular flow assessment is to reconstruct a complex flow field via CFD simulation, using vascular geometries and input/output functions derived from non-invasive imaging data (Kim et al., 2006;Papathanasopoulou et al., 2003;Steinman et al., 2002). The accuracy of conventional CFD routines hinges on many modeling assumptions that are not strictly true for in vivo vascular flow, including rigid vessel walls and uniform blood viscosity. Furthermore, conventional CFD routines typically employ a structured non-Cartesian finite-element mesh that conforms to the vessel geometry, which leads to relatively complex algorithms that incur significant computational costs. This chapter proposes two new approaches for MRI assessment of carotid flow. First, we introduce a rapid MRI method for fully-localized measurement of cardiovascular blood flow. The proposed method consists of combining slice-selective spiral imaging with Fourier velocity-encoded (FVE) imaging. The "spiral FVE" method provides intrinsically higher SNR than PC. Furthermore, it is not affected by partial-volume effects, as it measures the velocity distribution within each voxel. This makes this method particularly useful for assessment of flow in stenotic vessels. We show that spiral FVE measurements of carotid flow show good agreement with Doppler ultrasound measurements. We also introduce a mathematical model for deriving spiral FVE velocity distributions from PC velocity maps, and we use this model to show that spiral FVE provides good quantitative agreement with PC measurements, in an experiment with a pulsatile carotid flow phantom. Then, we show that it is possible to accurately estimate in vivo wall shear rate (WSR) and oscillatory shear index (OSI) at the carotid bifurcation using spiral FVE. Finally, we propose a hybrid MRI/CFD approach that integrates low-resolution PC flow measurements directly into the CFD solver. The feasibility of this MRI-driven CFD approach is demonstrated in the carotid bifurcation of a healthy volunteer. We show that MRI-driven CFD has a regularizing effect on the flow fields obtained with MRI alone, and produces flow patterns that are in better agreement with direct MRI measurements than CFD alone. This approach may provide improved estimates of clinically significant blood flow patterns and derived hemodynamic parameters, such as wall shear stress (WSS).

Magnetic resonance flow imaging 2.1 Basic principles of MRI
MRI is a modality uniquely capable of imaging all aspects of cardiovascular disease, and is a potential "one-stop shop" for cardiovascular health assessment. MRI can generate cross-sectional images in any plane (including oblique planes), and can also measure blood flow. The image acquisition is based on using strong magnetic fields and non-ionizing radiation in the radio frequency range, which are harmless to the patient. The main component of a MRI scanner is a strong magnetic field, called the B 0 field. This magnetic field is always on, even when the scanner is not being used. Typically, MR is used to image hydrogen nuclei, because of its abundance in the human body. Spinning charged particles (or "spins"), such as hydrogen nuclei, act like a tiny bar magnet, presenting a very small magnetic field, emanating from the south pole to the north pole. In normal conditions, each nucleus points to a random direction, resulting in a null net magnetization. However, in the presence of an external magnetic field (such as the B 0 field), they will line up with that field. However, they will not all line up in the same direction. Approximately half will point north, and half will point south. Slightly more than half of these spins (about one in 514 Fluid Dynamics, Computational Modeling and Applications www.intechopen.com Assessment of Carotid Flow Using Magnetic Resonance Imaging and Computational Fluid Dynamics 3 a million) will point north, creating a small net magnetization M 0 , which is strong enough to be detected. The net magnetization is proportional to the strength of the B 0 field, so MRI scanners with stronger magnetic fields (e.g., 3 Tesla) provide higher SNR. Other important components of the scanner are the gradient coils. There are typically three gradient coils (G x , G y , and G z ), that produce an intentional perturbation in the B 0 field when turned on ("played"). This perturbation varies linearly along each spatial direction (x, y and z), such that no perturbation exists at the iso-center of the magnet when these gradients are used. In the presence of an external magnetic field, the spins rotate about the axis of that field. B 0 is (approximately) spatially uniform, so all spins initially rotate at the same frequency (the Larmor frequency), ω = γB 0 , where γ is the gyromagnetic ratio (γ = 42.6 MHz/Tesla for hydrogen protons). However, when any of the gradients is played, the magnetic field becomes spatially varying, and so does the rotation frequency of the spins. Therefore, G x , G y , and G z are used to frequency-encode (or phase-encode) spatial position along the x, y and z directions, respectively. The final major component of the MR scanner is the radio-frequency (RF) coil. This is used to transmit an RF "excitation" pulse to the body, and also to receive the frequency-encoded signal from the "excited" portion of the body. In practice, independent coils may be used for transmission and reception. The RF pulse is typically modulated to the Larmor frequency. While B 0 is aligned with the z-axis (by definition), B 1 , which is a very weak magnetic field associated with the RF pulse, is aligned with the x-axis (also by definition). When the RF pulse is played, some of the spins which are in resonance with the RF pulse (i.e., rotating at the RF pulse's frequency) will now begin to rotate around the x-axis (thus the name magnetic resonance). This tilts the net magnetization towards the x-y plane, and the net magnetization will now have a component in the x-y plane,M xy . The RF pulse is typically designed to have a somewhat rectangular profile in Fourier domain, centered at the modulation frequency (e.g., a modulated windowed sinc). This implies that the RF pulse in fact contains a certain range of frequencies, thus all spins rotating within that range become "excited", or tilted towards the x-axis. So, by playing gradient(s) of appropriate amplitude, and designing the RF pulse accordingly, one can excite only a thin slice of the body, which correspond to the region containing all spins that are in resonance with the RF pulse's range of frequencies. Excitation profiles other than "slices" may also be obtained (e.g., a pencil beam, or cylindrical excitation (Hu et al., 1993)), by designing an appropriate gradient/pulse combination. When the RF pulse is turned off, M xy begins to rotate (at the Larmor frequency) around the z-axis, as the net magnetization begins to realign with B 0 . This rotating magnetization generates an oscillating signal, which can be detected by the receive coil. The frequency content of the received signal can be used to obtain spatial information about the excited portion of the body. In order to frequency-encode spatial information, gradients are also played during signal acquisition. These are called readout gradients. For imaging a slice perpendicular to the z-axis (an axial image), G z is played during excitation (for slice-selection), and G x and G y are played during acquisition. These can be switched, for acquiring sagittal or coronal images, or all three gradients may be used during both excitation and acquisition to image oblique planes. When the readout gradients are played, the acquired signal at a particular time instant corresponds to the sum of sinusoidal signals generated by spins located at different regions of the body, each rotating at different frequencies corresponding to their spatial locations. If 515 Assessment of Carotid Flow Using Magnetic Resonance Imaging and Computational Fluid Dynamics www.intechopen.com an axial slice is being acquired, for example, the demodulated signal value is equivalent to a sample of the Fourier transform M(k x , k y ) of the cross-sectional image m(x, y). In this case, by changing the amplitudes of G x and G y during acquisition, one may acquire different samples of M(k x , k y ). In fact, by playing G x and/or G y , one can move along the k x -k y plane (which is known in MRI as "k-space"), collecting samples of M(k x , k y ). When enough samples of M(k x , k y ) have been collected, an inverse Fourier transform produces m(x, y). The required coverage of k-space, and the number of samples, depend on the specified spatial resolution and field-of-view. For low spatial resolution imaging, only the central portion of k x -k y needs to be sampled. For higher spatial resolution, the periphery of k-space must also be covered. The field of view is associated with the spacing between samples. For a larger field-of-view, k-space needs to be more densely sampled, requiring an increased number of samples. If k-space is not sufficiently sampled, and the resulting field-of-view is not large enough to cover the entire object, overlap in spatial domain (aliasing) is observed. Because signal amplitude is lost as the net magnetization realigns with B 0 (this is called relaxation), multiple acquisitions (excitation + readout) may be needed in order to cover the entire k-space. Some trajectories are more efficient in covering k-space than others. For example, spiral imaging, which uses oscillating gradients to achieve spiral k-space trajectories (Figure 1b), are generally faster than 2DFT imaging, i.e., require fewer acquisitions. In 2DFT imaging, each acquisition readout acquires a single line of k-space, sampling k x -k y in a Cartesian fashion (Figure 1a). This is generally slower, but may be advantageous in some applications with respect to the nature of associated image artifacts. The full sequence of RF pulses and gradients is called a "pulse sequence". The time between acquisitions is called the pulse repetition time, or TR.

Mathematical formalism
As discussed on section 2.1, the acquired MR signal s(t) at a particular time instant corresponds to a sample of the Fourier transform M(k x , k y ) of the excited magnetization m(x, y): x y m(x, y) e −j2π(k x x+k y y) dx dy. (1) The Fourier coordinates k x and k y vary with time, according to the zeroth moment of the readout gradients G x and G y : These equations explain how the gradients can be used to "move" along k-space, as discussed in section 2.1. This formalism can be generalized for any combination of G x , G y , and G z gradients: where G r is the oblique gradient resulting from the combination of the G x , G y and G z gradients, and r is its corresponding axis along which the linear variation in magnetic field intensity is realized. Given a spatial position function r(t) and a magnetic field gradient G r (t), the magnetization phase is: For static spins, r(t) is constant ( r), and this becomes: as in the exponential in equation 4.

Principles of MR flow imaging
The basic principles of quantitative flow measurement using magnetic resonance were first proposed by Singer (1959) and Hahn (1960) in the late 1950s. However, clinical applications of MR flow quantification were not reported until the early 1980s (Moran et al., 1985;Nayler et al., 1986;Singer & Crooks, 1983;van Dijk, 1984 velocity times the first moment of the gradient waveform along the direction in which they are moving. For spins moving along the r-axis with a constant velocity v, and initial position r 0 , we can write r(t)= r 0 + vt. Rewriting equation 6, for t = t 0 : where M 0 and M 1 are the zeroth and first moments of the r-gradient waveform at the time of signal acquisitions ("echo time", or "time to echo" (TE)), respectively. Thus, if a gradient with null zeroth moment is used (e.g., a bipolar gradient, aligned with v), the phase accrued for a constant velocity spin is φ = γ v · M 1 . Therefore, if a bipolar gradient waveform is played between the excitation and the readout, the phase measured in a pixel of the acquired image is directly proportional to the velocity of the spins contained within its corresponding voxel. However, factors other than flow (such as inhomogeneities of the magnetic field) may cause additional phase shifts that would cause erroneous interpretation of the local velocity (Rebergen et al., 1993).

Phase contrast
The phase contrast method addresses the problem mentioned above by using two gradient-echo data acquisitions in which the first moment of the bipolar gradient waveform is varied between measurements (O'Donnell, 1985). The velocity in each voxel is measured as: where φ a (x, y) and φ b (x, y) are the phase images acquired in each acquisition, and M a 1 and M b 1 are the first moment of the bipolar gradients used in each acquisition.

Fourier velocity encoding
While phase contrast provides a single velocity measurement associated with each voxel, Fourier velocity encoding (Moran, 1982) provides a velocity histogram for each spatial location, which is a measurement of the velocity distribution within each voxel. FVE involves phase-encoding along a velocity dimension. Instead of only two acquisitions, as in phase contrast, multiple acquisitions are performed, and a bipolar gradient with a different amplitude (and first moment) is used in each acquisition. Equation 10 can be rewritten as: where k v is the velocity frequency variable associated with v, and is proportional to the first moment of G r (t): Each voxel of the two-dimensional image is associated with a distribution of velocities. This three-dimensional function m(x, y, v) is associated with a three-dimensional Fourier space Thus, an extra dimension is added to k-space, and multiple acquisitions are required to cover the entire k x -k y -k v space. In order to move along k v , a bipolar gradient with the appropriate amplitude (and first moment) is played before the k x -k y readout gradients, in each acquisition. Placing the bipolar gradient along the z-axis will encode through-plane velocities. Placing the bipolar gradient along x or y will encode in-plane velocities. Oblique flow can be encoded using a combination of bipolar gradients along the x, y and z axes. Each acquisition along k v is called a velocity encode. The number of required velocity encodes depends on the desired velocity resolution and velocity field-of-view (the maximum range of velocities measured without aliasing). For example, to obtain a 25 cm/s resolution over a 600 cm/s field-of-view, 24 velocity encodes are needed. The spatial-velocity distribution, m(x, y, v), is obtained by inverse Fourier transforming the acquired data, M(k x , k y , k v ). If cine imaging (Glover & Pelc, 1988) is used, measurements are also time resolved, resulting in a four-dimensional dataset: m(x, y, v, t).

FVE signal model
2DFT phase contrast provides two two-dimensional functions, m(x, y) and v o (x, y), the magnitude and velocity maps, respectively. If these maps are measured with sufficiently high spatial resolution, and flow is laminar, one can assume that each voxel contains only one velocity, and therefore the spatial-velocity distribution associated with the object is approximately: where δ(v) is the Dirac delta function. In 2DFT FVE, k-space data is truncated to a rectangular cuboid in k x -k y -k v space. The associated object domain spatial-velocity blurring can be modeled as a convolution of the true object distribution, s(x, y, v), with sinc(x/∆x), sinc(y/∆y), and sinc(v/∆v), where ∆x and ∆y are the spatial resolutions along the x and y axes, respectively, and ∆v is the velocity resolution, as follows: whereŝ(x, y, v) is the measured object distribution, and * denotes convolution. This is equivalent to: This approach for deriving FVE data from high-resolution velocity maps can be used for many simulation purposes (Carvalho et al., 2010).

Slice-selective FVE with spiral readouts (spiral FVE)
Phase contrast imaging is fast, but has limitations associated with partial-volume effects. On the other hand, 2DFT FVE addresses those limitations, but requires long scan times. Thus, we propose the use of slice-selective FVE MRI with spiral acquisitions. The proposed spiral FVE method is capable of acquiring fully localized, time-resolved velocity distributions in a short breath-hold. Scan-plane prescription is performed using classic protocols. We present practical implementations for measuring blood flow through the carotid arteries, and comparisons with Doppler ultrasound and high-resolution 2DFT phase contrast MRI. The proposed method is demonstrated in healthy volunteers. Subjects provided informed consent, and were imaged using a protocol approved by the institutional review board of the University of Southern California.

Pulse sequence
The spiral FVE imaging pulse sequence ( Figure 2) consists of a slice-selective excitation, a velocity-encoding bipolar gradient, a spiral readout, and refocusing and spoiling gradients.
The dataset corresponding to each temporal frame is a stack-of-spirals in k x -k y -k v space ( Figure 3). The bipolar gradient effectively phase-encodes in k v , while each spiral readout acquires one "disc" in k x -k y .

Spiral FVE signal model
As spiral FVE acquisitions follow a stack-of-spirals pattern in k x -k y -k v space (Figure 3), k-space data is truncated to a cylinder, i.e., a circle along k x -k y (with diameter 1/∆r), and a rect function along k v (with width 1/∆v), where ∆r and ∆v are the prescribed spatial and velocity resolutions, respectively. Using the same approach we used in section 2.6 for 2DFT FVE, the associated object domain spatial-velocity blurring in spiral FVE can be modeled as a convolution of the true object distribution, s(x, y, v), with jinc( x 2 + y 2 /∆r) and sinc(v/∆v), resulting in:

In vitro validation
An in vitro comparison of velocity distributions measured with spiral FVE with those derived from high-resolution 2DFT phase contrast -the current MR gold standard -was performed.
The signal model presented in section 3.2 was used to generate simulated FVE data based on high-resolution 2DFT phase contrast data.
The validation experiments were performed using a pulsatile carotid flow phantom (Phantoms by Design, Inc., Bothell, WA). A slice perpendicular to the phantom's carotid bifurcation was prescribed, and through-plane velocities were measured. A cine gradient-echo 2DFT phase contrast sequence with high spatial resolution and high SNR (0.33 mm resolution, 10 averages, 80 cm/s Venc) was used as a reference. Cine spiral FVE data with ∆r = 3m ma n d∆v = 10 cm/s was obtained from the same scan plane. Both acquisitions were prospectively gated, and used the same TR (11.6 ms), flip angle (30 • ), slice profile (3 mm), temporal resolution (23.2 ms), and pre-scan settings. The total scan time was 40 minutes for phase contrast, and 12 seconds for FVE. A simulated spiral FVE dataset was computed from the PC magnitude and velocity maps, using the convolution model described in Eq. 18. The PC-derived and FVE-measured data were registered by taking one magnitude image m(x, y) from each dataset, and then using the phase difference between their Fourier transforms M(k x , k y ) to estimate the spatial shift between the images. Amplitude scaling was performed by normalizing the ℓ 2 -norm of each FVE dataset. The difference between PC-derived and FVE-measured time-velocity distributions was calculated for select voxels, and the associated signal-to-error ratios were computed. This was used as a quantitative assessment of spiral FVE's accuracy. Figure 4 shows measured and PC-derived time-velocity FVE distributions from nine representative voxels, selected around the circumference of the vessel wall of the pulsatile carotid flow phantom's bifurcation. The signal-to-error ratio between measured and PC-derived time-velocity distributions was measured to be within 9.3-11.7 dB. Imperfect registration between the datasets, combined with spatial blurring due to off-resonance in the measured spiral FVE data, may have contributed to this moderate signal-to-error ratio. Nevertheless, the two datasets show good visual agreement, and no significant spatial variation was observed in terms of accuracy. These results show that velocity distributions

In vivo evaluation
The spiral FVE method was evaluated in vivo, aiming at quantifying flow through the common carotid artery of a healthy volunteer. Doppler ultrasound was used as a gold standard and was qualitatively compared with the proposed method. We also show experiments for evaluating the most appropriate view-ordering scheme for measuring carotid flow, and we demonstrate spiral FVE's ability to measure multiple flows in a single acquisition.

View ordering
Experiments were performed to determine the appropriate view-ordering scheme for carotid spiral FVE studies. Flow was measured through the carotid artery of a healthy volunteer, using a 4-interleaf spiral FVE acquisition. The measurement was performed twice, and in each acquisition, a different view-ordering scheme was used. In the first acquisition, one spiral interleaf was acquired per heartbeat, and two different k v levels were encoded throughout each R-R interval. In the second acquisition, two spiral interleaves were acquired per heartbeat, encoding one k v level per cardiac cycle. Reconstructed velocity histograms were compared with respect to data inconsistency artifacts related to view-ordering. The results are shown in Figure 5. Note that the ghosting artifacts that arise in the velocity histograms when acquiring two different velocity encode levels during the same heartbeat ( Figure 5a) do not appear when we used interleaf segmentation instead (Figure 5b). In contrast to view-sharing along k v , which causes ghosting in the velocity direction, view-sharing among spiral interleaves introduces swirling artifacts in image domain, reducing the effective unaliased spatial field-of-view by a factor of 2. However, only moving spins (flowing blood) experience these artifacts, and vessels on the same side of the neck are relatively close to each other. Therefore, the unaliased field-of-view is wide enough to enclose all vessels on one side of the neck, so that quantification of a vessel of interest will not be disturbed by flow from neighboring vessels (e.g. measurement of flow in the left carotid arteries will not be disturbed by flow in the left jugular vein). In order to suppress signal from the opposite side of the neck, we separately reconstruct data from the left and right neck-coil elements. This uses the receiver coil sensitivity profile to avoid spiral view-sharing artifacts. Based on these observations, we used spiral interleaf view-sharing in all subsequent carotid studies. When two or more k v levels are acquired during the same heartbeat (a), velocity distribution changes between consecutive TRs cause ghosting artifacts along the velocity axis (arrow). This artifact is not seen if, in the same heartbeat, different spiral interleaves, but only one k v encoding, are acquired (b).

Comparison with Doppler ultrasound
Representative in vivo results are compared with Doppler ultrasound in Figure 6. The MRI measured time-velocity histogram shows good agreement with the ultrasound measurement, as the peak velocity and the shape of the flow waveform were comparable to those observed in the ultrasound studies.

Wall shear rate estimation using spiral FVE
Arterial wall shear stress, the drag force acting on the endothelium as a result of blood flow, is widely believed to influence the formation and growth of atherosclerotic plaque, and may have prognostic value. According to Newton's law of viscosity, WSS can be estimated as the product of wall shear rate and blood viscosity (µ), where WSR is the radial gradient of blood flow velocity (dv/dr) at the vessel wall. Low WSS (Tang et al., 2008;Zarins et al., 1983) and highly oscillatory WSS (Ku et al., 1985) have been linked to the formation and growth of atherosclerotic plaques, and this link has been validated in vitro (Dai et al., 2004). High WSS has also been hypothesized as a factor responsible for the topography of atherosclerotic lesions (Thubrikar & Robicsek, 1995). Phase contrast MRI (PC-MRI) suffers from data inconsistency, partial-volume effects (Tang et al., 1993), intravoxel phase dispersion and inadequate SNR at high spatial resolutions ( Figure 8). PC-MRI is not currently capable of providing accurate absolute measurements of WSS (Boussel et al., 2009), and has been shown to underestimate blood flow velocities in the carotid bifurcation by 31-39% . As shown in section 3, spiral FVE produces accurate velocity histograms compared with those acquired with Doppler ultrasound. In addition, spiral FVE method is capable of rapid acquisition of fully localized, time-resolved velocity distributions. In this section, we propose High spatial resolution is typically associated with low SNR (a). Averaging multiple acquisitions improves SNR (b), but also increases the total scan time, and may cause loss of effective resolution due to subject motion. Scan parameters: 0.33x0.33x3 mm 3 spatial resolution, 37 ms temporal resolution, 30 • flip angle, 80 cm/s Venc, 2-minute scan (120 heartbeats) per acquisition. the use of spiral FVE to estimate WSR and the oscillatory shear index. This is made possible by the method proposed by Frayne & Rutt (1995) for FVE-based WSR estimation, which is discussed next.

FVE-based WSR estimation: the Frayne method
The Frayne method for FVE-based WSR estimation (Frayne & Rutt, 1995) involves obtaining the velocity distribution for a voxel spanning the blood/vessel wall interface, and then using this distribution to reconstruct the velocity profile across the voxel, with sub-voxel spatial resolution. Assuming that signal intensity is spatially and velocity invariant, the sum of the signals from all velocities within a voxel is proportional to the total volume of material within the voxel. Furthermore, two assumptions can be made about the shape of the velocity profile: (i) the fluid velocity at the vessel wall is approximately zero, and (ii) the velocity profile within a voxel is monotonically increasing or decreasing. Using these assumptions, a step-wise discrete approximationṽ(r) to the true velocity profile across a voxel can be obtained from its velocity distribution s(v) by inverting the discrete function r(v), which is constructed as follows: Note that, for each velocity bin v i , the intra-voxel position r is incremented by a fraction of the total radial extent of the voxel (∆r). This fraction is proportional to h(v i ), which is the volume fraction of the voxel that has velocity v = v i . The volume fractions are 525 Assessment of Carotid Flow Using Magnetic Resonance Imaging and Computational Fluid Dynamics www.intechopen.com calculated by normalizing the velocity distribution. This process is demonstrated graphically in Figure 9. Spatial variations in signal intensity due to radio-frequency saturation (i.e., inflow enhancement) and due to differences in 1 H density and relaxation properties between vessel wall tissue and blood may be compensated by adjusting s(v) accordingly, prior to calculating h(v) (Frayne & Rutt, 1995).  (Frayne & Rutt, 1995). WSR is estimated from FVE velocity distributions in voxels spanning the blood/vessel wall interface. First, a threshold is applied to the velocity histogram to reduce noise sensitivity (a). Then, the volume fraction within each velocity bin is converted into a radial position across the voxel, using Eq. 19. Finally, the velocity gradient is estimated from the reconstructed velocity profile (b), within a small velocity interval [v 0 ,v 1 ].
In order to reduce the effects of ringing and noise rectification due to the magnitude operation in Eq. 19, a threshold is applied to s(v) before normalization (Figure 9a). The appropriate threshold level must be determined by analyzing the signal intensities in the velocity distribution for a range of velocities outside the range of expected blood flow velocities. It is assumed that signal outside this expected range is exclusively due to rectified noise (Frayne & Rutt, 1995). In our implementation, only components that are below the specified threshold and outside this expected range of velocities are set to zero. The WSR can be estimated by prescribing a velocity interval [v 0 ,v 1 ] and then fitting a first-order polynomial to the points ofṽ(r) within this interval ( Figure 9b). Ideally, v 0 = 0 and v 1 = ∆v, because we wish to estimate the velocity derivative at the blood-wall interface. The SNR of shear rate estimates will increase as this velocity interval becomes larger, because of averaging across multiple velocity steps. However, as the interval becomes larger, the shear rate is averaged over a larger distance within the voxel and may deviate from the true local shear at the wall (Frayne & Rutt, 1995). Therefore, it is important to prescribe a reasonable [v 0 ,v 1 ] interval. In our implementation, v 0 = ∆v and v 1 ≈ 30 cm/s are used for an initial assessment, and then the interval is manually adjusted for selected voxels of interest. The same voxel-based approach is used with respect to the noise threshold discussed above, with a fixed threshold value being used for the initial assessment.

In vivo WSR measurement
The in vivo measurement of carotid WSR using spiral FVE acquisitions with Frayne's reconstruction is now demonstrated. Three healthy subjects were studied. For each volunteer, five 5 mm contiguous slices (2.5 cm coverage) were prescribed perpendicular to the left carotid bifurcation. Each slice was imaged independently (separate acquisitions). Localized gradient shimming was performed, and acquisitions were prospectively ECG-gated. Several cardiac phases were acquired, spanning the systolic portion of the cardiac cycle. A time-bandwidth product 2 RF pulse was used for excitation, and the flip angle was 30 • . Only through-plane velocities were measured, using 32 k v encoding steps over a 160 cm/s velocity field-of-view (5 cm/s resolution). The velocity field-of-view was shifted from the −80 to 80 cm/s range to the −40 to 120 cm/s range during reconstruction, in order to avoid aliasing and maximize field-of-view usage. Negative velocities are encoded in order to assess negative WSR values, and also to accommodate leakage and ringing due to k v truncation (i.e., finite velocity resolution). Spatial encoding was performed using eight 4 ms variable-density spiral interleaves (16∼4 cm field-of-view, 1.4 mm resolution). The temporal resolution was 24 ms (12 ms TR, 2 views per beat). Scan time was 128 heartbeats per slice, i.e., approximately 2 minutes per slice. The subjects provided informed consent, and were scanned using a protocol approved by the institutional review board of the University of Southern California. Figures 10 and 11 show two representative sets of in vivo results. The WSR values are shown for manually-segmented regions-of-interest. Figure 10 illustrates the variation in WSR along all three spatial dimensions near the carotid bifurcation of subject #1. These results correspond to the cardiac phase with the highest peak velocity. Figure 11 illustrates the temporal variation of pulsatile WSR in the common carotid artery of subject #2. The results show a circumferential variation in WSR around the wall of the common carotid artery. Markl et al. (2009) recently observed a similar variation using a PC-based approach.

Estimation of oscillatory shear index
The OSI is important for the evaluation of shear stress imposed by pulsatile flow. This index describes the shear stress acting in directions other than the direction of the temporal mean shear stress vector (He & Ku, 1996;Ku et al., 1985). It is defined as the relation between the time-integral of the shear stress component acting in the direction opposite to the main direction of flow and the time-integral of the absolute shear stress. In this work, OSI was calculated as defined by He & Ku (1996), and assuming spatially and temporally invariant blood viscosity. Measuring negative WSR with the proposed method requires voxels to be small enough to contain only reverse flow. Voxels containing both forward and reverse flow would violate the assumption of a monotonically increasing/decreasing velocity profile within the voxel. Under sufficient spatial resolution, negative WSR may be measured simply by using a negative [v 0 ,v 1 ] interval. Figure 12 presents a demonstration of in vivo OSI estimation, using the proposed method. Results are shown for eight representative voxels, selected around the circumference of the wall of the carotid bifurcation of subject #3. Measured velocity distributions, wall shear rates, and OSI values, corresponding to the systolic portion of the cardiac cycle, are shown for each voxel. High OSI values were observed at the wall corresponding to the carotid bulb. Non-zero OSI was also observed at the opposite wall. These findings are in agreement with PC-based OSI measurements recently reported by Markl et al. (2009). The results also suggest   Fig. 11. Temporal variation of pulsatile WSR measured in the common carotid artery of subject #2. Results correspond to the cardiac phases acquired 48-168 ms after the ECG trigger, and to a slice prescribed 10 mm below the carotid bifurcation.

528
Fluid Dynamics, Computational Modeling and Applications www.intechopen.com that higher spatial resolution is needed for accurately estimating OSI in some of the voxels. Notably, voxel (h) presents both positive and negative velocity components simultaneously during post-systolic deceleration. This indicates that the spatial resolution was insufficient, and the assumption of a monotonically decreasing velocity profile within the voxel was violated. OSI=0

Assessment of carotid flow using CFD (and MRI)
Computational fluid dynamics methods are concerned with the approximate solution of the fluid motion as well as with the interaction of the fluid with solid bodies. Fluid dynamic problems, essentially, are based on nonlinear systems of coupled partial differential equations. Attempting to solve those systems analytically, in general, is an impracticable task. Over the years, researchers have developed methodologies, codes, schemes and algorithms to find approximate solutions, focusing on accuracy and speed of convergence, for problems involving fluid dynamics and heat transfer. Historically, CFD started in the early 1970s, triggered by the availability of increasingly more powerful mainframes. In the beginning, the study was limited to high-technology engineering areas of aeronautics and astronautics. Nowadays, computational fluid dynamics methodologies are routinely employed in many fields such as: racing car design, ship design, meteorology, oil recovery, civil engineering, airspace engineering, and biomedical engineering. This section is dedicated to the use of CFD for evaluation of blood flow in the carotid arteries. In vivo 3D blood flow patterns can be either measured directly using PC-MRI, or obtained from model-based CFD calculations. PC-MRI is accurate, but suffers from long scan times and limited spatio-temporal resolution and SNR. CFD provides arbitrarily high resolution and reduced scan times, but its accuracy hinges on the model assumptions. A numerical framework for constructing a flow field that is influenced by both direct measurements and a fluid physics model is described.

CFD in biomedical engineering
Recently, CFD is playing an important role in the analysis of blood flow. This technique of flow visualization has been widely applied in problems involving arterial diseases with reconstruction of hemodynamics in realistic models based on images generated by standard in vivo medical visualization tools, such as MRI. CFD can be used to improve the data obtained using MRI for estimation of the flow properties of blood vessels. CFD can also be used to compare real data obtained from patients with simulated data models using realistic geometries. These geometries may even be constructed using MRI data. Boussel et al. (2009) compared a time-dependent 3D phase-contrast MRI sequence with patient-specific CFD models for patients who had intracranial aneurysms. The evolution of intracranial aneurysms is known to be related to hemodynamic forces such as wall shear stress and maximum shear stress. Harloff et al. (2010) used CFD image-based modeling as an option to improve the accuracy of MRI-based WSS and OSI estimation, with the purpose of correlating atherogenic low WSS and high OSI with the localization of aortic plaques. Steinman et al. (2002) used a novel approach for noninvasively reconstructing artery wall thickness and local hemodynamics at the human carotid bifurcation. Three-dimensional models of the lumen and wall boundaries, from which wall thickness can be measured, were reconstructed from black blood MRI. Along with time-varying inlet/outlet flow rates measured via PC-MRI, the lumen boundary was used as input for a CFD simulation of the subject-specific flow patterns and wall shear stress. Long et al. (2003) were concerned with the reproducibility of geometry reconstruction, one of the most crucial steps in the modeling process. Canstein et al. (2008) used rapid prototyping to transform aortic geometries as measured by contrast-enhanced MR angiography into realistic vascular models with large anatomical coverage. Visualization of characteristic 3D flow patterns and quantitative comparisons of the in vitro experiments with in vivo data and CFD simulations in identical vascular geometries were performed to evaluate the accuracy of vascular model systems. CFD can also be used with other imaging techniques such as tomography. Howell et al. (2007) used CFD to study the temporal and spatial variations in surface pressure and shear through the cardiac cycle on models of bifurcated stent-grafts derived from computed tomography in patients who had previously undergone endovascular repair of abdominal aortic aneurysm.

MRI-driven CFD
A numerical framework for constructing a flow field that is influenced by both direct measurements and a fluid physics model is now described. The PC-MRI signal is expressed as a linear function of the velocities on the underlying computational grid, and a tunable parameter controls the relative influence of direct measurements and the model assumptions. volunteer. The results show that this methodology produces flow fields that are in better agreement with direct PC-MRI measurements than CFD alone.

Proposed approach
Here, blood is modeled as an incompressible Newtonian fluid with constant viscosity µ. This model is widely used in in vivo CFD analysis, and underlies the commonly held definition of wall shear stress as being proportional to dv/dr, where v is the velocity tangential to the vessel wall, and r is the perpendicular distance from the wall. The task of a CFD routine is to solve the momentum and continuity equations that the flow field -subject to the assumption of Newtonian flow -must obey. Following the control-volume formulation introduced by Patankar (1980), the momentum equation is where ρ is the fluid density, v =(u, v, w) is the velocity vector field, p is the pressure field, ∇ is the gradient operator, and ∆ is the Laplacian operator. The flow field must also satisfy mass conservation (or continuity), which can be expressed as Equations (20) and (21) must be solved for the unknown scalar field variables u, v, w, and p. These equations are non-linear and coupled, and attempting to solve them directly in one step is a formidable, if not impossible, task. The solver was built on the SIMPLER algorithm developed by Patankar (1980), which is a well-known and established numerical routine for solving Eqs. (20) and (21). SIMPLER starts with an initial estimate for the velocity field, and updates this estimate in an iterative fashion. At each iteration i, the discretized equations are linearized using the velocity estimate v i−1 at the previous iteration, which produces a square system matrix A for each velocity component. For example, for the z velocity component, we can write where w is the z velocity in all voxels in the 3D calculation domain (stacked to form a 1D column vector), and b is a constant column vector. Note that the pressure field is updated periodically based on the current velocity field estimate, and hence does not appear explicitly in (22).
The key step in our approach is to add additional rows to A and b that incorporate MRI measurements of one or more velocity components. The underlying assumption we make here is that the velocity measured with PC-MRI is equal to the average velocity within the voxel, and can hence be expressed as a linear combination of the velocities on the underlying CFD calculation grid. For example, for the z (S/I) velocity component we have  (22) and (23), we have The weighting parameter s is an adjustable scalar that determines the influence of the MRI measurements on the solution. For example, s = 0 produces a conventional, unconstrained CFD solution. Here, we solve Eq. (24) in the least squares sense, using the conjugate gradient method on the normal equations. Hence, the solution w at each iteration represents a weighted least squares estimate, with relative contributions of fluid physics and direct MRI measurement being controlled by the parameters after multiple interactions. The velocity field v =(u, v, w) converges toward the solution for a particular point t. Time is then incremented by an amount ∆t, and the iterative procedure is repeated to obtain the solution at t + ∆t. This way, it is possible to predict a time-dependent flow field, as long as the true velocity field v at t = 0 is known. In this work, however, the exact velocity field is not known, since only a relatively low-resolution and noisy PC-MRI velocity field is available. Instead, we will calculate the steady flow v s subject to the measured inlet and outlet velocities. We will furthermore assume that the velocity field obtained with PC-MRI at one time-point (near peak flow) is representative of the steady flow given the inlet and outlet velocities at that time-point. We obtain v s by starting with an initial guess for v, and carrying the simulations forward in time until convergence.

In vivo demonstration
PC-MRI data were obtained from four time-resolved 3DFT FGRE image volumes in the carotid artery in one healthy volunteer (1×1×2.5 mm 3 voxel size; FOV 16×12×7.5 cm 3 ; TR 7.0 ms; flip angle 15 • ; temporal resolution 56 ms; Venc 1.6 m/s; 7 minutes per scan), on a GE Signa 3T EXCITE HD system (4 G/cm and 15 G/cm/ms maximum gradient amplitude and slew rate) with a 4-channel carotid receive coil array. The through-slab (z) axis was oriented along the S/I direction. The subject provided informed consent, and was scanned using a protocol approved by the institutional review board of the University of Southern California. PC-MRI velocity component maps u MR , v MR , and w MR were calculated using data from one receive coil. Residual linear velocity offsets in each velocity component map (e.g. due to eddy-currents) were removed by performing a linear fit to manually defined 3D regions containing only stationary tissue. The vessel lumen was segmented by manually outlining the vessel borders from a stack of 2D axial slices. The solver calculations assumed a blood viscosity of 0.0027 Pa·sec and a blood density of 1060 kg/m 3 (Reynolds number of order 1000), and no-slip boundary conditions. Calculations were performed on a Cartesian grid of 1 mm isotropic resolution, and all algorithms were implemented in Matlab. Low-resolution (truncated to 1×1×3mm 3 ) measurements of the S/I velocity component w MR near the time-point of peak flow were incorporated into the solver, and a steady velocity vector field v s was calculated as described above. Hence, w was partially constrained by the measured velocity component w MR , whereas u and v were determined solely from the fluid physics model. As in the original SIMPLER algorithm (Patankar, 1980), u, v, and w were defined on regular grids that were staggered by half a grid spacing (in different directions) with respect to the centered grid. This is done to avoid a checkerboard solution for the pressure and velocity fields (Patankar, 1980 Figure 13 compares flow fields in the carotid artery obtained with PC-MRI only, CFD only (s = 0), and the proposed combined solver algorithm with s = 5. These flow fields were calculated from 4, 1, or 2 time-resolved 3D image acquisitions, corresponding to a total scan time of 28, 7, and 14 minutes, as indicated in the figure. In the common carotid artery, flow is predominantly along z, and all methods produce comparable flow fields. In the bifurcation, the bulk flow pattern appears qualitatively similar for all three methods, which indicates that the underlying fluid physics model makes reasonable predictions regarding the transverse velocities. However, comparison of the CFD and PC-MRI results shows that CFD underestimates the velocities in the external carotid artery. The combined solver, which strikes a compromise between CFD and PC-MRI, brings the calculated flow field in closer agreement with the values measured with PC-MRI. These velocity fields were obtained from 4 (left), 1 (center), or 2 (right) MRI scans, corresponding to total scan times of 28, 7 and 14 minutes, respectively. Compared to CFD, the combined solver produces flow fields that are in better qualitative and quantitative agreement with PC-MRI.

Conclusion
In this chapter, we have introduced spiral FVE, a rapid MRI method for fully-localized measurement of cardiovascular blood flow. The proposed method was shown to be capable of measuring blood flow in the carotid arteries and estimating wall shear stress and oscillatory shear index at the carotid bifurcation. In addition, we proposed a combined CFD-MRI solver that integrates the non-linear coupled system of the fluid partial differential equations using MRI data as initial data. This methodology, which has a tunable parameter, is capable of The content of this book covers several up-to-date topics in fluid dynamics, computational modeling and its applications, and it is intended to serve as a general reference for scientists, engineers, and graduate students. The book is comprised of 30 chapters divided into 5 parts, which include: winds, building and risk prevention; multiphase flow, structures and gases; heat transfer, combustion and energy; medical and biomechanical applications; and other important themes. This book also provides a comprehensive overview of computational fluid dynamics and applications, without excluding experimental and theoretical aspects.

How to reference
In order to correctly reference this scholarly work, feel free to copy and paste the following: